帶正則化校正的TTI介質(zhì)qP波方程及其逆時偏移方法和應(yīng)用

劉文卿，王西文，王宇超，雍學(xué)善，王小衛(wèi)，張濤

1 成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都　610059　2　中國石油勘探開發(fā)研究院西北分院，蘭州　730020

?

帶正則化校正的TTI介質(zhì)qP波方程及其逆時偏移方法和應(yīng)用

劉文卿1,2，王西文2，王宇超2，雍學(xué)善2，王小衛(wèi)2，張濤2

1 成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都6100592中國石油勘探開發(fā)研究院西北分院，蘭州730020

摘要各向異性研究對地下介質(zhì)精確成像有著重要的意義,在當(dāng)前計算機硬件迅速發(fā)展及寬方位地震數(shù)據(jù)采集日益普遍的情況下，成像必須考慮介質(zhì)的各向異性.逆時偏移是基于雙程波動方程的較為精確的數(shù)值解的成像方法，所以相對于其他地震成像方法，它具有很大的優(yōu)勢,譬如不受反射界面的傾角限制、偏移速度結(jié)構(gòu)合適時能夠使回轉(zhuǎn)波及多次波正確成像.在各向同性介質(zhì)中，可使用標量波方程來模擬波場.而在各向異性介質(zhì)中，P波和SV波是相互耦合的，即不存在單純的標量波傳播，通常利用能代表耦合波場中P波分量運動學(xué)特征的擬聲波(qP波)進行偏移成像.本文中，我們推導(dǎo)出了TTI介質(zhì)下qP波控制方程.該方程可采用顯式有限差分格式進行求解.通過聲學(xué)近似，若沿對稱軸方向的剪切波速度為零，對于對稱軸方向不變且ε≥δ的模型來說，可得到穩(wěn)定的數(shù)值解.但對于TTI介質(zhì)來說，由于沿對稱軸方向各向異性參數(shù)是變化的，聲學(xué)近似會引起波場傳播及數(shù)值計算的不穩(wěn)定.因此，我們提出了正則化有限橫波的方法，很好地解決了這一問題.最后，給出了Foothill模型的測試結(jié)果及某探區(qū)實際資料試算結(jié)果，展示了采用這個方程進行復(fù)雜TTI模型正演和高質(zhì)量逆時偏移成像結(jié)果，證實了該方法的正確性和實際資料應(yīng)用中的有效性.

關(guān)鍵詞TTI介質(zhì)； 各向異性； 波場耦合； 有限橫波； 正則化； 逆時偏移

1引言

隨著寬方位地震勘探技術(shù)的發(fā)展及勘探精度的不斷提高，各向異性地震成像受到業(yè)界更多的關(guān)注.逆時偏移方法基于精確的雙程波動理論，在時間方向上對震源波場進行正向延拓、檢波點波場反向延拓(Baysal et al.，1983；McMechan，1983；Whitmore，1983)，正反向延拓的波場零相位相關(guān)產(chǎn)生成像結(jié)果.該方法具有無反射傾角限制、能適應(yīng)強橫向變速,能夠?qū)剞D(zhuǎn)波及多次反射波進行正確成像的特點.在各向同性介質(zhì)中，可以對標量聲波方程進行求解，得到逆時偏移成像結(jié)果.嚴格來說實際接收的地震波場更接近于彈性波方程所描述的波場特征.因此，一些學(xué)者考慮能否直接用彈性波方程外推成像(Sun et al.，1986；Chang and McMechan，1986).但利用彈性波方程進行波場外推成像存在很多困難，首先是多參數(shù)建模困難，其次是直接外推耦合的方程組計算量巨大，另外多分量數(shù)據(jù)波場分解及相關(guān)成像難度很大.因此，目前開展彈性波逆時偏移成像實用化的條件還不成熟.但并不能否認彈性波逆時偏移成像的重要性.

