重力數據LTHG 均衡邊界識別方法對比分析及其應用

彭莉紅，張偉盟，程莎莎，孫棟華，駱燕，陳偉

1 核工業航測遙感中心，河北 石家莊 050002

2 中核集團鈾資源地球物理勘查技術（重點實驗室），河北 石家莊 050002

3 河北省航空探測與遙感技術重點實驗室，河北 石家莊 050002

邊界識別方法在重力數據解釋中起著至關重要的作用，國內外學者對它進行了廣泛研究，提出多種邊界識別方法［1-16］。眾多邊界識別方法以導數分析類為主，而導數分析類中最基本的方法是總水平梯度（THG）［1-3］、垂向導數（VDR）［4-6］和解析信號（AS）［7-9］，但上述三種方法均存在異常幅度隨地質體埋深的增加而迅速衰減的情況，因此對于埋深較深、異常幅度較小的地質體，邊界識別能力稍顯不足。近年來，多種通過求取重力數據導數比值來均衡不同幅度異常的邊界識別方法相繼被提了出來。Miller 和Singh［17］提出了均衡不同幅度異常的邊界識別方法——傾斜角法（TDR），該方法對于場源體的深度不敏感，對埋深不同的多場源異常體的邊界識別效果較好；Verduzco 等［18］針對TDR 的處理結果峰值位于場源體中心，識別結果不直接的問題，提出了傾斜角的總水平導數法（THDR），該方法采用極大值識別場源邊界，結果相較于TDR 簡單、直接，且具有更高的分辨率，然而也存在隨場源深度的增加、異常幅度降低的情況［19］；Wijns 等［20］提出另一種均衡不同振幅的識別方法，利用解析信號與總水平導數的比值來進行地質體的邊界識別，即Theta圖法（TM），該方法相對于解析信號法的分辨能力增強，缺點是所識別出的邊界較模糊。Ferreira 等人［19］提出了總水平梯度的傾斜角法（TTHG），該方法能較好地識別出不同埋深的場源體邊界信息；Cooper［21］提出了另一種改進的傾斜角法，稱為解析信號的傾斜角法（TAS），雖然TAS 可以更好地識別出埋深較深的場源體邊界信息，但存在與AS 法同樣的不足，即隨埋深的增加，其極大值位置會向實際邊界的內部收斂。王萬銀等［22］提出了歸一化總水平導數垂向導數邊界識別方法（NVDR_THDR），該方法采用零閾值技術消除地質體邊界，結果直觀，易于識別，是均衡不同埋深場源邊界識別效果的有效方法。馬國慶等提出增強型均衡濾波器，該濾波器利用不同階導數之間的組合識別地質體邊界。該方法能有效識別深、淺部地質體邊界，但高階導數的引入使噪聲不可避免地被放大，從而影響整體邊界識別的穩定性［23-24］。Luan 等提出基于Logistic 函數及總水平導數的重力數據邊界增強方法（LTHG 法），認為該方法處理結果簡單、直接且對于場源邊界具有更高的分辨能力。LTHG 方法是目前最新提出的可用于重力數據的邊界識別方法，但還未有其與主流邊界識別方法效果對比的相關研究成果，也缺少實際應用方面的公開論述。

為此，本文基于理論模型及實測數據，對比分析了不同情況下LTHG 法與常用的THG、AS、TDR、THDR、TTHG、TAS、TM 和NVDR_THDR8 種方法對不同埋深場源體的均衡能力、分辨能力及抗噪能力，進而驗證LTHG 法的適用性及應用效果。

1 LTHG 方法簡介

LTHG 法是Luan 等（2019）新提出的一種基于Logistic 函數及總水平導數的位場數據均衡邊界識別方法，其核心原理是基于Logistic 函數，通過對總水平導數的一階垂向導數和水平導數的比值進行Logistic 函數運算來達到突出顯示不同深度的場源體邊界的效果。

Logistic 函數是機器學習領域中常用的一種激活函數，定義域（-∞，+∞），值域為（0，1），形式如下：

Logistic 函數（函數曲線形式見圖1）在定義域內單調連續，呈現出先緩慢增長，然后加速增長，最后逐漸穩定的趨勢，與常用于識別重力數據邊界的arctan 函數相似。

