王保利，高靜懷，陳文超，張喚蘭

1 西安交通大學(xué)電子與信息工程學(xué)院，西安 710049

2 中煤科工集團西安研究院，西安 710077

3 西安科技大學(xué)地質(zhì)與環(huán)境學(xué)院，西安 710054

1 引 言

20世紀80年代就提出的逆時偏移［1－2］在對復(fù)雜地下結(jié)構(gòu)進行成像時體現(xiàn)出了較強的優(yōu)勢，因此越來越受到人們的重視［3－10］.多路徑、多次波、大傾角、回轉(zhuǎn)波、棱鏡波以及保幅等關(guān)鍵問題上，逆時偏移幾乎都可以較好地解決.但準確的逆時偏移成像比其它成像方法（如單程波偏移、Kirchhoff偏移和beam偏移等）對速度模型的準確度更加敏感.為了給逆時偏移建立一個較準確的速度模型，需要一種迭代的方法來不斷地更新速度.而共成像點道集（CIGs）是疊前深度偏移中的很重要的輸出道集，它可以提供速度建模所需要的信息，以及反映地下巖性的振幅和相位等信息.正確的速度模型得到的共成像點道集中的同相軸是平的，而不正確的速度得到的是彎曲的，且彎曲的程度與速度校正量是有關(guān)系的，通過這種關(guān)系就可以不斷地更新速度模型［11］.

通常，有兩種共成像點道集，一種是偏移距域的（common－offset CIGs，簡稱為COCIGs），一種是角度域的（angle－domain CIGs，簡稱為 ADCIGs）.在地下反射體是水平層狀的情況，這兩種道集是可互換的.Kirchhoff偏移和beam偏移方法可以很方便快速地輸出COCIGs［12］.但多路徑問題會影響到COCIGs的質(zhì)量，進而會影響速度分析和AVO分析的精度.而單程波偏移和逆時偏移中的COCIGs是很難提取的.Save和 Formel（2003）［13］利用單程波偏移在二維情況下首先提取了地下局部偏移距域的CIGs，然后再轉(zhuǎn)換為 ADCIGs.Fomel（2004）［14］將這種方法推廣到三維情況，并指出三維情況下該方法非常復(fù)雜.Save和 Fomel（2006）［15］又提出一種新的方法，該方法需要先生成延時CIGs，然后再轉(zhuǎn)換為 ADCIGs.Xu等 （2010，2011）［16－17］提出一種直接從三維逆時偏移中提取ADCIGs的方法，該方法在頻率波數(shù)域求取傳播方向，進而計算出反射角和方位角信息，最后使用給定的成像條件直接輸出ADCIGs.這種方法可以得到質(zhì)量較好的ADCIGs，但需要保存所有時刻的炮點波場值和接收點波場值，并進行四維的傅里葉變換，計算量和存儲量都非常得大，不適合進行逆時偏移速度分析等處理.

Yoon和 Marfurt（2003）［5］用坡印廷矢量來確定傳播方向并用此修改成像條件來去除逆時偏移中的低頻噪聲，取得了較好的效果.本文借助于這種思想，用一階應(yīng)力速度方程（與常規(guī)標量應(yīng)力波動方程相比，增加了額外的存儲量，但減少了計算坡印廷矢量的耗時）進行波場延拓，并用坡印廷矢量求取成像點處的傳播方向，最后使用給定的成像條件輸出逆時偏移中的ADCIGs.

2 方法原理

2.1 波場傳播

逆時偏移是通過波動方程在時間上對采集到的地震數(shù)據(jù)進行反向傳播來實現(xiàn)的.其偏移過程使用的是雙程波波動方程，因此不受傾角限制，且能正確地對回轉(zhuǎn)波、多次波等進行成像.逆時偏移的原理本文不再贅述.

通常，逆時偏移采用的聲波波動方程為標量應(yīng)力波動方程［18－22，10］，其形式如下：

其中，υ為介質(zhì)縱波速度，P為質(zhì)點應(yīng)力波場.該波動方程不利于使用坡印廷矢量（原因見下節(jié)）.因此本文采用一階應(yīng)力－速度聲波波動方程：

