一種艦船高精度感應(yīng)磁場(chǎng)快速正演建模方法

何保委 孫兆龍 劉月林 周國(guó)華 唐烈崢

一種艦船高精度感應(yīng)磁場(chǎng)快速正演建模方法

何保委 孫兆龍 劉月林 周國(guó)華 唐烈崢

（海軍工程大學(xué)電氣工程學(xué)院 武漢 430033）

為了解決矢量積分方程法計(jì)算艦船感應(yīng)磁場(chǎng)時(shí)效率較低和計(jì)算量較大的問題，提出一種基于簡(jiǎn)化標(biāo)量磁位的積分方程正演模型，并引入多層自適應(yīng)交叉近似（MLACA）算法，不僅能夠保證磁場(chǎng)計(jì)算精度，還能大幅減少計(jì)算機(jī)的內(nèi)存需求和計(jì)算時(shí)間。數(shù)值仿真表明，使用基于MLACA的標(biāo)量積分方程法能夠快速獲取高精度的鐵磁物體感應(yīng)磁場(chǎng)。針對(duì)艦船鐵磁材料的磁性參數(shù)不易獲取的問題，基于磁場(chǎng)實(shí)測(cè)值和正演耦合模型，以磁場(chǎng)擬合度、磁化率先驗(yàn)分布和光滑約束條件為目標(biāo)函數(shù)建立磁化率反演模型，并采用模擬退火（SA）算法優(yōu)化得到等效磁化率的空間分布。球殼數(shù)值仿真和船模試驗(yàn)結(jié)果表明，磁化率反演優(yōu)化模型有效可行，其中球殼磁場(chǎng)的預(yù)測(cè)誤差為2.2%，船模磁場(chǎng)的預(yù)測(cè)誤差約為4.8%。

艦船磁場(chǎng) 簡(jiǎn)化標(biāo)量磁位 多層自適應(yīng)交叉近似 模擬退火法 等效磁化率

0 引言

大多數(shù)艦船由鐵磁材料制成，在地磁場(chǎng)的磁化作用和其他因素的影響下會(huì)產(chǎn)生艦船磁場(chǎng)，成為被攻擊或監(jiān)視的信號(hào)源，因此，全面、準(zhǔn)確地了解和掌握艦船磁場(chǎng)的空間分布特性是實(shí)施磁防護(hù)措施的重要前提，也是磁隱身技術(shù)的一項(xiàng)重要研究?jī)?nèi)容[1-2]。艦船磁場(chǎng)按照產(chǎn)生機(jī)理可以分為鐵磁性磁場(chǎng)、渦流磁場(chǎng)、雜散磁場(chǎng)和極低頻磁場(chǎng)等。而鐵磁性磁場(chǎng)又包含了感應(yīng)磁場(chǎng)和固定磁場(chǎng)兩類[3]。獲取艦船固定磁場(chǎng)和感應(yīng)磁場(chǎng)的一般做法是利用艦船下方的磁傳感器陣列直接測(cè)量艦船在不同外磁場(chǎng)作用下產(chǎn)生的磁場(chǎng)并進(jìn)行分解計(jì)算[3]，但是該方法要求艦船左右大致對(duì)稱，且在測(cè)量過程中艦船的位置不能出現(xiàn)較大偏移，測(cè)得的也僅為有限的磁場(chǎng)信息。若想得到艦船的空間磁場(chǎng)分布，還需要利用有限的測(cè)量點(diǎn)磁場(chǎng)值再次進(jìn)行延拓和換算。針對(duì)測(cè)量法這一弊端和通用性不夠的問題，需要研究數(shù)值建模方法來計(jì)算得到更加精確的艦船全空間磁場(chǎng)。艦船感應(yīng)磁場(chǎng)與地磁場(chǎng)的大小、方向以及材料磁導(dǎo)率等因素有關(guān)，對(duì)于非規(guī)則的鐵磁物體來說，可以使用等效源法[3-6]、邊界元法[7]、有限元法[8]和積分方程法[9]等數(shù)值方法進(jìn)行建模。因此，本文主要研究艦船感應(yīng)磁場(chǎng)的建模和計(jì)算方法，進(jìn)而為艦船磁性設(shè)備的位置選定、磁性處理、消磁系統(tǒng)設(shè)計(jì)等提供幫助。

