基于Chebyshev多項(xiàng)式的非對(duì)稱走時(shí)Kirchhoff疊前時(shí)間偏移角道集求取

劉志遠(yuǎn)，常 旭，劉伊克，王一博，杜向東

1 中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院工程地質(zhì)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100029

2 中海油研究總院，北京 100027

3 中石化勘探開發(fā)研究院油氣地球物理研究中心，北京 100083

1 引 言

偏移是地震成像技術(shù)中最重要的環(huán)節(jié)之一，Kirchhoff疊前時(shí)間偏移對(duì)觀測系統(tǒng)適應(yīng)性強(qiáng)、計(jì)算成本低，在生產(chǎn)中應(yīng)用廣泛.Kirchhoff疊前時(shí)間偏移的關(guān)鍵技術(shù)之一是計(jì)算地震波走時(shí)，為了使其獲得較高的精度，前人做過很多工作.目前疊前時(shí)間偏移走時(shí)方程中常見的有：二階雙曲NMO方程、時(shí)移雙曲 NMO 方程［1］、六階近似 NMO 方程［2］、Alkhalifah連分式逼近 NMO 方 程［3－4］、射 線追蹤時(shí)差方程［5］、分式展開二次方程［6］以及基于李群算法和擬微分算子象征理論的非對(duì)稱走時(shí)［7］方程等.切比雪夫多項(xiàng)式在逼近理論中有重要的應(yīng)用［8－10］，本文提出基于第一類切比雪夫多項(xiàng)式的非對(duì)稱走時(shí)方程，該方程能在對(duì)稱彎曲射線走時(shí)方程的基礎(chǔ)上添加切比雪夫非對(duì)稱項(xiàng)，達(dá)到計(jì)算量少，計(jì)算精度高的效果.

在油氣勘探開發(fā)的地震方法中，角度域共成像道集（ADCIG）是連接地震數(shù)據(jù)處理與儲(chǔ)層反演的重要紐帶.為抽取更準(zhǔn)確的時(shí)間域ADCIG，對(duì)走時(shí)計(jì)算加以改進(jìn)一直是研究的熱點(diǎn).自1990年de Bruin提出計(jì)算地震數(shù)據(jù)反射系數(shù)隨角度變化開始［11］，國內(nèi)外很多學(xué)者對(duì)ADCIG開展了一系列的研究工作，成為近年來的研究熱點(diǎn)之一.角道集主要有以下幾點(diǎn)應(yīng)用：其一，基于共成像道集的偏移速度分析［12－13］，在 ADCIG 上，當(dāng)偏移速度與實(shí)際速度一致時(shí)，道集的同相軸曲線不隨入射角的變化而變化，是拉平的，當(dāng)偏移速度與實(shí)際速度不一致時(shí)，道集的同相軸曲線就隨入射角增大上翹或者下彎；其二，在ADCIG上可以進(jìn)行AVA分析［14］；其三，可對(duì)角道集進(jìn)行適當(dāng)?shù)娜ピ搿⑹Ｓ嗲市Ｕ忍幚?，由于角度域共成像點(diǎn)道集能夠提供其它類型道集上不易表示的地震波信息，因此實(shí)施有針對(duì)性的數(shù)據(jù)處理后再疊加成像，可提高成像的質(zhì)量［15］.角道集研究用于疊前深度偏移（PSDM）的較多［16－17］，但 PSDM 速度建模困難，計(jì)算效率低，而疊前時(shí)間偏移（PSTM）角道集的獲取效率較高［18］，因此基于PSTM方法的角道集研究得以展開.近年來，針對(duì)PSTM的ADCIG研究方法有，基于直射線近似的Kirchhoff積分角度域成像方法［19］，基于射線追蹤的表驅(qū)Kirchhoff疊前時(shí)間偏移角度域算法以及非雙曲時(shí)差方位／角度域疊前時(shí)間偏移方法［5］，基于平面波分解的射線參數(shù)域成像道集方法［20］，基于李代數(shù)的非對(duì)稱走時(shí)角道集獲取方法［21］等.本文角道集的獲得采用基于切比雪夫多項(xiàng)式的非對(duì)稱彎曲射線走時(shí).

