位場垂向梯度最佳自比值的邊界檢測技術

王彥國，張鳳旭，劉 財，王祝文，張 瑾

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

2東華理工大學核工程與地球物理學院，南昌 330013

1 引 言

應用位場數據識別地質體邊界是地質解釋的一項重要工作.位場異常中包含有場源邊界的信息，但邊界信息的提取需要進行相關的數據處理.為此，地球物理工作者提出了許多基于位場梯度的邊界檢測方法.這些方法按邊界檢測的計算模式，可分為導數分析和數理統(tǒng)計兩類［1］.

導數分析法主要有垂向導數分析法［2］、水平總梯 度［3］、tilt angle［4］及 其 水 平 總 梯 度［5］、Theta map［6］等方法.垂向導數分析法是Evjen提出的，主要是利用垂向一、二階導數的零值線確定地質體邊界［2］；Cordell提出了水平梯度法，利用水平梯度模極大值檢測場源邊界［3］；這兩種方法有效地提高了對原異常信息的識別能力，但難以同時反映出深度不同的地質體邊界［7］.Miller等提出了tilt angle法，該法識別地質體邊界的原理與垂向導數類似［4］；Verduzco在tilt angle法基礎上提出了tilt angle水平總梯度法，用計算結果的極大值檢測地質體邊界，但在計算過程中采用了異常的二階導數，更易受高頻干擾的影響［5］.Wijns等采用水平總梯度與一階解析信號比值（Theta map）來分析磁源邊界，對低緯度磁性體有較好的邊界檢測能力，但檢測的邊界與深部地質體的實際邊界位置有偏差［6］.

數理統(tǒng)計法主要有Euler反褶積［8－10］、小子域濾波［11］、導數歸一化標準差［12］、均方差比歸一化垂向梯度［13］等方法.Euler反褶積法是將位場異常和其梯度聯(lián)系在一起的Euler齊次方程，該方法可用于地質體邊界的識別，但方程式靈敏度較高，導數計算結果稍有偏差，便會在計算結果中產生較大的偏離，從而影響反演的精度［14］.小子域濾波法是楊高印針對異常分離提出的，可用于邊界的識別，但只對原異常中以梯級帶形式顯示的邊界位置有效［11］；張鳳旭等［15－16］在小子域濾波基礎上，提出了三方向小子域濾波法，獲得了較為詳細的斷裂構造信息，但計算結果與實際地質體邊界也有一定的偏差.Cooper等［12］借助于異常導數的標準差，提出了可以檢測不同深度地質體邊界的歸一化標準差法；王彥國等［13］在2012年提出了均方差比歸一化垂向梯度的邊界檢測方法，對不同深度的地質體邊界也有較好的探測效果，但這兩種方法在計算的過程中同樣采用了導數計算，因此也會受噪聲等高頻干擾.

顯然，上述邊界檢測方法大都是以位場梯度計算為基礎的，方法運算簡單，物理意義也明確，但在實際資料處理中都存在著受高頻干擾的影響［14］.在處理這個問題上，如果直接采用低通濾波進行降噪，則會損失部分有用信號.

在利用位場異常進行場源邊界識別中，一階導數分析是目前常用的一類方法.不過筆者認為，導數的階次越高，對場源的邊界識別能力和定位精度就越強，然而高階導數易受噪聲等干擾的影響，而且導數的階次越高，高頻干擾的影響就越為嚴重.如果能采用某種數學處理手段，既能盡量減少導數異常中有效邊界信息的損失，又具有抗干擾能力的話，勢必可以提高場源的邊界識別能力.為此，本文提出了位場垂向梯度最佳自比值的邊界檢測方法.

2 方法原理分析

2.1 自比值的定義

對于位場異常數據f（i，j），設其垂向m階導數為fz（m），那么可以將導數異常fz（m）窗口內數據的均值與其均方根的比值定義為fz（m）的一次自比值κ，即

2.2 自比值數值分析

構造如下二次函數Fij（x）：

2.3 位場垂向梯度自比值的物理意義

設垂向導數fz（m）在所選定的計算窗口內數據平均值為：

2.4 最佳自比次數和最佳窗口長度的確定

