起伏地表棱柱波逆時(shí)偏移成像方法

冉 崎,陳 康,韓 嵩,馬 博,劉 暢,曲英銘

(1.中國(guó)石油天然氣股份有限公司西南油氣田分公司勘探開(kāi)發(fā)研究院,四川成都610021;2.中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院地球物理系,山東青島266580)

山前帶地區(qū)地表起伏劇烈,復(fù)雜的高陡構(gòu)造發(fā)育。劇烈變化的起伏地表給常規(guī)地震勘探工作帶來(lái)了極大的挑戰(zhàn)。常規(guī)起伏地表地震資料處理方法包括靜校正、基準(zhǔn)面校正以及直接成像法[1-4],但當(dāng)?shù)乇砜v橫向速度變化劇烈時(shí),利用常規(guī)靜校正技術(shù)和基準(zhǔn)面校正方法難以得到準(zhǔn)確的成像結(jié)果,因此,起伏地表疊前深度偏移方法得到了快速發(fā)展。

起伏地表深度偏移成像方法包括射線類、單程波動(dòng)方程類和雙程波動(dòng)方程類偏移成像三大類方法。射線類偏移方法主要包括Kirchhoff偏移[5-6]和束偏移[7-13]方法。該類偏移方法運(yùn)算效率高且對(duì)于起伏地表具有較好的適應(yīng)性。WIGGINS[5]提出了面向起伏地表地震數(shù)據(jù)的Kirchhoff偏移方法,在此基礎(chǔ)上,JAGER等[6]提出了保幅Kirchhoff起伏地表偏移方法。秦寧等[7]針對(duì)復(fù)雜山前帶成像難題提出了一種高斯束疊前深度偏移技術(shù)。GRAY[8]將束偏移方法引入Kirchhoff偏移方法,既保留了高計(jì)算效率,又提升了偏移剖面的信噪比。岳玉波等[9]基于平面波分解,提出面向起伏地表地震數(shù)據(jù)的保幅高斯束成像方法。黃建平等[10]將岳玉波等[9]提出的面向起伏地表地震數(shù)據(jù)的保幅高斯束成像方法推廣到了彈性介質(zhì)。楊繼東等[11]提出了一種適用于陸地復(fù)雜地表?xiàng)l件的疊前菲涅爾束偏移方法。徐秀剛等[12]基于“逐步累加”的“直接下延”法,將小束源與Fourier傳播算子相結(jié)合,提出了基于起伏地表的小波束疊前深度偏移成像技術(shù)。黃建平等[13]基于有效鄰域波場(chǎng)近似理論提出了一種成像精度更高且適用于復(fù)雜起伏地表?xiàng)l件的疊前保幅高斯束偏移方法。雖然射線理論憑借較高的運(yùn)算效率,在工業(yè)界已經(jīng)得到了廣泛應(yīng)用,但受方法本身的限制,射線類偏移方法只能反映運(yùn)動(dòng)學(xué)特征,不能反映動(dòng)力學(xué)特征,且無(wú)法刻畫(huà)高陡構(gòu)造等復(fù)雜地質(zhì)體。因此,基于單程波動(dòng)方程的疊前深度偏移方法得到了快速發(fā)展。陳林謙[14]推導(dǎo)并分析了帶誤差補(bǔ)償?shù)念l率空間域有限差分法單程波延拓方法,在此基礎(chǔ)上,提出了帶誤差補(bǔ)償?shù)钠鸱乇懑B前深度偏移技術(shù)流程。段心標(biāo)等[15]提出了基于GPU的三維起伏地表裂步傅里葉疊前深度偏移技術(shù)。張宇等[16]提出了真振幅全傾角單程波方程偏移方法。近些年基于單程波方程的起伏地表直接成像方法發(fā)展較為緩慢,因?yàn)榫统上窬榷?單程波成像方法不如逆時(shí)偏移方法;就計(jì)算效率而言,單程波成像方法不如射線類方法。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,基于雙程波方程的起伏地表疊前逆時(shí)偏移技術(shù)得到了發(fā)展。逆時(shí)偏移成像方法因無(wú)傾角限制,在高陡構(gòu)造成像中發(fā)揮了重要作用。WHITMORE等[17]、BAYSAL等[18]和LOEWENTHAL等[19]均提出了成像效果良好的逆時(shí)偏移成像技術(shù)。劉紅偉等[20]給出了針對(duì)起伏地表直接進(jìn)行疊前逆時(shí)偏移的實(shí)現(xiàn)過(guò)程。王延光等[21]提出利用逆時(shí)偏移實(shí)現(xiàn)基于起伏地表的波動(dòng)方程疊前深度偏移。魯雁翔等[22]推導(dǎo)出了彈性波正向傳播和逆時(shí)傳播的曲線網(wǎng)格差分格式,實(shí)現(xiàn)了起伏層狀介質(zhì)中的波場(chǎng)逆時(shí)偏移。蘭海強(qiáng)等[23]推導(dǎo)出了一種曲線坐標(biāo)系下穩(wěn)定的、顯式二階精度有限差分方法的彈性VTI方程。李慶洋等[24]提出了曲線坐標(biāo)系下基于仿真型有限差分的貼體全交錯(cuò)網(wǎng)格波場(chǎng)延拓算法。

