重磁勘探地質異常體邊界識別方法研究及對比

唐新功，侯帥，李丹丹，2

1.油氣資源與勘探技術教育部重點實驗室(長江大學)，湖北 武漢 430100

2.武漢大學中國南極測繪研究中心，湖北 武漢 430079

作為最成熟的地球物理位場勘探方法之一，重磁勘探在油氣勘探、固體礦產勘探、地下水、工程物探、軍事、大地測量、深部構造以及地震預報中都發揮著不可或缺的重要作用[1]。在重磁勘探中，實測重磁數據包含了地下不同規模、不同深度各種地質異常體組合產生的疊加異常，為了對目標地質異常體進行勘探，必須在疊加重磁異常中提取出目標異常。因此，重磁資料處理中異常分離是不可或缺的一步，其中使用重磁數據進行異常邊界識別也是位場資料處理的重要過程，能否分辨出地下斷裂構造或密度體邊界是衡量重磁勘探方法應用效果的重要指標。隨著勘探精度的不斷提高，直接從布格重力異常或磁異常中得到地質邊界信息的難度越來越大，因此邊界增強或突出技術顯得尤為重要。邊界增強方法的基本原理是利用導數的零值點或極值點位置來判定地質異常體的邊界，在地質異常體邊界附近，重磁異常的變化率或梯度較大。邊界識別方法眾多，主要可分為數理統計和數值計算兩大類。

數理統計類重磁邊界識別方法主要包括小子域濾波和歸一化標準差兩種。楊高印[2]最早提出了小子域濾波法，該方法是在滑動平均法的基礎上進行改進的，其基本原理是將滑動平均的窗口劃分成位于中心的不同側面的8個小子域，然后將均方差最小的子域的平均值作為窗口中心的值。段曉旭[3]在小子域濾波的基礎上提出了滑動子域濾波法，改進了拐點處的數據處理效果，但是存在的折線走樣現象并沒有得到很好的解決。另一種數理統計類方法是歸一化標準差方法，該方法通過計算一個滑動窗口內的垂向一階導數的標準差與X、Y、Z方向的標準差之和的比值作為窗口中心的值。數理統計類方法的優點是可以根據實際情況的不同，選擇不同的窗口大小來消除噪音干擾，提高邊界識別的分辨率。

相對于數理統計方法，數值計算類方法更多地被用于重磁資料的處理中。該類方法主要包括垂向導數法、總水平方向導數法、解析信號法、傾斜角法和Theta圖法等[4-9]。數值計算類方法的邊界識別原理主要是通過對重磁異常進行求取垂向導數、水平導數等處理之后再根據結果的極大值位置或者零值位置來識別出地質異常體的邊界。此外，數值計算類方法還衍生出了一些處理方法，主要包括解析信號傾斜角法、廣義導數算子法、增強型均衡濾波器法、平面全張量梯度角度法和改進二階傾斜角法等[10-14]。

上述邊界識別方法的適用條件不同，處理效果也不同，了解其適用條件及如何選擇哪種方法或哪些方法更加合理至關重要。目前，針對不同處理方法適用條件及其優缺點的詳細探討和對比的文獻較少。因此，筆者使用自行編制的重磁處理軟件對上述邊界識別方法的優缺點進行了詳細的研究和對比，對于幫助提高重磁勘探的應用效果具有重要的理論和實用價值。

1 邊界識別方法

1)一階導數類邊界識別方法。主要包括垂向一階導數法、總水平方向導數法、解析信號法、傾斜角法、Theta圖法、歸一化標準差法等幾種處理方法，也是重磁勘探中最基本和最常用的方法。其中，垂向一階導數法、總水平方向梯度法以及解析信號法是對原函數直接求取一階導數，而傾斜角法和Theta圖法則是在前3種方法的基礎上經過求取比值而獲得的[4-10]。

2)二階導數類邊界識別方法。主要包括了解析信號傾斜角法、廣義導數算子法、增強型均衡濾波器法、平面全張量梯度角度法、改進的二階傾斜角法等[11-15]。解析信號傾斜法與傾斜角法利用零值位置識別邊界不同的是，該方法利用極大值位置來識別地質異常體的邊界位置。

2 軟件編制