艦船感應(yīng)磁場(chǎng)的計(jì)算歸結(jié)為開域問題，而對(duì)于開域的靜磁場(chǎng)計(jì)算問題來說，由于積分方程法僅需對(duì)鐵磁區(qū)進(jìn)行離散，且不需要考慮邊界條件和可以求解連續(xù)函數(shù)等的優(yōu)點(diǎn)，受到廣泛關(guān)注和應(yīng)用。國(guó)內(nèi)，F(xiàn)an Mingwu等[10]提出了三維靜磁場(chǎng)的矢量積分方程法；倪振群等[11]將該方法應(yīng)用于艦船感應(yīng)磁場(chǎng)的計(jì)算；周國(guó)華等[12]、郭成豹等[13]、朱武兵等[14]使用非規(guī)則體單元、面單元對(duì)艦船幾何結(jié)構(gòu)進(jìn)行網(wǎng)格離散，研發(fā)了用于計(jì)算艦船感應(yīng)磁場(chǎng)的工具包；何保委等[15]使用標(biāo)量磁位建模提升了求解精度。針對(duì)在感應(yīng)磁場(chǎng)求解過程中鐵磁物體磁化率難以確定的問題，郭成豹等提出了利用試驗(yàn)樣品的磁性參數(shù)替代的方法[16]，周國(guó)華等則提出一種利用粒子群算法對(duì)磁化率進(jìn)行優(yōu)化的模型[17]，為艦船感應(yīng)磁場(chǎng)的正演計(jì)算奠定了基礎(chǔ)。國(guó)外，C. W. Trowbridge等[18]早在1970年間就將積分方程法用于求解物體的電磁場(chǎng)；O. Chadebec等[19]首次將積分方程法引入薄殼體艦船模型的閉環(huán)消磁算法建模當(dāng)中；S. Guerin[20]和Y. Vuillermet[21]研究對(duì)比了以磁化強(qiáng)度、點(diǎn)磁荷、切向磁偶極子為磁源的感應(yīng)磁場(chǎng)建模方法。近年來，法國(guó)Grenoble Alpes大學(xué)的研究團(tuán)隊(duì)陸續(xù)提出了新的關(guān)于積分方程法的加速改進(jìn)算法和耦合算法[22-24]，在減少網(wǎng)格數(shù)量的同時(shí)提升了磁場(chǎng)計(jì)算精度?？偟膩碚f，積分方程法已經(jīng)廣泛應(yīng)用于磁場(chǎng)計(jì)算領(lǐng)域，但是求解時(shí)間過長(zhǎng)和對(duì)計(jì)算機(jī)配置要求較高等問題還有待進(jìn)一步改善。

本文以艦船感應(yīng)磁場(chǎng)的計(jì)算為研究背景，提出一種改進(jìn)型標(biāo)量積分方程法。首先，通過選取三角形單元離散薄殼體，依據(jù)插值原理[25]將計(jì)算點(diǎn)由單元中心點(diǎn)轉(zhuǎn)移至節(jié)點(diǎn)，推導(dǎo)得到了基于標(biāo)量磁位的積分方程公式。其次，在求解稠密耦合系數(shù)矩陣的逆運(yùn)算時(shí)，為減少計(jì)算機(jī)的內(nèi)存需求和計(jì)算時(shí)間，引入了通用性較強(qiáng)的多層自適應(yīng)交叉近似（Multi- Level Adaptive Cross Approximation, MLACA）算法和廣義最小殘差（Generalized Minimal Residual, GMRES）法，將源點(diǎn)和場(chǎng)點(diǎn)按照相對(duì)位置進(jìn)行分組計(jì)算。由于感應(yīng)磁場(chǎng)正演需要已知鐵磁材料的磁性參數(shù)，考慮到臨近單元的磁化率應(yīng)當(dāng)平滑過渡。然后，基于磁場(chǎng)實(shí)測(cè)值和正演耦合算法建立了等效磁化率的反演模型，以磁場(chǎng)擬合度、磁化率先驗(yàn)分布和光滑約束條件構(gòu)建目標(biāo)函數(shù)，采用模擬退火（Simulated Annealing, SA）算法對(duì)磁化率尋優(yōu)求解，在具有磁場(chǎng)解析解的鐵質(zhì)球殼上進(jìn)行數(shù)值仿真對(duì)比驗(yàn)證。最后，將本文所提出的高精度感應(yīng)磁場(chǎng)建模方法應(yīng)用于一艘磁性參數(shù)未知的縮比船模進(jìn)行試驗(yàn)研究，通過將測(cè)量值和實(shí)測(cè)值進(jìn)行比較，最后得出結(jié)論。

1 基于標(biāo)量磁位的感應(yīng)磁場(chǎng)正演建模方法

1.1 標(biāo)量磁位積分公式

在艦船內(nèi)部，假設(shè)船體鐵磁材料各向同性，即

圖1 地磁場(chǎng)作用下的艦船離散單元模型和感應(yīng)磁場(chǎng)方向

基于疊加定理，空間任意一點(diǎn)的磁位為所有離散單元在該點(diǎn)產(chǎn)生的磁位的代數(shù)和。聯(lián)立式（3）～式（6），可以得到

1.2 離散單元的選取

（1）從離散單元來看，使用面單元要比體單元更加適應(yīng)薄殼體的剖分。

（2）對(duì)于同一艘艦船或者薄殼體，使用面單元離散后的單元數(shù)量相比較于體單元數(shù)量能夠大幅減少，從而提高求解速度、降低計(jì)算機(jī)的內(nèi)存需求。

不失一般性，采用如圖2所示的三棱柱作為離散單元，此時(shí)式（8）仍然適用，并且可以轉(zhuǎn)化為

圖2 三棱柱離散單元

考慮到薄殼體結(jié)構(gòu)的厚度e相比于長(zhǎng)度和寬度都小得多，因此可以使用三角形面單元對(duì)艦船的薄殼體結(jié)構(gòu)進(jìn)行離散，從而得到

（12）

1.3 三角形線性插值

