基于界面反演增強的位場邊緣識別方法

馮旭亮，魏澤坤

(1.西安石油大學 陜西省油氣成藏地質學重點實驗室，陜西 西安 710065； 2.西安石油大學 地球科學與工程學院，陜西 西安 710065)

0 引言

油氣資源、深部礦體等的形成與分布往往與斷裂、火成巖體等有著直接或間接的關系，因此利用地球物理方法探測斷裂或巖體的位置在油氣藏勘查和深部找礦中起著至關重要的作用[1-4]，這些斷裂或巖體的位置可統稱為地質邊界。重、磁方法在研究地質體橫向不均勻性、特別是地質邊界方面有獨特的優勢。利用重、磁位場資料識別地質體邊界的方法可分為數理統計、數值計算和其他方法[5]，其中數值計算類方法研究和使用最多，比較有代表性方法主要有垂向導數[6-7]、總水平導數[8-9]、解析信號振幅[10-12]、傾斜角(Tilt導數)[13-15]、θ圖[16]5種基本方法，以及在這5種基本方法上發展起來的其他方法，如傾斜角總水平導數[17]、改進Tilt梯度水平總梯度[18]等。

王萬銀等[5]概括了各類位場邊緣識別方法的研究現狀，重點對數值計算類中垂向導數、總水平導數、解析信號振幅、傾斜角、θ圖這5種基本方法及在其基礎上發展的各方法原理和應用效果進行了研究。Ekinci等[19]對以上5種常用的邊緣識別方法進行了系統研究，并處理了土耳其安納托利亞半島地區的實際資料，與已知構造的對比結果證實了以上邊緣識別方法的有效性。為進一步增強位場邊緣識別方法的處理效果，特別是對深部小構造的識別能力，一種措施為對邊緣識別方法處理結果進行濾波增強[20-23]。不同基本方法的組合也是常用的提高邊緣識別效果的措施，如傾斜角解析信號模法[24]、解析信號模傾斜角法[25]等。英高海等[26]對比了15種具有代表性的邊界增強方法，認為解析信號模垂向導數、總水平導數、解析信號模傾斜角法是對于磁異常適用性較強、應用效果較好的3種方法。

為提高弱異常(深部異常體)的邊緣識別效果，一些學者研究了均衡類方法。馬國慶等[27]利用水平與垂直導數的相關系數進行位場邊緣識別，并基于不同階水平導數構造了多個增強型的位場邊緣識別方法[28-29]。此外，王彥國等[30]利用不同階次的垂向導數的比值進行位場邊緣檢測，Du等[31]則利用垂向導數的平均值與標準差的相關系數進行位場邊緣檢測，于平等[32]將水平方向解析信號及其垂向導數與傳統的均衡邊界識別方法結合形成了位場邊界探測方法。郭燦文等[33]利用水平導數與垂向導數的標準差的相關系數識別了磁源邊界。Pham等[34]利用總水平導數構建邏輯函數并據此識別位場邊界，可明顯增強邊界識別的效果。

歸一化類方法也是增強深部場源弱異常識別能力的主要措施，如歸一化總水平導數法[35-37]、歸一化θ圖法[38]。Wang等[39]提出了歸一化總水平導數垂向導數(NVDR_THDR)邊緣識別技術，該技術不但具有邊緣探測功能，并且具有邊緣增強功能，在劃分斷裂[40-42]、解釋巖體邊界[43]等方面取得了較好的應用效果。我們曾在識別裂陷盆地基底不同區域界面形態時，嘗試將直接迭代密度界面反演方法[44-45]加入到NVDR_THDR位場邊緣識別技術之中，提高了基底深部小規模斷裂的識別能力[46]，但并未對該方法的細節及應用效果進行詳細研究。本文將詳細闡述這一技術的原理，并通過不同類型的模型試算及實際資料處理，以分析該方法的正確性及應用效果。

1 方法原理

基于界面反演增強的位場邊緣識別方法由界面反演技術和位場邊緣識別技術組成。首先利用密度界面反演方法對重力異常進行處理(若為磁異常，可利用磁性界面反演方法進行，或轉換為磁源重力異常)，之后利用位場邊緣識別技術對密度界面反演結果進行處理。其中密度界面反演技術采用直接迭代法計算，計算公式為[45,47]：

i=1,2,…,N

(1)