式中：ρ為介質(zhì)密度，υ為介質(zhì)縱波速度，v和w 分別為x和z軸上的質(zhì)點振動速度.本文對公式（2）采用交錯網(wǎng)格技術(shù)［18］來構(gòu)建波場傳播算子.考慮到逆時偏移的計算量和為避免數(shù)值頻散問題對逆時偏移結(jié)果的影響，在時間方向采用二階而在空間方向采用高階差分格式（≥8階）.

2.2 Poynting矢量求取反射角

坡印廷矢量也稱為能流密度，1884年由坡印廷提出用來測量單位時間內(nèi)穿過垂直于此矢量方向的單位表面的電磁場能量.Poynting矢量值等于電場強度與磁場強度的叉積，矢量方向代表電磁波的傳播方向.

Yoon和 Marfurt（2003）［5］給出了地震波場中的坡印廷矢量計算式：

圖1 坡印廷矢量快照圖（a）水平分量；（b）垂直分量.Fig.1 The snapshots of Poynting vector（a）Horizontal component；（b）Vertical component.

其中，v表示位移速度矢量，v＝（v，w），P為應(yīng)力波場.通過（3）式可知，如果使用（1）式的標量聲波波動方程，需計算應(yīng)力波場值的時間和空間導(dǎo)數(shù)，計算比較復(fù)雜.而如果采用（2）式的一階應(yīng)力速度方程，則應(yīng)力P和位移速度矢量v都是逆時偏移過程中已有的，不需任何額外計算.圖1是均勻介質(zhì)中的某時刻計算的坡印廷矢量快照圖，可以很明顯地判斷出波的大致傳播方向：在水平分量圖（a）中藍色和紅色分別代表向右傳和向左傳，而在垂直分量圖（b）中藍色和紅色則分別代表向下傳和向上傳.

用SS和SR分別表示炮點波場的坡印廷矢量和接收點波場的坡印廷矢量，則根據(jù)余弦定理，可以求出兩矢量的夾角θ，

2.3 提取ADCIGs

獲得反射角β后，便可以進行ADCIGs的提取.常規(guī)的逆時偏移成像條件是對所有角度道集的疊加，因此不能直接用來提取ADCIGs，需要進行修改.考慮到計算得到的反射角通常與我們給定的采樣角度（一般是等間隔的）不一致，我們用高斯采樣函數(shù)將常規(guī)的成像條件

其中S和R分別表示炮點波場和接收點波場（本文中為應(yīng)力波場值），而β表示當前時刻該成像點處計算得到的反射角，βk表示ADCIGs中的離散角，σ表示高斯函數(shù)的方差.根據(jù)（7）便能容易地提取出ADCGIs.

需特別說明的是，（5）式中，分母在很小甚至為零情況下，計算得到的角度β誤差較大，但此時，其能量貢獻也非常小，因此不會影響最終角道集的質(zhì)量.

具體的計算流程如圖2.

3 模型試驗

3.1 Marmousi模型的ADCIGs

用本文提出的方法提取出的Marmousi模型的ADCIGs如圖3所示.角度從0°～90°，角度間隔為2°，每100個CMP顯示一個ADCIG，一共23個ADCIGs.從圖上可以看出每個ADCIG上的同相軸都是拉平的，且低頻噪聲主要分布在大角度的道上.圖4為第500個CMP位置處的ADCIG.

3.2 低頻噪聲

逆時偏移的低頻噪聲主要是由透射波成像造成的，而入射波和透射波的夾角在模型比較光滑情況下接近于180°（傳播路徑一微小段范圍內(nèi)），因此可通過在角道集上對大角度的道進行衰減或直接置零來對逆時偏移成像剖面進行去噪.理論上這種去噪方法等同于拉普拉斯濾波［7，11］.圖5、圖6和圖7分別為用0°～40°、40°～70°和70°～90°等角道集疊加后的結(jié)果，可明顯看出在0°～40°內(nèi)幾乎沒有任何低頻噪聲，而在70°～90°的偏移剖面幾乎全部是低頻噪聲.