對于各向異性介質(zhì)，嚴格意義上的標量波是不存在的.但對于一般的偏移成像，與各向同性介質(zhì)相似，可通過聲學(xué)近似的方式，能夠找到一個能代表TI介質(zhì)中P波運動學(xué)特征的控制方程.利用這樣的方程就能進行P波波場的外推和偏移成像.Alkhalifah(1998)首先提出了擬聲波近似思想，通過聲學(xué)近似簡化耦合的頻散關(guān)系，引入了TI介質(zhì)擬聲波近似頻散關(guān)系并導(dǎo)出方便數(shù)值求解的qP波控制方程.盡管這個標量波場的頻散關(guān)系具有和矢量波場中P波相似的運動學(xué)特征，但它會產(chǎn)生偽橫波.對于VTI介質(zhì)，基于Alkhalifah的頻散關(guān)系式，許多學(xué)者已經(jīng)推導(dǎo)出了擬聲波雙程波動方程并用于VTI介質(zhì)逆時深度偏移成像中，Alkhalifah(2000)，Zhou等(2006)，Du等(2008)，F(xiàn)letcher等(2008)，Duveneck等(2008)利用不同輔助變量簡化了擬聲波方程的求解.同時在二維情況下通過解耦的P-SV關(guān)系導(dǎo)出‘純P’波方程，但‘純P波’方程的快速數(shù)值算法仍不成熟(Liu et al.，2009；Chu et al.,2011).

從Hooke定律和運動方程出發(fā)，通過聲學(xué)近似也可推導(dǎo)出VTI介質(zhì)擬聲波方程.從平面波傳播的角度，VTI介質(zhì)向TTI介質(zhì)轉(zhuǎn)化可以看成在旋轉(zhuǎn)后的坐標系下的平面波傳播.坐標系旋轉(zhuǎn)后，僅增加了數(shù)學(xué)運算的復(fù)雜度，并不增加物理量綱.通過設(shè)置對稱軸方向SV波速度為零，可推導(dǎo)出TTI介質(zhì)下的擬聲波方程(Fletcher et al.，2009；Duveneck and Bakker,2011).

本文將著重介紹一個新的TTI介質(zhì)擬聲波方程推導(dǎo)及實現(xiàn)策略，該方程從精確P-SV波頻散關(guān)系推導(dǎo)而來，并且使用Thomsen(1986)定義的各向異性參數(shù)ε和δ進行參數(shù)化.對于對稱軸方向不變且ε≥δ的介質(zhì)(VTI或HTI介質(zhì))來說，通過聲學(xué)近似，可得到穩(wěn)定的數(shù)值解，可以直接將該方程應(yīng)用于正演和偏移中去.然而，對于非均勻TTI介質(zhì)中的正演和逆時偏移，各向異性對稱軸方向的變化會引起數(shù)值計算的不穩(wěn)定問題.為了得到更加穩(wěn)定的TTI介質(zhì)中qP波方程，F(xiàn)lecter等(2009)通過在聲學(xué)近似方程中加入非零橫波速度得到有限橫波方程來增加方程求解穩(wěn)定性.更自然的解決辦法由Duveneck和Bakker(2011)首先提出，即對彈性波方程直接做聲學(xué)近似，避開從頻散關(guān)系推導(dǎo)qP波方程時對稱軸方向緩變的假設(shè)，導(dǎo)出了TTI介質(zhì)中與彈性波方程對應(yīng)的qP波方程，該方程可以穩(wěn)定模擬P波分量的傳播.Zhang Y和Zhang H Z(2009)從構(gòu)造自共軛微分算子角度導(dǎo)出了一種TTI介質(zhì)中的qP波方程.該方程可以看作是Duveneck和Bakker(2011)導(dǎo)出的qP波方程的近似方程，其忽略了部分對稱軸方向空間導(dǎo)數(shù)項以提高計算效率，但保留了Duveneck方程中微分算子的自共軛性來保證穩(wěn)定性(Zhang Y and Zhang H Z,2009；Zhang H Z and Zhang Y，2011).本文采用一種新的方式增強qP波傳播時的穩(wěn)定性，采用的是有限橫波正則化的方法，這樣就可以消除由傾斜對稱軸變化以及由ε<δ引起的傳播不穩(wěn)定問題，同時最大程度地減弱波場傳播中的偽橫波影響.在數(shù)值試算及實際資料試算部分，展示了采用所導(dǎo)出方程的正演和逆時偏移結(jié)果，證明了方法理論的正確性.

2TTI介質(zhì)qP波控制方程推導(dǎo)

