斜磁化條件下多磁源目標的邊界識別方法

李金朋， 范紅波， 張英堂， 李志寧， 尹剛

(1.陸軍工程大學石家莊校區 車輛與電氣工程系， 河北 石家莊 050003；2.中國空氣動力研究與發展中心 高速空氣動力研究所， 四川 綿陽 621000)

0 引言

邊界識別是解釋磁測數據的重要過程，其主要目的是通過獲得的磁測數據對地下或水下的未知磁源的物理邊緣進行可視化，在軍事偵察(未爆彈、潛艇和水雷等)、水文及工程地質勘探等領域有著廣泛的應用[1-5]。

傳統的邊界識別方法主要包括總水平導數法、傾斜角法以及Theta 圖(TM)法等[6-9]，但是這些傳統的邊界識別方法存在一定局限性：總水平導數法利用磁場沿z軸方向的分量數據Bz在水平方向上梯度的模值對水平邊界進行計算，該方法不能同時對強弱磁異常進行顯示；傾斜角法計算Bz分量沿z軸方向空間導數獲得的Bzz與總水平導數法模值的反正切，獲得的弧度值為磁性目標的二維邊界結果，該方法能夠有效地平衡不同強弱磁異常的振幅，但是磁源的邊界沒有非常明確的定義；TM法計算總水平導數法獲得的模值與Bz分量在3個方向上梯度模值比值的反余弦獲得弧度值，利用弧度值結果能夠同時描繪不同磁性強弱的磁源，但是計算精度難以滿足對磁源邊界進行清晰描述的要求。在傾斜磁化條件下，利用傳統邊界識別方法獲得的磁源邊界將會偏離磁源的真實位置[10-12]。此外，在噪聲影響下，磁源的識別精度會受到一定影響。

盡管目前已經開始將磁梯度張量信息引入到磁源的二維邊界反演計算中，但是多局限于對垂直磁化條件下的磁源進行計算[13]，而基于磁梯度張量斜磁化條件下的多磁源目標二維邊界反演技術仍然處于起步階段。且在多磁源目標邊界計算過程中，測量得到的磁分量信號受環境磁場、系統誤差以及人為操作不當等因素會對磁測數據引入噪聲，嚴重影響計算結果。

因此，針對多磁源目標二維邊界反演中存在的問題，本文提出基于磁梯度張量的多磁源邊界識別方法。該方法首先基于改進K奇異值分解(K-SVD)的方法對磁測信號進行降噪處理，以實現對磁測數據的預處理，然后利用互相關理論對不同計算區域內的單目標磁源磁化方向進行估計，以實現將磁測數據轉化為垂直磁化條件下的數據，最后計算磁梯度張量矩陣特征值的絕對值之和與垂直磁化條件下的磁梯度張量Bzz分量的比值，獲得磁源的水平邊界。

1 磁梯度張量理論

在笛卡爾坐標系下，磁梯度張量數據是磁場B=[Bx,By,Bz]T沿著x軸、y軸、z軸3個方向的空間導數,Bx、By和Bz為磁場三分量數據。假設x軸指向地理北向，y軸指向地理東向，z軸指向垂直向下的方向，在觀測平面(x,y,z)的磁梯度張量數據可表示為

(1)

式中：Γ為磁梯度張量矩陣；Bij(i,j=x,y,z)為磁梯度張量分量，磁梯度張量矩陣為對稱矩陣，其矩陣的跡為0. 因此，磁梯度張量具有Bxx、Bxy、Bxz、Byy、Byz這5個獨立分量，其他分量都可以從這5個獨立分量推導出來。磁梯度張量矩陣的特征值分解可以表示為

(2)

式中：λ1、λ2、λ3為磁梯度張量矩陣的特征值，按照單調遞減的順序依次排列為|λ1|≥|λ2|≥|λ3|≥0，v1、v2、v3分別為特征值λ1、λ2、λ3對應的特征向量.

