王璐琛，常旭，王一博

1 中國科學(xué)院頁巖氣與地質(zhì)工程重點實驗室，中國科學(xué)院地質(zhì)與地球物理研究所，北京　100029　2　中國科學(xué)院大學(xué)，北京　100049

?

TTI介質(zhì)的交錯網(wǎng)格偽P波正演方法

王璐琛1,2，常旭1*，王一博1

1 中國科學(xué)院頁巖氣與地質(zhì)工程重點實驗室，中國科學(xué)院地質(zhì)與地球物理研究所，北京1000292中國科學(xué)院大學(xué)，北京100049

摘要研究了三維弱各向異性近似下，利用偽P波(偽縱波)模擬彈性波場P分量在傾斜對稱軸的橫向各向同性(TTI)介質(zhì)中的傳播過程，并對比了分別基于彈性Hooke定律、彈性波投影和運動學(xué)色散方程所建立的三種二階差分偽P波方程的正演特點.目前這些偽P波方程數(shù)值計算主要采用規(guī)則網(wǎng)格差分，但是規(guī)則網(wǎng)格在TTI模擬中有低效率、低精度以及不穩(wěn)定的缺點.為了提高計算的精度，本文構(gòu)建出相應(yīng)方程的交錯網(wǎng)格有限差分格式.通過對比偽P波方程在三維TTI介質(zhì)中不同的數(shù)值模擬的表達形式，本文認(rèn)為基于色散方程所建立的偽P波方程在模擬彈性波中P波傳播的過程中具有最小的噪聲.本文分析不同的各向同性對稱軸空間角度的頻散特征，并引入適當(dāng)?shù)臋M波速度維持計算的穩(wěn)定.二維模型算例表明，本文提出的交錯網(wǎng)格正演算法可以得到穩(wěn)定光滑的偽P波正演波場.使用本文交錯網(wǎng)格算法對二維BP TTI模型的逆時偏移也具有較穩(wěn)定的偏移結(jié)果.

關(guān)鍵詞偽P波； 傾斜各向異性； 交錯網(wǎng)格

1引言

目前油氣開采中重點研究的頁巖氣主要聚集吸附在致密頁巖層或泥巖層中.研究表明頁巖的各向異性可達70%(Vernik and Nur,1992; Johnston and Christensen,1995)，因此使用傾斜對稱軸的橫向各向同性(TTI)介質(zhì)(Qin et al.,2003a; 2003b; 吳國忱等,2010; 黃翼堅等,2011)，能更好地描述地震波在頁巖中的傳播規(guī)律.頁巖氣的開采需要利用被動源微地震進行監(jiān)測.被動源地震的震源不同于主動源地震，其震源存在多種破裂形式.這些破裂形式受到所處地下環(huán)境的介質(zhì)各向異性所影響，尤其是在最近研究比較熱門的微地震領(lǐng)域，其震源基本都是處在很明顯的TTI沉積層介質(zhì)中.研究震源處介質(zhì)的各向異性有助于了解地震區(qū)域震源機制的構(gòu)成.

由于彈性波的傳播本身已經(jīng)比較復(fù)雜，介質(zhì)的各向異性又會增加非零的Hooke彈性參數(shù)從而加重這種復(fù)雜的狀況，會大幅度增加波場正演的計算時間以及計算量.同時考慮到縱波勘探應(yīng)用廣泛，技術(shù)也更加成熟，是比較常用的研究方法.因此，只研究類似P波的波場傳播(稱為偽P波)，無疑是在各向異性研究中一種比較折中的方式(Ekkehart,2010; Yoon et al.,2010).針對偽P波的波場傳播，可以通過限定不同的約束條件(Alkhalifah,1998; Zhou et al.,2006; Fletcher et al.,2008)或不同的分離方式從TTI彈性波方程中建立出偽P波的傳播方程(Fletcher et al.,2009; Zhan et al.,2011; Kang and Cheng,2011).

這些方程普遍都是復(fù)雜的二階微分方程，所以在數(shù)值實現(xiàn)上基本采用最直觀的離散方式，用規(guī)則網(wǎng)格差分將二階方程進行數(shù)值離散.然而規(guī)則網(wǎng)格由于差分精度有限，一方面不利于在有限的存儲與計算時間下提高計算精度；另一方面在各向異性參數(shù)變化較為劇烈的區(qū)域，由于差分網(wǎng)格點間距過大，容易造成計算的不穩(wěn)定性.同時有些規(guī)則網(wǎng)格算法會將各向異性參數(shù)的微商引入導(dǎo)數(shù)項的計算中，這樣容易因各向異性參數(shù)微商計算的精度不夠而造成數(shù)值計算發(fā)散.

