改進(jìn)二維變分模態(tài)分解的磁源分離

李金朋，任國全，張英堂，范紅波，李志寧

(陸軍工程大學(xué)石家莊校區(qū) 車輛與電氣工程系，河北 石家莊 050003)

1 引 言

磁測數(shù)據(jù)解釋主要應(yīng)用于軍事偵察(未爆彈、潛艇和水雷等)、缺陷檢測、剩磁估計(jì)及工程地質(zhì)勘探等領(lǐng)域[1-5]。在實(shí)際測量中獲得的磁測數(shù)據(jù)包含了從淺部到深部所有磁源共同引起的磁異常，包含了豐富了磁源信息。當(dāng)需要對某個(gè)目標(biāo)進(jìn)行解釋時(shí)(如對淺部未爆彈進(jìn)行探測)，需要對磁測數(shù)據(jù)進(jìn)行分離，將特定目標(biāo)引起的磁異常從疊加異常中分離出來，進(jìn)行下一步的磁測數(shù)據(jù)解釋。因此，對磁源分離方法進(jìn)行研究是當(dāng)前研究的熱點(diǎn)問題之一。

當(dāng)前，磁源的分離方法主要參考重力數(shù)據(jù)上的分離方法，向上延拓法是應(yīng)用較為廣泛的一種方法[6-7]，通過對觀測數(shù)據(jù)進(jìn)行向上延拓，作為異常目標(biāo)的區(qū)域異常，將觀測數(shù)據(jù)與區(qū)域異常數(shù)據(jù)進(jìn)行做差求解最終的局部異常。該方法需要已知最佳向上延拓高度，在向上延拓的過程中，淺源信號被壓制的同時(shí)，深源信號也會(huì)產(chǎn)生一部分損失。Zeng等[6]提出了向上延拓法的最佳延拓高度估計(jì)法，計(jì)算一系列的相鄰延拓面磁測數(shù)據(jù)的互相關(guān)系數(shù)，互相關(guān)系數(shù)曲線的最大偏差點(diǎn)對應(yīng)的高度為最佳延拓高度。Pawlowski[7]提出了優(yōu)選向上延拓法，該方法向上延拓時(shí)能夠在對淺源信號進(jìn)行壓制的同時(shí)對深源信號長波不衰減。Meng等[8]在優(yōu)選向上延拓法的基礎(chǔ)上，提出了優(yōu)選向上延拓的差值法，實(shí)現(xiàn)對某一特定頻段信號進(jìn)行提取。Guo[9-10]提出了基于格林等效層概念和維納濾波的優(yōu)化濾波方法，該方法不需要計(jì)算延拓高度。上述方法需要計(jì)算對實(shí)測數(shù)據(jù)的平均對數(shù)功率譜進(jìn)行分段擬合來確定頻段范圍，在實(shí)際應(yīng)用中存在一定難度。在實(shí)際測量中，當(dāng)區(qū)域異常較大時(shí)，對局部異常分離精度會(huì)有所降低，使得對磁性目標(biāo)的識別結(jié)果產(chǎn)生一定影響，因此，對不同深度磁源磁測數(shù)據(jù)的分離方法進(jìn)行研究仍處于起步階段。為提高磁源的分離能力，提出了利用改進(jìn)二維變分模態(tài)分解法對磁源進(jìn)行分離。本文在最佳延拓高度估計(jì)法基礎(chǔ)上，將延拓法獲得的淺源局部異常利用二維變分模態(tài)分解進(jìn)行二次分離，進(jìn)一步去除延拓法存留的區(qū)域異常數(shù)據(jù)，提高磁源局部異常數(shù)據(jù)的精度，并轉(zhuǎn)化為磁梯度張量數(shù)據(jù)實(shí)現(xiàn)對不同深度的磁性目標(biāo)進(jìn)行有效識別。

2 理論分析

2.1 基于最佳延拓高度的向上延拓方法

