林濤，曾昭發(fā)，于平，周帥，焦健

吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

重力勘探方法即根據(jù)地質(zhì)體的密度差異測(cè)量其在地表引起的重力異常，解釋地質(zhì)體場(chǎng)源幾何參數(shù)和物性參數(shù)的方法。重力勘探被廣泛應(yīng)用于區(qū)域地質(zhì)調(diào)查、礦區(qū)普查、工程地質(zhì)及水文勘察等領(lǐng)域[1-2]。隨著重力梯度儀的不斷發(fā)展，航空重力梯度測(cè)量技術(shù)由于其測(cè)量分辨率高、測(cè)量效率高等優(yōu)勢(shì)，在油氣資源勘查、礦產(chǎn)資源勘探等方面逐漸發(fā)揮更大的作用。綜合利用重力及其梯度數(shù)據(jù)可對(duì)地質(zhì)體的場(chǎng)源參數(shù)反演提供更多的信息，提高現(xiàn)有解釋方法的精度[3-4]。常規(guī)重力及其梯度異常反演是基于反演理論，設(shè)置目標(biāo)函數(shù)，使用非線性或線性優(yōu)化算法，使目標(biāo)函數(shù)達(dá)到最低的方法[5-10]。但重力及其梯度反演數(shù)學(xué)上呈病態(tài)，多解性強(qiáng)，需要較多的約束條件，雖然能夠?qū)Φ叵陆Y(jié)構(gòu)進(jìn)行定量刻畫，卻很難在實(shí)際應(yīng)用中取得良好效果[11]。

Patella[12]基于歸一化互相關(guān)系數(shù)對(duì)自然電位測(cè)量數(shù)據(jù)進(jìn)行地質(zhì)目標(biāo)體空間分布的快速成像。Mauriello和Patella[13-14]分別將這種方法推廣應(yīng)用于大地電磁和重力測(cè)量領(lǐng)域。此后，相關(guān)成像方法被廣泛應(yīng)用于地球物理數(shù)據(jù)處理中[15-18]。郭良輝等[19]基于重力梯度數(shù)據(jù)建立了相關(guān)成像公式，并采用異常分離方法提高了成像精度。陳召曦等[20]、侯振隆[21]實(shí)現(xiàn)了基于GPU并行的重力及其梯度相關(guān)成像，使得算法的運(yùn)算效率大為提高。近年來(lái)，相關(guān)成像與其他地球物理數(shù)據(jù)和處理方法的結(jié)合成為趨勢(shì)[22-26]。重力及其梯度相關(guān)成像技術(shù)主要用于定位目標(biāo)體的中心埋深，由于完全忽略先驗(yàn)信息，相關(guān)成像方法的縱向分辨率較差，一般用于為精細(xì)反演提供初始模型[27]。

層源位場(chǎng)是通過(guò)異常分離方法得到的地下某一深度層產(chǎn)生的位場(chǎng)。程方道等[28]首先提出了可以分離出重力層源位場(chǎng)的插值-切割分離法。汪炳柱等[29]使用這種方法，將勝利油田某地區(qū)磁異常劃分為淺、中、深3層的層源位場(chǎng)，并經(jīng)過(guò)化極垂向一次導(dǎo)數(shù)和垂向二次導(dǎo)數(shù)的驗(yàn)證，證明該方法得到的淺部異常的可信度高。余海龍[30]首次將分離得到的層源位場(chǎng)使用積分-迭代法向下延拓至地層頂層，并將每一層分別進(jìn)行反演，提高了反演的速度和精度。

筆者在前人工作的基礎(chǔ)上，提出基于層源位場(chǎng)的重力及其梯度數(shù)據(jù)聯(lián)合相關(guān)成像的方法，將重力及其梯度數(shù)據(jù)進(jìn)行異常分離，并采用向下延拓技術(shù)得到某一深度地層引起的重力及其梯度異常，再對(duì)每一層的數(shù)據(jù)分別進(jìn)行聯(lián)合相關(guān)成像，通過(guò)合成數(shù)據(jù)和實(shí)測(cè)Humble鹽丘重力數(shù)據(jù)的試驗(yàn)，驗(yàn)證了該方法能有效提升縱向分辨率。

1 方法原理

1.1 重力及其梯度聯(lián)合三維相關(guān)成像原理

