改進的均衡濾波器在位場數據邊界識別中的應用

馬國慶，黃大年，于 平，李麗麗

吉林大學地球探測科學與技術學院，長春 130021

1 引 言

位場主要是指重力和磁力場，重磁勘探具有快速、經濟、范圍大的優點，是進行構造劃分和異常圈定快速而有效的方法，但是原始異常與地質體邊界之間不存在良好的對應關系，所以直接根據重磁異常不易識別出地質體的邊界.重磁異常水平導數的極大值和垂直導數的零值與地質體邊界相對應，該性質常被用來進行地質體邊界的識別.垂直導數（1936）是最早被用來進行地質體邊界識別的方法［1-2］；后來人們發現水平導數的極值位于密度與磁化率發生突變的位置，該方法被證明是一種非常有效的邊界識別工具［3-4］；解析信號的極大值與磁性體的邊界也存在良好的對應關系［5］，Hsu等利用高階解析信號使磁性體的邊界更加清晰［6］，但上述邊界識別方法均存在不能同時顯示淺部和深部地質體邊界的缺點［7］，這是由于導數隨地質體埋深的衰減速度較快，因此更易凸顯出淺部地質體的效應.為了解決這一問題，人們開始致力于均衡濾波器的研究.Miller和Singh在1994年提出了第一個均衡邊界識別濾波器—tilt angle［8］，該方法能很好地均衡強弱異常之間的幅值，但不是一個很好的邊界識別工具；Rajagopalan和 Milligan在1995年提出利用自動增益控制來進行磁源體邊界的識別［9］，該方法的識別結果依賴于窗口的尺寸，靈活性較差；Verduzco等在2004年提出利用tilt angle的總水平導數（THDR）進行地質體邊界的識別［10］，取得不錯的效果；Wijns等在2005年利用總水平導數與解析信號的比值（theta map）來進行邊界識別［11］，取得了很好的效果；Cooper和Cowan在2006總結了多種不同形式的邊界識別工具［12］，并提出了其它形式的傾斜角進行異常體邊界的識別；Cooper和Cowan 2008年利用水平與垂直導數的均方差進行異常體邊界的識別［13］，該方法的結果與窗口的尺寸有關；后來人們利用異常的希爾伯特變換來進行地質體邊界的識別［7，14］，該方法可有效地降低噪聲的干擾，但所識別出的地質體邊界比較模糊；馬國慶等在2011年提出利用總水平導數與垂直導數的相關系數進行地質體邊界的識別［15］，取得了不錯的應用效果，但是該方法對小范圍異常體的識別效果不是很好.Ma和Li在2012年采用歸一化總水平導數法進行地質體邊界的識別［16］，有效地提高了輸出結果的穩定性，但識別出的邊界存在一定的發散.

本文提出增強型均衡濾波器來進行地質體邊界的識別，該濾波器是由不同階導數構建的，其能更加清晰和準確地識別出地質體的邊界.

2 增強型均衡濾波器

現有的邊界識別濾波器均是同階導數之間的組合，增強型濾波器是利用不同階導數之間的非線性組合，其輸出結果更加準確和清晰.本文給出了三種不同形式的增強型均衡濾波器.

2.1 增強型tilt angle

Miller和Singh在1994年提出采用tilt angle進行地質體邊界的識別，其表達式為：

其中，f為原始重力或磁異常.tilt angle能很好地均衡不同深度異常之間的幅度，但是該方法并不能很好地進行地質體邊界的識別［15］，因此Verduzco等在2004年提出利用tilt angle的總水平導數（THDR）進行地質體的邊界識別工作.

本文提出了另一種形式的tilt angle，稱之為增強型tilt angle（ETA），其表達式為：

其中，n表示垂直導數的階數，也以此表示濾波器的階數，

系數cn是為了使公式（3）滿足數學運算法則，并可均衡不同階導數之間的強度.增強型tilt angle為垂直導數與其水平導數的組合，其分子與分母是不同階的，這是其與常規tilt angle的主要區別，增強型tilt angle的最大值與地質體的邊界相對應.