圖2 Poynting矢量提取角道集的流程圖Fig.2 The flow chart of extracting ADCIGs using Poynting vector

圖4 Marmousi模型第500個CMP位置處的ADCIGFig.4 A ADCIG located at 500th CMP of Marmousi model

3.3 AVA

角道集的另一個用途就是用于AVA 分析［23－25］，本文使用單層模型分析了角道集上振幅的特征.圖8為用振幅隨角度遞增的模型參數(shù)（上下層速度分別為1500m／s和1800m／s，密度均為2.1g／cm3）得到的角道集AVA曲線與理論AVA曲線的對比.圖9為用振幅隨角度遞減的模型參數(shù)（上下層速度分別為1500m／s和1300m／s，密度均為2.1g／cm3）得到的角道集AVA曲線與理論AVA曲線的對比.對比可以看出，兩種模型得到的角道集的AVA曲線和理論AVA曲線基本一致.

4 實際數(shù)據(jù)試驗

圖10為從某地區(qū)的一條測線中提取出的部分ADCIGs.可以看出1.3～1.8km段中同相軸基本是平的，而在1.8～2.3km斷層同相軸是向下彎曲的，說明偏移過程中所用的速度偏高.通過不斷迭代地進行速度調(diào)整即可實現(xiàn)逆時偏移速度分析.

圖10 實際數(shù)據(jù)中提取的部分角道集Fig.10 Part of the ADCIGs from field data

5 結(jié) 論

本文在逆時偏移過程中采用一階應(yīng)力速度波動方程進行波場傳播，并用Poynting矢量來確定傳播方向并計算反射角，進而通過給定的成像條件進行逆時偏移角道集的提取.整個計算過程中幾乎不額外增加計算量和存儲量，提取的共角度道集可用于進行速度分析和AVA分析等.該方法可用于三維逆時偏移中角道集等的提?。òǚ轿唤堑兰?，實現(xiàn)簡單，易于進行并行加速（如GPU等）.

后續(xù)工作將在此工作基礎(chǔ)上，使用提取的角道集進行逆時偏移速度分析，以得到較好的偏移速度剖面.

（References）

［1］Whitmore N D.Iterative depth migration by backward time propagation.53rdAnnual International Meeting，SEG Expanded Abstracts，1983，2：827－830.

［2］Baysal E，Kosloff D D，Sherwood J W C.Reverse time migration.Geophysics，1983，48（11）：1514－1524.

［3］張美根，王妙月.各向異性彈性波有限元疊前逆時偏移.地球物理學(xué)報，2001，44（5）：711－719.Zhang M G Wang M Y.Prestack finite element reverse－time migration for anisotropic elastic waves.Chinese J.Geophys.（in Chinese），2001，44（5）：711－719.

［4］Yoon K，Shin C，Hong S.3Dreverse－time migration using acoustic wave equation.An experience with the SEG／EAGE data set.The Leading Edge，2003，22（1）：38－41.

［5］Yoon K，Marfurt Kurt J，Starr W.Challenges in reversetime migration.74thAnnual International Meeting，SEG Expanded Abstracts，2003：1057－1060.

［6］Bednar J B，Bednar C J，Shin C.Two－way versus one－way：A case study style comparison.76thAnnual International Meeting，SEG Expanded Abstracts，2006，25：2343－2347.

［7］Zhang Yu，Xu S，Bleistein N，et al.True－amplitude，angledomain，common－image gathers from one－way wave－equation migrations.Geophysics，2007，72（1）：S49－S58.

［8］Zhang Y，James S. Practical issues of reverse time migration：True amplitude gathers，noise removal and harmonic－source encoding.First Break，2009，27（1）：53－60.

［9］Liu F Q，Zhang G Q，Morton S A，et al.Reverse time migration using one way wavefield imaging condition.77thAnnual International Meeting，SEG Expanded Abstracts，2007，26（1）：2170－2174