由于規(guī)則網(wǎng)格差分在TTI模擬中存在如上弊端，構(gòu)造出結(jié)構(gòu)簡約的交錯網(wǎng)格差分格式就顯得尤為必要.首先，交錯網(wǎng)格本身的計算精度就要高于同樣差分階數(shù)的規(guī)則網(wǎng)格，另外由于交錯網(wǎng)格是在半網(wǎng)格點上進行差分的，能更好的適應(yīng)各向異性的復(fù)雜變化，所以數(shù)值計算的穩(wěn)定性也能提高.

本文主要研究在弱各向異性(Thomsen,1986)近似下，地震波傳播交錯網(wǎng)格差分的實現(xiàn)方法，分析了TTI介質(zhì)中交錯網(wǎng)格相比于規(guī)則網(wǎng)格所具有的高效、低存儲等計算優(yōu)勢，建立了統(tǒng)一的交錯網(wǎng)格偽P波差分方式.本文對比了三種主流的偽P波正演方程，并分別建立了交錯網(wǎng)格的離散公式.本文的交錯網(wǎng)格差分方法將各向異性的旋轉(zhuǎn)特性加載在空間導(dǎo)數(shù)上，這樣一方面保留了各向異性的物性參數(shù)，便于提取物理屬性，方便研究地層的巖性(杜啟振和秦童,2009)，另一方面避免了交叉導(dǎo)數(shù)項的出現(xiàn)，簡化了計算，提高了計算穩(wěn)定性.本文在研究TTI介質(zhì)的性質(zhì)時是先假定TTI傾角(θ)為零，然后再推廣到任意的空間轉(zhuǎn)角，因此本文交錯網(wǎng)格差分方程的適應(yīng)性比較好.

2三維TTI介質(zhì)的波場傳播理論

地震波在各向異性介質(zhì)中的傳播是一個復(fù)雜的彈性波傳播過程，各向異性會使得Hooke彈性系數(shù)矩陣變得更加復(fù)雜，同時也不利于波場模擬.由于彈性各向異性介質(zhì)波動方程的復(fù)雜，所以大多數(shù)各向異性正演和偏移方法通常都是在Thomsen弱各向異性的近似下使用偽P波方程來實現(xiàn).

2.1TTI介質(zhì)的弱各向異性

彈性波波動方程的一般數(shù)學(xué)表示為

(1)

其中，u通常代表質(zhì)點的位移場，ρ和F分別表示介質(zhì)的密度和所受到的外力.cijkl表示4階Hooke定律中的剛度矩陣，表示應(yīng)力與應(yīng)變的對應(yīng)關(guān)系.由于應(yīng)力的對稱性，所以空間上9個應(yīng)力分量中只有Pxx,Pyy,Pzz,Pxy,Pyz,Pzx，這6個是獨立的，分別表示3個正應(yīng)力以及3個剪切應(yīng)力.又由于彈性系數(shù)張量的對稱性，所以剛度矩陣中一共只有21個是獨立的變量.剛度矩陣具體的形式如下：

(2)

如果是橫向各向同性的情況，彈性剛度矩陣則化簡如公式(3)所示：

(3)

式(3)中的彈性參數(shù)，無法直接對應(yīng)物理量，因此按照Thomsen引入的5個弱各向異性彈性參數(shù)(Thomsen,1986)進行如下替換，可賦予各參數(shù)可測的物理意義：

(4)

(5)

2.2弱各向異性下的P波運動學(xué)色散方程

弱各向異性下的偽P波方程，并不是簡單的在彈性波動方程中令橫波波速為0.因為各向異性本身就是對彈性的一種描述，不相等的Thomsen弱各向異性參數(shù)(例如ε≠δ，此時稱為非橢圓各向異性)(Thomsen,1986)會導(dǎo)致方程在不同方向的傳播上出現(xiàn)不同的橫波假象.所以如設(shè)定了橫波波速為0，還是會產(chǎn)生偽橫波的出現(xiàn)(Tsvankin et al.,2001).同時根據(jù)Alkhalifah提出的基于橫波速度為0的縱波各向異性方程(Alkhalifah,1998; 2000)，被假定為0的也僅僅是沿著TTI介質(zhì)中各向同性對稱軸方向的橫波速度.雖然縱波的傳播在運動學(xué)上基本滿足各向異性傳播性質(zhì)，但是不同的推導(dǎo)方式還是會有不同的振幅以及頻散表現(xiàn)，尤其是在不同的TTI對稱軸的空間角度上，表現(xiàn)會更加劇烈.

