單志超，周家新，劉賢忠，楊 磊，趙海英

（1.海軍航空大學(xué)，山東煙臺(tái)264001；2.92351部隊(duì)，海南三亞572000；3.92474部隊(duì)，海南三亞572018）

鐵磁性目標(biāo)在地球磁場(chǎng)的磁化下產(chǎn)生的磁場(chǎng)與其周圍的地磁場(chǎng)疊加產(chǎn)生磁異常[1-2]，航空磁探儀通過(guò)探測(cè)磁異常信號(hào)來(lái)發(fā)現(xiàn)目標(biāo)[3]。為提高航空磁探儀的檢測(cè)能力[4]，進(jìn)一步對(duì)鐵磁性目標(biāo)進(jìn)行定位和識(shí)別[5]，需要建立鐵磁性目標(biāo)磁場(chǎng)預(yù)測(cè)模型生成目標(biāo)的高空磁場(chǎng)分布為航空磁探儀提供目標(biāo)信號(hào)參考[6]。高精度的目標(biāo)磁場(chǎng)預(yù)測(cè)模型不但能夠用于輔助航空磁探[7-9]，還可以應(yīng)用于磁性導(dǎo)航[10]、艦船消磁[11]等領(lǐng)域。

目前，建立水下鐵磁性目標(biāo)磁場(chǎng)預(yù)測(cè)模型的方法主要有邊界元法[12]、積分方程法[13]、有限元法[14]和磁體模擬法[15-18]。邊界積分法和曲面積分法對(duì)鐵磁性目標(biāo)磁場(chǎng)建模的原理與邊界元法相同[12]。磁體模擬法主要有磁偶極子模擬[15]、旋轉(zhuǎn)橢球體模擬[19]以及面磁荷模擬等方法。有限元法通過(guò)計(jì)算有限區(qū)域內(nèi)的標(biāo)量磁位來(lái)求解空間場(chǎng)點(diǎn)的總磁場(chǎng)[14]。在近場(chǎng)范圍內(nèi)，有限元需要剖分的場(chǎng)域小，能夠進(jìn)行精細(xì)剖分，計(jì)算精確。當(dāng)計(jì)算場(chǎng)域較大時(shí)，有限元法需要進(jìn)行大量的剖分，實(shí)時(shí)性和計(jì)算精度下降。積分方程法將鐵磁目標(biāo)源區(qū)進(jìn)行離散化，分解為一個(gè)個(gè)鐵磁性目標(biāo)源，將空間場(chǎng)點(diǎn)處的磁場(chǎng)值認(rèn)為是源區(qū)所有目標(biāo)源的疊加。三維積分方程法考慮目標(biāo)源的三維特性，減少假設(shè)條件，使得計(jì)算精度較傳統(tǒng)的二維積分方程法更高。

1 三維積分方程法的基本原理

空間磁場(chǎng)預(yù)測(cè)點(diǎn)P(x,y,z)處的磁場(chǎng)可以認(rèn)為是源區(qū)的鐵磁體以及空間自由電流的矢量疊加，如圖1所示，則點(diǎn)P處的磁場(chǎng)為：

式（1）中：HJ(P)為源區(qū)中空間電流在點(diǎn)P處產(chǎn)生的磁場(chǎng)；HM(P)為源區(qū)中鐵磁性物體產(chǎn)生的磁場(chǎng)。

通過(guò)求解空間中的標(biāo)量磁位分布φ(P)，可以解得源區(qū)中鐵磁性物體產(chǎn)生的磁場(chǎng)：

式（2）、（3）中：grad為梯度算子；BM(P)為源區(qū)中鐵磁性物體產(chǎn)生磁場(chǎng)的磁感應(yīng)強(qiáng)度；μ0=4π×10-7H/m為真空磁導(dǎo)率。

建立如圖2所示的柱坐標(biāo)，由矢量格林公式可以得到，空間場(chǎng)點(diǎn)P處的標(biāo)量磁位為：

式（4）中：div為散度算子；M(N)為空間場(chǎng)點(diǎn)N處的磁化強(qiáng)度；rNP為空間點(diǎn)P到場(chǎng)點(diǎn)N的矢徑；Vm為場(chǎng)點(diǎn)N處磁源的體積；n是邊界Sm的外法相矢量。

在柱坐標(biāo)系中，可以得到：

當(dāng)Mr和Mz為常數(shù)時(shí)，可以得到：

式（7）中，curl為旋度算子。