從平面波傳播的角度，從VTI介質(zhì)向TTI介質(zhì)轉(zhuǎn)化可看作是在旋轉(zhuǎn)坐標系下的平面波傳播.首先，基于Christoffel方程從相速度/頻散關(guān)系出發(fā)，可導(dǎo)出VTI介質(zhì)耦合頻散關(guān)系：

(1)

(2)

則(1)式可以重新寫為：

(3)

聯(lián)立(2)、(3)式,并將其反變換到時空域，得到VTI介質(zhì)中的有限橫波方程(Zhou et al.,2006)：

×(p-q).

(4)

同理,通過求解非均勻TTI介質(zhì)的Christoffel方程,可推導(dǎo)出旋轉(zhuǎn)坐標系下TTI介質(zhì)P-SV波耦合頻散關(guān)系:

(5)

由觀測坐標與旋轉(zhuǎn)坐標系下波數(shù)矢量的關(guān)系可得：

(6)

(7)

將(4)式中的二階微分算子換成(7)式旋轉(zhuǎn)坐標系下的形式，即可得TTI介質(zhì)中的有限橫波方程：

(8)

將(7)式對應(yīng)的空間微分算子用H1、H2重寫可得：

(9)

則有限橫波方程(8)可表示為：

(10)

式中，θ是傾角，φ是對稱軸的方位角.注意對于VTI介質(zhì)，微分算子H1只作用于垂直方向，而H2只作用于水平方向.兩個算子都不包含空間交叉導(dǎo)數(shù)項.在各向同性介質(zhì)中，式(10)中的兩個二階偏微分方程是等價的，波場p等價于波場q.

方程(10)可以采用偽譜法求解，求解中要用到波數(shù)域空間導(dǎo)數(shù)計算(Fletcher,2009).然而，由于偽譜法需要多次傅里葉變換，不利于大型數(shù)據(jù)的波傳播并行計算.因此，采用時空域高階有限差分法求解(10)式.盡管波場p和q都是四階偏微分方程組的解，但是我們將波場p看作壓力場，將q看作輔助波場時,式(10)可看作是兩個二階偏微分方程,與各向同性聲波方程類似，可通過高階差分方法求解(董良國，2000；劉紅偉，2010；李博，2010；劉文卿，2012,2013)，本文不再討論其解法.

3帶正則化校正的TTI介質(zhì)qP波方程推導(dǎo)

方程(10)即為Fletcher等(2009)提出的所謂的有限橫波方程計算(即剪切波速度不為0).從方程(10)出發(fā)，根據(jù)Alkhalifah的思想，將對稱軸方向剪切波速度設(shè)為0，可得到一個更為簡化的TTI介質(zhì)擬聲波方程，這個方程更容易求解，但數(shù)值求解也更加不穩(wěn)定(Fletcher et al.,2009)，Zhang Y和Zhang H Z(2009)采用偽譜法計算時，也發(fā)現(xiàn)了類似的不穩(wěn)定現(xiàn)象.當(dāng)然，我們可以對模型做平滑處理能使波場穩(wěn)定傳播，但這種方式同時會顯著地改變波的運動學(xué)特征，影響成像精度.值得注意的是，有限橫波方程也并非完全穩(wěn)定，這種不穩(wěn)定性通常開始于模型中對稱軸傾角存在劇烈變化的地方.本文提出了一種帶正則化項的有限橫波方法，數(shù)值試驗表明相較于有限橫波方程，這種方程能使波場更加精確穩(wěn)定地傳播下去.下面討論帶正則化項的有限橫波方程實現(xiàn)過程.

所謂帶正則化項方程，即是設(shè)計一個濾波器，對波場傳播中不穩(wěn)定的高波數(shù)成分做適當(dāng)衰減，以達到減少頻散或使波場計算穩(wěn)定的目的.首先來具體推導(dǎo)各向同性介質(zhì)中帶正則化項平面波解的形式，以形象說明正則化項的作用.對于構(gòu)建的帶正則化項的方程：

(11)

方程(11)兩邊同時對時間/空間做傅氏變化，經(jīng)過整理并帶入平面波解通解形式，可得到各向同性介質(zhì)下帶正則化項聲波方程平面波解形式：

(12)