根據重磁正演理論方法，使用C#語言編寫了三維可視化重磁正演程序。該程序能夠計算球體和直立長方體的重磁正演異常，也可以進行重磁數據的初步處理。在程序啟動界面可以選擇計算重力異常或磁異常，然后設置測網大小和點距。

為了驗證本軟件的正確性，設計并計算了一個球體和一個長方體模型的磁異常響應。球體的半徑是200 m，球心坐標X是300 m，Y是0 m，Z是500 m，磁化強度1 A/m，磁化傾角45°，磁化偏角40°。長方體的長為500 m，寬為300 m，高為200 m，中心坐標X是-300 m，Y是200 m，Z是500 m，磁化強度是1 A/m，磁化傾角45°，磁化偏角20°。測量區域網格大小X方向坐標是-1 000～1 000 m，Y方向的坐標是-1 000～1 000 m，網格點數X方向是201個，Y方向是201個，計算結果如圖1(a)所示。為了進一步驗證本軟件的正確性，使用RGIS商業軟件的重磁模型庫模塊計算了相同模型的磁異常正演，并對比了兩者的計算結果，誤差如圖1(b)所示。可以看出，兩個軟件的計算結果非常接近，誤差極小，表明了本軟件的正確性。

圖1 作者編制軟件的磁異常正演計算結果及其與RGIS商業軟件計算結果的誤差

3 不同地質異常體邊界識別方法的計算與對比

為了檢驗上述不同地質異常體邊界識別方法在不同地質條件下對邊界增強效果的優劣及不同特征，設置了3種理論重磁模型(地下單一地質異常體模型，地下兩個規模相同但埋深不同的地質異常體模型，添加隨機噪音的單一地質異常體模型)，分別使用不同地質異常體邊界識別方法對得到的重磁響應進行處理和對比。

3.1 地下單一地質異常體模型

在均勻半空間中設置一個長、寬、高分別為800 m、800 m和400 m的三維地質異常體模型，頂面埋深100 m、剩余密度為300 kg/m3；磁化傾角45°，磁化偏角0°，測量區域網格大小X方向坐標為-1 000～1 000 m，Y方向坐標為-1 000～1 000 m，網格點數X方向為201個，Y方向為201個(見圖2(a))。通過計算其重磁響應以及使用不同地質異常體邊界識別方法處理后的效果對比不同方法的邊界識別能力的優劣。

圖2 均勻半空間中地下單一地質異常體模型產生的重力異常及使用不同地質異常體識別方法處理后的結果圖

圖2是均勻半空間中地下單一地質異常體模型產生的重力異常及使用不同地質異常體識別方法處理后的結果圖。由圖2可知，垂向一階導數法利用其零值位置范圍可以確定出模型的邊界位置，只是比模型的范圍略大一點(見圖2(b))；解析信號法的極大值位置集中在模型的中部，基本能夠反映模型的位置，但是邊界識別不清晰(見圖2(c))；總水平方向導數法的極大值位置與模型邊界位置基本重合，基本可以識別出地質異常體邊界，但在角落處誤差較大(見圖2(d))；傾斜角法和垂向一階導數法一樣都利用零值位置識別邊界，其效果與垂向一階導數法類似，但是識別出的邊界有偏差，邊界識別不清晰(見圖2(e))；Theta圖法利用極大值位置識別邊界，其極大值的位置和垂向一階導數法零值位置重合，但是極大值的峰值較寬，對邊界識別的分辨率不足(見圖2(f))；歸一化標準差法與總水平方向導數法結果相似，其極大值與模型邊界基本重合，對邊界識別的分辨率較好(見圖3(g))；解析信號傾斜角法的識別結果與解析信號法基本一致(見圖3(h))；廣義導數算子法計算時使用的仰角是90°此時方位角未發揮作用，利用其零值位置識別的邊界效果和垂向一階導數法的結果相同(見圖2(i))；增強型傾斜角法的極大值位置與模型邊界重合，且極大值邊界尖銳，有利于識別邊界(見圖2(j))；平面全張量傾斜角法的極大值位置較模型邊界更大一些，識別結果沒有總水平方向導數法好(見圖2(k))；改進二階傾斜角法的計算結果非常銳利清晰，能夠最佳地突出地質異常體的邊界位置(見圖2(l))。

