姚振寧，劉大明，周?chē)?guó)華，喻洲

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

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

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

1 基于混合等效源的艦船磁場(chǎng)建模方法

1.1 邊界等效源

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

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

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

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

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

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

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

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

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

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

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

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

1.2 磁體等效源

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

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

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

1.3 基于混合等效源的逆問(wèn)題模型

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

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

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

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

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

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

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

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

2.1 正則化原理

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

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

Tikhonov正則化解可表示為

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

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

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

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

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

3 船模實(shí)驗(yàn)

為了檢驗(yàn)本文所述艦船磁場(chǎng)建模方法的有效性，設(shè)計(jì)了一個(gè)船模磁場(chǎng)測(cè)量實(shí)驗(yàn)。

圖4 實(shí)驗(yàn)示意圖Fig.4 The sketch map of experiment

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

式中:Emax為各點(diǎn)絕對(duì)誤差的最大值，Bm，max為各點(diǎn)理論測(cè)量值的最大值，Bc為計(jì)算值向量，Bm為理論測(cè)量值向量。

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

圖6分別給出了船模磁場(chǎng)在左舷、龍骨和右舷下方測(cè)量點(diǎn)的計(jì)算值與2種測(cè)量值的對(duì)比曲線?？梢钥闯?由于船模左舷中間區(qū)域有磁性目標(biāo)的存在，左舷下方測(cè)量點(diǎn)的船模磁場(chǎng)干擾測(cè)量值與理論測(cè)量值的差別很大，其x、y、z分量的最大相對(duì)誤差分別為15.38%、13.67%、12.78%，相對(duì)殘差分別為7.81%、9.35%、7.76%;而采用本文建模方法進(jìn)行計(jì)算，左舷下方測(cè)量點(diǎn)的船模磁場(chǎng)計(jì)算值與理論測(cè)量值吻合較好，其x、y、z分量的最大相對(duì)誤差分別為7.23%、4.07%、2.69%，相對(duì)殘差分別為5.76%、4.04%、2.17%;另外，測(cè)量點(diǎn)越遠(yuǎn)離磁性目標(biāo)，干擾測(cè)量值與理論測(cè)量值差別越不明顯，因而右舷下方測(cè)量點(diǎn)差別最不明顯，其計(jì)算值和兩種測(cè)量值3種曲線幾乎重合。

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

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

4 結(jié)論

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

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

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

［7］晁立東，仵杰，王仲奕.工程電磁場(chǎng)基礎(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ì)的時(shí)域逆散射方法［J］.物理學(xué)報(bào)，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.