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

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

?

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

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

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

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

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

1引言

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

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

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

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

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

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

(1)

(2)

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

(3)

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

×(p-q).

(4)

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

(5)

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

(6)

(7)

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

(8)

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

(9)

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

(10)

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

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

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

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

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

(11)

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

(12)

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

(13)

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

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

(14)

其初始條件為：

(15)

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

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

(16)

其邊值條件為：

(17)

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

5數(shù)值試算

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

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

6實(shí)際資料試算與分析

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

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

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

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

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

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

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

7結(jié)論

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

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

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

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

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

致謝感謝同濟(jì)大學(xué)王華忠老師、周陽(yáng)博士在研究中給予的指導(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.

附中文參考文獻(xiàn)

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

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

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

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

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

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

(本文編輯何燕)

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

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

doi:10.6038/cjg20160326 中圖分類號(hào)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波方程及其逆時(shí)偏移方法和應(yīng)用.地球物理學(xué)報(bào)，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.