觀察式（12）的被積函數(shù)，需要對(duì)未知的標(biāo)量磁位求梯度，因此需要消去該梯度算子?；谌切螁卧木€性插值原理[25]，可以將單元內(nèi)部的標(biāo)量磁位用插值函數(shù)和節(jié)點(diǎn)標(biāo)量磁位來表示，即

聯(lián)立式（12）～式（15），可以得到

式（5）和式（18）組成了基于面單元求解艦船感應(yīng)磁場(chǎng)的標(biāo)量磁位積分方程法公式。首先以地磁場(chǎng)作為左端的輸入項(xiàng)，通過求解該代數(shù)方程組可以得到每個(gè)離散單元的標(biāo)量磁位，然后再求解空間任意場(chǎng)點(diǎn)的磁場(chǎng)值。

1.4 計(jì)算耦合系數(shù)矩陣

式中，為局部坐標(biāo)系的軸與全局坐標(biāo)系的軸的夾角，=，=。全局坐標(biāo)系和局部坐標(biāo)系的位置關(guān)系如圖4所示。

至此，基于標(biāo)量磁位的薄殼體感應(yīng)磁場(chǎng)正演建模方法和公式已經(jīng)給出。為了體現(xiàn)標(biāo)量建模方法相比于矢量建模方法的優(yōu)越性，參照式（18），直接給出以磁感應(yīng)強(qiáng)度為待求量的積分方程作為后續(xù)對(duì)比研究，其表達(dá)式為

2 基于MLACA的加速求解方法

2.1 MLACA算法原理

MLACA為一種矩陣壓縮算法，它將離散單元?jiǎng)澐譃榻鼒?chǎng)組和遠(yuǎn)場(chǎng)組兩類，對(duì)于近場(chǎng)組的單元采用直接計(jì)算方法，對(duì)于遠(yuǎn)場(chǎng)組的單元采用兩個(gè)滿秩矩陣的乘積來近似表示，其原理為

為了使得近似矩陣盡量逼近精確矩陣，定義殘差計(jì)算公式為

2.2 基于MLACA的感應(yīng)磁場(chǎng)正演流程

MLACA算法的優(yōu)點(diǎn)在于它是純數(shù)學(xué)方法，不會(huì)涉及積分方程核函數(shù)的展開和變換，所以適用于各種不同核函數(shù)的積分方程組，通用性較強(qiáng)，直接給出MLACA算法應(yīng)用于標(biāo)量磁位積分方程法的磁場(chǎng)正演流程[29]。

1）網(wǎng)格離散。對(duì)鐵磁目標(biāo)進(jìn)行幾何建模和網(wǎng)格剖分，獲取單元和節(jié)點(diǎn)的信息。

2）建立八叉樹。將鐵磁目標(biāo)所占空間按照八叉樹模式劃分為多個(gè)小區(qū)域，以層數(shù)=2為例，鐵磁區(qū)域空間一共劃分成了64個(gè)小立方體，每個(gè)立方體包含了一組場(chǎng)點(diǎn)或源點(diǎn)的集合，若64組點(diǎn)的集合相互作用，得到的是一個(gè)64×64的矩陣塊。在該矩陣中，近場(chǎng)區(qū)為位于對(duì)角線的矩陣塊，是滿秩矩陣；遠(yuǎn)場(chǎng)區(qū)為非對(duì)角線的矩陣塊，是低秩矩陣。按照近場(chǎng)區(qū)和遠(yuǎn)場(chǎng)區(qū)的劃分，分別采用直接求解方法和調(diào)用MLACA算法求解系數(shù)矩陣。

3）對(duì)遠(yuǎn)場(chǎng)區(qū)調(diào)用MLACA算法。

4）采用廣義最小殘差（GMRES）法求解式（18），輸出值為單元節(jié)點(diǎn)的標(biāo)量磁位。

2.3 鐵質(zhì)球殼模型數(shù)值仿真

為了檢驗(yàn)標(biāo)量磁位積分方程法結(jié)合MLACA算法的應(yīng)用效果，本文給出了具有磁場(chǎng)解析解的鐵質(zhì)球殼算例進(jìn)行驗(yàn)證。為了排除網(wǎng)格疏密帶來的影響，分別采用不同的網(wǎng)格密度對(duì)球殼模型進(jìn)行剖分，用于檢驗(yàn)所提出算法的計(jì)算精度、運(yùn)行速度以及對(duì)計(jì)算機(jī)所消耗內(nèi)存的壓縮程度。

鐵質(zhì)球殼模型及測(cè)量點(diǎn)示意圖如圖5所示，外徑為10 m，厚度為0.01 m，材料磁導(dǎo)率為61。外磁場(chǎng)大小ez=36 000 nT，方向豎直向下。以球殼模型的中心點(diǎn)為原點(diǎn)建立三維直角坐標(biāo)系（豎直向下為軸正方向），測(cè)量點(diǎn)的范圍為∈[-30, 30 ] m，步長(zhǎng)為1 m；∈[-10, 10] m，步長(zhǎng)為10 m；=15 m。

球殼模型在垂向外磁場(chǎng)的作用下，三分量感應(yīng)磁場(chǎng)的解析結(jié)果可以表示為

圖5 鐵質(zhì)球殼模型及測(cè)量點(diǎn)示意圖

