黏聲波方程波場模擬的通用格式自適應系數頻域有限差分方法

徐文豪,巴晶*,J.M.Carcione,朱鶴松

1 河海大學地球科學與工程學院,南京 211100 2 National Institute of Oceanography and Applied Geophysics-OGS,Borgo Grotta Gigante 42/c,Trieste,Sgonico 34010,Italy

0 引言

由于地下介質的黏彈性,地震波在地下傳播時存在著衰減效應.黏聲波方程可以有效表征地震波傳播的衰減特征,因此基于黏聲波方程的波場數值模擬在地震資料處理、反演研究中具有重要應用,如基于黏聲波方程的全波形反演(Malinowski et al.,2011;李海山等,2018;Chen et al.,2021)和逆時偏移(Yang and Zhu,2018;周彤等,2018;Qu et al.,2020;谷丙洛等,2023)等.黏聲波模擬的主要方法包括時域有限差分方法(趙建國等,2014;Guo et al.,2016)、頻域有限差分(FDFD)方法(張立彬和王華忠,2010;Aghamiry et al.,2022)、有限元法(Zyserman and Gauzellino,2004;Groby and Tsogka,2006)、偽譜法(Zhu et al.,2014;吳玉等,2017)、近似解析離散法(宋國杰等,2019)、傅里葉有限差分法(張金海等,2007,2008)等.在這些方法中,頻域有限差分方法具有可靈活處理衰減效應、高效模擬多炮或窄帶地震資料、避免時間頻散、并行效率高等優點,因而在黏聲波模擬中具有較多的應用案例(Arntsen et al.,1998;李添才等,2014;Keating and Innanen,2019).用FDFD方法模擬黏聲波方程的主要挑戰在于其需要求解大型的稀疏線性方程組,而大型稀疏線性方程組的求解往往需要花費較多的計算量和內存.克服這一挑戰的一個常用方法是發展所需一個波長內的網格點數(number of points per wavelength,NPPW)更少的FDFD方法,從而減少波場模擬所需的網格點數,進而減少相應大型稀疏線性方程組的階數.

作為使用FDFD方法求解黏聲波方程的特例以及發展黏聲波FDFD方法的基礎,關于聲波方程的FDFD方法,前人已經進行了較多的研究.針對二維聲波方程,Pratt和Worthington(1990)提出了經典的五點FDFD方法,在相速度誤差不超過1%的精度要求下,該方法所需的NPPW為13.Liu(2014)提出了一個優化的五點FDFD方法,其精度明顯優于傳統的五點方法.Jo等(1996)提出一個適用于正方形網格的優化的九點FDFD方法,當x方向和z方向的空間采樣間隔相等的時候,該方法能將所需的NPPW減少到4.通過引入平均導數方法(average-derivative method,ADM),Chen(2012)將Jo等(1996)的FDFD方法推廣到長方形網格情況下.劉璐等(2013)提出一個優化的十五點FDFD方法,將所需的NPPW減少到2.97,同時相應的稀疏系數矩陣具有可接受的帶寬.Shin和Sohn(1998)提出一個優化二十五點FDFD方法,其所需的NPPW為2.5.基于ADM方法,張衡等(2014)將Shin和Sohn(1998)的方法推廣到了長方形網格情況.Fan等(2017)提出一個優化通用二十五點FDFD方法,該方法將所需的NPPW減少到2.13,而且很多常用的二維FDFD格式都可以被視作該方法的特殊情況.Xu和Gao(2018)提出一個自適應系數的九點FDFD方法,其在正方形網格和非正方形網格情況下所需的NPPW分別為2.12和2.44,同時所需的計算量和內存需求與優化九點FDFD方法相當.針對三維聲波方程,在相速度誤差不超過1%的精度要求下,傳統的二階七點FDFD方法所需的NPPW為13(Chen,2014).Operto等(2007)發展了一個適用于立方體網格的優化二十七點方法,將所需的NPPW減少到了4.Chen(2014)發展了一個優化ADM二十七點方法,該方法對一般的長方體網格所需的NPPW為4.Xu等(2021a)發展了自適應系數的三維二十七點FDFD方法,并給出了一種更容易使用和推廣的自適應系數FDFD通用格式.自適應系數的二十七點FDFD方法在立方體網格和非立方體網格下所需的NPPW分別為2.2和2.7.