在航空磁異常探測(cè)中，鐵磁性目標(biāo)的高磁場(chǎng)無(wú)法直接使用式（4）、（5）得到，一般認(rèn)為源區(qū)的空間自由電流J=0，則可以得到點(diǎn)P處的磁感應(yīng)強(qiáng)度為：

在柱坐標(biāo)系中，磁場(chǎng)強(qiáng)度可分解為Hr和Hz兩個(gè)分量。航空磁探中鐵磁性目標(biāo)的磁體材料一般可以認(rèn)為是均勻，將磁源剖分成N1個(gè)單元，并假設(shè)第n個(gè)單元處磁化強(qiáng)度滿足

式中，n=1,2,…,N1為剖分單元的編號(hào)。

則可以得到：

式（10）、（11）中：|H|表示磁場(chǎng)強(qiáng)度的絕對(duì)值；χ(|H|)表示材料的磁化率；是由第n個(gè)磁源位置以及空間觀測(cè)點(diǎn)P決定的幾何系數(shù)，其物理含義是第n個(gè)磁源的單位磁化強(qiáng)度在空間觀測(cè)點(diǎn)P處產(chǎn)生的磁場(chǎng)強(qiáng)度。

式（12）～（16）中：rN1、rN2、zN1、zN2分別為柱坐標(biāo)系中第n個(gè)單元的邊界，rN1＜rN2＜zN1＜zN2。

2 水下目標(biāo)磁場(chǎng)預(yù)測(cè)模型

根據(jù)三維積分方程法的基本原理，求解N1個(gè)離散單元內(nèi)部的磁感應(yīng)強(qiáng)度。

式中：μri為第i個(gè)單元在地磁場(chǎng)磁化下的磁導(dǎo)率；vi為第i個(gè)單元的中心；Qi為第i個(gè)單元的中心。

由式（17）即可實(shí)現(xiàn)水下目標(biāo)在空間任意場(chǎng)點(diǎn)P處磁場(chǎng)的預(yù)測(cè)。直接計(jì)算式（17）中耦合系數(shù)體積分形式相當(dāng)困難，不妨先將其轉(zhuǎn)化為等價(jià)面積分形式。對(duì)線性均勻磁介質(zhì)或均勻磁化單元，假設(shè)源區(qū)內(nèi)的自由電流為0，由矢量恒等式和高斯散度定理可以得到：

式（18）中：si為第i個(gè)單元的表面積；ni為第i個(gè)單元的外法線方向單位矢量。

對(duì)式（18）進(jìn)行負(fù)梯度運(yùn)算可得：

進(jìn)而有：

將面積分中矢量展開(kāi)成分量形式，可得式（20）的矩陣形式：

由式（21）可以看出，該單元耦合系數(shù)各元素為面積分形式，較體積分形式更方便求解。由于計(jì)算場(chǎng)點(diǎn)取在單元中心，而面積分在單元表面進(jìn)行，即使當(dāng)i=j時(shí)，面積分耦合系數(shù)的計(jì)算也不存在積分奇異性問(wèn)題，耦合系數(shù)具體求解過(guò)程可以參見(jiàn)文獻(xiàn)[17]。

3 模型精度理論驗(yàn)證

為檢驗(yàn)本文的水下目標(biāo)磁場(chǎng)預(yù)測(cè)模型精度，使用鐵磁性長(zhǎng)方體對(duì)模型精度進(jìn)行驗(yàn)證。圖3所示為鐵磁性長(zhǎng)方體，磁化率為99，長(zhǎng)1 m，寬0.5 m，高0.3 m，放置在地磁場(chǎng)垂向分量為34 500 nT的外部磁場(chǎng)中。鐵磁性長(zhǎng)方體的空間磁場(chǎng)分布計(jì)算可以參見(jiàn)文獻(xiàn)[17]。

以鐵磁性長(zhǎng)方體中心為原點(diǎn)，長(zhǎng)方體的長(zhǎng)軸指向磁東向，建立圖4坐標(biāo)系。其中：X軸平行于水平面，以指向地磁東向?yàn)檎?，稱縱軸；Y軸平行于水平面，以指向地磁南向?yàn)檎?，稱橫軸；Z軸垂直于水平面，以向下為正，稱垂軸。Bx、By、Bz分別是水下目標(biāo)空間磁場(chǎng)BP在X軸、Y軸、Z軸的投影，分別稱縱向、橫向、垂向分量。空間磁場(chǎng)預(yù)測(cè)點(diǎn)P取在長(zhǎng)方體中心正上方，x∈[-10m,10m]，y=0 m，z=-3 m?？臻g磁場(chǎng)預(yù)測(cè)值與理論值的對(duì)比如圖5所示。

