基于LSQR算法的二維聲波方程頻率域正演模擬與數值實現

張鑫磊， 陳建宇

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

0 引言

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

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

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

1 基本原理

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

17點差分格式與傳統的最優9點差分格式類似，將各項同性介質的頻率域波動方程表示為0°和45°坐標系下的加權，用最優化方法計算得到加權系數[16]。

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

(1)

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

(2)

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

(3)

其中:Δ=Δx=Δz為網格間距，對四階9點差分格式在45°坐標系下進行擴展。則構造出的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 頻散特征

加權系數可以通過頻散方程采用最優化方法得到[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，為求出系數最優解，令k*=1/G，并使式(6)的值接近“1”，定義相速度的誤差為：

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

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

(7)

則目標函數為：

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

(8)

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

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

1.3 完全匹配層

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

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

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

(9)

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

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

(10)

則在x=0處的反射系數為：

(11)

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

(12)

(13)

(14)

其中：

(15)

(16)

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

(17)

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

(18)

結合式(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點網格厚度；(b)10點網格厚度

1.4 LSQR算法求解波場

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

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

(20)

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

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

(21)

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

(22)

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

(23)

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

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

Am×mxm×1=bm×1

(24)

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

1)初始化。

(25)

2)迭代循環。

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

圖6 均勻層狀模型Fig.6 Uniform layer model

(26)

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

3)參數修改。

(27)

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

2 模型試算

2.1 PML厚度取值

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

2.2 均勻層狀模型

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

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

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

圖10 Overthrust模型Fig.910 Overthrust model

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

2.3 Overthrust模型

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

表1 均勻層狀模型計算效率對比

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

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

3 結論

本次正演模擬采用優化17點差分算子格式，運用稀疏矩陣LSQR算法，并使用PML完全匹配層吸收邊界反射。優化17點格式在大網格間距的條件下表現出很小的頻散特性，可以節省大量的系統資源。在簡單層狀模型和Overthrust模型上的數值實驗也證明了LSQR在快速和準確的求解頻率域正演波場的同時只消耗很小的內存，高于直接解法的計算效率，因此在全波形反演領域將具有很好的應用前景。

表2 Overthrust模型計算效率對比

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