張春麗，張偉

1. 中國科學技術大學地球和空間科學學院，安徽 合肥 230026

2. 南方科技大學深圳市深遠海油氣勘探技術重點實驗室，廣東 深圳 518055

3. 南方科技大學地球與空間科學系，廣東 深圳 518055

1 背 景

地震波數(shù)值模擬是復雜介質逆時偏移成像和全波形反演的基礎。有限差分法（FDM）由于易于理解和實施，網(wǎng)格生成簡單，易于進行大規(guī)模并行計算等原因，是地震波成像和反演最主要的地震波正演方法。提高地震波有限差分模擬的計算效率，可以顯著節(jié)省地震波成像和反演工作的成本。有限差分模擬的計算資源需求與格點數(shù)量和時間步數(shù)量近似成正比，對確定的計算區(qū)域和時間窗長度，網(wǎng)格步長越小則網(wǎng)格數(shù)越多；時間步長越小則時間步數(shù)越多。為壓制數(shù)值頻散和耗散誤差，網(wǎng)格步長需要在所有的網(wǎng)格點上滿足所用格式的每波長網(wǎng)格點數(shù)要求，如果使用均勻網(wǎng)格離散模型，則網(wǎng)格步長由整體速度模型中的最低速度決定。另外，離散幾何形狀復雜尖銳的散射體，例如地質尖滅，也需要使用較小的網(wǎng)格步長。而在速度值高的區(qū)域和速度變化平緩的區(qū)域，全局細網(wǎng)格會造成空間過采樣（圖1（a））。此外，由于顯式時間積分格式穩(wěn)定條件限制，在高速區(qū)采用較小的空間步長時，需要同時采用較小的時間步長。綜上所述，對于速度差異大的復雜介質模型，采用均勻網(wǎng)格進行地震波模擬會導致空間和時間過采樣，浪費計算資源，增加計算時間。

克服均勻網(wǎng)格計算效率低的問題的一種方法是采用可變網(wǎng)格?？勺兙W(wǎng)格可按照網(wǎng)格步長是否連續(xù)變化分為兩類。文獻［1］最早提出了一種簡單的網(wǎng)格步長連續(xù)變化的方案［2-4］，如圖1（b）所示。由于和均勻網(wǎng)格相比，網(wǎng)格的拓撲結構沒有變化，所以在粗細網(wǎng)格的波場間不需要進行變量信息交換；然而，由于細網(wǎng)格區(qū)域會不可避免地從目標細化區(qū)延伸至計算區(qū)域邊界，所以這種方案只是部分減少了網(wǎng)格數(shù)，并且對于時間步長過小問題沒有顯著改善。另外一種可變網(wǎng)格是間斷網(wǎng)格，即網(wǎng)格步長在變化邊界突變整數(shù)倍，如圖1（c）所示。同級網(wǎng)格覆蓋區(qū)域呈現(xiàn)層狀或者塊狀［5-8］。與連續(xù)變化的可變網(wǎng)格相比，間斷網(wǎng)格更為靈活，減少了計算的網(wǎng)格點數(shù)。但在低速體形狀不規(guī)則時，受網(wǎng)格區(qū)域簡單幾何形狀的限制，間斷網(wǎng)格仍然會細化部分非目標區(qū)域造成浪費。

自適應網(wǎng)格（AMR，adaptive mesh refinement）由一系列邊界形狀不規(guī)則的多級網(wǎng)格塊組成，如圖1（d）所示，能夠靈活地覆蓋不規(guī)則的目標細化區(qū)域。從圖1中的不同類型網(wǎng)格的格點數(shù)目可以看出，自適應網(wǎng)格可以最大化減少不必要的網(wǎng)格數(shù)，是比間斷網(wǎng)格更為高效的網(wǎng)格劃分方案。文獻［9］最早提出自適應網(wǎng)格技術，同有限差分法結合，用于求解雙曲型偏微分方程；文獻［10］提出了一種格點聚類算法（point clustering）用于自適應網(wǎng)格的自動生成。在彈性動力學領域，許多研究者將自適應網(wǎng)格算法應用到有限體積法（finitevolume method）或者間斷伽遼金方法（discontinuous Galerkin method）中：文獻［11］提出使用高階間斷伽遼金和自適應網(wǎng)格模擬固液耦合介質中波場的傳播。文獻［12］在二維聲波方程有限差分模擬中實現(xiàn)了自適應網(wǎng)格，顯著提高了復雜速度模型中的聲波模擬效率。由于彈性波方程比聲波更為復雜，在相同網(wǎng)格點需要定義更多的變量，導致實現(xiàn)自適應網(wǎng)格難度更大，目前尚未見到基于自適應網(wǎng)格的彈性波有限差分模擬算法，本文將首次實現(xiàn)二維彈性波模擬的自適應網(wǎng)格有限差分算法。