其中k=(kx,ky,kz)，x=(x,y,z)分別為波數(shù)以及空間坐標向量,σ為正則化系數(shù)(通常取值(0～1)).與不帶正則化項的聲波方程平面波解相比，(12)式多出了一項振幅衰減項，其物理意義為，對平面波波場中傳播時間較長的高波數(shù)成分做衰減.利用這一思想，假設(shè)波場中不穩(wěn)定成分主要為高波數(shù)情況下，構(gòu)造TTI介質(zhì)中帶正則化項的有限橫波方程，使得波場數(shù)值計算更加穩(wěn)定.對于方程(10)可以構(gòu)建成帶正則化項的有限橫波方程：

(13)

4帶正則化校正的TTI介質(zhì)qP波方程的逆時偏移成像方法

帶正則化項的方程(13)看起來似乎要比原始的有限橫波方程(10)復(fù)雜許多，但在實際計算中，其增加的計算量和儲存量并不多.原因在于正則化項中旋轉(zhuǎn)坐標系下波場微分無需重復(fù)計算，僅需多儲存兩個波場，記錄上一時刻所有正則化項波場之和對應(yīng)的兩個波場，然后對時間求微分即可.實驗表明選擇合理的σ值，即使傾斜同相軸的方位及傾角變化較大，依然能夠保證波場穩(wěn)定地傳播.這也是這種方式實際資料計算中比較可行的原因.

(14)

其初始條件為：

(15)

其中s代表震源位置，f(t)為模擬震源函數(shù).在此RTM實現(xiàn)方案中，首先將邊界層最里側(cè)的波場記錄下來.二維情況下是四個邊界面處的波場記錄，三維情況下是六個邊界面處的三維波場.視機器內(nèi)存大小，可存儲在內(nèi)存中，也可存儲在硬盤上.

其次，在求解炮波場正向外推基礎(chǔ)上再求解波場沿時間逆向外推定解問題，如式(16)—(17).

(16)

其邊值條件為：

(17)

其中Bp(x,t)、Bq(x,t)分別是兩個波場在邊界處的值.利用每個時間層上的邊界波場，解上述邊值問題，沿時間逆時外推可實現(xiàn)炮波場逆向外推.此時，與檢波點波場的逆向外推相同步.因此在每個時間步，利用成像條件進行炮波場與檢波點波場零延遲相關(guān)成像.

5數(shù)值試算

為了測試采用帶正則化校正的TTI介質(zhì)qP波方程波場傳播的穩(wěn)定性,我們選擇了FoothillTTI模型中參數(shù)變化最為劇烈的部分作為測試模型.圖1展示了該模型的參數(shù)數(shù)據(jù)，模型中包含了速度、傾角及各向異性參數(shù).圖2展示的是不同控制方程不同時刻的波場快照,圖2a在正演方程中采取不帶正則化項的有限橫波方程，在模型中存在對稱軸傾角劇烈變化的地方，波場傳播不穩(wěn)定.模型正演中采用25Hz的雷克子波作為震源，震源放在(3550，0)m處,計算網(wǎng)格在兩個方向上均為10m.利用空間10階、時間2階有限差分去逼近控制方程中相應(yīng)導(dǎo)數(shù)項.圖2c、2d分別為采用帶正則化校正的TTI介質(zhì)qP波方程2000ms、3000ms的波場快照，與圖2a、2b同一時刻的波場快照對比可以看出，此時波的傳播是穩(wěn)定的.通過數(shù)值試算可看到，直接由頻散關(guān)系導(dǎo)出的TTI介質(zhì)中的有限橫波qP波方程在忽略參數(shù)的空間的導(dǎo)數(shù)時，由于受對稱軸傾角變化的影響仍然會產(chǎn)生數(shù)值求解不穩(wěn)定的現(xiàn)象，我們通過正則化的方式衰減高波數(shù)成分，最終使得整個波場計算穩(wěn)定.

