徐世剛 包乾宗 任志明

(長安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，陜西西安 710054)

0 引言

地下介質(zhì)普遍具有非彈性和各向異性特征，對(duì)地震波的傳播速度、振幅和相位等產(chǎn)生影響，因此研究介質(zhì)的各向異性和黏滯性對(duì)于地震波場(chǎng)數(shù)值模擬、偏移成像和波形反演具有重要意義。

TI(Transversely Isotropic，橫向各向同性)介質(zhì)是一種常用的地球介質(zhì)近似模型，按照對(duì)稱軸的方向可以劃分為VTI(垂直對(duì)稱軸)、HTI(水平對(duì)稱軸)和TTI(傾斜對(duì)稱軸)三種類型。由于各向異性彈性波方程涉及波場(chǎng)變量多、計(jì)算效率低、多波信息復(fù)雜等，因此目前主要采用偽聲波方程進(jìn)行相關(guān)計(jì)算，該類方程是直接將TI介質(zhì)耦合頻散關(guān)系中的橫波速度置為0簡化而來[1]?；谠撍枷耄瑢?dǎo)出了多種改進(jìn)版本的各向異性偽聲波方程[2-6]。但是當(dāng)介質(zhì)參數(shù)不滿足橢圓假設(shè)時(shí)，大多數(shù)偽聲波方程模擬容易出現(xiàn)噪聲干擾及數(shù)值不穩(wěn)定。發(fā)展各向異性純聲波方程能夠較好地解決偽聲波方程存在的問題，其核心思想是從原始的TI介質(zhì)P-SV波頻散關(guān)系中解耦出純P波頻散關(guān)系，獲得對(duì)應(yīng)的純聲波方程，采用不同數(shù)值算法求解，以獲得純聲波響應(yīng)[7-10]。相比傳統(tǒng)偽聲波方程，純聲波方程能夠徹底消除偽橫波干擾及數(shù)值不穩(wěn)定，近年來在各向異性波場(chǎng)模擬及成像中受到較多關(guān)注[11-15]。目前關(guān)于各向異性純聲波方程的研究主要集中于模擬精度和計(jì)算效率的相互平衡，以及如何進(jìn)一步的推廣應(yīng)用。

在考慮地下介質(zhì)各向異性特征的同時(shí)，地層的黏滯性容易引起地震波的振幅衰減和相位畸變。在地震頻帶內(nèi)，廣義標(biāo)準(zhǔn)線性體和常Q模型被廣泛用于描述介質(zhì)的衰減特性[16-19]。多位學(xué)者利用差分法求解了基于標(biāo)準(zhǔn)線性體推導(dǎo)的黏聲或黏彈方程[20-23]。杜啟振等[24]采用有限元法求解黏彈各向異性波動(dòng)方程； Zhu等[25]采用差分法求解了基于標(biāo)準(zhǔn)線性體的黏聲和黏彈波動(dòng)方程，并與常Q理論計(jì)算的參考解進(jìn)行對(duì)比分析。相比于標(biāo)準(zhǔn)線性體模型，常Q模型對(duì)介質(zhì)衰減描述更準(zhǔn)確，但是相應(yīng)的理論推導(dǎo)和數(shù)值計(jì)算更復(fù)雜[25]。Carcione[26]采用偽譜法求解了均勻介質(zhì)中常Q黏聲波方程； Zhu等[27]利用偽譜法模擬了振幅和衰減項(xiàng)相互獨(dú)立的常Q黏聲波方程。為了同時(shí)兼顧常Q黏滯聲波方程的計(jì)算精度和效率，不同學(xué)者提出了相應(yīng)的解決策略[28-29]。

基于標(biāo)準(zhǔn)線性體理論，本文首先介紹了各向同性介質(zhì)中能表征振幅衰減和相位畸變的三維黏滯聲波方程； 其次，將泊松算子和有限差分聯(lián)合用于求解三維TTI介質(zhì)純聲波方程； 最后，結(jié)合各向同性黏滯聲波方程和TTI介質(zhì)純聲波方程，推導(dǎo)了適用于各向異性介質(zhì)、同時(shí)考慮振幅衰減和相位畸變的三維TTI黏滯介質(zhì)純聲波方程。數(shù)值算例表明本文方法能夠同時(shí)考慮地震波的振幅衰減和相位畸變，模擬精度高，適用于復(fù)雜模型的P波模擬。

