基于重力異常增強的改進小子域濾波方法

王彥國 黃遠生 張 瑾

(東華理工大學放射性地質與勘探技術國防重點學科實驗室，江西南昌 330013)

0 引言

運用位場數(shù)據(jù)識別地質體邊界是對地球物理資料進行地質解釋的一項重要內(nèi)容。近年來，在利用位場資料推算斷裂構造的水平位置及深度方面已出現(xiàn)了許多方法。其中基于位場梯度的邊界識別方法算法[1]簡單、易于實現(xiàn)，是一類常用的方法(如水平總梯度[2-3]、 垂向導數(shù)[4-5]、 Theta map[6]、 Tilt angle[7]、歸一化標準差[8]等)。但是這類方法對區(qū)域邊界信息不敏感，且易受隨機干擾的影響。Euler反褶積法[9-11]也可以用于邊界識別，該方法是通過計算位場一階導數(shù)和點位坐標實現(xiàn)的，導數(shù)計算結果或計算點位稍有偏差，便會在反演中產(chǎn)生較大的偏離，從而極大地影響邊界反演精度。利用重力歸一化總梯度及其相位同樣可以確定場源深度和邊界水平位置[12],但此法受諧波數(shù)、向下延拓的影響，場源深度和水平位置的估計受到了一定限制，因此精度不高。

小子域濾波法是楊高印[13]提出的一項位場數(shù)據(jù)處理技術，是對滑動平均法的一種改進，該方法能夠增強區(qū)域異常的界限特征。夏玲燕等[14]將小子域濾波和Tilt梯度相結合，許海紅等[15]將小子域濾波和水平總梯度相結合，均獲得了較良好的邊界識別效果。馬濤等[16]針對小子域濾波存在偏心問題對梯度保持濾波法進行改進，使濾波結果更加穩(wěn)定可靠； 肖鋒等[17]采用“田”字型子域劃分模式，使計算更便捷； 馬國慶等[18]針對小子域濾波存在異常曲線扭曲的缺陷，提出了優(yōu)化小子域濾波，更好地展示了區(qū)域異常間的界限。段曉旭[19]采用加密點距的方式解決小子域存在異常曲線扭曲的問題，并提出點位移動法增強異常梯度帶。蔡鐘等[20]針對小子域濾波存在邊界異常失真的現(xiàn)象，提出了五/六邊形子域濾波法，在一定程度上改善了小子域濾波的效果。許海紅等[21]對上述多種小子域濾波方法進行了對比分析，給出了不同方法的優(yōu)缺點及適用條件。但是上述改進方法僅能對具有明顯梯級帶特征的斷裂構造進行識別，而無法檢測具有弱梯級帶或非梯級帶特征的場源邊界。為此，張鳳旭等[22-23]提出了三方向小子域濾波進行斷裂構造識別，提高了邊界識別能力； Jiang等[24]對小子域濾波的剖分模式進行了改進，并實現(xiàn)了基于梯度張量的小子域濾波，同樣獲得了豐富的邊界信息； Tai等[25]在三方向小子域濾波基礎上，提出了加權小子域濾波法，使斷裂構造的劃分更快捷、更精確。但這三種改進方法均采用了差分或導數(shù)計算，因此計算更易受到干擾的影響。

針對上述方法的優(yōu)缺點，本文提出了新的子域剖分模式，使其可以更準確地檢測出不同方向的異常界限; 同時為了消除干擾對小子域濾波的影響和提高異常分辨率，給出了一種新的濾波技術——基于迭代差分的穩(wěn)定增強濾波。

1 改進小子域濾波

1.1 傳統(tǒng)小子域濾波

傳統(tǒng)小子域濾波[13]是將滑動窗口按中心點的不同側面劃分出8個子域，并將8個子域均方差最小子域的平均值作為濾波輸出，圖1為傳統(tǒng)小子域濾波子域剖分模式的示意圖。該方法可以較好地保留異常間的界限，但還存在著以下不足： ①雖然對位場異常進行了低通濾波處理，然而不同子域對噪聲的敏感程度不同，從而導致濾波結果中出現(xiàn)明顯的“虛假”異常[26]； ②8個子域的重心均偏離窗口中心點，易造成異常形態(tài)不規(guī)則扭曲，從而產(chǎn)生虛假信息[15,18]； ③不能較好地檢測出梯級帶不明顯或非梯級帶特征的異常界限，且檢測結果存在較大偏差[22]。

