一階速度—壓力常分?jǐn)?shù)階黏滯聲波方程及其數(shù)值模擬

陳漢明 汪燚林 周 輝*

(①中國石油大學(xué)(北京)地球物理學(xué)院，北京 102249； ②油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室，北京102249； ③CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京102249； ④同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院，上海200092)

0 引言

地震波在地下傳播會表現(xiàn)出依賴頻率的衰減，同時(shí)其相位也會因介質(zhì)的黏滯性而發(fā)生變化?；谒p介質(zhì)假設(shè)建立相應(yīng)的波動方程并研發(fā)相應(yīng)的波動方程偏移和反演算法，是目前眾多學(xué)者的研究內(nèi)容[1-3]。通常使用常Q模型描述地震波的衰減[4]，在頻率域常Q模型可利用復(fù)速度表達(dá)，但數(shù)值求解頻率域的波動方程涉及大量的矩陣和向量運(yùn)算，效率低。

相對而言，時(shí)間域的模擬方法更高效，但需要近似應(yīng)力—應(yīng)變之間的卷積關(guān)系，且難以準(zhǔn)確匹配常Q模型。孫成禹等[5]實(shí)現(xiàn)了利用廣義線性體構(gòu)建常Q模型的優(yōu)化方法，并證實(shí)至少需要三個(gè)線性體元在給定的頻率范圍內(nèi)擬合常Q曲線，才能保證大炮檢距地震波衰減和頻散的模擬精度。由此所發(fā)展的廣義線性體黏滯波動方程所包含的表征參數(shù)較多(每個(gè)線性體元包含兩個(gè)弛豫時(shí)間參數(shù))，方程個(gè)數(shù)也較多，將其作為正演模擬的工具，會增加后續(xù)地震偏移和反演的難度。杜啟振[6]將廣義線性體模型進(jìn)一步推廣至各向異性介質(zhì)，推導(dǎo)了更加復(fù)雜的黏彈性各向異性介質(zhì)波動方程，并采用偽譜法進(jìn)行數(shù)值模擬。蔡瑞乾等[7]提出了變線性體元數(shù)量的數(shù)值模擬方法，該方法在不同的區(qū)域采用不同個(gè)數(shù)的線性體元擬合常Q曲線，提高了廣義線性體黏滯波動方程數(shù)值模擬的計(jì)算效率。苑春方等[8]、李曉波等[9]基于Kelvin-Voigt模型和孫成禹等[10]基于Maxwell模型建立了非常Q黏滯波動方程并在時(shí)間域進(jìn)行地震波場模擬。

與傳統(tǒng)整數(shù)階導(dǎo)數(shù)的局部性質(zhì)不同，有學(xué)者利用分?jǐn)?shù)階導(dǎo)數(shù)的記憶性質(zhì)描述固體介質(zhì)的蠕變過程，由此發(fā)展了分?jǐn)?shù)階黏彈性波動方程，用于描述地震波隨頻率呈指數(shù)衰減的物理現(xiàn)象[11]。Carcione等[12]首先引入時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)建立符合常Q模型的黏滯聲波方程，能準(zhǔn)確描述常Q模型，但時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的全局性使數(shù)值模擬需要消耗大量的內(nèi)存存儲所有歷史時(shí)刻波場。孫成禹等[13]利用加權(quán)窗函數(shù)提高了時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的計(jì)算精度，同時(shí)利用自適應(yīng)的窗函數(shù)寬度提高了黏滯聲波和黏彈性波場數(shù)值模擬的計(jì)算效率。劉財(cái)?shù)萚14]基于VIT—雙相介質(zhì)的本構(gòu)關(guān)系和常Q模型推導(dǎo)了時(shí)間域分?jǐn)?shù)階黏彈性波動方程，用于模擬復(fù)雜介質(zhì)中的波場傳播。