2 方法原理

2.1 切比雪夫多項(xiàng)式非對(duì)稱走時(shí)計(jì)算方法

2.1.1 第一類切比雪夫多項(xiàng)式基本概念

切比雪夫多項(xiàng)式與棣美弗定理有關(guān)，是以遞歸方式定義的一系列正交多項(xiàng)式序列.第n階第一類切比雪夫多項(xiàng)式以符號(hào)Tn表示.第一類切比雪夫多項(xiàng)式有自變量w介于－1和1之間的均勻分布的n個(gè)零點(diǎn)和n＋1個(gè)極值點(diǎn)，其最大誤差在各種級(jí)數(shù)中是最小的，它的根可以用于多項(xiàng)式插值，相應(yīng)的插值多項(xiàng)式能最大限度地降低龍格現(xiàn)象［8－10］，并且提供多項(xiàng)式連續(xù)函數(shù)的最佳一致逼近.其遞推公式為

文中用到的第一類切比雪夫多項(xiàng)式的前幾項(xiàng)可以寫為

相對(duì)應(yīng)的，

2.1.2 對(duì)稱彎曲射線泰勒走時(shí)方程分析

Slotnick于1959年［22］根據(jù)Snell定律，給出了射線參數(shù)表征的地震反射波時(shí)距曲線方程（4）和（5）：

t（p）為射線走時(shí)，x（p）為偏移距，p為射線參數(shù)，vk為第k層地震波波速，Δtk為第k層垂直射線走時(shí).

求解此方程的方法以對(duì)稱走時(shí)理論方法居多.Taner和Koehler于1969年將方程（4）、（5）進(jìn)行Taylor展開［23］，得到了對(duì)稱彎曲射線走時(shí)方程（6）.通常情況下，為求得高精度的走時(shí)計(jì)算式，有兩類方法：一種是對(duì)泰勒展開項(xiàng)的系數(shù)進(jìn)行優(yōu)化，一種是采用高階走時(shí)計(jì)算式［4］.具體而言，優(yōu)化六階的走時(shí)精度是目前精度比較高的.不過，兩種方法各有不足：采用高階走時(shí)計(jì)算式會(huì)增加計(jì)算量，尤其在計(jì)算大規(guī)模地震數(shù)據(jù)時(shí)所需計(jì)算時(shí)間較長；而采用優(yōu)化泰勒系數(shù)的方法時(shí)，需要先對(duì)走時(shí)計(jì)算式進(jìn)行優(yōu)化，目前優(yōu)化方法主要是六階優(yōu)化方程方法［2］和四、六階調(diào)節(jié)系數(shù)走時(shí)方程方法［4］，因?yàn)閮烧叩碾A數(shù)比較高，所以并沒有減少計(jì)算量.以上情況限制了高精度計(jì)算走時(shí)在實(shí)際生產(chǎn)中的應(yīng)用.而本文提出切比雪夫三階優(yōu)化走時(shí)，既能節(jié)省計(jì)算量，又能保證較高的精度.

其中各個(gè)系數(shù)為

dk為第k層厚度.

2.1.3 切比雪夫多項(xiàng)式非對(duì)稱走時(shí)方程

在泰勒展開式中，偏移距x屬于（k1，k2）區(qū)間，設(shè)w＝ （x－q）／p，q＝ （k1＋k2）／2，p＝ （k2－k1）／2，即可將f（x）等價(jià)變化到f（w），w 屬于（－1，1），于是，C′i＝f（ai，p，q），將四階泰勒展開走時(shí)方程等價(jià)變化為

這時(shí)，用第一類切比雪夫多項(xiàng)式對(duì)f（w）進(jìn)行轉(zhuǎn)化，得到：

式中Mi＝f（Cj），截去高階項(xiàng)T4，再反變換回f（x），最終得到：

