磁性目標(biāo)干擾下的艦船磁場建模

姚振寧，劉大明，周國華，喻洲

(海軍工程大學(xué)電氣工程學(xué)院，湖北武漢430033)

由于艦船形狀不規(guī)則、磁化不均勻，用解析方法很難計算其磁場。目前，用于艦船磁場計算的方法有磁體模擬法、邊界元法和積分方程法等［1］。文獻［2-4］利用磁場測量數(shù)據(jù)，通過求解正演問題中的線性方程組反演得到模型參數(shù)，即建立了艦船磁場模型。文獻［5］利用某磁化狀態(tài)下鐵磁物體的磁場測量數(shù)據(jù)，反演出鐵磁物體的等效磁化率而得到感應(yīng)磁場計算模型，這種方法也可以用于艦船感應(yīng)磁場的計算。但上述幾種方法只能對單獨存在的艦船進行建模，對周圍還有磁性目標(biāo)存在的艦船難以建立數(shù)學(xué)模型，從而無法由原始磁場測量數(shù)據(jù)計算艦船磁場。

通過把邊界元法和磁體模擬法相結(jié)合，本文提出了一種在磁性目標(biāo)干擾下建立艦船磁場模型的方法，首先用邊界等效源和磁體等效源組成的混合等效源等效艦船和磁性目標(biāo)，然后通過Tikhonov正則化方法反演求解出混合等效源參數(shù)，進而由獲得的邊界等效源參數(shù)計算出艦船磁場。

1 基于混合等效源的艦船磁場建模方法

1.1 邊界等效源

如圖1所示，S為包圍磁源的閉合邊界面，S∞為延伸到無窮遠處的無限大閉合邊界面，S之外、S∞之內(nèi)的空間V為所研究場域。若場域V內(nèi)無其他的磁源(鐵磁物質(zhì)或電流)存在，則標(biāo)量磁位φ在V內(nèi)滿足拉普拉斯方程，記φ*=，由格林第二公式可得［1，6］

式中:φM為場域V內(nèi)任意場點M的磁位，r為邊界面上的點到場點M的距離，n為方向指向場域外的邊界面法線，φ和φ＊為邊界面上的磁位和磁位的法向?qū)?shù)。cM取值如下:若場點M位于場域V內(nèi)，則cM=1;若場點M位于光滑邊界面上，則cM=1/2;若場點M位于場域V外，則cM=0。

圖1 磁源與邊界面示意圖Fig.1 The sketch map of magnetic source and boundary surface

當(dāng)r→∞時，式(1)中的被積函數(shù)將隨1/r3減小，在S∞面上的積分為零，則式(1)化為

從式(2)看出，磁源在區(qū)域V內(nèi)任意點產(chǎn)生的效應(yīng)可用邊界面S上φ和φ＊的面積分等效表示。因此，把φ和φ＊作為邊界等效源參數(shù)。

從廣義敘述學(xué)的視角來看，無論是新聞還是電影，都是一種由符號組合構(gòu)成的敘述文本。如果從一般意義的敘述分類來說，以事實報道為主的新聞體裁屬于紀(jì)實型敘述，以反映劇情為主的電影體裁(非紀(jì)錄片)屬于虛構(gòu)型敘述。

由式(2)可以計算場域V內(nèi)任意場點M的磁位φM(cM=1)，然后求其負梯度，即可得到該點處的磁感應(yīng)強度:

式中:μ為該點處磁導(dǎo)率。

由式(2)可以導(dǎo)出邊界等效源參數(shù)φ和φ＊滿足的關(guān)系(cM=1/2):

另外，邊界等效源參數(shù)φ＊還應(yīng)該滿足磁通連續(xù)性原理［7］，即

如果在場域V內(nèi)通過測量得到一組磁感應(yīng)強度數(shù)據(jù)，則根據(jù)式(3)～(5)能夠反演出邊界等效源參數(shù)φ和φ＊，進而根據(jù)式(3)可計算場域V內(nèi)其他場點的磁場。