進(jìn)行向上延拓計(jì)算時(shí)，需要對最佳延拓高度進(jìn)行估計(jì)，通過計(jì)算相鄰延拓面上的互相關(guān)系數(shù)，繪制互相關(guān)系數(shù)與延拓高度的互相關(guān)圖像，其中磁源的最佳向上延拓高度位于互相關(guān)圖像的最大偏轉(zhuǎn)處[6]。

相鄰延拓面之間的互相關(guān)公式C為：

(1)

其中：M和N分別為沿x和y軸的觀測點(diǎn)數(shù)，B1和B2分別為兩個(gè)相鄰延拓面上的磁測數(shù)據(jù)。

2.2 二維變分模態(tài)分解算法

二維變分模態(tài)分解(2D-Variational Mode Decomposition,TVMD)理論是由一維變分模態(tài)分解理論推導(dǎo)而來[11-12]，主要應(yīng)用于生物化學(xué)和納米科學(xué)領(lǐng)域。該方法從經(jīng)驗(yàn)?zāi)B(tài)分解方法(Empirical Mode Decomposition,EMD)發(fā)展而來[13]，通過對構(gòu)造的變分模型進(jìn)行尋優(yōu)獲得不同的子信號，使得子信號的帶寬之和最小。在二維解析信號中，需要將頻域一半的平面設(shè)置為零，該半平面是相對于矢量選擇的，為wk。二維解析信號的頻率域定義如下:

uAS,k(w)=(1+sgn(〈w,wk〉))uk(w),

(2)

其中：{wk}={w1,w2,…,wK}為模態(tài)分量uk的中心頻率，k=1,2,…K，K為模態(tài)分量個(gè)數(shù)。

具有上述傅里葉特性的二維解析信號為[14]:

uAS,k(η)=uk(η)*(δ(〈η,wk〉)+

j/π〈η,wk〉)δ(〈η,wk,⊥〉),

(3)

式中*代表卷積且變換是可分的。

二維變分模態(tài)分解的函數(shù)最小化公式為:

(4)

其中：αk為Tikhonov正則化參數(shù)，f(x)為觀測數(shù)據(jù)。

利用二次懲罰項(xiàng)以及拉格朗日乘數(shù)計(jì)算約束問題，用乘法交替方向算法(Alternating Direction Method of Multipliers,ADMM)[15,16]進(jìn)行優(yōu)化計(jì)算。對模態(tài)分量uk進(jìn)行優(yōu)化計(jì)算:

(5)

式中λ為拉格朗日乘數(shù)項(xiàng)。求解uk會(huì)產(chǎn)生維納濾波的更新結(jié)果:

?w∈Ωk:Ωk={w|〈w,wk〉≥0}.

(6)

對中心頻率wk進(jìn)行優(yōu)化計(jì)算:

(7)

因此，在半平面Ωk模型的功率譜|uk(w)|2表示為:

(8)

通過公式(6)以及公式(8)，能夠?qū)ΧS觀測信號做變分模態(tài)分解得到K個(gè)固有模態(tài)分分量。在本文計(jì)算過程中，由于需要假設(shè)信號是由深源和淺源組合而成，因此，需要定義K=2。

利用二維變分模態(tài)分解計(jì)算磁源分量，能夠自適應(yīng)的對兩個(gè)頻帶內(nèi)的觀測數(shù)據(jù)進(jìn)行分離，不需要對觀測數(shù)據(jù)進(jìn)行人為的進(jìn)行分段擬合，提高了自適應(yīng)性。通過利用二維變分模態(tài)分解對延拓法分離出來淺源局部異常進(jìn)行計(jì)算，能夠有效保留區(qū)域異常的磁源信息，提高局部異常的分離精度。

2.3 磁梯度張量理論及應(yīng)用