1 理論方法

1.1 各向同性介質(zhì)三維黏滯聲波方程

地球介質(zhì)對(duì)地震波的衰減作用不僅體現(xiàn)在對(duì)振幅的衰減，同樣會(huì)引起相位的畸變[16-23]。在各向同性介質(zhì)中，基于標(biāo)準(zhǔn)線性體理論的黏滯聲波方程，Wang等[20]推導(dǎo)了一種同時(shí)包含振幅衰減和相位畸變的簡化黏滯聲波方程，其三維形式為

(1)

式中：P為縱波波場(chǎng)；v為縱波速度； ?2為三維拉普拉斯算子；Q為品質(zhì)因子；τ為與應(yīng)力和應(yīng)變松弛時(shí)間τσ、τε相關(guān)的參數(shù)

(2)

其中ω0為參考角頻率。

將式(1)分解為以下三個(gè)方程[20]

(3)

(4)

(5)

式(3)為傳統(tǒng)的聲波方程，式(4)為僅包含相位畸變項(xiàng)的黏滯聲波方程，式(5)是僅包含振幅衰減項(xiàng)的黏滯聲波方程。

(6)

式中：F和F-1分別表示傅里葉正、反變換；k為波數(shù)矢量；h為空間采樣間隔； Δt為時(shí)間采樣間隔；M為有限差分算子長度的一半；cI為有限差分系數(shù)。式(3)～式(5)可以采用類似的遞推格式求解。

采用三維各向同性層狀模型對(duì)上述方程進(jìn)行測(cè)試。模型尺寸為2km×2km×2km，上層速度為2400m/s，品質(zhì)因子為20，厚度為1km； 下層速度為3000m/s，品質(zhì)因子為40。震源是主頻為26Hz的Ricker子波，位于(1.0km，1.0km，0.9km)。圖1為聲波(式(3))、僅含相位畸變黏滯聲波(式(4))、僅含振幅衰減黏滯聲波(式(5))和同時(shí)考慮相位畸變和振幅衰減黏滯聲波方程(式(1))模擬的波場(chǎng)切片，圖2為四個(gè)方程模擬的單道波形對(duì)比。由圖1和圖2可以看出：與聲波方程的計(jì)算結(jié)果相比，僅含相位畸變黏滯聲波方程僅能模擬相移而無明顯振幅差異； 僅含振幅損失黏滯聲波方程僅能模擬振幅衰減而無明顯相移； 同時(shí)考慮相位畸變和振幅損失的黏滯聲波方程能夠同時(shí)模擬振幅衰減和相位畸變，更精確地描述了黏滯介質(zhì)的吸收衰減效應(yīng)。

圖1 三維各向同性層狀模型四種方程模擬的切片對(duì)比(a)xz平面，y=0.5km； (b)xy平面，z=0.5km； (c)yz平面，x=1km。由左往右依次為式(3)、式(4)、式(5)、式(1)

圖2 四個(gè)方程模擬的(x=1.0km，y=0.5km)處波形對(duì)比

1.2 基于泊松算子和有限差分相結(jié)合的三維TTI介質(zhì)純聲波方程

迄今為止，TI介質(zhì)中常用的聲波方程主要是基于偽聲波近似推導(dǎo)的，該類方程在各向異性參數(shù)不滿足近似條件時(shí)模擬結(jié)果會(huì)出現(xiàn)偽橫波干擾，甚至模擬過程不穩(wěn)定[1]，尤其是TTI介質(zhì)。多種改進(jìn)版本的TI偽聲波方程僅能在一定程度上壓制干擾，提高穩(wěn)定性，但不能從根本上解決問題[2-6]。眾多研究表明，發(fā)展各向異性純聲波方程能夠有效消除偽橫波污染和不穩(wěn)定現(xiàn)象[7-15]。相比于顯式表達(dá)的偽聲波方程，純聲波方程通常包含復(fù)雜的偏微分算子，常規(guī)差分方法很難直接求解?？紤]到計(jì)算精度與效率，本文應(yīng)用泊松算子和優(yōu)化有限差分相結(jié)合的方法求解三維各向異性純聲波方程[10]。在TI介質(zhì)中，為了便于計(jì)算，通常采用一階泰勒展開對(duì)解耦后的純聲波頻散關(guān)系進(jìn)行近似，三維TTI介質(zhì)近似純聲波頻散關(guān)系可以簡化為[7-12]