為了獲得磁梯度張量數據，設計了一個由4個三軸磁通門組成的十字形磁梯度張量系統，系統如圖1所示。兩個磁力計在同一方向上的基線距離為2d，d為傳感器測量軸到磁梯度張量系統中心之間的距離。十字形磁梯度張量系統的測量數據可以定義為

(3)

式中：Bij(i=1,2,3,4,j=x,y,z)為第i個傳感器j軸的分量數據。

圖1 十字形磁梯度張量系統Fig.1 Cross magnetic gradient tensor system

該測量系統的結構誤差[14]可以表示為

(4)

根據(4)式可以看出，構造的磁梯度張量系統測量值與理論值間的誤差較小，可以在實際測量中忽略不計。

2 多磁源二維邊界反演方法

2.1 基于改進K-SVD的磁測數據預處理方法

在對磁源的探測與識別過程中，常常需要考慮噪聲的影響，因此在計算過程中要降低磁場數據中的噪聲。在降噪過程中，為了保證不同方向上磁場分量之間的相互關系，將傳感器上不同的磁場分量利用矢量合成原理轉化為總幅值異常以及方向余弦進行降噪處理：

(5)

式中：B為磁總場幅值異常；l、m、n分別為沿坐標軸3個方向的方向余弦。

K-SVD降噪模型[15]可以表示為

(6)

式中：αe為最終的稀疏向量解；y是測試數據，y=b+q，b是理論數據，q是均值為0的高斯白噪聲；α∈Rs為K-SVD字典矩陣的稀疏向量；D∈Ra×s(s≥a)為字典矩陣，其中s為字典大小，a為測試數據個數；T為容差極限。在實際應用中，噪聲為標準高斯白噪聲的情況很少，因此，K-SVD降噪方法在應用中存在一定的局限性。因此，本文提出了一種改進K-SVD(IK-SVD)降噪模型來重建磁數據。

(7)

式中：yi是測試數據中的第i列數據；Mi是觀測數據去均值化后的第i列數據。

(8)

利用IK-SVD降噪方法對磁總場幅值異常B及方向余弦(l,m,n)進行降噪處理，利用(5)式對信號進行逆向重構[16]：

(9)

式中：B′x、B′y、B′z分別為降噪后的磁測數據三分量；B′為降噪后的磁總場幅值異常；(l′,m′,n′)為降噪后的方向余弦。

2.2 磁梯度張量數據轉化理論

在斜磁化條件下，對磁性目標進行邊界識別計算時會受到磁化方向的影響而對計算結果產生干擾。為了解決這一問題，本文利用基于歸一化磁源強度(NSS)和化極(RTP)互相關的方法對磁化方向進行估計，進而利用獲得的磁化方向對降噪后的磁測數據進行RTP計算，最終獲得垂直磁化條件下磁源的磁梯度張量數據。NSS是由磁梯度張量矩陣推導出來的一個旋轉不變量[17]，可以表示為

(10)

定義實測區域內的磁異常RTP結果ΔTr與NSS數值μ的互相關系數C[18]為

(11)

磁總場異常TMI RTP公式[18]為

(12)

磁場三分量的傅里葉變換[19]可以表示為

(13)

2.3 磁源邊界識別理論

在對多目標磁源進行邊界計算過程中，利用磁梯度張量矩陣獲得的特征值的和(SAE)對磁源的分布區域進行劃分,這是因為SAE數據的主正極值與磁源的實際分布位置相一致，受磁化方向影響較小。同時，為了獲得磁源的邊界，利用特征值邊緣檢測的方法來進行計算，邊界計算結果Es的具體計算公式如下：

(14)

在垂直磁化條件下，磁梯度張量Bzz分量的正負數據分布與磁源的邊界一致，因此，磁源的邊界位于Es的正負數據交界處，即磁源的邊界可以用Es的零輪廓線來進行描述。

3 仿真分析