在笛卡爾坐標(biāo)系下，假設(shè)x軸正方向?yàn)楸毕?，y軸正方向?yàn)闁|向，z軸正方向?yàn)榇怪毕蛳?。磁梯度張量?shù)據(jù)是由磁場矢量B=[Bx,By,Bz]T沿坐標(biāo)軸3個(gè)正交方向的空間變化率，磁梯度張量可以由1個(gè)二階矩陣表示為[17-19]:

(9)

其中：G為磁梯度張量；Bi(i=x,y,z)為磁場三分量；Bαβ(α,β=x,y,z)為磁梯度張量的9個(gè)分量。磁梯度張量矩陣G為對稱矩陣，因此包含5個(gè)獨(dú)立分量。

利用磁梯度張量矩陣Bzx以及Bzy分量能夠?qū)Υ旁吹乃椒植歼M(jìn)行計(jì)算，計(jì)算公式為：

(10)

在計(jì)算過程中，利用Bz異常獲得的磁梯度張量分量的計(jì)算公式為[20]:

(11)

2.4 計(jì)算步驟

根據(jù)以上分析，可以利用本文方法對磁源進(jìn)行分離，并獲得磁梯度張量數(shù)據(jù)，計(jì)算過程如圖1所示，計(jì)算要點(diǎn)可歸納為：

(1)計(jì)算磁異常Bz分量數(shù)據(jù)的最佳向上延拓高度h，利用延拓法對數(shù)據(jù)進(jìn)行分離，獲得局部異常BJ和區(qū)域異常BQ；

(2)對局部異常BJ數(shù)據(jù)進(jìn)行變分模態(tài)分解計(jì)算，獲得局部異常BJQ和區(qū)域異常BJJ；

(3)利用分離出來局部異常以及區(qū)域異常，通過頻率域計(jì)算方法獲得磁梯度張量分量數(shù)據(jù)并計(jì)算磁源分布位置。

圖1 磁測數(shù)據(jù)分離步驟Fig.1 Separation steps of magnetic data

3 仿真分析

為了證明本文方法的有效性，建立5個(gè)長方體模型，模型形狀、位置及物性參數(shù)如表1所示。假設(shè)觀測表面的深度為0 m，觀測點(diǎn)的間距為1 m，觀測點(diǎn)網(wǎng)格為66×66。理論模型具有兩個(gè)深度層模型A和B，深度層模型A產(chǎn)生的磁場模擬區(qū)域異常數(shù)據(jù)，深度層模型B產(chǎn)生的磁場模擬局部異常數(shù)據(jù)，多目標(biāo)磁源的空間分布如圖2所示。對不同深度層模型以及所有深度層模型的磁異常Bz分量數(shù)據(jù)進(jìn)行仿真計(jì)算，獲得的計(jì)算結(jié)果如圖3所示。黑線表示深度層模型A的水平分布位置，白線表示深度層模型B的水平分布位置，紅色虛線代表測線位置。根據(jù)圖3可以看出，區(qū)域異常具有較高的磁場大小。對觀測面上所有組合模型的磁異常Bz分量數(shù)據(jù)進(jìn)行觀察可知，深度層模型B的磁測數(shù)據(jù)由于磁性較弱，幾乎被區(qū)域異常數(shù)據(jù)所淹沒。

表1 模型位置及物理參數(shù)Tab.1 Location and physical parameters of models

圖2 多目標(biāo)磁源空間分布圖Fig.2 Spatial distribution of multi-target magnetic sources

圖3 磁源磁異常Bz分量理論數(shù)據(jù)Fig.3 Theoretical data of magnetic anomaly Bz component

利用公式(11)對本文方法得到的區(qū)域異常與局部異常的磁梯度張量Bxz與Byz分量進(jìn)行計(jì)算，獲得的磁梯度張量分量的區(qū)域異常數(shù)據(jù)與局部異常數(shù)據(jù)如圖6(a)～圖6(f)所示。通過對比分離前后的磁測數(shù)據(jù)可知，通過本文的分離方法，能夠?qū)⒉煌疃却旁吹拇艤y數(shù)據(jù)進(jìn)行有效地分離。對分離得到的局部異常與區(qū)域異常數(shù)據(jù)的THDR進(jìn)行計(jì)算，獲得的計(jì)算結(jié)果如圖6(g)～圖6(i)所示。