我們將帶正則化校正的TTI介質(zhì)qP波方程方法應(yīng)用到TTI逆時偏移中,Foothill模型作為該測試所用模型，模型參數(shù)如圖1所示，波場傳播的計算均是采用高階有限差分法來求解,我們采用的是時間方向2階差分，空間方向10階差分的有限差分算法.圖3展示了逆時偏移成像結(jié)果，反射層都正確歸位，陡傾角地層都能成像.但在地層傾角較大、各向異性較強的區(qū)域，其下覆地層影響較大，各向同性及VTI各向異性逆時偏移成像效果都不理想，而TTI各向異性逆時偏移能使陡傾角地層及下覆地層準確成像.因此，對于各向異性介質(zhì)的地震數(shù)據(jù),采用TTI各向異性算法進行成像可得到更為精確的結(jié)果.

6實際資料試算與分析

本文所采用的實際數(shù)據(jù)是我國西部某油田的資料.該區(qū)發(fā)育大規(guī)模巖溶縫洞型儲層，非均質(zhì)性及各向異性都比較明顯.以往由于各向異性的影響，人們認為窄方位角資料效果好，如今方位各向異性現(xiàn)象逐步被人們認識，它可以提供更多的裂縫信息和儲層信息.因此，常規(guī)各向同性地震成像方法由于沒有考慮地層各向異性特性，很難準確刻畫縫洞的空間 位置.傳統(tǒng)的偏移距道集是反射點所有信息的綜合，并沒有準確包含傾角信息和方位角信息，無法將角度信息從道集中提取出來、也無法精確描述與方位有關(guān)的各向異性(王西文，2010).因此，只有全方位的角度道集才能獲得準確的全方位速度信息、傾角信息及與方位有關(guān)的各向異性參數(shù)信息，可以基于全方位共反射角道集進行各向異性速度建模,建模流程如圖4.

圖1　Foothill模型參數(shù)(a) vP速度； (b) 傾角模型； (c) epsilon模型； (d) delta模型.

圖2　不同方程傳播的波場快照(a) 有限橫波方程2 s波場； (b) 有限橫波方程3 s波場； (c) 正則化有限橫波方程2 s波場； (d) 正則化有限橫波方程3 s波場.

圖3　Foothill 模型不同方法逆時偏移結(jié)果對比(a) 各向同性逆時偏移結(jié)果； (b) VTI各向異性逆時偏移結(jié)果； (c) TTI各向異性逆時偏移結(jié)果.

圖5為傳統(tǒng)Kirchhoff方法獲得的共反射點(CRP)道集與全方位反射角道集對比.從道集特征可看出，常規(guī)的CRP道集是反射點所有信息的綜合，掩蓋了不同方位的速度變化特征及各向異性特征.全方位反射角道集定義了一個新的數(shù)據(jù)結(jié)構(gòu)(KorenandRavve,2011)，其角度為5°、10°、15°、20°、25°、30°、35°，每一個角度包含全方位的信息.從道集特征可看出，隨反射角度增大，相同反射角,不同方位道集不平直，具有不同的速度，其各向異性特征更加明顯.利用全方位反射角道集不同方位的剩余時差(圖5中綠線)可獲得有效的各向異性參數(shù)，通過網(wǎng)格層析成像方法建立精確的各向異性速度模型,如圖6.

圖7a為各向同性逆時偏移成像結(jié)果，圖7c為TTI各向異性逆時偏移成像結(jié)果，由于TTI各向異性逆時偏移成像考慮整個波場信息、地層的各向異性參數(shù)及傾角和方位角等空間屬性參數(shù)，所以串珠成像更聚焦、深度誤差更小，成像精度更高.從成像結(jié)果及鉆井誤差對比分析得出，各向同性逆時偏移結(jié)果深度誤差為190m(圖7a)，而各向異性逆時偏移結(jié)果深度誤差為30m(圖7b,圖7c),對于TTI各向異性逆時偏移成像由于上覆地層傾角不大，垂向消除深度誤差的效果與VTI介質(zhì)逆時偏移基本相 同，但在橫向上串珠位置相差12.5m，與鉆井更加吻合.同時深層成像TTI介質(zhì)各向異性逆時偏移成像相對VTI介質(zhì)及各向同性介質(zhì)的成像有較大的改善.

圖4　全方位角度域網(wǎng)格層析各向異性速度建模流程

7結(jié)論