(7)

(8)

其中：Kx、Ky和Kz為正交波數(shù)；θ為對(duì)稱軸傾角；φ為對(duì)稱軸方位角。

將式(7)變換到時(shí)間—空間域，有

(9)

式中算子Rx、Ry和Rz的形式為

(10)

(11)

(12)

引入輔助波場(chǎng)H，式(9)可分解為兩個(gè)方程[8]

(13)

(14)

式(13)為隱式方程，其運(yùn)算過程需要求解線性方程組。為了提高效率，本文結(jié)合泊松算子實(shí)現(xiàn)[10]。

將式(13)和式(14)整理為

?2H=P

(15)

(16)

式(15)和式(16)組成新的TTI介質(zhì)三維純聲波方程，后者可采用優(yōu)化有限差分求解。前者為泊松方程，基于前人的研究，可以通過泊松算子求解，其實(shí)現(xiàn)流程可分為以下步驟[10]。

(1)基于五點(diǎn)差分法離散式(15)

6Hi,j,k-Hi+1,j,k-Hi-1,j,k-Hi,j+1,k-Hi,j-1,k-

Hi,j,k+1-Hi,j,k-1=-h2Pi,j,k

(17)

式中：i、j、k分別為x、y、z方向網(wǎng)格點(diǎn)序號(hào)。

(18)

式中：Nx、Ny和Nz為三個(gè)方向的網(wǎng)格點(diǎn)數(shù)；l、m、n為波數(shù)域三個(gè)方向的網(wǎng)格點(diǎn)序號(hào)。

(19)

(20)

通過以上四步，可實(shí)現(xiàn)泊松方程式(15)的求解。采用優(yōu)化有限差分可實(shí)現(xiàn)顯式方程式(16)的求解。上述流程能夠保證純聲波方程式(15)和式(16)的高效精確求解。

采用三維各向異性層狀模型對(duì)橫波速度為0的偽聲波方程[3]、限制性橫波速度的偽聲波方程[4]和本文的純聲波方程進(jìn)行測(cè)試。模型上層介質(zhì)參數(shù)為： 厚度為1.2km，vP0=2400m/s，ε=0.10，δ=0.05，θ=0°，φ=0°； 下層介質(zhì)參數(shù)為：厚度為0.8km，vP0=3000m/s，ε=0.10，δ=-0.12，θ=45°，φ=60°。震源是主頻為26Hz的Ricker子波，位于(1.0km，1.0km，0.9km)。圖3為不同方程計(jì)算的波場(chǎng)切片，可知，基于橫波速度為0的偽聲波方程和限制性橫波速度的偽聲波方程不可避免地會(huì)出現(xiàn)偽橫波干擾，而純聲波方程能夠獲得高精度P波響應(yīng)。

圖3 三維TTI層狀模型中三種方程模擬的切片對(duì)比(a)基于橫波速度為0的偽聲波方程； (b)基于限制性橫波速度的偽聲波方程； (c)本文純聲波方程從左到右依次為：xz平面，y=1.0km； xy平面，z=0.9km； yz平面，x=1.0km

1.3 簡化的三維TTI黏滯介質(zhì)純聲波方程

現(xiàn)有黏滯聲波方程主要是針對(duì)各向同性介質(zhì)設(shè)計(jì)的。對(duì)于各向異性介質(zhì)，尤其是各向異性純聲波方程，相應(yīng)的黏滯聲波方程研究較少。本文將三維各向同性黏滯聲波方程與TTI純聲波方程進(jìn)行結(jié)合，推導(dǎo)了適用于三維TTI介質(zhì)的黏滯聲波方程。

純聲波方程式(16)可寫為

(21)

式中算子S可表述為

S(P)=(1+2ε)(Rx+Ry)P+RzP-

2(ε-δ)(RxRz+RyRz)H

(22)

