雙程聲波方程逆時(shí)深度偏移

何兵壽,張會(huì)星,范國苗

(中國海洋大學(xué) 海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島 266100)

0 前言

地震波偏移是地震勘探領(lǐng)域的核心技術(shù),工業(yè)界對(duì)地震波場疊前偏移的要求主要有二方面:

(1)復(fù)雜構(gòu)造的精確成像。

(2)提供疊前反演所需的道集。

在理論上,基于雙程波方程的地震資料,疊前逆時(shí)深度偏移技術(shù)能夠同時(shí)達(dá)到以上目標(biāo),因此該項(xiàng)技術(shù)在提出之初就展現(xiàn)出良好的工業(yè)應(yīng)用前景。近年來,隨著油氣生產(chǎn)對(duì)地震勘探精度要求的不斷提高,勘探目標(biāo)構(gòu)造與巖性復(fù)雜程度的不斷增加,業(yè)界開始逐漸研究并采用實(shí)用的地震波疊前逆時(shí)深度偏移技術(shù)來解決問題,以達(dá)到提高勘探精度與實(shí)現(xiàn)巖性勘探等目標(biāo)。

多年來,國內(nèi)、外地球物理學(xué)者在逆時(shí)偏移領(lǐng)域進(jìn)行了大量的研究工作,并取得了許多成果[1～5]。但目前為止,基于雙程波方程的逆時(shí)偏移技術(shù)在工業(yè)生產(chǎn)資料的實(shí)際應(yīng)用方面還少有成功實(shí)例,限制該項(xiàng)技術(shù)進(jìn)入工業(yè)應(yīng)用領(lǐng)域的主要原因在于:

(1)地震波方程逆時(shí)延拓本身還有許多技術(shù)問題需要進(jìn)一步完善,如計(jì)算精度的進(jìn)一步提高,截?cái)噙吔鐐畏瓷涞膲褐频取?/p>

(2)雙程波方程逆時(shí)偏移的成像準(zhǔn)則與成像條件的計(jì)算方法需進(jìn)一步深入研究,在地震波方程逆時(shí)延拓過程中,地下各點(diǎn)在何時(shí)成像,如何成像是決定偏移效果的又一關(guān)鍵因素。而目前常用的,基于時(shí)間一致性準(zhǔn)則的激發(fā)時(shí)間成像條件,無法補(bǔ)償?shù)卣鸩鏁r(shí)延拓過程中的層間能量反射損失,難以達(dá)到理想的偏移效果。

(3)在雙程波方程逆時(shí)延拓過程中產(chǎn)生的層間反射干擾,需要新的方法來壓制,逆時(shí)偏移中的層間反射會(huì)使淺層產(chǎn)生強(qiáng)振幅干擾,惡化剖面淺層部份的成像結(jié)果。

作者在本文引入了波動(dòng)方程正演中的交錯(cuò)網(wǎng)格高階差分技術(shù)[6]和最佳匹配層邊界吸收技術(shù)[7],以提高聲波方程逆時(shí)延拓的精度。對(duì)于疊后逆時(shí)偏移,作者采用了基于爆炸反射界面的零時(shí)間成像條件進(jìn)行波場成像。對(duì)于疊前逆時(shí)深度偏移,應(yīng)用上行、下行波歸一化互相關(guān)成像條件,壓制了部份層間反射,并補(bǔ)償了深層能量,最終實(shí)現(xiàn)了聲波方程的逆時(shí)偏移。

1 雙程聲波方程的逆時(shí)延拓算法

1.1 雙程聲波方程逆時(shí)延拓的高階差分格式

在二維情況下,各向同性介質(zhì)中的雙程聲波方程為:

其中 vx、vz分別為質(zhì)點(diǎn)在x與z方向的振動(dòng)速度,p為位移;v為縱波速度;ρ為密度;x、z分別為空間坐標(biāo);t為時(shí)間。

借鑒一階聲波方程的正演方法算法,在圖1所示的交錯(cuò)網(wǎng)格空間中,對(duì)式(1)進(jìn)行高階差分離散,可得反射縱波逆時(shí)延拓的高階差分格式:

其中:Δx、Δz分別為x與z方向的空間離散步長;Δt為時(shí)間離散步長;N為差分階數(shù)的一半;i、j為空間離散點(diǎn)序號(hào);n為時(shí)間離散點(diǎn)序號(hào);為差分系數(shù),其值詳見參考文獻(xiàn)[6]。