為求解一般的黏聲波方程,傳統的二階五點FDFD方法所需NPPW較多(Amini and Javaherian,2011).高階FDFD方法可以有效減少所需的NPPW(Xu et al.,2021b),但其對應的大型稀疏線性方程組含有較大的帶寬和較多的非零元素,因此所需的計算量和內存較大(Takekawa and Mikada,2016).Arntsen等(1998)將Holberg(1987)采用的十字形優化有限差分系數應用于黏聲波方程的FDFD求解,并針對均勻模型應用解析解驗證了其有效性,十字形優化FDFD相對高階FDFD可以減少所需空間差分的長度,但相對緊湊格式的FDFD方法所需計算量和內存仍然較大,同時十字形優化FDFD在低波數區域的精度相對高階FDFD會有所減少(Zhang and Yao,2013).通過將聲波方程FDFD中的實波數替換為黏聲波方程的復波數,Amini和Javaherian(2011)將Jo等(1996)關于聲波方程的優化九點FDFD直接應用于黏聲波方程,并給出了在一般黏聲波模型中的測試,然而Jo等(1996)的FDFD只適用于正方形網格,同時Amini和Javaherian(2011)并沒有通過黏聲波方程的解析解或參考解驗證相應FDFD方法在一般黏聲波模型中的有效性.Chen(2012,2014)在提出聲波方程ADM優化FDFD方法的同時,給出了在黏聲波方程中的相應推廣,但并未給出一般黏聲波模型的數值算例.在將自適應系數FDFD應用到聲波全波形反演的同時,Aghamiry等(2022)給出了自適應系數FDFD方法在黏聲波方程的推廣,但同樣未給出一般黏聲波模型的數值算例,同時Aghamiry等(2022)所采用的基于旋轉坐標系的自適應系數FDFD方法只適用于空間采樣間隔相等的情形(Chen,2012).

為了減少不同空間采樣間隔之比下的黏聲波波場模擬所需的NPPW,本文針對黏聲波波場模擬發展了一種通用格式自適應系數FDFD方法.同時,為驗證自適應系數FDFD對一般黏聲波模型的有效性,本文針對三個典型的黏聲波模型,分別使用解析解和高階FDFD參考解驗證所提方法的有效性.在本方法中,所需的自適應系數FDFD格式可通過在傳統的二階FDFD格式的基礎上加入相關的系數隨NPPW和空間采樣間隔之比變化的校正項得到,其中校正項按網格點與中心點的距離進行分類選取.本文首先給出了基于黏聲波方程的通用格式自適應系數FDFD方法,然后分別針對均勻黏聲波模型、層狀黏聲波模型、Marmousi黏聲波模型,將所提FDFD與二階五點FDFD、ADM優化九點FDFD的波場模擬結果進行對比,并分別通過解析解和基于高階FDFD的參考解驗證所提方法的有效性.

1 通用格式自適應系數FDFD方法

1.1 自適應系數FDFD的通用格式

頻率域二維無源黏聲波方程可表示為

(1)

其中u表示頻率域波場,ω表示角頻率,vc表示復速度.關于復速度vc,本文采用如下常用的Kolsky-Futterman模型進行描述(Kolsky,1956; Futterman,1962):

(2)

其中vr表示在參考角頻率ωr下的參考波傳播速度,Q表示品質因子,sgn表示符號函數,i表示虛數單位.據Toverud和Ursin (2005)以及Aghamiry等(2020),本文將參考角頻率對應的頻率值設為50 Hz.需指出的是,公式(2)對應的一維平面波解的形式是u0e-i(ω/vc)x,如果對應的一維平面波解的形式是u0ei(ω/vc)x,則式(2)虛部的符號將改變.

(3a)

r=Δx/Δz,

(3b)

(4)

其中um,n=u(mΔx,nΔz),(vc)m,n=vc(mΔx,nΔz),ηi(i=1,2,3)表示通用格式中的自適應FDFD系數.

1.2 自適應FDFD系數的確定

據Operto等(2007)、Chen (2012,2014)、Aghamiry等(2022),黏聲波FDFD的優化或自適應FDFD系數可通過先假設品質因子趨于無窮大,再極小化相應平面波解的代入誤差得到.參考已有的關于黏聲波方程的FDFD方法(Chen,2012; Aghamiry et al.,2022),本文假設黏聲波FDFD的誤差主要來源于數值頻散,因此在確定自適應FDFD系數的過程中,只考慮品質因子趨于無窮大的情況(即Q→∞).此時黏聲波方程的平面波解與聲波方程的平面波解具有相同形式,可表示為