式(1)中的算子?2和式(21)中的算子S均是由波場(chǎng)相關(guān)的不同偏導(dǎo)數(shù)組合而成，因此將各向同性黏滯聲波方程式(1)與TTI介質(zhì)純聲波方程結(jié)合，即用算子S替換?2，可以得到三維TTI介質(zhì)中簡化的黏滯純聲波方程

(23)

類似地，可將式(23)分解為三個(gè)方程

(24)

(25)

(26)

式(24)為三維TTI純聲波方程，式(25)為僅包含相位畸變的TTI黏滯純聲波方程，式(26)是僅包含振幅衰減的TTI黏滯純聲波方程。利用式(23)即可實(shí)現(xiàn)三維TTI黏滯純聲波方程數(shù)值模擬。相應(yīng)的離散方案與前文所述類似。

應(yīng)用三維TTI層狀黏滯模型對(duì)比式(23)～式(26)數(shù)值模擬結(jié)果(圖4、圖5)。模型上層介質(zhì)參數(shù)為：厚度為1.2km，vP0=2400m/s，ε=0.10，δ=0.05，θ=0°，φ=0°，Q=20； 下層介質(zhì)參數(shù)為：vP0=3000m/s，ε=0.10，δ=-0.12，θ=45°，φ=60°，Q=40。震源是主頻為26Hz的Ricker子波，位于(1.0km，1.0km，0.9km)。由圖4、圖5可以看出，與TTI純聲波方程(式(24))的模擬結(jié)果相比，僅含相位畸變TTI黏滯純聲波方程(式(25))僅能模擬時(shí)移，僅含振幅損失TTI黏滯純聲波方程(式(26))僅能模擬振幅衰減，同時(shí)考慮相位畸變和振幅損失TTI黏滯純聲波方程(式(23))能夠同時(shí)模擬振幅衰減和相位畸變，更精確地描述了三維TTI介質(zhì)中黏滯純聲波的吸收衰減屬性。

將三維TTI層狀黏滯模型上下層品質(zhì)因子由(20，40)變?yōu)?40，60)和(60，80)，分析品質(zhì)因子對(duì)波場(chǎng)模擬結(jié)果的影響(圖6、圖7)。由圖6、圖7可見，相比于純聲波方程的計(jì)算結(jié)果(圖4、圖5左)，隨著品質(zhì)因子的增大，吸收衰減效應(yīng)導(dǎo)致地震波的振幅減弱、相位變化的程度降低(圖5右、圖6)。

圖4 三維TTI黏滯層狀模型中四個(gè)方程模擬的(x=1.0km，y=1.0km)處波形曲線對(duì)比

圖5 三維TTI黏滯層狀模型不同方程模擬的波場(chǎng)快照對(duì)比(a)xz平面，y=1.0km； (b)xy平面，z=0.9km； (c)yz平面，x=1.0km。由左往右依次為式(24)、式(25)、式(26)、式(23)

圖6 三維TTI黏滯層狀不同Q值模型的簡化黏滯純聲波方程模擬的波場(chǎng)快照(a)xz平面，y=1.0km； (b)xy平面，z=0.9km； (c)yz平面，x=1.0km。左：上下層介質(zhì)Q為(40，60)； 右：上下層介質(zhì)Q為(60，80)

2 數(shù)值模擬

2.1 楔狀體模型

首先采用修改的三維各向異性楔狀體模型進(jìn)行測(cè)試。模型尺寸為1.5m×1.5m×1.5m，圖8為該模型xz平面示意圖，沿y方向直接延拓可得其三維形式。模型A區(qū)的參數(shù)為：vP0=1500m/s，ε=0，δ=0，θ=-25°，φ=45°，Q=34.2； B區(qū)的參數(shù)為：vP0=1500m/s，ε=0.30，δ=0.02，θ=-25°，φ=45°，Q=34.2； C區(qū)的參數(shù)為：vP0=4500m/s，ε=0.10，δ=0.01，θ=60°，φ=45°，Q=382.9； D區(qū)的參數(shù)為：

圖7 三維TTI層狀不同Q值模型簡化黏滯純聲波方程模擬的在(x=1.0km，y=1.0km)處波形對(duì)比(a)上層介質(zhì)； (b)下層介質(zhì)

