陳漢明 汪燚林 周 輝*

(①中國石油大學(北京)地球物理學院，北京 102249； ②油氣資源與探測國家重點實驗室，北京102249； ③CNPC物探重點實驗室，北京102249； ④同濟大學海洋與地球科學學院，上海200092)

0 引言

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

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

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

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

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

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

1 常分數(shù)階拉普拉斯算子黏滯聲波方程

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

Chen等[21]推導的位移形式的常分數(shù)階黏滯聲波方程為

(1)

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

(2)

且有

(3)

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

(4)

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

(5)

式中

(6)

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

(7)

(8)

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

(9)

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

(10)

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

(11)

上式為頻率域的相乘關(guān)系，返回到時間域?qū)木矸e關(guān)系為

(12)

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

(13)

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

(14)

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

(15)

其中αx、κx、dx為

(16)

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

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

(17)

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

(18)

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

(19)

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

(20a)

(20b)

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

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

2 數(shù)值算例

2.1 均勻模型

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

2.2 穩(wěn)定性分析

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

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

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

圖3 BP鹽丘速度模型

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

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

圖5 單炮模擬記錄

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

3 結(jié)論

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

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

附錄A

由Zhu等[17]推導的分數(shù)階拉普拉斯算子黏滯聲波方程為

(A-1)

式中

(A-2)

當介質(zhì)均勻時，分數(shù)階拉普拉斯算子可采用快速傅里葉變換(FFT)進行計算，即

(A-3)

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

(A-4)

(A-5)

(A-6)

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