式中pk(xi,yi)和pk-1(xi,yi)分別為第i個計算點處第k次和第k-1次的反演結果；G為萬有引力常量；Δρ(xi,yi)為第i個計算點處的密度差；go(xi,yi)和gcal(xi,yi)分別為第i個點處實測的重力異常和第k-1次迭代時該點的正演擬合重力異常。計算擬合重力異常時將密度界面與參考界面之間的物質剖分為垂直并置的棱柱體，重力計算點位于棱柱體正上方，其個數與計算點數相同，水平尺寸與計算點間距一致，計算公式為：

(2)

式中非線性函數fi(pj)用來計算第j個棱柱體在第i點處引起的重力異常，其表達式為：

(3)

式中：(dξj,dηj,dζj)為剖分棱柱體內任一微元的坐標，ξj1～ξj2、ηj1～ηj2、ζj1～ζj2分別為第j個棱柱體在3個方向的坐標范圍。

利用式(1)反演密度界面深度時，通常令p0為零開始迭代，并利用go(xi,yi)與gcal(xi,yi)的均方差小于重力數據的噪聲水平而終止迭代計算。然而，本文重點將其應用于深部場源體異常的增強，并不要求反演結果有較高的準確性，因此一般迭代3～5次即可。需要說明的是，對于多個孤立地質體引起的異常，從幾何形態而言，并非密度界面，但其引起的異常也可利用該方法進行增強。

本文選用歸一化總水平導數垂向導數(NVDR_THDR)技術作為位場邊緣識別方法，用該技術直接處理密度界面的反演結果。結合本文處理流程，其方法原理[39]如下：

計算密度界面反演結果總水平導數(THDR)：

(4)

式中：f(x,y,z)為密度界面反演結果。

計算總水平導數的n階垂向導數：

(5)

式中：n為垂向導數的階次，階次越高，計算結果的水平分辨率越高，但計算結果不穩定，階次一般選1或2即可。

計算總水平導數的峰值：

PTHDR(x,y,z)=

(6)

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

VDR_THDR(x,y,z)=

(7)

計算總水平導數垂向導數的極大值并利用其對總水平導數垂向導數歸一化：

(8)

2 模型試算

2.1 裂陷盆地基底模型

裂陷盆地模型如圖1a所示，該盆地由7個大小不一的次盆組成(編號為A～G)，基底最大深度為6 km。假設沉積層與基底的密度差為-0.4×103kg/m3，將盆地的沉積層剖分為水平尺寸為1 km×1 km垂直并置的棱柱體計算其重力異常，并加入均值為0、標準差為0.1×10-5m/s2的高斯白噪聲，作為盆地模型的實測重力異常，如圖1b所示。理論重力異?；灸芊从掣鞔闻璧奈恢眉靶螒B，但距離較近的幾個次盆的邊界(實質為斷裂)無法區分。

利用本文方法對圖1b所示的重力異常進行處理，其中在利用密度界面反演技術時，給定密度差為-0.4×103kg/m3，迭代5次的結果如圖2a所示。與原始重力異常(圖1b)相比，圖2a的反演結果中可以明顯區分7個次盆的位置及形態，邊界較為清晰，利用NVDR_THDR對密度界面反演結果做進一步處理，結果如圖2b所示，其極大值位置與理論模型中次盆的邊界較為吻合，并且也能在一定程度上區分距離較近的幾個次盆的范圍，如次盆B、C、D以及次盆E、F。利用NVDR_THDR方法直接對重力異常進行處理，結果如圖2c所示。與圖2b相比，圖2c的結果中次盆B的右邊界不完整，無法呈現E和F的南北邊界及分界線。為提高小規模構造的識別能力，在NVDR_THDR技術中垂向導數的階次選為2重新處理，結果如圖2d所示，與垂向一階導數結果(圖2c)相比，垂向導數階次為2時計算結果幅值更窄，但其仍然無法有效識別上述次盆的邊界。

本文方法的實質是將重力異常轉換為一種與真實構造更為接近的“異?！睆亩岣哌吘壸R別能力，因此若能找到一種與理論構造更為接近的處理結果或對重力異常局部增強結果并用NVDR_THDR進行邊緣識別，是否也能得到理想的效果？圖2e～h分別為采用位場分離、穩定向下延拓、傾斜角和解析信號振幅對重力異常處理并利用NVDR_THDR對異常進行提取的結果。其中，剩余重力異常的NVDR_THDR結果(圖2e)極大值幅值較窄，也能清晰地呈現次盆B的東部邊界，但無法識別次盆D、E和F的南北邊界。采用穩定向下延拓技術向下延拓3 km重力異常的NVDR_THDR結果(圖2f)與剩余異常結果較為相似，但其仍然無法反映E和F的南北邊界。傾斜角的NVDR_THDR結果(圖2g)與垂向導數階次為2的NVDR_THDR結果(圖2d)較為類似，但其幅值更寬一些。解析信號振幅的NVDR_THDR結果(圖2h)對于兩個規模較大的次盆A和G識別較好，但對于規模較小且距離較近的次盆無法有效識別其邊界。圖2的結果表明，本文提出的方法更有利于盆地基底斷裂的識別，尤其對于小規模斷裂構造更為有效。