1.2 改進型小子域濾波

針對傳統(tǒng)小子域濾波存在的問題，本文對子域剖分模式進行改進，改進型小子域濾波與傳統(tǒng)方法的區(qū)別在于剖分子域的個數(shù)不同。以5×5窗口為例，傳統(tǒng)剖分是按中心點的不同側面劃分出8個子域(圖1)，而本文方法的剖分是將5×5窗口剖分出21個子域，每個子域包含9個數(shù)據(jù)點(圖2)。從圖2中可以看出，改進后的子域劃分方式可以檢測任意方向上的異常變化，能更好地體現(xiàn)出不同走向上的構造，獲得更為準確、豐富的異常界限。該方法的計算流程與傳統(tǒng)方法一致，即將窗口內(nèi)21個子域均方差最小的子域數(shù)據(jù)平均值作為窗口中心點的濾波結果，以緊縮梯級帶作為場源邊界的識別依據(jù)。

圖1 傳統(tǒng)小子域濾波5×5滑動窗口及分解成的8個子域[13]

圖2 改進型小子域濾波5×5滑動窗口及分解成的21個子域

2 異常增強濾波技術

傳統(tǒng)濾噪技術一般是低通濾波。低通濾波在消除噪聲時，對低頻有效信號同樣也有一定的消減，從而導致異常分辨率降低。因此，在原異常中表現(xiàn)為非梯度帶的異常界限，小子域濾波仍不能有效識別。鑒于此，本文提出了一種穩(wěn)定增強濾波技術，該方法在濾除噪聲的同時，可以提高對異常的分辨能力。

2.1 基本原理

設位場異常u的波譜為U，低通濾波算子為φ，0≤φ≤1，濾波后的波譜可以表示為

Uφ=Uφ

(1)

雖然上述低通濾波可以濾除隨機干擾，但對低頻有效信息也存在一定的消減作用。因此，采用迭代補償模式完成對低頻成分的彌補，迭代過程如下：

(1)將U與Uφ的差值利用濾波算子φ進行修正，修正結果為

(2)

=Uφ[1+(1-φ)+(1-φ)2]

(3)

(3)按照如上迭代補償模式，第n次修正結果為

(4)

根據(jù)式(2)～式(4)，得到迭代通式

(1-φ)n]=U[1-(1-φ)n+1]

(5)

為了提高異常分辨率，對迭代結果進行如下處理：

=F-1{U(1-φ)n+1[1-(1-φ)Δn]}

(6)

式中F-1表示反傅里葉變換。

=F-1{U(1-φ)n+1[1-(1-φ)Δn]2}

(7)

(3)仿照上述流程，第m階迭代差分結果為

(8)

由式(8)可知，迭代差分m次的濾波因子Φ=(1-φ)n+1[1-(1-φ)Δn]m事實上是一個低通濾波ΦL=[1-(1-φ)Δn]m與高通濾波ΦH=(1-φ)n+1的組合。

2.2 濾波特性

令濾波算子φ=exp(-wh)，m=1， Δn=10，其中w表示波數(shù)，h=10Δx， Δx=0.1km。圖3a給出了上文提及的濾波算子的濾波曲線?？梢钥闯?，低通濾波φ(向上延拓算子)雖然對高頻成分具有較好的壓制作用，但對中、低頻的有用信號一定程度上也有消減作用；相對于低通濾波φ，增強濾波器中的低通濾波ΦL可以更好地保留低頻成分；穩(wěn)定增強濾波Φ在中、低頻則趨近于高通濾波算子ΦH，而在高頻處則趨于低通濾波ΦL，即濾波算子Φ為一帶通濾波器，對中頻成分起到了放大作用，而對高頻成分起到了壓制作用。也就是說，穩(wěn)定增強濾波后的位場異常在理論上不僅具有較強的穩(wěn)定性，而且還可以提高原異常的分辨能力，因此，將這種迭代差分濾波稱為穩(wěn)定增強濾波。

