磁性目标干扰下的舰船磁场建模

姚振宁，刘大明，周国华，喻洲

(海军工程大学电气工程学院，湖北武汉430033)

由于舰船形状不规则、磁化不均匀，用解析方法很难计算其磁场。目前，用于舰船磁场计算的方法有磁体模拟法、边界元法和积分方程法等［1］。文献［2-4］利用磁场测量数据，通过求解正演问题中的线性方程组反演得到模型参数，即建立了舰船磁场模型。文献［5］利用某磁化状态下铁磁物体的磁场测量数据，反演出铁磁物体的等效磁化率而得到感应磁场计算模型，这种方法也可以用于舰船感应磁场的计算。但上述几种方法只能对单独存在的舰船进行建模，对周围还有磁性目标存在的舰船难以建立数学模型，从而无法由原始磁场测量数据计算舰船磁场。

通过把边界元法和磁体模拟法相结合，本文提出了一种在磁性目标干扰下建立舰船磁场模型的方法，首先用边界等效源和磁体等效源组成的混合等效源等效舰船和磁性目标，然后通过Tikhonov正则化方法反演求解出混合等效源参数，进而由获得的边界等效源参数计算出舰船磁场。

1 基于混合等效源的舰船磁场建模方法

1.1 边界等效源

如图1所示，S为包围磁源的闭合边界面，S∞为延伸到无穷远处的无限大闭合边界面，S之外、S∞之内的空间V为所研究场域。若场域V内无其他的磁源(铁磁物质或电流)存在，则标量磁位φ在V内满足拉普拉斯方程，记φ*=，由格林第二公式可得［1，6］

式中:φM为场域V内任意场点M的磁位，r为边界面上的点到场点M的距离，n为方向指向场域外的边界面法线，φ和φ＊为边界面上的磁位和磁位的法向导数。cM取值如下:若场点M位于场域V内，则cM=1;若场点M位于光滑边界面上，则cM=1/2;若场点M位于场域V外，则cM=0。

图1 磁源与边界面示意图Fig.1 The sketch map of magnetic source and boundary surface

当r→∞时，式(1)中的被积函数将随1/r3减小，在S∞面上的积分为零，则式(1)化为

从式(2)看出，磁源在区域V内任意点产生的效应可用边界面S上φ和φ＊的面积分等效表示。因此，把φ和φ＊作为边界等效源参数。

从广义叙述学的视角来看，无论是新闻还是电影，都是一种由符号组合构成的叙述文本。如果从一般意义的叙述分类来说，以事实报道为主的新闻体裁属于纪实型叙述，以反映剧情为主的电影体裁(非纪录片)属于虚构型叙述。

由式(2)可以计算场域V内任意场点M的磁位φM(cM=1)，然后求其负梯度，即可得到该点处的磁感应强度:

式中:μ为该点处磁导率。

由式(2)可以导出边界等效源参数φ和φ＊满足的关系(cM=1/2):

另外，边界等效源参数φ＊还应该满足磁通连续性原理［7］，即

如果在场域V内通过测量得到一组磁感应强度数据，则根据式(3)～(5)能够反演出边界等效源参数φ和φ＊，进而根据式(3)可计算场域V内其他场点的磁场。

采用解析方法求解以上积分方程是十分困难的，但可以将边界面S离散成一系列小的面积元，也即为边界元。边界元采用常数单元，即把每个边界元上的φ和φ＊值都设定为相应的常数，且等于该边界元中心点上的值。设将边界面S离散成N个边界元，用φj和φj＊分别表示第j个边界元上的φ和φ＊的值，则式(3)～(5)化为式中:φi为第i个边界元上的φ值，rMj为场点M与第j个边界元的距离，rij为第i个边界元与第j个边界元的距离，ASj为第j个边界元的面积。式中积分项与边界几何形状有关，一般由数值积分算出。

1.2 磁体等效源

在磁体模拟法中，一般用若干个具有特定磁矩的磁性模拟体(磁偶极子或椭球体)所产生的磁场去模拟不规则物体磁场。因此，把这里的磁性模拟体作为磁体等效源，其磁矩也即为磁体等效源参数。如图2所示，不规则物体磁场采用n个椭球体组成的椭球阵列产生的磁场来模拟。大椭球体的几何中心与不规则物体的中心重合，其长轴、短轴与不规则物体长度、最大宽度相当;其余n-1个小椭球体均匀分布在大椭球体的x轴线上。设不规则物体在任意场点M产生的磁感应强度3个分量为BxM、ByM、BzM，则有方程组

式中:mxi、myi、mzi为第 i个椭球的 x、y、z方向磁矩分量;fxMi、fyMi、fzMi、gxMi、gyMi、gzMi、exMi、eyMi、ezMi是与场点M的坐标、椭球位置和椭球焦距有关的常数。如果在场域内通过测量得到一组磁感应强度数据，则根据上式能够反演出磁体等效源参数(椭球磁矩)，进而可计算场域其他场点的磁场。

图2 椭球阵列模型Fig.2 The ellipsoidal sphere array model

1.3 基于混合等效源的逆问题模型

