徐蔚亞 朱成宏 魏哲楓 張曉語(yǔ) 杜啟振

(①頁(yè)巖油氣富集機(jī)理與有效開(kāi)發(fā)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100083；②中國(guó)石化彈性波理論與探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100083；③中國(guó)石化石油勘探開(kāi)發(fā)研究院，北京 100083；④齊魯工業(yè)大學(xué)(山東省科學(xué)院)山東省科學(xué)院海洋儀器儀表研究院，山東青島 266061；⑤中國(guó)石油大學(xué)(華東)深層油氣重點(diǎn)實(shí)驗(yàn)室，山東青島 266580)

0 引言

逆時(shí)偏移(RTM)方法[1-4]于20世紀(jì)80年代提出，該方法基于波動(dòng)方程實(shí)現(xiàn)波場(chǎng)延拓，以其無(wú)傾角限制、能對(duì)多種波(反射波、多次波、棱柱波)成像而受到廣泛關(guān)注[5-6]。彈性波逆時(shí)偏移(ERTM)算法[7-8]由McMechan研究組提出，采用更接近真實(shí)的彈性介質(zhì)和彈性波場(chǎng)，與傳統(tǒng)的聲波RTM成像方法相比，理論假設(shè)與實(shí)際更接近，并能同時(shí)實(shí)現(xiàn)縱波和橫波的深度域成像。

地震波波場(chǎng)分離是ERTM方法中的重要環(huán)節(jié)，也是得到純波模式成像結(jié)果的關(guān)鍵一步。除了基于Helmholtz定理的波場(chǎng)分離方法[9-12]以外，目前常用的是解耦延拓方程法[13-16]。解耦延拓方程法可在波場(chǎng)延拓過(guò)程中實(shí)現(xiàn)均勻介質(zhì)的縱、橫波質(zhì)點(diǎn)位移的解耦或者質(zhì)點(diǎn)振動(dòng)速度的矢量解耦，但在應(yīng)力解耦時(shí)，橫波應(yīng)力存在能量串?dāng)_問(wèn)題[17-19]。為解決橫波應(yīng)力中的縱波串?dāng)_問(wèn)題，Du等[16-17]提出了偽橫波應(yīng)力構(gòu)建方法，實(shí)現(xiàn)了縱波和橫波應(yīng)力的解耦。

共成像點(diǎn)道集(CIG)作為疊前深度偏移的重要輸出道集，可以用于改善成像質(zhì)量，提供反映地下巖性的振幅和相位等信息，同時(shí)能用于速度建模。CIG[20]主要分為共炮檢距道集(Common-offset CIG，COCIG)和角度域道集(Angle-domain CIG，ADCIG)。Sava等[21]利用單程波偏移提取地下局部COCIG，通過(guò)轉(zhuǎn)換獲得ADCIG。三維情況下，Xu等[22-23]在頻率—波數(shù)域通過(guò)傳播方向求取反射角信息和方位角信息，利用合適的成像條件提取角道集，優(yōu)點(diǎn)是角道集質(zhì)量較好，缺點(diǎn)是計(jì)算量和存儲(chǔ)量均較大。王保利等[24]在實(shí)現(xiàn)過(guò)程中利用波場(chǎng)延拓獲得的質(zhì)點(diǎn)振動(dòng)速度和應(yīng)力，提取了成像點(diǎn)處的Poynting矢量，從而確定入射波和反射波的傳播方向，最后通過(guò)構(gòu)建成像條件提取了ADCIG，該方法簡(jiǎn)單易行；汪天池等[25]將該方法拓展到傾角域進(jìn)行反射波與繞射波分離以實(shí)現(xiàn)繞射波RTM。Zhao等[26]通過(guò)提取地下COCIG進(jìn)而提取了PP波和PS波的ADCIG，并分析了二者的運(yùn)動(dòng)學(xué)特征。Wang等[27]通過(guò)構(gòu)建Poynting矢量提取了PP波和PS波的ADCIG，但該方法沒(méi)有考慮分離得到的橫波應(yīng)力存在的縱波能量串?dāng)_問(wèn)題。

針對(duì)橫波應(yīng)力串?dāng)_問(wèn)題，本文基于解耦延拓方程及偽橫波應(yīng)力構(gòu)建方法，開(kāi)展了基于Poynting矢量提取橫波反射角的研究，構(gòu)建了基于ERTM的PP波和PS波ADCIG求取方法；針對(duì)PP波ADCIG不同角度的成像特性，提出了把ADCIG按角度衰減的疊加成像策略。最后應(yīng)用模型數(shù)據(jù)驗(yàn)證了本文方法的正確性。