u(x,z)=u0e-ikr[sin(θ)x+cos(θ)z],

(5)

(6a)

(6b)

(6c)

(6d)

-1],

(6e)

其中η=[η1,η2,η3].極小化平面波解的代入誤差,所需的自適應FDFD系數可通過求解如下優化問題得到:

(7)

式(7)中的優化問題可通過Polak-Ribière-Polyak共軛梯度法(Polyak,1969)快速求解,

(8a)

ηk+1=ηk+αkdk,

(8c)

1.3 FDFD系數查找表

表1 基于聲波數值頻散關系的黏聲波方程自適應系數九點FDFD方法的部分系數查找表(r=3)

需指出的是,對于通用格式自適應系數九點FDFD方法,相應的查找表只是一個3×500的實數雙精度矩陣,因此幾乎不會增加波場模擬所需的內存需求.此外,對網格尺寸為nx×nz的波場模擬,構造自適應系數FDFD系數矩陣相對構造常系數FDFD系數矩陣只需額外增加3nxnz次線性插值,其計算量相對計算FDFD系數矩陣的各分量以及求解FDFD線性方程組而言很小.

1.4 數值頻散分析

使用FDFD方法數值求解波動方程時容易導致數值頻散,即FDFD方法對應的相速度與模型真實速度有偏差.這里給出通用格式自適應系數九點FDFD方法的數值頻散分析結果.令vph表示自適應FDFD方法對應的數值相速度.則自適應系數FDFD方法的平面波解(和波場方程的平面波解不同)可表示為

(9)

為得到通用格式自適應系數九點FDFD的數值相速度vph和模型真實速度vr的偏差,將自適應系數FDFD的平面波解代入自適應九點FDFD格式,即式(4),得到

(10a)

(10b)

(10c)

-1],

(10e)

圖1和圖2分別給出空間采樣間隔比率r=1和r=3時自適應系數九點FDFD和二階五點FDFD、優化ADM九點FDFD的數值頻散曲線對比.圖1和圖2結果表明,相對于二階五點FDFD和ADM優化九點FDFD,通用格式自適應系數九點FDFD可以更為有效的壓制正方形網格和非正方形網格下的數值頻散.以數值頻散誤差不超過1%為標準,相對于二階五點FDFD和ADM優化九點FDFD,新方法可將黏聲波方程一個波長內所需的網格點數從4減少到2.5.

圖1 r=1時不同FDFD方法的數值頻散曲線對比

圖2 r=3時不同FDFD方法的數值頻散曲線對比

1.5 吸收邊界條件

本文采用完美匹配層(perfectly matched layer,PML)吸收邊界條件(Liu and Tao,1997; Ren and Liu,2013;高靜懷等,2018)壓制人工邊界反射.帶有PML吸收邊界和震源項的二維黏聲波方程可表示為

-f(ω)δ(x-xs)δ(z-zs),

(11)

(12)

其中(xm,zn)=(mΔx,nΔz),δm,ms和δn,ns表示克羅內克函數,(ms,ns)表示震源的網格點下標.

2 數值算例

為驗證所提方法的有效性,針對三個典型的黏聲波模型,將新方法的模擬結果分別與二階五點FDFD以及ADM優化九點FDFD的模擬結果進行對比.同時,對均勻模型采用解析解作為參考解,對非均勻模型采用階數足夠高的高階FDFD作為參考解,高階FDFD采用的具體階數參考Zhang和Yao(2013)關于使用不同階數的有限差分近似空間偏導的波數域誤差的研究.對于所有算例,FDFD生成的大型稀疏線性方程組均通過PARDISO直接求解器(Petra et al.,2014)求解,計算平臺均為超算平臺的30個節點(AMD EPYC 7452 @2.35GHz),每個節點有64核和256 GB的內存,不同頻率的波場通過消息傳遞接口(message passing interface,MPI)多節點并行計算.

2.1 均勻模型算例