2.2 增強型theta map

Wijns等在2005年提出theta map來進行地質體邊界的識別，其具體公式為：

增強型theta map（ETM）的具體表達式為：

增強型theta map的最大值也與地質體邊界相對應.

2.3 增強型cosine function

第三種是利用不同階導數之間比值的余弦函數來進行異常體邊界的識別，稱之為增強型cosine function（ECF），其具體表達式為：

其中，R表示取實部運算，這是因為反余弦函數的運算結果中會出現虛數.增強型cosine function的最大值與地質體邊界相對應.為了顯示增強型濾波器的可行性，給出一系列重力異常剖面（圖1）.

圖1a為模型的原始異常，2個模型埋深分別為12m和15m.圖1b表示原始異常的THDR計算結果，該方法不能清晰地給出深部地質體的邊界，且會產生多余的干擾異常；圖1c表示異常的theta map計算結果，該方法能同時顯示淺部與深部異常的邊界，但是所識別出的邊界較發散；圖1d—1f分別表示ECF1、ETA1和ETM1的計算結果；圖1g—1i分別表示ECF2、ETA2和ETM2的計算結果；一階和二階增強型均衡濾波器均能很清晰地識別出地質體的邊界，邊界非常集中.圖1j—1l分別表示ECF3、ETA3和ETM3的計算結果，三階增強型均衡濾波器的識別結果產生了多余的峰值，不利于識別出地質體的邊界，無法對異常進行解釋.因此僅可以利用一階和二階增強型均衡濾波器來進行地質體邊界的識別.

從增強型均衡濾波器的表達式中可以看出，二階增強型濾波器需要計算二階垂直導數，二階垂直導數的計算會明顯增大噪聲的干擾.為了減小噪聲的干擾，采用Laplace方程［17］來完成二階垂直導數的計算：

圖1 不同方法邊界識別結果的剖面圖（a）原始重力異常；（b）異常的THDR結果；（c）異常的theta map結果；（d）異常的ETA1 結果；（e）異常的ETM1 結果；（f）異常的ECF1 結果；（g）異常的ETA2 結果；（h）異常的ETM2 結果；（i）異常的ECF2 結果；（j）異常的ETA3 結果；（k）異常的ETM3 結果；（l）異常的ECF3 結果.Fig.1 Profiles showing different edge identifications from original gravity data（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）ETA1of the data in（a）；（e）ETM1of thedata in（a）；（f）ECF1of the data in（a）；（g）ETA2of the data in（a）；（h）ETM2of the data in（a）；（i）ECF2of the data in（a）；（j）ETA3of the data in（a）；（k）ETM3of the data in（a）；（l）ECF3of the data in（a）.

從（7）式中可以看出，二階垂直導數可以采用兩個二階水平導數之和來完成，水平導數采用空間域方法來計算得到，不會增大噪聲的干擾.公式中所涉及的其它低階垂直導數的計算采用傅里葉變換來完成.

3 理論模型試驗

為了試驗增強型均衡濾波器的可行性，建立如下地質模型：地下布設了2個中心埋深分別為15m和20m、與圍巖密度差為ρ＝2g／cm3、半邊長為10m的棱柱體.圖2a為模型理論重力異常，白虛線標識地質體的真實水平位置.分別利用上述的方法對重力異常進行處理（圖2）.

圖2b表示異常THDR的結果，THDR法不能很清晰地給出地質體的邊界信息；圖2c表示異常theta map的結果，theta map能夠清晰地識別出地質體的邊界，但是邊界比較發散，導致兩個地質體的邊界已經相連.高階導數對地質體的邊界具有更高的分辨率，因此可以通過增加導數的階次來使地質體邊界更加清晰，考慮采用一階垂直導數構建theta map濾波器，其具體表達式為：