1 方法原理

1.1 解耦延拓方程

ERTM是利用彈性波動(dòng)方程進(jìn)行波場(chǎng)延拓，并結(jié)合波場(chǎng)分離方法得到解耦的震源端和檢波端的縱/橫波波場(chǎng)，然后運(yùn)用彈性波成像條件得到最終的PP波和PS波成像結(jié)果。解耦延拓方程[14-19]能實(shí)現(xiàn)彈性波場(chǎng)的構(gòu)建和解耦。其中，一階彈性波波場(chǎng)方程[18]為

(1)

縱波方程為

(2)

橫波波場(chǎng)可表示為

(3)

式中：vx、vz分別為總彈性波場(chǎng)中x和z方向的質(zhì)點(diǎn)振動(dòng)速度分量；vP,x、vP,z為縱波波場(chǎng)的兩個(gè)方向速度分量；vS,x、vS,z為橫波波場(chǎng)的兩個(gè)方向速度分量；τxx、τzz、τxz為總彈性波場(chǎng)應(yīng)力分量；τP為縱波體應(yīng)力；λ、μ和ρ分別為彈性介質(zhì)的拉梅常數(shù)和密度。

總彈性波場(chǎng)應(yīng)力分量減去縱波體應(yīng)力τP，可以構(gòu)建橫波應(yīng)力分量[15，17]

(4)

式中τS,xx、τS,zz和τS,xz為橫波應(yīng)力分量。對(duì)上式進(jìn)行時(shí)間積分，可得由彈性應(yīng)變表示的橫波應(yīng)力

τS,ij=2μ(εij-bδij)

(5)

式中：εij為彈性應(yīng)變分量；b=εxx+εzz，為彈性體應(yīng)變。由于彈性應(yīng)變既包含縱波應(yīng)變?chǔ)臥,ij，又包含橫波應(yīng)變?chǔ)臩,ij，則有

εij=εP,ij+εS,ij

(6)

同時(shí)橫波應(yīng)力可以改寫為

τS,ij=2μ[(εP,ij+εS,ij)-(bP+bS)δij]

(7)

式中：bP=εP,xx+εP,zz，為縱波體應(yīng)變；bS=εS,xx+εS,zz，為橫波體應(yīng)變。可見(jiàn)橫波應(yīng)力除了橫波應(yīng)變作用2μ(εS,ij-bSδij)外，還有縱波應(yīng)變作用2μ(εP,ij-bPδij)，因此，橫波應(yīng)力中含有縱波串?dāng)_。

為壓制縱波能量在橫波應(yīng)力中的影響，Du等[16-17]基于橫波質(zhì)點(diǎn)振動(dòng)速度提出了偽橫波應(yīng)力構(gòu)建方程

(8)

式中：τqS,xx和τqS,zz分別為x方向和z方向的偽橫波正應(yīng)力分量；τqS,xz為偽橫波切應(yīng)力?；跈M波質(zhì)點(diǎn)振動(dòng)速度構(gòu)建的偽橫波應(yīng)力與橫波應(yīng)力在形式上保持一致。不同于橫波應(yīng)力，偽橫波應(yīng)力只與橫波質(zhì)點(diǎn)振動(dòng)速度相關(guān)，避免了縱波應(yīng)變對(duì)橫波應(yīng)力的影響。由于橫波質(zhì)點(diǎn)振動(dòng)速度滿足?vS,z/?z+?vS,x/?x=0，則偽橫波應(yīng)力滿足τqS,xx=-τqS,zz。因此，基于解耦延拓方程和偽橫波應(yīng)力構(gòu)建方法可以實(shí)現(xiàn)縱橫波應(yīng)力場(chǎng)的解耦，進(jìn)而可以實(shí)現(xiàn)彈性應(yīng)力和質(zhì)點(diǎn)振動(dòng)速度的解耦。

1.2 基于Poynting矢量提取PP波和PS波反射角

Poynting將能流密度矢量用于描述電磁能傳播方向，即通過(guò)電場(chǎng)強(qiáng)度與磁場(chǎng)強(qiáng)度點(diǎn)乘得到電磁波的傳播方向，后人把能流密度矢量稱為Poynting矢量。Yoon等[28]給出了地震波場(chǎng)Poynting矢量表達(dá)式

S=-v·τ

(9)