由圖5的鐵磁性長(zhǎng)方體空間磁場(chǎng)理論值與磁場(chǎng)預(yù)測(cè)模型計(jì)算得到的預(yù)測(cè)值的對(duì)比，可以發(fā)現(xiàn)，本文建立的水下目標(biāo)磁場(chǎng)預(yù)測(cè)模型在磁場(chǎng)的變化趨勢(shì)、幅度上擬合較準(zhǔn)確。

4 實(shí)驗(yàn)驗(yàn)證

實(shí)驗(yàn)使用長(zhǎng)半軸為1 m短半軸為0.35 m的鐵磁性長(zhǎng)旋轉(zhuǎn)橢球體作為水下目標(biāo)，以橢球體的中心作為坐標(biāo)原點(diǎn)，測(cè)量橢球體正上方x∈[-10 m,10 m]，y=0 m，z=-5 m處的磁感應(yīng)強(qiáng)度值。實(shí)驗(yàn)中，磁場(chǎng)測(cè)量使用G858銫光泵磁力儀，分辨率可達(dá)0.05 nT。根據(jù)磁場(chǎng)預(yù)測(cè)模型，計(jì)算空間場(chǎng)點(diǎn)的磁場(chǎng)三分量預(yù)測(cè)值，將三分量預(yù)測(cè)值投影到地磁場(chǎng)方向與光泵測(cè)量值進(jìn)行比較，結(jié)果如圖6所示。

圖6為磁感應(yīng)強(qiáng)度預(yù)測(cè)值BP和實(shí)際測(cè)量值BM，平均絕對(duì)誤差A(yù)T和平均相對(duì)誤差A(yù)R定義如下：

可計(jì)算得：AT=0.320 5 nT，AR=12.368%。

5 結(jié)論

本文使用三維積分方程法對(duì)水下目標(biāo)的空間磁場(chǎng)進(jìn)行建模，通過(guò)鐵磁性長(zhǎng)方體對(duì)模型進(jìn)行理論驗(yàn)證，使用銫光泵磁力儀測(cè)量鐵磁性長(zhǎng)旋轉(zhuǎn)橢球體的高空磁場(chǎng)；然后，用本文的預(yù)測(cè)模型推算鐵磁性長(zhǎng)旋轉(zhuǎn)橢球體的磁場(chǎng)，根據(jù)實(shí)際測(cè)量的磁場(chǎng)數(shù)據(jù)和預(yù)測(cè)值進(jìn)行比較。結(jié)果表明，預(yù)測(cè)模型的推算精度較高，平均絕對(duì)誤差為0.320 5 nT，平均相對(duì)誤差為12.368%。

[1]林春生，龔沈光.艦船物理場(chǎng)[M].北京：兵器工業(yè)出版社，2007：45-49.LIN CHUNSHENG，GONG SHENGUANG.Physical field of ship[M].Beijing：Weapon Industry Press，2007：45-49.（in Chinese）

[2]崔漢國(guó).艦船磁場(chǎng)可視化技術(shù)的研究與實(shí)現(xiàn)[J].計(jì)算機(jī)仿真，2004，21（4）：22-24.CUI HANGUO.Research and application of visualization technology for vessel’s magnetic field[J].Computer Simulation，2004，21（4）：22-24.（in Chinese）

[3]SHEINKER A，F(xiàn)RUMKIS L，GINZBURG B.Magnetic anomaly detection using a three-axis magnetometer[J].IEEE Transactions on Magnetics，2009，45（1）：160-167.

[4]熊雄，楊日杰，王鴻吉.海浪磁噪聲背景中動(dòng)目標(biāo)航空磁異常檢測(cè)算法[J].華中科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2015，43（5）：100-105.XIONG XIONG，YANG RIJIE，WANG HONGJI.Airborne magnetic anomaly detection algorithm for moving target under ocean wave generated magnetic noise[J].Journal of Huazhong University of Science and Technology：Natural Science Edition，2015，43（5）：100-105.（in Chinese）

[5]林春生.艦船磁場(chǎng)信號(hào)檢測(cè)與磁性目標(biāo)定位[D].武漢：海軍工程大學(xué)，1996.LIN CHUNSHENG.Signal detection of ship magnetic field and magnetic target location[D].Wuhan：Naval Engineering Institute，1996.（in Chinese）

[6]翁行泰，曹梅芬.磁異探潛中潛艇的數(shù)學(xué)模型[J].上海交通大學(xué)學(xué)報(bào)，1995，29（3）：27-32.WENG XINGTAI，CAO MEIFEN.Submarine math model of magnetic abnormity detection[J].Journal of Shanghai Jiaotong University，1995，29（3）：27-32.（in Chinese）