圖3 均勻半空間中單一地質異常體模型產生的磁異常及使用不同地質異常體識別方法處理后的結果圖

通過圖2不同地質異常體識別方法的處理效果可以看出，對目標體邊界位置識別最準確和最清晰的是改進二階傾斜角法和增強型傾斜角法；其次是總水平方向導數法與歸一化標準差法，二者基本能夠識別出模型的邊界位置；平面全張量傾斜角法相比總水平方向導數法的邊界位置略微向外偏移，對邊界識別效果不佳；垂向一階導數法、傾斜角法、Theta圖法和廣義導數算子法的邊界識別結果相似，雖然基本能夠識別模型的邊界位置，但是得到的邊界位置都與模型邊界有偏差，識別效果一般；解析信號法和解析信號傾斜角法對邊界的識別效果最差。

圖3是均勻半空間中單一地質異常體模型產生的磁異常及使用不同地質異常體識別方法處理后的結果圖。為了考察傾斜磁化的影響，這里均不做化極處理。由圖3可知，垂向一階導數的零值位置與模型邊界位置相差較大，不利于邊界識別(見圖3(b))；解析信號法的極大值與模型邊界位置基本重合，能夠較為準確地指示邊界位置，但是不夠完整，只能識別與磁化方向垂直的邊界(見圖3(c))；總水平導數法的極大值存在雙峰現象，也能在一定程度上識別出模型邊界位置(見圖3(d))；傾斜角法的零值位置和垂向一階導數法的零值位置重合，同樣與模型邊界位置相差較大(見圖3(e))；Theta圖法的極大值位置與垂向一階導數法的零值位置基本相同，識別效果較好，但在邊界識別中具有一定誤差(見圖3(f))；歸一化標準差法的極大值位置識別出的邊界相對完整，與磁化方向平行的邊界位置能夠準確地識別出來，但與磁化方向垂直的邊界識別會出現雙峰現象(見圖3(g))；解析信號傾斜角法的極大值位置與模型邊界位置基本重合，能夠準確且較完整地識別出邊界(見圖3(h))；廣義導數算子法的識別結果與垂向一階導數基本相同，無法準確指示出模型邊界位置(見圖3(i))；增強型傾斜角法的邊界識別結果比較準確，但是存在較多的假邊界(見圖3(j))；平面全張量傾斜角法也能準確的指示出部分邊界的位置，但是存在較多的假邊界(見圖3(k))；改進二階傾斜角法的極大值位置清晰，邊界識別結果也較準確，但是同樣也存在多余的假邊界現象(見圖3(l))。

綜上所述，使用磁異常數據進行地下地質異常體邊界識別時，解析信號傾斜角法的邊界識別效果最好；其次是解析信號法和Theta圖法，受到傾斜磁化的影響也較小；總水平導數法的識別結果較好，能夠顯示邊界位置，但是由于傾斜磁化的影響，對邊界的位置識別得不太準確，且會出現雙峰現象；歸一化標準差法、增強型傾斜角法和改進二階傾斜角法都可以識別出邊界，但是在正常場區域會存在較多的假邊界；其他方法受到斜磁化的影響較大，邊界識別的誤差也較大，不太適合直接用于邊界的識別。

3.2 地下兩個規模相同但埋深不同的地質異常體模型

為了考察上述不同方法對相鄰多個規模相同但埋深不同地質異常體模型的分辨能力，在均勻半空間中設置了兩個規模相同(長寬高規模均分別為500 m、500 m和400 m)、剩余密度相同(300 kg/m3)、但埋深不同(左右兩個地質異常體的頂面埋深分別為100 m和300 m)的地質異常體模型。磁化強度1 A/m，磁化傾角45°，磁化偏角0°。左邊長方體中心位置為坐標X=-400 m，Y=0 m，Z=300 m；右邊長方體中心位置為坐標X=400 m，Y=0 m，Z=500 m。

