基于LSQR算法的二維聲波方程頻率域正演模擬與數(shù)值實現(xiàn)

張鑫磊， 陳建宇

(中國石油大學(北京) 地球科學學院，北京 102200)

0 引言

近幾十年人們對地震波數(shù)值模擬進行了許多研究，取得了不少成就。大多數(shù)的波場模擬都是在時間域進行的，時間域計算方法是按時間片遞推計算，由前一個時間片的波場值推出下一個時間片的波場值。因此前一時刻的計算誤差會累加到下一時刻，如果計算的時間切片數(shù)目較多，計算誤差將會累計，導致計算結(jié)果信噪比下降、精度太低而無法為接下來的波形反演提供研究基礎(chǔ)。近三十年來，學者做了大量工作，頻率域正演模擬最早由Lysmer等[1]提出；Shin等[2]進一步發(fā)展了這種方法，將頻率域有限元法應用于地震波形反演; Jo等[3]提出了最優(yōu)化9點差分格式，極大地減小了數(shù)值頻散效應;Stekl等[4]在Jo的基礎(chǔ)上，將最優(yōu)化9點差分格式推廣到變密度聲波方程和黏彈性介質(zhì)的彈性波方程; Min等[5]提出了一種頻率域標量波動方程25點差分格式，通過計算最優(yōu)化系數(shù)來減小數(shù)值頻散，減少了空間采樣點數(shù);Hustedt[6]研究了交錯網(wǎng)格上導數(shù)算子的四階精度有限差分格式，給求解方程組帶來更大帶寬。在國內(nèi)，近些年來吳國忱教授[7-8]做了大量研究，將25點優(yōu)化差分算子理論應用于VIT介質(zhì)和TTI介質(zhì)的正演模擬，取得了較好的模擬效果;王華忠教授[9]也在這方面做了大量研究，利用最優(yōu)化9點差分算子對標量波方程做了數(shù)值模擬。

頻率域正演模擬涉及到大型稀疏矩陣的求解問題，求解的方法主要分為直接法(高斯消元法和LU分解法等)與迭代法(LSCG最小二乘共軛梯度法和LSQR最小二乘正交分解法等)。直接法通常占用內(nèi)存大，資源消耗嚴重，因此筆者采用LSQR(Least Square QR-factorization)最小二乘正交分解法，LSQR算法是Sanders和Paige[10-11]提出的，它是利用Lanczos法求解最小二乘問題的一種方法。劉伊克等[12-13]將其運用到地震層析成像處理中；隨后楊文采[14]利用該算法完善了地球物理反演理論；劉勁松等[15]又對該算法進行了并行改進，使其計算效率大大增加。LSQR具有占用儲存空間少、計算速度快等優(yōu)點，可以利用矩陣稀疏性簡化運算，適合用來求解大型稀疏矩陣。與其他迭代方法相比，LSQR具有更快的收斂速度及更精確的結(jié)果。

筆者采用的差分格式為優(yōu)化17點差分算子格式[15]，適合高精度正演模擬。該差分格式是在四階精度9點差分格式的基礎(chǔ)上，結(jié)合0°系和45°坐標系構(gòu)造而來，具有計算精度高，時間短，矩陣帶寬小等優(yōu)點。在17點差分格式基礎(chǔ)上同時采用稀疏矩陣LSQR算法計算波場，并使用PML完全匹配層吸收邊界反射，最終達到提高正演模擬效率的目的。

1 基本原理

1.1 頻率域17點有限差分格式

17點差分格式與傳統(tǒng)的最優(yōu)9點差分格式類似，將各項同性介質(zhì)的頻率域波動方程表示為0°和45°坐標系下的加權(quán)，用最優(yōu)化方法計算得到加權(quán)系數(shù)[16]。

首先構(gòu)造四階精度的9點差分格式，時間域二維聲波方程對時間t做傅里葉變換可以得到頻率域聲波方程:

(1)

地震波場上任意位置拉普拉斯算子可以寫成加權(quán)和的形式：

(2)

對式(1)采用4階精度有限差分格式可得0°坐標系9點差分格式：

(3)

其中:Δ=Δx=Δz為網(wǎng)格間距，對四階9點差分格式在45°坐標系下進行擴展。則構(gòu)造出的17點差分格式可表示為：

圖1 17點差分格式頻散曲線Fig.1 Dispersion curves of 17-point scheme

F(ω)=β1Um-2,n-2+β2Um-2,n+β1Um-2,n+2+

β3Um-1,n-1+β4Um-1,n+β3Um-1,n+1+

β2Um,n-2+β4Um,n-1+γUm,n+β4Um,n+1+

β2Um,n+2+β3Um-1,n+1+β4Um+1,n+

β3Um+1,n+1+β1Um+2,n-2+β2Um+2,n+

β1Um+2,n+2

(4)

1.2 頻散特征