以地面記錄作為邊值條件,利用式(2)即可求解地震波的逆時(shí)傳播問題。

圖1 交錯(cuò)網(wǎng)格示意圖Fig.1 The diagram staggered grids

1.2 穩(wěn)定性條件與吸收邊界條件

經(jīng)推導(dǎo),式(2)～式(4)的穩(wěn)定性條件為:

采用最佳匹配層(PML)吸收邊界條件,解決截?cái)噙吔绲膫畏瓷鋯栴}。依據(jù)PML的方程分裂思路[7],可得:

其中 p⊥和p‖分別為p分量在x與z方向的分裂算子;d(x)和d(z)分別為x和z方向的衰減因子。對(duì)于左邊界和右邊界,有d(z)=0;對(duì)于上邊界和下邊界,有d(x)=0。

同樣在交錯(cuò)網(wǎng)格空間中對(duì)式(4)進(jìn)行差分離散,可得到適用于雙程聲波方程逆時(shí)偏移的PML邊界差分格式,詳細(xì)推導(dǎo)過程這里不再贅述。

2 雙程聲波方程逆時(shí)偏移的成像條件

2.1 疊后逆時(shí)偏移的成像條件

疊后反射縱波剖面可等效為在各巖性分界面上,以不同強(qiáng)度同時(shí)激發(fā)的地震波,以速度v/2(v為地層速度)傳播至地表并被記錄的結(jié)果,這就是目前疊后偏移中常用的爆炸反射界面原理。這一原理同樣對(duì)聲波方程的疊后逆時(shí)偏移有效,故疊后逆時(shí)偏移的成像條件為:其中 I mage(x,z)為對(duì)應(yīng)網(wǎng)格點(diǎn)的偏移成像結(jié)果;為接收點(diǎn)記錄逆時(shí)延拓至t=0時(shí)刻地下各點(diǎn)的波場值;x,z為空間坐標(biāo)。

2.2 疊前逆時(shí)偏移的成像條件

在縱波方程的疊前逆時(shí)深度偏移領(lǐng)域,Claerbout[8]提出的上行、下行波互相關(guān)成像條件可概括為:其中 I mage(x,z)為對(duì)應(yīng)網(wǎng)格點(diǎn)的偏移成像結(jié)果;為接收點(diǎn)記錄逆時(shí)延拓波場值;x、z為空間坐標(biāo);tmax為記錄長度;為炮點(diǎn)子波正向延拓波場值。

式(8)實(shí)質(zhì)上是炮點(diǎn)正向延拓波場與檢波點(diǎn)逆時(shí)延拓波場的互相關(guān)運(yùn)算,對(duì)于深部地層來說,由于在檢波點(diǎn)波場逆時(shí)延拓過程中,同樣存在能量的衰減,但這種衰減不是由介質(zhì)的彈性參數(shù)所引起,而是由逆時(shí)延拓算法本身引起的,其本質(zhì)是一種誤差,且深度越大,這種誤差對(duì)偏移結(jié)果的影響也越大。為克服這一局限,可利用炮點(diǎn)波場值對(duì)式(8)進(jìn)行歸一,得到新的成像條件:

式(9)中的分母項(xiàng)與地震波正向傳播過程中的能量損失有關(guān),能量損失越大,其值越大。因此,式(9)對(duì)深層地震波能量有較強(qiáng)的補(bǔ)償作用,便于深部地層的成像,作者在本文采用這種歸一化的互相關(guān)成像條件,來實(shí)現(xiàn)雙程聲波方程的疊前逆時(shí)深度偏移。

3 雙程聲波方程逆時(shí)偏移的脈沖響應(yīng)

3.1 均勻介質(zhì)中的脈沖響應(yīng)

如圖2所示,脈沖主頻為50 Hz且各脈沖能量相同,均勻介質(zhì)模型的縱波速度為3 000 m/s。

圖2 脈沖響應(yīng)試驗(yàn)數(shù)據(jù)Fig.2 seismic data for impulse response test