圖4 測線上不同方法的延拓結(jié)果Fig.4 Continuation results of different methods on survey line

圖5 磁源磁異常Bz分量數(shù)據(jù)分離結(jié)果Fig.5 Separation results of magnetic source magnetic anomaly Bz component

表2 不同分離方法的分離結(jié)果
Tab.2 Separation results of different separation methods

方法磁測數(shù)據(jù)無噪聲RMSE(×104)/nT互相關(guān)有噪聲RMSE(×104)/nT互相關(guān)延拓法區(qū)域異常1.170 40.997 81.172 00.997 7局部異常1.170 40.493 01.172 00.491 4變分模態(tài)分解區(qū)域異常0.764 20.703 70.764 50.703 5局部異常0.764 20.150 20.764 50.149 9本文方法區(qū)域異常0.048 90.998 40.067 90.997 9局部異常0.048 90.926 60.067 90.859 6

圖6 磁梯度張量Bxz, Byz以及THDRFig.6 Magnetic gradient tensors Bxz, Byz and THDR

根據(jù)計(jì)算的結(jié)果可以看出，在分離之前由于區(qū)域異常數(shù)據(jù)的磁源強(qiáng)度較大，導(dǎo)致局部異常數(shù)據(jù)的分布被淹沒，難以得到局部數(shù)據(jù)實(shí)際的分布位置。而利用本文分離數(shù)據(jù)進(jìn)行計(jì)算，可以有效獲得不同深度磁源的實(shí)際分布范圍。同時(shí)，在區(qū)域異常較強(qiáng)的情況下，能夠?qū)植慨惓５姆植挤秶M(jìn)行有效求解。

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

在河北省石家莊某地進(jìn)行實(shí)測實(shí)驗(yàn)測試，測試裝置如圖7所示。該測試裝置系統(tǒng)主要包括：英國Bartingtong公司制造的Mag-03MS傳感器，數(shù)字采集模塊、軟件操作終端以及無磁實(shí)驗(yàn)臺架。在實(shí)驗(yàn)過程中，傳感器固定在無磁實(shí)驗(yàn)臺架上，利用掃描的方式對待測區(qū)域的每一點(diǎn)進(jìn)行測量，其中傳感器的量程為-70 000～70 000 nT之間，傳感器分辨率為0.01 nT。測區(qū)為2.1×2.1 m的區(qū)域，測量間隔為0.1 m。在測區(qū)內(nèi)放置一個(gè)南北走向的水平圓柱體鐵管以及一個(gè)南北走向水平放置的未爆彈。其中，水平圓柱體的軸線距離觀測面為0.58 m作為區(qū)域異常，未爆彈的軸線距離觀測面為0.32 m作為局部異常。分別對單個(gè)磁源Bz分量數(shù)據(jù)以及組合磁源的Bz分量之和進(jìn)行測試，獲得的數(shù)據(jù)如圖8所示。在圖8中，模型1代表水平圓柱體，黑色實(shí)線代表模型1的水平分布位置；模型2代表未爆彈，白色實(shí)線代表模型2的水平分布位置，紅色虛線代表測線位置(彩圖見期刊電子版)。

圖7 實(shí)驗(yàn)裝置Fig.7 Experimental equipment

