矩形旋轉交錯網格PML邊界條件的研究與實現

張 亮，孫盡堯

（武漢大學 物理科學與技術學院，湖北 武漢430072）

矩形旋轉交錯網格PML邊界條件的研究與實現

張 亮，孫盡堯

（武漢大學 物理科學與技術學院，湖北 武漢430072）

現有的頻率域波全形反演多是默認基于正方形旋轉交錯網格，而基于矩形旋轉交錯網格的全波形反演有更大的靈活性，能更好地適應實際地下介質結構，有更好的內存使用率和計算效率。為了解決原始的基于正方形網格的PML（完全匹配層）吸收邊界形式不再適用的問題，文中將通過聲波方程完全匹配層吸收邊界的經典方法，構建并推導出PML吸收邊界條件在矩形旋轉交錯網格中的具體形式，通過Marmousi模型正演和反演實驗，得出本文所構建的PML邊界獲得良好的吸收效果，對反演結果有很好的改善。

頻率域全波形反演；矩形旋轉交錯網格；正演；完全匹配層

頻率域全波形反演是在時域波形反演的基礎上發展起來的，經過Marfurt[1]、Shin[2]、Worthington[3]、Pratt[4-5]等進一步地深入研究與完善，利用頻率域全波形反演的方法來重建地下介質結構，越來越受到人們的關注。與時域全波形反演相比，頻域全波形反演計算效率更高，能處理更大的地下介質模型，有更好的實用性。

交錯網格最早由Madariage[6]提出，Pratt改進了經典五點有限差分法求解波動方程，Jo[7]等人提出了九點旋轉交錯網格有限差分法來提高拉普拉斯域和質量加速項的精度[8-9]。而這些方法是基于正方形網格的，這樣造成的影響是在一些領域無法適應地下結構，同時可能增加不必要的內存和計算量，這對于大模型全波形反演是一個難點。而實際地下介質結構，橫向速度變化較小，縱向速度變化較大，基于矩形網格的全波形反演有更大的靈活性，能更好地適應實際地質結構，同時可減少網格點數目，有更好的內存使用率和計算效率。

在實際地震勘探中，地震波在地下是無限傳播的，而在我們波場模擬的過程中，由于計算機存儲的限制，在我們所需研究的區域外，加入了人為的截斷邊界。地震波會在人為邊界上產生反射，反射回的波會產生新的波場，對我們預期的波場產生干擾。因此，需要引入吸收邊界條件，消除人為截斷邊界的影響?，F在應用最廣泛的人工邊界條件是Berenger提出的完全匹配層法[10]，有著比較理想的吸收效果。其基本思想[11]是，在模擬區域周圍添加完全匹配層，在完全匹配層中加入吸收衰減因子，吸收人為截斷邊界差生的反射波。

在矩形旋轉交錯網格頻率域全波形反演中，隨著矩形網格寬長比的變化，原始的基于正方形網格的PML邊界形式已經不再適用。因此，文中從聲波方程完全匹配層吸收邊界的經典方法出發，構建并推導出PML邊界條件在矩形旋轉交錯網格中的形式。在原矩形旋轉交錯網格頻域全波形反演算法的基礎上，添加文中所構建的矩形旋轉交錯網格中的PML邊界條件，以Marmousi模型為目標模型，通過正演和反演實驗，表明文中所提方法的有效性。

1 方法原理

1.1 一般網格PML層聲波方程的離散形式

二維線性標量聲波方程形式如下[12]：

式中P（x，z，t）為模型中（x，z）點處的聲波波場，v（x，z）為模型中（x，z）點處的介質速度，S（x，z，t）為時域震源，?2為拉普拉斯算子。

為簡單起見，不考慮PML層中震源S的作用。按照經典PML邊界的理論[13]，將波場P分裂為兩項，即P=Px+Pz，通過引入中間變量、加入衰減因子，可構造如下PML吸收邊界條件：

式中Q為中間變量，γ為PML層衰減系數。再將將式（2）變換回頻域、差分離散、消去中間變量、合并，可推出一般網格PML層聲波方程的離散形式為：

式中Δx為橫向網格間距，Δz為縱向網格間距，ξ=1+iγ/ω。

1.2 矩形旋轉交錯網格中PML邊界的構建

矩形旋轉交錯網格的基本思路如圖1所示，對任意點和其8鄰域點進行有限差分，將微分形式近似轉換為差分形式，對聲波方程進行離散，然后求解離散后的聲波方程。如圖2所示，里面包含兩種坐標系，一種是直角坐標系，另外一種是旋轉后的坐標系。