圖3(a)為采用式(8)得到的雙程聲波方程逆時(shí)偏移的脈沖響應(yīng),此脈沖響應(yīng)表明雙程聲波方程逆時(shí)偏移技術(shù),具有良好的傾角適應(yīng)性。由于雙程聲波方程在推導(dǎo)過程中未做任何傾角假設(shè),故以其為理論基礎(chǔ)的逆時(shí)偏移算法,可對(duì)任意傾角的地層進(jìn)行成像。通過分析圖3(a)也可得到互相關(guān)成像條件的主要缺陷,表現(xiàn)為:深層成像結(jié)果中地震波振幅衰減嚴(yán)重,且其衰減機(jī)理與地震波正向傳播時(shí)的衰減機(jī)理相同。而事實(shí)上,震源激發(fā)的地震波,在下行傳播至反射界面過程中會(huì)產(chǎn)生能量衰減。當(dāng)下行波遇到反射界面后,產(chǎn)生的反射波在上行至地表接收點(diǎn)的過程中,也會(huì)產(chǎn)生能量的衰減,其衰減程度與波的傳播路徑和所經(jīng)過介質(zhì)的彈性性質(zhì)有關(guān)。這說明,對(duì)于地表接收到的一次反射波來說,地震波傳播過程中能量的衰減必然發(fā)生二次(上行一次,下行一次),也只會(huì)發(fā)生二次。而利用雙程聲波方程進(jìn)行逆時(shí)延拓時(shí),地震波的能量就會(huì)產(chǎn)生第三次衰減,這顯然與一次反射波的傳播機(jī)理相悖。同時(shí),地震波逆時(shí)偏移的目標(biāo),是將來自地下各反射點(diǎn)的一次反射波在成像時(shí)刻,全部聚焦到對(duì)應(yīng)的反射點(diǎn)位置處,這就要求接收點(diǎn)波場在逆時(shí)延拓過程中不會(huì)產(chǎn)生能量的損失,而基于求解雙程波方程的逆時(shí)延拓算法是無法做到這一點(diǎn)的。

圖3 均勻介質(zhì)中雙程聲波方程逆時(shí)偏移的脈沖響應(yīng)Fig.3 Reverse-time migration impulse response of two-way acoustic equation in homogeneous media

作者在本文采用式(9)作為逆時(shí)偏移的成像條件,以部份補(bǔ)償?shù)卣鸩ㄔ谀鏁r(shí)延拓過程中發(fā)生的第三次能量損失。圖3(b)為基于歸一化互相關(guān)成像條件得到的脈沖響應(yīng)。分析對(duì)比圖3(a)與圖3(b),可以得以下結(jié)論:

(1)在二種成像條件下,各波組內(nèi)部各道的能量關(guān)系未發(fā)生變化,這表明式(9)未改變逆時(shí)偏移中各成像點(diǎn)地震波動(dòng)力學(xué)特征的相對(duì)關(guān)系。

(2)在圖3(b)中,各脈沖的響應(yīng)在能量上保持一致,這表明算法補(bǔ)償了地震波逆時(shí)延拓過程中的第三次能量衰減。在本例中,由于模型為均勻各向同性介質(zhì),此時(shí)能量的衰減只可能由地震波的幾何擴(kuò)散引起,說明在均勻介質(zhì)情況下,式(9)可補(bǔ)償?shù)卣鸩ㄔ谀鏁r(shí)延拓過程中,由于球面擴(kuò)散造成的地震波能量衰減。因此,歸一化互相關(guān)成像條件更有利于深部地層的成像。在復(fù)雜模型條件下,由各種因素所造成能量衰減的補(bǔ)償方法,是今后逆時(shí)偏移領(lǐng)域的研究重點(diǎn)。作者在本文后續(xù)的疊前偏移算例中,均以式(9)作為成像條件。

3.2 層狀介質(zhì)中的脈沖響應(yīng)

研究均勻介質(zhì)情況下的偏移脈沖響應(yīng)只有理論意義,層狀介質(zhì)在油氣勘探中最為常見,為分析問題方便,作者將在本文中研究水平層狀介質(zhì)的脈沖響應(yīng)。假設(shè)地下共有四層介質(zhì),其速度分別為3 000 m/s、3 400 m/s、3 800 m/s、4 200 m/s;三個(gè)速度分界面埋深分別為500 m、1 000 m、1 500 m。圖4是圖2所示數(shù)據(jù)的疊前逆時(shí)深度偏移結(jié)果,脈沖響應(yīng)同樣顯示出了逆時(shí)偏移良好的陡傾角適應(yīng)性,但水平層狀介質(zhì)的逆時(shí)偏移脈沖響應(yīng),受到了層間反射的干擾(圖4中箭頭所指處)。層間反射是由檢波點(diǎn)波場逆時(shí)延拓過程中在速度分界面處產(chǎn)生的反射波所引起的,其本質(zhì)是一種偏移噪聲,在偏移結(jié)果中主要表現(xiàn)為近地表處出現(xiàn)許多低頻強(qiáng)振幅干擾。目前,地球物理學(xué)界針對(duì)逆時(shí)偏移中層間反射的壓制問題,進(jìn)行了許多研究工作[9～12],取得了一些成果,但該問題的根本解決還需要結(jié)合逆時(shí)偏移算法本身,或采用單程波方程來實(shí)現(xiàn)。有關(guān)層間反射的詳細(xì)壓制方法,作者暫不討論。