［10］Liu F Q，Zhang G Q，Morton S A，et al.An effective imaging condition for reverse－time migration using wavefield decomposition.Geophysics，2011，76（1）：S29－S39.

［11］Zhou H B，Gray S H，Young J，et al.Tomographic residual curvature analysis：The process and its components.73rdAnnual International Meeting，SEG Expanded Abstracts，2003，22：666－669.

［12］Xu S，Chauris H，Lambaré G，et al.Common－angle migration：A strategy for imaging complex media.Geophysics，2001，66（6）：1877－1894.

［13］Sava P C，F(xiàn)ormel S.Angle－domain common－image gathers by wavefield continuation methods.Geophysics，2003，68（3）：1065－1074.

［14］Fomel S.Theory of 3－D angle gathers in wave－equation imaging.74thAnnual International Meeting，SEG Expanded Abstracts，2004，1（1）：1053－1056.

［15］Sava P，F(xiàn)ormel S.Time－shift imaging condition in seismic migration.Geophysics，2006，71（6）：S209－S217.

［16］Xu S，Zhang Y，Tang B.3Dcommon image gathers from Reverse time migration.80thAnnual International Meeting，SEG Expanded Abstracts，2010：3257－3260.

［17］Xu S，Zhang Y，Tang B.3Dangle gathers from Reverse time migration.Geophysics，2011，76（2）：S77－S92.

［18］李博，劉國峰，劉洪.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法.地球物理學(xué)報，2009，52（1）：245－252.Li B，Liu G F，Liu H.A method of using GPU to accelerate seismic pre－stack time migration.Chinese J.Geophys.（in Chinese），2009，52（1）：245－252.

［19］李博，劉紅偉，劉國峰等.地震疊前逆時偏移算法的CPU／GPU實施對策.地球物理學(xué)報，2010，53（12）：2938－2943.Li B，Liu H W，Liu G F，et al.Computational strategy of seismic pre－stack reverse time migration on CPU／GPU.Chinese J.Geophys.（in Chinese），2010，53（12）：2938－2943.

［20］劉紅偉，李博，劉洪等.地震疊前逆時偏移高階有限差分算法及 GPU實現(xiàn).地球物理學(xué)報，2010，53（7）：1725－1733.Liu H W，Li B，Liu H，et al.The algorithm of high order finite difference pre－stack reverse time migration and GPU implementation.Chinese J.Geophys.（in Chinese），2010，53（7）：1725－1733.

［21］劉紅偉，劉洪，鄒振等.地震疊前逆時偏移中的去噪與存儲.地球物理學(xué)報，2010，53（9）：2171－2180.Liu H W，Liu H，Zou Z，et al.The problems of denoise and storage in seismic reverse time migration.Chinese J.Geophys.（in Chinese），2010，53（7）：2171－2180.

［22］劉紅偉，劉洪，李博等.起伏地表疊前逆時偏移理論及GPU加速技術(shù).地球物理學(xué)報，2011，54（7）：1883－1892.Liu H W，Liu H，Li B，et al.Pre－stack reverse time migration for rugged topography and GPU acceleration technology.Chinese J.Geophys.（in Chinese），2011，54（7）：1883－1892.

［23］楊午陽.模擬退火疊前AVA同步反演方法.地球物理學(xué)進展，2010，25（1）：219－224.Yang W Y.Pre－stack simultaneous AVA inversion algorithm with simulated annealing.Progress in Geophys.（in Chinese），2010，25（1）：219－224.

［24］王保利.彈性阻抗反演及應(yīng)用研究.地球物理學(xué)進展，2005，20（1）：89－92.Wang B L.Elastic impedance inversion and its application.Progress in Geophys.（in Chinese），2005，20（1）：89－92.

［25］喻岳鈺，楊長春，王彥飛等.疊前彈性阻抗反演及其在含氣儲層預(yù)測中的應(yīng)用.地球物理學(xué)進展，2009，24（2）：574－580.Yu Y Y，Yang C C，Wang Y F，et al.Application of prestack seismic elastic impedance inversion to gas reservoir.Progress in Geophys.（in Chinese），2009，24（2）：574－580.