首先研究VTI介質(zhì)的情況.在各向異性的情況下，波的相速度能夠描述波在介質(zhì)中傳播的特點，因此可以從相速度的角度建立波場的運動學(xué)傳播公式.先求解Christoffel方程可以得到精確的P-SV波和SH波相速度公式(Dellinger and Etgen,1990; Tsvankin.,1996;2001)，其中P-SV波的相速度公式如下：

(6)

(7)

2.3偽P波方程組構(gòu)造

(1)依據(jù)彈性Hooke定律：由公式(3)和(5)所表示的彈性Hooke方程進行變量代換(Zhang et al.,2011)，設(shè)

(8)

其中，ux，uy，uz表示三維空間三個方向上的位移.將公式(8)代入TTI彈性Hooke方程就可以得到如下的偽P波時域方程組：

(9)

(10)

(11)

那么就可以得到如下的運動學(xué)偽P波時域方程組(Zhang et al.,2011; Fletcher et al.,2009; Zhan et al.,2011; 張巖和吳國忱,2013)：

(12)

這里要強調(diào)的是，之所以稱是偽P波，因為簡單的令沿著TTI對稱軸的橫波速度為0并不能將介質(zhì)其他部分的橫波相速度消除(Grechka et al.,2004).當(dāng)然可以通過輔助波場壓制偽橫波傳播，或者通過震源處理(Duveneck et al.,2008)減小偽橫波的影響.

(13)

(14)

這樣就得到了TTI的偽P波方程.其中包含一個通常意義上的聲波壓力場和一個計算引入的輔助場，本質(zhì)上是對TTI彈性波方程的一種計算上的簡化.

2.4偽P波方程交錯網(wǎng)格差分格式

公式(9)、(10)、(12)均給出的是二階位移方程，按照公式格式進行數(shù)值離散即可得到相應(yīng)的VTI數(shù)值方程，進一步借由公式(14)旋轉(zhuǎn)得到TTI數(shù)值方程，這樣就可以實現(xiàn)公式的規(guī)則網(wǎng)格差分格式.這樣處理的一個好處是，從VTI過度到TTI所需的空間旋轉(zhuǎn)計算被加載在了變量的空間二階微分上，避免了對各向異性參數(shù)的空間旋轉(zhuǎn)插值，提高了計算的穩(wěn)定性以及精確性，缺點就是二階微分中包含了交叉導(dǎo)數(shù)項，不利于高精度的計算.為了提高計算精度，現(xiàn)在考慮如何實現(xiàn)公式(9)、(10)、(12)的應(yīng)力、速度交錯網(wǎng)格差分格式.這里以二維情況的公式(12)為例說明.二維方程與空間導(dǎo)數(shù)的旋轉(zhuǎn)形式如公式(15)和(16)所示：

(15)

(16)

對比聲波方程交錯網(wǎng)格格式，記速度場vpx,vpz,vrx,vrz滿足

(17)

(18)

(19)

(20)

(21)

速度場v rx,v rz空間導(dǎo)數(shù)計算同理為

(22)

代入公式(22)到公式(19)中，得旋轉(zhuǎn)后的速度場空間導(dǎo)數(shù)差分計算為

(23)

所以公式(18)所對應(yīng)的應(yīng)力方程差分格式由上式代入可得

(24)

公式(20)和公式(24)合在一起就是公式(12)所對應(yīng)的方程的一階速度應(yīng)力方程交錯網(wǎng)格差分格式.我們很容易發(fā)現(xiàn)，公式(12)的方程形式僅僅體現(xiàn)在交錯網(wǎng)格差分格式的應(yīng)力方程(24)中，這樣的目的在于可以最大限度的簡化速度場方程的差分形式，同時也易于擴展到其他形式的偽P波方程.因此，公式(9)和(10)所對應(yīng)方程的交錯網(wǎng)格速度方程同公式(20)完全一樣，應(yīng)力方程同理公式(24)其格式分別如下:

(25)

和

(26)