根據(jù)式（18）和式（22），分別使用基于MLACA的標(biāo)量積分方程法（Scalar Integral Equation Method, SIEM）和矢量積分方程法（Vector Integral Equation Method, VIEM）求解不同網(wǎng)格密度球殼下方測(cè)量點(diǎn)的感應(yīng)磁場(chǎng)。以磁場(chǎng)解析解為理論值，評(píng)估不同方法計(jì)算得到的磁場(chǎng)的平均相對(duì)誤差。

為保證使用不同算法進(jìn)行比較的公平性，本文使用的編程軟件統(tǒng)一為Matlab 2020，選用的計(jì)算機(jī)平臺(tái)為搭載12核心的Core i7 10 750H型處理器，主頻為2.6 GHz，RAM為32 GB的筆記本，計(jì)算結(jié)果見表1，當(dāng)網(wǎng)格密度模式為3號(hào)時(shí)，球殼模型的磁場(chǎng)三分量結(jié)果如圖6所示。

對(duì)于球殼模型來說，根據(jù)表1和圖6的數(shù)據(jù)可以得出以下結(jié)論：

（1）網(wǎng)格密度越高，積分方程法的磁場(chǎng)計(jì)算誤差越小，說明適當(dāng)增加離散單元的數(shù)量可以得到高精度的鐵磁物體感應(yīng)磁場(chǎng)，但是當(dāng)網(wǎng)格密度增加到一定程度時(shí)，計(jì)算精度提高的不明顯。

表1 使用積分方程法計(jì)算不同網(wǎng)格密度球殼模型磁場(chǎng)的結(jié)果

Tab.1 The calculation results of spherical shell models with different mesh densities utilizing integral equation method

（2）本文使用了MLACA算法壓縮內(nèi)存（層數(shù)為6）和加速計(jì)算后，SIEM+MLACA相比較VIEM節(jié)省了內(nèi)存消耗約97%，節(jié)省了計(jì)算時(shí)間約65%。

圖6 鐵質(zhì)球殼磁場(chǎng)計(jì)算結(jié)果比較

（3）在相同網(wǎng)格密度的情況下，基于MLACA的SIEM在保證計(jì)算精度不變的前提下，能夠大幅減少計(jì)算時(shí)間和降低計(jì)算機(jī)的內(nèi)存需求，驗(yàn)證了MLACA算法的有效性。

3 等效磁化率反演優(yōu)化方法

3.1 等效磁化率反演模型

根據(jù)積分方程法原理可知，計(jì)算物體感應(yīng)磁場(chǎng)的前提是獲得材料的磁化率（或磁導(dǎo)率）。對(duì)于一般的鐵磁性材料來說，可以通過查閱磁化率曲線或者采取樣品試驗(yàn)的方法來獲取磁化率[17]，但由于艦船鐵磁材料的多樣性以及船體需要經(jīng)過多次的高溫鍛造和沖擊試驗(yàn)，意味著艦船磁化率不是一個(gè)常數(shù)，而是一個(gè)與空間位置相關(guān)的變量。

首先假設(shè)在某一磁化狀態(tài)下個(gè)測(cè)量點(diǎn)獲取的艦船三分量感應(yīng)磁場(chǎng)為

在反演模型中，通常根據(jù)磁場(chǎng)計(jì)算值和測(cè)量值的擬合程度來尋找合適的磁化率分布，令擬合函數(shù)為

根據(jù)式（27）和式（28）可以得到眾多符合截止條件的空間磁化率分布，然而對(duì)于物體“真實(shí)”的磁化率來說，該方法得到的解可能是眾多局部最優(yōu)解的其中一個(gè)。為確保得到擬合度和穩(wěn)定性較好的等效磁化率，加入兩個(gè)約束條件來減少解空間的維度，從而改善反演結(jié)果的質(zhì)量。

（2）約束條件2：等效磁化率的分布應(yīng)當(dāng)符合物理客觀規(guī)律和滿足平滑性準(zhǔn)則，即不允許某一單元的磁化率無限大，也不允許相鄰單元的磁化率發(fā)生突變。具體實(shí)現(xiàn)辦法是使用有限差分格式[30]使得相鄰單元的磁化率過渡更加平滑，具體表達(dá)式為

3.2 基于SA算法的等效磁化率求解方法

模擬退火算法借鑒于高溫固體的退火過程，當(dāng)固體的溫度較高時(shí)，由于內(nèi)能較大，固體內(nèi)部的粒子處于快速且無序的運(yùn)動(dòng)狀態(tài)[31]。在溫度緩慢降低的過程中，由于內(nèi)能減小，粒子緩慢趨于有序。最終，在固體處于常溫時(shí)，內(nèi)能達(dá)到最小，此時(shí)粒子狀態(tài)最為穩(wěn)定。SA算法由兩層循環(huán)組成，內(nèi)外層循環(huán)分別是同一溫度下尋找局部最優(yōu)解和溫度冷卻過程中尋找全局最優(yōu)解，它是通過在內(nèi)循環(huán)搜索解的過程中賦予一定的概率來跳出局部最優(yōu)解從而達(dá)到全局尋優(yōu)的目的[32]。下面直接給出由SA算法對(duì)磁化率進(jìn)行反演建模的流程，如圖7所示。