為了證明本文方法的有效性，在測量區域內放置兩個斜磁化條件下的模型進行仿真分析，如圖2所示。地磁背景場的磁化方向為Ib=60°、Db=-20°. 模型1為多邊形薄板，磁矩為50 A/m，磁化方向為Ia1=60°、Da1=-20°，厚度為0.1 m，中心深度為0.25 m. 模型2為長方形薄板，磁矩為50 A/m，磁化方向為Ia2=45°、Da2=55°，厚度為0.1 m,中心深度為0.35 m.

圖2 組合模型Fig.2 Combined model

圖3 觀測面上磁測數據Fig.3 Magnetic measurement data on the observation surface

分別向4個傳感器每個軸上加入信噪比為30%的高斯白噪聲，模擬環境因素帶來的誤差，并向觀測數據中隨機加入0.5×rand(max(B))的沖擊噪聲，模擬實際測量中由于操作誤差所帶入的沖擊誤差。觀測面上獲得的磁測數據如圖3所示，圖3(h)中的白線為模型的理論水平位置，圖3(i)中的白線為劃分的計算區域。為了證明本文降噪方法的有效性，分別利用形態學濾波器、低通濾波器、小波變換方法、EMD分解方法以及K-SVD方法與本文的IK-SVD方法進行對比，其中形態學濾波器通過構造開- 閉和閉- 開組合的組合濾波器并利用三角形結構元素進行濾波運算，低通濾波器采用巴特沃斯濾波器確定截止頻率，EMD濾波器與小波濾波器采用軟閾值的濾波方法[16]，小波方法中的小波基為coif5，分解層數為5層。K-SVD計算方法中，塊大小為4，訓練字典大小為512. 本文方法的窗口大小為4,移動步長為1. 分別利用信噪比SNR以及均方根誤差RMSE作為評判降噪能力的指標，其中

表1 不同計算方法的降噪結果Tab.1 Denoised results of different calculation methods nT/m

圖4 去噪聲后觀測面上磁測數據Fig.4 Denoised magnetic measurement data on the observation surface

利用IK-SVD方法降噪后的磁測數據如圖4所示。利用(11)式估計兩個計算窗口內模型的磁化方向，其中模型1的磁化方向估計值為Is1=56°、Ds1=-18°，模型2的磁化方向估計值為Is2=48°、Ds2=59°. 圖5為磁梯度張量轉化前后不同計算區域內的Bzz數據。利用(13)式獲得垂直條件下的磁梯度張量Bzz分量，如圖5(c)和圖5(d)所示，白色實線為磁源的理論水平分布區域。與圖5(a)和圖5(b)對比可以看出：在垂直磁化條件下，實線的位置與Bzz主正極值的分布更加一致；而在斜磁化條件下，Bzz主正極值與實際輪廓位置存在一定的偏移。因此，利用本文的相關性來估計磁化方向的方法可以有效地降低斜磁化的影響。

圖6 不同方法獲得的邊界識別結果Fig.6 Edge detection results obtained by different methods

圖7 降噪及磁梯度張量轉化前后的Es輪廓Fig.7 Es values before and after noise reduction and magnetic gradient tensor conversion

4 試驗驗證

為了驗證上述方法的有效性，進行實測試驗。測區位于石家莊某地，探測儀器及磁源位置分布如圖8所示。試驗所用探頭是利用英國Bartington公司生產的Mag-03三軸磁通門傳感器搭建的十字形磁梯度張量探頭，將傳感器固定在塑料支架上，兩個磁傳感器在同一方向上的基線距離為0.4 m，在試驗過程中，探頭固定在無磁試驗臺架上，利用掃描的方式對待測區域的每一點進行測量，利用非線性最小二乘方法對傳感器安裝誤差及軟磁干擾進行校正及補償[20]。測區的背景磁化方向為Ib=55°、Db=-16°，測點的間距為0.1 m，在測區內對兩個長方形鐵塊的磁測數據進行測量，獲得磁測數據如圖9所示。