首先考慮網格點數為201×201的均勻模型.模型在50 Hz參考頻率下的速度和品質因子分別為2000 m·s-1和50.模型的采樣間隔為Δx=Δz=16.7 m,各側PML吸收邊界層數均為20,震源選為主頻20 Hz的雷克子波,震源位置在模型中心.檢波器布置在地表所有網格點處.波場模擬的頻率范圍是0.5～60 Hz且頻率采樣間隔為0.5 Hz.從頻域地震記錄轉換為時域地震記錄的時間范圍是0到2 s且時間采樣間隔為1 ms.圖3給出本算例使用的部分自適應FDFD系數(r=1),表明所使用的自適應FDFD系數是較光滑的.圖4給出了0.9 s時二階五點FDFD、ADM優化九點FDFD、自適應系數九點FDFD的波場快照對比.圖5給出了各檢波器接收的地震記錄對比.圖6給出了在(x,z)=(1670 m,0 m)的地表中心處上述三種FDFD接收的地震記錄與解析解(Carcione et al.,1988)的對比.圖4、圖5、圖6表明對這個均勻黏聲波模型,相對于二階五點FDFD和ADM優化九點FDFD,本文所提的自適應系數九點FDFD可以更有效地減少黏聲波模擬誤差.同時,圖4、圖5、圖6也給出了相應速度模型的聲波模擬結果,圖4和圖5體現了黏聲波傳播相對聲波傳播主頻降低的特征,圖6體現了黏聲波傳播相對聲波傳播振幅衰減的特征.此外上述三種FDFD在超算30個節點上的計算時間分別為1 s、1 s、1 s,所需內存分別為0.22 GB、0.25 GB、0.25 GB,所需的計算時間和內存相當.

圖3 均勻黏聲波模型使用的r=1時的自適應FDFD系數

圖4 均勻黏聲波模型(Q=50)在0.9 s時不同方法得到的波場快照對比

圖5 均勻黏聲波模型(Q=50)不同方法得到的接收地震記錄對比

圖6 均勻黏聲波模型(Q=50)中聲波方程解析解、黏聲波方程二階五點FDFD、黏聲波方程ADM優化九點FDFD、黏聲波方程通用格式自適應系數九點FDFD波場模擬結果與解析解單道對比

為驗證通用格式自適應系數FDFD的誤差主要來源于數值頻散,我們在上述均勻模型的基礎上,進一步考慮品質因子Q值的變化對模擬效果的影響.圖7給出了Q值分別取250、50、10的三種情形下,在(x,z)=(1670 m,0 m)地表中心處通用格式自適應系數FDFD接收的地震記錄與解析解的對比.圖7表明在不同的Q值情況下,通用格式自適應系數FDFD均能與解析解較好符合.

圖7 對取不同Q值的均勻黏聲波模型,通用格式自適應系數九點FDFD波場模擬結果與解析解單道對比

2.2 層狀模型算例

其次考慮一個網格點數為401×401的兩層模型.在50 Hz參考頻率下,模型在第一層和第二層的速度分別為2000 m·s-1和4000 m·s-1,品質因子分別為50和200.模型的采樣間隔為Δx=16.7 m、Δz=8.35 m,沿x和z方向的單側PML吸收邊界層數分別為20和40,震源選為主頻20 Hz的雷克子波,震源位置為(xs,zs)=(3340 m,83.5 m).檢波器布置在地表的所有網格點處.波場模擬的頻率范圍是0.25～60 Hz且頻率采樣間隔為0.25 Hz.從頻域地震記錄轉換為時域地震記錄的時間范圍是0到4 s且時間采樣間隔為1 ms.本算例采用沿x方向72階、沿z方向12階的高階FDFD作為參考解,該高階FDFD為傳統二階五點FDFD的高階版本(Xu et al.,2021b).圖8給出了本算例使用的部分自適應FDFD系數(r=2).圖9給出了1.5 s時上述四種FDFD的波場快照對比.圖10給出了各檢波器接收的地震記錄的對比.圖11給出了在(x,z)=(1670 m,0 m)處上述四種FDFD接收的地震記錄的對比.圖9、圖10、圖11表明,對這個層狀黏聲波模型,相對于二階五點FDFD和ADM優化九點FDFD,自適應系數九點FDFD可更為有效的減少黏聲波模擬的誤差.同時,圖9、圖10、圖11也給出了相應速度模型的聲波模擬結果,圖9和圖10體現了黏聲波傳播相對聲波傳播主頻降低的特征,圖11體現了黏聲波傳播相對聲波傳播振幅衰減的特征.此外上述四種FDFD在超算30個節點上的計算時間分別為3 s、4 s、4 s、216 s,所需內存分別為0.78 GB、0.91 GB、0.91 GB、16.5 GB.通用格式自適應系數九點FDFD的計算時間和所需內存與二階五點FDFD、ADM優化九點FDFD相當,均遠小于參考解所需的計算時間和內存.