圖7 模擬退火法用于反演磁化率的流程

將SIEM+MLACA方法獲得的鐵質(zhì)球殼磁場(chǎng)三分量作為輸入，使用SA算法對(duì)鐵磁球殼的磁化率進(jìn)行反演優(yōu)化，使用和不使用約束條件前后的磁化率光滑度對(duì)比結(jié)果如圖8所示?？梢钥闯?，增加了5%的光滑度約束之后，相鄰單元的磁化率過渡更加平滑，說明了約束方程的作用。

圖8 優(yōu)化前后的磁化率曲線光滑度對(duì)比

基于反演優(yōu)化得到的等效磁化率對(duì)球殼在25 m深度的五條測(cè)線的磁場(chǎng)進(jìn)行預(yù)測(cè)，結(jié)果如圖9所示。與“真實(shí)磁化率”計(jì)算得到的磁場(chǎng)解析值進(jìn)行對(duì)比，預(yù)測(cè)磁場(chǎng)的平均相對(duì)誤差僅為2.2%，證明所獲取的等效磁化率是相對(duì)準(zhǔn)確的。

圖9 基于等效磁化率的球殼模型磁場(chǎng)預(yù)測(cè)值和解析值比較（左舷-龍骨-右舷）

4 船模試驗(yàn)研究

4.1 試驗(yàn)方案設(shè)計(jì)

為驗(yàn)證所提出的等效磁化率反演優(yōu)化方法在工程上的實(shí)用，以一磁性參數(shù)未知的縮比船模（見圖1）為試驗(yàn)研究對(duì)象，通過測(cè)量其感應(yīng)磁場(chǎng)求解得到等效磁化率分布，并將其用于預(yù)測(cè)其他深度的感應(yīng)磁場(chǎng)，從而驗(yàn)證等效磁化率的穩(wěn)定性。

船模磁場(chǎng)測(cè)量示意圖如圖10所示，船模長(zhǎng)度為=2.75 m，寬度為=0.5 m，厚度為=0.8 mm，在船模下方設(shè)置上層磁傳感器陣列1和下層磁傳感器陣列2，分別作為磁場(chǎng)測(cè)量平面和磁場(chǎng)預(yù)測(cè)平面。此外，為了排除環(huán)境磁場(chǎng)的影響，在遠(yuǎn)離船模的地方設(shè)置了單個(gè)傳感器用于監(jiān)測(cè)環(huán)境磁場(chǎng)的變化，所使用的傳感器為高精度三分量磁通門傳感器，測(cè)量分辨率為1 nT。船模所在試驗(yàn)環(huán)境的磁場(chǎng)水平分量為34 200 nT，垂向分量為35 000 nT，通過逐點(diǎn)測(cè)量法獲取船模下方的垂向感應(yīng)磁場(chǎng)，測(cè)量點(diǎn)的范圍為∈[-,]，步長(zhǎng)為/60；∈[-,]，步長(zhǎng)為；1=1.3，2=2。

圖10 磁性船模測(cè)量磁場(chǎng)測(cè)量示意圖

4.2 船模試驗(yàn)結(jié)果

將船模在1平面測(cè)量的磁場(chǎng)值作為輸入，通過反演優(yōu)化模型求解得到其等效磁化率分布，根據(jù)磁場(chǎng)擬合值判斷等效磁化率是否合理，然后將其用于預(yù)測(cè)2平面的船模磁場(chǎng)。船模在1平面的磁場(chǎng)擬合效果和2平面的磁場(chǎng)預(yù)測(cè)值和實(shí)測(cè)值的對(duì)比結(jié)果如圖11所示。

由圖11看出，磁場(chǎng)擬合值和預(yù)測(cè)值都和實(shí)測(cè)值比較吻合，其誤差分別為5.1%和4.8%，表明本文所提出方法得到的等效磁化率是有效的，可用于艦船感應(yīng)磁場(chǎng)的高精度建模。從整個(gè)測(cè)量和計(jì)算的過程來看，試驗(yàn)結(jié)果的誤差主要來源以下方面：

（1）測(cè)量誤差：主要包括測(cè)傳感器自身的測(cè)量誤差、使用逐點(diǎn)測(cè)量法測(cè)量船磁時(shí)的人為操作引入的誤差以及試驗(yàn)環(huán)境變化帶來的偶然誤差。

（2）計(jì)算誤差：考慮到計(jì)算時(shí)間不宜過長(zhǎng)，船模的剖分單元不能無限多，但理論上只有離散單元足夠小、足夠多時(shí)，每個(gè)單元才可以被看作均勻磁化體，但是這種做法會(huì)導(dǎo)致計(jì)算量急劇增加。

（3）模型誤差：船模受地磁場(chǎng)磁化時(shí)，實(shí)際磁化率仍會(huì)有微小變化，但是本文為了求解等效磁化率，認(rèn)為磁化率在弱磁環(huán)境中基本不變，這種線性化處理方式引入了誤差。

5 結(jié)論