傳統(tǒng)重力數(shù)據(jù)相關(guān)成像方法的基本原理是：對(duì)地下空間采用均勻網(wǎng)格進(jìn)行剖分，對(duì)每個(gè)剖分的網(wǎng)格節(jié)點(diǎn)采用正演計(jì)算公式計(jì)算其引起的重力異常，將計(jì)算的重力異常與觀測(cè)異常進(jìn)行歸一化互相關(guān)系數(shù)的計(jì)算，基于相關(guān)系數(shù)三維成像結(jié)果刻畫地下場(chǎng)源的空間分布特征。

假設(shè)剖分網(wǎng)格的中心坐標(biāo)為(xq,yq,zq)，體積為vq，剩余密度為σq的某點(diǎn)q在一個(gè)測(cè)點(diǎn)(x,y,z)處的重力異常Δgq(x,y,z)可表示為：

(1)

式中：G為萬(wàn)有引力常數(shù)，G=6.67×10-11m3/(kg·s2)。

計(jì)算第q網(wǎng)格引起的重力異常與實(shí)測(cè)重力異常的歸一化互相關(guān)系數(shù)[19]：

(2)

式中：(xi,yi,zi)為第i個(gè)觀測(cè)點(diǎn)的坐標(biāo)，Δg(xi,yi,zi)為第i個(gè)觀測(cè)點(diǎn)的重力異常，Δgq(xi,yi,zi)為q網(wǎng)格在第i個(gè)觀測(cè)點(diǎn)位置引起的重力異常。

將式(1)代入式(2)中，有：

(3)

其中，

(4)

式中：Bq為地下q點(diǎn)的質(zhì)量對(duì)于觀測(cè)點(diǎn)(xi,yi,zi)的基函數(shù)，Cq取值范圍在-1和1之間。相關(guān)系數(shù)Cq代表了觀測(cè)的重力異常與某一地下網(wǎng)格節(jié)點(diǎn)的質(zhì)量基函數(shù)的相關(guān)程度。

郭良輝等[19]將互相關(guān)系數(shù)推廣到重力梯度，得到重力梯度的相關(guān)系數(shù)計(jì)算方法：

(5)

其中，

(6)

(7)

根據(jù)該理論，可以將重力梯度不同分量相關(guān)系數(shù)進(jìn)行聯(lián)合成像。

1.2 重力異常分離原理

重力及重力梯度測(cè)量在水平方向的分辨率較高，為提高相關(guān)成像的縱向分辨率，筆者引入徐世浙等[32]使用的插值-切割分離方法進(jìn)行不同深度重力場(chǎng)的分離，該方法計(jì)算速度快，而且能將異常近似分離為不同地層深度引起的異常，其基本思路如下。

設(shè)位場(chǎng)異常Z(x,y)由區(qū)域場(chǎng)R(x,y)和局部場(chǎng)L(x,y)組成，即：

Z(x,y)=R(x,y)+L(x,y)

(8)

令A(yù)(x,y)代表與點(diǎn)(x,y)相距r的4個(gè)測(cè)點(diǎn)測(cè)得異常的平均值：

Z(x,y+r)+Z(x,y-r)]

(9)

式中：r為切割半徑。

由式(9)可以得出區(qū)域場(chǎng)是Z(x,y)和A(x,y)的加權(quán)平均：

R(x,y)=aA(x,y)+bZ(x,y)

(10)

式中：a和b是加權(quán)系數(shù)，且a+b=1。

令：

(11)

(12)

ΔBx=Z(x+r,y)-Z(x-r,y)

(13)

ΔBy=Z(x,y+r)-Z(x,y-r)

(14)

(15)

(16)

e=c+d(0≤e≤2)

(17)

取加權(quán)系數(shù)：

(18)

由此得到一次切割場(chǎng)R1(x,y)。對(duì)R1(x,y)重復(fù)以上步驟，不斷迭代下去，有：

(19)

于是有：

(20)

式中：R(x,y)為切割半徑為r的區(qū)域場(chǎng)。由式(8)可以得到最終局部場(chǎng)：

L(x,y)=Z(x,y)-R(x,y)

(21)