圖3b～圖3d是不同參數(shù)的增強濾波響應曲線?？梢钥闯?，隨著初始迭代次數(shù)n的增加(圖3b)，濾波算子的主峰向高頻移動，但幅值卻隨之降低；隨著迭代次數(shù)間隔Δn的增加(圖3c)，增強濾波算子頻帶變寬，主峰同時向高頻移動；隨著迭代差分階數(shù)m的增加(圖3d)，濾波算子頻帶變窄且主峰向低頻移動。因此選擇較少的初始迭代次數(shù)n、較小的迭代次數(shù)間隔Δn和較大的迭代差分次數(shù)m，會使濾波結果更穩(wěn)定。

圖3 穩(wěn)定增強濾波算子濾波響應曲線

理論上，穩(wěn)定增強濾波器中的低通濾波ΦL、迭代初始次數(shù)n、迭代間隔Δn及迭代差分階數(shù)m的選擇具有隨意性，但為確保計算結果的穩(wěn)定性及異常分辨能力，建議φ選用向上延拓不少于5倍點距的延拓算子，選取n=1、Δn<50及m< 5。后文中的模型及實例部分低通濾波算子均是選取向上延拓10倍的延拓算子，m=n=1，Δn=10。

3 理論模型分析

為了驗證本文方法的有效性和優(yōu)越性，建立了由4個埋深不同、尺度不同、密度不同的重力異常體(表1)組成的復雜模型進行試驗。同時，在模型的正演重力異常中添加了隨機干擾，以檢驗方法的穩(wěn)定性。圖4是該模型產(chǎn)生的重力異常和分別添加1%、 3%噪聲的重力異常及利用本文方法對圖4c進行穩(wěn)定增強濾波處理的結果。可以看出，理論重力異常(圖4a)可以通過梯度帶大致識別異常體A、B、C的邊界位置，但梯度帶較寬，邊界位置不明確；模型體D的異常受疊加場影響，圈閉極值與模型體中心不對應，在異常圖主要以等值線同向扭曲為主。在重力異常含1%隨機噪聲時(圖4b)，異常體D所對應的異常更加模糊。當重力異常含3%噪聲時(圖4c)，難以通過異常梯度帶識別所有模型邊界。穩(wěn)定增強濾波后的重力異常(圖4d)不僅等值線相對于圖4c更加圓滑，且異常分辨率也得到明顯提高，尤其對規(guī)模小、埋深大的異常體D。

表1 理論模型異常體的參數(shù)

圖4 疊加模型產(chǎn)生的理論重力異常(a)、含1%噪聲的重力異常(b)、含3%噪聲的重力異常(c)及其穩(wěn)定增強濾波結果(d)

圖5是無噪數(shù)據(jù)在不同窗口下的傳統(tǒng)小子域濾波與本文改進方法濾波處理結果?？梢钥闯觯斢嬎愦翱谳^小時(如3×3)，傳統(tǒng)小子域(圖5a上)與改進小子域濾波(圖5a下)僅對埋深較淺的異常體C的邊界異常進行了有效緊縮；隨著窗口尺度的增加，無論是傳統(tǒng)小子域濾波還是改進小子域濾波，地質體邊界異常被緊縮得更加明顯，如9×9窗口時，除異常體D缺失部分邊界，其他三個異常體的邊界異常均得到了有效增強。由圖5可以明顯發(fā)現(xiàn)，傳統(tǒng)小子域濾波結果在模型體角點處存在明顯的等值線扭曲現(xiàn)象，而改進型小子域濾波的處理結果并不存在此現(xiàn)象，能更好地識別異常體的邊界。