本文首先提出了一種結(jié)合多層自適應(yīng)交叉近似算法的改進(jìn)型標(biāo)量積分方程法，用于快速高效獲取高精度的艦船感應(yīng)磁場(chǎng)。其次通過建立以磁場(chǎng)擬合度、磁化率先驗(yàn)分布和平滑約束條件為目標(biāo)函數(shù)的等效磁化率反演模型，并使用模擬退火算法進(jìn)行優(yōu)化求解，解決了工程中鐵磁物體磁化率不易獲取的難題，得出以下結(jié)論：

1）球殼數(shù)值仿真表明，所提出的鐵磁物體感應(yīng)磁場(chǎng)正演模型高效準(zhǔn)確。當(dāng)幾何模型的網(wǎng)格密度較小時(shí)，基于MLACA法的標(biāo)量積分方程法微弱提高了磁場(chǎng)計(jì)算精度；隨著網(wǎng)格密度增大，該方法在保持計(jì)算精度的同時(shí)，大幅減少了計(jì)算時(shí)間和計(jì)算機(jī)的內(nèi)存消耗。

2）球殼數(shù)值仿真和縮比船模試驗(yàn)研究結(jié)果表明，等效磁化率反演模型有效可行，其中船模試驗(yàn)結(jié)果的磁場(chǎng)擬合誤差和預(yù)測(cè)誤差約為5.0%，表明所提出的感應(yīng)磁場(chǎng)正演模型和等效磁化率反演優(yōu)化模型可以用于大型艦船的高精度磁場(chǎng)建模，能夠?yàn)榕灤F磁設(shè)備安裝位置的確定和消磁繞組的設(shè)計(jì)提供 依據(jù)。

[1] Holmes J J. Reduction of a ship's magnetic field signatures[M]. Maryland, USA: Morgan & Claypool Publishers, 2008.

[2] 劉大明. 艦船消磁理論與方法[M]. 北京: 國(guó)防工業(yè)出版社, 2011.

[3] 郭成豹, 殷琦琦. 艦船磁場(chǎng)磁單極子陣列法建模技術(shù)[J]. 物理學(xué)報(bào), 2019, 68(11): 106-115.

Guo Chengbao, Yin Qiqi. Magnetic monopole array model for modeling ship magnetic signatures[J]. Acta Physica Sinica, 2019, 68(11): 106-115.

[4] 戴忠華, 周穗華, 張曉兵. 多目標(biāo)優(yōu)化的艦船磁場(chǎng)建模方法[J]. 物理學(xué)報(bào), 2021, 70(16): 147-159.

Dai Zhonghua, Zhou Suihua, Zhang Xiaobing. Multi- objective optimization of ship magnetic field mode- ling method[J]. Acta Physica Sinica, 2021, 70(16): 147-159.

[5] 劉琪, 孫兆龍, 姜潤(rùn)翔, 等. 一種艦船下方磁場(chǎng)的信號(hào)重構(gòu)及換算方法[J]. 電工技術(shù)學(xué)報(bào), 2022, 37(15): 3723-3732.

Liu Qi, Sun Zhaolong, Jiang Runxiang, et al. A signal reconstruction and conversion method of magnetic field under ship[J]. Transactions of China Electro- technical Society, 2022, 37(15): 3723-3732.

[6] 劉芙妍, 顏冰. 磁偶極子陣列模型的適用性研究與優(yōu)化分析[J]. 物理學(xué)報(bào), 2022, 71(12): 36-48.

Liu Fuyan, Yan Bing. Applicability and optimization analysis of magnetic dipole array model[J]. Acta Physica Sinica, 2022, 71(12): 36-48.

[7] 徐慶林, 王向軍, 張建春, 等. 921A鋼板腐蝕電場(chǎng)的Frumkin修正[J]. 電工技術(shù)學(xué)報(bào), 2020, 35(14): 2951-2958.

Xu Qinglin, Wang Xiangjun, Zhang Jianchun, et al. Frumkin correction of corrosion electric field gen- erated by 921A steel[J]. Transactions of China Electrotechnical Society, 2020, 35(14): 2951-2958.

[8] 饒凡, 吳旭升, 高嵬, 等. 兩極永磁電機(jī)靜態(tài)內(nèi)外磁場(chǎng)研究[J]. 電工技術(shù)學(xué)報(bào), 2021, 36(14): 2936- 2944.

Rao Fan, Wu Xusheng, Gao Wei, et al. Study on internal and external magnetic field of static two-pole permanent magnet motor[J]. Transactions of China Electrotechnical Society, 2021, 36(14): 2936-2944.

[9] Chadebec O, Coulomb J L, Janet F. A review of magnetostatic moment method[J]. IEEE Transactions on Magnetics, 2006, 42(4): 515-520.

[10] Fan Mingwu, Shao Hanguang, Wang Jingguo. Some experiences of using integral equation method to calculate magnetostatic fields[J]. IEEE Transactions on Magnetics, 1985, 21(6): 2185-2187.

[11] 倪振群, 蔡雪祥, 翁行泰. 艦船主甲板模型感應(yīng)磁場(chǎng)的計(jì)算[J]. 上海交通大學(xué)學(xué)報(bào), 1996, 30(7): 83-88.

Ni Zhenqun, Cai Xuexiang, Weng Xingtai. Com- putation of inductive magnetic field induced by vessel’s main deck model[J]. Journal of Shanghai Jiaotong University, 1996, 30(7): 83-88.