對于復(fù)雜構(gòu)造成像，疊前逆時偏移是一個有力的工具.寬方位數(shù)據(jù)采集的普及使得各向異性逆時偏移變得必不可少.本文著重討論了TI介質(zhì)下的聲學(xué)近似、基于頻散關(guān)系的控制方程推導(dǎo)及高階有限差分數(shù)值解法.基于Christoffel方程從相速度/頻散關(guān)系出發(fā)推導(dǎo)出的qP方程，通過強制方程中沿對稱軸方向的剪切波速度為0，可以得到一個簡 化的波動方程.這個退化的方程可以提高求解的效率，并且在ε-δ≥0的VTI介質(zhì)中是穩(wěn)定的.若將其推廣到TTI介質(zhì)，在速度模型復(fù)雜時存在數(shù)值計 算的不穩(wěn)定問題.我們采用一個新的TTI介質(zhì)擬聲波方程及其解決方案，在有限橫波基礎(chǔ)上提出了一種帶正則化項的有限橫波方法，保證波場能精確穩(wěn)定地傳播下去，同時最大程度地減弱SV反射系數(shù)的影響.

圖5　全方位反射角道集(a)與傳統(tǒng)CRP道集(b)對比

圖6　各向異性參數(shù)剖面(a) epslion參數(shù)剖面； (b) delta參數(shù)剖面； (c)層速度參數(shù)剖面； (d) 地層傾角參數(shù)剖面； (e)方位角參數(shù)剖面.

圖7　各向同性逆時偏移(a)、VTI各向異性逆時偏移(b)與TTI各向異性逆時偏移(c)對比

FoothillTTI模型的正演和逆時偏移試驗表明，在波場傳播步數(shù)較多時，常規(guī)的有限橫波方程也開始出現(xiàn)計算不穩(wěn)定的現(xiàn)象，通過正則化的方式衰減高波數(shù)成分，最終使得整個波場計算穩(wěn)定，同時成像效果也最好.我們提出的方案在3DTTI介質(zhì)實際資料中也能得到非常好的效果.在模型參數(shù)劇烈變化的區(qū)域，上述方法仍然會產(chǎn)生SV波.為了進行P波偏移，這些SV波被看作是噪聲.但在實際資料偏移中，我們并沒有在最終偏移的剖面上發(fā)現(xiàn)來自這些SV能量的污染，可能是SV能量相對于P波能量比較弱的緣故，這些SV波能量是否在RTM角道集上有反映，還有待于進一步研究.從具有實際物理含義及數(shù)值計算穩(wěn)定角度出發(fā)，我們應(yīng)在彈性波框架下做近似，即TI介質(zhì)中qP控制方程的推導(dǎo)應(yīng)在時空域進行，這些是下一步研究的方向.

致謝感謝同濟大學(xué)王華忠老師、周陽博士在研究中給予的指導(dǎo)與幫助.

References

Alkhalifah T.1998.Acoustic approximations for processing in transversely isotropic media.Geophysics,63(2): 623-631.

Alkhalifah T.2000.An acoustic wave equation for anisotropic media.Geophysics,65(4): 1239-1250.

Baysal E,Kosloff D D,Sherwood J W C.1983.Reverse time migration.Geophysics,48(11): 1514-1524.

Chang W F,McMechan G A.1986.Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition.Geophysics,51(1): 67-84.

Chu C L,Nacy B K,Anno P D.2011.An accurate and stable wave equation for pure acoustic TTI modeling.∥ 81st Ann.Internat.Mtg.,Soc.Expl.Geophys.,Expanded Abstracts,179-184.

Dong L G,Ma Z T,Cao J Z,et al.2000.A staggered-grid high-order difference method of one-order elastic wave equation.Chinese J.Geophys.(in Chinese),43(3): 411-419.

Du X,Fletcher R P,Fowler P J.2008.A new pseudo-acoustic wave equation for VTI media.∥ 70thEAGE Conference & Exhibition.Rome,Italy.Duveneck E,Milcik P,Bakker P M,et al.2008.Acoustic VTI wave equations and their application for anisotropic reverse-time migration.∥ 78th Ann.Internat.Mtg.,Soc.Expl.Geophys.,Expanded Abstracts,2186-2190.

Duveneck E,Bakker P M.2011.Stable P-wave modeling for reverse-time migration in tilted TI media.Geophysics,76(2): S65-S75.