圖1 不同類型的網(wǎng)格和格點數(shù)（灰度值表示波速，黃線表示不同網(wǎng)格步長區(qū)域的邊界。）Fig.1 Different grids and corresponding numbers of grid points(The gray values indicate the velocities of geological bodies,and the yellow lines represent the boundary between patches of different grid spacings.)

2 二維彈性波控制方程與同位網(wǎng)格有限差分算法

本文求解二維一階速度應力方程組形式的PSV波方程，該方程組的矢量形式為

其中U為速度-應力矢量：

其中vi代表i方向的速度分量，σij代表應力分量。震源項

其中ρ代表密度，fi表示體力，Mij，t表示地震矩張量關于時間的導數(shù)。系數(shù)A和B分別是

其中λ和μ是Lamé常量。

空間域采用中心差分格式對方程進行求解。以x方向為例，同位網(wǎng)格空間精度為2N階的中心差分格式為

其中下標表示空間索引，Lx表示x方向的空間差分，Δx表示空間步長。我們使用八階精度的差分格式［13］，N= 4，系數(shù)分別為a1= 4/5，a2= -1/5，a3= 4/105，a4= -1/280。

時間域采用如下的四階Runge-Kutta 時間積分格式，

其中Un表示時間步為n時波場各分量的值，L(U) =ALx(U) +BLz(U)，Δt表示時間步長，系數(shù)分別為α2=α3= 0.5，α4= 1，β1=β4= 1/6，β2=β3= 1/3.

為了解決中心差分格式的奇偶失聯(lián)導致的不穩(wěn)定問題，我們還進行了顯式濾波處理［14］。求解得到下一時刻的波場值Un+1后，對其應用模板長度為2N+1 個格點的濾波算子進行濾波處理。濾波算子可以表示為（以變量q為例）

其中

3 自適應網(wǎng)格有限差分法求解彈性波方程的實現(xiàn)

在有限差分模擬方法中實現(xiàn)自適應網(wǎng)格，需要選擇或實現(xiàn)以下幾個關鍵部分：數(shù)據(jù)結構，多級網(wǎng)格生成，不同級別網(wǎng)格間的波場交換，以及施加邊界條件和震源。

3.1 數(shù)據(jù)結構

用于實現(xiàn)自適應網(wǎng)格的數(shù)據(jù)結構方案主要有兩種：樹狀結構的自適應網(wǎng)格［15］和塊狀結構的自適應網(wǎng)格［16］。樹狀結構的自適應網(wǎng)格使用四叉樹（二維）或八叉樹（三維）的數(shù)據(jù)結構。塊狀結構的自適應網(wǎng)格使用的是常規(guī)的線性表結構，將細化的目標區(qū)域使用一系列形狀類似，包含格點數(shù)量相近的矩形網(wǎng)格區(qū)進行離散［10］，負載均衡將在這些矩形區(qū)域的基礎上進行。和樹形結構相比，塊狀結構的自適應網(wǎng)格實現(xiàn)方式較簡單，能夠訪問離當前格點較遠的波場值，與本文采用的八階格式的長差分模板適應，因此本文采用塊狀自適應網(wǎng)格結構。我們采用Berkeley Lab 開發(fā)的自適應網(wǎng)格框架AMReX［16］解決自適應網(wǎng)格復雜的數(shù)據(jù)結構帶來的網(wǎng)格生成和管理問題，生成網(wǎng)格的具體流程在下一節(jié)中詳細介紹。

3.2 網(wǎng)格生成

