基于磁力梯度全張量特征值的均衡邊界識別方法

鄭 強 郭 華 張貴賓 韓 松 王 明 劉浩軍

(①中國地質大學(北京)地球物理與信息技術學院，北京 100083； ②中國自然資源航空物探遙感中心，北京 100083)

0 引言

邊界識別是位場(重磁)數據處理與解釋的重要內容，其對圈定異常范圍，劃分斷裂構造等具有重要意義。基于重磁異常的傳統導數類邊界識別方法，如垂向導數法[1-3]、總水平導數法[4-5]、解析信號法[6-9]等，均存在不能清晰識別深部地質體邊界的缺點。為了提高對深部地質體邊界的識別效果，地球物理學家提出了均衡邊界識別方法。Miller等[10]提出了第一種均衡邊界識別方法——傾斜角法，該方法分析垂向導數與總水平導數的比值，基于其零等值線推斷地質體邊界。傾斜角法能很好地均衡深、淺部地質體邊界，但識別精度較低。Verduzco等[11]利用傾斜角的總水平導數進行地質體的邊界識別，該方法對淺部地質體的邊界識別效果較好，且其受地質體傾斜側面的影響較小，但對邊界的增強效果不明顯。Wijns等[12]提出Theta法，該方法計算總水平導數與解析信號的比值，根據最大值的分布定義地質體邊界。Theta法受場源埋深影響較小，但識別精度不高，且輸出結果穩定性較差。Cooper等[13]提出了水平傾斜角法，根據水平傾斜角的極大值解釋地質體邊界。該方法能同時識別深、淺部異常體的邊界，但也存在識別精度不高的問題，且對噪聲極其敏感。馬國慶等[14]提出了基于總水平導數與垂向導數的相關系數法，利用相關系數的極小值識別地質體邊界。該方法具有一定的抗噪性，但識別效果受滑動窗口影響較大。Ma等[15]提出利用歸一化總水平導數的極大值識別地質體邊界。該方法有效提高了輸出結果的穩定性，但識別結果同樣依賴于滑動窗口的選擇。馬國慶等[16]提出增強型均衡濾波器，該濾波器利用不同階導數之間的組合識別地質體邊界。該方法能有效識別深、淺部地質體邊界，但高階導數的引入使噪聲不可避免地被放大，從而影響整體邊界識別的穩定性?？傊卮女惓５膫鹘y導數類邊界識別方法大多存在識別精度低和輸出結果穩定性差等缺點。

隨著重磁勘探技術的不斷發展，目前主要進行重磁梯度全張量測量。重磁梯度全張量數據具有高精度、高分辨率、多參量的優點[17]，成為近年來位場數據處理與解釋的研究熱點[18-21]。大量學者對重磁梯度全張量的不變量[22-23]、特征值[24]以及總模值[25]進行了研究，并提出了相應的邊界識別方法，一定程度上提高了地質體邊界的識別效果，但還沒有一種既能夠提高識別精度、又能提高結果穩定性的邊界識別方法。

基于磁力梯度全張量的最大特征值，提出了一種新的均衡邊界識別方法。該方法不僅能夠同時識別深淺模型體的邊界，而且具有較高的識別精度、一定的抗傾斜磁化能力和較強的穩定性，是一種能夠同時提高識別精度和輸出結果穩定性的邊界識別新方法。

1 基本理論

1.1 磁力梯度全張量特征值

磁力梯度全張量G可表示為

(1)

式中：U表示磁位；Bx、By、Bz為磁場矢量的三個分量；Bmn(m,n=x,y,z)是磁場張量的9個分量；為了方便下文公式推導，將磁場張量的9個分量記為Tij(i,j=1,2,3)。

利用下式求取式(1)的特征值λ

=(T11-λ)[(T22-λ)(T33-λ)-T23T32]-

T12[T21(T33-λ)-T23T31]+

T13[T21T32-T31(T22-λ)]=0

(2)

式中： det(·)表示求取矩陣行列式；E表示單位矩陣。

T11T22T33-T12T23T31-T13T21T32=0

(3)

式(3)可等效為如下形式的一元三次方程

λ3+pλ+q=0

(4)

式中

(5)

T12T23T31-T13T21T32

(6)

利用卡爾丹公式，可得式(4)的三個解分別為

式中

(10)

(11)

磁力梯度全張量矩陣的最大特征值λ1對應地質體邊界[24]，但由于其主要構成元素為磁場張量，而張量信號隨著地質體埋深的增加會急劇衰減，因此該方法對深部地質體的分辨率不高。

1.2 均衡邊界識別方法

磁力梯度全張量的總模值M為

(12)

M含有全部9個張量元素的信息，其最大值對應著地質體的邊界[25]，但是其識別精度較低。為了增強有用信號的幅值，本文提出將磁力梯度全張量矩陣的最大特征值與總模值相乘

S=λ1M

(13)

相乘的優點是使淺部邊界處異常幅值較大的有用信號增強，缺點是使深部邊界處異常幅值較小的有用信號減弱，因此S能夠提高對淺部地質體邊界的識別精度，但無法有效識別深部地質體的邊界，為此，對其進行了均衡處理