Carcione[15]從常Q模型的頻散關(guān)系出發(fā)，引入分?jǐn)?shù)階拉普拉斯算子描述地震波的衰減和頻散。與時(shí)間分?jǐn)?shù)階黏滯聲波方程相比，分?jǐn)?shù)階拉普拉斯算子黏滯聲波的時(shí)間導(dǎo)數(shù)為二階，可利用傳統(tǒng)的有限差分法近似，同時(shí)分?jǐn)?shù)階拉普拉斯算子可使用快速傅里葉變換(FFT)計(jì)算，因此具有較高的計(jì)算效率。Zhang等[16]基于正則化方法推導(dǎo)了解耦的分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程，由于推導(dǎo)過程不是基于波動方程的頻散關(guān)系，其物理意義不明確。Zhu等[17]及吳玉等[18]受Treeby等[19]在醫(yī)學(xué)成像領(lǐng)域研究的啟發(fā)，基于常Q模型和平面波的頻散關(guān)系推導(dǎo)了新的解耦分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程，該方程形式簡單，模擬常Q模型的精度高，同時(shí)在形式上顯式地分離了控制振幅衰減和相位畸變的算子。這種解耦的特點(diǎn)有利于發(fā)展穩(wěn)定的Q補(bǔ)償逆時(shí)偏移技術(shù)[20]。

數(shù)值求解分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程的難點(diǎn)在于處理空間可變階數(shù)的拉普拉斯算子。由于階數(shù)隨空間變化，分?jǐn)?shù)階拉普拉斯算子的譜響應(yīng)也隨空間變化，對于空間每個(gè)不同的階數(shù)，需要單獨(dú)使用FFT計(jì)算分?jǐn)?shù)階拉普拉斯算子的譜響應(yīng)，計(jì)算效率低。Zhu等[17]直接取空間平均的階數(shù)代替空間可變階數(shù)，其模擬的結(jié)果存在較大的誤差。Chen等[21]提出兩種方法近似空間可變階數(shù)的拉普拉斯算子：其一為利用泰勒近似將空間可變階數(shù)的拉普拉斯算子轉(zhuǎn)化為若干個(gè)常分?jǐn)?shù)階拉普拉斯算子，由此發(fā)展了常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程，簡化了數(shù)值計(jì)算；另一種方法將空間可變分?jǐn)?shù)階拉普拉斯算子的譜響應(yīng)視為波數(shù)—空間的混合域矩陣，采用矩陣低秩分解方法近似。Chen等[22]證實(shí)，基于矩陣低秩分解的時(shí)間外推方法比傳統(tǒng)偽譜法有更高的時(shí)間近似精度和穩(wěn)定性。Li等[23]提出最小二乘法近似空間可變階數(shù)的拉普拉斯算子，并發(fā)展了相應(yīng)的衰減補(bǔ)償逆時(shí)偏移算法。最小二乘法與泰勒近似法類似，都將變階數(shù)轉(zhuǎn)化為固定階數(shù)，但由最小二乘法導(dǎo)出的黏滯波動方程的表征參數(shù)與速度和Q沒有顯式的表達(dá)關(guān)系，不利于相應(yīng)的全波形反演。Wang等[24]將常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程推廣到了黏彈介質(zhì)； Wang等[25]基于常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程發(fā)展了自適應(yīng)的衰減補(bǔ)償逆時(shí)偏移方法； Chen等[26]基于常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程發(fā)展了全波形反演方法； Chen等[27]發(fā)展了一種矩陣轉(zhuǎn)換法用于求解分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程，該方法是傳統(tǒng)有限差分法的推廣，與偽譜法相比，矩陣轉(zhuǎn)換法處理各種復(fù)雜邊界條件更靈活，但計(jì)算量更大。

本文基于位移形式的常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程，推導(dǎo)了一階速度—壓力形式常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程，并考慮了密度的空間變化；推導(dǎo)了相應(yīng)的完全匹配層(PML)邊界條件以壓制虛假反射；采用交錯(cuò)網(wǎng)格偽譜法進(jìn)行數(shù)值模擬。借助數(shù)值實(shí)驗(yàn)，分析Q對交錯(cuò)網(wǎng)格偽譜法穩(wěn)定性條件的影響，并驗(yàn)證新的一階速度—壓力黏滯波動方程模擬常Q模型具有較高的精度。

1 常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程

1.1 一階速度—壓力形式常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程的建立

Chen等[21]推導(dǎo)的位移形式的常分?jǐn)?shù)階黏滯聲波方程為

(1)

式中：u為位移；c0為定義在參考角頻率ω0處的地震波速度；Q為地震品質(zhì)因子； 算子D和L為