自適應網(wǎng)格框架AMReX 提供了自適應網(wǎng)格生成的函數(shù)接口，使用時只需以數(shù)組形式給出每個網(wǎng)格點的某個物理量取值，以及以該物理量取值為標準的細化閾值，則AMReX 函數(shù)自動將需要細化的網(wǎng)格進行標記（設置特定取值）。圖2 展示了三級自適應網(wǎng)格生成過程。首先我們要計算基網(wǎng)格，也就是最大的網(wǎng)格步長，使用基網(wǎng)格劃分整個計算模型（圖2（a））；然后按照某種標準對需要細化的目標網(wǎng)格區(qū)域進行標記（圖2（b）），其中黃色點劃線表示標記出的在基網(wǎng)格基礎上需要細化的區(qū)域，使用集群算法將被標記的所有點整合，劃分整理成組成下一級網(wǎng)格的網(wǎng)格塊（圖2（c））。每一級網(wǎng)格都重復此步驟進行生成（圖2（d）），直到達到需要的最小網(wǎng)格步長。圖2（d）中覆蓋表層高速區(qū)的細網(wǎng)格將在后文自由表面條件實現(xiàn)中進行介紹。需要并行計算時，程序將按照每個網(wǎng)格塊的點數(shù)，將各網(wǎng)格塊分配給各個進程，進行負載均衡。在本文中，細化標準選擇為每級網(wǎng)格的步長將保證每波長格點數(shù)（PPW）滿足減少數(shù)值頻散的精度需求。考慮到使用的八階精度的中心差分格式和波傳播路徑的長度，本文采用PPW=6以兼顧精度與效率。

圖2 三級自適應網(wǎng)格生成過程（灰度圖表示速度模型，較淺色塊表示低速的需要細化的區(qū)域）Fig.2 Scheme of AMR grid generation with three levels

對于給定的速度模型，最小網(wǎng)格步長hmin可以按照模型最小速度vmin，震源最大頻率fmax和PPW來計算，即

本文采用固定的相鄰級別網(wǎng)格步長比為2。繼續(xù)計算其他各級粗網(wǎng)格的步長，理論基網(wǎng)格步長hmax可以使用模型最大速度vmax來估計

實際使用的基網(wǎng)格步長應該小于等于這個值。本文計算基網(wǎng)格步長時采用的速度值略小于模型最大速度，目的是保證基網(wǎng)格作為最高級別網(wǎng)格的覆蓋范圍不會太小。

3.3 粗細網(wǎng)格信息交換

不同級別網(wǎng)格間的波場交換分為兩個部分：從粗網(wǎng)格波場到細網(wǎng)格波場的插值操作，目的是計算細網(wǎng)格區(qū)域邊緣點差分模板缺失的格點波場值；從細網(wǎng)格波場到粗網(wǎng)格波場的降采樣操作，目的是使計算穩(wěn)定。

對于粗網(wǎng)格到細網(wǎng)格的波場傳遞，本文采用之前相關研究中普遍采用的雙線性插值的方法［5，17-18］：

其中Wc和Wf分別是粗網(wǎng)格和細網(wǎng)格上的物理量，F(xiàn)(i)和N(i)用于計算細網(wǎng)格點周圍空間距離最近的粗網(wǎng)格點的索引，

Ai是插值系數(shù)，可以由粗細格點的空間位置計算得到［19］：A1= 9/16，A2=A3= 4/16，A4= 1/16.

對于細網(wǎng)格到粗網(wǎng)格的波場傳遞，本文采用以下公式由細網(wǎng)格波場獲得粗網(wǎng)格點的波場值：

其中系數(shù)取值A1=A2=A3=A4= 1/4.

粗細網(wǎng)格邊界處各類點的空間分布見圖3。

圖3 自適應網(wǎng)格波場空間分布Fig.3 Wavefield variable locations in AMR

3.4 邊界條件和震源施加

本文采用海綿層吸收邊界［20］吸收來自人工截斷邊界的虛假反射。在不同位置采用相同物理長度的吸收層。采用應力鏡像法施加自由表面邊界條件。為確保穩(wěn)定性，在接近地表的位置速度和應力計算均采用降階處理。由于自適應網(wǎng)格中不同等級網(wǎng)格點空間位置不會相互重合，如果地表層的介質（按照波速）需要不同等級的網(wǎng)格進行離散，地表格點的深度將無法一致，所以我們將自適應網(wǎng)格的地表層處理成同級的、地表最小速度需要的網(wǎng)格等級（如圖2（d））。