山前帶地區(qū)高陡構(gòu)造的成像難度極大,逆時(shí)偏移成像方法雖可以改善高陡構(gòu)造的成像精度,但因觀測(cè)系統(tǒng)的局限性,高陡構(gòu)造的成像能量依舊較弱。傳統(tǒng)的地震勘探中,通常將棱柱波視為噪聲而將其壓制[25]。事實(shí)上,棱柱波攜帶了許多一次波所不包含的地下高陡構(gòu)造的信息。因此,棱柱波可以被用來(lái)提升照明度和高陡構(gòu)造的成像效果[26-30]。CAVALCA等[31]闡明了棱柱波包含兩個(gè)反射點(diǎn)和3段反射路徑。MARMALYEVSKYY等[32]提出了一種基于克?；舴虻睦庵ǔ上穹椒▉?lái)描述高陡鹽丘側(cè)翼。MALCOLM等[33]采用一種迭代單程波算法逐步對(duì)一次波和多次波成像。BROTO等[34]通過(guò)有限差分正演分析了棱柱波的運(yùn)動(dòng)學(xué)特征并提出了一種鑒別棱柱波的方法。LI等[35]提出了針對(duì)復(fù)雜鹽丘構(gòu)造的棱柱波成像流程。DAI[36]提出了一種無(wú)需提前拾取反射界面且對(duì)初始速度模型的精度要求較低的棱柱波逆時(shí)偏移成像方法。MAKSYM等[37]提出了最小二乘棱柱波逆時(shí)偏移方法;YANG等[38]提出了一次波與棱柱波聯(lián)合最小二乘逆時(shí)偏移方法。

本文分別介紹了曲坐標(biāo)系下的聲波方程、棱柱波逆時(shí)偏移方程以及曲坐標(biāo)系棱柱波逆時(shí)偏移方程,提出了適用于山前帶高陡構(gòu)造成像的起伏地表棱柱波逆時(shí)偏移方法,并利用兩個(gè)典型模型分別對(duì)該方法進(jìn)行了測(cè)試,證明了方法的正確性與有效性。

1 方法原理

1.1 曲坐標(biāo)系下的聲波方程

二維聲波介質(zhì)一階應(yīng)力-速度方程表示為:

(1)

式中:σ為應(yīng)力;vx和vz分別為質(zhì)點(diǎn)在x方向和z方向的速度分量;ρ為密度;vP為縱波傳播速度。

本文采用蔣麗麗等[39]提出的正交貼體網(wǎng)格技術(shù)對(duì)速度模型進(jìn)行剖分,通過(guò)調(diào)節(jié)網(wǎng)格線與邊界的間距、網(wǎng)格線與邊界的夾角,縮小了實(shí)際值與期望值的差距?；诓此煞匠痰馁N體網(wǎng)格表示為:

(2)

式中:P(x,z),Q(x,z)為不同的“源項(xiàng)”,其作用是調(diào)節(jié)網(wǎng)格間距與網(wǎng)格形狀;x,z為原始笛卡爾坐標(biāo)系下不同網(wǎng)格節(jié)點(diǎn)的坐標(biāo)值;ξ,η為轉(zhuǎn)換后的曲坐標(biāo)系下的不同網(wǎng)格節(jié)點(diǎn)的坐標(biāo)值;下角標(biāo)代表求二階導(dǎo)數(shù)。