[12] 周國(guó)華, 肖昌漢, 劉勝道, 等. 基于六面體單元表面磁場(chǎng)積分法求解三維靜磁場(chǎng)[J]. 電工技術(shù)學(xué)報(bào), 2009, 24(3): 1-7.

Zhou Guohua, Xiao Changhan, Liu Shengdao, et al. 3D magnetostatic field computation with hexahedral surface integral equation method[J]. Transactions of China Electrotechnical Society, 2009, 24(3): 1-7.

[13] 郭成豹, 劉大明. 薄鋼殼物體磁特征建模研究[J]. 兵工學(xué)報(bào), 2012, 33(8): 912-915.

Guo Chengbao, Liu Daming. Modeling of magnetic signatures of thin-sheet objects[J]. Acta Armamentarii, 2012, 33(8): 912-915.

[14] 朱武兵, 莊勁武, 趙文春, 等. 載體感應(yīng)磁場(chǎng)的改進(jìn)積分方程法[J]. 國(guó)防科技大學(xué)學(xué)報(bào), 2018, 40(3): 101-106.

Zhu Wubing, Zhuang Jinwu, Zhao Wenchun, et al. Modified integral equation method for carrier’s induced magnetic field solving[J]. Journal of National University of Defense Technology, 2018, 40(3): 101- 106.

[15] 何保委, 周國(guó)華, 劉勝道, 等. 標(biāo)量磁位積分方程法計(jì)算艦艇三維靜磁場(chǎng)[J/OL]. 兵工學(xué)報(bào), 2022, http://kns.cnki.net/kcms/detail/11.2176.TJ.20221027.1405.004.html.

He Baowei, Zhou Guohua, Liu Shengdao, et al. Ship’s 3D magnetostatic field computation utilizing integral equation method of scalar magnetic potential[J]. Acta Armamentarii, 2022, http://kns.cnki.net/kcms/detail/ 11.2176.TJ.20221027.1405.004.html.