圖2 不同方法的邊界識別結果（a）原始重力異常；（b）異常的THDR結果；（c）異常的theta map結果；（d）異常的一階垂直導數的theta map結果；（e）異常的ETA1 結果；（f）異常的ETM1 結果；（g）異常的ECF1 結果；（h）異常的ETA2 結果；（i）異常的ETM2 結果；（j）異常的ECF2 結果.Fig.2 Identifications of edges by different methods（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）theta map of the first vertical derivative of the data in（a）；（e）ETA1of the data in（a）；（f）ETM1of the data in（a）；（g）ECF1of the data in（a）；（h）ETA2of the data in（a）；（i）ETM2of the data in（a）；（j）ECF2of the data in（a）.

圖2d表示一階垂直導數的theta map（Theta2），該方法使邊界更加收斂，但是在高值外側存在明顯的低值，為異常的解釋工作帶來了困難.一階垂直導數與其水平導數組成的邊界識別濾波器（圖2e，2f，2g）能很好地完成邊界識別工作，且未產生干擾異常，識別出的邊界很清晰，采用本文方法構建的邊界識別工具避免了采用高階導數直接構建邊界識別濾波器所帶來的干擾.二階垂直導數與其水平導數組成的邊界識別工具（圖2h，2i，2j）能更加準確地描述地質體的水平位置，且也不存在干擾異常.從該試驗中可以看出，增強型均衡濾波器能準確地識別出地質體的水平位置，識別出的邊界非常清晰，能很好地完成異常的解釋工作.

為了進一步試驗方法的有效性，建立如下較為復雜的模型：重力異常由四個長方體組成，埋深分別為20m（1），30m（2），30m（3），5m（4）.圖3a為理論重力異常，圖中白虛線表示地質體的真實水平位置，分別利用上述邊界識別方法對該異常進行處理（圖3）.

圖3b表示異常THDR的結果，該方法能大致地給出地質體的邊界信息，但是埋深較深的地質體邊界不是很清晰；圖3c表示異常theta map的結果，該方法能給出較大的地質體邊界，而不能給出較小地質體邊界，且邊界比較發散.一階增強型均衡濾波器（圖3d，3e，3f）能很清晰地給出較大的地質體的邊界，而對于較小的地質體的邊界比較模糊.二階增強型均衡濾波器（圖3g，3h，3i）能更加準確地給出地質體的邊界信息，且較小的地質體的邊界反映也十分清晰.從圖3中可以看出，本文方法能更好地完成邊界識別工作.從圖3d，3e，3f與3g，3h，3i的對應中可以看出，隨著所使用的導數階次的增加能識別出更多的細節信息.

為了試驗方法的穩定性，在圖3a所示的重力異常中加入信噪比為50的高斯白噪聲，圖4a為加入噪聲后的重力異常，利用邊界識別方法對該異常進行處理（圖4）.

圖3 不同方法的邊界識別結果（a）原始重力異常；（b）異常的THDR結果；（c）異常的theta map結果；（d）異常的ETA1 結果；（e）異常的ETM1 結果；（f）異常的ECF1 結果；（g）異常的ETA2 結果；（h）異常的ETM2結果；（i）異常的ECF2 結果.Fig.3 Identifications of edges by different methods（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）ETA1of the data in（a）；（e）ETM1of the data in（a）；（f）ECF1of the data in（a）；（g）ETA2of the data in（a）；（h）ETM2of the data in（a）；（i）ECF2of the data in（a）.

圖4 加入噪聲后不同方法的邊界識別結果（a）原始重力異常；（b）異常的THDR結果；（c）異常的theta map結果；（d）異常的ETA1 結果；（e）異常的ETM1 結果；（f）異常的ECF1 結果；（g）由方程（7）計算的異常的ETA2結果；（h）由方程（7）計算的異常的ETM2結果；（i）由方程（7）計算的異常的ECF2結果；（j）由Fourier變換計算的異常的ETA2結果；（k）由Fourier變換計算的異常的ETM2結果；（l）由Fourier變換計算的異常的ECF2結果.Fig.4 Identifications of edges by different methods after adding noise（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）ETA1of the data in（a）；（e）ETM1of the data in（a）；（f）ECF1of the data in（a）；（g）ETA2of the data in（a）computed by the Eq.（7）；（h）ETM2of the data in（a）computed by the Eq.（7）；（i）ECF2of the data in（a）computed by the Eq.（7）；（j）ETA2of the data in（a）computed by Fourier transformation；（k）ETM2of the data in（a）computed by Fourier transformation；（l）ECF2of the data in（a）computed by Fourier transformation.