這樣就構(gòu)建出了三種方法所對應(yīng)方程(公式(9)(10)(12))的交錯網(wǎng)格差分格式.其中對TTI對稱軸的空間旋轉(zhuǎn)效應(yīng)的計算僅加載在了應(yīng)力方程中的速度場空間導(dǎo)數(shù)差分上.

3偽P波在TTI介質(zhì)中的數(shù)值計算

3.1偽P波的正演

先考察TTI介質(zhì)中不同正演差分策略的計算效率.以二維均勻模型正演計算為例，計算區(qū)域網(wǎng)格數(shù)為400×400，正演步數(shù)為3000步，對偽P波交錯網(wǎng)格、偽P波規(guī)則網(wǎng)格兩種正演策略各進行40次計算，結(jié)果顯示于表1.

作為對比，對同樣大小的速度模型采用聲波交錯網(wǎng)格和聲波規(guī)則網(wǎng)格兩種正演策略各進行相同的40次正演計算，結(jié)果也顯示于表1.

表1　二維模型下不同正演方式的內(nèi)存需求和計算效率

從表1可以看出，對于偽P波的正演模擬，交錯網(wǎng)格方法無論是在內(nèi)存需求上還是計算時間上都要少于規(guī)則網(wǎng)格算法，這一點同傳統(tǒng)聲波方程的有限差分是不同的.對于聲波方程而言，交錯網(wǎng)格的計算精度比同階差分的規(guī)則網(wǎng)格精度要高，但是計算效率要慢，這是因為一階速度應(yīng)力方程比二階聲波方程需要更多的差分計算.由于TTI偽P波的規(guī)則網(wǎng)格正演方程需要計算二階空間微分的交叉導(dǎo)數(shù)，因此對每個網(wǎng)格點需要額外的一次循環(huán)計算(先計算一個方向的一階差分，再計算另一個方向上的差分)，并需要引入相應(yīng)的交叉導(dǎo)數(shù)臨時變量.這樣即增加了額外計算量導(dǎo)致總的計算效率下降，以至于平均計算時間反而超過TTI偽P波交錯網(wǎng)格，又增大了內(nèi)存需求(三維情況下計算的交叉導(dǎo)數(shù)項更多，對應(yīng)需要更多的臨時變量).另外由于交錯網(wǎng)格差分精度比同階的規(guī)則網(wǎng)格精度要高，因此TTI偽P波交錯網(wǎng)格算法相比于規(guī)則網(wǎng)格算法有著更優(yōu)越的應(yīng)用效果.

再考察TTI不同正演方程的模擬效果.通常情況下降階得出的二元二階耦合方程中，輔助波場可以不具有實質(zhì)上的物理含義，也就是說在此情況下并不會影響到P波運動學(xué)的特性.而如公式(10)的方程(Kang and Cheng,2011)，又是從各向異性彈性波運動方程中推導(dǎo)的P波方程，分別表示兩個應(yīng)力方向上的波場傳播(二維情況下為沿著傾斜對稱軸的方向和垂直于傾斜對稱軸的方向)，其波場快照如圖1a所示(使用了聲波近似條件，即vSz=0).其中P波速度vPz=2500 m·s-1，各向異性參數(shù)和TTI對稱軸的空間角度分別是ε=0.24，δ=0.1，θ=45°.將其二者相加組合就是完整的P波波場.

可以看到，在圖1中各個分量波場中((a1—c1)和(a2—c2)圖)內(nèi)部尖角所表示的偽S波正好呈現(xiàn)振幅相反的特點，相加后的波場中(a3—c3圖)，最外層的P波成分變得更加完整清晰，同時內(nèi)部的偽S波尖角也得到了衰減.所以這里就將耦合公式中的兩個分量進行相加來表示模擬的偽P波波場.

將此處理應(yīng)用到公式(9)，公式(12)所表示的三種偽P波耦合方程，二維情況下的各自的波場快照對比結(jié)果如圖1B和圖1C，其中P波速度vPz=2500 m·s-1，各向異性參數(shù)和TTI對稱軸的空間角度分別是ε=0.24，δ=0.1，θ=45°，同樣使用聲波近似條件，即令vSz=0.

對比圖1B、1C和圖1A，不難看出，不管采用何種公式，輔助波場r都體現(xiàn)著更多的內(nèi)部偽S波成分，對于外圈的P波成分，或者沒有實際意義(圖1B)，或者起到補償消減的作用(圖1C).所以說明，將兩個波場組合相加，可以提高偽P波波場的精確性.由于兩個波場所表示的主體不同，輔助場主要體現(xiàn)的是次生的偽S波場成分，所以對這兩個波場就可以進行各自的濾波處理，從而進一步壓制偽S波的成分(張巖和吳國忱,2013).