圖4是在均勻半空間中放置了兩個埋深與規模相同的地質異常體產生的重力異常及使用不同地質異常體識別方法處理的結果圖，可以看出，不同的處理方法對于多個地質異常體邊界的刻畫能力與單個地質異常體情況很相似，對邊界位置刻畫最準確、分辨率最高的仍是改進二階傾斜角法和增強型傾斜角法；總水平方向導數法、歸一化標準差法的極值均與模型邊界基本重合，識別效果較好，但總水平方向導數法受到相鄰模型的影響較大，導致相鄰兩個模型之間的邊界突出不夠；平面全張量傾斜角法受相鄰模型的影響也較大，甚至在二者之間產生了虛假異常；垂向一階導數法對邊界識別能力較好，只是邊界略微偏大一點；傾斜角法、Theta圖法和廣義導數算子法的邊界識別結果與模型邊界均有偏差，分辨精度不足；解析信號法和解析信號傾斜角法的極值大致能與兩個模型位置相對應，并指示出模型的中心位置，但是清晰展現邊界的能力稍顯不足。

圖4 均勻半空間中兩個埋深不同地質異常體模型產生的重力異常及使用不同地質異常體識別方法處理后的結果圖

圖5是兩個埋深不同的長方體模型產生磁異常圖和不同地質異常體識別方法的處理結果圖。同樣為了考察傾斜磁化的影響，所有圖件也均不做化極處理。由圖5可知，垂向一階導數法的處理結果由于受到傾斜磁化和不同埋深的地質異常體的影響，其零值位置對比復雜，埋深較大的地質異常體反映不明顯，不能準確指示地質異常體邊界位置(見圖5(b))；解析信號法能夠準確地識別淺部地質異常體的部分邊界位置，但是卻無法顯示深部地質異常體的邊界位置(見圖5(c))；總水平導數法同樣受到不同埋深的影響，無法有效地識別埋深較大地質異常體的邊界(見圖5(d))；傾斜角法受埋深的影響較小，能夠平衡不同埋深的異常，但是其零值位置與垂向一階導數的零值位置一樣，存在較大誤差，對識別邊界位置不利(見圖5(e))；Theta圖法也能夠平衡不同埋深的異常，其邊界識別結果同樣受到傾斜磁化的影響，存在較大誤差(見圖5(f))；歸一化標準差法受地質異常體埋深不同的影響較小，能夠平衡不同埋深的異常值，淺部地質異常體邊界識別結果較準確，深部的識別結果存在一定誤差(見圖5(g))；解析信號傾斜角法不易受到埋深的影響，邊界識別結果較準確，出現的假邊界也較少(見圖5(h))；廣義導數算子法能夠平衡不同埋深的異常，但是其零值位置與地質異常體邊界位置相差較大(見圖5(i))；增強型傾斜角法的識別結果對比準確，深部地質異常體的邊界識別存在一定誤差，而且存在假邊界(見圖5(j))；平面全張量傾斜角法能夠平衡不同埋深的異常，埋深較淺的地質異常體邊界識別結果較好，深部識別誤差較大(見圖5(k))；改進二階傾斜角法的結果與增強型傾斜角類似，但是對邊界識別的結果更加清晰(見圖5(l))。

圖5 均勻半空間中兩個埋深不同地質異常體模型產生的磁異常及使用不同地質異常體識別方法處理后的結果圖

綜上所述，解析信號傾斜角法、歸一化標準差法和傾斜角法等使用了比值原理進行歸一化的方法受埋深影響較小，能夠較好地同時識別出不同埋深的地質異常體邊界；沒有進行歸一化的方法由于受埋深影響，對深部地質異常體邊界的識別效果較差。

3.3 添加隨機噪音的單一地質異常體模型

實測的重磁異常數據往往包含著各種噪音，在進行識別之前要經過濾波處理以消除噪音干擾。但是濾波等處理通常僅能夠消除部分噪音，剩余的噪音往往會影響邊界識別的效果。為了測試上述地質異常體識別方法在存在噪音情況下對邊界的分辨能力，設置了與3.1完全相同的地質異常體模型，在正演的重磁異常值中添加了5%的隨機噪音，以此研究上述方法在未處理干凈所有噪音情況下的邊界識別效果。