本文施加震源時，采用同時在震源點周圍的高斯分布上的點同時加源的平滑處理，避免震源處空間差分的奇異性：

其中(xs，zs)為震源點空間坐標，(x，z)為周圍點坐標（xs-x≤dx），f(t)是震源時間函數(shù)。

本文使用自適應網(wǎng)格進行彈性波有限差分模擬的主要步驟為：①輸入計算參數(shù)；②生成自適應網(wǎng)格；③從基網(wǎng)格到最高級網(wǎng)格，依次遍歷每級網(wǎng)格，求解彈性波方程；④使用粗網(wǎng)格點的波場值計算細網(wǎng)格邊緣模板缺失的值；⑤使用細網(wǎng)格點的波場值重新計算與細網(wǎng)格重疊的粗網(wǎng)格的部分波場值；⑥重復步驟3 到步驟5 直到最終的時間步；⑦輸出結果。

4 數(shù)值算例

我們使用4個模型對本文提出的自適應網(wǎng)格方法的準確性和穩(wěn)定性進行測試：均勻半空間模型、層狀模型、盆地模型和SEAM 模型。自適應網(wǎng)格中的網(wǎng)格級別從0開始計數(shù)，使用h0表示基網(wǎng)格空間步長，使用hn表示第n級網(wǎng)格的空間步長。

4.1 均勻半空間模型和層狀模型

我們首先使用自適應網(wǎng)格模擬均勻半空間模型中彈性波的傳播，并將結果和常規(guī)的使用均勻細網(wǎng)格的結果作比較。半空間均勻模型中只有地表是不連續(xù)面，可以排除其他因素的影響，通過簡單的震相直觀地顯示波場在不同級別的網(wǎng)格傳播時產(chǎn)生的誤差。模型縱波速度為6 900 m/s，橫波速度為4 000 m/s，密度為2 600 kg/m3。使用爆炸源，震源時間函數(shù)為主頻為2.67 Hz 的雷克子波。使用兩級的自適應網(wǎng)格對模型進行劃分，由于震源最高頻率是6.67 Hz（約2.5倍主頻），我們設置基網(wǎng)格步長h0= 100 m（PPW=6），h1= 0.5h0=50 m 呈層狀分布，粗細網(wǎng)格交界處位于z=8 100 m（圖4）。設置時間步dt=0.002 s，使其滿足穩(wěn)定性條件。使用厚度為800 m 的指數(shù)吸收邊界吸收來自計算區(qū)域的截斷邊界的虛假反射。我們使用全局細網(wǎng)格（h= 50 m）的計算結果作為參考解。

圖4 半空間均勻速度模型Fig.4 Homogeneous half-space velocity model

計算得到的vx分量在1.2 s和2 s時的波場快照如圖5。由圖5可知，在自適應網(wǎng)格的波場快照中，在粗細網(wǎng)格邊界處沒有明顯的由波場交換引起的虛假反射。和全局細網(wǎng)格相比，只有在吸收邊界處有明顯差異。由于不同位置粗細網(wǎng)格的吸收邊界物理長度一致，細網(wǎng)格處吸收層數(shù)多，衰減程度高于粗網(wǎng)格。圖6展示了自適應網(wǎng)格和全局細網(wǎng)格得到的vx和vz分量（接收點坐標見圖4），模擬記錄的最高頻率是震源子波主頻的2.5倍。在每個接收點處，二者波形無明顯差異。我們計算了各個接收點的兩組波形從0 s 到6 s 的單值包絡誤差EM（single-valued envelope misfit）和單值相位誤差PM（single-valued phase misfit）［21-22］。這種誤差標準基于連續(xù)小波變換，可以區(qū)分包絡誤差和相位誤差。EM和PM表達式分別為

圖5 均勻半空間模型彈性波模擬的波場快照Fig.5 Snapshots at 1.2 s and 2 s for the homogeneous half-space model

