彈性波逆時偏移的角度域共成像點道集提取方法

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

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

0 引言

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

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

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

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

1 方法原理

1.1 解耦延拓方程

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

(1)

縱波方程為

(2)

橫波波場可表示為

(3)

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

總彈性波場應力分量減去縱波體應力τ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，為橫波體應變。可見橫波應力除了橫波應變作用2μ(εS,ij-bSδij)外，還有縱波應變作用2μ(εP,ij-bPδij)，因此，橫波應力中含有縱波串擾。

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

(8)

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

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

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

S=-v·τ

(9)

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

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

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

(10)

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

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

(11)

式中

(12)

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

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

(13)

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

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

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

IPP_ADCIG(x,βk)

(14)

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

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

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

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

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

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

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

IPS_ADCIG(x,βk)

(15)

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

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

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

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

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

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

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

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

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

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

(16)

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

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

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

2 模型驗證

2.1 水平層狀模型

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

圖3 水平層狀模型

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

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

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

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

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

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

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

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

2.2 Salt模型

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

圖8 Salt縱波速度模型

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

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

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

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

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

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

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

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

3 結論和認識

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

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

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

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