根據(jù)該方法在實(shí)際應(yīng)用中的經(jīng)驗(yàn)，可以近似地認(rèn)為，切割半徑為r時(shí)得到的局部場(chǎng)是深度r以上的地層產(chǎn)生的異常場(chǎng)。在此基礎(chǔ)上，余海龍[30]提出了層切割技術(shù)，即用不同的兩個(gè)切割半徑對(duì)原始異常進(jìn)行切割，用切割結(jié)果的局部場(chǎng)做減法，其結(jié)果便可近似地表征為地下某具有上下底面的層位在地表引起的層源異常。利用這個(gè)技術(shù)，將重力數(shù)據(jù)劃分為不同深度的地層引起的多組重力異常，從而在數(shù)據(jù)處理過(guò)程中提高縱向分辨率。

1.3 重力異常向下延拓法原理

用切割法分離異常之后所得的異常場(chǎng)是地面測(cè)得的地下不同深度層的異常場(chǎng)。若直接采用分層相關(guān)成像的方法，深部底層距測(cè)點(diǎn)仍較遠(yuǎn)，會(huì)產(chǎn)生不必要的誤差。所以在切割分離之后，采用向下延拓的方法，將測(cè)量面下延至切割層的頂層，此時(shí)再采用相關(guān)成像的方法，可以更加精確地分辨層內(nèi)密度特征。

本文采用的向下延拓方法為徐世浙[33]提出的迭代法，其原理如下：

假設(shè)存在兩個(gè)平面ΓA和ΓB，用u(x,y,0)代表平面的位，用U(kx,ky,0)=F[u(x,y,0)]表示u(x,y,0)的Fourier變換，則z=h(ΓB)平面的位可表示為：

(22)

式中：Γ-1是Fourier反變換。

(23)

式中：s是步長(zhǎng)，一般為1。

重復(fù)進(jìn)行以上過(guò)程，可以得到迭代公式：

(24)

一般迭代20～50次即可。這種方法將不穩(wěn)定的向下延拓問(wèn)題轉(zhuǎn)換為穩(wěn)定的向上延拓問(wèn)題，使得延拓的結(jié)果更加穩(wěn)定，并且能延拓更大的深度。經(jīng)過(guò)異常分離和向下延拓后，可以將某一深度范圍內(nèi)產(chǎn)生的重力及重力梯度異常進(jìn)行聯(lián)合相關(guān)成像，得到更為精確和易于解釋的結(jié)果。

2 模型驗(yàn)證

2.1 傳統(tǒng)聯(lián)合相關(guān)成像

設(shè)計(jì)由11個(gè)不同大小、不同密度的棱柱體組成的地質(zhì)模型，其具體參數(shù)如表1所示。選擇y=8 000 m測(cè)線為研究對(duì)象，地質(zhì)體大致分為3層，A1、A2層埋深為1 500～3 000 m；B1、B2、B3層埋深為400～800 m；C1、C2層為100～200 m。其俯瞰圖和y=8 000 m截面圖見(jiàn)圖1。

表1 復(fù)合模型的幾何參數(shù)和剩余密度

對(duì)組合體的Δg、Gxz、Gyz、Gzz分別進(jìn)行網(wǎng)格為20 000 m×20 000 m×3 000 m，步長(zhǎng)為m的三維相關(guān)成像，選取y=8 000 m測(cè)線作剖面圖。為使圖像具有一致性，文中的相關(guān)系數(shù)結(jié)果均經(jīng)過(guò)歸一化處理。由圖2可以看出，關(guān)于重力異常的相關(guān)成像可以突出底部較大地質(zhì)體的異常A1、A2，而關(guān)于重力梯度的相關(guān)成像方法對(duì)中等深度的地質(zhì)體B1、B2、B3有異常反應(yīng)，但它們對(duì)淺層地質(zhì)體C1、C2的信號(hào)都不明顯。

對(duì)Gxz、Gyz、Gzz三個(gè)分量進(jìn)行聯(lián)合相關(guān)成像，得到圖3。圖像明顯收斂，B3產(chǎn)生的異常雖然較小，但精確地反映了其位置，仍存在以下問(wèn)題：

①A1的質(zhì)心比A2要高，推測(cè)是測(cè)線在A1的中間部分而在A2的邊緣部分，使得A2的異常值降低更大。

②受A1、A2影響比較大的B1、B2的異常消失了。