依據文獻［17－19］，可定義垂向m階導數的第n次自比值與第n＋1次自比值的歸一化互相關系數R為：

其中M、N分別為測線數和每條測線的測點數.

事實上，在選定的窗口下，計算的自比值互相關系數與自比次數的關系曲線中存在一個明顯極大值（后文將進一步闡述），現將這個極大值點對應的自比次數定義為最佳次數，將該次數對應的自比值稱為最佳自比值.

由于在選定的初始窗口下，最佳自比值未必是確定地質體邊界位置的最佳結果，因此需要選擇最佳窗口的長度.最佳窗口長度的選取方案是：首先設定一個置信度R0（其取值一般介于0.96和0.99之間），如果在初始窗口下，自比值互相關系數的極大值max（R）大于給定的置信度R0，則認為最佳自比值去噪能力強，結果是可信的，那么初始窗口即為最佳窗口；若max（R）不大于R0，則認為初始窗口下的自比值結果可信度較低，需調整計算窗口的長度，重新計算自比值互相關系數與自比次數的關系，直到某一窗口下滿足max（R）＞R0時，終止計算，這時所對應的窗口大小便是最佳計算窗口長度.

2.5 最佳自比值計算步驟

由以上分析知，本文提出的垂向梯度最佳自比值法在給定置信度R0之后，可以自適應地尋找到最佳窗口下的最佳自比次數，然后通過最佳自比值進行場源邊界的識別.為了方便起見，在此給出算法的步驟：

（1）采用快速傅里葉變換計算位場數據的垂向梯度fz（m）；

（2）給定初始窗口長度（一般選取D＝5），利用公式（1）和（2）計算fz（m）的n次自比值

（4）給定置信度R0，將初始窗口下的互相關系數極大值max（R）與R0進行比較.若max（R）＞R0，則輸出最佳自比值；若max（R）≤R0，則調整窗口長度，重新計算（2）、（3）步驟，直到滿足max（R）＞R0為止.

值得說明的是：為了進一步提高地質體邊界的定位精度，可采用楊高印提出的小子域濾波［11］對最佳自比值進行濾波輸出.

3 理論模型分析

為了驗證本文方法的邊界識別能力，設計了一個包含4個不同參數的長方體組合模型（模型體參數見表1，位置見圖1a），其產生的理論重力異常見圖1b.同時為了驗證方法消除噪聲的能力，在疊加模型體產生的重力異常中加入了隨機噪聲（圖1c）.

從模型重力異常（圖1b）中可以看出，地質體A由于埋深和規(guī)模較大，其產生的異常在邊界位置表現出的梯級帶范圍較寬，因此利用該異常所顯示的邊界特征是無法精確確定地質體邊界位置的；地質體B、C、D規(guī)模較小，受地質體A的影響，其產生的異常在疊加異常中表現為等值線同形扭曲，顯然只依靠這種異常特征難以確定出這三個地質體的邊界；而且在加入噪聲后，異常圖（圖1c）在三個地質體位置所顯示的等值線同形扭曲特征變的模糊不清，這進一步增加了邊界識別的難度.

表1 模型參數Table 1 The model parameters

在模型試驗中，筆者先用tilt angle、水平總梯度、Theta map、導數歸一化標準差、小子域濾波和均方差比歸一化垂向梯度6種已有的邊界識別方法對含噪聲的重力異常分別進行了計算.各種方法計算的結果見圖2，從中可以看出，由于受噪聲和異常疊加的影響，6種方法的計算結果均無法清晰地顯示模型體A的邊界，更無法有效地識別其它3個較試驗表明，該極值點就是確定最佳窗口大小和自比次數的重要依據.

在計算最佳窗口和最佳自比次數（圖4）時，首先選擇初始窗口D＝5，由于該窗口條件下垂向一階和二階導數的自比值互相關系數極大值均大于可信度R0（文中選取R0＝0.98），因此可確定垂向一階、二階導數的最佳窗口與初始窗口相同（D＝5），相對應的最佳自比次數分別為n＝3和n＝4.對于垂向三階導數，由于在初始窗口下的互相關系數極大值max（Rz（3））＝0.956，小于R0，因此將計算窗口長度增加至D＝7，此窗口下的max（Rz（3））＝0.994，大于R0，即垂向三階導數自比值的最佳窗口為D＝7，其對應的最佳自比次數為n＝4.