Fletcher R P,Du X,Fowler P J.2008.A new pseudo-acoustic wave equation for TI media.∥ 78thAnnual.Internationaln Meeting,SEG,Expanded Abstracts,2082-2086.

Fletcher R P,Du X,Fowler P J.2009.Reverse time migration in tilted transversely isotropic (TTI) media.Geophysics,74(6): WCA179-WCA187.

Koren Z,Ravve L.2011.Full-azimuth subsurface angle domain wavefield decomposition and imaging Part I: Directional and reflection image gathers.Geophysics,76(1): S1-S13.

Li B,Liu H W,Liu G F,et al.2010.Computational strategy of seismic pre-stack reverse time migration on CPU/GPU.Chinese J.Geophys.(in Chinese),53(12): 2938-2943,doi: 10.3969/j.issn.0001-5733.2010.12.017.

Liu F Q,Morton S A,Jiang S S,et al.2009.Decoupled wave equations for P and SV waves in an acoustic VTI media.∥ 79th Annual International Meeting,SEG,Expanded Abstracts,2844-2848.

Liu H W,Li B,Liu H,et al.2010.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation.Chinese J.Geophys.(in Chinese),53(7): 1725-1733,doi: 10.3969/j.issn.0001-5733.2010.07.024.

Liu W Q,Wang Y C,Yong X S,et al.2012.Prestack reverse time migration on GPU/CPU co-parallel computation.Oil Geophysical Prospecting (in Chinese),47(5): 712-716.

Liu W Q,Wang X W,Liu H,et al.2013.Application of velocity modeling and reverse time migration to subsalt structure.Chinese J.Geophys.(in Chinese),56(2): 612-625,doi: 10.6038/cjg20130225.McMechan G A.1983.Migration by extrapolation of time-dependent boundary values.Geophysical Prospecting,31(3): 413-420.

Sun R,McMechan G A,Lee C S,et al.2006.Prestack scalar reverse-time depth migration of 3D elastic seismic data.Geophysics,71(5): S199-S207.Thomsen L.1986.Weak elastic anisotropy.Geophysics,51(10): 1954-1966.

Wang X W,Liu W Q,Wang Y C,et al.2010.Research and application of common reflection angle domain pre-stack migration.Chinese J.Geophys.(in Chinese),53(11): 2732-2738,doi: 10.3969/j.issn.0001-5733.2010.11.021.

Whitmore N D.1983.Iterative depth migration by backward time propagation.∥53rdAnnual.International Meeting,SEG,Expanded Abstracts,382-385.

Zhang H Z,Zhang Y.2011.Reverse time migration in vertical and tilted orthorhombic media.∥ 81st Ann.Internat.Mtg.,Soc.Expl.Geophys.,Expanded Abstracts,185-189.

Zhang Y,Zhang H Z.2009.A stable TTI reverse time migration and its implementation.∥ 79th Ann.Internat.Mtg.,Soc.Expl.Geophys.,Expanded Abstracts,2794-2798.

Zhou H,Zhang G,Bloor R.2006.An anisotropic acoustic wave equation for VTI media.∥ 68th EAGE Conference & Exhibition.SPE,EAGE.

附中文參考文獻

董良國,馬在田,曹景忠等.2000.一階彈性波方程交錯網(wǎng)格高階差分解法.地球物理學(xué)報,43(3): 411-419.

李博,劉紅偉,劉國峰等.2010.地震疊前逆時偏移算法的CPU/GPU實施對策.地球物理學(xué)報,53(12): 2938-2943,doi: 10. 3969/j.issn.0001-5733.2010.12.017.

劉紅偉,李博,劉洪等.2010.地震疊前逆時偏移高階有限差分算 法及GPU實現(xiàn).地球物理學(xué)報,53(7): 1725-1733,doi: 10.3969/j.issn.0001-5733.2010.07.024.

劉文卿,王宇超,雍學(xué)善等.2012.基于GPU/CPU疊前逆時偏移研究及應(yīng)用.石油地球物理勘探,47(5): 712-716.

劉文卿,王西文,劉洪等.2013.鹽下構(gòu)造速度建模與逆時偏移成像研究及應(yīng)用.地球物理學(xué)報,56(2): 612-625,doi: 10.6038/cjg20130225.