(14)

式中k是均衡深淺部異常的調節系數，一般取值為0～1，其大小與地質體異常信號的幅度有關，需要解釋人員通過模型和實際數據調試。

1.3 對比邊界識別方法

為了說明本文提出的邊界識別方法的有效性，選取三種基于磁力梯度全張量的邊界識別方法進行對比，即不同分量組合法(HG)、總水平解析信號法(TA)和方向解析信號法(BDA)，其參數分別定義如下

式中

(18)

分別為方向解析信號。

2 模型試驗

為了檢驗本文方法的邊界識別效果，分別建立了簡單模型和復雜模型進行計算分析。

2.1 簡單模型

首先建立由兩個立方體組成的簡單模型(圖1)。兩個立方體的邊長均為200m，中心點坐標(x，y，z)分別為(-200m，200m，120m)和(200m，-200m，180m)，頂面埋深分別為20、80m。

下面分別考慮垂直磁化、傾斜磁化和加噪垂直磁化三種情形分析本文方法對異常體邊界的識別效果。

2.1.1 垂直磁化

地質體產生的重磁異常會隨著地質體埋深的增加而迅速衰減，對深部地質體產生的弱異常信號進行有效提取是邊界識別的主要問題，因此埋深是影響邊界識別效果的主要因素之一。

圖1 簡單模型俯視圖(a)和平視圖(b)

首先分析埋深對邊界識別效果的影響。對圖1所示模型設置其地磁傾角為90°，地磁偏角為0°，磁化傾角為90°，磁化偏角為0°，磁化強度為1A/m。

圖2是垂直磁化下磁力梯度全張量6個元素的平面等值線圖(由于張量矩陣的對稱性，有3組分量完全相同，因此只展示每組中的一個，下同)。可以看出，較淺的異常體的異常信息較豐富，所反映的模型體邊界很清晰；反之，較深異常體的異常信息比較弱，所反映的異常體邊界較模糊。

圖3是垂直磁化情況下各邊界識別參數或識別方法計算結果。從圖3a～圖3c可以看出，λ1和M的最大值可以識別淺部異常體邊界，但識別精度較低，且對深部異常體邊界無法識別。相比于λ1和M，S對淺部異常體的識別精度更高，但仍無法有效識別深部異常體邊界。從圖3d～圖3f可以看出，相比于參數S，均衡邊界識別方法的參數BS能夠識別深部異常體邊界，調節系數k起到均衡深淺部異常的作用：大的k值對深部模型邊界的識別均衡能力強，但識別精度低；小的k值對深部模型邊界的均衡能力弱，但識別精度高。模型試驗表明，k取0.001時可兼顧深部和淺部異常體的邊界識別效果。從圖3g～圖3i可以看出： 參數HG的總體識別效果一般，且邊界處異常分布發散； 參數TA對淺部模型的識別效果較好，但無法有效識別深部模型體的邊界； 參數BDA能夠識別深部異常體邊界，但總體識別精度不高。

圖2 垂直磁化情況下磁力梯度全張量6個元素的平面等值線圖

圖3 垂直磁化情況下不同參數的平面分布圖

由上分析可知，本文均衡邊界識別方法在垂直磁化情況下，能夠很好地識別深部和淺部異常體邊界，調節系數k取0.001時，有著較高的識別精度，可兼顧深部和淺部異常體。

2.1.2 傾斜磁化

實際情況中，地質體所處的磁場環境往往都是傾斜磁化的，而傾斜磁化會使地質體產生的異常形態變得復雜，不再具有對稱性，因此傾斜磁化也是影響邊界識別效果的因素之一。本節分析傾斜磁化對邊界識別效果的影響。對圖1所示模型，設地磁傾角為60°，地磁偏角為5°，磁化傾角為60°，磁化偏角為5°，磁化強度為1A/m。

圖4是傾斜磁化情況下磁力梯度全張量6個元素的平面等值線圖。相比于圖2，圖4中各張量元素所體現的異常形態不再具有對稱性，與模型體邊界的實際位置也存在一定的偏移。

圖4 傾斜磁化情況下磁力梯度全張量6個元素的平面等值線圖

圖5是傾斜磁化情況下各邊界識別參數或識別方法計算結果的平面等值線圖。相比于圖3，圖5中λ1、M以及S受傾斜磁化影響較小，但總體來說識別的邊界向南移動；均衡邊界識別方法的參數BS也體現了相似特征，即具有一定的抗傾斜磁化能力；參數HG、TA及BDA受傾斜磁化的影響較大，對深部和淺部異常體邊界均不能有效識別，尤其是BDA受傾斜磁化的影響最嚴重，產生了錯誤的邊界信息。

總體上看，本文提出的均衡邊界識別方法在傾斜磁化情況下，具有一定的抗干擾能力，對傾斜磁化敏感度較低，可對深部和淺部異常體邊界進行識別。

2.1.3 加噪垂直磁化

