基于边界积分法的舰艇磁场模型参数确定

基于边界积分法的舰艇磁场模型参数确定

0重点难以确定的磁场范围难以延拓的归位方法或原理消磁是提高潜水员耐磁性素质的有效手段。评价消磁效果或者评估潜艇的磁性防护能力,一般都是针对潜艇水线以下的某种特定深度上某些特定点的磁感应强度或潜艇上方的特定高度上的磁信号来进行的。受实际测量条件(水深,风浪和高度等)的限制,在实际测量中往往是很难直接获得这类特定的磁场信息,所以需要从已知场点的磁场值来推算考核点上的磁场值。目前,用于舰艇磁场延拓的方法主要有有限元法、积分方程法和磁体模拟法。有限元和积分方程法都需要对舰船进行剖分,其剖分的合理性和疏密程度严重影响着延拓结果。对于形体复杂且庞大的舰船而言,剖分是非常复杂困难的,并且得到的方程组也将非常庞大,运算时间较长。磁体模拟法常用的磁性模拟体主要有磁偶极子和磁椭球体,该方法易受人为经验的影响,稳定性和精度难以保证。人们在20世纪70年代将边界积分方程用于对位场的研究,由于该方法只对求解区域边界剖分,因而使问题的维数降低,这就使代数方程组的维数大为减少,节约了计算成本。1磁体扩展数学模型1.1各区域v内磁场的各磁感应强度计算如图1所示,闭合面S将场域空间分成2部分V1和V,所有的磁源(包括铁磁物质和传导电流)都被限制在V1内。V可以看成是由2个闭合面限定的区域(从内部限定的闭合面S和伸展到无限远的S∞),该区域内没有电流分布(J=0),标量磁位φ满足拉普拉斯方程∇∇2φ=0.由于磁位φP在无穷远处定义为0,则直接边界积分方程可以写成:cPφP=∮S14πr∂φ∂ndS-14πr∂φ∂ndS−∮Sφ∂∂nφ∂∂n14πr14πrdS,(1)其中:cP对V内的点为1,对光滑边界S上的点为1/2,对区域V1内的点为0;r是边界面上任意点到场点的矢径;S面上的外法线n从区域V指向V1.记φ*=∂φ∂nφ∗=∂φ∂n,则从(1)式可以看出:V1内的磁源对V中任意位置磁位φP的贡献可以用边界上的面积分来表示。若已知边界面S上的等效源磁性参数(φ和φ*)后,可以由该式计算V中任何一点的标量磁位φP.在得到该点的φP后,通过求取该点处的负梯度,并结合该点处的磁导率μ,即可得到该点处的磁感应强度。BP=-μ∇∇´φΡ=-μ4π´φP=−μ4π∮S∇∇′1r∂φ∂n-φ∂∂n1r∂φ∂n−φ∂∂n1r1rdS,(2)其中,法线方向的偏微分是对边界源点坐标进行的,而梯度∇∇′是对场点坐标进行的。将边界面S划分为N个面积单元,则两式可离散化为这就是边界上的等效源磁性参数与无磁区域V中的磁感应强度的正演关系。如果在区域V内测量得到一组潜艇磁场感应强度,就可以利用(4)式反演出边界上等效的磁性参数,再利用这组等效的磁性参数可以计算V区域内任意位置的磁感应强度。在反演的过程中,还需要考虑一些其他信息,如边界面磁通连续,φ和φ*的平衡,以及两者在边界上的连续性。1.2单元中心处的磁位平衡从测量值反演等效磁性参数时,约束条件为:1)φ和φ*的平衡;2)磁通连续。如果将场点取在各个单元的中心,则cp=1/2(为了计算方便,一般取光滑的边界),φ和φ*在边界上的平衡,可以描述为φi为第i个单元中心处的磁位。由于所有磁源都被限制在闭合面S内,进出S的磁力线条数相等。由高斯定理∮SBdS=0,经推导并离散化后得到Ν∑j=1∫Sj∑j=1N∫Sj∂φj∂ndSj=0.(6)∂φj∂ndSj=0.(6)1.3型磁管磁模型的正演分析边界形状和边界单元的划分严重影响着正演和反演计算的速度和精度。出于单元划分和积分计算简单考虑,一般都选取较规则的边界(如长方体表面,图2).由于边界S上没有实际电荷和磁源,所以边界上的φ和φ*连续;同时边界S存在拐点(如长方体的棱),在这些拐点处φ*的方向未定。线性单元适合于φ和φ*在边界上连续的场合,但是如果将边界全部使用线性单元分割,解在拐点处的误差就大。为了改善解的精度,本文将边界元视为混合单元,即在每个单元内φ是连续线性变化的,而φ*是定常的,此时(4)式和(5)式相应的变为式中:Q为单元节点数;ωkiki和ωkjkj为加权系数;φkiki和φkjkj为单元i和j的第k个节点处的磁位。ωkiki和Q由边界面的具体形状和分割方式决定,一般采用简单的四边形元素(本文中为长方形,见图2),它由4个节点组成,此时Q=4,ωkiki可按下式计算:{ω1i=(1-ε1)(1-ε2)4,ω2i=(1+ε1)(1-ε2)4,ε3i=(1+ε1)(1+ε2)4,ε4i=(1-ε1)(1+ε2)4,-1≤ε1,ε2≤1.(9)⎧⎩⎨⎪⎪⎪⎪⎪⎪ω1i=(1−ε1)(1−ε2)4,ω2i=(1+ε1)(1−ε2)4,ε3i=(1+ε1)(1+ε2)4,ε4i=(1−ε1)(1+ε2)4,−1≤ε1,ε2≤1.(9)(7)式和(8)式中的积分项与场点Pi和边界单元Sj有关,积分方法与单元的具体形状密切相关,一般难以求得解析表达式,都采用数值积分方法。积分项的具体表达式可以参阅文献,这里不再赘述。对于(8)式中,应当注意的是:当i=j时,源点与节点在同一单元上,出现积分奇异,两者的积分结果如下:(11)式只适用于图2中的长方形边界单元,其中l和d为长方形的长和宽。至此,作为延拓潜艇磁感应强度所需的第1个步骤,我们完成了正演问题的分析。延拓的具体过程可以描述为:在潜艇周围空间测量得到m个不同位置的潜艇磁场三分量,则综合(6)～(8)式形成方程组,写成矩阵的形式为:Ax=b;求解方程组Ax=b反演得到边界上一组等效的磁性参数后,再利用(7)式计算舰船在其他位置的磁场感应强度。2模型求解和正则化参数的选取由上节模型得到的方程组都呈现严重的病态特性,磁感应强度测量值的微小扰动会导致解的很大变化。而在求解方程组时,通常按线性最小二乘来求解,即求解法方程组ATAx=ATb.(12)若方程本身是病态的,则法方程会更加的病态。国内外很多学者都对这类病态问题进行了研究,这些研究表明TSVD和Tikhonov正则化等方法可以较好的处理逆问题中由系数矩阵病态特性而造成结果的不稳定性。本文将采用Tikhonov正则化方法对方程组进行求解。2.1Tikhonov正则化基本原理不妨假设实际磁场测量过程中测量误差为ξ,理想情况下磁场测量值为c,则有b=c+ξ,Ax=c+ξ.对系数矩阵进行奇异值分解A=USVT,则(12)式的解为x=k∑i=1∑i=1k(uΤic)σivi+(uTic)σivi+k∑i=1∑i=1k(uΤiξ)σivi,(13)其中ui和vi分别为U和V矩阵的第i列,σ1≥σ2≥…≥σk为A的奇异值,k=rank(A).分析(13)式可知,噪声对真解的污染主要反映在右端第2项,当奇异值σi异常小时,该项将被严重放大,甚至掩盖真解(右端第1项)。TSVD方法把容易造成不稳定的较小奇异值直接截去,反映在物理实际中就是将该单元内的磁化参数置为0,加剧了磁化参数的不连续性。而Tikhonov方法的处理思路:在奇异值上叠加一个极小的正实数α(正则化参数),抑制异常小的奇异值σi对噪声的放大作用,达到正则滤波的目的。其相应的解为xα=k∑i=1σi(uΤic)σ2i+αvi+k∑i=1σi(uΤiξ)σ2i+αvi.(14)从(14)式可以看出,若正则化参数α选取较大(过正则化),则xα与真解差别较大,模型的精确性降低;若正则化参数α选取较小(欠正则化),则无法抑制对噪声的放大作用,模型的稳定性下降。所以恰当的正则化参数需要在模型的精确性和稳定性之间做出一个合适的折中。2.2正则化参数的选取正则化参数的选取方法分为两类:先验选取法和后验选取法。先验选取方法有Morozov偏差原理,该类方法需要知道与测量数据相关的先验信息(如测量数据的误差水平),然而这些先验信息在潜艇磁场反演建模中非常少。后验选取方法主要有L曲线法,交叉校验(CV)和广义交叉校验(GCV).这3种方法直接从测量数据出发,通过预解方程组来确定正则化参数,而不需要先验信息,适合本文中逆问题的求解。1)L曲线法如果以‖xa‖为横坐标,‖Axa-b‖为纵坐标绘成的曲线为L形状,称为L曲线。L曲线反映了逆问题模型的估计参数的精确性和稳定性随着正则化参数α的变化规律。一般选取L曲线上曲率最大的点(都为L曲线的拐点)所对应的α作为正则化参数,但是磁场延拓模型实际的L曲线在拐点附近集聚了大量的α,最优的α值难以确定。所以本文借助极小化如下辅助函数来选取α.ρ(α)=lg(‖Axα-b‖)+lg(‖xα‖),(15)磁场延拓模型的L曲线和辅助函数ρ(α)曲线见图3和图4.2)CV和GCV定义如下的CV和GCV函数:CVα=1nn∑i=1bi-ˆbi1-Ηii(α)2,(16)GCVα=1nn∑i=1bi-ˆbi1-tr(Η(α))/n2,(17)式中:Hα=A(ATA+αI)-1AT;■=Hαb;tr(H(α))为H(α)的迹;Hii(α)是矩阵H(α)的对角元。使函数CV或GCV取得最小值的α就是由交叉校验或广义交叉校验所求得的正则化参数。潜艇磁场延拓模型的CV函数和GCV函数曲线见图5和图6.从图4～6可以看出,CV函数确定的α最小,其次是GCV函数,最大的是L曲线。大量分析发现:CV函数确定的模型精度较高,但鲁棒性较差;L曲线确定的模型鲁棒性较强,但精度较差;综合比较之下,GCV函数选取的正则化参数比较恰当。本文在船模实验中,采用GCV函数确定正则化参数,取得了较为满意的效果。3船模右病z1-12模以一个长308cm、宽34cm的铁质结构船模为实验对象,采用三分量磁通门传感器阵列进行测量。在4m×0.4m的测量范围内,以等间距分别测量得到21×5点阵列上的磁感应强度,如图7所示(船模中心在测量平面上的投影与平面中心重合,测量平面高度分别为z1=57.8cm和z2=72cm).闭合边界面取包含船模的长方体表面,测量平面较其在水平面的投影稍大,整个边界共划分为280个单元。分别利用z1(z2)平面上的磁感应强度三分量的测量数据,建立磁场模型,并由模型计算得到了z2(z1)平面磁感应强度的三分量值。定义相对残差:‖Bm-Bc‖2/‖Bm,max‖2,推算误差:Bci/Bm,max×100%,其中Bm和Bc分别为磁感应强度的测量值向量和推算值向量,Bm,max为磁感应强度测量值的最大值,Bci为第i点磁感应强度的推算值。从表1可以看出:X分量的相对残差最小,Z分量其次,Y分量最大,这是因为在Y方向的磁场信息比X方向和Z方向少造成的。图8和图9分别给出了近场到远场及远场到近场的推算值与测量值的对比曲线(由于篇幅的原因,只画出潜艇龙骨下的磁场比较图,其他位置有同样的性质)。从图中可以看出,由近场推算远场时,3个分量的推算误差均可控制在6%以内,而从远场推算近场模型误差稍大,但均可控制在9%以内。图10和图11分别给出了经正则化处理和未经正则化处理2种情况下,从近场到远场和从远场到近场,船模右舷z分量的推算值与测量值的对比曲线(其他位置和其他分量具有同样的性质)。从图中可以看出,采用GCV函数正则