x1對(duì)公式的誤差影響較大，且有礙于優(yōu)化，去掉此項(xiàng)，并賦予x3一調(diào)節(jié)系數(shù)m，用于調(diào)節(jié)截取的誤差，賦予其初始值m＝1，得到：

本文以模擬退火法對(duì)m進(jìn)行優(yōu)化，用此式逼近彎曲射線優(yōu)化六階泰勒走時(shí)的精度，即可將m優(yōu)化得到值m′，最終的切比雪夫多項(xiàng)式非對(duì)稱走時(shí)方程如下：

需要注意的是，在優(yōu)化m的時(shí)候，會(huì)增加一定的復(fù)雜性：一般在偏移孔徑內(nèi)，優(yōu)化一個(gè)m′即可，但如果數(shù)據(jù)要求進(jìn)行長偏移距的走時(shí)計(jì)算，則需要進(jìn)行分段優(yōu)化［10］，使得不同偏移距范圍具有不同的m′值；另一方面，與泰勒展開走時(shí)計(jì)算等常規(guī)方法相比，優(yōu)化m本身雖然多了一步工作量，不過優(yōu)化步驟并不會(huì)消耗大量的計(jì)算時(shí)間.

設(shè)計(jì)一個(gè)三層水平層狀模型進(jìn)行走時(shí)分析與驗(yàn)證，如圖1所示，網(wǎng)格大小為400×800，網(wǎng)格間距比例為10m，即模型X方向8000m，Z方向4000m.在對(duì)模型正演時(shí)，采樣間隔4ms，采樣時(shí)長6s，觀測系統(tǒng)采用全局接收，道間距20m，炮點(diǎn)位置選擇在模型網(wǎng)格點(diǎn)坐標(biāo)（1，1）處.用單炮數(shù)據(jù)對(duì)m進(jìn)行優(yōu)化.

在圖2的走時(shí)精度分析與圖3走時(shí)誤差對(duì)比中，橫坐標(biāo)X／Dmax為偏移距X與最大深度Dmax的比值.在走時(shí)計(jì)算式中，四階泰勒走時(shí)具有計(jì)算量低的優(yōu)點(diǎn)，六階優(yōu)化泰勒走時(shí)具有精度高的優(yōu)點(diǎn)，本文將優(yōu)化三階切比雪夫走時(shí)與上述兩者進(jìn)行比較，如圖2和圖3，證明優(yōu)化后的三階切比雪夫走時(shí)比高階的四階泰勒走時(shí)精度高，接近更高精度的優(yōu)化六階泰勒走時(shí).并且由走時(shí)計(jì)算公式可以很明顯的發(fā)現(xiàn)，優(yōu)化三階切比雪夫走時(shí)的計(jì)算量在三者中是最小的.

圖1 水平三層速度模型Fig.1 3layers velocity model

2.2 角道集的獲取

可以輸入任一未經(jīng)偏移的地震道集，見圖4.對(duì)于輸入的每一道地震數(shù)據(jù)，掃描地下成像點(diǎn)，入射線走時(shí)為ts，反射線走時(shí)為tr.角度計(jì)算采用公式（12）、（13），此種計(jì)算方法考慮射線彎曲效應(yīng)，計(jì)算地下散射點(diǎn)處真實(shí)的地震射線角度［21］.

如圖4，當(dāng)滿足ts＋tr＝ti時(shí)，將地震道上對(duì)應(yīng)的ti時(shí)刻的能量e加權(quán)歸位到成像點(diǎn)對(duì)應(yīng)的i z時(shí)刻上，加權(quán)采用保幅權(quán)因子［24］公式（14）：

把所有地震道依次這樣處理后將同一成像點(diǎn)相同角度能量疊加，不同角度橫向排列，便可得到i x成像點(diǎn)的角度域共成像點(diǎn)道集.依次類推，便可以求得各個(gè)成像點(diǎn)的角道集.這樣求取角道集的方法能克服水平層狀介質(zhì)限制，考慮地層傾斜影響，可以求得真實(shí)入射角度，使能量正確歸位，并且可以由原始的地震數(shù)據(jù)中直接求得角道集.