圖4 水平層狀介質(zhì)中雙程聲波方程逆時(shí)偏移的脈沖響應(yīng)Fig.4 Reverse-time migration impulse response of acoustic equation in horizontally layered elastic medium

3.3 連續(xù)速度介質(zhì)中的脈沖響應(yīng)

連續(xù)介質(zhì)的速度為:

式中 z為深度。

圖5是在該模型條件下的雙程聲波方程逆時(shí)偏移脈沖響應(yīng)。由圖5可見,在連續(xù)速度介質(zhì)中,逆時(shí)偏移層間反射干擾遠(yuǎn)低于存在較大速度分界面的層狀介質(zhì)情況,這是由于連續(xù)速度介質(zhì)中不存在大的波阻抗差,層間反射波的能量極弱所致。這說明對(duì)存在明顯速度分界面的層速度模型進(jìn)行適當(dāng)平滑,可壓制層間反射干擾。

圖5 連續(xù)速度介質(zhì)中雙程聲波方程逆時(shí)偏移的脈沖響應(yīng)Fig.5 Reverse-time migration impulse response of acoustic equation in continuous medium

4 模型算例

4.1 疊后逆時(shí)深度偏移模型算例

利用sigsbee_2b縱波速度模型(見圖6),檢驗(yàn)雙程聲波方程疊后逆時(shí)偏移對(duì)復(fù)雜構(gòu)造的成像能力。該模型中存在一個(gè)高速異常體、一系列斷層、二組繞射點(diǎn)(圖6中虛箭頭所指處)和若干套地層。在高速體上方8 km和12 km位置處,存在一些高速的零碎沉積物。疊后逆時(shí)偏移所用的零偏移距理論數(shù)據(jù),可以通過有限差分求解波動(dòng)方程獲得,正演所用的震源為主頻50 Hz的雷克子波。

圖7為正演得到的零偏移距剖面。在圖7中存在大量繞射波,各波組十分復(fù)雜,無法用此疊加剖面直接進(jìn)行構(gòu)造解釋。

圖6 sigsbee_2b速度模型Fig.6 Sigsbee_2b compression wave velocity model

圖7 sigsbee_2b模型正演零偏移距剖面Fig.7 Zero offset seismogram of sigsbee_2b model

圖8 圖7的疊后逆時(shí)深度偏移結(jié)果Fig.8 Post stack reverse-time migration result of Fig.7

圖8為采用本文的疊后逆時(shí)偏移方法得到的成像結(jié)果。在圖8中,海底地形、高速體頂?shù)捉缑?、高速體上部地層和模型左部的地層均得到了很好的成像,但高速體底部地層由于受高速體屏蔽的影響難以準(zhǔn)確成像,這是疊后偏移的局限,利用疊前逆時(shí)深度偏移技術(shù)可解決高速體下部地層的成像問題。

4.2 疊前逆時(shí)深度偏移算例

對(duì)sigsbee_2b模型的理論單炮數(shù)據(jù)進(jìn)行疊前逆時(shí)深度偏移,單炮記錄按以下觀測系統(tǒng),通過雙程聲波方程有限差分正演獲得:縱波源激發(fā),震源為主頻50 Hz的雷克子波,中間放炮,炮點(diǎn)二側(cè)共501道接收,測線二端固定排列不動(dòng),只移動(dòng)炮點(diǎn)位置,道間距0.005 km,炮間距0.05 km,記錄長度7.5 s,共得321炮合成記錄。

圖9(見下頁)為各單炮記錄疊前逆時(shí)深度偏移后的共成像點(diǎn)疊加剖面,高速體下部地層的成像質(zhì)量得到的明顯改善,說明基于雙程聲波方程的疊前逆時(shí)深度偏移算法,能夠?qū)?fù)雜構(gòu)造進(jìn)行成像。