(2)

且有

(3)

為了考慮密度的空間可變以及方便加載PML吸收邊界條件，將式(1)等價(jià)變換為一階速度—壓力形式，為此引入中間變量

(4)

將中間變量代入式(1)，可得

(5)

式中

(6)

進(jìn)一步將密度項(xiàng)ρ加入到一階波動方程(式(5))中，由此得到變密度的常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程

(7)

(8)

式中： F和F-1分別表示傅里葉正、反變換；kx為x方向的波數(shù)；hx為x方向的網(wǎng)格間距。同樣，分?jǐn)?shù)階拉普拉斯算子也采用FFT計(jì)算，如

(9)

1.2 卷積型完全匹配層吸收邊界

(10)

式中：dx為PML區(qū)域的衰減函數(shù)；κx為改善以掠射形式入射到PML層的地震波吸收效果的系數(shù)；αx為改善對低頻成分吸收的系數(shù)。根據(jù)鏈?zhǔn)角髮?dǎo)法則，有

(11)

上式為頻率域的相乘關(guān)系，返回到時(shí)間域?qū)?yīng)的卷積關(guān)系為

(12)

根據(jù)式(10)中sx的表達(dá)式可知

(13)

式中：δ(t)是狄拉克函數(shù)；H(t)是單位階躍函數(shù)。將式(13)代入式(12)，可得

(14)

式中φx(t)為輔助變量

(15)

其中αx、κx、dx為

(16)

式中：w為PML厚度；x為PML區(qū)域計(jì)算點(diǎn)到PML內(nèi)邊界的距離；cmax為最大速度；R是理論反射系數(shù)，一般取R=10-5；κmax可取不同值，為簡單起見本文取為0；αmax=πfd，其中fd為震源的主頻。

式(15)中φx(t)的計(jì)算可采用迭代更新

(17)

式中Δt為波場外推的時(shí)間步長，其他參數(shù)定義為[29]

(18)

式(14)和式(17)表明CPML中的卷積項(xiàng)可以利用遞推公式計(jì)算，且只需要存儲前一個(gè)時(shí)刻的輔助變量值。

(19)

由交錯(cuò)網(wǎng)格的配置方式可知輔助變量φvx,x和φvz,z配置在整數(shù)時(shí)間網(wǎng)格節(jié)點(diǎn)上，φp,x和φp,z配置在半時(shí)間網(wǎng)格節(jié)點(diǎn)上，因此其迭代更新的公式為

(20a)

(20b)

值得注意的是，用于更新φp,x的變量bx和ax落在x方向的半網(wǎng)格節(jié)點(diǎn)上，用于更新φp,z的bz和az落在z方向的半網(wǎng)格節(jié)點(diǎn)上。

需要指出的是，式(19)中加載CPML吸收邊界的方法并不完全嚴(yán)格，因?yàn)榭刂祁l散的算子Dv和控制衰減的算子L中仍然包含空間導(dǎo)數(shù)項(xiàng)，理論上這些導(dǎo)數(shù)項(xiàng)也需要作相應(yīng)的復(fù)坐標(biāo)擴(kuò)展，但此項(xiàng)工作相當(dāng)復(fù)雜，本文不作深入的研究。數(shù)值實(shí)驗(yàn)將證實(shí)，本文所采用的近似CPML方法仍然能取得較好的吸收效果。

2 數(shù)值算例

2.1 均勻模型

首先模擬網(wǎng)格數(shù)為1100×1100、網(wǎng)格間距為10m的均勻介質(zhì)中的波場，并與格林函數(shù)計(jì)算的解析解[30]進(jìn)行對比。模擬的時(shí)間長度為2.5s，時(shí)間步長為1ms，采用主頻為20Hz的雷克子波作為震源并置于模型的中心，距離震源4km設(shè)置一個(gè)檢波點(diǎn)。定義參考頻率20Hz的地震波的速度為2000m/s，密度為2000kg/m3。圖1為接收點(diǎn)處數(shù)值模擬的地震道與解析解的對比，其中Q分別為5、10、20、50、80、100。由圖1可知，當(dāng)Q=5時(shí)，數(shù)值解的相位明顯滯后于解析解，說明本文所推導(dǎo)的常分?jǐn)?shù)階黏滯聲波方程在Q較小時(shí)存在較明顯的誤差； 隨著Q逐漸增大，數(shù)值解的相位與解析解更吻合，但在波峰和波谷處仍存在微小的振幅差異，這可能是由于傳播過程中數(shù)值頻散誤差累積造成的(偽譜法采用二階差分近似時(shí)間微分算子，時(shí)間近似精度只有二階)。圖1中數(shù)值解與解析解的對比證實(shí)了本文所推導(dǎo)的黏滯聲波方程能較好地匹配常Q模型，模擬精度高。

