一階速度—壓力常分數階黏滯聲波方程及其數值模擬

陳漢明 汪燚林 周 輝*

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

0 引言

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

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

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

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

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

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

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

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

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

(1)

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

(2)

且有

(3)

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

(4)

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

(5)

式中

(6)

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

(7)

(8)

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

(9)

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

(10)

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

(11)

上式為頻率域的相乘關系，返回到時間域對應的卷積關系為

(12)

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

(13)

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

(14)

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

(15)

其中αx、κx、dx為

(16)

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

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

(17)

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

(18)

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

(19)

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

(20a)

(20b)

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

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

2 數值算例

2.1 均勻模型

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

2.2 穩定性分析

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

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

2.3 非均勻介質模型

圖3 BP鹽丘速度模型

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

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

圖5 單炮模擬記錄

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

3 結論

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

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

附錄A

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

(A-1)

式中

(A-2)

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

(A-3)

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

(A-4)

(A-5)

(A-6)

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