采用解析方法求解以上積分方程是十分困難的，但可以將邊界面S離散成一系列小的面積元，也即為邊界元。邊界元采用常數(shù)單元，即把每個邊界元上的φ和φ＊值都設(shè)定為相應(yīng)的常數(shù)，且等于該邊界元中心點上的值。設(shè)將邊界面S離散成N個邊界元，用φj和φj＊分別表示第j個邊界元上的φ和φ＊的值，則式(3)～(5)化為式中:φi為第i個邊界元上的φ值，rMj為場點M與第j個邊界元的距離，rij為第i個邊界元與第j個邊界元的距離，ASj為第j個邊界元的面積。式中積分項與邊界幾何形狀有關(guān)，一般由數(shù)值積分算出。

1.2 磁體等效源

在磁體模擬法中，一般用若干個具有特定磁矩的磁性模擬體(磁偶極子或橢球體)所產(chǎn)生的磁場去模擬不規(guī)則物體磁場。因此，把這里的磁性模擬體作為磁體等效源，其磁矩也即為磁體等效源參數(shù)。如圖2所示，不規(guī)則物體磁場采用n個橢球體組成的橢球陣列產(chǎn)生的磁場來模擬。大橢球體的幾何中心與不規(guī)則物體的中心重合，其長軸、短軸與不規(guī)則物體長度、最大寬度相當(dāng);其余n-1個小橢球體均勻分布在大橢球體的x軸線上。設(shè)不規(guī)則物體在任意場點M產(chǎn)生的磁感應(yīng)強度3個分量為BxM、ByM、BzM，則有方程組

式中:mxi、myi、mzi為第 i個橢球的 x、y、z方向磁矩分量;fxMi、fyMi、fzMi、gxMi、gyMi、gzMi、exMi、eyMi、ezMi是與場點M的坐標(biāo)、橢球位置和橢球焦距有關(guān)的常數(shù)。如果在場域內(nèi)通過測量得到一組磁感應(yīng)強度數(shù)據(jù)，則根據(jù)上式能夠反演出磁體等效源參數(shù)(橢球磁矩)，進而可計算場域其他場點的磁場。

圖2 橢球陣列模型Fig.2 The ellipsoidal sphere array model

1.3 基于混合等效源的逆問題模型

當(dāng)艦船周圍有磁性目標(biāo)存在時，磁傳感器的原始磁場測量數(shù)據(jù)是艦船磁場和磁性目標(biāo)磁場之和，如果用單一的等效源去等效艦船與磁性目標(biāo)，則相互耦合比較嚴(yán)重，此時需要用邊界等效源和磁體等效源組成的混合等效源進行等效。

艦船形狀極不規(guī)則，其磁場用參數(shù)較多的邊界等效源模擬，艦船的邊界面原則上可以是任意形狀的閉合面，但為了計算簡單，一般選擇比較規(guī)則的長方體表面。與艦船相比，磁性目標(biāo)一般較小，其磁場用參數(shù)較少的磁體等效源模擬，磁性模擬體選擇橢球陣列。圖3是艦船的邊界面和磁性目標(biāo)的橢球陣列示意圖。

圖3 艦船的邊界面和磁性目標(biāo)的橢球陣列示意圖Fig.3 The sketch map of boundary surface and ellipsoidal sphere array

由于獲得的原始磁場測量數(shù)據(jù)是艦船磁場和磁性目標(biāo)磁場之和，因而把式(6)與式(9)相加合并為一個等式并與式(7)、(8)統(tǒng)一寫成矩陣形式

式中:φ=［φ1φ2…φN］T和 φ*=［…φN*］T為待求的兩類邊界等效源參數(shù)向量，m=［mx1mx2…mxnmy1my2… mynmz1mz2… mzn］T為待求的磁體等效源參數(shù)向量;BM為已知的原始磁場測量數(shù)據(jù)向量，通過磁傳感器測量獲得;FM和GM是將原始磁場測量數(shù)據(jù)和邊界等效源聯(lián)系起來的觀測矩陣，F(xiàn)和G是描述兩類邊界等效源之間的關(guān)系矩陣，As為邊界元面積向量，a是將原始磁場測量數(shù)據(jù)和磁體等效源聯(lián)系起來的系數(shù)矩陣，以上矩陣均可從式(6)～(9)得到，將式(10)進一步合并成方程組