除了異常體埋深和傾斜磁化會對邊界識別產生影響，測量環境和測量儀器等帶來的各種噪聲也會對邊界識別結果產生較大干擾，尤其是基于重磁異常的傳統導數類邊界識別方法會大幅度放大高頻噪聲，無法有效地識別邊界。

因此，本節討論噪聲對邊界識別效果的影響。對垂直磁化情況下的9個張量元素，分別加入各自異常最大幅值10%的高斯噪聲。

圖6是含噪垂直磁化情況下磁力梯度全張量6個元素的平面等值線圖。相比于圖2，圖6中各張量元素的異常形態由于噪聲的干擾變得十分紊亂，無法反映模型體的邊界位置。

圖5 傾斜磁化情況下不同參數的平面分布圖

圖6 含噪垂直磁化情況下磁力梯度全張量6個元素的平面等值線圖

圖7是含噪垂直磁化情況下不同邊界識別參數或識別方法計算結果。可以看出，參數λ1、M的識別效果受到了噪聲的影響，但參數S受噪聲的影響較小； 參數BS對于調節系數k取值不同情況下受噪聲影響的大小不同：k值越大，受噪聲的影響越大，k值越小，受噪聲的影響越??；基于磁力梯度全張量的邊界識別方法中，參數HG具有一定的抗噪性，而參數TA和BDA對噪聲極其敏感，尤其是BDA受噪聲的影響最為嚴重。

圖7 含噪垂直磁化情況下不同參數的平面分布圖

(a)λ1； (b)M最大值； (c)S； (d)BS：k=0.100； (e)BS：k=0.010； (f)BS：k=0.001； (g)HG； (h)TA； (i)BDA

因此，總體上來看，本文提出的均衡邊界識別方法在加噪垂直磁化情況下，當調節系數k取0.001時，受噪聲的影響較小，具有較強的魯棒性。

2.2 復雜模型

為了進一步說明本文方法的有效性并分析如何選取恰當的k值，下面建立一個包括四個立方體的復雜模型(表1和圖8)進行計算。

對復雜模型僅考慮垂直磁化不加噪的情況，并選取不同的k值，對比分析邊界識別效果，從而得到最佳k值。

圖9是復雜模型垂直磁化情況下調節系數k不同取值時利用本文均衡邊界識別方法得到的平面等值線圖。可以看出，隨著k取值由大到小，邊界識別的效果也在變化： 圖9a中k=0.1000，此時對異常體的邊界識別精度不高； 圖9b中k=0.0100，識別的邊界精度較圖9a而有所提高； 圖9c中k=0.0010時，邊界識別精度較圖9b有所提高； 圖9d中當k=0.0001時，幾乎沒有異常顯示，即無法識別異常體邊界。

表1 復雜模型參數表

圖8 復雜模型俯視圖(a)和平視圖(b)

圖9 復雜模型垂直磁化情況下的BS平面等值線圖

綜上，可以得出： 利用本文提出的均衡邊界識別方法能夠完成對復雜模型體邊界的識別，且當調節系數k取0.0010時，取得的邊界識別效果最好。本文提出的方法理論及k值選取的合理性得到了進一步驗證。

3 實際應用

由于磁力梯度全張量測量對測量儀器的要求很高，目前很少能夠在地面直接測量較高精度的磁力梯度全張量數據，但是航磁三分量的測量技術已發展成熟，國內以中國自然資源航空物探遙感中心為代表，已獲得了高精度的航磁三分量數據，而對航磁三分量數據求取三維空間各方向上的一階導數，便能獲得航磁梯度全張量數據。

為了檢驗本文均衡邊界識別方法的實際應用效果，截取了中國北部某地區的實測航磁三分量數據，并對其進行求導，進而得到航磁梯度全張量數據(圖10)。

對轉換后的航磁梯度全張量數據(圖10)進行了邊界識別方法的計算，結果見圖11。可以看出，λ1識別的邊界信息不夠清晰；M識別出的邊界只包含淺部地質體的信息；S識別結果精度有所提高，但仍無法顯示深部地質體的邊界信息； 參數BS隨著調節系數k取值的變化，識別結果越來越清晰，當k=0.001時識別出精度較高的地質體邊界；而基于磁力梯度全張量的三種邊界識別方法中，參數HG識別的地質體邊界不清晰；TA和BDA識別出的邊界十分模糊，對異常體的刻畫不夠清晰。

圖10 研究區磁力梯度全張量6個元素的平面等值線圖

圖11 研究區不同參數的平面分布圖

綜上，本文提出的均衡邊界識別方法對于實際航磁數據的應用具有良好的邊界識別效果，能夠獲得較清晰的邊界信息。

4 結論

本文基于磁力梯度全張量的最大特征值和總模值，提出了一種新的均衡邊界識別方法。模型試驗和實際應用證明，利用本文提出的均衡邊界識別方法，當調節系數k=0.0010時，能夠有效識別地質體邊界，并具一定的抗傾斜磁化能力和較強的抗噪性，是一種既能夠提高識別精度、又能提高計算穩定性的磁場數據邊界識別新方法。