圖4b表示異常THDR的結果，由于噪聲的干擾，該方法已經無法給出地質體的邊界；圖4c表示異常theta map的結果，該方法依舊可以給出地質體的邊界.一階增強型均衡濾波器（圖4d，4e，4f）能很清晰地給出地質體的邊界，受噪聲影響較小.圖4g—4i為采用方程（7）計算得到的異常的ETA2、ETM2和ECF2結果，該方法計算出的二階增強型均衡濾波器仍能很清晰地給出地質體的邊界.圖4j—4l為采用常規Fourier變換計算得到的ETA2、ETM2和ECF2結果，但是采用該方法計算出的二階增強型均衡濾波器根本無法給出地質體的邊界，邊界已經被損壞.因此采用方程（7）進行高階垂直導數的計算能有效地提高輸出結果的穩定性.從對比中可以看出，本文提出的增強型均衡濾波器的穩定性不比現有的僅由一階導數組成的邊界識別濾波器差，即使在存在噪聲的情況下依舊能得到穩定的結果.

在應用本文方法進行磁異常處理前，要對磁異常進行化極處理，因為磁數據及其導數均受磁化方向的干擾［18］，因此采用化極異常進行解釋所獲得的結果將更加準確.

4 實際應用效果

為了驗證方法在實際中的應用效果，對四川地區重力異常和朱日和地區磁異常進行處理.圖5a為四川地區重力異常，重力數據采自國家測繪總局編繪的1∶100萬布格重力異常圖，并用黑虛線在圖中標識出已知斷裂的水平位置.利用上述方法對重力數據進行處理，來得到該地區的線性構造和地層之間的界線（圖5）.

圖5 同圖3，但為四川地區重力異常邊界識別結果Fig.5 Same as Fig.3，but for edge identification of gravity anomalies in Sichuan

圖6 同圖3，但為朱日和地區磁力異常邊界識別結果Fig.6 Same as Fig.3，but for edge identification of magnetic anomalies in the Zhurihe area

THDR法（圖5b）與theta map法（圖5c）劃分斷裂的能力較差，僅能給出四川地區部分已知斷裂的位置且不是十分清晰.一階增強型均衡濾波器（圖5d，5e，5f）能清晰地給出斷裂的分布特征，根據該結果可以很容易地劃分出斷裂的位置及走勢，且與已知斷裂對應較好.二階增強型均衡濾波器的結果（圖5g，5h，5i）與已知斷裂也存在較好的對應，此外還發現了更多的小的斷裂及小范圍的異常，能得到更多的細節信息.

下面試驗本文方法在磁異常數據處理中的作用，對中國西北部朱日和地區的磁異常進行處理（圖6）.圖6a為該地區化極后磁異常，并根據前期的研究工作劃定出已知成礦帶的大致范圍.

從圖6中可以看出，增強型均衡濾波器的邊界識別效果要明顯好于其它邊界識別技術，能更加清晰地顯示出成礦帶的界線，且與前期劃定的成礦帶范圍擬合較好，為進一步開發提供了有力的保證.二階增強型均衡濾波器能發現更多小范圍的異常，很好地描述了異常的細節特征.

5 結 論