圖8 各向異性楔狀體模型xz平面切片

vP0=3500m/s，ε=0.05，δ=0.01，θ=10°，φ=45°，Q=220.2。震源是主頻為20Hz的Ricker子波，位于(0.75km，0.75km，0.15km)。圖9對(duì)比了偽聲波和純聲波方程的模擬結(jié)果，可以看出，相比于偽聲波方程的模擬結(jié)果，純聲波方程能夠較好地消除偽橫波干擾。圖10為黏滯純聲波方程模擬波場(chǎng)快照，與純聲波方程的模擬結(jié)果(圖9c)相比，能量明顯減弱、分辨率降低，有明顯的振幅和相位差異(圖11)。

圖9 三維TTI楔狀體模型不同方程模擬的波場(chǎng)快照(a)基于橫波速度為0的偽聲波方程； (b)基于限制性橫波速度的偽聲波方程； (c)本文純聲波方程。從左到右依次為：xz平面，y=0.75km； xy平面，z=0.40km； yz平面，x=0.75km

圖10 三維TTI楔狀體模型黏滯純聲波模擬的波場(chǎng)快照(a)xz平面，y=0.75km； (b)xy平面，z=0.40km； (c)yz平面，x=0.75km

圖11 三維TTI楔狀體模型的純聲波和黏滯純聲波方程模擬結(jié)果的(x=0.75km，y=0.75km)處波形對(duì)比

2.2 復(fù)雜模型

采用修改的三維Marmousi TTI模型驗(yàn)證本文方法對(duì)復(fù)雜模型的適應(yīng)性。模型尺寸為1.5km×1.5km×2.1km，圖12為該改進(jìn)模型的xz平面示意圖，沿y方向直接延拓可得三維形式，對(duì)稱軸方位角固定為30°。震源選用主頻為20Hz的Ricker子波，置于(0.74km，0.74km，0.18km)。圖13和圖14分別為純聲波和黏滯純聲波方程模擬的波場(chǎng)快照，圖15為兩種方程模擬波形曲線對(duì)比，可以看出，與純聲波方程的模擬結(jié)果相比，基于黏滯純聲波方程的模擬結(jié)果能量顯著減弱，出現(xiàn)明顯的振幅降低和相位畸變。該算例表明，本文的TTI黏滯純聲波方程可以適應(yīng)復(fù)雜各向異性黏滯介質(zhì)模型的P波波場(chǎng)模擬。

圖12 修改的Marmousi TTI模型的xz切片(a)vP0； (b)ε； (c)δ； (d)θ； (e)Q

圖13 三維Marmousi TTI模型純聲波方程模擬的波場(chǎng)快照左：xz平面，y=0.74km； 中：xy平面，z=0.18km； 右：yz平面，x=0.74km

圖14 三維Marmousi TTI模型黏滯純聲波方程模擬的波場(chǎng)快照左：xz平面，y=0.74km； 中：xy平面，z=0.18km； 右：yz平面，x=0.74km

圖15 三維Marmousi TTI模型純聲波和黏滯純聲波方程模擬的(x=0.74km,y=0.74km)處波形對(duì)比

3 結(jié)論

本文首先介紹了一種三維各向同性介質(zhì)中同時(shí)考慮振幅衰減和相位畸變的黏滯聲波方程。針對(duì)各向異性介質(zhì)傳統(tǒng)偽聲波方程模擬容易產(chǎn)生偽橫波干擾及出現(xiàn)數(shù)值不穩(wěn)定，結(jié)合泊松算子和有限差分法給出了一種適用于三維TTI介質(zhì)的純聲波求解算法，有效消除了偽聲波方程模擬結(jié)果存在的干擾及數(shù)值不穩(wěn)定?；诟飨蛲责暡ǚ匠毯蚑TI介質(zhì)純聲波方程，推導(dǎo)了一種針對(duì)三維TTI介質(zhì)的黏滯純聲波方程，該方程能夠同時(shí)考慮黏滯各向異性介質(zhì)中的地震波振幅衰減和相位畸變屬性。

應(yīng)用三維TTI層狀模型模型驗(yàn)證了本文方法的模擬精度； TTI楔狀體和Marmousi模型驗(yàn)證了方法的有效性和適用性。