物理域(笛卡爾坐標(biāo)系)與計(jì)算域(曲坐標(biāo)系)關(guān)系表示為:

(3)

(3)式兩邊分別對(duì)x和z求導(dǎo),表示為:

(4)

(5)

整理(4)式和(5)式得到的坐標(biāo)變換關(guān)系表示為:

(6)

式中:J為雅可比行列式。

聯(lián)立(1)式和(6)式并結(jié)合鏈?zhǔn)椒▌t得到曲坐標(biāo)系下的波動(dòng)方程為:

(7)

1.2 棱柱波逆時(shí)偏移

棱柱波作為一種多次波,攜帶了更多的地質(zhì)體信息,特別是來(lái)自地下高陡構(gòu)造的信息。在描述高陡構(gòu)造方面,棱柱波逆時(shí)偏移可以發(fā)揮很大的作用。根據(jù)棱柱波的產(chǎn)生機(jī)理可將其分為兩種類型,第一種棱柱波是地震波先在低傾角反射界面上發(fā)生反射,隨后又在高陡體側(cè)翼發(fā)生反射,傳播至地面被檢波器接收;第二種棱柱波是地震波先在高陡體側(cè)翼發(fā)生反射,隨后又在低傾角沉積界面上發(fā)生反射,傳播至地面被檢波器接收。

棱柱波逆時(shí)偏移是針對(duì)高陡構(gòu)造而發(fā)展出的一種多次波成像技術(shù)。該方法包括兩個(gè)部分:①首先輸入偏移速度場(chǎng)以及炮記錄,應(yīng)用常規(guī)逆時(shí)偏移成像得到低傾角反射層位置;②對(duì)第①步得到的成像結(jié)果繼續(xù)應(yīng)用棱柱波逆時(shí)偏移及其特定成像條件,得到棱柱波逆時(shí)偏移結(jié)果。

頻率域共炮域逆時(shí)偏移公式可表示為[24]:

(8)

式中:mmig(x|xs)為逆時(shí)偏移波場(chǎng)值;W*(ω)為震源共軛頻譜;xg為檢波點(diǎn)位置;G*為共軛格林函數(shù);d(xg|xs)為炮記錄;xs為震源點(diǎn)位置。

對(duì)于給定的偏移速度場(chǎng),可以求得震源波場(chǎng)的格林函數(shù)G0,進(jìn)而可求得反射波場(chǎng)的格林函數(shù)G1,即[40]:

(9)

式中:m1(x′)為水平反射層在常規(guī)逆時(shí)偏移處理后得到的反射系數(shù)模型。

由此可得笛卡爾坐標(biāo)系下的頻率域棱柱波逆時(shí)偏移公式,可表示為:

(10)

式中:d1(xg|xs)為一階散射反射波;d2(xg|xs)為二階散射棱柱波;P1,Q0可由如下公式求得:

(11)

式中:ω為頻率;v0(x)為偏移速度;P0(x)為震源波場(chǎng);P1(x)為反射波場(chǎng);W(ω)為震源頻譜;Q0(x)為檢波點(diǎn)波場(chǎng);m1(x)為水平反射層經(jīng)常規(guī)逆時(shí)偏移處理后得到的反射系數(shù)模型。

1.3 曲坐標(biāo)系棱柱波逆時(shí)偏移

(7)式還可以寫(xiě)成如下形式,即:

(12)

速度擾動(dòng)可表示為:

vP=vP0+δvP

(13)

式中:vP0為背景速度;δvP為擾動(dòng)速度。

應(yīng)力場(chǎng)p0可通過(guò)下式求解:

(14)

由泰勒公式可知,背景速度與擾動(dòng)速度的關(guān)系可表示為:

(15)

將(15)式代入(12)式中,并減去(14)式,可得:

(16)

反射系數(shù)模型可表示為:

(17)

將(17)式代入(16)式,可得:

(18)

本文中,擾動(dòng)場(chǎng)(δp,δvx,δvz)為反射波場(chǎng)(p1,vx1,vz1)。采用常規(guī)逆時(shí)偏移成像結(jié)果I代替m(vP),(18)式可表示為:

(19)

曲坐標(biāo)系棱柱波逆時(shí)偏移的實(shí)現(xiàn)主要包括如下9個(gè)步驟:①輸入笛卡爾坐標(biāo)系下的地震炮記錄和偏移速度模型;②計(jì)算曲網(wǎng)格并將偏移速度轉(zhuǎn)換到曲坐標(biāo)系下;③計(jì)算正向延拓的震源波場(chǎng)和逆時(shí)延拓的檢波點(diǎn)波場(chǎng);④得到曲坐標(biāo)系下的常規(guī)逆時(shí)偏移;⑤在步驟③的正向延拓震源波場(chǎng)的基礎(chǔ)上,再一次進(jìn)行波場(chǎng)正傳計(jì)算以正向延拓反射波場(chǎng);⑥在步驟③的逆時(shí)延拓檢波點(diǎn)波場(chǎng)的基礎(chǔ)上,再一次進(jìn)行波場(chǎng)逆?zhèn)饔?jì)算逆時(shí)延拓的反射波場(chǎng);⑦得到曲坐標(biāo)系下的棱柱波逆時(shí)偏移成像結(jié)果;⑧將棱柱波逆時(shí)偏移的成像結(jié)果反變換到笛卡爾坐標(biāo)系下;⑨輸出棱柱波成像結(jié)果。

2 模型試算

2.1 簡(jiǎn)單起伏地表凹陷模型

笛卡爾坐標(biāo)系下的簡(jiǎn)單起伏地表凹陷速度模型如圖1a所示,3層的速度分別為3000m/s,4000m/s,4500m/s,曲坐標(biāo)系下的速度模型如圖1b所示。該速度模型x方向與z方向均為2400m,網(wǎng)格數(shù)目均為301,網(wǎng)格間距均為8m,笛卡爾坐標(biāo)系下的貼體網(wǎng)格模型如圖2a所示,圖2b為其局部放大圖。

圖1 簡(jiǎn)單起伏地表凹陷速度模型a 笛卡爾坐標(biāo)系;b 曲坐標(biāo)系

圖2 貼體網(wǎng)格模型(a)及其局部放大結(jié)果(b)

采用主頻為25Hz的雷克子波,在地表布置30個(gè)均勻分布的震源,炮間距為80m,每炮由301道接收,道間距為8m,總記錄時(shí)間為2.4s,每道采樣間隔為0.6ms。第15炮的單炮記錄如圖3所示。

對(duì)圖3所示的地震數(shù)據(jù)分別進(jìn)行常規(guī)逆時(shí)偏移和棱柱波逆時(shí)偏移處理。圖4為540ms時(shí)正向延拓的波場(chǎng)快照,圖4a和圖4c分別為540ms的曲坐標(biāo)系下正向延拓的p0和p1波場(chǎng)快照,笛卡爾坐標(biāo)系下的波場(chǎng)快照分別如圖4b和圖4d所示,可以看出,p0波場(chǎng)快照中來(lái)自水平層的地震波能量(圖4a和圖4b 空心箭頭所示)比來(lái)自高陡構(gòu)造的地震波能量(圖4a和圖4b實(shí)線箭頭處)強(qiáng)得多;而p1波場(chǎng)快照中兩者的能量更均衡(圖4c和圖4d)。由此可知,在棱柱波逆時(shí)偏移過(guò)程中應(yīng)用p1信息進(jìn)行成像,相比于傳統(tǒng)逆時(shí)偏移方法應(yīng)用p1信息進(jìn)行成像,高陡構(gòu)造的能量分布更均衡。圖5a和圖5b分別為常規(guī)逆時(shí)偏移和棱柱波逆時(shí)偏移的成像結(jié)果,從圖5中可以看出,采用常規(guī)逆時(shí)偏移方法得到的高陡構(gòu)造成像效果較差,采用棱柱波逆時(shí)偏移方法得到的水平構(gòu)造成像效果雖然不及常規(guī)逆時(shí)偏移方法效果,但高陡構(gòu)造的成像效果得到了明顯改善(圖5白色箭頭處)。

圖3 用于常規(guī)逆時(shí)偏移和棱柱波逆時(shí)偏移的單炮記錄