王西文,劉文卿,王宇超等.2010.共反射角疊前偏移成像研究及應(yīng)用.地球物理學(xué)報,53(11): 2732-2738,doi: 10.3969/j.issn.0001-5733.2010.11.021.

(本文編輯何燕)

基金項目國家油氣專項“天然氣地球物理烴類檢測評價技術(shù)及應(yīng)用(2016ZX05007-006)資助.

作者簡介劉文卿，男，1977年生，在讀博士研究生，2000年于成都理工大學(xué)獲應(yīng)用地球物理學(xué)士學(xué)位.主要從事地震速度建模與地震成像方法研究工作.E-mail：liuwq@petrochina.com.cn

doi:10.6038/cjg20160326 中圖分類號P631

收稿日期2014-11-30，2015-12-28收修定稿

A regularized qP-wave equation for TTI media and its application to reverse time migration

LIU Wen-Qing1,2,WANG Xi-Wen2,WANG Yu-Chao2,YONG Xue-Shan2,WANG Xiao-Wei2,ZHANG Tao2

1StateKeyLaboratoryofOilandGasReservoirGeologyandExploitation,ChengduUniversityofTechnology，Chengdu610059,China2ResearchInstituteofPetroleumExploration&Development-Northwest,Petrochina,Lanzhou730020,China

AbstractThe research of anisotropy is significant for subsurface precise imaging.With rapid development of computation capability and gradual generalization of seismic acquisition with wide azimuth,anisotropy should be taken into consideration.The reverse time migration is based on the accurate solution of two-way wave equation,therefore compared with other methods,it has many advantages,such as no dip limitation and the ability to image turning waves and multiples.Scalar wave equation could be used to calculate the wave field in isotropic media.In anisotropic media,however,P and SV waves are coupled and there is no pure scalar wave mode.Consequently,pseudo acoustic wave (qP-waves) which can represent the kinematic property of P-component in the coupled wave field are usually used for imaging in anisotropic media.

In this paper,starting with the P-SV coupled dispersion relation for VTI media,we derive qP wave equation for TTI media,which can be solved using an explicit finite difference format.With the help of acoustic approximation,stable numerical solutions can be obtained for a model with vertical symmetric axis,ε≥δ and,if the velocity of shear waves becomes zero along the symmetric axis.However,acoustic approximation will cause the instability of wave field propagation and numerical computation as anisotropic parameters vary in the direction of symmetric axis in TTI media.To solve this problem,we proposed to regularize the finite shear wave equation by adding a regularized term to the original equation.The idea of our method is similar to designing a filter which could attenuate the high wave number components appropriately.We give out the specific form of regularized TTI qP wave equation.Then,the TTI reverse time migration (RTM) method using the regularized wave equation was discussed.In numerical implementation,the computing complexity of our regularized wave equation only increases by one wave field differential and the storage needed only increases by two more wave field.Numerical tests show that by choosing appropriate σ,we can obtain stable wave field even with sharp variation of azimuth and dip of tilted axis.

For numerical examples,we give out the modelling and RTM results on the Foothill TTI model and the imaging result of field data from a carbonate rock area.When prorogating time of the wavefield is much longer,the traditional finite shear wave equation shows wavefield instability.In contrast,our method can obtain wavefield without any instabilities.In areas with abrupt changes of model parameters,our method may still produce weak SV wave which are viewed as artificial noise in reverse time migration.But for real seismic data,we didn’t find any of these SV wave energy in the final migration profile which may be caused by the quite weaker SV energy compared with P wave energy.The RTM imaging with high quality for field seismic data demonstrates that our method works well and can be applied to field data effectively.

KeywordsTTI media; Anisotropy; Wave field coupling; Finite shear wave; Regularization; RTM

劉文卿，王西文，王宇超等.2016.帶正則化校正的TTI介質(zhì)qP波方程及其逆時偏移方法和應(yīng)用.地球物理學(xué)報，59(3):1059-1069,doi:10.6038/cjg20160326.

Liu W Q,Wang X W,Wang Y C,et al.2016.A regularized qP-wave equation for TTI media and its application to reverse time migration.Chinese J.Geophys.(in Chinese),59(3):1059-1069,doi:10.6038/cjg20160326.