其中WREF(t，f)是參考解進行連續(xù)小波變換的結果，ΔE和ΔP分別是自適應網(wǎng)格計算得到的波形經(jīng)連續(xù)小波變換的結果W(t，f)與WREF(t，f)的包絡誤差和相位誤差，具體公式見文獻［21］。文獻［22］中將同時滿足|EM|＜0.2 和|PM|＜0.22 的兩組波形的擬合程度評價為“excellent”。我們將每個接收點的包絡誤差和相位誤差分別標記在圖6中。可以看到，因為2 號和5 號檢波點的vx分量能量微弱，吸收邊界沒有完全吸收的能量占主導，自適應網(wǎng)格吸收邊界的吸收效果和全局細網(wǎng)格的吸收效果有差異，所以EM 和PM 值較大。與均勻細網(wǎng)格不同，自適應網(wǎng)格的震源位于粗網(wǎng)格中。5 號檢波器距離震源較近，所以vz分量誤差較大。其他4 個檢波點的vx分量及全部檢波點的vz分量的誤差都遠小于“excellent”等級的誤差標準。由此說明，使用自適應網(wǎng)格模擬彈性波傳播時，模擬結果比較準確。

圖6 自適應網(wǎng)格和全局細網(wǎng)格的地震記錄（接收點位于圖4）Fig.6 Comparison of the velocity components calculated with the AMR grid and those with the uniform grid at receivers in Fig.4

另外，我們將自適應網(wǎng)格震源位于細網(wǎng)格（震源坐標為10 025 m 和5 025 m）和粗網(wǎng)格（震源坐標為10 050 m 和12 050 m）情況下的波場計算至第10萬步（時間=200 s），將各接收點的分量記錄展示在圖7 中。圖7 顯示在長時程模擬的過程中波形是穩(wěn)定的，表明了自適應網(wǎng)格算法的穩(wěn)定性。

圖7 自適應網(wǎng)格均勻半空間模型長時模擬結果Fig.7 The results of long time simulation using AMR for homogeneous half-space model

自適應網(wǎng)格提高地震波模擬計算效率，一方面是通過減少網(wǎng)格數(shù)，更重要的一方面是自適應網(wǎng)格可以在滿足穩(wěn)定性的條件下，顯著增大了時間步長。我們通過層狀速度模型來展示自適應網(wǎng)格對時間步長的增加效果。所用的三層速度模型由淺到深的橫波速度分別是800 m/s，1 600 m/s，3 200 m/s，縱波速度分別是1 386 m/s，2 771 m/s和5 543 m/s，密度分別是2 300 kg/m3，2 500 kg/m3和2 600 kg/m3。爆炸源的震源時間函數(shù)為主頻為2 Hz 的雷克子波。使用三級自適應網(wǎng)格，由淺到深的基網(wǎng)格，一級網(wǎng)格和二級網(wǎng)格的空間步長分別是100 m，50 m 和25 m。圖8 展示了速度模型、網(wǎng)格分布和觀測系統(tǒng)。圖8中黑色實線表示速度界面，黑色虛線表示不同級別間網(wǎng)格的邊界。覆蓋從淺到深的層的最細網(wǎng)格分別是二級網(wǎng)格，一級網(wǎng)格和基網(wǎng)格。由于穩(wěn)定性條件，如果采用全局均勻細網(wǎng)格計算，穩(wěn)定的時間步長為dt=0.005 s。而自適應網(wǎng)格由于在最高速區(qū)域的網(wǎng)格步長更大，穩(wěn)定的時間步長為dt=0.01 s。

圖9 比較了接收點處（見圖8）自適應網(wǎng)格和全局細網(wǎng)格計算得到的速度分量?？梢钥吹綌M合程度較高。時間效率方面，與全局細網(wǎng)格相比，自適應網(wǎng)格將需要進行計算的網(wǎng)格點減少到約26%，由于時間步長大，計算時間減少到約15%。該測試表明自適應網(wǎng)格可以從網(wǎng)格數(shù)和時間步長兩個方面提高模擬的效率。

圖8 層狀模型及網(wǎng)格分布Fig.8 The layer model and the refinement scheme

圖9 層狀模型自適應網(wǎng)格和均勻細網(wǎng)格計算結果比較Fig.9 Comparisons of the velocities in x-and z-direction calculated with the AMR and uniform finer grid for the layer model

