基于界面反演增強(qiáng)的位場邊緣識(shí)別方法

馮旭亮，魏澤坤

(1.西安石油大學(xué) 陜西省油氣成藏地質(zhì)學(xué)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710065； 2.西安石油大學(xué) 地球科學(xué)與工程學(xué)院，陜西 西安 710065)

0 引言

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

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

為提高弱異常(深部異常體)的邊緣識(shí)別效果，一些學(xué)者研究了均衡類方法。馬國慶等[27]利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進(jìn)行位場邊緣識(shí)別，并基于不同階水平導(dǎo)數(shù)構(gòu)造了多個(gè)增強(qiáng)型的位場邊緣識(shí)別方法[28-29]。此外，王彥國等[30]利用不同階次的垂向?qū)?shù)的比值進(jìn)行位場邊緣檢測，Du等[31]則利用垂向?qū)?shù)的平均值與標(biāo)準(zhǔn)差的相關(guān)系數(shù)進(jìn)行位場邊緣檢測，于平等[32]將水平方向解析信號(hào)及其垂向?qū)?shù)與傳統(tǒng)的均衡邊界識(shí)別方法結(jié)合形成了位場邊界探測方法。郭燦文等[33]利用水平導(dǎo)數(shù)與垂向?qū)?shù)的標(biāo)準(zhǔn)差的相關(guān)系數(shù)識(shí)別了磁源邊界。Pham等[34]利用總水平導(dǎo)數(shù)構(gòu)建邏輯函數(shù)并據(jù)此識(shí)別位場邊界，可明顯增強(qiáng)邊界識(shí)別的效果。

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

1 方法原理

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

i=1,2,…,N

(1)

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

(2)

式中非線性函數(shù)fi(pj)用來計(jì)算第j個(gè)棱柱體在第i點(diǎn)處引起的重力異常，其表達(dá)式為：

(3)

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

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

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

計(jì)算密度界面反演結(jié)果總水平導(dǎo)數(shù)(THDR)：

(4)

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

計(jì)算總水平導(dǎo)數(shù)的n階垂向?qū)?shù)：

(5)

式中：n為垂向?qū)?shù)的階次，階次越高，計(jì)算結(jié)果的水平分辨率越高，但計(jì)算結(jié)果不穩(wěn)定，階次一般選1或2即可。

計(jì)算總水平導(dǎo)數(shù)的峰值：

PTHDR(x,y,z)=

(6)

計(jì)算總水平導(dǎo)數(shù)峰值與總水平導(dǎo)數(shù)的比值：

VDR_THDR(x,y,z)=

(7)

計(jì)算總水平導(dǎo)數(shù)垂向?qū)?shù)的極大值并利用其對(duì)總水平導(dǎo)數(shù)垂向?qū)?shù)歸一化：

(8)

2 模型試算

2.1 裂陷盆地基底模型

裂陷盆地模型如圖1a所示，該盆地由7個(gè)大小不一的次盆組成(編號(hào)為A～G)，基底最大深度為6 km。假設(shè)沉積層與基底的密度差為-0.4×103kg/m3，將盆地的沉積層剖分為水平尺寸為1 km×1 km垂直并置的棱柱體計(jì)算其重力異常，并加入均值為0、標(biāo)準(zhǔn)差為0.1×10-5m/s2的高斯白噪聲，作為盆地模型的實(shí)測重力異常，如圖1b所示。理論重力異常基本能反映各次盆的位置及形態(tài)，但距離較近的幾個(gè)次盆的邊界(實(shí)質(zhì)為斷裂)無法區(qū)分。

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

本文方法的實(shí)質(zhì)是將重力異常轉(zhuǎn)換為一種與真實(shí)構(gòu)造更為接近的“異?！睆亩岣哌吘壸R(shí)別能力，因此若能找到一種與理論構(gòu)造更為接近的處理結(jié)果或?qū)χ亓Ξ惓＞植吭鰪?qiáng)結(jié)果并用NVDR_THDR進(jìn)行邊緣識(shí)別，是否也能得到理想的效果？圖2e～h分別為采用位場分離、穩(wěn)定向下延拓、傾斜角和解析信號(hào)振幅對(duì)重力異常處理并利用NVDR_THDR對(duì)異常進(jìn)行提取的結(jié)果。其中，剩余重力異常的NVDR_THDR結(jié)果(圖2e)極大值幅值較窄，也能清晰地呈現(xiàn)次盆B的東部邊界，但無法識(shí)別次盆D、E和F的南北邊界。采用穩(wěn)定向下延拓技術(shù)向下延拓3 km重力異常的NVDR_THDR結(jié)果(圖2f)與剩余異常結(jié)果較為相似，但其仍然無法反映E和F的南北邊界。傾斜角的NVDR_THDR結(jié)果(圖2g)與垂向?qū)?shù)階次為2的NVDR_THDR結(jié)果(圖2d)較為類似，但其幅值更寬一些。解析信號(hào)振幅的NVDR_THDR結(jié)果(圖2h)對(duì)于兩個(gè)規(guī)模較大的次盆A和G識(shí)別較好，但對(duì)于規(guī)模較小且距離較近的次盆無法有效識(shí)別其邊界。圖2的結(jié)果表明，本文提出的方法更有利于盆地基底斷裂的識(shí)別，尤其對(duì)于小規(guī)模斷裂構(gòu)造更為有效。

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

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

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

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

2.2 孤立形體組合模型

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

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

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

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

3 實(shí)際資料處理

為測試本文方法在實(shí)際資料處理中的效果，選擇鄂爾多斯盆地北部重力數(shù)據(jù)(圖6a)進(jìn)行計(jì)算，該數(shù)據(jù)來源于區(qū)域1∶50萬實(shí)測布格重力異常。圖6b為利用NVDR_THDR技術(shù)直接對(duì)布格重力異常的處理結(jié)果，研究區(qū)北部NVDR_THDR極大值以近EW向?yàn)橹?，其次是NE向和NEE向，極大值分布較為密集。研究區(qū)中部和南部NVDR_THDR極大值以NW向、EW向和NE向?yàn)橹鳎捎谠搮^(qū)域重力異常變化相對(duì)較為平緩，因此NVDR_THDR的極大值分布較為稀疏。研究區(qū)西部NVDR_THDR極大值密集分布，走向以近SN向?yàn)橹鳌?/p>

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

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

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

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

圖7 鄂爾多斯盆地北部位場邊緣識(shí)別結(jié)果與地震剖面對(duì)比(地震剖面據(jù)文獻(xiàn)[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 結(jié)論

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

基于密度界面增強(qiáng)的位場邊緣識(shí)別方法可以直接用于重力異常， 但不能直接用于磁力異常處理。對(duì)于磁異常，可將其轉(zhuǎn)化為磁源重力異常進(jìn)行邊緣識(shí)別，或?qū)⑦@一思路推廣至磁力異常處理之中，將密度界面反演改為磁性界面反演，之后利用NVDR_THDR進(jìn)行邊緣提取和增強(qiáng)。