王璐琛，常旭，王一博

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

?

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

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

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

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

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

1引言

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

由于彈性波的傳播本身已經(jīng)比較復(fù)雜，介質(zhì)的各向異性又會(huì)增加非零的Hooke彈性參數(shù)從而加重這種復(fù)雜的狀況，會(huì)大幅度增加波場(chǎng)正演的計(jì)算時(shí)間以及計(jì)算量.同時(shí)考慮到縱波勘探應(yīng)用廣泛，技術(shù)也更加成熟，是比較常用的研究方法.因此，只研究類似P波的波場(chǎng)傳播(稱為偽P波)，無疑是在各向異性研究中一種比較折中的方式(Ekkehart,2010; Yoon et al.,2010).針對(duì)偽P波的波場(chǎng)傳播，可以通過限定不同的約束條件(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ù)值實(shí)現(xiàn)上基本采用最直觀的離散方式，用規(guī)則網(wǎng)格差分將二階方程進(jìn)行數(shù)值離散.然而規(guī)則網(wǎng)格由于差分精度有限，一方面不利于在有限的存儲(chǔ)與計(jì)算時(shí)間下提高計(jì)算精度；另一方面在各向異性參數(shù)變化較為劇烈的區(qū)域，由于差分網(wǎng)格點(diǎn)間距過大，容易造成計(jì)算的不穩(wěn)定性.同時(shí)有些規(guī)則網(wǎng)格算法會(huì)將各向異性參數(shù)的微商引入導(dǎo)數(shù)項(xiàng)的計(jì)算中，這樣容易因各向異性參數(shù)微商計(jì)算的精度不夠而造成數(shù)值計(jì)算發(fā)散.

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

(5)

2.2弱各向異性下的P波運(yùn)動(dòng)學(xué)色散方程

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

首先研究VTI介質(zhì)的情況.在各向異性的情況下，波的相速度能夠描述波在介質(zhì)中傳播的特點(diǎn)，因此可以從相速度的角度建立波場(chǎng)的運(yùn)動(dòng)學(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方程進(jìn)行變量代換(Zhang et al.,2011)，設(shè)

(8)

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

(9)

(10)

(11)

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

(12)

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

(13)

(14)

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

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

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

(15)

(16)

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

(17)

(18)

(19)

(20)

(21)

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

(22)

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

(23)

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

(24)

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

(25)

和

(26)

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

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

3.1偽P波的正演

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4　三維空間中對(duì)稱軸不同轉(zhuǎn)角下的垂向波場(chǎng)(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ā)的不同時(shí)刻地震波的波場(chǎng)快照如圖6所示.其震源位置為橫向2000 m，深度1000 m處.圖6中從左至右的波場(chǎng)時(shí)刻依次是0.2、0.3、0.4、0.6 s.

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

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

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

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

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

圖7　單震源炮記錄

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

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

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

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

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

4結(jié)論

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

致謝感謝國家自然科學(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.

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

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

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

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

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

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

(本文編輯汪海英)

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

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

doi:10.6038/cjg20160325 中圖分類號(hào)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ì)的交錯(cuò)網(wǎng)格偽P波正演方法.地球物理學(xué)報(bào)，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.