二維彈性波自適應網格有限差分模擬方法*

張春麗，張偉

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

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

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

1 背 景

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

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

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

圖1 不同類型的網格和格點數（灰度值表示波速，黃線表示不同網格步長區域的邊界。）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 二維彈性波控制方程與同位網格有限差分算法

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

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

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

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

其中λ和μ是Lamé常量。

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

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

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

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

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

其中

3 自適應網格有限差分法求解彈性波方程的實現

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

3.1 數據結構

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

3.2 網格生成

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

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

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

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

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

3.3 粗細網格信息交換

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

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

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

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

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

其中系數取值A1=A2=A3=A4= 1/4.

粗細網格邊界處各類點的空間分布見圖3。

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

3.4 邊界條件和震源施加

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

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

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

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

4 數值算例

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

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

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

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

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

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

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

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

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

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

自適應網格提高地震波模擬計算效率，一方面是通過減少網格數，更重要的一方面是自適應網格可以在滿足穩定性的條件下，顯著增大了時間步長。我們通過層狀速度模型來展示自適應網格對時間步長的增加效果。所用的三層速度模型由淺到深的橫波速度分別是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。爆炸源的震源時間函數為主頻為2 Hz 的雷克子波。使用三級自適應網格，由淺到深的基網格，一級網格和二級網格的空間步長分別是100 m，50 m 和25 m。圖8 展示了速度模型、網格分布和觀測系統。圖8中黑色實線表示速度界面，黑色虛線表示不同級別間網格的邊界。覆蓋從淺到深的層的最細網格分別是二級網格，一級網格和基網格。由于穩定性條件，如果采用全局均勻細網格計算，穩定的時間步長為dt=0.005 s。而自適應網格由于在最高速區域的網格步長更大，穩定的時間步長為dt=0.01 s。

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

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

圖9 層狀模型自適應網格和均勻細網格計算結果比較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 盆地模型

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

用于測試的速度模型中含有兩個低速盆地，形狀分別為直角梯形和矩形，如圖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）各級網格分布Fig.10 （a）The basin model，（b）Diagram of grid refinement

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

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

圖11 盆地模型彈性波模擬的波場快照（同一時刻的波場快照的色標的最大值調整至全局細網格波場快照振幅的最大值）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 使用自適應網格和全局細網格模擬彈性波在盆地模型中的傳播得到的地震記錄，接收點見圖10（a）Fig.12 The velocity components calculated with theAMR grid and those with the uniform grid at receivers in Fig.10(a).

圖13 使用自適應網格模擬彈性波在盆地模型中的傳播得到的地表測線的結果，和全局細網格得到的參考解做比較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

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

4.3 SEAM模型

為測試自適應網格在更復雜的模型中的有效性，采用的速度模型來自國際勘探地球物理學家學會的先進模擬計劃第一期（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中偽彩色圖表示橫波速度，灰線、紅線和綠線分別表示基網格、一級網格和二級網格塊的邊界。使用爆炸源，震源時間函數是主頻為4 Hz的雷克子波。

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

圖14 SEAM模型及網格Fig.14 SEAM model and diagram of grid refinement

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

圖15 1～7號檢波點處（如圖14）自適應網格模擬結果和參考解的比較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）自適應網格模擬結果和參考解的比較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 結 論

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

我們采用4個模型對方法的準確性和效率進行了驗證：均勻模型、層狀模型、盆地模型和復雜的SEAM 模型。使用自適應網格的結果和使用全局均勻細網格的結果進行比較。基于數值實驗結果得到了如下結論：

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

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

3）自適應網格通過降低總網格數，和通過增大高速區網格大小提高滿足穩定性的時間步長，使總體計算效率提升。

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