当舰船周围有磁性目标存在时，磁传感器的原始磁场测量数据是舰船磁场和磁性目标磁场之和，如果用单一的等效源去等效舰船与磁性目标，则相互耦合比较严重，此时需要用边界等效源和磁体等效源组成的混合等效源进行等效。

舰船形状极不规则，其磁场用参数较多的边界等效源模拟，舰船的边界面原则上可以是任意形状的闭合面，但为了计算简单，一般选择比较规则的长方体表面。与舰船相比，磁性目标一般较小，其磁场用参数较少的磁体等效源模拟，磁性模拟体选择椭球阵列。图3是舰船的边界面和磁性目标的椭球阵列示意图。

图3 舰船的边界面和磁性目标的椭球阵列示意图Fig.3 The sketch map of boundary surface and ellipsoidal sphere array

由于获得的原始磁场测量数据是舰船磁场和磁性目标磁场之和，因而把式(6)与式(9)相加合并为一个等式并与式(7)、(8)统一写成矩阵形式

式中:φ=［φ1φ2…φN］T和 φ*=［…φN*］T为待求的两类边界等效源参数向量，m=［mx1mx2…mxnmy1my2… mynmz1mz2… mzn］T为待求的磁体等效源参数向量;BM为已知的原始磁场测量数据向量，通过磁传感器测量获得;FM和GM是将原始磁场测量数据和边界等效源联系起来的观测矩阵，F和G是描述两类边界等效源之间的关系矩阵，As为边界元面积向量，a是将原始磁场测量数据和磁体等效源联系起来的系数矩阵，以上矩阵均可从式(6)～(9)得到，将式(10)进一步合并成方程组

式中:x=［φφ*m］T为混合等效源参数。由式(11)反演求解出混合等效源参数，即建立了基于混合等效源的逆问题模型，再根据式(6)由获得的边界等效源参数计算出待求的舰船磁场。

2 Tikhonov正则化求解混合等效源参数

在边界元法中，为了保证计算精度，一般将边界划分较多的边界单元，边界元的增多，会使得边界元方程组系数矩阵的列相关性增强，系数矩阵的条件数增大，导致方程组呈现严重的病态特性。另外，舰船与磁性目标的距离越近，使边界等效源与磁体等效源的模型参数相关性越强，即使求出方程组的最小二乘解，也不能真实反映舰船与磁性目标的磁性。要使病态方程有稳定的解，必须采用正则化方法，Tikhonov正则化是处理方程组病态问题的常用正则化方法［8-9］。

2.1 正则化原理

将式(11)的系数矩阵进行奇异值分解，即A=USVT，则其最小二乘解可表示为

式中:ui和 vi分别为正交矩阵 U和 V的第 i列，σ1≥σ2≥…≥σk＞0为A的奇异值，k为A的秩。可以看出，当奇异值σi很小时，b中很小的噪声误差都将在解x中产生很大的误差分量。

Tikhonov正则化解可表示为

式中:λ为正则化参数，且满足σ1≥λ≥σk。可以看出，Tikhonov正则化解其实是在最小二乘解基础上加入了一个滤波因子fi=/(σi2+λ2)，滤波因子随着奇异值减小而减小，从而抑制极小奇异值对噪声的放大作用。

2.2 正则化参数的选取

如果正则化参数λ选择合理，正则化解将很稳定，如果正则化参数λ选择过大或过小，会造成有用信号的过滤或无法抑制噪声。

L曲线法可以很好地选取正则化参数［10-12］，即以正则化参数λ为参变量绘出正则化解范数‖xλ‖对残差范数‖Axλ-b‖变化的L形状曲线。L曲线上弯曲最厉害(曲率最大)的点所对应的λ值作为正则化参数，但是一般实际的L曲线拐角处聚集了大量的λ，正则化参数仍然难以确定。因此，可以寻找辅助函数 lg(‖Axλ-b‖)+lg(‖xλ‖)的最小值，其所对应的λ值作为正则化参数。

另外，磁性目标的磁场一般较小，而原始磁场测量数据中一般存在较大的人为测量误差、传感器固有误差，这些误差将严重影响磁体等效源对磁性目标的建模，进而影响边界等效源对舰船的建模，此时直接利用L曲线法选取正则化参数将造成过正则化(选取的正则化参数大于最优正则化参数)。因此，要对原始磁场测量数据进行一次预处理，即先通过磁体模拟法求出原始磁场测量数据的拟合数据，然后通过L曲线法进行正则化参数的选取，由得到的正则化参数再对原方程组进行正则化处理。

3 船模实验

为了检验本文所述舰船磁场建模方法的有效性，设计了一个船模磁场测量实验。

图4 实验示意图Fig.4 The sketch map of experiment

如图4所示，铁质船模长4.8 m、最大宽0.54 m，磁性目标是长0.8 m、宽0.11 m的长方体形状的铁磁物体，平行放置在船模左舷外侧，用5×49三分量磁通门传感器阵列在船模下方平面进行磁场测量。当磁性目标存在时，测量船模下方磁场，其测量值为干扰测量值(即建模用的原始磁场测量数据);当磁性目标不存在时，测量船模下方磁场，把其测量值当作船模磁场的理论测量值。船模的闭合边界面被剖分为690个边界单元，磁性目标用椭球阵列等效。反演等效源参数首先需要通过L曲线法选取正则化参数，图5给出了本实验的辅助函数随正则化参数变化的曲线，当辅助函数最小时，选取的正则化参数值为0.003。为了分析数据，定义最大相对误差(REmax)和相对残差(RRE)分别为