分析圖1以及對應(yīng)的公式(公式9、10、12)，可 以看出，圖1A(公式10)是按照波場運動的方向進行波場分解的，兩個分量場表示偽P波波場兩個方向的運動分量.圖1B(公式9)是按照動力學(xué)應(yīng)力場的方向分量進行波場分解的，兩個分量場表示偽P波波場的兩個方向的應(yīng)力分量.其波場反映著動力學(xué)的振幅信息.圖1C(公式12)是按照波場的P、SV成分進行分解的，兩個分量場分別表示P波場和SV波場.圖2表示的是三個偽P波公式(公式9、10、12)正演合成波場的振幅和頻散對比：

圖1　三種正演公式模擬的波場快照(A)公式(10)： (B)公式(9)； (C)公式(12)； (a1—c1)p成分; (a2—c2)r成分; (a3—c3)p+r成分

從圖2A中可以明顯看出公式(10)與公式(12)對應(yīng)的偽P波波場是比較相似的，而基于應(yīng)力的公式(9)對應(yīng)的綠色實線同另外兩個相差較大，無論是波形的形狀還是到時都存在較明顯的差異.它反映了TTI介質(zhì)中質(zhì)點的應(yīng)力場同速度場之間的差異.圖2B中可看出，頻散主要集中在波場中間的偽橫波尖角處.并且公式(9)的頻散明顯要強于另外兩個公式正演的波場.所以本文之后的模型正演部分，采用的是頻散較小的公式(12)構(gòu)造的偽P波方程.同時采用有約束的橫波速度(Tsvankin et al.,2001)(不是簡單設(shè)為0)，以及特征投影濾波的方式(Zhou and Zhang,2009)壓制正演中的偽橫波噪聲，可以進一步降低數(shù)值頻散，提高計算的穩(wěn)定性.設(shè)vPz=2500 m·s-1，ε=0.24，δ=0.1，圖3展示了偽橫波數(shù)值頻散的壓制效果.其中外圈橢圓為正演的偽P波波前，內(nèi)部波場為偽橫波噪聲.

除了不同降階公式本身的影響，TTI介質(zhì)中偽P波的傳播還會受到對稱軸角度的影響.這種影響主要是源于數(shù)值計算上產(chǎn)生的頻散現(xiàn)象，角度越傾斜，頻散的現(xiàn)象越嚴(yán)重，尤其是在三維模擬的情況下.令θ表示TTI對稱軸相對于豎直的傾角，φ表示TTI對稱軸的方位角，圖4所示的是對稱軸在不同空間轉(zhuǎn)角下的波場垂直切片.

可以看到，在空間轉(zhuǎn)角為0°的時候，數(shù)值頻散現(xiàn)象是非常弱的，但是隨著轉(zhuǎn)角的出現(xiàn)，在內(nèi)部偽S波尖角的地方還是出現(xiàn)了頻散現(xiàn)象，右圖的空間角度更加復(fù)雜，頻散的現(xiàn)象也比中間的要更加復(fù)雜.

3.2模型正演與偏移結(jié)果

首先考察二維情況下本文所述方法的正演效果.這里作為簡化，我們只研究固定TTI對稱軸傾角的情況.所用模型如圖5所示.

模型縱向深2000 m，橫向長4000 m，網(wǎng)格大小為4 m，其中TTI各向異性僅僅存在于中間的兩層，也就是說最上層和最下層是普通的各向同性介質(zhì)，中間兩層中TTI的對稱軸角度設(shè)定為固定值 45°.此模型表示一個斜插入的各向異性層，并且該層內(nèi)還具有構(gòu)造起伏的特征.

圖2　不同公式合成波場的道記錄(A)及不同公式合成波場的頻散(B)(a1) 第101道;(a2) 第52道; (b1) 公式(9);(b2) 公式(10); (b3) 公式(12).

圖3　二維TTI正演波場中的偽橫波(a) 橫波速度vPz為0; (b) 橫波速度vPz不為0; (c) 橫波速度vPz不為0并使用濾波.

圖4　三維空間中對稱軸不同轉(zhuǎn)角下的垂向波場(a) 轉(zhuǎn)角:θ=0°,φ=0°; (b) 轉(zhuǎn)角:θ=45°,φ=0°; (c) 轉(zhuǎn)角:θ=45°,φ=30°.