圖1 裂陷盆地模型及其重力異常Fig.1 A rifted basin model and its gravity anomaly

利用密度界面反演方法增強位場邊緣識別效果時，界面上下的密度差是一個關鍵參數。圖2a為密度差是真實值的結果，若密度差取值不準確，是否會對位場邊緣識別結果有影響？為此，分別給定密度差為-0.3×103kg/m3和-0.5×103kg/m3進行反演計算，結果分別如圖3a和3c所示，相應地NVDR_THDR處理結果如圖3b和3d所示。當密度差不準確時，直接迭代法均能得到密度界面的起伏變化，盡管其與理論模型的真實深度差別較大，但相對起伏變化是一致的。在此情形下，密度界面反演結果的NVDR_THDR計算結果與準確密度差時的處理結果幾乎是一致的，均能較清晰地反映理論模型中次盆的邊界位置?？梢?，密度差的選擇幾乎不會影響最終的位場邊緣識別結果。

圖2 裂陷盆地模型位場邊緣識別結果Fig.2 The edge recognition results of the rifted basin model

圖3 密度差不同時裂陷盆地模型位場邊緣識別結果Fig.3 The edge recognition results of the rifted basin model with different density contrasts

2.2 孤立形體組合模型

孤立形體組合模型如圖4所示，該組合模型由7個大小不一、埋深不同的立方體組成，組合模型的觀測面為z=0 km的平面，x和y方向的范圍均為0～20 km，觀測面網格間距為0.25 km×0.25 km，共81×81個觀測點。正演計算該組合模型引起的重力異常，并加入均值為0、標準差為0.05×10-5m/s2的高斯白噪聲，最終得到的理論重力異常如圖4c所示。

利用本文方法進行位場邊緣識別時，需要將理論重力異常視為密度界面引起的重力異常并反演該“界面”的起伏變化，反演結果取決于界面上下密度差及界面平均深度。為保證反演的密度界面起伏變化與重力異常大小變化呈正相關以利于直觀分析，密度差應取為負值。盡管盆地模型試驗結果表明了密度差幾乎不影響最終的位場邊緣識別結果，但應盡可能地保證界面反演結果較淺以使其橫向分辨力較高。實際資料處理時可根據巖石物性特征選擇一個符合實際情況的密度差即可，這里選為-0.25×103kg/m3，界面平均深度取為2 km，反演結果如圖5a所示，相應的NVDR_THDR處理結果見圖5b。為便于對比，本文亦分別計算了模型理論重力異常的NVDR_THDR垂向導數階次為1和2的結果，見圖5c和5d。此外，裂陷盆地基底模型試驗中，剩余重力異常和向下延拓重力異常的NVDR_THDR結果較佳，因此，本文也對孤立模型進行相同的處理，結果分別見圖5e和5f。

與理論重力異常(圖4c)相比，密度界面反演結果(圖5a)的形態更接近理論模型的位置和形態，利用NVDR_THDR技術處理的結果(圖5b)中，除模型E的左邊界之外，淺部的5個模型C～G的邊界較為清晰，深部模型B的邊界位置也基本能識別出來。直接利用NVDR_THDR技術對理論重力異常的處理結果(圖5c)中，深部模型B的位置有所顯示，但其右側邊界偏差較大；淺部模型中，僅能完整的識別F和G的邊界。垂向導數階次為2時NVDR_THDR處理結果(圖5d)中，異常的幅值較窄，更容易確定模型邊界，且深部模型B的位置也更準確一些，但仍然無法完整識別淺部模型C、D、E的邊界。剩余重力異常NVDR_THDR結果(圖5e)與本文方法計算結果(圖5b)較為相似，但異常極大值稍寬一些，模型C的邊界也不太清晰。向下延拓重力異常NVDR_THDR結果(圖5f)可以清晰地反映淺部模型C～G的邊界， 甚至模型E的左邊界也能顯示出來，但對于深部模型B則無法有效反映，NVDR_THDR的幅值與噪聲相當。不難理解，重力異常向下延拓之后，可以突出淺部異常、壓制深部異常，因此，處理結果中只能突出淺部異常的邊界?？梢?，對于孤立形體的組合模型，本文方法的優勢不如在盆地基底斷裂識別中明顯，但其在一定程度上可平衡深部和淺部異常的識別。