式中：v=(vx，vz)為質(zhì)點(diǎn)振動(dòng)速度矢量；τ為應(yīng)力張量。

利用解耦的質(zhì)點(diǎn)振動(dòng)速度及應(yīng)力，可以確定縱波和橫波的Poynting矢量，進(jìn)而確定縱波和橫波能量傳播方向。

圖1 縱波Poynting矢量與反射角的關(guān)系示意圖

(10)

由于縱波反射角等于入射角，因此開(kāi)角的一半即為縱波反射角。根據(jù)Snell定律，縱波入射角β和橫波反射角α之間滿足sinβ/VP=sinα/VS，其中VP、VS為縱、橫波傳播速度。

圖2 縱波入射角β與橫波反射角α關(guān)系示意圖

(11)

式中

(12)

由VP、VS和開(kāi)角θPS就可以確定PS波反射角α。

由于Poynting矢量的非穩(wěn)健性，往往采用在一定空間窗上進(jìn)行平滑的方式求取θPP或θPS，即

(13)

式中Ω表示平滑的空間窗口。

1.3 PP波和PS波角度域道集提取方法

獲得縱波入射角β和橫波反射角α后，便可以提取ADCIG以對(duì)成像過(guò)程進(jìn)行質(zhì)控，從而進(jìn)一步提高成像精度。常規(guī)成像條件是對(duì)角道集所有角度的疊加，故不能直接用于提取ADCIG。為此，需對(duì)常規(guī)成像條件進(jìn)行修改。借鑒王保利等[24]的方法，引入高斯采樣函數(shù)，將常規(guī)成像條件拓展為PP波角度域共成像點(diǎn)成像條件

IPP_ADCIG(x,βk)

(14)

利用成像域PP波成像條件可以得到PP波的ADCIG，其中ERTM中PP波ADCIG的計(jì)算流程如下：

(1)在炮域構(gòu)建震源波場(chǎng)，計(jì)算每一時(shí)刻的縱波質(zhì)點(diǎn)振動(dòng)速度場(chǎng)，并將其存儲(chǔ)至硬盤；

(2)在炮域構(gòu)建檢波點(diǎn)波場(chǎng)，計(jì)算每一時(shí)刻的縱波質(zhì)點(diǎn)振動(dòng)速度場(chǎng)及應(yīng)力場(chǎng)；

(3)從硬盤讀取縱波質(zhì)點(diǎn)振動(dòng)速度場(chǎng)，利用式(9)計(jì)算當(dāng)前時(shí)刻的Poynting矢量，并利用式(10)計(jì)算縱波的開(kāi)角和入射角；

(4)利用式(13)進(jìn)行角度平滑；

(5)以離散角作為角道集輸出的角度變量，利用式(14)輸出PP波ADCIG。

類似地，利用得到的橫波反射角α，本文給出了PS波ADCIG成像條件

IPS_ADCIG(x,βk)

(15)

利用成像域PS波成像條件，可以得到PS波的ADCIG。其中ERTM中PS波ADCIG的計(jì)算流程如下：

(1)在炮域構(gòu)建震源波場(chǎng)，計(jì)算每一時(shí)刻的縱波質(zhì)點(diǎn)振動(dòng)速度場(chǎng)，并將其存儲(chǔ)至硬盤；

(2)在炮域構(gòu)建檢波點(diǎn)波場(chǎng)，計(jì)算每一時(shí)刻的橫波質(zhì)點(diǎn)振動(dòng)速度場(chǎng)及應(yīng)力場(chǎng)；

(3)從硬盤讀取縱波質(zhì)點(diǎn)振動(dòng)速度場(chǎng)，利用式(9)計(jì)算當(dāng)前時(shí)刻的Poynting矢量，并求取縱波入射及橫波反射之間的開(kāi)角θPS；

(4)結(jié)合Snell定律及波場(chǎng)傳播關(guān)系，利用式(11)計(jì)算PS波反射角；

(5)應(yīng)用式(13)進(jìn)行角度平滑；

(6)以離散角作為角道集輸出的角度變量，利用式(15)輸出PS波ADCIG。

1.4 PP波ADCIG角度衰減+Laplace濾波的組合疊加成像策略

受雙程波動(dòng)方程影響，ERTM方法進(jìn)行PP波成像過(guò)程中引入了低波數(shù)噪聲。低波數(shù)噪聲反射角接近90°，主要存在于角道集的大角度區(qū)域，因此對(duì)角道集進(jìn)行疊加時(shí)，可通過(guò)對(duì)大角度切除獲得消除低波數(shù)噪聲的RTM結(jié)果。但是，道集中大角度道也含有部分有效信息。