2.2 穩(wěn)定性分析

圖1 炮檢距為4km處數(shù)值解與解析解的對比

圖2 Q值不同時(shí)黏滯聲波方程交錯(cuò)網(wǎng)格偽譜法模擬所允許的最大網(wǎng)格比

2.3 非均勻介質(zhì)模型

圖3 BP鹽丘速度模型

圖4a和4b為模擬的黏滯聲波質(zhì)點(diǎn)振動速度波場切片，時(shí)刻分別為1.0和1.5s，圖4c和4d分別為相同時(shí)刻的完全彈性聲波波場切片。由圖可知，介質(zhì)的黏滯性引起地震波振幅顯著衰減。此外，由于使用了CPML吸收邊界條件，所模擬的波場中未見明顯的邊界反射，證實(shí)了本文加載CPML方法的可行性。圖5a為地表接收的完全彈性聲波共炮點(diǎn)記錄，圖5b為對應(yīng)的黏滯聲波記錄，可以看出，黏滯性導(dǎo)致1.5s以下反射波振幅明顯衰減。為了更清楚地展示介質(zhì)黏滯性對地震波振幅和相位的影響，抽取(2400m，0)處的地震道進(jìn)行對比，其中所使用的Q模型包括上述由經(jīng)驗(yàn)公式生成的非均勻Q模型(近地表最小Q值約為20)以及Q=200的均勻Q模型。圖6a為直達(dá)波波形對比，可見介質(zhì)的黏滯性造成地震波相位的明顯滯后，Q越小，相位滯后越明顯，振幅衰減也越強(qiáng)；圖6b為主要的反射波形對比，可觀察到相同的現(xiàn)象。

圖4 BP模型模擬的波場切片

圖5 單炮模擬記錄

圖6 BP模型(2400m，0)處不同Q值模擬結(jié)果單道波形對比

3 結(jié)論

本文將位移形式的常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程整理推導(dǎo)至一階速度—壓力形式，考慮了空間可變的密度，并仿照聲波方程討論了其卷積型完全匹配層吸收邊界條件的加載方法。雖然該完全匹配層吸收邊界方法并未完全遵循完全匹配層的理論推導(dǎo)，但數(shù)值實(shí)驗(yàn)仍證實(shí)了其較好的吸收效果。

均勻介質(zhì)模擬的數(shù)值解與解析解的對比證實(shí)，本文所推導(dǎo)的常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程匹配常Q模型的精度較高。數(shù)值模擬實(shí)驗(yàn)表明，當(dāng)Q值越小時(shí)，其穩(wěn)定性條件越苛刻，當(dāng)Q值大于40時(shí)，其穩(wěn)定性條件與模擬傳統(tǒng)聲波方程相同。BP鹽丘模型波場數(shù)值模擬結(jié)果證實(shí)了本文的黏滯聲波方程數(shù)值模擬方法對復(fù)雜介質(zhì)的適用性。

附錄A

由Zhu等[17]推導(dǎo)的分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程為

(A-1)

式中

(A-2)

當(dāng)介質(zhì)均勻時(shí)，分?jǐn)?shù)階拉普拉斯算子可采用快速傅里葉變換(FFT)進(jìn)行計(jì)算，即

(A-3)

將式(A-1)變換到波數(shù)域并經(jīng)過整理，可得

(A-4)

(A-5)

(A-6)

同樣的近似方法可用于式(A-1)右端的第二項(xiàng)，Chen等[21]通過數(shù)值分析證實(shí)當(dāng)ε=1/16時(shí)，式(A-6)具有很高的精度。經(jīng)過上述近似之后，空間可變階數(shù)的黏滯聲波方程(式(A-1))簡化為空間常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程(式(1))。