3 模型和實(shí)際數(shù)據(jù)測試

3.1 模型數(shù)據(jù)應(yīng)用

采用數(shù)據(jù)模型，見圖5，橫向1000點(diǎn)縱向400點(diǎn)，樣點(diǎn)間隔均為10m，采樣間隔4ms，采樣時(shí)長6s，觀測系統(tǒng)采用全局接收，炮間距40m，道間距20m，共250炮.角度分辨率為1°.我們計(jì)算如圖水平方向280號(hào)CMP坡度點(diǎn)處的ADCIG.

圖5 四層速度模型Fig.5 4layers model

由圖6所示可以發(fā)現(xiàn)，經(jīng)切比雪夫改進(jìn)的三階非對(duì)稱走時(shí)計(jì)算方法計(jì)算出的角道集，其精度與圖6c優(yōu)化六階泰勒走時(shí)展開的精度十分接近；與圖6a四階泰勒走時(shí)計(jì)算方法相比，其第二層大角度處的同相軸有所延長，且更加連續(xù)、平整，能量也得以較為準(zhǔn)確的歸位.在表1中列出了只改變走時(shí)計(jì)算方法的情況下三種計(jì)算方法所需時(shí)間.表1給出的數(shù)據(jù)說明，三階切比雪夫非對(duì)稱走時(shí)方法在時(shí)間上比泰勒四階走時(shí)快8.4%，比優(yōu)化六階泰勒走時(shí)快24.2%.即證明，經(jīng)切比雪夫多項(xiàng)式改進(jìn)求得的非對(duì)稱走時(shí)同時(shí)具備計(jì)算量較少和精度較高的優(yōu)勢(shì).

3.2 實(shí)際數(shù)據(jù)應(yīng)用

取某一實(shí)際數(shù)據(jù)為例，offset范圍（100m，4875m），道間隔25m，角度分辨率為1°，見圖7，文中顯示了從1300ms到3000ms的三種不同計(jì)算走時(shí)得出的同一成像點(diǎn)的角道集.

如圖7b中三階切比雪夫非對(duì)稱走時(shí)的計(jì)算方法精度較高，其ADCIG同相軸與圖7c優(yōu)化六階泰勒對(duì)稱走時(shí)展開的結(jié)果精度幾乎一樣，比圖7a四階泰勒走時(shí)的結(jié)果同相軸更加平整.在表1中列出了在計(jì)算該成像點(diǎn)時(shí)，三種計(jì)算方法計(jì)算時(shí)間的對(duì)比.表1給出的數(shù)據(jù)說明，求取此成像點(diǎn)的ADCIG采用經(jīng)切比雪夫改進(jìn)的三階非對(duì)稱走時(shí)方法比四階泰勒走時(shí)的計(jì)算結(jié)果相對(duì)縮短了約3.6%.圖8列舉了三種走時(shí)方法計(jì)算出的角道集經(jīng)過疊加后的最終剖面的一部分，可以看出，圈出部分的同相軸的連續(xù)性和清晰度由不好到較好到很好的變化.由此得到與模型數(shù)據(jù)相同的結(jié)果，即，經(jīng)切比雪夫多項(xiàng)式改進(jìn)求得的非對(duì)稱走時(shí)同時(shí)具備計(jì)算量較少和精度較高的優(yōu)勢(shì).

圖8 四階泰勒展開走時(shí)的剖面（a）；改進(jìn)的切比雪夫三階非對(duì)稱走時(shí)的剖面（b）；優(yōu)化六階泰勒展開走時(shí)的剖面（c）Fig.8 Shows the 4－order taylor method profile result（a），Shows the optimized 3－order asymmetric Chebyshev method profile result（b）；shows the optimized 6－order taylor profile result（c）

表1 走時(shí)計(jì)算量統(tǒng)計(jì)與對(duì)比Table 1 Contrast of the 3travel－time calculation methods

4 結(jié)論與討論

