雙程聲波方程逆時深度偏移

何兵壽,張會星,范國苗

(中國海洋大學 海底科學與探測技術教育部重點實驗室,山東青島 266100)

0 前言

地震波偏移是地震勘探領域的核心技術,工業界對地震波場疊前偏移的要求主要有二方面:

(1)復雜構造的精確成像。

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

在理論上,基于雙程波方程的地震資料,疊前逆時深度偏移技術能夠同時達到以上目標,因此該項技術在提出之初就展現出良好的工業應用前景。近年來,隨著油氣生產對地震勘探精度要求的不斷提高,勘探目標構造與巖性復雜程度的不斷增加,業界開始逐漸研究并采用實用的地震波疊前逆時深度偏移技術來解決問題,以達到提高勘探精度與實現巖性勘探等目標。

多年來,國內、外地球物理學者在逆時偏移領域進行了大量的研究工作,并取得了許多成果[1～5]。但目前為止,基于雙程波方程的逆時偏移技術在工業生產資料的實際應用方面還少有成功實例,限制該項技術進入工業應用領域的主要原因在于:

(1)地震波方程逆時延拓本身還有許多技術問題需要進一步完善,如計算精度的進一步提高,截斷邊界偽反射的壓制等。

(2)雙程波方程逆時偏移的成像準則與成像條件的計算方法需進一步深入研究,在地震波方程逆時延拓過程中,地下各點在何時成像,如何成像是決定偏移效果的又一關鍵因素。而目前常用的,基于時間一致性準則的激發時間成像條件,無法補償地震波逆時延拓過程中的層間能量反射損失,難以達到理想的偏移效果。

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

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

1 雙程聲波方程的逆時延拓算法

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

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

其中 vx、vz分別為質點在x與z方向的振動速度,p為位移;v為縱波速度;ρ為密度;x、z分別為空間坐標;t為時間。

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

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

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

圖1 交錯網格示意圖Fig.1 The diagram staggered grids

1.2 穩定性條件與吸收邊界條件

經推導,式(2)～式(4)的穩定性條件為:

采用最佳匹配層(PML)吸收邊界條件,解決截斷邊界的偽反射問題。依據PML的方程分裂思路[7],可得:

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

同樣在交錯網格空間中對式(4)進行差分離散,可得到適用于雙程聲波方程逆時偏移的PML邊界差分格式,詳細推導過程這里不再贅述。

2 雙程聲波方程逆時偏移的成像條件

2.1 疊后逆時偏移的成像條件

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

2.2 疊前逆時偏移的成像條件

在縱波方程的疊前逆時深度偏移領域,Claerbout[8]提出的上行、下行波互相關成像條件可概括為:其中 I mage(x,z)為對應網格點的偏移成像結果;為接收點記錄逆時延拓波場值;x、z為空間坐標;tmax為記錄長度;為炮點子波正向延拓波場值。

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

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

3 雙程聲波方程逆時偏移的脈沖響應

3.1 均勻介質中的脈沖響應

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

圖2 脈沖響應試驗數據Fig.2 seismic data for impulse response test

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

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

作者在本文采用式(9)作為逆時偏移的成像條件,以部份補償地震波在逆時延拓過程中發生的第三次能量損失。圖3(b)為基于歸一化互相關成像條件得到的脈沖響應。分析對比圖3(a)與圖3(b),可以得以下結論:

(1)在二種成像條件下,各波組內部各道的能量關系未發生變化,這表明式(9)未改變逆時偏移中各成像點地震波動力學特征的相對關系。

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

3.2 層狀介質中的脈沖響應

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

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

3.3 連續速度介質中的脈沖響應

連續介質的速度為:

式中 z為深度。

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

圖5 連續速度介質中雙程聲波方程逆時偏移的脈沖響應Fig.5 Reverse-time migration impulse response of acoustic equation in continuous medium

4 模型算例

4.1 疊后逆時深度偏移模型算例

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

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

圖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的疊后逆時深度偏移結果Fig.8 Post stack reverse-time migration result of Fig.7

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

4.2 疊前逆時深度偏移算例

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

圖9(見下頁)為各單炮記錄疊前逆時深度偏移后的共成像點疊加剖面,高速體下部地層的成像質量得到的明顯改善,說明基于雙程聲波方程的疊前逆時深度偏移算法,能夠對復雜構造進行成像。

5 結論及展望

波動方程逆時偏移技術具有不受傾角限制與高保真性等優點,是解決復雜構造成像問題的理想工具。作者在本文推導的基于交錯網格的雙程聲波方程高階逆時延拓算子,與相應的PML吸收邊界條件可實現地表波場的高精度逆時延拓。并且歸一化互相關成像條件,可有效補償地震波逆時延拓過程中的能量損失。綜合應用這二項技術,可實現復雜構造的精確成像。

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

圖9 Sigsbee_2b模型的疊前逆時深度偏移結果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] 董良國,馬在田,曹景忠,等.一階彈性波方程交錯網格高階差分解法[J].地球物理學報,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.