4.2 盆地模型

沉積盆地內部經(jīng)常存在較厚的低速風化層，地震波傳播進入盆地后，其能量會陷入在盆地中，導致地震動幅值增大和地震動持續(xù)時間增加的現(xiàn)象。準確模擬地震波在盆地模型中的傳播是強地面運動模擬的關鍵。使用數(shù)值方法模擬地震波時，由于盆地內速度往往較低，需要使用小網(wǎng)格，導致計算量增加。已有研究表明，即使只計算低頻地震波，不使用等效介質參數(shù)處理時，粗網(wǎng)格計算的低頻地震動和細網(wǎng)格計算結果也會存在較大差異［23］。自適應網(wǎng)格可以幫助解決這類問題：在低速區(qū)使用較小空間步長，在高速區(qū)使用較大空間步長，由此達到準確模擬波場的同時，節(jié)省計算資源的目的。本例中，我們使用自適應網(wǎng)格模擬彈性波在盆地中的傳播，將結果和常規(guī)的使用均勻細網(wǎng)格的結果作比較。

用于測試的速度模型中含有兩個低速盆地，形狀分別為直角梯形和矩形，如圖10 所示。兩個盆地的介質的縱波波速為2 500 m/s，橫波波速為650 m/s，密度為2 300 kg/m3；圍巖的縱波波速4 500 m/s， 橫 波 波 速 為2 600 m/s， 密 度 為2 600 kg/m3。

圖10 （a）盆地模型，（b）各級網(wǎng)格分布Fig.10 （a）The basin model，（b）Diagram of grid refinement

爆炸源采用震源時間函數(shù)為主頻4.34 Hz 的雷克子波。最大與最小的橫波速度比值為4，本次試驗采用三級的自適應網(wǎng)格進行計算?？紤]到震源的最大頻率和模型速度值，我們設置基網(wǎng)格步長h0=40 m 覆蓋整個計算區(qū)域，第2 級網(wǎng)格（h2=10 m）對兩個低速盆地進行劃分，基網(wǎng)格和第2級網(wǎng)格間有第1級網(wǎng)格過渡，以滿足相鄰網(wǎng)格級別連續(xù)的要求。具體做法是設置過渡層的最小寬度，每個網(wǎng)格塊的范圍由自適應網(wǎng)格框架自動生成。第二級網(wǎng)格包含部分高速區(qū)，避免低速區(qū)部分格點的差分模板包含粗網(wǎng)格點。根據(jù)速度模型自動生成的各級網(wǎng)格的空間分布如圖10（b），深灰線、紅線和綠線分別表示基網(wǎng)格，1 級和2 級網(wǎng)格塊的邊界。采用這種網(wǎng)格步長設置方式，滿足壓制頻散誤差的每波長格點PPW≥6 的要求。設置時間步長dt=0.002 s 滿足穩(wěn)定性要求。震源點和檢波點的位置如圖10（a），另外，我們在地表設置測線。我們將自適應網(wǎng)格的模擬結果和全局細網(wǎng)格（空間步長h=h2= 10 m）得到的參考解進行比較。