圖8 兩層黏聲波模型使用的r=2時的自適應FDFD系數

圖9 兩層黏聲波模型在1.5 s時不同方法得到的波場快照對比

圖10 兩層黏聲波模型不同方法得到的接收地震記錄對比

圖11 兩層黏聲波模型中聲波方程高階FDFD(沿x方向72階、沿z方向12階)、黏聲波方程二階五點FDFD、黏聲波方程ADM優化九點FDFD、黏聲波方程通用格式自適應系數九點FDFD波場模擬結果與參考解(黏聲波方程沿x方向72階、沿z方向12階的高階FDFD)的單道對比

2.3 Marmousi模型算例

本算例考慮一個網格點數為737×750的Marmousi模型.在50 Hz參考頻率下,模型速度如圖12a所示,模型的品質因子(圖12b)是基于李氏經驗公式通過速度模型換算得到,即Q≈3.516×v2.2×10-6(李慶忠,1993).模型的采樣間隔為Δx=12.5 m、Δz=4 m,沿x和z方向的單側PML吸收邊界層數分別為20和40,震源選為主頻20 Hz的雷克子波,震源位置為(xs,zs)=(4600 m,40 m).檢波器布置在地表所有網格點處.波場模擬頻率范圍是0.125～60 Hz且頻率采樣間隔為0.125 Hz.從頻域地震記錄轉換為時域地震記錄的時間范圍是0～8 s且時間采樣間隔為1 ms.本算例采用沿x方向72階、沿z方向8階的高階FDFD作為參考解,該高階FDFD為傳統二階五點FDFD的高階版本(Xu et al.,2021b).圖13給出了本算例使用的部分自適應FDFD系數(r=3.125).圖14給出了各檢波器接收的地震記錄的對比.圖15給出了在(x,z)=(3750 m,0 m)處上述四種FDFD接受的地震記錄的對比.圖14和圖15的結果表明,對該Marmousi黏聲波模型,相對于二階五點FDFD和ADM優化九點FDFD,自適應系數九點FDFD可以更為有效的減少模擬誤差.同時,圖14和圖15也給出了相應速度模型的聲波模擬結果,圖14體現了黏聲波傳播相對聲波傳播主頻降低的特征,圖15體現了黏聲波傳播相對聲波傳播振幅衰減的特征.此外上述四種FDFD在超算30個節點的計算時間分別為20 s、23 s、24 s、974 s,所需內存分別為2.4 GB、3.0 GB、3.0 GB、47.5 GB,通用格式自適應系數九點FDFD的計算時間和所需內存與另兩種方法相當,遠小于參考解所需時間和內存.

圖12 Marmousi黏聲波模型

圖13 Marmousi黏聲波模型使用的r=3.125時的自適應FDFD系數

圖14 Marmousi黏聲波模型不同方法得到的接收地震記錄對比

圖15 Marmousi黏聲波模型中聲波方程高階FDFD(沿x方向72階、沿z方向8階)、黏聲波方程二階五點FDFD、黏聲波方程ADM優化九點FDFD、黏聲波方程通用格式自適應系數九點FDFD波場模擬結果與參考解(黏聲波方程沿x方向72階、沿z方向8階的高階FDFD)的單道對比

3 結論

本文針對黏聲波方程給出了一種適用于不同空間采樣間隔之比的通用格式自適應系數FDFD方法,通過假設黏聲波FDFD的誤差主要來源于數值頻散,所需的自適應FDFD系數根據聲波方程的數值頻散關系和查找表快速給出,并用于黏聲波方程FDFD的數值頻散分析和數值模擬中.針對三個典型的黏聲波模型,分別采用解析解和基于高階FDFD的參考解驗證了所提出方法的有效性.數值頻散分析表明,在空間采樣間隔相等或不等的情況下,以相速度誤差不超過1%為標準,通用格式自適應系數FDFD方法所需的一個波長內的采樣點數均小于2.5.數值實驗表明,在不同的空間采樣間隔之比的較粗網格情況下,相對于二階五點FDFD和ADM優化九點FDFD,新方法均可在相似的計算量和內存需求下,有效提高黏聲波數值模擬的精度.