圖1 Logistic 函數變化曲線Fig. 1 The Logistic function curve

基于Logistic 函數，LTHG 法的計算公式如下：

式（2）中： THG—重力數據的總水平導數；a—常數。LTHG 法求解后得到一個比值，是無量綱，其平面峰值表征了地質體的邊界。

圖2 為二維情況下，采用LTHG 法a分別取值為0.5、1.0、5.0、10、20、50 和100 時直立長方體密度模型的邊界識別結果，可以看出a 取不同參數時，均可以較好的識別出模型的邊界位置，且邊界位置處異常幅值基本相同，但不同取值計算的結果還是存在一定的差異，當a＞10 時，計算結果當中，邊界位置處異常峰型會變的更加尖銳，但當a取值大于10 時，處理結果中會因存在畸變點而出現空值的情況，當a＜1 時，模型邊界位置處異常峰型平緩，1≤a＜10 時LTHG 法應用效果較好。

圖2 不同參數LTHG 法邊界識別結果Fig. 2 The LTHG filter of gravity anomaly due to the 2D block with different a.

2 邊界識別方法對比

為說明LTHG 法識別邊界的有效性，選取八種常用的重力數據的邊界識別方法進行對比，即總水平梯度（THG）、解析信號（AS）、傾斜角法（TDR）、傾斜角的總水平導數法（THDR）、Theta 圖法（TM）、總水平梯度傾斜角法（TTHG）、解析信號傾斜角法（TAS）和歸一化總水平導數垂向導數（NVDR_THDR）法，各方法定義如下：

THG 法［25］：

AS 法［25］：

TDR 法［25］：

THDR 法［25］：

TM 法［25］：

TAS 法［25］：

TTHG 法［25］：

NVDR_THDR 法［22］：

計算總水平導數（THG）的n階垂向導數（VDRn），n—垂向導數階數，n=1，2，3…。階數越大，橫向分辨力越高，通常對于重磁異常，階數取2 時較為合適；

使用閾值大于0 來計算總水平導數峰值（PTHDR）：

計算總水平導數峰值及總水平導數的比值：

計算總水平導數垂向導數的最大值（VDR - THDRmax），并使用最大值進行總水平導數垂向導數歸一化，最終得到歸一化總水平導數垂向導數（NVDR-THDR）：

建立了由5 個不同規模、不同埋深的長方體組成的組合模型，模型三維立體圖見圖3，參數見表1，正演計算了3 種情況下組合模型的理論重力異常。情況1 為所有長方體的剩余密度均為正值，長方體物體1 和2 剩余密度為0.5×10-3kg·m-3，3、4、5 體的剩余密度為0.4×10-3kg·m-3，計算的重力異常如圖4A 所示。情況2 為長方體1、4 剩余密度為負值，體1 剩余密度為-0.5×10-3kg·m-3，體4 剩余密度為-0.4×10-3kg·m-3，體2、3、5 剩余密度與第一種情況保存一致，計算的理論重力異常如圖5A 所示。情況3 為在情況2 數據中加入了2 %的高斯噪聲，結果如圖6A 所示。選用包括THG、AS、TDR、THDR、TM、TAS、TTHG、NVDR_THDR 及LTHG在內的九種方法對以上三種情況下的理論重力異常數據進行處理，對比分析各方法的不同振幅重力異常的均衡能力、邊界分辨能力及抗干擾能力。

表1 組合模型參數Table 1 Parameters of the synthetic model

圖3 模型三維立體圖Fig.3 A 3D view of the synthetic model

圖4 模型理論重力異常（重力正異常）及邊界識別結果Fig. 4 Synthetic gravity anomaly of the first model and the edges detected by different methods

圖5 模型理論重力異常（重力正、負異常）及邊界識別結果Fig. 5 Synthetic gravity anomaly of the second model and the edges detected by different methods

圖6 加入噪聲的模型理論重力異常（重力正、負異常）及邊界識別結果Fig. 6 Synthetic gravity anomaly of the third model and the edges detected by different methods