5 結(jié)論及展望

波動(dòng)方程逆時(shí)偏移技術(shù)具有不受傾角限制與高保真性等優(yōu)點(diǎn),是解決復(fù)雜構(gòu)造成像問題的理想工具。作者在本文推導(dǎo)的基于交錯(cuò)網(wǎng)格的雙程聲波方程高階逆時(shí)延拓算子,與相應(yīng)的PML吸收邊界條件可實(shí)現(xiàn)地表波場的高精度逆時(shí)延拓。并且歸一化互相關(guān)成像條件,可有效補(bǔ)償?shù)卣鸩鏁r(shí)延拓過程中的能量損失。綜合應(yīng)用這二項(xiàng)技術(shù),可實(shí)現(xiàn)復(fù)雜構(gòu)造的精確成像。

層間反射的壓制技術(shù)與層速度建模技術(shù),是波動(dòng)方程逆時(shí)偏移領(lǐng)域的難點(diǎn),目前的技術(shù)無法在不損傷有效信號(hào)的前提下完全去除層間反射干擾。層間反射的壓制問題,已成為制約基于雙程波方程的逆時(shí)偏移技術(shù)進(jìn)入工業(yè)應(yīng)用領(lǐng)域的主要因素之一。該問題的最終解決,將會(huì)有力推動(dòng)逆時(shí)偏移技術(shù)進(jìn)入工業(yè)化生產(chǎn)階段,并進(jìn)一步提高地震資料成像處理的精度。層速度建模方面,在不考慮地層各向異性的前提下,逆時(shí)偏移所需的層速度,可參考其它疊前深度偏移技術(shù)的建模方法,并綜合利用地質(zhì)、測井,及其它已知資料來獲得。更為精確的速度建模技術(shù),需要針對(duì)逆時(shí)偏移算法本身的特點(diǎn)來研究。除相關(guān)理論與方法需進(jìn)一步完善外,巨大的計(jì)算量也是限制逆時(shí)偏移技術(shù)進(jìn)入工業(yè)應(yīng)用領(lǐng)域的一大障礙。但隨著計(jì)算機(jī)硬件的高速發(fā)展,特別是微機(jī)群的普及與GPU的推廣,逆時(shí)偏移技術(shù)必將在未來的油氣勘探中發(fā)揮巨大作用。

圖9 Sigsbee_2b模型的疊前逆時(shí)深度偏移結(jié)果Fig.9 Pre-stack reverse-time migration stack profile of sigsbee_2b data

[1] WHITEMORE N D.Iterative depth migration by backward time propagation[C].Expanded Abstracts of 53rdSEG Annual Int Mtg,1983.

[2] MCMECHAN GA.Migration by extrapolation of timedependent boundary values[J].Geophysical Prospecting,1983,31(2):413.

[3] CHATTOPADHYAY S,MCMECHANG A. Imaging conditions for prestack reverse-time migration[J].Geophysics,2008,73(3):S81.

[4] STEVE T H.Reverse-t ime depth migration: Impedance imaging condition[J].Geophysics,1987,52(8):1060.

[5] SCHLEI CHER J,COSTA J C.A comparison of imaging conditions for wave-equation shot-profile migration[J].Geophysics,2008,73(6):S219.

[6] 董良國,馬在田,曹景忠,等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法[J].地球物理學(xué)報(bào),2000,43(3):411.

[7] BERFENGER J P.A perfectly matched layer for the absorption of electro magnetics waves[J].Journal Computation Physics,1994,114(2):185.

[8] CLAERBOUT J F.Toward a unified theory of reflector mapping[J].Geophysics,1971,36(3):467.

[9] FLETCHER R P,FOWLER P J,KITCHENS IDE P,et al.Suppressing unwanted internal reflections in prestack reverse-time migration[J].Geophysics,2006,71(6):E79.

[10]YOON K,MARFURT K J.Revsrse-time migration using the pointing vector[J].Exploration Geophysics,2006,59(1):102.

[11]GUITTON A,KAELIN B.Least-square attenuation of reverse time migration artifacts[C].Expanded Abstracts of 76rdSEG Annual Int Mtg,2006:2348.

[12]GUITTON A,VALENCIANO A,BEVC D,et al.Smoothing imaging condition for shot-profile migration[J].Geophysics,2007,72(3):S149.