圖5　TTI各向異性參數(shù)模型(a) P波速度vPz； (b) δ； (c) ε； (d) θ.

采用偽P波模擬的單一地下震源激發(fā)的不同時刻地震波的波場快照如圖6所示.其震源位置為橫向2000 m，深度1000 m處.圖6中從左至右的波場時刻依次是0.2、0.3、0.4、0.6 s.

從圖6B中可以發(fā)現(xiàn)，在非各向異性介質(zhì)中，輔助波場直接衰減為0，從另一方面也說明了本文中使用的偽P波方程的輔助場描述的是偽S波的產(chǎn)生與傳播.因此也就是說，合成波場在地表的檢波器處同偽P波波場是相同的，波場合成并不影響最終的道記錄.如果說存在影響的話，影響也僅會發(fā)生在TTI介質(zhì)的井中接收情況中.此地下震源對應(yīng)的偽P波單炮記錄如圖7所示，從圖中可見模擬的地震記錄波形穩(wěn)定，并且沒有計算頻散與發(fā)散.說明我們的正演效果較好.

對于三維情況，速度模型與各向異性參數(shù)模型如圖8所示(只顯示非恒定值的量)，采用前面圖5的二維4層速度模型在Y方向上平推而得，因此三維模型的中間兩層TTI的對稱軸角度為固定值θ=45°，φ=0°.地面觀測系統(tǒng)如圖9a所示，用水平地表“米”字型共4條測線進行觀測.模型大小為X方向4000 m，Y方向800 m，深度2000 m，網(wǎng)格大小4 m.

假定震源位置為X方向2000 m，Y方向400 m，深度1000 m處.此震源事件的單炮記錄如圖9b所示.從圖9b中可見，三維TTI情況下，本文正演的波場也具有良好的數(shù)值穩(wěn)定性.

最后再用標(biāo)準(zhǔn)TTI模型的逆時偏移(RTM)來進一步測試本文的正演算法.RTM偏移采用炮點沿時間正向延拓震源子波，檢波點同時沿時間反向延拓地震記錄的方式，得到炮點與檢波點波場.接著利用成像條件與這兩個波場對地下反射點進行成 像.其中，波形正確、穩(wěn)定、頻散小的正演差分求解，是RTM偏移成像的基礎(chǔ).BP TTI模型是國際上普遍采用的二維各向異性模型，其模型網(wǎng)格大小為6.25 m×6.25 m，橫向距離78.7 km，最大深度11.25 km，我們用模型橫向0～51 km的部分進行計算.其速度的模型與各向異性參數(shù)模型如圖10所示，使用發(fā)布的數(shù)據(jù)以及本文所述的正演方式對其進行偏移的結(jié)果如圖11所示.

圖6　(A) 單一震源的合成波場(p+r)和(B)輔助波場(r)在不同時刻的波場從左至右時間依次是：0.2 s，0.3 s，0.4 s，0.6 s.

圖7　單震源炮記錄

圖8　(a) 三維P波vPz速度模型； (b) TTI各向異性模型δ； (c) TTI各向異性模型ε

圖9　三維TTI模型的地表觀測系統(tǒng)(a)及對應(yīng)的地震事件的炮記錄(b)

圖10　BP TTI模型的速度與各向異性參數(shù)(a) vPz； (b) ε； (c) δ； (d) θ.

圖11　使用本文正演方式的逆時偏移(RTM)結(jié)果

圖11偏移結(jié)果同相軸連續(xù)性較好，并且偏移結(jié)果的鹽丘位置同給定模型中的鹽丘位置能相互匹配，說明本文的正演算法能正確的模擬TTI介質(zhì)中地震波波形.另外，成像結(jié)果無發(fā)散，說明本文的正演差分算法具有良好的穩(wěn)定性.

4結(jié)論