圖11 比較了自適應網(wǎng)格和全局細網(wǎng)格在時刻t=1.2 s 和時刻t=2 s 的波場快照。為排除吸收邊界的影響，波場快照內沒有包含吸收層區(qū)域。由圖可知，除吸收邊界外，使用全局細網(wǎng)格和自適應網(wǎng)格得到的結果相差不大；在t=2 s 時，可以看到地震波在盆地內的能量很強，而盆地周圍的波場能量微弱。圖12 比較了幾個接收點處（見圖10）自適應網(wǎng)格得到的速度分量和全局細網(wǎng)格得到的速度分量。1～3 號檢波點位于盆地內（2 級細網(wǎng)格），4～6 號檢波點位于高速圍巖區(qū)域（粗網(wǎng)格）。自適應網(wǎng)格和全局細網(wǎng)格的指數(shù)吸收邊界的吸收效果很難做到完全一致，這種差異在波形上有所體現(xiàn)。整體擬合程度較高。為了定量地評估兩種方法得到的波形的擬合程度，我們使用從0～24 s的各接收點的x和z方向的速度分量計算了EM 和PM值，結果也展示在圖12 中。所有檢波點的EM 和PM 值均小于“excellent”水平的誤差閾值。盆地內部1～3 號和6 號檢波點的vx的EM 值較大：6 號檢波點可能是因為地震波振幅較小，受吸收邊界吸收效果差異的影響比較大。1～3 號檢波點vx分量的主要能量比較集中在波場能量在盆地內多次反射的時段，對整體誤差影響大，而4～6號的vx分量主要能量集中在直達波，自由表面反射和前幾次盆地下界面反射的時段；vz分量的主要能量出現(xiàn)的時段在1～6號檢波點相似。自適應網(wǎng)格波場的誤差隨傳播逐漸積累，所以相比4～6 號檢波點，1～3 號檢波點的vx誤差明顯偏大，vz誤差接近。圖13是在地表測線處使用全局細網(wǎng)格和自適應網(wǎng)格的計算得到的水平速度分量和垂直速度分量，可以看到，自適應網(wǎng)格計算結果和全局細網(wǎng)格的計算結果無明顯差異；另外，由于地震波能量被盆地捕捉，地表測線的各檢波點（除在盆地范圍外，橫坐標10 km左右）在很長時間內一直存在地震波能量。

圖11 盆地模型彈性波模擬的波場快照（同一時刻的波場快照的色標的最大值調整至全局細網(wǎng)格波場快照振幅的最大值）Fig.11 Snapshots at 1.2 s and 2 s for the basin model(The maximum values in grayscale bars of the results at the same moment have been adjusted to the maximum value for the corresponding result using the uniform grid.)

圖12 使用自適應網(wǎng)格和全局細網(wǎng)格模擬彈性波在盆地模型中的傳播得到的地震記錄，接收點見圖10（a）Fig.12 The velocity components calculated with theAMR grid and those with the uniform grid at receivers in Fig.10(a).

圖13 使用自適應網(wǎng)格模擬彈性波在盆地模型中的傳播得到的地表測線的結果，和全局細網(wǎng)格得到的參考解做比較Fig.13 Comparison of the velocity components for the receivers at the survey line calculated using the AMR grid and those using the uniform grid

時間效率方面，與全局細網(wǎng)格相比，自適應網(wǎng)格將需要進行計算的網(wǎng)格點減少到18%，計算時間也減少到約18%。自適應網(wǎng)格中除求解彈性波控制方程之外消耗計算時間的來源，主要包括粗細網(wǎng)格波場交換和自適應網(wǎng)格復雜的存儲結構（可能破壞內存訪問的局部性，影響內存訪問效率）。這些因素在本次測試中幾乎沒有影響。本次測試的結果證實了自適應網(wǎng)格模擬的準確性和高效性，以及在強地面運動模擬場景下的適用性。

4.3 SEAM模型

為測試自適應網(wǎng)格在更復雜的模型中的有效性，采用的速度模型來自國際勘探地球物理學家學會的先進模擬計劃第一期（SEAM Phase I，SEG advanced modeling program）［24］。SEAM 一期的鹽下三維模型（SEAM Phase I subsalt earth model）取自墨西哥灣的一個深水鹽域，本次測試使用的是這個模型的一個二維剖面，截取的是北坐標等于23 900 m 的位置。本次測試在原模型的基礎上進行了部分調整：將原模型的橫波速度為0的海水部分改成橫波速度600 m/s。SEAM 模型的內部速度變化復雜（見圖14）。圖14中偽彩色圖表示橫波速度，灰線、紅線和綠線分別表示基網(wǎng)格、一級網(wǎng)格和二級網(wǎng)格塊的邊界。使用爆炸源，震源時間函數(shù)是主頻為4 Hz的雷克子波。