式中:x=［φφ*m］T為混合等效源參數(shù)。由式(11)反演求解出混合等效源參數(shù)，即建立了基于混合等效源的逆問題模型，再根據(jù)式(6)由獲得的邊界等效源參數(shù)計算出待求的艦船磁場。

2 Tikhonov正則化求解混合等效源參數(shù)

在邊界元法中，為了保證計算精度，一般將邊界劃分較多的邊界單元，邊界元的增多，會使得邊界元方程組系數(shù)矩陣的列相關(guān)性增強，系數(shù)矩陣的條件數(shù)增大，導(dǎo)致方程組呈現(xiàn)嚴(yán)重的病態(tài)特性。另外，艦船與磁性目標(biāo)的距離越近，使邊界等效源與磁體等效源的模型參數(shù)相關(guān)性越強，即使求出方程組的最小二乘解，也不能真實反映艦船與磁性目標(biāo)的磁性。要使病態(tài)方程有穩(wěn)定的解，必須采用正則化方法，Tikhonov正則化是處理方程組病態(tài)問題的常用正則化方法［8-9］。

2.1 正則化原理

將式(11)的系數(shù)矩陣進行奇異值分解，即A=USVT，則其最小二乘解可表示為

式中:ui和 vi分別為正交矩陣 U和 V的第 i列，σ1≥σ2≥…≥σk＞0為A的奇異值，k為A的秩?？梢钥闯觯?dāng)奇異值σi很小時，b中很小的噪聲誤差都將在解x中產(chǎn)生很大的誤差分量。

Tikhonov正則化解可表示為

式中:λ為正則化參數(shù)，且滿足σ1≥λ≥σk?？梢钥闯?，Tikhonov正則化解其實是在最小二乘解基礎(chǔ)上加入了一個濾波因子fi=/(σi2+λ2)，濾波因子隨著奇異值減小而減小，從而抑制極小奇異值對噪聲的放大作用。

2.2 正則化參數(shù)的選取

如果正則化參數(shù)λ選擇合理，正則化解將很穩(wěn)定，如果正則化參數(shù)λ選擇過大或過小，會造成有用信號的過濾或無法抑制噪聲。

L曲線法可以很好地選取正則化參數(shù)［10-12］，即以正則化參數(shù)λ為參變量繪出正則化解范數(shù)‖xλ‖對殘差范數(shù)‖Axλ-b‖變化的L形狀曲線。L曲線上彎曲最厲害(曲率最大)的點所對應(yīng)的λ值作為正則化參數(shù)，但是一般實際的L曲線拐角處聚集了大量的λ，正則化參數(shù)仍然難以確定。因此，可以尋找輔助函數(shù) lg(‖Axλ-b‖)+lg(‖xλ‖)的最小值，其所對應(yīng)的λ值作為正則化參數(shù)。

另外，磁性目標(biāo)的磁場一般較小，而原始磁場測量數(shù)據(jù)中一般存在較大的人為測量誤差、傳感器固有誤差，這些誤差將嚴(yán)重影響磁體等效源對磁性目標(biāo)的建模，進而影響邊界等效源對艦船的建模，此時直接利用L曲線法選取正則化參數(shù)將造成過正則化(選取的正則化參數(shù)大于最優(yōu)正則化參數(shù))。因此，要對原始磁場測量數(shù)據(jù)進行一次預(yù)處理，即先通過磁體模擬法求出原始磁場測量數(shù)據(jù)的擬合數(shù)據(jù)，然后通過L曲線法進行正則化參數(shù)的選取，由得到的正則化參數(shù)再對原方程組進行正則化處理。

3 船模實驗

為了檢驗本文所述艦船磁場建模方法的有效性，設(shè)計了一個船模磁場測量實驗。

