瑞雷波數值模擬中的數值頻散校正策略及實例分析

豐赟 ，周竹生，沙椿

(1. 中南大學 地球科學與信息物理學院，湖南 長沙，410083；(2. 四川中水成勘院 工程勘察有限責任公司，四川 成都，610072)

瑞雷波勘探作為一種新物探方法，近年來發展迅速，其具有經濟、快捷、高效的特點，被廣泛應用于工程地質勘察和工程無損檢測。目前，人們對瑞雷波的理論研究主要集中在頻散曲線的正演和反演方面。為更準確、全面地了解瑞雷波的傳播特性，近年來，人們對基于波動方程的瑞雷波數值模擬進行了研究，如：周竹生等[1]采用交錯網格有限差分法對水平層狀各向同性彈性介質進行了全波場模擬，并對瑞雷波的傳播特征進行了探討，取得了較好的模擬效果，但差分精度較低導致模擬的波場出現了較明顯的數值頻散現象；Xu等[2]在模擬中采用AEA法對自由邊界進行處理，由于采用空間二階差分精度，為滿足精度要求，必須減小網格間距以增加最小波長中的網格數，導致計算量增加；Bohlen[3]對普通交錯網格和旋轉交錯網格進行了討論，發現由于差分階數較低，要達到一定的計算精度同樣需要增加最小波長中的網格數；熊章強等[4]的研究主要集中在邊界處理方面，采用PML吸收邊界條件以消除邊界反射，取得了較好的效果，但其采用2×12階高精度交錯網格有限差分法必然導致計算效率降低。以上研究大多采用有限差分方法對瑞雷波進行數值模擬，而數值(網格)頻散是有限差分模擬不可回避的問題，究其原因，是因為基于波動方程的有限差分求解過程通常是利用1個離散化的差分方程去逼近波動方程，從而使得相速度變成離散空間間隔的函數。因此，當每一波長內空間采樣太少(即粗網格)時，就會產生數值頻散，其實質是一種因離散化求解波動方程而產生偽波動[5]。為了壓制數值頻散，Dablain[6]對聲波二階空間差分的數值模擬進行了分析，指出網格大小和傳播方向是影響頻散的重要因素；蔡其新等[7]就高階差分、優化差分參數和通量校正傳輸(FCT)等方面對壓制頻散進行了討論，提出了最小頻散算法；吳國忱等[8]從空間離散、時間離散和算子近似3方面對TI介質中波場模擬產生的頻散進行分析，并對壓制數值頻散的方法進行討論；孫成禹等[9]就空間采樣間隔、傳播距離及速度等因素對頻散的影響進行了定量分析，并指出了由頻散導致的假頻現象及其主要特點。在此，本文作者針對瑞雷波有限差分模擬中的數值頻散現象，通過通量校正傳輸方法對其進行壓制，分析和解釋波場特征，并對計算效率進行討論。

1 二維一階彈性波波動方程及其差分格式

1.1 波動方程

在二維均勻各向同性彈性介質中，一階速度-應力標量彈性波波動方程組為[10]：

其中：vx和vz分別為質點振動速度在x和z方向的分量；τxx和 τzz分別為質點在 x和 z方向的正應力；τxz為質點在xz平面內的剪應力；ρ為介質的密度；μ為介質的剪切摸量，與λ一并統稱為拉梅常數。

1.2 差分格式

交錯網格是指把速度v和應力τ分別定義于2套不同網格上的網格系統，其優點在于能夠很好地處理一階波動方程中速度和應力的耦合關系。圖1所示即為其網格剖分方案。

圖1 交錯網格剖分示意圖Fig.1 Scheme of staggered

本文采用時間二階和空間四階的差分格式，差分格式如下：

其中：Dt，Dx和 Dz分別為對 t，x和 z的差分離散算子的集合；a1=1.125；對波動方程(1a～1e)進行差分離散，離散方程組如下：

其中：U和V分別為vx和vz的離散量；R，T和H分別為 τxx，τzz和 τxz的離散量；B 為 ρ-1的離散值。

2 邊界條件及穩定性條件

2.1 自由邊界

在進行瑞雷波數值模擬時，自由邊界能否得到很好處理將直接決定著模擬效果。本文在自由邊界采用空間四階差分精度，用鏡像法和零速度法進行處理[11]。

對于二維模型，自由邊界條件(z=0)為τxz=0，τzz=0。當采用高階差分格式時,為了保證在自由邊界條件(z=0)成立，則 τzz和 τxz在高階差分算子的長度上一定是關于z=0反對稱的，這就是鏡像原理。設定vx，τxx和τzz位于自由界面，并在自由界面上設置2層真空層(vs和vp為0 m/s)，此時，自由界面存在一反射界面，從而使得球面波在界面反射產生瑞雷波。速度分量U和V計算時滿足：

2.2 吸收邊界

吸收邊界條件廣泛應用于無邊界介質中的波場模擬，以減少來自人工邊界的反射。關于吸收邊界條件已有較多的討論，完全匹配層(PML)吸收邊界條件是消除邊界反射的理想方法之一，其實質是理論上假設存在一種各項異性有耗介質(PML層)，適當選取參數可使任意頻率、任意極化的波以任意入射角通過它與自由空間的交界面時相速度和特征阻抗均不變，從而達到與自由空間完全匹配，入射波進入 PML層后迅速衰減，理論上不產生任何反射。本文采用文獻[4]中的方法對人工邊界進行處理。

2.3 穩定性條件

二維均勻各向同性彈性介質，差分精度為O(Δt2+Δx4)，Δx=Δz，其穩定性條件為[12]：

3 FCT方法

FCT(Flux-corrected transport)方法又稱通量校正傳輸，由 Boris等[13]提出。該方法有效壓制了粗網格情況下差分計算產生的數值頻散和數值震蕩。FCT方法是基于守恒型方程的差分格式，要求網格間交界面處守恒物理量的通量相互抵消，即交界面兩側一邊流出和另一邊流出的通量相等。其基本原理是：假設所有的極值點都是由數值頻散引起的，對所有網格點進行漫射校正處理，達到消除數值頻散的效果；再對非局部極值點進行補償的反漫射校正，使任何不需要漫射的地方得到恢復，反漫射校正是非線性的。Tong等[14]將 FCT方法應用于基于波動方程的波場數值模擬，很好地壓制了數值頻散；楊頂輝等[5,15]將FCT方法應用于求解各向異性波動方程組；鄭海山等[16]采用有限差分法，結合FCT方法和幅值限制器對橫向各向同性介質的二維非線性彈性波進行數值模擬，均有效地壓制了傳統有限差分數值模擬中的數值頻散,取得了較好的模擬效果。

對于一階速度-應力波動方程，在應用FCT方法對差分數值解進行校正的過程中，只需對和進行計算，因為在完成對的校正后，通過方程3(c)～3(e)對進行了校正。FCT方法計算步驟如下(它適用于二維一階速度-應力波動方程的有限差分計算)。

(1) 利用有限差分公式 3(a)和 3(b)計算所對應的差分數值解

(2) 計算第k-1時間層的漫射通量：

其中：η1為漫射參數，它為常量或為1個線性函數，其取值取決于有限差分階段頻散誤差，在數值試驗中，控制0.008≤η1≤0.100(當η1小于0.008時，已失去壓制頻散的效果；當其大于0.100時，真實波場已被平滑模糊不清，分辨率顯著下降)；0.01≤η1≤0.05段比較合適，模擬效果較好，本文中η1=0.010。

(3) 對第k+1時間層的反漫射進行校正：

η2的取值與η1的取值不同。這是由于傳統有限差分運算引起的頻散，又有人為加入的漫射。反漫射運算不僅要補償人為加入射引起的損失，而且要補償傳統有限差分運算帶來的振幅損失。因此，η2的取值應比η1大10%～15%。本文中，η1=0.012。

4 數值模擬算例

為更加接近實際地模擬淺層各向同性彈性介質的全波場，選取較小的網格間距，給定較低的縱橫波速，但低速地層更容易產生數值頻散和假頻[11]，因此，采取有效措施以消除數值頻散很有必要。設計模型寬×深為200 m×150 m，縱橫向空間采樣間隔均為1 m，時間采樣間隔為0.2 ms，震源采用主頻30 Hz的Ricker子波，置于第100號檢波點處，埋深0 m。

4.1 模型1——均勻半空間介質模型

均勻介質物性參數分別為：密度ρ=2.0 g/cm3，縱波速度vp=1.5 km/s，橫波速度vs=800 m/s。圖2(a)和(b)所示是未進行FCT方法校正的快照圖和單炮記錄，圖3(a)和(b)所示是經FCT方法處理后的快照圖和單炮記錄(圖中，D，SW，S，P，P-S和PP分別代表直達波、面波、橫波、縱波、P-S頭波和縱波反射波，記錄和快照顯示的都是vz分量，下同)。

圖2 未加入FCT方法的60 ms時的快照(a)和模擬單炮記錄Fig.2 Snaps at 60 ms and synthetic single-shot record without using FCT method

對比圖2和圖3可以看到：傳統有限差分模擬造成比較嚴重的頻散情況，而加入FCT方法后，數值頻散得到有效壓制，更接近于真實波場。圖3顯示：在均勻半空間介質中，瑞雷波不發生頻散現象。

4.2 模型2——含低速夾層的介質模型

表1所示為含低速夾層的層狀介質模型參數，圖4(a)和(b)所示是未進行FCT方法處理的快照圖和單炮記錄，圖5(a)和(b)所示是經FCT方法處理后的快照圖和單炮記錄。

圖3 加入FCT方法后的60 ms時的快照和模擬單炮記錄Fig.3 Snaps at 60 ms and synthetic single-shot record by using FCT method

表1 含低速夾層的層狀介質的模型參數Table 1 Model parameters of layered media with a low velocity interlayer

圖4 未加入FCT方法的100 ms時刻的快照(a)和模擬單炮記錄Fig.4 Snaps at 100 ms and synthetic single-shot record without using FCT method

不同頻率的波在介質中以不同的速度傳播的現象稱為波的頻散，頻率ω和波數k之間的關系ω(k)稱為頻散關系。波的頻散主要包括數值(網格)頻散和物理頻散。數值頻散是由于網格離散的截斷誤差導致群速度和相速度依賴于波的頻率，它是一種偽波動，表現為波在傳播過程中波前形狀發生變化,并且逐漸散開；物理頻散是指由于物理介質不同，波傳播相速度隨波數發生變化的現象。瑞雷波在層狀介質中傳播發生的頻散屬于物理頻散。

對比圖 4(a)和圖 5(a)可見：對于圖 4(a)中 2個矩形區域，由于各層之間的多次反射造成振蕩和網格頻散相混淆，不能清晰地分辨真實波場，經FCT方法處理后，圖 5(a)中各反射波層次比較清楚，網格頻散被壓制，而真實的反射波場被保留。

圖5 加入FCT方法的100 ms時的快照和模擬單炮記錄Fig.5 Snaps at 100 ms and synthetic single-shot record by using FCT method

對比圖4(b)和圖5(b)可見：圖4(b)中三角形區域網格頻散現象突出，干擾了真實波場；圖5(b)中，由于消除了相應區域的網格頻散，使得瑞雷波的頻散特征更加清晰。

圖 5(a)和圖5(b)中，波場快照中顯示的各類波在單炮記錄里有良好的反映，印證了經過FCT方法后，網格頻散被壓制，真實波場被凸顯。這是因為網格頻散主要是以接近于 Nyquist頻率的高頻成分引起的，而真實波場主要是低頻成分。FCT方法往往會壓制那些接近于 Nyquist頻率的數值頻散，同時恢復大部分真實波場[14]。因此，FCT方法不會破壞真正的波場振蕩，它壓制的是數值頻散，而不是在層狀介質中瑞雷波產生的頻散。

5 計算效率分析

模型3和模型4的寬×深為400 m×300 m，均勻介質物性參數分別為：密度 ρ=2.0 g/cm3，縱波速度vp=1.2 km/s，橫波速度vs=650 m/s。模型3是粗網格下的FCT模型，模型4是細化網格模型，其具體參數和計算效率見表2。為壓制頻散，模型3計算時間為473 s，模型4計算時間為1 096 s, 使用FCT方法所得計算效率比細化網格來壓制頻散的計算效率提高了 1.3倍，具有明顯的優勢；另外，在效果方面，圖 6(a)中還存在部分頻散現象，而6(b)沒有頻散現象，FCT方法壓制頻散效果比細化網格方法的優。

圖6 模型3和模型4在160 ms時的快照Fig.6 Snap at 160 ms of models 3 and 4

表2 2個模型的參數及計算時間Table 2 Parameters and computing time of models 3 and 4

6 結論

(1) 將FCT方法應用于瑞雷波的數值模擬中。該方法很好地壓制了離散網格帶來的數值頻散現象，同時保證了真實的震波場震蕩，使瑞雷波在層狀介質中傳播產生的頻散現象被凸顯，能更準確、直觀地了解瑞雷波的傳播特征。

(2) FCT方法的引入雖然增加了計算量，但在達到同一精度時，與細化網格相比，FCT方法的計算效率提高了1.3倍，為解決粗網格下瑞雷波的有限差分模擬中出現的數值頻散問題提供了一種快速有效的方法。此外，在FCT方法中漫射與反漫射步驟可以引入代碼并行計算，計算效率還將進一步提高。

(3) 建議將模擬中得到的波場數據進行處理，提取頻散曲線，反演模型地層參數，驗證數值模擬的正確性和頻散壓制的可靠性。

[1] 周竹生, 劉喜亮, 熊孝雨. 彈性介質中瑞雷面波在有限差分正演模擬[J]. 地球物理學報, 2007, 50(2): 567-573.ZHOU Zhu-sheng, LIU Xi-liang, XIONG Xiao-yu.Finite-difference modeling of Rayleigh surface wave in elastic media[J]. Chinese Journal of Geophysics, 2007, 50(2): 567-573.

[2] XU Yi-xian, XIA Jiang-hai, Miller R D. Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach[J]. Geophysics, 2007, 72(5): 147-153.

[3] Bohlen S R. Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves[J]. Geophysics,2006, 71(4): T109-T115.

[4] 熊章強, 張大洲, 秦臻, 等. 瑞雷波數值模擬中的邊界條件及模擬實例分析[J]. 中南大學學報: 自然科學版, 2008, 39(4):824-830.XIONG Zhang-qiang, ZHANG Da-zhou, QIN Zheng, et al.Boundary conditions and case analysis of numerical modeling of Rayleigh wave[J]. Journal of Central South University:Science and Technology, 2008, 39(4): 824-830.

[5] 楊頂輝, 騰吉文. 各向異性介質中三分量地震記錄的 FCT有限差分模擬[J]. 石油地球物理勘探, 1997, 32(2): 181-189.YANG Ding-hui, TENG Ji-wen. FCT finite difference modeling of three-component Seismic records in anisotropic medium[J].Oil Geophysical Prospecting, 1997, 32(2): 181-190.

[6] Dablain M A. The application of high-order differencing to scalar wave equation[J]. Geophysics, 1986, 51(1): 54-66.

[7] 蔡其新, 何佩軍, 秦廣勝, 等. 有限差分數值模擬的最小頻散算法及其應用[J]. 石油地球物理勘探, 2003, 38(3): 247-251,262.CAI Qi-xin, HE Pei-jun, QIN Guang-sheng, et al. Minimum dispersion algorithm by finite-difference numeric modeling and its applications[J]. Oil Geophysical Prospecting, 2003, 38(3):247-251, 262.

[8] 吳國忱, 王華忠. 波場模擬中的數值頻散分析與校正策略[J].地球物理學進展, 2005, 20(1): 58-65.WU Guo-chen, WANG Hua-zhong. Analysis of numerical dispersion in wave-field simulation[J]. Progress in Geophysics,2005, 20(1): 58-65.

[9] 孫成禹, 宮同舉, 張玉亮, 等. 波動方程有限差分法中的頻散與假頻分析[J]. 石油地球物理勘探, 2009, 44(1): 43-48.SUN Chen-yu, GONG Tong-ju, ZHANG Yu-liang, et al.Analysis on dispersion and alias in finite-difference solution of wave equation[J]. Oil Geophysical Prospecting, 2009, 44(1):43-48.

[10] 董良國, 馬在田, 曹景忠, 等. 一階彈性波方程交錯網格高階差分解法[J]. 地球物理學報, 2000, 43(3): 411-419.DONG Liang-guo, MA Zai-tian, CAO Jin-zhong, et al. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese Journal of Geophysics, 2000, 43(3):411-419.

[11] 裴正林. 任意起伏地表彈性波方程交錯網格高階有限差分法數值模擬[J]. 石油地球物理勘探, 2004, 39(6): 629-634.PEI Zheng-lin. Numerical modeling using staggered-grid high order finite-difference of elastic wave equation on arbitrary relief surface[J]. Oil Geophysical Prospecting, 2004, 39(6): 629-634.

[12] 董良國, 馬在田, 曹景忠. 一階彈性波方程交錯網格高階差分解法穩定性研究[J]. 地球物理學報, 2000, 43(6): 856-864.DONG Liang-guo, MA Zai-tian, CAO Jin-zhong. A study on stability of the staggered-gird high-order difference method of first-order elastic wave equation[J]. Chinese Journal of Geophysics, 2000, 43(6): 856-864.

[13] Boris J P, Book D L. Flux-corrected transportⅠ: SHASTA, a fluid transport algorithm that works[J]. Journal of Computational Physics, 1973, 11(1): 38-69.

[14] Tong Fei, Larner K. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport[J]. Geophysics, 1995, 60(6): 1830-1842.

[15] Yang D H, Liu E, Zhang Z J, et al. Finite-difference modeling in two-dimensional anisotropic media using a flux-corrected transport technique[J]. Geophysical Journal International, 2002,148(2): 320-328.

[16] 鄭海山, 張中杰. 橫向各向同性(VTI)介質中非線性地震波場模擬[J]. 地球物理學報, 2005, 48(3): 660-671.ZHENG Hai-shan, ZHANG Zhong-Jie. Synthetic seismograms of nonlinear seismic waves in anisotropic (VTI) media[J].Chinese Journal of Geophysics, 2005, 48(3): 660-671.