本文引進(jìn)第一類切比雪夫多項(xiàng)式，對(duì)泰勒展開對(duì)稱走時(shí)進(jìn)行改進(jìn)，推導(dǎo)出非對(duì)稱走時(shí)奇數(shù)項(xiàng)，得到低階切比雪夫非對(duì)稱走時(shí)計(jì)算公式，通過在數(shù)值模型和實(shí)際數(shù)據(jù)中抽取角道集的實(shí)驗(yàn)，取得了較好的結(jié)果，驗(yàn)證了該方法具有高階走時(shí)的精度，同時(shí)比實(shí)際生產(chǎn)中較為廣泛應(yīng)用的低階高精度走時(shí)消耗更少的計(jì)算量，而且應(yīng)用這種走時(shí)直接求得的角道集，能較真實(shí)的反映地震射線在地下層位的入射角度.

切比雪夫非對(duì)稱走時(shí)與基于李代數(shù)理論的非對(duì)稱走時(shí)［7，25］具有不同的理論計(jì)算公式，切比雪夫非對(duì)稱走時(shí)的計(jì)算量較小，在實(shí)際應(yīng)用中有一定的優(yōu)勢(shì).

在切比雪夫非對(duì)稱走時(shí)計(jì)算中，優(yōu)化算子m′的取值較為關(guān)鍵.針對(duì)不同的地震數(shù)據(jù)，k1、k2值各有不同，優(yōu)化方法也可以根據(jù)情況擇優(yōu)選取，因此優(yōu)化得到的m′值也不盡一致，但一般情況下都可以達(dá)到計(jì)算量較小和精度較高的計(jì)算效果.

（References）

［1］Castle R J.A theory of normal moveout.Geophysics，1994，59（6）：983－999.

［2］Sun C W，Wang H W，Nartinez R D.Optimized 6th order NMO correction for long－offset seismic data.72nd Annual International Meeting，SEG，Expanded Abstracts，2002.

［3］Alkhalifah T，Tsvankin I.Velocity analysis for transversely isotropic media.Geophysics，1995，60（5）：1550－1566.

［4］尤建軍，常旭，劉伊克.VTI介質(zhì)長偏移距非雙曲動(dòng)校正公式優(yōu)化.地球物理學(xué)報(bào)，2006，49（6）：1770－1778.You J J，Chang X，Liu Y K.Optimization of nonhyperbolic moveout correction equation of long－offset seismic data in VTI media.Chinese J.Geophys.（in Chinese），2006，49（6）：1770－1778.

［5］程玖兵，王楠，馬在田.表驅(qū)三維角度域Kirchhoff疊前時(shí)間偏移成像方法.地球物理學(xué)報(bào)，2009，52（3）：792－800.Cheng J B，Wang N，Ma Z T.Table－driven 3－D angledomain imaging approach for Kirchhoff prestack time migration.Chinese J.Geophys.（in Chinese），2009，52（3）：792－800.

［6］劉洋.反射波分式展開時(shí)距方程及其精度分析.石油物探，2003，42（4）：441－447.Liu Y.Fraction expansion time－distance equation of reflection wave and its accuracy analysis.Geophysical Progress for Petroleum （in Chinese），2003，42（4）：441－447.

［7］劉洪，劉國鋒，李博等.基于橫向?qū)?shù)的走時(shí)計(jì)算方法及其在疊前時(shí)間偏移中的應(yīng)用.石油物探，2009，48（1）：3－10.Liu H，Liu G F，Li B，et al.The travel time calculation method via lateral derivative of velocity and its application in pre－stack time migration.Geophysical Progress for Petroleum（in Chinese），2009，48（1）：3－10.

［8］Zhang J H，Wang W M，Wang S Q.Optimized Chebyshev Fourier migration：A wide－angle dual－domain method for media with strong velocity contrasts.Geophysics，2010，75（2）：S23－S24.

［9］丁帆，張金海，姚振興.長偏移距地震資料的優(yōu)化契比雪夫動(dòng)校正方法.地球物理學(xué)進(jìn)展，2011，26（3）：836－842.Ding F，Zhang J H，Yao Z X.Optimized Chebyshev method for normal moveout of long－offset seismic data.Progress in Geophysics （in Chinese），2011，26（3）：836－842.