圖6是不同噪聲水平的重力異常在不同窗口下的傳統(tǒng)小子域濾波與改進小子域濾波處理結果對比?？梢钥闯觯捎陔S機噪聲對不同方向子域的影響程度不同，導致了兩種方法的濾波結果中均存在著多條不連續(xù)分布的異常緊縮帶，且無法確定哪些緊縮帶是有效的；而且，噪聲水平越高，小子域濾波處理結果中的虛假信息就越多，處理結果的可靠性越差。

圖7是含3%噪聲重力異常經(jīng)過穩(wěn)定增強濾波后的傳統(tǒng)小子域濾波和改進型小子域濾波處理結果對比。相對于圖6來說，穩(wěn)定增強濾波后的小子域濾波更加穩(wěn)定，不但不包含明顯的虛假緊縮帶，更能有效地檢測出地質體邊界，尤其異常體D的邊界也有明顯的展示。

圖5 不同大小窗口下的傳統(tǒng)小子域濾波(上)和改進小子濾波(下)結果對比

圖6 含1%及3%噪聲時不同窗口下的傳統(tǒng)小子域濾波(上)與本文濾波(下)結果對比

圖7 含3%噪聲重力異常經(jīng)穩(wěn)定增強濾波后的傳統(tǒng)小子域濾波和改進小子域濾波結果對比

4 實際重力資料應用

為了驗證穩(wěn)定增強異常改進型小子域濾波對實際資料的處理效果，選取吉林省南部鴨綠江盆地實測重力數(shù)據(jù)(點距和線距均為1km)進行試驗。鴨綠江盆地位于中朝板塊東北緣，二級大地構造單元隸屬于遼東臺隆區(qū)，盆地主體為太子河—渾江坳陷，盆地內(nèi)除志留系、泥盆系地層缺失外，其他時代的巖石均有出露(圖8)。研究區(qū)存在一系列逆沖推覆構造，使得地層關系十分復雜、連續(xù)性極差及地層倒置現(xiàn)象非常普遍[28]，這給地質體邊界的確定甚至地質解釋都帶來了較大難度。研究區(qū)內(nèi)各個時代的巖石出露較全，奧陶系與石炭系間存在明顯的密度差(表2)，這為方法試驗的可靠性提供了較好的物性基礎。

從盆地內(nèi)布格重力異常圖(圖9a)可以看出，反映地質體邊界位置或斷裂構造的梯級帶異常梯度較平緩，因而地質體邊界不易直接根據(jù)重力異常進行精確解釋。較小型地質體的邊界受區(qū)域異常的影響較大，在異常圖中主要表現(xiàn)為異常等值線突然變寬或變窄以及同形扭曲等非梯級帶特征，地質體的邊界位置確定難度較大。從圖9b可以看出，經(jīng)穩(wěn)定增強處理的重力異常圖的分辨率獲得了明顯提高，不僅異常梯級帶連續(xù)性增強，且異常呈現(xiàn)更明顯的NE走向，與地質圖(圖8)中的巖石分布呈現(xiàn)出了很強的一致性。圖10 是布格重力異常的傳統(tǒng)小子域和改進型小子域濾波結果，可以看出，對于原始布格重力異常，傳統(tǒng)小子域濾波(圖10上)對大型異常梯度帶進行了較好的緊縮，然而緊縮帶的連續(xù)性較差，且存在明顯的無規(guī)律性彎曲；而改進型小子域濾波(圖10下)同樣也反映了大型斷裂的位置，且邊界異常信息更豐富，緊縮帶的連續(xù)性也更強。對于穩(wěn)定增強異常，無論是經(jīng)傳統(tǒng)小子域濾波(圖11上)還是改進型小子域濾波(圖11下)處理的重力數(shù)據(jù)對邊界的識別能力顯然都獲得了大幅度的提升，邊界異常信息更豐富，而改進型小子域濾波結果顯示出更豐富的構造邊界信息，且緊縮帶的走向和連續(xù)性也更好。需要指出的是，5×5窗口的小子域濾波(圖11a)數(shù)據(jù)中的一些緊縮帶在9×9窗口的小子域濾波中力圖(圖11b)上沒有顯示。如在大安鎮(zhèn)的東南側，5×5窗口的小子域濾波數(shù)據(jù)中存在兩條NE走向的異常緊縮帶，但9×9窗口的小子域濾波數(shù)據(jù)中并未被有效識別出來，且檢測出的緊縮帶連續(xù)性也較差。這說明在實際資料處理中，需要結合不同窗口的小子域濾波數(shù)據(jù)進行地質解釋。另外，新方法獲得的異常緊縮帶與不同時代的巖性界限(圖8)具有良好的對應性；增強異常的改進小子域濾波數(shù)據(jù)的負異常分布也較好地反映了研究區(qū)中—新生代地層和燕山期花崗巖等低密度巖石的分布。