圖8 試驗測試儀器及磁源分布Fig.8 Experimental equipment and magnetic source distribution

圖9(h)中鐵塊1尺寸為沿南北方向長度為0.36 m，沿東西方向寬度為0.23 m，中心位于(1.15 m，0.65 m)，頂部深度為0.2 m，厚度為0.1 m. 圖9(h)中鐵塊2尺寸為沿南北方向長度為0.36 m，沿東西方向寬度為0.234 m，中心為(0.45 m，1.5 m)。頂部深度為0.15 m，厚度為0.05 m. 分別利用形態學濾波器、低通濾波器、小波變換方法、EMD分解方法以及K-SVD方法與本文的IK-SVD方法進行對比，計算結果如圖10所示。由圖10可以看出，通過對獲得的磁梯度張量5個獨立分量進行對比，可以看出本文降噪方法具有最優的降噪能力。將SAE圖劃分為覆蓋各個鐵塊的計算區域，如圖9(i)中的實線所示。

利用IK-SVD方法降噪后的磁測數據如圖11所示。利用(11)式估計兩個計算窗口內鐵塊的磁化方向，其中鐵塊1的磁化方向估計值為Ie1=12°、De1=-5°，鐵塊2的磁化方向估計值為Ie2=3°、De2=-13°. 圖12為磁梯度張量轉化前后不同計算區域內的Bzz數據。利用(13)式獲得垂直條件下的磁梯度張量Bzz分量，如圖12(c)和圖12(d)所示，白色實線為磁源的理論水平分布區域。與圖12(a)和圖12(b)對比可以看出：在垂直磁化條件下，實線的位置與Bzz主正極值的分布更加一致；而在斜磁化條件下，Bzz主正極值與實際輪廓位置存在一定的偏移。圖13為具有噪聲的斜磁化條件下不同計算方法識別的邊界輪廓。圖14為及磁梯度張量轉化前后的Es輪廓。由圖13和圖14可以看出，本文方法能夠有效提高多磁源目標在噪聲及斜磁化條件下的邊界識別能力，進一步證明了本文方法的有效性。

圖9 觀測面上磁測數據Fig.9 Magnetic measurement data on the observation surface

圖10 不同計算方法的降噪結果Fig.10 Denoised results of different calculation methods

圖11 去噪聲后觀測面上磁測數據Fig.11 Denoised magnetic measurement data on the observation surface

圖12 磁梯度張量轉化前后不同計算區域內的Bzz數據Fig.12 Bzz data in different computing areas before and after magnetic gradient tensor conversion

圖13 不同方法獲得的邊界識別結果Fig.13 Edge detection results obtained by different methods

圖14 降噪及磁梯度張量轉化前后的Es輪廓Fig.14 Es values before and after noise reduction and magnetic gradient tensor conversion

5 結論

本文針對磁源在二維邊界反演過程中存在的問題，提出了斜磁化條件下多磁源目標二維邊界識別方法，該方法能有效地在噪聲以及斜磁化條件下對多磁源目標的二維邊界進行識別。通過仿真和試驗證明了本文方法具有如下優勢：

1)提出了利用IK-SVD方法對磁測數據總幅值異常以及方向余弦進行預處理，能夠在預處理過程中能夠有效保證磁測數據三分量之間的內在聯系，且具有較強的抗噪性。

2)利用NSS與RTP互相關的方法估計不同計算區域內磁源的磁化方向，將磁測數據轉化為垂直磁化條件下的磁測數據，能夠有效消除斜磁化對計算結果帶來的影響。

3) 通過計算SAE與垂直條件下的磁梯度張量Bzz分量的比值零線位置，能夠準確對磁源的水平邊界進行求解。

本文方法適用的探測范圍可以從米到千米，具有一定的國防和工業應用背景，未來的工作旨在檢測剩磁以提高其準確性，該成果可應用于未爆彈探測、礦產勘探等領域。