圖4 540ms時(shí)正向延拓的波場(chǎng)快照a 曲坐標(biāo)系下的p0波場(chǎng);b 笛卡爾坐標(biāo)系下的p0波場(chǎng);c 曲坐標(biāo)系下的p1波場(chǎng);d 笛卡爾坐標(biāo)系下的p1波場(chǎng)

圖5 簡(jiǎn)單起伏地表凹陷模型偏移結(jié)果a 常規(guī)逆時(shí)偏移;b 棱柱波逆時(shí)偏移

2.2 起伏地表火山巖模型

本文根據(jù)川西北部龍門(mén)山前帶地區(qū)實(shí)際模型建立了起伏地表火山巖速度模型。該模型中廣泛發(fā)育具有高陡傾角的火山通道構(gòu)造。笛卡爾坐標(biāo)系下的起伏地表火山巖速度模型如圖6a所示,速度范圍為2800～6400m/s,曲坐標(biāo)系下的速度模型如圖6b所示。該速度模型x方向?yàn)?0km,z方向?yàn)?km,x方向網(wǎng)格點(diǎn)數(shù)為1001,z方向網(wǎng)格點(diǎn)數(shù)為501,網(wǎng)格間距均為10m。采用主頻為25Hz的雷克子波,在地表布置125個(gè)均勻分布的震源,炮間距為80m,每炮由1251道接收,總記錄時(shí)間為4.4s,每道采樣間隔為0.5ms。圖7a和圖7b分別為采用常規(guī)逆時(shí)偏移和采用棱柱波逆時(shí)偏移得到的成像結(jié)果。從圖7a中可以看出,起伏地表?xiàng)l件下常規(guī)逆時(shí)偏移可對(duì)低傾角的火山巖體及沉積層清晰成像,但對(duì)于地震、地質(zhì)解釋非常重要的火山通道構(gòu)造并沒(méi)有得到較好的刻畫(huà)。對(duì)比圖7a和圖7b發(fā)現(xiàn),起伏地表棱柱波逆時(shí)偏移方法在成像時(shí)增強(qiáng)了棱柱波的能量,因此可以對(duì)含高陡傾角的火山巖通道清晰成像,為地震、地質(zhì)解釋提供了直觀的可靠依據(jù)。

圖6 起伏地表火山巖速度模型a 笛卡爾坐標(biāo)系;b 曲坐標(biāo)系

圖7 起伏地表火山巖模型偏移結(jié)果a 常規(guī)逆時(shí)偏移;b 棱柱波逆時(shí)偏移

3 討論

利用本文方法進(jìn)行實(shí)際數(shù)據(jù)處理時(shí),不需要進(jìn)行高程靜校正,只需根據(jù)地表高程函數(shù)生成貼體網(wǎng)格,本文方法與常規(guī)逆時(shí)偏移一樣都需要較為準(zhǔn)確的偏移速度場(chǎng)。在應(yīng)用本文方法進(jìn)行偏移前,首先采用常規(guī)偏移方法得到低傾角沉積層成像結(jié)果,輸入的數(shù)據(jù)包括炮域地震數(shù)據(jù)、偏移速度場(chǎng)、貼體網(wǎng)格文件、觀測(cè)系統(tǒng)文件和常規(guī)成像結(jié)果。雖然本文方法能夠改善高陡構(gòu)造的成像結(jié)果,但損失了低傾角沉積層的成像效果,因此,在地震、地質(zhì)解釋時(shí),需要與傳統(tǒng)成像方法的成像結(jié)果進(jìn)行聯(lián)合分析。

4 結(jié)論與認(rèn)識(shí)

地表起伏劇烈、高陡構(gòu)造廣泛發(fā)育導(dǎo)致山前帶地區(qū)成像困難,針對(duì)這一難題,將基于貼體網(wǎng)格的曲坐標(biāo)系聲波方程與棱柱波逆時(shí)偏移技術(shù)結(jié)合,提出并實(shí)現(xiàn)了基于貼體網(wǎng)格剖分的曲坐標(biāo)系棱柱波逆時(shí)偏移成像方法,該方法克服了山前帶地表劇烈起伏對(duì)地震成像的影響,同時(shí)充分利用包含大量高陡構(gòu)造信息的棱柱波進(jìn)行成像,改善了山前帶高陡構(gòu)造的成像效果。