改進的均衡濾波器在位場數(shù)據(jù)邊界識別中的應(yīng)用

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

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130021

1 引 言

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

本文提出增強型均衡濾波器來進行地質(zhì)體邊界的識別，該濾波器是由不同階導(dǎo)數(shù)構(gòu)建的，其能更加清晰和準(zhǔn)確地識別出地質(zhì)體的邊界.

2 增強型均衡濾波器

現(xiàn)有的邊界識別濾波器均是同階導(dǎo)數(shù)之間的組合，增強型濾波器是利用不同階導(dǎo)數(shù)之間的非線性組合，其輸出結(jié)果更加準(zhǔn)確和清晰.本文給出了三種不同形式的增強型均衡濾波器.

2.1 增強型tilt angle

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

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

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

其中，n表示垂直導(dǎo)數(shù)的階數(shù)，也以此表示濾波器的階數(shù)，

系數(shù)cn是為了使公式（3）滿足數(shù)學(xué)運算法則，并可均衡不同階導(dǎo)數(shù)之間的強度.增強型tilt angle為垂直導(dǎo)數(shù)與其水平導(dǎo)數(shù)的組合，其分子與分母是不同階的，這是其與常規(guī)tilt angle的主要區(qū)別，增強型tilt angle的最大值與地質(zhì)體的邊界相對應(yīng).

2.2 增強型theta map

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

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

增強型theta map的最大值也與地質(zhì)體邊界相對應(yīng).

2.3 增強型cosine function

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

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

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

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

圖1 不同方法邊界識別結(jié)果的剖面圖（a）原始重力異常；（b）異常的THDR結(jié)果；（c）異常的theta map結(jié)果；（d）異常的ETA1 結(jié)果；（e）異常的ETM1 結(jié)果；（f）異常的ECF1 結(jié)果；（g）異常的ETA2 結(jié)果；（h）異常的ETM2 結(jié)果；（i）異常的ECF2 結(jié)果；（j）異常的ETA3 結(jié)果；（k）異常的ETM3 結(jié)果；（l）異常的ECF3 結(jié)果.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）式中可以看出，二階垂直導(dǎo)數(shù)可以采用兩個二階水平導(dǎo)數(shù)之和來完成，水平導(dǎo)數(shù)采用空間域方法來計算得到，不會增大噪聲的干擾.公式中所涉及的其它低階垂直導(dǎo)數(shù)的計算采用傅里葉變換來完成.

3 理論模型試驗

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

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

圖2 不同方法的邊界識別結(jié)果（a）原始重力異常；（b）異常的THDR結(jié)果；（c）異常的theta map結(jié)果；（d）異常的一階垂直導(dǎo)數(shù)的theta map結(jié)果；（e）異常的ETA1 結(jié)果；（f）異常的ETM1 結(jié)果；（g）異常的ECF1 結(jié)果；（h）異常的ETA2 結(jié)果；（i）異常的ETM2 結(jié)果；（j）異常的ECF2 結(jié)果.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表示一階垂直導(dǎo)數(shù)的theta map（Theta2），該方法使邊界更加收斂，但是在高值外側(cè)存在明顯的低值，為異常的解釋工作帶來了困難.一階垂直導(dǎo)數(shù)與其水平導(dǎo)數(shù)組成的邊界識別濾波器（圖2e，2f，2g）能很好地完成邊界識別工作，且未產(chǎn)生干擾異常，識別出的邊界很清晰，采用本文方法構(gòu)建的邊界識別工具避免了采用高階導(dǎo)數(shù)直接構(gòu)建邊界識別濾波器所帶來的干擾.二階垂直導(dǎo)數(shù)與其水平導(dǎo)數(shù)組成的邊界識別工具（圖2h，2i，2j）能更加準(zhǔn)確地描述地質(zhì)體的水平位置，且也不存在干擾異常.從該試驗中可以看出，增強型均衡濾波器能準(zhǔn)確地識別出地質(zhì)體的水平位置，識別出的邊界非常清晰，能很好地完成異常的解釋工作.

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

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

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

圖3 不同方法的邊界識別結(jié)果（a）原始重力異常；（b）異常的THDR結(jié)果；（c）異常的theta map結(jié)果；（d）異常的ETA1 結(jié)果；（e）異常的ETM1 結(jié)果；（f）異常的ECF1 結(jié)果；（g）異常的ETA2 結(jié)果；（h）異常的ETM2結(jié)果；（i）異常的ECF2 結(jié)果.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 加入噪聲后不同方法的邊界識別結(jié)果（a）原始重力異常；（b）異常的THDR結(jié)果；（c）異常的theta map結(jié)果；（d）異常的ETA1 結(jié)果；（e）異常的ETM1 結(jié)果；（f）異常的ECF1 結(jié)果；（g）由方程（7）計算的異常的ETA2結(jié)果；（h）由方程（7）計算的異常的ETM2結(jié)果；（i）由方程（7）計算的異常的ECF2結(jié)果；（j）由Fourier變換計算的異常的ETA2結(jié)果；（k）由Fourier變換計算的異常的ETM2結(jié)果；（l）由Fourier變換計算的異常的ECF2結(jié)果.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的結(jié)果，由于噪聲的干擾，該方法已經(jīng)無法給出地質(zhì)體的邊界；圖4c表示異常theta map的結(jié)果，該方法依舊可以給出地質(zhì)體的邊界.一階增強型均衡濾波器（圖4d，4e，4f）能很清晰地給出地質(zhì)體的邊界，受噪聲影響較小.圖4g—4i為采用方程（7）計算得到的異常的ETA2、ETM2和ECF2結(jié)果，該方法計算出的二階增強型均衡濾波器仍能很清晰地給出地質(zhì)體的邊界.圖4j—4l為采用常規(guī)Fourier變換計算得到的ETA2、ETM2和ECF2結(jié)果，但是采用該方法計算出的二階增強型均衡濾波器根本無法給出地質(zhì)體的邊界，邊界已經(jīng)被損壞.因此采用方程（7）進行高階垂直導(dǎo)數(shù)的計算能有效地提高輸出結(jié)果的穩(wěn)定性.從對比中可以看出，本文提出的增強型均衡濾波器的穩(wěn)定性不比現(xiàn)有的僅由一階導(dǎo)數(shù)組成的邊界識別濾波器差，即使在存在噪聲的情況下依舊能得到穩(wěn)定的結(jié)果.

在應(yīng)用本文方法進行磁異常處理前，要對磁異常進行化極處理，因為磁數(shù)據(jù)及其導(dǎo)數(shù)均受磁化方向的干擾［18］，因此采用化極異常進行解釋所獲得的結(jié)果將更加準(zhǔn)確.

4 實際應(yīng)用效果

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

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

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

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

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

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

5 結(jié) 論

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

（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，F(xiàn)airhead 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］駱遙，王明，羅鋒等.重磁場二維希爾伯特變換——直接解析信號解釋方法.地球物理學(xué)報，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］馬國慶，杜曉娟，李麗麗.利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進行位場數(shù)據(jù)的邊界識別.吉林大學(xué)學(xué)報（地球科學(xué)版），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，F(xiàn)lorio 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.