圖4B、C 為THG、AS 的處理結果，可以看出當存在多個埋深不同的場源體時，THG 和AS不能均衡不同埋深的場源體邊界信息，隨著深度增加，場源體邊界變得越來越不清晰，同時雖然薄板1 和2 埋深較淺，但AS 仍無法很好的顯示出它們的邊界信息。圖4D～I 分別為TDR、THDR、TM、TAS、TTHG 和NVDR_THDR的處理結果，雖然TDR 均衡了不同埋深的場源體的邊界信息，但處理結果的峰值位于場源體中心，處理結果不直接。THDR 法克服了TDR法的不足，利用極大值來識別場源體的邊界，但對于埋深較深的場源體，邊界分辨能力仍稍顯不足。TM 法識別的邊界位置不夠準確，且存在假的邊界信息。TAS 法同樣利用極大值來識別場源體邊界，然而對于較深的場源，其峰值會從實際的邊界位置向內偏移，使識別到的邊界位于真實邊界內；對于較薄場源，TAS 的峰值直接位于場源體中心。TTHG 法較好地均衡了不同埋深的場源體邊界信息，識別出的邊界信息接近真實的場源體邊界，但隨著場源埋深的增加，識別的邊界信息也存在一定程度的擴散。NVDR_THDR 法可直接用峰值識別出邊界，效果比上述各方法更好。圖4J 為 LTHG 法的處理結果，可以看出LTHG 很好地克服了導數計算隨深度增加而迅速衰減的問題，均衡了不同埋深場源體的邊界信息，結果與NVDR_THDR 相似，但LTHG 方法對邊界的分辨能力要更強。對比分析結果表明：LTHG 不僅能更清晰、準確地識別出場源邊界，而且能給出較高的邊界分辨率。

為對比分析上述各方法同時識別正負異常邊界的能力，選用包括THG、AS、TDR、THDR、TM、TAS、TTHG、NVDR_THDR 及LTHG在內的9 種方法對第2 種情況下的理論重力異常數據進行了處理。圖5B、C 為THG、AS 的處理結果，同樣THG 法和AS 法不能同時均衡振幅不同的異常。圖5D～I 分別為TDR、THDR、TM、TAS、TTHG 和NVDR_THDR 的處理結果，由圖5D 可知，當同時存在正負異常時，TDR 法無法準確識別場源的邊界，且結果中出現了假邊界信息。THDR 法也存在假邊界信息。TM法除不能很好地識別出埋深較深或薄的場源體邊界外假邊界信息更加顯著。TAS 法能避免假邊界信息的產生，但識別到的邊界位于真實邊界內。TTHG 能避免假邊界信息的產生，但識別的邊界信息隨著場源埋深的增加有一定程度的擴散。NVDR_THDR 法識別出的邊界信息輪廓清晰、位置準確，無虛假邊界信息。圖5J 為 LTHG 法的處理結果，可以看出LTHG 法識別的邊界位置接近真實場源邊界，無假邊界信息，相較于NVDR_THDR，LTHG 方法對邊界的分辨能力要更強。

為對比分析各方法的抗噪能力，在第二種情況的數據中加入幅度為異常振幅2 %的高斯噪聲。圖6B～J 分別為THG、AS、TDR、THDR、TM、TAS、TTHG、NVDR_THDR 及LTHG 的處理結果。THG、AS 法對淺源體3 和4 的邊界識別效果相對較好，但對深源體5 的邊界識別效果卻相對模糊，同時對于薄板1 和2 ，AS 法也不能很好地發揮作用。TDR 和TM 依然無法對場源體邊界進行準確的識別，受噪聲的影響THDR、TAS 法識別的場源邊界變的模糊，分辨率降低，TTHG、NVDR_THDR 和LTHG 三種方法依然可以較準確的識別出場源體邊界位置，且NVDR_THDR 和LTHG 受深度變化的影響小，未出現明顯的邊界向內或外側偏移的情況，相對而言LTHG 法識別的位置更清晰，分辨率更高。

3 應用實例

埃塞俄比亞歐加登盆地是非洲東海岸油氣勘探前景最好的盆地之一，是發育在古生代結晶基底之上的一個古內陸裂谷—被動大陸邊緣盆地，根據基底構造特征，可進一步劃分為6 個次一級構造單元：北部隆起、北部斜坡、西部斜坡、中部坳陷、南部斜坡和南部隆起（圖7），目前在中部坳陷帶發現了C、H 和 D 三個氣田，展示了盆地良好的勘探潛力［26-30］。但盆地的勘探還處于早期評價階段，對斷裂的整體空間延伸及展布規律還認識不清。對重力數據及其轉換參數進行分析，可以為研究盆地斷裂構造發育特征及其與油氣成藏的關系提供重要依據［31-33］。