［10］劉璐，梁光河，符超等.基于Chebyshev多項(xiàng)式的彎曲射線Kirchhoff疊前時(shí)間偏移.地球物理學(xué)報(bào)，2011，54（10）：2665－2672.Liu L，Liang G H，F(xiàn)u C，et al.Bend－ray Kirchhoff pre－stack time migration based on Chebyshev polynomial.Chinese J.Geophys.（in Chinese），2011，54（10）：2665－2672.

［11］de Bruin，C G M，Wapenaar C P A，Berkhout A J.Angledependent reflectivity by means of prestack migration.Geophysics，1990，55（9）：1223－1234.

［12］Biondo B， William W S. Angle－domain common－image gathers for migration velocity analysis by wavefieldcontinuation imaging.Geophysics，2004，69（5）：1283－1298.

［13］劉奇琳，劉伊克，常旭.雙平方根波動(dòng)方程偏移速度分析.地球物理學(xué)報(bào)，2009，52（7）：1891－1898.Liu Q L，Liu Y K，Chang X.Wave－equation migration velocity analysis by double square root method.Chinese J.Geophys.（in Chinese），2009，52（7）：1891－1898.

［14］Brandsberg－Dahl S，de Hoop M V，Ursin B.Focusing in dip and AVA compensation on scattering－angle／azimuth common image gathers.Geophysics，2003，68（1）：232－254.

［15］Lee S，King D，Lin S.Efficient true－amplitude weights in Kirchhoff time migration.74th Annual International Meeting，SEG，Expanded Abstracts，2004：1089－1092.

［16］Rickett J，Sava P.Offset and angle－domain common imagepoint gathers for shot－profile migration.Geophysics，2002，67（3）：883－889.

［17］陳凌，吳如山，王偉君.基于Gabor－Daubechies小波束疊前深度偏移的角度域共成像道集.地球物理學(xué)報(bào)，2004，47（5）：876－885.Chen L，Wu R S，Wang W J.Common angle image gathers obtained from Gabor－daubechies beamlet prestack depth migration.Chinese J.Geophys.（in Chinese），2004，47（5）：876－885.

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

［19］Perez G， Marfurt K J.Improving lateral and vertical resolution of seismic images by correcting for wavelet stretch in common－angle migration.Geophysics，2007，72（6）：94－104.

［20］王棣，王華忠，馬在田等.疊前時(shí)間偏移方法綜述.勘探地球物理進(jìn)展，2004，27（5）：313－316.Wang D，Wang H Z，Ma Z T，et al.Review of prestack time migration methods.Progress in Exploration Geophysics （in Chinese），2004，27（5）：313－316.

［21］鄒振，劉洪，劉紅偉.Kirchhoff疊前時(shí)間偏移角度道集.地球物理學(xué)報(bào)，2010，53（5）：1207－1214.Zou Z，Liu H，Liu H W.Common－angle gathers based on Kirchhoff pre－stack time migration.Chinese J.Geophys.（in Chinese）.2010，53（5）：1207－1214.

［22］Slotnick M.Lessons in Seismic Computing：A Memorial to the Author.Lawrence：Society of Exploration Geophysicists，1959.

［23］Taner M T，Koehler F.Velocity spectra－digital computer derivation applications of velocity functions.Geophysics，1969，34（6）：859－881.

［24］Dellinger J A，Gray S H，Murphy G E，et al.Efficient 2.5－Dtrue－amplitude migration.Geophysics，2000，65（3）：943－950.

［25］張廉萍，劉洪.適于Kirchhoff疊前深度偏移的地震走時(shí)李代數(shù)積分算法.地球物理學(xué)報(bào)，2010，53（8）：1893－1901.Zhang L P，Liu H.Lie algebra integral algorithm of traveltime calculation for pre－stack Kirchhoff depth migration.Chinese J.Geophys.（in Chinese），2010，53（8）：1893－1901.