徐蔚亞 朱成宏 魏哲楓 張曉語 杜啟振

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

0 引言

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

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

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

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

1 方法原理

1.1 解耦延拓方程

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

(1)

縱波方程為

(2)

橫波波場可表示為

(3)

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

總彈性波場應力分量減去縱波體應力τP，可以構建橫波應力分量[15，17]

(4)

式中τS,xx、τS,zz和τS,xz為橫波應力分量。對上式進行時間積分，可得由彈性應變表示的橫波應力

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

(5)

式中：εij為彈性應變分量；b=εxx+εzz，為彈性體應變。由于彈性應變既包含縱波應變εP,ij，又包含橫波應變εS,ij，則有

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

(6)

同時橫波應力可以改寫為

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

(7)

式中：bP=εP,xx+εP,zz，為縱波體應變；bS=εS,xx+εS,zz，為橫波體應變?？梢姍M波應力除了橫波應變作用2μ(εS,ij-bSδij)外，還有縱波應變作用2μ(εP,ij-bPδij)，因此，橫波應力中含有縱波串擾。

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

(8)

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

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

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

S=-v·τ

(9)

式中：v=(vx，vz)為質(zhì)點振動速度矢量；τ為應力張量。

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

圖1 縱波Poynting矢量與反射角的關系示意圖

(10)

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

圖2 縱波入射角β與橫波反射角α關系示意圖

(11)

式中

(12)

由VP、VS和開角θPS就可以確定PS波反射角α。

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

(13)

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

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

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

IPP_ADCIG(x,βk)

(14)

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

(1)在炮域構建震源波場，計算每一時刻的縱波質(zhì)點振動速度場，并將其存儲至硬盤；

(2)在炮域構建檢波點波場，計算每一時刻的縱波質(zhì)點振動速度場及應力場；

(3)從硬盤讀取縱波質(zhì)點振動速度場，利用式(9)計算當前時刻的Poynting矢量，并利用式(10)計算縱波的開角和入射角；

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

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

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

IPS_ADCIG(x,βk)

(15)

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

(1)在炮域構建震源波場，計算每一時刻的縱波質(zhì)點振動速度場，并將其存儲至硬盤；

(2)在炮域構建檢波點波場，計算每一時刻的橫波質(zhì)點振動速度場及應力場；

(3)從硬盤讀取縱波質(zhì)點振動速度場，利用式(9)計算當前時刻的Poynting矢量，并求取縱波入射及橫波反射之間的開角θPS；

(4)結合Snell定律及波場傳播關系，利用式(11)計算PS波反射角；

(5)應用式(13)進行角度平滑；

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

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

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

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

(16)

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

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

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

2 模型驗證

2.1 水平層狀模型

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

圖3 水平層狀模型

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

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

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

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

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

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

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

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

2.2 Salt模型

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

圖8 Salt縱波速度模型

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

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

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

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

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

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

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

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

3 結論和認識

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

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

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

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