[16] 郭成豹, 張曉鋒, 肖昌漢, 等. 薄鋼板船舶磁性的物理模型和數(shù)值模擬混合建模[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2008, 29(3): 247-250.

Guo Chengbao, Zhang Xiaofeng, Xiao Changhan, et al. Magnetization determination of ships made of thin steel sheets by combining the physical model and numerical calculations[J]. Journal of Harbin Engin- eering University, 2008, 29(3): 247-250.

[17] 周國(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 a weak magnetic field[J]. Journal of Harbin Engineering University, 2009, 30(1): 91-95.

[18] Trowbridge C W. Progress in magnet design by computer[C]//Proceedings of 4th Conference Mag- netostatic Technology, NewYork, USA, 1972: 555- 565.

[19] Chadebec O. Modelling of magnetic field induced by shells-identification of the magnetisation application to the closed loop degaussing of a ferromagnetic hull[D]. Grenoble, France: Institut National Polyte- chnique de Grenoble, 2001.

[20] Guerin S. Magnetic source identification measurement robustness and optimization- application to ship magnetization reconstruction[D]. Grenoble, France: Institut National Polytechnique de Grenoble, 2005.

[21] Vuillermet Y. Closed loop degaussing applied to double hull submarine Magnetization identification from near magnetic field measurements[D]. Grenoble, France: Institut National Polytechnique de Grenoble, 2008.

[22] Le-Van V. Developement of magnetostatic volume integral formulations[D]. Grenoble, France: Univer- site Grenoble Alpes, 2015.

[23] Chavin-Collin G, Bannwarth B, Cavallera D, et al. An integral face formulation for thin non-conductive magnetic regions[J]. IEEE Transactions on Magnetics, 2019, 55(6): 1-4.

[24] Chavin-Collin G. Identification de l’aimantation de t?les ferromagnétiques soumises à des champs magnétiques faibles et à de fortes contraintes mécaniques: application à l’immunisation magnétique des sous-marins en boucle fermée[D]. Grenoble, France: Universite Grenoble Alpes, 2020.

[25] 謝德馨, 楊仕友. 工程電磁場(chǎng)數(shù)值分析與優(yōu)化設(shè)計(jì)[M]. 2版. 北京: 機(jī)械工業(yè)出版社, 2017.

[26] Le-Van V, Bannwarth B, Carpentier A, et al. The adaptive cross approximation technique for a volume integral equation method applied to nonlinear magnetostatic problems[J]. IEEE Transactions on Magnetics, 2014, 50(2): 445-448.

[27] 周耀忠, 張國(guó)友. 艦船磁場(chǎng)分析計(jì)算[M]. 北京: 國(guó)防工業(yè)出版社, 2004.

[28] 陳杰, 魯習(xí)文. 艦船感應(yīng)磁場(chǎng)預(yù)測(cè)的一種新方法[J]. 物理學(xué)報(bào), 2010, 59(1): 239-245.

Chen Jie, Lu Xiwen. A new method for predicting the induced magnetic field of naval vessels[J]. Acta Physica Sinica, 2010, 59(1): 239-245.

[29] 裴亞杰. 電磁場(chǎng)積分方程自適應(yīng)交叉近似算法的研究[D]. 南京: 南京郵電大學(xué), 2011.

[30] Vijn A R P J. Development of a closed-loop degaussing system towards magnetic unobservable vessels[D]. Delft, Netherlands: Delft University of Technology, 2021.

[31] 何保委, 劉勝道, 周國(guó)華, 等. 基于多目標(biāo)模擬退火法的垂向檢測(cè)線圈改進(jìn)[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2021, 49(12): 46-50.

He Baowei, Liu Shengdao, Zhou Guohua, et al. Improvement of vertical detecting coil based on multi-objective simulated annealing method[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2021, 49(12): 46-50.

[32] 陳龍, 易瓊洋, 賁彤, 等. 全局優(yōu)化算法在Preisach磁滯模型參數(shù)辨識(shí)問題中的應(yīng)用與性能對(duì)比[J]. 電工技術(shù)學(xué)報(bào), 2021, 36(12): 2585-2593, 2606.

Chen Long, Yi Qiongyang, Ben Tong, et al. Appli- cation and performance comparison of global optimization algorithms in the parameter identi- fication problems of the Preisach hysteresis model[J]. Transactions of China Electrotechnical Society, 2021, 36(12): 2585-2593, 2606.

A Fast Forward Modeling Method for High Precision Induced Magnetic Field of Ships

（School of Electrical Engineering Naval University of Engineering Wuhan 430033 China）

Mastering the distribution of a ship’s induced magnetic field is an important issue for implementing magnetic stealth technology, and the integral equation method is one of the main methods to calculate the ship’s induced magnetic field. The integral equation method only needs to discretize the ferromagnetic region and does not need to consider the boundary conditions, so it has been widely concerned and applied. The traditional vector integral equation for modeling induced magnetic fields has the problems of low efficiency and an enormous computational burden. Considering that there are many discrete elements of large ferromagnetic objects such as ships, the coupling coefficient between elements forms a huge asymmetric dense matrix. As a result, the computing time and memory requirement will increase sharply with the increase of the number of elements. Therefore, a scalar magnetic potential integral equation method based on surface elements is proposed. Since the integral equation method needs to obtain the magnetic susceptibility of ferromagnetic materials, the equivalent magnetic susceptibility inversion model is established based on a multi-level adaptive cross approximation (MLACA) algorithm.

Firstly, the scalar magnetic potential integral formula based on triangular surface elements is derived. According to the principle of linear interpolation, the scalar magnetic potential of the center point of discrete elements is expressed by interpolation function and node scalar magnetic potential, and the elements’ coupling coefficient matrix is obtained by establishing a local coordinate system. Therefore, the scalar integral method of surface elements is obtained for solving the ship’s induced magnetic field. Secondly, the MLACA algorithm is introduced to guarantee the accuracy of magnetic field calculation and reduce the memory requirement and computing time of the computer. Finally, aiming at the problem that the magnetic parameters of ferromagnetic materials of ships are not easy to obtain, a magnetic susceptibility inversion model is established based on the measured magnetic field values and the forward coupling model. The magnetic field fitting degree, prior distribution of magnetic susceptibility, and smooth constraint are the objective function. The spatial distribution of equivalent magnetic susceptibility is optimized by simulated annealing (SA) algorithms.

A numerical simulation of the spherical iron shell shows that the proposed scalar magnetic potential coupling forward modeling method can efficiently obtain the ship’s induced magnetic field with high precision. For the same discrete elements, compared with the vector method, the scalar method can save about 97% of the memory consumption and 65% of the computing time, which verifies the effectiveness of the scalar magnetic potential integral equation method based on the MLACA algorithm. According to the inversion model, the equivalent magnetic susceptibility of the spherical iron shell is optimized by the SA algorithm, which is utilized to predict the spherical shell of other positions, and the average relative error is only 2.2%. After using of smooth constraint condition, the equivalent magnetic susceptibility obtained is smoother.

In order to verify the practicability of the proposed forward modeling method and magnetic susceptibility inversion method in engineering, an experimental scheme was designed for a reduced-scale ship model with unknown magnetic parameters. Firstly, the equivalent magnetic susceptibility distribution of the ship model is obtained by measuring the magnetic field of the ship model at the1 plane, and it is used to predict the induced magnetic field at the2 plane. The magnetic field fitting and prediction errors of the ship are about 5.0%, indicating that the proposed forward modeling model of induced magnetic field and the inverse optimization model of equivalent magnetic susceptibility can be used for induced magnetic field modeling of large ships with high precision and can provide support for the implementation of magnetic stealth technology on ships.

Ship magnetic field, reduced scalar magnetic potential, multi-level adaptive approximate cross, simulated annealing method, equivalent magnetic susceptibility

10.19595/j.cnki.1000-6753.tces.222204

TM153

國(guó)家自然科學(xué)基金資助項(xiàng)目（51377165, 52207020）。

2022-11-22

2022-12-29

何保委 男，1995年生，博士研究生，研究方向?yàn)殡姶怒h(huán)境與防護(hù)技術(shù)。E-mail: hebaowei188@163.com

劉月林 女，1989年生，碩士，講師，研究方向?yàn)殡姶怒h(huán)境與防護(hù)技術(shù)。E-mail: mooncake1024@163.com（通信作者）

（編輯 陳 誠(chéng)）