圖3 含噪聲重力異常垂向導數（a）垂向一階導數；（b）垂向二階導數；（c）垂向三階導數.Fig.3 Vertical derivatives of gravity anomaly including random noise（a）First－order vertical derivative；（b）Second－order vertical derivative；（c）Third－order vertical derivative.

圖4 自比值互相關系數與自比次數的關系曲線Fig.4 The relationship of cross－correlation coefficient of two successive auto－ratios and the number of auto－ratio

以上對比分析充分證明，本文提出的邊界檢測方法不但具有較高的精度，而且計算過程更為穩(wěn)定.

4 實際重力資料1）應用

為了驗證垂向梯度最佳自比值對實際資料的處理效果，選取了吉林省南部鴨綠江盆地2）實測重力數據進行試驗.

鴨綠江盆地位于中朝板塊東北緣，二級大地構造單元隸屬于遼東臺隆區(qū)，主體位于太子河—渾江坳陷內（圖6）.在盆地內，除志留紀、泥盆紀和早石提供了較好的對比依據.

從研究區(qū)布格重力異常（圖8a）中可以看出，反映地質體邊界位置或斷裂構造的梯級帶主要分布在二道江—新立屯—板石鎮(zhèn)一線、四道江—六道江—白山—孫家堡子一線、孤砬子—紅土崖—三道湖—石人一線、公益村—大路村一線以及螞蟻河—鬧枝溝屯圈閉區(qū)等位置，這幾組異常梯級帶可能是大型地質體邊界（或規(guī)模較大的斷裂構造）的反映.但這些異常梯度帶較平緩，邊界位置難以直接從重力異常圖中精確定位.另外，較小型地質體邊界受區(qū)域異常影響較大，在異常圖中主要表現為異常等值線突然變寬或變窄以及同形扭曲等非梯級帶特征，這類地質體邊界位置確定難度更大.

圖8b為結合區(qū)域地質和前人地質－地球物理工作成果勾劃的構造分區(qū)圖3）.圖中給出的大型斷裂位置、燕山期花崗巖出露區(qū)、中新生代凹陷分布范圍以及渾江煤田工作區(qū)等信息，可為文中方法的有效性提供佐證.

圖8c為研究區(qū)垂向三階導數自比值互相關系數與自比次數的關系曲線，仿照圖4做法，可檢測出最佳窗口長度D＝5和最佳自比次數n＝4，進而獲得了垂向三階導數最佳自比值的計算結果（圖8d）.從中可以看出，最佳自比值不僅清晰地反映出龍崗隆起和渾江坳陷之間以及渾江坳陷和老嶺推覆區(qū)之間的大型斷裂，還較好地顯示出了不同巖性間的接觸帶（圖7）.另外，自比值圈定的負異常較好地反映出了研究區(qū)中－新生代地層和燕山期花崗巖等相對低密度巖石的分布，也客觀地反映出了渾江煤田的工作范圍.這些結果與已知地質資料以及前人的工作成果符合較好，有力地佐證了方法的有效性.

圖7 鴨綠江盆地地質圖Fig.7 Geology map of Yalujiang basin

5 結 論

有效地利用位場相應階次的導數異常，能夠提高地質構造解釋的分辨率，但高階導數換算對干擾的放大作用一直困擾著地球物理工作者.本文在前人工作的基礎上，另辟蹊徑，提出了位場垂向梯度最佳自比值的邊界檢測方法.在文中，筆者對方法進行了數值分析，并闡述了該算法檢測地質體邊界的物理機制.模型試驗和實際算例表明，垂向梯度最佳自比值算法不但能夠有效地消除導數中的干擾成分，而且能夠精細地識別出地質體的邊界.

本文實現了高階導數在位場邊界識別中的應用.由于自比值可以有效地壓制干擾，因此可以提供更豐富、精確的邊界信息，為高階導數在位場數據處理中的有效應用提供了新的思路.