分別利用本文分離方法、延拓法及變分模態(tài)分解法對測線上磁源的Bz分量數(shù)據(jù)進(jìn)行分離計(jì)算，獲得的分離結(jié)果如圖9所示。圖10為不同分離方法獲得分離結(jié)果的二維數(shù)據(jù)。對分離數(shù)據(jù)計(jì)算得到的RMSE以及互相關(guān)系數(shù)C如表3所示。根據(jù)圖9、圖10以及表3可以看出，本文分離方法計(jì)算精度最好，延拓法的計(jì)算精度其次，變分模態(tài)分解結(jié)果最差。根據(jù)計(jì)算結(jié)果可以看出，在延拓高度不高的情況下，延拓法獲得的磁源深度層數(shù)據(jù)損失較小，因此計(jì)算精度有所提高。本文計(jì)算方法獲得的計(jì)算結(jié)果互相關(guān)系數(shù)不低于0.966 4，RMSE不高于7 506 nT，優(yōu)于其他兩種方法。利用公式(11)對本文方法得到的水平圓柱體與未爆彈的磁梯度張量Bxz與Byz分量進(jìn)行計(jì)算，獲得的磁梯度張量分量的區(qū)域異常數(shù)據(jù)與局部異常數(shù)據(jù)如圖11(a)～圖11(f)所示。計(jì)算局部異常與區(qū)域異常數(shù)據(jù)的THDR，計(jì)算結(jié)果如圖11(g)～圖11(i)所示。根據(jù)計(jì)算結(jié)果可以看出，在分離之前由于水平圓柱體與未爆彈數(shù)據(jù)混合在一起，導(dǎo)致在計(jì)算結(jié)果中難以區(qū)分未爆彈的分布位置。通過利用本文分離數(shù)據(jù)進(jìn)行計(jì)算，可以有效獲得不同深度磁源的實(shí)際分布范圍。實(shí)驗(yàn)表明，本文方法能夠?qū)⒔M合磁源中的未爆彈進(jìn)行有效分離，在實(shí)踐中能夠應(yīng)用于未爆彈的探測。同時(shí)，該方法的適用范圍可以從幾米擴(kuò)展到幾千米，因此對于礦物以及石油勘探也同樣適用。

圖8 磁源磁異常Bz分量數(shù)據(jù)Fig.8 Magnetic anomaly Bz component of sources

圖9 測線上不同方法的延拓結(jié)果Fig.9 Continuation results of different methods on survey line

圖10 磁源磁異常Bz分量數(shù)據(jù)分離結(jié)果Fig.10 Separation results of magnetic source magnetic anomaly Bz component

表3 不同分離方法的分離結(jié)果

Tab.3 Separation results of different separation methods

方法磁測數(shù)據(jù)RMSE(×104)互相關(guān)延拓法區(qū)域異常1.675 40.9887局部異常1.675 40.9613變分模態(tài)分解區(qū)域異常1.870 60.970 6局部異常1.870 60.944 3本文方法區(qū)域異常0.750 60.990 9局部異常0.750 60.966 4

圖11 磁梯度張量Bxz, Byz以及THDRFig.11 Magnetic gradient tensors Bxz, Byz and THDR

4 結(jié) 論

在區(qū)域異常較大的條件下，局部異常數(shù)據(jù)不易被發(fā)現(xiàn)，利用改進(jìn)的二維變分模態(tài)分解理論對局部異常具有較好的分離效果。

(1) 利用最佳延拓高度估計(jì)法對磁源進(jìn)行一次分離，該過程能夠直接對磁測數(shù)據(jù)進(jìn)行分離并獲得局部異常與區(qū)域異常數(shù)據(jù)；

(2)對延拓后的局部異常進(jìn)行二維變分模態(tài)分解，能夠有效提高區(qū)域異常與局部異常的計(jì)算精度；

(3)實(shí)驗(yàn)結(jié)果表明：本文方法能夠?qū)^(qū)域異常較大條件下的局部異常數(shù)據(jù)進(jìn)行有效分離，對小尺度磁性體組合磁源(高度差為26 cm)的分離數(shù)據(jù)與單目標(biāo)觀測數(shù)據(jù)的互相關(guān)系數(shù)在0.966 4以上。該方法為后期不同深度磁源的定位、邊界識別以及三維反演提供了理論基礎(chǔ)與工程經(jīng)驗(yàn)。