圖8 鴨綠江盆地地質圖及工區(qū)位置

表2 鴨綠江盆地地層/巖石密度統(tǒng)計表 單位： g·cm-3

圖9 鴨綠江盆地布格重力異常(a)及穩(wěn)定增強重力異常(b)

圖10 鴨綠江盆地布格重力異常的傳統(tǒng)小子域(上)和改進型小子域(下)濾波結果對比

圖12a是結合區(qū)域地質及前人工作成果繪制的構造分區(qū)圖。圖中給出了大型斷裂的位置、低密度巖石的分布區(qū)(燕山期花崗巖及石炭—白堊紀地層)以及渾江煤田工作區(qū)等信息，以此驗證本文方法的可靠性。圖12b是基于圖11下解釋的3條大型斷裂和17條中小型斷裂以及在渾江坳陷內(nèi)圈定的負異常區(qū)域，其中F1為二道江—江源斷裂，是龍崗隆起與渾江坳陷的界限斷裂； F2為石人—大青溝斷裂，是渾江坳陷與老嶺隆起的分界斷裂； F3為鴨綠江斷裂，在研究區(qū)東部則是老嶺隆起與四道溝坳陷的界限。這三條大型斷裂在重力異常中均以大型梯級帶為異常特征，其余17條中、小型斷裂在重力異常圖上表現(xiàn)為弱異常梯級帶或等值線突然變寬或變窄以及同形扭曲等異常特征。需要指出的是，在渾江坳陷內(nèi)，地表出露的低密度巖石分布區(qū)僅占30%左右的面積，而本文方法識別出的負異常分布范圍可達50%以上，也就是說渾江坳陷內(nèi)有大面積的低密度體隱伏于高密度巖石之下，其中勘探程度較高的渾江煤田工作區(qū)中的一些飛來峰下的煤層(主要分布在石炭、二疊、侏羅紀地層中)已被鉆探或采掘所證實[28-29]，這也再一次證實了本文方法的實用性，同時該方法也為研究區(qū)后續(xù)煤田勘查有利區(qū)(尤其勘探程度較低的六道溝—三道湖—石人鎮(zhèn)一帶)和地質—地球物理綜合解釋的有效手段之一。

圖11 鴨綠江盆地穩(wěn)定增強重力異常的傳統(tǒng)小子域(上)和改進型小子域(下)濾波結果對比

圖12 鴨綠江盆地構造分區(qū)圖(a)和斷裂構造分布及渾江坳陷負異常分布區(qū)(b)

5 結論

針對傳統(tǒng)小子域濾波存在異常曲線扭曲的缺陷，本文提出了一種包含更多子域的新子域剖分模式，可以更準確地檢測不同走向構造的界限。針對小子域濾波存在易受干擾影響和難以識別非梯級帶特征的斷裂構造的缺陷，提出了基于迭代差分的穩(wěn)定增強濾波技術。該技術在濾除噪聲干擾的同時，還可以提高異常的分辨率。模型試驗和實例應用表明，改進型小子域濾波比傳統(tǒng)算法效果更佳，異常的界限可以得到更加準確、精細的刻畫；穩(wěn)定增強濾波后的重力異常不僅穩(wěn)定性強，且異常分辨率明顯提高；增強異常改進型小子域濾波則能夠更加準確、精細地刻畫密度異常體的界限，同時還可以通過重力正、負異常的分布進一步提高方法的解釋能力。