圖1 矩形網格示意圖

圖2 旋轉坐標系示意圖

矩形旋轉交錯網格PML邊界的構建分為3個部分：一是旋轉坐標系，推導出聲波方程在新旋轉坐標系下的具體形式；二是參照經典PML邊界的構建方法，推導出新旋轉坐標系下的離散形式；三是將原正交坐標系和新旋轉坐標系下的PML邊界混合加權，融合生成矩形旋轉交錯網格下的PML邊界。

旋轉后坐標系x′z′下的坐標與原正交坐標系xz下的坐標關系為：

對（4）式兩邊求二階微分并整理可得：

由（1）（4）式可以得出旋轉后新坐標系x′z′下的聲波方程為：

與構建一般網格PML邊界同樣的過程，可推導出旋轉新坐標系x′z′下PML層聲波方程的離散形式為：

設原正交坐標系xz下的PML層所占比重為η，新旋轉坐標系x′z′下的PML層所占比重為1-η，構建權重因子如下：

將式（3）、（6）分別乘以η和1-η，即可將兩個坐標系下PML層聲波方程按權重融合即可得到矩形旋轉交錯網格中PML層的離散形式。至此，矩形旋轉交錯網格中PML邊界的離散形式構建完畢。

2 算法數值實驗

為了驗證文中提出的基于矩形旋轉交錯網格PML邊界的有效性，文中將以Marmousi模型進行基于矩形旋轉交錯網格的頻域全波形反演實驗，從正演和反演結果兩個方面進行說明。在算法實現過程中，軟件平臺使用Visual Studio 2010，運行環境為64位windows7操作系統，編程語言是C++，在處理線性方程組時調用MUMPS[14]庫來進行求解。

Marmousi模型橫向距離范圍為0～9 200米，橫坐標網格間距為12.5米，共737個網格；縱向深度0～3 000米，縱坐標網格間距12.5米，共240個網格。為將Marmousi按矩形網格劃分，文中對Marmousi模型的縱向距離進行線性插值，在不改變整個模型大小的情況下，將模型縱坐標網格間距減少到10米，縱向網格數目變為360。炮點共240個，分布于地表，炮點覆蓋范圍分布于模型的橫坐標從3 000米到8 975米處，炮間距25米。每個炮點的接收點均分布于該炮點左側，最大炮檢距（炮點和接收點之間的距離）為2 575米，最小炮檢距200米，道間距25米，每個炮點有96個接收點，即96道。

2.1 正演實驗

正演是利用初始速度模型得到模擬波場記錄，是反演的前提和基礎。文中初始速度模型選取勻速遞增模型。正演的采樣時間為3 s，采樣間隔為4 ms，共750個采樣點。震源為Ricker子波，主頻為8 Hz。

以Marmousi模型第一炮的數據為例，第一個炮點位于地表3 000米處，96個接收點分布在425～2 800米間，即接收道的第1個接收點位于425米處，第96個接收點位于2 800米，中間94個檢波點平均分布于此區域內。分別記錄無PML層和添加本文所構建的PML層（20層網格厚度）時，第一個炮點正演的波形記錄，其中橫坐標表示96個接收點的道號，縱坐標表示采樣時間。

圖3（a）顯示，當原原矩形旋轉交錯網格頻域全波形反演算法不添加PML吸收層時，在邊界處反射現象明顯。由于正演是反演的前提和基礎，正演結果的精度將在很大程度上影響反演結果，所以邊界處產生的反射必將對反演過程帶來很大的干擾。圖3（b）顯示，加入文中所構建的PML層后反射現象得到明顯抑制。由此可見，文中所構建的PML邊界取得了良好的吸收效果，這也將會對反演結果帶來積極影響。

2.2 反演實驗

反演是建立理論波場與模擬波場之間殘差的目標函數，通過使目標函數極小來更新速度模型，是一個對目標函數最優化的問題。文中使用的是預處理梯度算法[15]。

圖3 Marmousi模型第一個炮點正演模擬的波形記錄

圖4 Marmous模型的頻率域全波形反演