加權(quán)系數(shù)可以通過頻散方程采用最優(yōu)化方法得到[17-18]，將平面簡諧波表達式代入式(4)，可得到17點格式頻散方程：

4c-4d-4e)D)]/{(3aΔ2[30-16A+

B]+3Δ2(1-a)[15-16C+D])}

(5)

相速度的定義為：vρh=ω/k，代入式(5)可得到相速度的表達式：

2(1-b-4c-4d-4e)D)}/{(3aΔ2[30-

16A+B]+3Δ2(1-a)[15-16C+D])}

(6)

傳播角度的取值范圍為[0,π/4]，1/G的最大取值為0.4，為求出系數(shù)最優(yōu)解，令k*=1/G，并使式(6)的值接近“1”，定義相速度的誤差為：

ε(θ,k*,a,b,c,d,e)= [1-(vρh(θ,k*,a,

b,c,d,e)/v)]2

(7)

則目標函數(shù)為：

c,d,e)/v]2dk*dθ

(8)

依據(jù)式(8)可以解出最終的優(yōu)化系數(shù)：a=1.067 3、b=0.887 5、c=0.025 1、d=0.023 7、e=-0.020 4，圖2給出了 17點格式的頻散曲線，可以看出相速度誤差基本可控制在1%以內(nèi)。

圖2 完全匹配層邊界示意圖Fig.2 PML border diagram

1.3 完全匹配層

為了消除邊界反射的影響必須引入人工邊界條件，筆者采用頻率-空間域的完全匹配層PML邊界條件，在最外層網(wǎng)格區(qū)域的四個方向上構(gòu)建PML[19-20]。

設(shè)計算區(qū)域從x= -d開始，則x=0處為PML層邊界，PML層的厚度為0至 -d，設(shè)為L，對平面波方程加入衰減系數(shù)，則在邊界處入射波經(jīng)過了兩次PML層的衰減，可表示為：

U(kx,kz,ω)=U(kx,kz,ω)

(9)

在x=0處的反射波方程為：

Ur(kx,kz,ω)=RUr(kx,kz,ω)

(10)

則在x=0處的反射系數(shù)為：

(11)

我們可以通過設(shè)置計算區(qū)域衰減因子為零值，邊界外匹配層衰減因子為非零值而達到消除邊界反射的目的[21-22]。對式(11)兩邊做傅里葉反變換，并求x的二階導數(shù)，可得：

(12)

(13)

(14)

其中：

(15)

(16)

同理可得縱向匹配層關(guān)系式為式(17)。

(17)

角點區(qū)域與邊界區(qū)域的處理方式不同，如圖2所示，在四個角點區(qū)域?qū)Σ▓龅乃p是對x方向和z方向衰減的疊加。則PML條件下的聲波方程可寫為式(18)。

(18)

結(jié)合式(4)可得PML條件17點有限差分格式為式(19)

(19)

圖3 正演模擬的計算流程Fig.3 Calculation flow of forward modeling(a)直接法正演流程；(b)LSQR算法正演流程

圖4 PML效果圖Fig.4 PML effect diagram(a)5點網(wǎng)格厚度；(b)10點網(wǎng)格厚度

1.4 LSQR算法求解波場

將波場空間劃分為Nx×Nz的網(wǎng)格區(qū)域，可以令l=Nx×Nz，將式(17)改寫為：

F(ω)=S(v,ω)U(ω)

(20)

其中：S為復阻抗矩陣，大小為l×l；U為頻率域的離散波場值，是l×1的列向量；F為震源項，大小也是l×1。首先將復系數(shù)矩陣轉(zhuǎn)化為實系數(shù)分塊矩陣，則式(18)可改寫為：

(SRe+iSIm)(URe+iUIm)=FRe+iFIm

(21)

式中：SRe和SIm分別為阻抗矩陣的實部和虛部；URe和UIm分別為波場矩陣的實部和虛部，i為虛數(shù)單位。將(19)式展開可以得到一個方程組為式(22)。

(22)

寫成分塊矩陣為式(23)。

(23)

將復阻抗矩陣S按照式(21)排列，只存儲非零元，然后采用LSQR法求解方程，就可以快速高效的求解出各網(wǎng)格點的波場值。

將要求解的式(21)簡寫為式(24)。

Am×mxm×1=bm×1

(24)

設(shè)x0為式(21)要求解的波場x的初始值，β0為式(21)已知震源信息的二范數(shù)，則LSQR算法步驟總結(jié)如下：

1)初始化。

(25)

2)迭代循環(huán)。

圖5 PML效果圖(18點網(wǎng)格厚度)Fig.5 PML effect diagram (18 point mesh thickness)

圖6 均勻?qū)訝钅Ｐ虵ig.6 Uniform layer model

(26)

圖7 180 ms波場快照Fig.7 Snapshot of 180 ms wave field(a)頻率域LSQR法波場；(b)常規(guī)時間域波有限差分波場