圖7 歐加登盆地構造單元劃分及研究區位置（據參考文獻［30］）Fig. 7 Location and simplified geological maps of the study area

核工業航測遙感中心和中國石油集團東方地球物理有限責任公司于2018 年在歐加登盆地中部共同開展了1：10 萬航空重力測量，獲得一套高精度航空重力數據。本文為驗證LTHG 法實際應用效果，選取歐加登盆地中部坳陷帶作為研究區（圖7），對研究區的航空布格重力數據進行LTHG 處理，并與THG、AS、TDR、THDR、TM、TAS、TTHG 和NVDR_THDR在內的幾種方法進行對比分析，為降低數據噪聲的影響，對布格重力異常數據進行了向上延拓，延拓高度為2 km，經延拓處理后識別的場源體邊界會更加光滑，但場源體的形態不會發生改變；LTHG 法a取值為5（a為5 時，識別結果足夠尖銳且無空值現象）。從各方法的識別結果可以看出，研究區斷裂走向整體呈NE 和NEE 向，與地質資料相吻合，但各方法的識別結果及識別精度上存在一定差異。由圖8A、B可見，THG 法和AS 法對幅值較大的異常識別效果較好，對幅值較低的識別效果稍顯不足；圖8C 顯示，TDR 雖然均衡了不同埋深的場源體的邊界信息，但無法產生銳利的邊界信息；圖8D、E 和F 顯示THDR、TM 和TAS 識別效果不甚理想，假邊界信息較多。圖8G、H 和I 顯示TTHG、NVDR_THDR 和LTHG 三種方法的識別效果較好，處理得到的重力異常線性特征更加明顯，同時又顯示更多的地質細節，但相對而言LTHG 識別的結果更加直觀，分辨率更高。

圖8 研究區航空布格重力異常邊界識別結果Fig.8 The edge recognition results of airborne gravity data in the study area

結合研究區的區域地質資料和其他物探資料，以LTHG 方法為主，TTHG、NVDR_THDR等方法為輔助，推斷了研究區內的斷裂共10 條（表2），以近NE 向和近NW 向為主，其中，NE 向斷裂5 條，分別為 F1、F2、F3、F4、F5，NW 向斷裂5 條，為 F6、F7、F8、F9、F10。 將本次的邊界識別結果與前人文獻［30］進行對比分析，如圖9 所示。本次推斷的NE 向斷裂與前人研究成果基本一致，同時新推斷了一組NW 向斷裂，從布格重力數據及邊界識別結果來看，NW 向斷裂基本錯切了NE 向斷裂；前人研究成果中，中部凹陷、中央凸起及西部凹陷北界斷裂與本次推斷的F2 斷裂一致，南界斷裂與本次推斷F3 斷裂一致。

表2 研究區斷裂推斷成果表Table 2 The results of inference fault in the study area

圖9 研究區推斷斷裂與前人解釋成果對比圖Fig.9 Comparison of interpretation faults and previous interpretation results in the study area

4 結 論

1）理論模型的合成重力數據處理結果顯示，LTHG 法能很好的識別出不同埋深的場源體邊界，提供更豐富的場源體信息。與其他常用方法進行比，LTHG 法的識別的結果準確、清晰，分辨率更高，且具有一定的抗干擾能力；

2）理論模型試驗表明，采用LTHG 法進行邊界識別時，邊界位置不受a取值的影響，但當a＞10 時，識別結果會出現數據缺失的情況，a＜1 時，模型邊界位置處異常峰型平緩，因此1≤a＜10 時，LTHG 法應用效果較好；

3）對實測航空重力數據處理得到的重力異?？傮w上與研究區的斷裂有較好的相關性，同時顯示了更多的地質細節，最終以LTHG 法處理的結果為主，并綜合其他方法對研究區的斷裂進行了推斷解釋，推斷NE 向斷裂與前人研究成果基本一致，新推斷一組NW 向斷裂，且NW 向斷裂錯切了NE 向斷裂。