本文提出一種新的邊界識別濾波器，稱為增強型均衡濾波器，其利用不同階導數之間的組合來進行地質體邊界的識別.通過理論模型試驗證明增強型均衡濾波器相對于現有的邊界識別方法具有更高的分辨率，能更加清晰和準確地識別出地質體的邊界，且隨著增強型均衡濾波器階次的增加會得到更多的細節信息.在計算過程中為了降低二階垂直導數計算所帶來的噪聲干擾，引入了一種計算二階垂直導數的穩定算法，使濾波器的輸出結果更加穩定.最后將其應用于四川盆地重力異常及朱日和地區磁異常的解釋中，獲得了區域地質構造的水平位置，并發現了更多的小范圍構造.增強型均衡濾波器是一種非常有效的邊界識別工具，具有良好的應用前景.

（References）

［1］Evjen H M.The place of the vertical gradient in gravitational interpretations.Geophysics，1936，1（1）：127－136.

［2］Hood P J，Teskey D J.Aeromagnetic gradiometer program of the Geological Survey of Canada.Geophysics，1989，54（8）：1012－1022.

［3］Cordell L.Gravimetric expression of Graben faulting in Santa Fe country and the Espanola Basin，New Mexico.New Mexico Geol.Soc.Guidebook，30thField Conf.，1979：59－64.

［4］Cordell L，Grauch V J S.Mapping basement magnetization zones from aeromagnetic data in the San Juan basin，New Mexico.／／Hinzc W J ed.The Utility of Regional Gravity and Magnetic Anomaly.Society of Exploration Geophysicists，1985：181－197.

［5］Roest W R，Verhoef J，Pilkington M.Magnetic interpretation using the 3－D analytic signal.Geophysics，1992，57（1）：116－125.

［6］Hsu S，Sibuet J C，Shyu C.High－resolution detection of geologic boundaries from potential field anomalies：an enhanced analytic signal technique.Geophysics，1996，61（2）：373－386.

［7］Cooper G R J.Balancing images of potential－field data.Geophysics，2009，74（3）：L17－L20.

［8］Miller H G，Singh V.Potential field tilt—a new concept for location of potential field sources.Journal of Applied Geophysics，1994，32（2－3）：213－217.

［9］Rajagopalan S，Milligan P.Image enhancement of aeromagnetic data using automatic gain control.Exploration Geophysics，1995，25（4）：173－178.

［10］Verduzco B，Fairhead J D，Green C M，et al.New insights into magnetic derivatives for structural mapping.The Leading Edge，2004，23（2）：116－119.

［11］Wijns C，Perez C，Kowalczyk P.Theta map：edge detection in magnetic data.Geophysics，2005，70（4）：L39－L43.

［12］Cooper G R J，Cowan D R.Enhancing potential field data using filters based on the local phase.Computers ＆Geosciences，2006，32（10）：1585－1591.

［13］Cooper G R J，Cowan D R.Edge enhancement of potentialfield data using normalized statistics.Geophysics，2008，73（3）：H1－H4.

［14］駱遙，王明，羅鋒等.重磁場二維希爾伯特變換——直接解析信號解釋方法.地球物理學報，2011，54（7）：1912－1920.Luo Y，Wang M，Luo F，et al.Direct analytic signal interpretation of potential field data using 2－D Hilbert transform.Chinese J.Geophys.（in Chinese），2011，54（7）：1912－1920.

［15］馬國慶，杜曉娟，李麗麗.利用水平與垂直導數的相關系數進行位場數據的邊界識別.吉林大學學報（地球科學版），2011，41（S1）：345－348.Ma G Q，Du X J，Li L L.Edge detection of potential field data using correlation coefficients of horizontal and vertical derivatives.Journal of Jilin University （Earth Science Edition）（in Chinese），2011，41（S1）：345－348.

［16］Ma G，Li L.Edge detection in potential fields with the normalized total horizontal derivative.Computers ＆Geosciences，2012，41：83－87.

［17］Fedi M，Florio G.Detection of potential fields source boundaries by enhanced horizontal derivative method.Geophysical Prospecting，2001，49（1）：40－58.

［18］Li X.Understanding 3Danalytic signal amplitude.Geophysics，2006，71（2）：L13－L16.