圖2 y=8 000 m測(cè)線單分量相關(guān)成像Fig.2 Single-component probability tomography of y=8 000 m line

圖3 y=8 000 m測(cè)線聯(lián)合相關(guān)成像Fig.3 Joint probability tomography of y=8 000 m line

可見(jiàn)聯(lián)合相關(guān)成像能夠突出主要異常，但會(huì)忽略較小地質(zhì)體的異常。為突出淺部較小地質(zhì)體的異常，可以先進(jìn)行異常分離，將分離后的數(shù)據(jù)再進(jìn)行相關(guān)成像。圖4是進(jìn)行3×3均值濾波分離異常之后的多分量相關(guān)成像結(jié)果，中層B1、B2、B3和淺層C1、C2的地質(zhì)體異常更加明顯，雖然能根據(jù)最大值大致判斷其質(zhì)心的橫縱向位置關(guān)系，但對(duì)中淺層地質(zhì)體的深度分布卻很模糊，縱向分辨率并不高。

此時(shí)再對(duì)Gxz、Gyz、Gzz三個(gè)分量進(jìn)行聯(lián)合成像，所得到的結(jié)果如圖5所示，雖然淺層地質(zhì)體C1、C2的異常值減弱了，但通過(guò)最大值位置的判斷，可以清晰地分辨出B1、B2、B3以及C1、C2的質(zhì)心深度的相對(duì)關(guān)系，并能推斷其形態(tài)，但仍存在如下問(wèn)題：

①對(duì)C1、C2的質(zhì)心深度判斷并不準(zhǔn)確，推測(cè)是由于測(cè)點(diǎn)距為200 m，而C層位的分布深度在100～200 m導(dǎo)致精度不夠。

圖4 y=8 000 m測(cè)線基于均值分離的單分量相關(guān)成像Fig.4 Single-component probability tomography of y=8 000 m line based on mean separation

圖5 y=8 000 m測(cè)線基于均值分離的聯(lián)合相關(guān)成像Fig.5 Joint probability tomography of y=8 000 m line based on mean separation

②C2的位置比實(shí)際向下，推測(cè)是B2和A2的干擾所致。

盡管如此，聯(lián)合相關(guān)成像的精度還是顯著提高了。

根據(jù)聯(lián)合相關(guān)成像在組合體模型中的實(shí)驗(yàn)，可以看出，相對(duì)于傳統(tǒng)相關(guān)成像方法，聯(lián)合相關(guān)成像的方法可以壓制干擾，使得成像結(jié)果更加收斂，更加清楚地分辨地質(zhì)體的質(zhì)心深度。同時(shí)，該方法也會(huì)壓制一些實(shí)際存在的較小異常。

2.2 基于層源位場(chǎng)的聯(lián)合相關(guān)成像

對(duì)上述組合體模型采用切割分離的方法，將位場(chǎng)異常分離為不同深度的地質(zhì)體引起的異常，再向下延拓至地層頂面后進(jìn)行相關(guān)成像。研究切割厚度較大時(shí)的圖像特征，將地層切割為兩層，0～1 000 m為一層，1 000～2 000 m為一層，分別對(duì)其進(jìn)行傳統(tǒng)相關(guān)成像和聯(lián)合相關(guān)成像，得到圖6和圖7。

可以看出，在切割范圍較大的情況下，傳統(tǒng)成像方法受干擾較大，雖然其異常范圍可以與模型相對(duì)應(yīng)，但假影較多，難以解釋。而聯(lián)合相關(guān)成像方法的結(jié)果清晰。圖6中可以看到B層的埋深在600 m±，且其存在于0～1 000 m地層的下部，A1、A2的埋深在1 500 m和2 000 m，分布于1 000～2 000 m地層的下部，與模型較為符合。

圖6 0～1 000 m地層y=8 000 m測(cè)線相關(guān)成像Fig.6 Probability tomography of y=8 000 m line within 0～1 000 m strata

圖7 1 000～2 000 m地層y=8 000 m測(cè)線相關(guān)成像Fig.7 Probability tomography of y=8 000 m line within 1 000～2 000 m strata

圖8 y=8 000 m測(cè)線多地層聯(lián)合相關(guān)成像Fig.8 Multi-strata joint probability tomography of y=8 000 m line