在現(xiàn)階段對復(fù)雜地質(zhì)的地震資料的處理和解釋過程中，忽略各向異性將會導(dǎo)致很大的偏差，無法得到精度更高的結(jié)果.介質(zhì)的各向異性，尤其是更貼近實際情況的TTI介質(zhì)，能夠顯著地改善資料處理的質(zhì)量.對于TTI介質(zhì)的研究，首要的任務(wù)是估計正確的各向異性參數(shù)，然后是使用恰當(dāng)?shù)姆绞侥M波場的傳播.本文所使用的偽P波方程能夠簡化TTI波場的復(fù)雜性，并能較快捷的應(yīng)用到三維的成像、偏移中.這其中要注意的是，不同的偽P波方程具有不盡相同的物理意義，同時也有不同的頻散現(xiàn)象.本文對比了基于位移場和應(yīng)力場的偽P波方程各自的特點，推導(dǎo)并實現(xiàn)了正演模擬三維各向異性介質(zhì)中P波成分的偽P波交錯網(wǎng)格有限差分算法.從正演的波場上看，基于色散方程降階分離出的偽P波波場頻散較小，正演模擬效果最好.最后，由于頻散在TTI介質(zhì)中的影響會遠遠大于各向同性介質(zhì)，因此如何有效的去除頻散，認(rèn)識頻散在不同方程條件下的物理意義，仍將是后續(xù)研究的重點.

致謝感謝國家自然科學(xué)基金資助(41230317、41274112)，感謝英國石油公司(BP)提供的二維TTI理論模型與其合成數(shù)據(jù).

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.

Dellinger J,Etgen J.1990.Wave-field separation in two-dimensional anisotropic media.Geophysics,55(7): 914-919.Du Q Z,Qin T.2009.Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium.Chinese J.Geophys.(in Chinese),52(3): 801-807.

Duveneck E,Milcik P,Bakker P M,et al.2008.Acoustic VTI wave equations and their application for anisotropic reverse-time migration.∥ 78th Annual International Meeting,SEG,Expanded Abstracts,2186-2190.

Ekkehart T.2010.Reverse-time migration in TTI media.∥ 80th Annual International Meeting,SEG,Expanded Abstracts,3193-3197.

Fletcher R,Du X,Fowler P J.2008.A new pseudo-acoustic wave equation for TI media.∥ 78th Annual International 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.

Fowler P J,Du X,Fletcher R P.2010.Coupled equations for reverse time migration in transversely isotropic media.Geophysics,75(1): S11-S22.Grechka V,Zhang L B,Rector III J W.2004.Shear waves in acoustic anisotropic media.Geophysics,69(2): 576-582.

Hao C T,Yao C,Wang X.2006.The characteristics of velocities with azimuth variation for arbitrary spatial orientation TI media.Progress in Geophysics (in Chinese),21(2): 524-530.

Huang Y J,Zhu G M,Liu C Y.2011.An approximate acoustic wave equation for VTI media.Chinese J.Geophys.(in Chinese),54(8): 2117-2123,doi: 10.3969/j.issn.0001-5733.2011.08.019.

Johnston J E,Christensen N I.1995.Seismic anisotropy of shales.J.Geophys.Res.,100(B4): 5991-6003.

Kang W,Cheng J B.2011.Prestack scalar reverse time migration of elastic seismic data in TI media.∥ 81th Annual International Meeting,SEG,Expanded Abstracts,3367-3371.

Qin Y L,Zhang Z J,Li S L.2003a.CDP mapping in tilted transversely isotropic (TTI) media.Part I: Method and effectiveness.Geophysical Prospecting,51(4): 315-324.

Qin Y L,Zhang Z J,Xu S H.2003b.CDP mapping in tilted transversely isotropic (TTI) media.Part II: Velocity analysis by combining CDP mapping with a genetic algorithm.Geophysical Prospecting,51(4): 325-332.

Thomsen L.1986.Weak elastic anisotropy.Geophysics,51(10): 1954-1966.

Tsvankin I.1996.P-wave signatures and notation for transversely isotropic media: An overview.Geophysics,61(2): 467-483.

Tsvankin I.2001.Seismic Signatures and Analysis of Reflection Data in Anisotropic Media.Amsterdam: Pergamon.

Vernik L,Nur A.1992.Ultrasonic velocity and anisotropy of hydrocarbon source rocks.Geophysics,57(5): 727-735.

Wu G C,Liang K,Yin X Y.2010.The analysis of phase velocity and polarization feature for elastic wave in TTI media.Chinese J.Geophys.(in Chinese),53(8): 1914-1923,doi: 10.3969/j.issn.0001-5733.2010.08.017.

Yoon K,Sang S,Jean J,et al.2010.Stability and speedup issues in TTI RTM implementation.∥ 80th Annual Meeting,SEG,Expanded Abstracts,3221-3225.