圖6是均勻半空間中單個地質異常體產生的重力響應添加5%隨機噪音后再濾波的結果及不同地質異常體邊界識別方法處理后的結果圖，可以看出，垂向一階導數法、解析信號法和總水平方向導數法受噪音影響最小，非異常區域產生的假異常幅值與異常區域幅值差異明顯，能夠分辨出地質異常體的邊界；傾斜角法、廣義導數算子法和平面全張量傾斜角法受噪音影響相對較小，可以分辨出原始地質異常體的位置和形狀；Theta圖法容易受到噪音的影響，產生了較多虛假異常，雖可大致看出地質異常體邊界，但邊界不夠清晰；歸一化標準差法通過選取窗口計算標準差，雖在一定程度上放大了噪音的干擾，產生了較多虛假異常，但對邊界識別的效果尚可；解析信號傾斜角法、增強型傾斜角法和改進二階傾斜角法在計算中均使用了高階導數，可以提高邊界識別的分辨率，但同時也變得對噪音更加敏感，產生了較多的虛假邊界，對識別出邊界的形狀反而不利。

圖6 均勻半空間中單個地質異常體產生的重力響應添加5%隨機噪音后再濾波的結果及不同地質異常體邊界識別方法處理后的結果圖

圖7是均勻半空間中單個磁異常模型正演數據添加5%噪音并濾波后的結果及不同地質異常體邊界識別方法處理后的結果圖。同樣地均未做化極處理。由圖7可知，垂向一階導數法、解析信號法和總水平導數法受噪音影響較小，處理結果與未添加噪音的結果相差不大，極值位置的幅值與其他區域的幅值相差明顯，能夠較清晰地識別邊界的位置；傾斜角法、Theta圖法、歸一化標準差法、解析信號傾斜角法、廣義導數算子法和平面全張量傾斜角法均對噪音較敏感，邊界位置出現了變形，但也可以基本識別出邊界；增強型傾斜角法和改進二階傾斜角法受噪音影響更嚴重，非邊界位置產生了許多虛假干擾，無法對邊界進行準確識別。

圖7 均勻半空間中單個磁異常模型正演數據添加5%噪音并濾波后的結果及不同地質異常體邊界識別方法處理后的結果圖

綜上所述，在存在噪音的情況下，總水平方向導數法和平面全張量傾斜角法受到的影響較小，能夠較好地識別出地質異常體的邊界；垂向一階導數法、傾斜角法和廣義導數算子法對噪音不太敏感，抗噪能力較強，雖然其零值位置比邊界略大，但也基本能夠標識出地質異常體邊界；其他方法在存在明顯噪音情況下則難以準確識別出地質異常體的邊界。

此外，根據上述研究還可以看出，在重磁異常邊界識別中，總水平導數法雖然定位邊界的位置對比準確，但卻不利于識別埋深較大的地質異常體的邊界；歸一化標準差法和增強型傾斜角法都能夠平衡不同埋深的異常，識別出更多的邊界，其中增強型傾斜角的識別結果更加清晰。在磁異常邊界的識別中，總水平導數法和解析信號法不利于識別埋深較大的地質異常體的邊界；歸一化標準差法能夠識別出更多的邊界位置，解析信號傾斜角法的效果最好，不僅能夠定位更多的邊界位置，而且結果也更加清晰。

4 結束語

本文通過理論模型計算，對重磁勘探中多種地質異常體邊界識別方法進行了詳細的分析和對比。通過本文的對比研究發現，垂向一階導數法、傾斜角法和廣義導數算子法的極小值位置基本相同，與Theta圖法的極大值位置也基本相同。因此上述基于相似數學原理的處理方法對邊界的識別效果也相似，在實際應用中只需要選取一種方法即可。其中廣義導數算子法相較于其他方法，既能夠平衡不同深度的異常，而且受噪音的影響也更小，相對來說處理效果更好；解析信號傾斜角法、增強型傾斜角法和改進二階傾斜角法使用了高階導數的方法，雖然分辨率更高，但同時對噪音也更加敏感，使用之前應該首先對原始數據進行充分的濾波等去噪處理；總水平方向導數法的邊界識別結果較準確，受噪音影響也較小，但其也存在無法平衡不同埋深異常的缺點。通過本文研究，總體上認為，總水平方向導數法、增強型傾斜角法和歸一化標準差法相較于其他地質異常體邊界識別方法效果更好。