在此基礎(chǔ)上，為得到更精確的層源信息，將地層進(jìn)一步精細(xì)切割，并對(duì)切割的結(jié)果進(jìn)行聯(lián)合相關(guān)成像(圖8)。由此可推斷出地質(zhì)體的頂?shù)酌媛裆睿?/p>

①0～200 m地層中，只有C層位引起異常。由于頂界面埋深小于測(cè)點(diǎn)距，無(wú)法得出C層地質(zhì)體在該地層的分布位置。記C層埋深為0～200 m。

②200～400 m地層中，異常主要由B層引起，可以明顯看出該層異常存在于地層下部，可記其頂面埋深為300 m。

③400～800 m地層中，異常主要由B層地質(zhì)體引起。

④800～1 200 m地層中，沒(méi)有主要異常，故可將B層地質(zhì)體的底面埋深記為800 m。

⑤1 200～1 600 m的地層中，A1的異常出現(xiàn)，記A1的頂面埋深為1 200 m。

⑥1 600～2 400 m的地層中，A2的異常出現(xiàn)，記其頂面埋深為1 600 m。

由于切割成像和向下延拓法均有深度限制，為保證切割數(shù)據(jù)的精確性，只將地層切割至2 400 m的深度。

由此可以推斷出地質(zhì)體頂?shù)酌媛裆睿?/p>

淺層：0～200 m(實(shí)際100～200 m)。

中層：300～800 m(實(shí)際400～800 m)。

深層(頂部埋深)：A1為1 200 m；A2為1 600 m(實(shí)際均為1 500 m)。

由于在數(shù)據(jù)處理過(guò)程中明確了層位信息，基于層源位場(chǎng)的重力及其梯度聯(lián)合成像方法在模型實(shí)驗(yàn)中，精度高于傳統(tǒng)相關(guān)成像方法。由于異常分離和延拓的過(guò)程產(chǎn)生噪音，單獨(dú)數(shù)據(jù)的相關(guān)反演方法并不能在層源位場(chǎng)數(shù)據(jù)取得較好效果。

3 Humble鹽丘實(shí)測(cè)重力數(shù)據(jù)應(yīng)用

將本文方法應(yīng)用到美國(guó)Houston地區(qū)Humble鹽丘數(shù)據(jù)中。圖9是鹽丘的重力及重力梯度異常平面圖，采樣間隔為100 m。為研究鹽丘的縱向分布特征，對(duì)該測(cè)區(qū)的重力梯度數(shù)據(jù)分別按1 000～3 000 m，3 000～5 000 m，5 000～7 000 m進(jìn)行異常分離與延拓，將3個(gè)地層數(shù)據(jù)分別進(jìn)行步長(zhǎng)為100×100×100的多分量聯(lián)合相關(guān)成像，并畫出y=6 000 m截面的剖面圖(圖10)?？梢钥闯?，該鹽丘的中心埋深為4 km ，頂面在3 km。表2為本方法得到的鹽丘埋深與其他方法的對(duì)比。可以看出，該方法得到的埋深與其他方法相近，說(shuō)明重力梯度多分量聯(lián)合相關(guān)成像的反演結(jié)果是可靠的。

4 結(jié)論

(1)相較于單分量相關(guān)成像，聯(lián)合相關(guān)成像可以有效壓制假影，使成像結(jié)果更為聚焦，同時(shí)也會(huì)壓制一些真實(shí)影像。

表2 重力及其梯度聯(lián)合相關(guān)成像與其他方法的對(duì)比

圖9 Humble鹽丘重力及重力梯度異常平面圖Fig.9 Gravity and gravity gradiometry anomalies in Humble salt dome

圖10 Humble鹽丘測(cè)區(qū)y=6 000 m測(cè)線多地層聯(lián)合相關(guān)成像剖面圖Fig.10 Multi-strata joint probability tomography of y=6 000 m line within Humble salt dome survey area

(2)單分量相關(guān)成像在處理層源位場(chǎng)數(shù)據(jù)時(shí)會(huì)產(chǎn)生大量假影，聯(lián)合相關(guān)成像依然聚焦。

(3)對(duì)層源位場(chǎng)數(shù)據(jù)進(jìn)行聯(lián)合相關(guān)成像，可以在提升縱向分辨率的同時(shí)，抑制假影的產(chǎn)生，得到更準(zhǔn)確的結(jié)果。