[7]DAMES P M，SCHWAGER M，SCHWAGER D.Active magnetic anomaly detection using multiple micro aerial vehicles[J].IEEE Robotics and Automation Letters，2016，1（1）：153-160.

[8]SHEINKER A，MOLDWIN M B.Magnetic anomaly detection（MAD）of ferromagnetic pipelines using principal component analysis（PCA）[J].Measurement Science and Technology，2016，27（4）：1-7.

[9]SONG L，BILLINGS S，PASION L，et al.Transient electromagnetic scattering of a metallic object buried in underwater sediments[J].IEEE Transactions on Geoscience and Remote Sensing，2016，54（2）：1091-1102.

[10]LIU Y，ZHANG Y，YI H.The new magnetic survey method for underwater pipeline detection[J].Applied Mechanics and Materials，2013，239（2）：338-343.

[11]高俊吉，劉大明，姚瓊薈.艦船固定磁性分解方法研究[J].哈爾濱工程大學(xué)學(xué)報(bào)，2007，28（10）：1164-1170.GAO JUNJI，LIU DAMING，YAO QIONGHUI.Decomposition of the permanent magnetic field of ships[J].Journal of Harbin Engineering University，2007，28（10）：1164-1170.（in Chinese）

[12]楊明明，劉大明，劉勝道.采用邊界積分方程和Tikhonov正則化方法延拓潛艇磁場(chǎng)[J].兵工學(xué)報(bào)，2010，31（9）：1216-1221.YANG MINGMING，LIU DAMING，LIU SHENGDAO.Submarine’s magnetic field extrapolation using boundary element method and tikhonov regularization[J].Acta Armamentarii，2010，31（9）：1216-1221.（in Chinese）

[13]郭成豹，肖昌漢，劉大明.基于積分方程法和奇異值分解的磁性目標(biāo)磁場(chǎng)延拓技術(shù)研究[J].物理學(xué)報(bào)，2008，57（7）：4182-4188.GUO CHENGBAO，XIAO CHANGHAN，LIU DAMING.Research on the continuations of magnetic field of magnetic object based on integral equation method and singular value decomposition[J].Acta Physica Sinica，2008，57（7）：4182-4188.（in Chinese）

[14]翁行泰，曹梅芬.潛艇感應(yīng)磁場(chǎng)的三維有限元法計(jì)算研究[J].上海交通大學(xué)學(xué)報(bào)，1994，28（5）：69-76.WENG XINGTAI，CAO MEIFEN.A research on the calculation of the induced magnetic field of a submarine[J].Journal of Shanghai Jiaotong University，1994，28（5）：69-76.（in Chinese）

[15]牛龍飛，李斌，張靜.三磁偶極子磁體艦船磁場(chǎng)仿真[J].計(jì)算機(jī)仿真，2009，26（3）：12-15.NIU LONGFEI，LI BIN，ZHANG JING.Simulation of warship magnetic field based on three magnetic dipole magnetic bodies[J].Computer Simulation，2009，26（3）：12-15.（in Chinese）

[16]周家新，陳建勇，單志超，等.航空磁探中潛艇磁場(chǎng)建模方法分析[J].海軍航空工程學(xué)院學(xué)報(bào)，2017，32（1）：143-148.ZHOU JIAXIN，CHEN JIANYONG，SHAN ZHICHAO，et al.Analysis of submarine magnetic field modeling method for aerial magnetic detection[J].Journal of Naval Aeronautical and Astronautical University，2017，32（1）：143-148.（in Chinese）

[17]周耀忠，張國(guó)友.艦船磁場(chǎng)分析計(jì)算[M].北京：國(guó)防工業(yè)出版社，2004：104-200.ZHOU YAOZHONG，ZHANG GUOYOU.Ship magnetic field analyze and compute[M].Beijing：National Defend Industry Press，2004：104-200.（in Chinese）

[18]OLIVIER C，COULOMB J，BONGIRAUD J.Recent improvements for solving inverse magnetostatic problem applied to thin shells[J].IEEE Transactions on Magnetics，2002，38（2）：1005-1008.

[19]劉勝道，劉大明，肖昌漢，等.基于遺傳算法的磁性目標(biāo)磁模型[J].武漢理工大學(xué)學(xué)報(bào)，2008，32（6）：1017-1020.LIU SHENGDAO，LIU DAMING，XIAO CHANGHAN，et al.Magnetic model of the magnetic object based on genetic algorithm[J].Journal of Wuhan University of Technology，2008，32（6）：1017-1020.（in Chinese）