3)參數(shù)修改。

(27)

其中:式(25)中的xi+1為最終輸出的波場虛實分解的列向量，將其虛實結(jié)合重構(gòu)即可求出最終波場值，圖3為正演模擬的計算流程，左為直接法正演流程，右為LSQR算法正演流程，主要區(qū)別在于求解頻率域波場的方法上，LSQR算法采用阻抗矩陣虛實分解和稀疏矩陣存儲方法，具有更快的計算速度和更少的資源占用。

2 模型試算

2.1 PML厚度取值

PML厚度不同取值效果見圖4、圖5，均為均質(zhì)模型180 ms波場快照，由圖4、圖5可以看出當取值達到18點網(wǎng)格厚度時已經(jīng)有較好的吸收效果，所以最終取值為18點網(wǎng)格厚度。

2.2 均勻?qū)訝钅Ｐ?/h3>

試算模型為100×100網(wǎng)格的均勻?qū)訝钅Ｐ?圖6)，上層速度為1 500 m/s，下層速度3 000 m/s，網(wǎng)格間距Δ=Δx=Δz=10 m，震源采用Rick子波，主頻30 Hz，位于(21，51)點，采樣間隔1 ms，記錄長度1 024 ms，檢波器接收位置位于x=20，每隔一個網(wǎng)格點放置一個，邊界采用18網(wǎng)格厚度的PML邊界，求解波場使用LSQR算法(使用不完全LU分解算法作計算效率對比)，同時使用二階時間精度四階空間精度常規(guī)時域有限差分法取同時刻波場做結(jié)果對比，參數(shù)一致[24-27]。圖7、圖8為傅里葉反變換后時間域的波場快照與時間域有限差分波場快照對比，圖9為炮集記錄。

圖8 180 ms波場快照局部Fig.8 Snapshot of 180 ms local wave field(a)頻率域LSQR法波場；(b)常規(guī)時間域波有限差分波場

圖9 共炮點道集記錄Fig.9 Common shot point gather

圖10 Overthrust模型Fig.910 Overthrust model

從圖7(a)與圖7(b)對比以及圖8(a)與圖8(b)對比圖中可以看到，頻率域LSQR差分波場和時間域差分波場差別很小，PML對邊界反射波的吸收效果很好，沒有數(shù)值頻散，共炮點道集也顯示清晰的直達波同相軸和雙曲線形態(tài)的反射波同相軸，表1為不完全LU分解、LSQR和常規(guī)時間域有限差分法的計算效率對比，可以看出LSQR占用內(nèi)存最小，計算時間最短，節(jié)省了大量系統(tǒng)資源。

2.3 Overthrust模型

模型采用較為復雜的Overthrust模型中的部分區(qū)域(圖10)，網(wǎng)格范圍為120×120，最大速度為6 000 m/s，最小速度為2 674 m/s，平均速度為4 292 m/s，網(wǎng)格間距取Δx=Δz=25 m，震源位置網(wǎng)格點為(27，27)，檢波器位置x=22，其他參數(shù)與上個模型相同，表2為計算效率對比，圖11、圖12為波場快照對比，圖13為共炮點道集記錄，在波場對比可以發(fā)現(xiàn)頻率域LSQR有限差分波場與常規(guī)時間域波場只有細微差異，產(chǎn)生差異主要原因是在大網(wǎng)格間距下，四階精度時間域差分波場會產(chǎn)生輕微數(shù)值頻散，影響其模擬精度，而優(yōu)化17點差分波場擁有更好的頻散特性，因此模擬精度更高[16]。炮集記錄可以看到清晰的直達波和雜亂的反射波曲線，旅行時也與計算模型相吻合，證明了正演模擬的正確性。

表1 均勻?qū)訝钅Ｐ陀嬎阈蕦Ρ?/p>

圖11 250 ms波場快照Fig.11 Snapshot of 250 ms wave field(a)頻率域LSQR法波場；(b)常規(guī)時間域波有限差分波場

圖12 250 ms波場快照局部Fig.12 Snapshot of 250 ms local wave field (a)頻率域LSQR法波場；(b)常規(guī)時間域波有限差分波場

3 結(jié)論

本次正演模擬采用優(yōu)化17點差分算子格式，運用稀疏矩陣LSQR算法，并使用PML完全匹配層吸收邊界反射。優(yōu)化17點格式在大網(wǎng)格間距的條件下表現(xiàn)出很小的頻散特性，可以節(jié)省大量的系統(tǒng)資源。在簡單層狀模型和Overthrust模型上的數(shù)值實驗也證明了LSQR在快速和準確的求解頻率域正演波場的同時只消耗很小的內(nèi)存，高于直接解法的計算效率，因此在全波形反演領(lǐng)域?qū)⒕哂泻芎玫膽们熬啊?/p>

表2 Overthrust模型計算效率對比

圖13 共炮點道集記錄Fig.13 Common shot point gather