圖4 實驗示意圖Fig.4 The sketch map of experiment

如圖4所示，鐵質(zhì)船模長4.8 m、最大寬0.54 m，磁性目標(biāo)是長0.8 m、寬0.11 m的長方體形狀的鐵磁物體，平行放置在船模左舷外側(cè)，用5×49三分量磁通門傳感器陣列在船模下方平面進行磁場測量。當(dāng)磁性目標(biāo)存在時，測量船模下方磁場，其測量值為干擾測量值(即建模用的原始磁場測量數(shù)據(jù));當(dāng)磁性目標(biāo)不存在時，測量船模下方磁場，把其測量值當(dāng)作船模磁場的理論測量值。船模的閉合邊界面被剖分為690個邊界單元，磁性目標(biāo)用橢球陣列等效。反演等效源參數(shù)首先需要通過L曲線法選取正則化參數(shù)，圖5給出了本實驗的輔助函數(shù)隨正則化參數(shù)變化的曲線，當(dāng)輔助函數(shù)最小時，選取的正則化參數(shù)值為0.003。為了分析數(shù)據(jù)，定義最大相對誤差(REmax)和相對殘差(RRE)分別為

式中:Emax為各點絕對誤差的最大值，Bm，max為各點理論測量值的最大值，Bc為計算值向量，Bm為理論測量值向量。

圖5 輔助函數(shù)曲線Fig.5 The curve of accessorial function

圖6分別給出了船模磁場在左舷、龍骨和右舷下方測量點的計算值與2種測量值的對比曲線?？梢钥闯?由于船模左舷中間區(qū)域有磁性目標(biāo)的存在，左舷下方測量點的船模磁場干擾測量值與理論測量值的差別很大，其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%;另外，測量點越遠離磁性目標(biāo)，干擾測量值與理論測量值差別越不明顯，因而右舷下方測量點差別最不明顯，其計算值和兩種測量值3種曲線幾乎重合。

圖6 輸入和輸出的自功率譜Fig.6 Auto-power spectrum of input and output

由圖6可知，計算值的最大相對誤差和相對殘差均遠小于干擾測量值，所建模型能夠由磁性目標(biāo)干擾下的原始磁場測量數(shù)據(jù)有效計算出船模磁場。

4 結(jié)論

本文提出了在磁性目標(biāo)干擾下建立艦船磁場數(shù)學(xué)模型的方法，該方法用邊界等效源等效艦船以及磁體等效源等效磁性目標(biāo)，根據(jù)原始磁場測量數(shù)據(jù)建立基于混合等效源的逆問題模型，用Tikhonov正則化方法反演求解模型的混合等效源參數(shù)，進而由獲得的邊界等效源參數(shù)計算出艦船磁場。實驗結(jié)果表明，基于混合等效源的艦船磁場建模方法是可行的，能夠由原始磁場測量數(shù)據(jù)有效計算出艦船磁場，且計算結(jié)果誤差較小，能夠滿足實際工程需求。

［1］周耀忠，張國友.艦船磁場分析計算［M］.北京:國防工業(yè)出版社，2004:127-200.

［2］王金根，龔沈光，劉勝道.磁性目標(biāo)的高精度建模方法［J］.海軍工程大學(xué)學(xué)報，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］.兵工學(xué)報，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］.兵工學(xué)報，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］周國華，肖昌漢，閆輝，等.一種弱磁作用下鐵磁物體感應(yīng)磁場的計算方法［J］.哈爾濱工程大學(xué)學(xué)報，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］倪光正，楊仕友，錢秀英，等.工程電磁場數(shù)值計算［M］.北京:機械工業(yè)出版社，2004:238-242.

［7］晁立東，仵杰，王仲奕.工程電磁場基礎(chǔ)［M］.西安:西北工業(yè)大學(xué)出版社，2001:18-19.

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

［9］劉光東，張業(yè)榮.二維有耗色散介質(zhì)的時域逆散射方法［J］.物理學(xué)報，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.