Zhan G,Pestana R C,Stoffa P L.2011.An acoustic wave equation for pure P wave in 2D TTI media.∥ 81th Annual International Meeting,SEG,Expanded Abstracts,168-173.

Zhang Y,Wu G C.2013.Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration.Chinese J.Geophys.(in Chinese),56(6): 2065-2076,doi: 10.6038/cjg20130627.

Zhang Y,Zhang H Z,Zhang G Q.2011.A stable TTI reverse time migration and its implementation.Geophysics,76(3): WA3-WA11.

Zhou H B,Bloor R,Zhang G Q.2006.An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.∥ 76th Annual International Meeting,SEG,Expanded Abstracts,194-198.

Zhou H Z,Zhang G Q,Zhang Y.2009.Removing S-wave noise in TTI reverse time migration.∥ 79th Annual International Meeting,SEG,Expanded Abstracts,2849-2853.

附中文參考文獻

杜啟振,秦童.2009.橫向各向同性介質(zhì)彈性波多分量疊前逆時偏移.地球物理學(xué)報,52(3): 801-807.

郝重濤,姚陳,王迅.2006.任意空間取向TI介質(zhì)中速度隨方位變化特征.地球物理學(xué)進展,21(2): 524-530.

黃翼堅,朱光明,劉池洋.2011.一個近似的VTI介質(zhì)聲波方程.地球物理學(xué)報,54(8): 2117-2123,doi: 10.3969/j.issn.0001-5733.2011.08.019.

吳國忱,梁鍇,印興耀.2010.TTI介質(zhì)彈性波相速度與偏振特征分析.地球物理學(xué)報,53(8): 1914-1923,doi: 10.3969/j.issn.0001-5733.2010.08.017.

張巖,吳國忱.2013.TTI介質(zhì)qP波逆時偏移中偽橫波噪聲壓制方法.地球物理學(xué)報,56(6): 2065-2076,doi: 10.6038/cjg20130627.

(本文編輯汪海英)

基金項目國家自然科學(xué)基金(41230317、41274112)資助.

作者簡介王璐琛，男，1989年生，在讀博士，主要從事地震波成像，微地震的反演和震源機制研究.E-mail：wangluchen@mail.iggcas.ac.cn *通訊作者常旭，研究員，研究方向為地震波傳播、地震成像、地球物理監(jiān)測.E-mail：changxu@mail.iggcas.ac.cn

doi:10.6038/cjg20160325 中圖分類號P631

收稿日期2014-12-03，2016-01-15收修定稿

Forward modeling of pseudo P waves in TTI medium using staggered grid

WANG Lu-Chen1,2，CHANG Xu1*，WANG Yi-Bo1

1KeyLaboratoryofShaleGasandGeoengineeringInstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcadamyofSciences,Beijing100049,China

AbstractA tilted transversely isotropic (TTI) medium is a good approximation for anisotropic problems.We discuss the property of a 3D tilted transversely isotropic (TTI) medium and use pseudo P waves to simulate the P component of the elastic wave field in the weak elastic anisotropy approximation.We compare forward features of 3 second-order difference equations of pseudo P waves based on Hooke′s law,elastic wave projection and dispersion equation,respectively.Now numerical calculation of these equations uses regular grid finite difference method,which is of low efficiency,low accuracy and is unstable in TTI simulation.In order to improve the calculating accuracy,this paper builds the staggered grid finite difference format for each forwarding simulation.Comparing three different numerical expressions of pseudo P wave equation in a 3D TTI medium,we find that the method based on dispersion equation is more convenient to simulate elastic P wave propagation and has lowest numerical noise.By analyzing the frequency dispersion characteristics based on different angles of isotropic symmetry axis,we introduce appropriate SV wave velocity in calculation to keep calculation stable.An example of a 2D simple synthetic model shows that the staggered grid forwarding method allow us to obtain a smooth and stable pseudo P wave field.With this strategy,2D BP TTI reverse time migration can also achieve basically stable migration results.

KeywordsPseudo P wave; Tilted transversely isotropic; Staggered grid

王璐琛，常旭，王一博.2016.TTI介質(zhì)的交錯網(wǎng)格偽P波正演方法.地球物理學(xué)報，59(3):1046-1058,doi:10.6038/cjg20160325.

Wang L C,Chang X,Wang Y B.2016.Forward modeling of pseudo P waves in TTI medium using staggered grid.Chinese J.Geophys.(in Chinese),59(3):1046-1058,doi:10.6038/cjg20160325.