借鑒Laplace濾波的物理機(jī)制，本文提出了在PP波角道集中小角度按縱波入射角度余弦cosβk的方式進(jìn)行疊加成像、大角度實(shí)施Laplace濾波的組合疊加成像的策略，即

(16)

式中：Nβ2為離散角度個(gè)數(shù)；Nβ1為小角度與大角度的分界點(diǎn)序號(hào)。分界角取值范圍一般為30°～40°。

根據(jù)式(16)可知，本文的ADCIG疊加成像的策略利用了角道集中所有的角度成像，小角度cosβk值更大，因此最終成像結(jié)果中既突出了角度道集中小角度道的貢獻(xiàn)，又避免了大角度道有效信息的損失。

需要指出的是，對(duì)于PS波成像，由于震源端P波與檢波端S波傳播路徑不同，且振動(dòng)方向也不一致，因此PS波成像不存在低波數(shù)噪聲。為此，PS波角道集疊加成像時(shí)可以忽略低波數(shù)噪聲影響，不需要大角度道集噪聲壓制。

2 模型驗(yàn)證

2.1 水平層狀模型

構(gòu)建水平層狀模型驗(yàn)證本文的疊加成像策略的正確性。模型如圖3所示，網(wǎng)格數(shù)為400×350，尺寸為10m×10m，縱橫波速度比為1.73。加載主頻為15Hz的Ricker子波作為爆炸源，道間距為10m，炮間距為100m，在地表全覆蓋，共接收40炮；時(shí)間采樣間隔為1ms，記錄長(zhǎng)度為3.2s。采用時(shí)間二階、空間十階的有限差分算法進(jìn)行數(shù)值模擬。

圖3 水平層狀模型

用式(14)和式(15)分別提取出水平層狀模型的PP波角道集和PS波角道集，如圖4所示。其中，每個(gè)角道集角度范圍為0°～90°，間隔為1°；每隔50個(gè)CMP抽取一個(gè)角道集顯示。從PP波成像道集可以看出，同相軸水平，埋深1500m、2500m界面在ADCIG中的同相軸連續(xù)性較好，成像噪聲較少。由紅色箭頭處可以看出，PP低波數(shù)噪聲主要位于角道集的淺層大角度的位置。從PS波成像角道集可看出，雖然有多次波干擾，但是埋深1500m界面在ADCIG中同相軸是水平的，而且連續(xù)性較好。這證明了PP波角道集和PS波角道集提取方法的正確性。

圖4 水平層狀模型的PP波(a)與PS波角道集(b)

為便于AVO/AVA分析，利用Zoeppritz方程[29]計(jì)算深度1500m界面的PP波和PS波的理論反射系數(shù)曲線。提取深度1500m、CMP200處PP波角道集和PS波角道集振幅曲線并歸一化，并其與理論反射系數(shù)對(duì)比(圖5)顯示，提取的振幅隨角度變化曲線與理論反射系數(shù)曲線基本一致，說(shuō)明本文角道集提取方法基本可行。

圖5 1500m處界面PP波(a)和PS波(b)角道集歸一化振幅曲線與理論反射系數(shù)曲線對(duì)比

分別抽取水平層狀模型PP波角道集中0°～40°、41°～70°及71°～90°的道進(jìn)行疊加成像，分析角度對(duì)成像結(jié)果的影響(圖6)。0°～40°的疊加成像結(jié)果(圖6a)中幾乎沒(méi)有任何低波數(shù)噪聲，41°～70°疊加成像結(jié)果(圖6b)中低波數(shù)噪聲較少，而71°～90°的成像結(jié)果(圖6c)幾乎全部是低頻噪聲，但有部分有效成像信息。

圖6 水平層狀模型PP波不同角度范圍的疊加成像結(jié)果

采用本文方法提出的PP波角度衰減+Laplace濾波的組合疊加成像策略(式(16))，取分界角度為40°，得到組合疊加成像結(jié)果(圖7a)，并抽取CMP200處的0°～40°和40°～90°角道集(圖7b、圖7c)，可見(jiàn)低波數(shù)噪聲得以壓制，且同相軸連續(xù)。

圖7 水平層狀模型PP波組合疊加成像結(jié)果(a)及CMP200處

2.2 Salt模型