式中:Emax为各点绝对误差的最大值，Bm，max为各点理论测量值的最大值，Bc为计算值向量，Bm为理论测量值向量。

图5 辅助函数曲线Fig.5 The curve of accessorial function

图6分别给出了船模磁场在左舷、龙骨和右舷下方测量点的计算值与2种测量值的对比曲线。可以看出:由于船模左舷中间区域有磁性目标的存在，左舷下方测量点的船模磁场干扰测量值与理论测量值的差别很大，其x、y、z分量的最大相对误差分别为15.38%、13.67%、12.78%，相对残差分别为7.81%、9.35%、7.76%;而采用本文建模方法进行计算，左舷下方测量点的船模磁场计算值与理论测量值吻合较好，其x、y、z分量的最大相对误差分别为7.23%、4.07%、2.69%，相对残差分别为5.76%、4.04%、2.17%;另外，测量点越远离磁性目标，干扰测量值与理论测量值差别越不明显，因而右舷下方测量点差别最不明显，其计算值和两种测量值3种曲线几乎重合。

图6 输入和输出的自功率谱Fig.6 Auto-power spectrum of input and output

由图6可知，计算值的最大相对误差和相对残差均远小于干扰测量值，所建模型能够由磁性目标干扰下的原始磁场测量数据有效计算出船模磁场。

4 结论

本文提出了在磁性目标干扰下建立舰船磁场数学模型的方法，该方法用边界等效源等效舰船以及磁体等效源等效磁性目标，根据原始磁场测量数据建立基于混合等效源的逆问题模型，用Tikhonov正则化方法反演求解模型的混合等效源参数，进而由获得的边界等效源参数计算出舰船磁场。实验结果表明，基于混合等效源的舰船磁场建模方法是可行的，能够由原始磁场测量数据有效计算出舰船磁场，且计算结果误差较小，能够满足实际工程需求。

［1］周耀忠，张国友.舰船磁场分析计算［M］.北京:国防工业出版社，2004:127-200.

［2］王金根，龚沈光，刘胜道.磁性目标的高精度建模方法［J］.海军工程大学学报，2001，13(3):49-52.WANG Jingen，GONG Shenguang，LIU Shengdao.High accuracy method for modeling magnetic objects［J］.Journal of Naval University of Engineering，2001，13(3):49-52.

［3］杨明明，刘大明，刘胜道，等.采用边界积分方程和Tikhonov正则化方法延拓潜艇磁场［J］.兵工学报，2010，31(9):1216-1221.YANG Mingming，LIU Daming，LIU Shengdao，et al.Submarine’s magnetic field extrapolation using boundary elementmethod and Tikhonov regularization ［J］.Acta Armamentarii，2010，31(9):1216-1221.

［4］周国华，肖昌汉，刘胜道，等.一种舰艇磁隐身中磁场推算的新方法［J］.兵工学报，2009，30(7):951-957.ZHOU Guohua，XIAO Changhan，LIU Shengdao，et al.A new magnetic field extrapolation method in magnetic silencing of ships［J］.Acta Armamentarii，2009，30(7):951-957.

［5］周国华，肖昌汉，闫辉，等.一种弱磁作用下铁磁物体感应磁场的计算方法［J］.哈尔滨工程大学学报，2009，30(1):91-95.ZHOU Guohua，XIAO Changhan，YAN Hui，et al.A method to calculate the induced magnetic field of ferromagnetic objects in aweakmagnetic field［J］.Journal of Harbin Engineering University，2009，30(1):91-95.

［6］倪光正，杨仕友，钱秀英，等.工程电磁场数值计算［M］.北京:机械工业出版社，2004:238-242.

［7］晁立东，仵杰，王仲奕.工程电磁场基础［M］.西安:西北工业大学出版社，2001:18-19.

［8］WINKLER JR.Tikhonov regularisation in standard form for polynomial basis conversion［J］.Appl Math Modelling，1997，21:651-662.

［9］刘光东，张业荣.二维有耗色散介质的时域逆散射方法［J］.物理学报，2010，59(10):6969-6979.LIU Guangdong，ZHANG Yerong.Time-domain inverse scattering problem for two-dimensional frequency-dispersive lossy media［J］.Acta Physica Sinica，2010，59(10):6969-6979.

［10］HANSEN P C.Analysis of discrete ill-posed problems by means of the L-curve［J］.SIAM Review，1992，34(4):561-580.

［11］MORIGIS，SGALLARIF.A regularizing L-curve Lanczos method for under determined linear systems［J］.Appl Math Comput，2001，121:55-73.

［12］CALVETTID，MORIGIS，REICHEL L，et al.Tikhonov regularization and the L-curve for large discrete ill-posed problems［J］.JComput Appl Math，2000，123:423-446.