圖4 孤立形體組合及其重力異常Fig.4 The combination model and its gravity anomaly

3 實際資料處理

為測試本文方法在實際資料處理中的效果，選擇鄂爾多斯盆地北部重力數據(圖6a)進行計算，該數據來源于區域1∶50萬實測布格重力異常。圖6b為利用NVDR_THDR技術直接對布格重力異常的處理結果，研究區北部NVDR_THDR極大值以近EW向為主，其次是NE向和NEE向，極大值分布較為密集。研究區中部和南部NVDR_THDR極大值以NW向、EW向和NE向為主，由于該區域重力異常變化相對較為平緩，因此NVDR_THDR的極大值分布較為稀疏。研究區西部NVDR_THDR極大值密集分布，走向以近SN向為主。

利用本文方法對研究區重力異常進行處理，反演時為得到深度較合理的密度界面，首先對布格重力異常添加常數的背景重力值，將其調整為(-120～0)×10-5m/s2之間。之后，根據整個鄂爾多斯盆地的地層密度資料，并為保證反演界面盡量較淺以利于突出局部異常，給定密度差為-0.5×103kg/m3，得到界面如圖6c所示。利用NVDR_THDR技術對該界面處理結果如圖6d所示，圖中NVDR_THDR的極大值分布更密集一些，尤其在研究區中部和南部更為明顯。研究區北部和西部的極大值分布與圖6b中的較為相似，但異常幅值更窄，局部異常得到增強。

為進一步對比位場邊緣識別結果，沿研究區的一條地震剖面[48](圖6d中黑色實線)提取了圖6b和6d的處理結果，如圖7所示。位場邊緣識別結果僅能反映斷裂或巖性邊界的平面位置，因此NVDR_THDR的極大值對應于斷裂帶中點偏上的某一個位置。圖7中，斷距較大的斷裂附近重力異常NVDR_THDR和密度界面NVDR_THDR均為極大值，但重力異常NVDR_THDR無法有效識別斷距較小的斷裂。密度界面NVDR_THDR結果中，對于斷距較小的斷裂也能識別，例如剖面長度110、140、160 km處，這一特征尤為明顯?？梢?，本文提出的基于界面反演增強的位場邊緣識別方法在實際資料處理也較為有效。需要說明的是，理論上，當斷裂埋深較淺、斷距較大時，其引起的重力異常梯度也越大，相應的NVDR_THDR的值也越大。但實際地質情況較為復雜，NVDR_THDR的大小與斷裂的規模、埋深以及斷裂兩側的地層密度差等均有關，因此NVDR_THDR值的大小并不能完全反映斷裂的埋深、規模等特征。

圖5 孤立形體組合模型位場邊緣識別結果Fig.5 The edge recognition results of the combination model

圖6 鄂爾多斯盆地北部重力異常及位場邊緣識別結果Fig.6 Gravity anomaly in northern Ordos Basin and the edge recognition results

圖7 鄂爾多斯盆地北部位場邊緣識別結果與地震剖面對比(地震剖面據文獻[48])Fig.7 The comparison between the edge recognition results and the seismic profile in northern Ordos Basin (The seismic profile is from reference [48])

4 結論

為增強位場邊緣識別方法對于深部小規模地質體的識別能力，本文提出了一種基于界面反演增強技術，首先利用密度界面反演方法對重力異常進行處理，之后將NVDR_THDR技術作為一種邊緣提取和增強技術對密度界面反演結果進行提取，其比直接對重力異常進行NVDR_THDR計算的效果更佳。裂陷盆地模型試驗證明了該方法能較好地用于盆地基底斷裂的識別，密度差的大小幾乎不影響最終的邊緣識別結果。將該方法應用于孤立形體組合模型試驗，處理結果也與理論模型較為吻合，與其他增強方法相比，這一方法具有一定的優勢。鄂爾多斯盆地實際資料處理結果也表明了本文方法可以更好地識別盆地基底小規模斷裂，實際應用效果較好。

基于密度界面增強的位場邊緣識別方法可以直接用于重力異常， 但不能直接用于磁力異常處理。對于磁異常，可將其轉化為磁源重力異常進行邊緣識別，或將這一思路推廣至磁力異常處理之中，將密度界面反演改為磁性界面反演，之后利用NVDR_THDR進行邊緣提取和增強。