利用Salt模型(圖8)測(cè)試本文算法對(duì)復(fù)雜模型的成像精度和效果。模型網(wǎng)格數(shù)為1500×700，網(wǎng)格尺寸為10m×10m，縱橫波速度比為1.73，密度取常數(shù)2.0g/cm3。數(shù)值模擬時(shí)，震源采用主頻為30Hz的Ricker子波以爆炸形式加載至P波應(yīng)力分量上；

圖8 Salt縱波速度模型

時(shí)間采樣間隔為1ms；炮點(diǎn)間隔為100m，共150炮；在地表以10m為間隔均勻分布的檢波點(diǎn)進(jìn)行接收。采用在時(shí)間二階、空間十階的有限差分算法進(jìn)行正向和反向波場(chǎng)延拓。

提取的Salt模型的PP波的ADCIG如圖9a所示，角度范圍為0°～90°，角度間隔為2°，每隔100個(gè)CMP抽取一個(gè)角道集進(jìn)行展示。利用本文提出的組合疊加策略得到的PP波成像結(jié)果如圖9b所示。由圖9a可以看出，每個(gè)ADCIG的同相軸都是水平的，且存在低波數(shù)噪聲。采用組合疊加策略后，成像結(jié)果的低波數(shù)噪聲得到有效壓制。

圖9 Salt模型PP波共成像點(diǎn)道集(a)和成像結(jié)果(b)

用本文提出的基于偽橫波應(yīng)力角道集方法提取的Salt模型PS波ADCIG如圖10a所示。對(duì)0°～90°角道集進(jìn)行疊加得到的PS成像結(jié)果如圖10b所示。由圖可以看出，每個(gè)ADCIG上的同相軸都是水平的，驗(yàn)證了本文PS波角道集提取方法的正確性。在淺層PS波成像結(jié)果優(yōu)于PP波，在深層PP波成像結(jié)果優(yōu)于PS波。

圖10 Salt模型PS波共成像點(diǎn)道集(a)和成像結(jié)果(b)

圖11a和圖11b分別為基于偽橫波應(yīng)力(本文方法)提取的PS波部分CMP角道集及0°～45°疊加成像結(jié)果；圖11c和圖11d分別為基于橫波應(yīng)力提取的PS波部分CMP角道集及0°～45°的疊加成像結(jié)果。圖11c中同相軸不連續(xù)，圖11d中含有低波數(shù)噪聲；圖11a中同相軸較連續(xù)而且圖11b中不存在低波數(shù)噪聲。原因在于圖11c和圖11d中構(gòu)造的橫波應(yīng)力(式(4))中含有縱波成分，導(dǎo)致提取的PS波角道集也含有縱波成分，因此，基于橫波應(yīng)力提取的PS波角道集中存在縱波能量干擾和低波數(shù)噪聲；相比之下，圖11a中基于偽橫波應(yīng)力(式(8))提取的PS波角道集較好壓制了能量串?dāng)_。

圖11 基于偽橫波應(yīng)力的PS波角道集(a)及其0°～45°疊加成像結(jié)果(b)、基于橫波應(yīng)力提取的PS波角道集(c)及其0°～45°疊加成像結(jié)果(d)

此外，對(duì)比PS波0°～45°角道集疊加結(jié)果(圖11b)與0°～90°全部角道集疊加的PS波成像(圖10b)結(jié)果發(fā)現(xiàn)，后者更優(yōu)。

3 結(jié)論和認(rèn)識(shí)

(1)基于解耦延拓方程及偽橫波應(yīng)力實(shí)現(xiàn)的彈性波逆時(shí)偏移方法，可以構(gòu)建得到完全解耦的質(zhì)點(diǎn)振動(dòng)速度和應(yīng)力；

(2)根據(jù)余弦定理可以求得震源波場(chǎng)的Poyn-ting矢量與檢波點(diǎn)波場(chǎng)的Poynting矢量的夾角，進(jìn)而由縱、橫波背景速度和開(kāi)角確定PS波反射角；

(3)通過(guò)引入時(shí)間窗、空間窗以及高斯函數(shù)構(gòu)建的PP波和PS波共成像點(diǎn)道集成像條件，可穩(wěn)健地實(shí)現(xiàn)PP波和PS波角道集的提??；

(4)模型試驗(yàn)結(jié)果表明，組合疊加成像策略構(gòu)建的PP波角道集可以保護(hù)部分大角度有效信息，并實(shí)現(xiàn)低波數(shù)噪聲的有效壓制。