逆時(shí)偏移中用Poynting矢量高效地提取角道集

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

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

2 中煤科工集團(tuán)西安研究院，西安 710077

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

1 引 言

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

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

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

2 方法原理

2.1 波場(chǎng)傳播

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

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

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

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

2.2 Poynting矢量求取反射角

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

Yoon和 Marfurt（2003）［5］給出了地震波場(chǎng)中的坡印廷矢量計(jì)算式：

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

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

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

2.3 提取ADCIGs

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

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

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

具體的計(jì)算流程如圖2.

3 模型試驗(yàn)

3.1 Marmousi模型的ADCIGs

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

3.2 低頻噪聲

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

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

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

3.3 AVA

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

4 實(shí)際數(shù)據(jù)試驗(yàn)

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

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

5 結(jié) 論

本文在逆時(shí)偏移過(guò)程中采用一階應(yīng)力速度波動(dòng)方程進(jìn)行波場(chǎng)傳播，并用Poynting矢量來(lái)確定傳播方向并計(jì)算反射角，進(jìn)而通過(guò)給定的成像條件進(jìn)行逆時(shí)偏移角道集的提取.整個(gè)計(jì)算過(guò)程中幾乎不額外增加計(jì)算量和存儲(chǔ)量，提取的共角度道集可用于進(jìn)行速度分析和AVA分析等.該方法可用于三維逆時(shí)偏移中角道集等的提取（包括方位角道集），實(shí)現(xiàn)簡(jiǎn)單，易于進(jìn)行并行加速（如GPU等）.

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

（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］張美根，王妙月.各向異性彈性波有限元疊前逆時(shí)偏移.地球物理學(xué)報(bào)，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］李博，劉國(guó)峰，劉洪.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法.地球物理學(xué)報(bào)，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］李博，劉紅偉，劉國(guó)峰等.地震疊前逆時(shí)偏移算法的CPU／GPU實(shí)施對(duì)策.地球物理學(xué)報(bào)，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］劉紅偉，李博，劉洪等.地震疊前逆時(shí)偏移高階有限差分算法及 GPU實(shí)現(xiàn).地球物理學(xué)報(bào)，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］劉紅偉，劉洪，鄒振等.地震疊前逆時(shí)偏移中的去噪與存儲(chǔ).地球物理學(xué)報(bào)，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］劉紅偉，劉洪，李博等.起伏地表疊前逆時(shí)偏移理論及GPU加速技術(shù).地球物理學(xué)報(bào)，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］楊午陽(yáng).模擬退火疊前AVA同步反演方法.地球物理學(xué)進(jìn)展，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é)進(jìn)展，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］喻岳鈺，楊長(zhǎng)春，王彥飛等.疊前彈性阻抗反演及其在含氣儲(chǔ)層預(yù)測(cè)中的應(yīng)用.地球物理學(xué)進(jìn)展，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.