圖8 （a）鴨綠江盆地布格重力異常（單位10－5 m·s－2）；（b）研究區(qū)構造分區(qū)圖；（c）垂向三階導數自比值互相關系數與自比值次數關系曲線；（d）垂向三階導數最佳自比值Fig.8 （a）Bouguer gravity anomaly of Yalujiang basin（unit：10－5 m·s－2）；（b）The tectonic division of research area；（c）The relationship of cross－correlation coefficient of two successive auto－ratios of third－order vertical derivative and the number of auto－ratio；（d）Optimal auto－ratio of third－order vertical derivative

［1］ 王萬銀，邱之云，楊永等.位場邊緣識別方法研究進展.地球物理學進展，2010，25（1）：196－210.Wang W Y，Qiu Z Y，Yang Y，et al.Some advances in the edge recognition of the potential field.ProgressinGeophys.（in Chinese），2010，25（1）：196－210.

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

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

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

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

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

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

［8］ Thompson D T.EULDPH：A new technique for making computer－assisted depth estimates from magnetic data.Geophysics，1982，47（1）：31－37.

［9］ Reid A B，Allsop J M，Granser H，et al.Magnetic interpretation in three dimensions using Euler deconvolution.Geophysics，1990，55（1）：80－91.

［10］ Hansen R O，Laura S.Multiple－source Euler deconvolution.Geophysics，2002，67（2）：525－535.

［11］ 楊高印.位場數據處理的一項新技術——小子域濾波法.石油地球物理勘探，1995，30（2）：240－244.

Yang G Y.A new technique for potential－field data processing：small subdomain filtering.OilGeophysical Prospecting（in Chinese），1995，30（2）：240－244.

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

［13］ 王彥國，王祝文，張鳳旭等.基于均方差比歸一化垂向梯度法的位場邊界檢測.中國石油大學學報（自然科學版），2012，36（2）：86－90.

Wang Y G，Wang Z W，Zhang F X，et al.Edge detection of potential field based on normalized vertical gradient of mean square error ratio.JournalofChinaUniversityof Petroleum（EditionofNaturalScience），2012，36（2）：86－90.

［14］ 張恒磊，劉天佑，楊宇山.各向異性標準化方差計算重磁源邊界.地球物理學報，2011，54（7）：1921－1927.

Zhang H L，Liu T Y，Yang Y S.Calculation of gravity and magnetic source boundary based on anisotropy normalized variance.ChineseJ.Geophys.（in Chinese），2011，54（7）：1921－1927.

［15］ 張鳳旭，張鳳琴，劉財等.斷裂構造精細解釋技術——三方向小子域濾波.地球物理學報，2007，50（5）：1543－1550.

Zhang F X，Zhang F Q，Liu C，et al.A technique for elaborate explanation of faulted structures：three－directional small subdomain filtering.ChineseJ.Geophys.（in Chinese），2007，50（5）：1543－1550.

［16］ 張鳳旭，張興洲，張鳳琴等.中國東北地區(qū)重力場研究——利用改進的三方向小子域濾波劃分主構造線及大地構造單元.地球物理學報，2010，53（6）：1475－1485.

Zhang F X，Zhang X Z，Zhang F Q，et al.Study of gravity field in Northeastern China area：Classification of main structure lines and tectonic units using the improved threedirectional small subdomain filtering.ChineseJ.Geophys.（in Chinese），2010，53（6）：1475－1485.

［17］ 高光珠，李忠武，余理富等.歸一化互相關系數在圖像序列目標檢測中的應用.計算機工程與科學，2005，27（3）：38－40.

Gao G Z，Li Z W，Yu L F，et al.Application of the normalized cross correlation coefficient in image sequence object detection.ComputerEngineeringandScience（in Chinese），2005，27（3）：38－40.

［18］ 曾華霖，許德樹.最佳向上延拓高度的估計.地學前緣，2002，9（2）：499－504.

Zeng H L，Xu D S.Estimation of optimum upward continuation height.EarthScienceFrontiers（in Chinese），2002，9（2）：499－504.

［19］ 孟小紅，劉國鋒，陳召曦等.基于剩余異常相關成像的重磁物性反演方法.地球物理學報，2012，55（1）：304－309.

Meng X H，Liu G F，Chen Z X，et al.3－D gravity and magnetic inversion for physical properties based on residual anomaly correlation.ChineseJ.Geophys.（in Chinese），2012，55（1）：304－309.