圖4 （c）和4（d）分別為不添加PML層和添加文中所構建PML層反演結果的反演結果，可以看出添加本文所構建的PML層后，反演結果整體上結構更清晰。圖4（e）給出了在水平距離左邊界2.50 km處，兩組結果的對比，可以看出添加本文構建PML層的反演結果，速度分布更接近理論模型。這是由于PML層很好地吸收了地震波在人工邊界處的反射，減少了反射波的干擾，實驗結果更精確，從而更接近理論模型，這與正演的實驗結果也相符合，由此可以驗證文中所提算法的有效性。

3 結束語

文中構建并推導出PML邊界在矩形旋轉交錯網格中的形式，通過正演和反演實驗，驗證其對人為邊界處的反射波有很好的吸收效果，對最終的反演結果也有很好的提升。另外，由于本文在選取PML層衰減因子時，沿用了正方形網格PML層衰減因子的形式。而隨著矩形網格寬長比的變化，構建新的衰減因子也是一個值得研究的方向。

[1]Marfurt K J.Accuracy of finite difference and finite element modeling of the scalar and elastic waveequations[J].Geophysics，1984，49（5）:533-549.

[2]Shin C S.Nonlinear elastic wave inversion by blocky parameterization [D].Tulsa，Oklahoma，United States:University of Tulsa，1988.

[3]Pratt R G，Worthington M H.Inverse theory applied to multi-source cross-hole tomography.Part I:Acoustic wave-equation method[J].Geophys Prosp，1990，38（3）:287-310.

[4]Pratt R G.Inverse theory applied to multi-source cross-hole tomography， PartII:Elastic waveequation method[J].Geophys Prosp，1990，38（3）: 287-310.

[5]Pratt R G.Frequency-domain elastic wave modeling by finite differences:A tool for crosshole seismic imaging[J].Geophysics，1990，55（5）:626-632.

[6]Madariaga R.Dynamics of expanding circular fault [J].Bulletin of the Seismological Society of America，1976，66（3）:639-666.

[7]Jo CH，Shin C，Suh JH.An optimal 9-point，finite-difference，frequency-space，2-D scalar wave extrapolator[J].Geophysics，1996（61）：529-537.

[8]Erik H S，Norbert G，Serge A S.Modeling the propagation of elastic waves using a moodified finitedifference grid[J].Wave Motion，2000（31）:7-792.

[9]Reeshidev B，Mrinal K S.Finite-difference modeling of S-wave splitting in anisotropic media[J]. Geophysical Prospecting，2008（56）:293-312.

[10]Berenger J A.perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics，1994（114）:185-200.

[11]Hastings F D，Schneider J B，Broschat S L. Application of the perfectly matched layer（PML）absorbing boundary condition to elastic wave propagation[J].Journal of the Acoustical Society of America，1996，100（100）:3061-3069.

[12]Bjorn Engquist.Computationalhigh frequency wave propagation[J].Acta Numerica，2003:181-266.

[13]Collino F，Tsogka C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics，2012，66（1）:294-307.

[14]Solver.MUltifrontal Massively ParallelSolver（MUMPS 4.10.0）Users guide[R].French:MUMPS，2011.

[15]Operto S，Virieux J，Dessa J X，et al.Crustal imaging from multifold ocean bottom seismometers data by frequency-domain fullwaveform tomography: Application to the eastern ankaiTrough[J].Journal of Geophysical Research，2006，111（B9）：171.

Research and implementation of the PML boundary conditions based on rotated rectangular-grid modeling method

ZHANG Liang，SUN Jin-yao
（School of Physics and Technology，Wuhan University，Wuhan 430072，China）

The present Frequency-domain full waveform inversion（FWI）is mostly default to a rotated square-grid，but the FWI based on rotated rectangular-grid is more flexible，and can better adapt to the actualstructure ofunderground media，and have greateradventage in memory utilization and computational efficiency.To overcome the problem that the original square-grid PML（Perfectly Latched Layer）absorbing boundary conditions no longer apply，the paper give the new form of PML absorbing boundary conditions based on rotated rectangular-grid through the classic building strategy.Using the forward modeling and inversion experiment with Marmousi model，we show that the proposed PML boundary in this article can obtain satisfactory absorbing effect，and improve the inversion result.

frequency-domain full waveform inversion；rotate rectangular-grid；forward modeling；PML

TN91

：A

：1674-6236（2017）02-0163-05

2016-04-12稿件編號：201604110

大型油氣田及煤層氣開發國家科技重大專項課題（2011ZX05003-003）

張 亮（1991—），男，湖北棗陽人，碩士。研究方向：二維地震波頻率域全波形反演。