自適應網(wǎng)格基網(wǎng)格步長h0= 40 m，1級網(wǎng)格步長h1= 0.5h0= 20 m，2 級網(wǎng)格步長h2= 0.5h1=10 m，1 級和2 級網(wǎng)格的劃分標準是橫波波速分別小于2 400 m/s 和1 200 m/s，使得每波長網(wǎng)格點數(shù)滿足壓制數(shù)值頻散的要求。基網(wǎng)格和2級網(wǎng)格中間存在1 級網(wǎng)格進行過渡。圖14 展示了三級自適應網(wǎng)格的空間分布，可以看到淺層的海水和低速沉積層由第2級網(wǎng)格覆蓋，鹽丘和高速的沉積層則由基網(wǎng)格進行劃分，可以看到網(wǎng)格步長變化和高速區(qū)的不規(guī)則邊界較為貼合，體現(xiàn)了自適應網(wǎng)格的靈活性。另外，為了避免網(wǎng)格變化產(chǎn)生的誤差與正常的反射界面產(chǎn)生的反射波耦合，也避免細網(wǎng)格差分模板的部分點延伸到粗網(wǎng)格范圍內，細網(wǎng)格覆蓋范圍略微超出了低速區(qū)，覆蓋了周邊的一部分高速區(qū)。為滿足計算穩(wěn)定性要求，設置時間步長dt=0.000 8 s。震源位于（8 905 m，0 m），4～9號檢波點位于鹽丘深部高速區(qū)（基網(wǎng)格），12 號檢波點位于鹽丘淺部（基網(wǎng)格），10、11、13和14號檢波點位于地表低速區(qū)。本測試使用全局細網(wǎng)格（h=10 m）計算得到的波形為參考解。

圖14 SEAM模型及網(wǎng)格Fig.14 SEAM model and diagram of grid refinement

圖15 和16 比較了接收點處（見圖14）自適應網(wǎng)格和全局細網(wǎng)格計算得到的速度分量?？梢钥吹綌M合程度較高。為定量評估自適應網(wǎng)格的結果和參考解的擬合程度，分別計算了EM 和PM，結果見圖15 和16。與其他檢波點相比，5～8 號檢波器vx分量的主要能量集中在記錄后段，受記錄后段的誤差影響更大；vz分量的能量分布沒有明顯區(qū)別。自適應網(wǎng)格誤差隨傳播過程逐漸積累，記錄后段誤差更大。所以5～8號檢波點記錄vx分量誤差更大。所有檢波點自適應網(wǎng)格和均勻細網(wǎng)格的波形擬合程度均可以達到“excellent”。時間效率方面，與全局細網(wǎng)格相比，自適應網(wǎng)格將需要進行計算的網(wǎng)格點減少到約50%，計算時間也減少到約50%。本次測試的結果證實了自適應網(wǎng)格模擬的準確性和高效性，以及在復雜模型地震波模擬中的適用性。

圖15 1～7號檢波點處（如圖14）自適應網(wǎng)格模擬結果和參考解的比較Fig.15 Comparison of the velocity components calculated with the AMR grid and those with the uniform grid at 1-7 receivers in Fig.14

圖16 SEAM模型在8～14號檢波點處（如圖14）自適應網(wǎng)格模擬結果和參考解的比較Fig.16 Comparison of the velocity components calculated with the AMR grid and those with the uniform grid at 8-14 receivers in Fig.14

5 結 論

為了降低網(wǎng)格在模型內高速區(qū)域的過采樣，提高計算效率，我們將自適應網(wǎng)格和同位網(wǎng)格有限差分法相結合進行彈性波數(shù)值模擬。

我們采用4個模型對方法的準確性和效率進行了驗證：均勻模型、層狀模型、盆地模型和復雜的SEAM 模型。使用自適應網(wǎng)格的結果和使用全局均勻細網(wǎng)格的結果進行比較?；跀?shù)值實驗結果得到了如下結論：

1）使用自適應網(wǎng)格的計算結果滿足精度要求，相位誤差和幅值誤差在10%左右；

2）使用自適應網(wǎng)格計算過程中復雜的數(shù)據(jù)結構冗余較小，計算時間大致與網(wǎng)格點數(shù)成正比；

3）自適應網(wǎng)格通過降低總網(wǎng)格數(shù)，和通過增大高速區(qū)網(wǎng)格大小提高滿足穩(wěn)定性的時間步長，使總體計算效率提升。

此外，施加自由表面邊界條件時，由于不同等級的網(wǎng)格點不重合，需要將地表全部用細網(wǎng)格覆蓋。今后的研究將改進這一點，減少自由表面處不必要的細網(wǎng)格點。本文所述的自適應網(wǎng)格方案可以推廣至震源動力學模擬。