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

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

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

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

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

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

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

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

1 邊界識別方法

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

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

2 軟件編制

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

為了驗證本軟件的正確性，設(shè)計并計算了一個球體和一個長方體模型的磁異常響應(yīng)。球體的半徑是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°。測量區(qū)域網(wǎng)格大小X方向坐標是-1 000～1 000 m，Y方向的坐標是-1 000～1 000 m，網(wǎng)格點數(shù)X方向是201個，Y方向是201個，計算結(jié)果如圖1(a)所示。為了進一步驗證本軟件的正確性，使用RGIS商業(yè)軟件的重磁模型庫模塊計算了相同模型的磁異常正演，并對比了兩者的計算結(jié)果，誤差如圖1(b)所示。可以看出，兩個軟件的計算結(jié)果非常接近，誤差極小，表明了本軟件的正確性。

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

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

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

3.1 地下單一地質(zhì)異常體模型

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

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

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

圖3 均勻半空間中單一地質(zhì)異常體模型產(chǎn)生的磁異常及使用不同地質(zhì)異常體識別方法處理后的結(jié)果圖

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

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

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

3.2 地下兩個規(guī)模相同但埋深不同的地質(zhì)異常體模型

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

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

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

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

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

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

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

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

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

圖6 均勻半空間中單個地質(zhì)異常體產(chǎn)生的重力響應(yīng)添加5%隨機噪音后再濾波的結(jié)果及不同地質(zhì)異常體邊界識別方法處理后的結(jié)果圖

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

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

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

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

4 結(jié)束語

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