復(fù)雜地表起伏多重變加密網(wǎng)格有限差分法波動(dòng)方程數(shù)值模擬

摘要：

常規(guī)變加密網(wǎng)格有限差分波動(dòng)方程數(shù)值模擬方法采用水平分層加密網(wǎng)格，該網(wǎng)格剖分策略在適應(yīng)地形起伏和近地表速度結(jié)構(gòu)變化特征方面效果較差。針對(duì)該問(wèn)題，提出一種起伏多重變加密網(wǎng)格有限差分波動(dòng)方程數(shù)值模擬方法。該方法根據(jù)地形起伏和近地表低速層到高速層的速度分布特征進(jìn)行網(wǎng)格加密；采用不同網(wǎng)格中的變系數(shù)差分格式離散聲波方程，在保證波場(chǎng)模擬精度的同時(shí)兼顧計(jì)算效率；同時(shí)，為了進(jìn)一步保障地表附近波場(chǎng)模擬的精度，地表附近最細(xì)網(wǎng)格中的差分格式不做降階處理，針對(duì)位于地表以上的虛像點(diǎn)的波場(chǎng)值，提出一種融入自由地表邊界條件的法向虛像外推法。算例分析驗(yàn)證該算法對(duì)速度模型不同區(qū)域進(jìn)行的網(wǎng)格多次加密顯著提高了計(jì)算效率，以黃土塬實(shí)際模型為例，耗時(shí)為常規(guī)1 m×1 m網(wǎng)格耗時(shí)的43.3%，并可達(dá)到和細(xì)網(wǎng)格基本一致的模擬精度，模擬誤差控制在10-12范圍內(nèi)，同時(shí)表現(xiàn)出很好的近地表散射壓制和邊界吸收效果，且算法能穩(wěn)定地適應(yīng)實(shí)際復(fù)雜地表介質(zhì)。

關(guān)鍵詞：

復(fù)雜地表；波動(dòng)方程數(shù)值模擬；起伏多重變加密網(wǎng)格；法向虛像外推法；有限差分法

doi：10.13278/j.cnki.jjuese.20220185

中圖分類號(hào)：P631.4

文獻(xiàn)標(biāo)志碼：A

收稿日期：2022-06-23

作者簡(jiǎn)介： 張志佳（1998—），男，碩士研究生，主要從事地震波場(chǎng)正演模擬、地震波傳播理論與成像技術(shù)等方面的研究， E-mail： zhangzj20@mails.jlu.edu.cn

通信作者：韓復(fù)興（1981—），男，教授，博士，主要從事地震波傳播與成像、計(jì)算地球物理方面的教學(xué)與研究，E-mail： hanfx@jlu.edu.cn

基金項(xiàng)目：廣西重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2021AB30011）；國(guó)家自然科學(xué)基金項(xiàng)目（42074150，41404085）

Supported by the Key Ramp;D Program of Guangxi Zhuang Autonomous Region" （2021AB30011） and the National Natural Science Foundation of China （42074150， 41404085）

Numerical Simulation of Wave Equation Using Finite Difference Method with Rugged Variable Refined Multigrid in Complex Surface

Zhang Zhijia1， Sun Zhangqing1， Wang Ruihu2， Han Fuxing1，

Lu Jipu3， Liu Yaguang1， Cen Wenpan2， Gao Zhenghui1

1. College of" GeoExploration Science and Technology， Jilin University， Changchun 130026， China

2. Geological Survey Institute of Guangxi Zhuang Autonomous Region， Nanning 530031， China

3. Bureau of Geology and Mineral Exploration and

Development of Guangxi Zhuang Autonomous Region， Nanning 530023， China

Abstract：

The horizontal layered refined grid is adopted in the conventional finite difference method of variable refined grid for the wave equation based numerical simulation. But the method cannot well adapt to the characteristics of rugged topography and the variation of velocity structure in near-surface region. To solve this problem， a wave equation based numerical simulation using finite difference method with rugged variable refined multigrid is proposed in this paper. This method encrypts the grid based on rugged topography and velocity distribution characteristics of low-velocity layer to high-speed layer near the surface and adopts the variable-coefficient difference schemes in different meshes to discretize the acoustic wave equation， in order to ensure the accuracy of wave field simulation and calculation efficiency. Furthermore， in order to further guarantee the precision of wave field simulation near the surface， the difference schemes in the smallest grids near the earth's surface are not processed by order reduction. For the wave field of the virtual ghost points above the surface， a normal virtual ghost extrapolation method incorporating free surface boundary conditions is proposed. The numerical evaluation and example show that the algorithm significantly improves the computational efficiency by encrypting the grids in different regions of the velocity model. Taking the actual model of loess plateau as an example， the ratio of time consumption to conventional 1 m×1 m grid is 43.3% and the simulation accuracy is basically consistent with the fine grid， the simulation error is controlled within 10-12. The algorithm shows good near-surface scattering suppression and boundary absorption effect， and can stably adapt to the actual complex near-surface media.

Key words：

complex surface; wave equation based numerical simulation; rugged variable refined multigrid; normal virtual ghost extrapolation method; finite difference method

0 引言

有限差分法作為當(dāng)前數(shù)值模擬方法中最成熟、最有效的算法被廣泛使用［1-2］，其中的變加密網(wǎng)格有限差分方法更是具有兼顧精度和效率的特性。早在20世紀(jì)90年代，Jastram等［3-4］開(kāi)始對(duì)聲波方程、彈性波方程的變網(wǎng)格有限差分方法展開(kāi)研究。隨后，Moczo［5］引入插值函數(shù)在直角網(wǎng)格中實(shí)現(xiàn)網(wǎng)格間距的逐漸過(guò)渡變化。Wang等［6］、張劍鋒［7］均基于一階應(yīng)力-速度彈性波方程提出不規(guī)則網(wǎng)格方法并擴(kuò)展到黏彈介質(zhì)。Pitarka［8］嘗試在交錯(cuò)網(wǎng)格中使用該類型網(wǎng)格，模擬地表存在低速盆地的模型。雖然這種方法可以達(dá)到網(wǎng)格局部加密的要求，但由于在網(wǎng)格點(diǎn)排列的某個(gè)方向上總格點(diǎn)數(shù)固定，網(wǎng)格間距的變化沿著某一方向變化，隨著網(wǎng)格的連續(xù)，這種變化方式會(huì)橫向延伸至計(jì)算邊界，使得其他不需要細(xì)化的區(qū)域也被細(xì)化。局部的網(wǎng)格調(diào)整會(huì)導(dǎo)致網(wǎng)格結(jié)構(gòu)的變化，需要重新生成網(wǎng)格，使得網(wǎng)格對(duì)介質(zhì)變化有一定的要求，不能隨著介質(zhì)的變化靈活地變化。而后，李勝軍等［9］重新推導(dǎo)了聲波方程變網(wǎng)格有限差分格式，并將傳統(tǒng)的有限差分法與變網(wǎng)格算法從計(jì)算速度、內(nèi)存需求等方面做了對(duì)比，認(rèn)為變網(wǎng)格算法無(wú)論從模擬效率還是內(nèi)存需求方面均優(yōu)于傳統(tǒng)規(guī)則網(wǎng)格。黃超等［10-11］陸續(xù)提出一種空間網(wǎng)格大小與時(shí)間步長(zhǎng)均可任意變化的雙變網(wǎng)格高階有限差分?jǐn)?shù)值模擬方法，以及一種交錯(cuò)雙變網(wǎng)格高階有限差分?jǐn)?shù)值模擬方法。孫成禹等［12］對(duì)插值法、變差分系數(shù)法和變時(shí)間步長(zhǎng)法的精度進(jìn)行了分析。Li等［13］提出貼體網(wǎng)格上的過(guò)渡，在細(xì)網(wǎng)格和粗網(wǎng)格之間的過(guò)渡區(qū)域采用雙線性插值從粗網(wǎng)格中取值，采用高斯濾波抑制在粗網(wǎng)格中更新波動(dòng)方程長(zhǎng)時(shí)間不穩(wěn)定的人工數(shù)值噪聲。Fan等［14］改變了插值類變網(wǎng)格方法過(guò)渡帶的差分格式，根據(jù)以往的中心差分格式推導(dǎo)出一種針對(duì)過(guò)渡帶的旋轉(zhuǎn)差分格式以提高其數(shù)值模擬精度。李世中等［15］針對(duì)交錯(cuò)網(wǎng)格下變差分系數(shù)法的差分系數(shù)計(jì)算方法進(jìn)行了優(yōu)化，以提高其計(jì)算精度。姜占東等［16］將可變交錯(cuò)網(wǎng)格變差分系數(shù)法運(yùn)用于海洋地震勘探中鬼波的數(shù)值模擬。

基于起伏地形的變加密網(wǎng)格有限差分方法大部分采用水平分層加密的方式，這種加密方式不易適應(yīng)劇烈起伏地形、近地表低速風(fēng)化層、海水與海底的固液界面等地下介質(zhì)劇烈變化的分界面［17］，同時(shí)還易增加內(nèi)存和計(jì)算量。為了更好地解決這個(gè)問(wèn)題，本文提出了一種適應(yīng)復(fù)雜近地表介質(zhì)特征的起伏多重變加密網(wǎng)格方法。首先，利用法向虛像外推法對(duì)起伏自由地表邊界進(jìn)行處理，該邊界處理方法可以在地表加密網(wǎng)格區(qū)域穩(wěn)定實(shí)施，吸收邊界條件使用CPML（convolutional perfect matched layer）

吸收邊界條件［18-21］，并給出法向虛像外推法的實(shí)現(xiàn)過(guò)程；其次，對(duì)起伏多重變加密網(wǎng)格剖分方法進(jìn)行闡述并對(duì)過(guò)渡區(qū)域網(wǎng)格中的差分格式進(jìn)行推導(dǎo)；最后，對(duì)算法進(jìn)行各方面分析，并給出計(jì)算實(shí)例。

1 網(wǎng)格剖分

傳統(tǒng)多重變加密網(wǎng)格在地表部分采用細(xì)網(wǎng)格進(jìn)行剖分，隨著介質(zhì)速度遞增，網(wǎng)格間距逐漸增加（圖1a）。其網(wǎng)格剖分策略將由細(xì)網(wǎng)格過(guò)渡到粗網(wǎng)格區(qū)分成三部分：細(xì)網(wǎng)格區(qū)域、過(guò)渡區(qū)域（紅框區(qū)域）以及粗網(wǎng)格區(qū)域。這種剖分策略可以實(shí)現(xiàn)模型的逐漸過(guò)渡，避免因插值影響帶來(lái)數(shù)值不穩(wěn)定的問(wèn)題。從波動(dòng)方程高階空間差分離散公式可以看出，離散點(diǎn)處高階空間差分只需要單一方向上的離散點(diǎn)值即可，在數(shù)值模擬計(jì)算上更容易實(shí)現(xiàn)，計(jì)算效率得到提升，但在面對(duì)起伏地形和近地表復(fù)雜介質(zhì)時(shí)不具備靈活適應(yīng)性。

為了更好地描述起伏地形和復(fù)雜近地表特征，本文在多重變加密網(wǎng)格的基礎(chǔ)上，對(duì)其過(guò)渡帶區(qū)域進(jìn)行插值優(yōu)化（圖1b），應(yīng)用于起伏地表，并將其稱為起伏多重變加密網(wǎng)格。起伏多重變加密網(wǎng)格方法相比多重變加密網(wǎng)格方法增加了靈活性，保證模擬精度的同時(shí)提高了計(jì)算效率。以時(shí)間二階空間四階有限差分為例，加入起伏過(guò)渡帶后，只是在過(guò)渡帶邊界的處理上發(fā)生變化。圖1b上綠色點(diǎn)是小網(wǎng)格邊界處無(wú)法計(jì)算的點(diǎn)，計(jì)算綠色點(diǎn)需要先求得紫色點(diǎn)的波場(chǎng)值。為了計(jì)算這些點(diǎn)，在空間位置上對(duì)大網(wǎng)格進(jìn)行縱向內(nèi)插即可得到紫色點(diǎn)的波場(chǎng)值。由于僅采用縱向多點(diǎn)內(nèi)插，插值帶來(lái)的誤差很小。同時(shí)，起伏過(guò)渡帶并未因插值對(duì)計(jì)算精度帶來(lái)過(guò)多影響。

2 基于起伏地表的法向虛像外推法

Moczo［5］采用在直角網(wǎng)格中網(wǎng)格間距逐漸過(guò)渡變化的方式來(lái)處理起伏地表的問(wèn)題，用階梯狀網(wǎng)格近似地面起伏，需要足夠密集的網(wǎng)格點(diǎn)來(lái)壓制階梯狀網(wǎng)格所產(chǎn)生的角點(diǎn)散射。在地震波數(shù)值模擬計(jì)算中，地表的幾何形態(tài)會(huì)對(duì)地震波傳播的模擬結(jié)果產(chǎn)生顯著影響。常規(guī)差分格式為四階雙向五點(diǎn)差分格式，當(dāng)采取規(guī)則矩形網(wǎng)格離散建模時(shí)，起伏地表邊界會(huì)做降階

處理。考慮到復(fù)雜地表對(duì)波場(chǎng)數(shù)值模擬帶來(lái)的影響，在地表邊界的處理上不做降階處理，進(jìn)而引入虛像點(diǎn)［22］。為了表示虛像點(diǎn)的波場(chǎng)并實(shí)現(xiàn)自由地表邊界條件，在建立地表附近的差分格式時(shí)提出法向虛像外推法。

一般情況下，正交的規(guī)則矩形網(wǎng)格與需要處理的起伏自由地表邊界是不重合的，計(jì)算節(jié)點(diǎn)處于整數(shù)網(wǎng)格上而起伏自由地表邊界往往處于非整數(shù)網(wǎng)格上。因此，需要對(duì)起伏自由地表邊界處進(jìn)行一定的處理以滿足自由地表邊界條件。基于聲波方程的有

限差分法數(shù)值模擬，本文提出的法向虛像外推法首先根據(jù)空間離散階數(shù)確定待求虛像點(diǎn)數(shù)，然后對(duì)待求虛像點(diǎn)值進(jìn)行外推計(jì)算，最后根據(jù)外推虛像點(diǎn)實(shí)現(xiàn)自由地表邊界條件。

空間四階有限差分格式求取中心點(diǎn)（圖2a中的點(diǎn)M）時(shí)要求其周圍4個(gè)方向的點(diǎn)M1、M2、M3、M4、M5、M6、M7、G參與計(jì)算。其中，當(dāng)中心點(diǎn)靠近地表時(shí)，會(huì)有一些點(diǎn)延伸到自由地表之上，如圖2中的點(diǎn)G，我們稱之為虛像點(diǎn)。為了不對(duì)差分格式做降階處理，繼續(xù)在地表附近采用高階差分格式，對(duì)于地表以上的點(diǎn)，采用波場(chǎng)外推的方法用地表以下點(diǎn)的波場(chǎng)值來(lái)表示。

如圖2b所示，法向虛像外推法的實(shí)現(xiàn)原理為沿虛像點(diǎn)與地形切線的垂直方向（即地表法線方向）進(jìn)行多項(xiàng)式外推，以求取虛像點(diǎn)的波場(chǎng)表達(dá)式。首先建立局部坐標(biāo)系，以虛像點(diǎn)地表法線與切線的交點(diǎn)X0為原點(diǎn)，x方向與網(wǎng)格水平方向的3個(gè)交點(diǎn)分別為X1、X2、X3，它們的波場(chǎng)值分別為pX1、 pX2、pX3，可以由地表以下規(guī)則網(wǎng)格節(jié)點(diǎn)G1、G2、G3、G4、G5、G6處的波場(chǎng)值pG1、pG2、pG3、pG4、pG5、pG6表示：

pX1=pG1lG2X1Δh+pG2lG1X1Δh；pX2=pG3lG4X2Δh+pG4lG3X2Δh；pX3=pG5lG6X3Δh+pG6lG5X3Δh。（1）

式中：lG1X1、lG2X1、lG3X2、lG4X2、lG5X3、lG6X3分別為G1與X1、G2與X1、G3與X2、G4與X2、G5與X3、G6與X3的距離；SymbolDA@h為網(wǎng)格間距。根據(jù)自由地表邊界條件，地表處的壓力值為0，即pX0=0。

法向虛像外推法的核心思想是在X0、X1、X2和X3點(diǎn)的波場(chǎng)值已知的情況下，假設(shè)地表附近的波場(chǎng)值滿足如下局部三次函數(shù)：

f（x）=ax3+bx2+cx+d。（2）

式中，a、b、c、d為三次多項(xiàng)式的系數(shù)。

如圖2b所示，法向虛像外推的最終結(jié)果是G點(diǎn)的波場(chǎng)表達(dá)式可由多個(gè)地表以下規(guī)則網(wǎng)格節(jié)點(diǎn)（紅色點(diǎn)）的波場(chǎng)表示。法向虛像外推法處理自由地表主要有以下優(yōu)點(diǎn)：1）虛像點(diǎn)外推不需要近似于地表的幾何形態(tài)，且在地表處不做降階處理，有效抑制了模型地表階梯狀截?cái)嗨a(chǎn)生的散射和繞射等現(xiàn)象；2）該方法沿地形法線方向進(jìn)行外推，差分格式可由地下多個(gè)規(guī)則網(wǎng)格節(jié)點(diǎn)的線性組合來(lái)控制，更加符合物理意義，并可通過(guò)地下介質(zhì)中的大量網(wǎng)格點(diǎn)進(jìn)一步控制精度。

該變加密網(wǎng)格算法在模型中多次使用，使其形成網(wǎng)格間距逐漸增加的趨勢(shì)，整個(gè)過(guò)渡過(guò)程采取式（9）—式（11）格式進(jìn)行離散并將其應(yīng)用在起伏地表，稱之為起伏多重變加密網(wǎng)格算法。該算法可以很好地刻畫(huà)地表形態(tài)并描述地下介質(zhì)的分布特征，在過(guò)渡帶區(qū)域橫向上個(gè)別點(diǎn)采用了多次內(nèi)插，在保證計(jì)算穩(wěn)定性的同時(shí)，增強(qiáng)了在不同介質(zhì)分布情況下進(jìn)行有限差分計(jì)算時(shí)網(wǎng)格剖分的靈活性。在滿足穩(wěn)定性的條件下，可以根據(jù)介質(zhì)速度選擇最優(yōu)化的網(wǎng)格進(jìn)行差分計(jì)算，有效提高了計(jì)算效率。

4 算法分析

4.1 計(jì)算精度及誤差分析

為了分析算法的計(jì)算精度，采用如圖4a所示的模型，分別用細(xì)網(wǎng)格、粗網(wǎng)格、起伏多重變加密網(wǎng)格進(jìn)行剖分，并對(duì)它們的數(shù)值模擬結(jié)果進(jìn)行對(duì)比分析。對(duì)比分析3種不同類型網(wǎng)格剖分下0.4 s時(shí)的波場(chǎng)快照（圖4 b—d）及地震記錄（圖4 e—g）可知：細(xì)網(wǎng)格與變加密網(wǎng)格的波場(chǎng)特征和地震記錄基本一致，而粗網(wǎng)格中存在嚴(yán)重的數(shù)值頻散。

網(wǎng)格變化所引起的誤差是一個(gè)需要關(guān)注的問(wèn)題，因此，進(jìn)一步對(duì)單道地震記錄進(jìn)行對(duì)比分析。在均勻介質(zhì)模型（大小2 000 m×2 000 m，速度2 000 m/s）的情況下，在同一位置處利用不同網(wǎng)格剖分的單道記錄如圖5a所示。對(duì)比均勻介質(zhì)模型細(xì)網(wǎng)格情況下和起伏多重變加密網(wǎng)格情況下的單道記錄可見(jiàn)，圖5a紅框內(nèi)起伏多重變加密網(wǎng)格所引起的反射可以定義為虛假反射。解闖等［23］對(duì)聲波方程變網(wǎng)格有限差分?jǐn)?shù)值模擬的虛假反射進(jìn)行了分析，認(rèn)為當(dāng)網(wǎng)格比不為1 時(shí)，網(wǎng)格比越大，反射率越大，虛假反射現(xiàn)象越明顯。虛假反射現(xiàn)象在網(wǎng)格發(fā)生變化時(shí)難以避免，圖5b給出了該均質(zhì)模型細(xì)網(wǎng)格和起伏多重變加密網(wǎng)格兩種情況下同一位置單道記錄的差值，可以看出，由起伏多重變加密網(wǎng)格引起的誤差控制在10-12范圍以

內(nèi)。通過(guò)對(duì)上述均勻介質(zhì)模型的誤差分析可知，網(wǎng)格間距發(fā)生變化所引起的誤差難以避免，只能將誤差控制得盡量小。以圖4a所示的起伏層狀速度模型為例，在同一位置選取3種不同網(wǎng)格情況下地震記錄中的單道進(jìn)行對(duì)比分析，由其分析結(jié)果（圖5c）可知：起伏多重變加密網(wǎng)格方法的地震記錄較粗網(wǎng)格所計(jì)算的地震記錄很大程度上壓制了數(shù)值頻散，其計(jì)算精度與細(xì)網(wǎng)格計(jì)算結(jié)果基本一致。

4.2 計(jì)算效率分析

為了分析算法的計(jì)算效率，表1給出了3種不同大小（500 m×500 m、1 000 m×1 000 m、13 000 m× 6 600 m）的模型采取不同網(wǎng)格尺寸進(jìn)行有限差分?jǐn)?shù)值模擬的計(jì)算耗時(shí)對(duì)比。其中13 000 m×6 600 m模型為黃土塬模型，根據(jù)鄂爾多斯某工區(qū)的實(shí)際地質(zhì)情況構(gòu)建，包含起伏地表、近地表薄風(fēng)化層、近地表黃土低速層、小斷層和背斜等構(gòu)造特征。計(jì)算的硬件設(shè)備參數(shù)為：Intel（R） Xeon（R） CPU E5-2620 @ 2.00 GHz，32 G運(yùn)行內(nèi)存。本文將起伏多重變加密網(wǎng)格統(tǒng)一控制為兩次過(guò)渡，模型網(wǎng)格剖分為三部分，即為起伏三重變加密網(wǎng)格，網(wǎng)格大小分別為1 m×2 m、 1 m×4 m、1 m×8 m。由表1可知：1）在模型的網(wǎng)格剖分中，隨著起伏多重變加密網(wǎng)格加密重?cái)?shù)的增加，計(jì)算耗時(shí)逐漸減少；2）對(duì)于13 000 m×6 600 m的實(shí)際黃土塬模型而言，起伏三重變加密網(wǎng)格（1 m×8 m）的計(jì)算耗時(shí)為12 708 s，僅為常規(guī)網(wǎng)格1 m×1 m的43.3%，與常規(guī)網(wǎng)格2 m×2 m的計(jì)算耗時(shí)相當(dāng)。

綜上所述，通過(guò)上述計(jì)算精度和效率對(duì)比可知：本文提出的起伏多重變加密網(wǎng)格方法，通過(guò)多次對(duì)網(wǎng)格進(jìn)行不同程度上的加密使得該算法在保證計(jì)算模擬精度的前提下，對(duì)于不同的模型在計(jì)算效率上得到一定程度的提升，且算法能很好的壓制近地表和地下復(fù)雜介質(zhì)中的數(shù)值頻散。

5 計(jì)算實(shí)例

分別采用起伏層狀模型、SEG（society of exploration geophysicists）的Canadian Foothills模型以及鄂爾多斯某實(shí)際工區(qū)的黃土塬模型驗(yàn)證算法在復(fù)雜介質(zhì)中的適應(yīng)性，通過(guò)數(shù)值模擬對(duì)比上述模型常規(guī)網(wǎng)格有限差分方法與起伏多重變加密網(wǎng)格差分方法計(jì)算所得的波場(chǎng)特征，驗(yàn)證新算法的可靠性。

5.1 起伏地表層狀介質(zhì)模型

起伏層狀模型如圖6a所示，模型空白區(qū)域?yàn)榈乇硪陨希扇∽杂傻乇磉吔鐥l件后，地震波不會(huì)向上傳播，故設(shè)置空白區(qū)域速度為零。主要測(cè)試法向虛像外推法處理起伏自由地表對(duì)起伏多重變加密網(wǎng)格方法所帶來(lái)的影響，通過(guò)測(cè)試可知，得到的細(xì)網(wǎng)格（圖6b）和起伏多重變加密網(wǎng)格（圖6c）的地震記錄基本一致。對(duì)采用起伏多重變加密網(wǎng)格計(jì)算的波場(chǎng)結(jié)構(gòu)進(jìn)行分析，采用法向虛像外推法處理的自由邊界未產(chǎn)生散射現(xiàn)象，由層狀介質(zhì)引起的反射波明顯，反射波經(jīng)自由地表多次反射產(chǎn)生的多次波發(fā)育，邊界吸收效果良好。

5.2 Canadian Foothills模型

圖7a為SEG的Canadian Foothills模型。該模型為典型的起伏地表模型，主要測(cè)試起伏多重變加密網(wǎng)格法在起伏地表地下介質(zhì)速度突變時(shí)是否會(huì)造成數(shù)值不穩(wěn)定現(xiàn)象。通過(guò)測(cè)試可知，得到的細(xì)網(wǎng)格（圖 7b）和起伏多重變加密網(wǎng)格（圖 7c）的地震記錄基本一致。對(duì)采用起伏多重變加密網(wǎng)格計(jì)算的波場(chǎng)結(jié)構(gòu)進(jìn)行分析，可以明顯看出地形起伏劇烈造成的初至嚴(yán)重扭曲和能量分配不均勻，并且對(duì)速度發(fā)生變化區(qū)域所引起的反射波顯示清晰，邊界吸收效果顯著。

5.3 實(shí)際復(fù)雜地表模型

圖8a為根據(jù)鄂爾多斯某實(shí)際工區(qū)的地質(zhì)條件建立的實(shí)際黃土塬區(qū)地球物理模型。該模型地形劇烈起伏，近地表黃土層含起伏分層的風(fēng)化層、低速層、降速層等復(fù)雜近地表介質(zhì)［24-25］，模型大小為13 000 m×

6 600 m，計(jì)算時(shí)分別采用1 m×1 m網(wǎng)格和1 m×1 m多重過(guò)渡到1 m×8 m的起伏多重變加密網(wǎng)格，圖8b、c分別為二者的數(shù)值模擬結(jié)果。對(duì)黃土塬地區(qū)的波場(chǎng)特征進(jìn)行分析：首先，從圖8b、c可以看出，由于表層巨厚黃土低速層具有很強(qiáng)的能量束縛作用，它會(huì)使得大部分能量聚集在近地表的巨厚低速層內(nèi)，進(jìn)而導(dǎo)致近道數(shù)據(jù)能量強(qiáng)信噪比極低、遠(yuǎn)道數(shù)據(jù)能量很弱的特點(diǎn)。其次，由于自由地表起伏劇烈，其初至波形態(tài)扭曲嚴(yán)重。最后，對(duì)比常規(guī)1 m×1 m網(wǎng)格與起伏多重變加密網(wǎng)格的模擬結(jié)果可知，兩種算法皆能很好地適應(yīng)該復(fù)雜黃土塬區(qū)模型，且二者數(shù)值模擬的結(jié)果基本一致。但是，本文提出的多重變加密網(wǎng)格方法計(jì)算耗時(shí)12 708 s，僅為常規(guī)1 m×1 m網(wǎng)格情況下耗時(shí)29 376 s的43.3%。

6 結(jié)論

為了適應(yīng)近地表地形起伏和復(fù)雜近地表介質(zhì)分布特征，且兼顧波動(dòng)方程數(shù)值模擬的精度和效率，本文提出一種起伏多重變加密網(wǎng)格有限差分法，對(duì)所提出方法的實(shí)現(xiàn)原理以及實(shí)現(xiàn)步驟等內(nèi)容展開(kāi)研究并給予實(shí)際算例進(jìn)行分析。基于這些研究工作，可以得出如下結(jié)論：

1）在起伏自由地表上采取法向虛像外推法，可有效規(guī)避差分計(jì)算時(shí)近地表處的降階處理，地表成階梯狀所產(chǎn)生的地表散射以及繞射等現(xiàn)象得到有效抑制。

2）起伏多重變加密網(wǎng)格可以更精細(xì)地刻畫(huà)起伏地表形態(tài)和地下介質(zhì)分布特征。算法通過(guò)縱向多點(diǎn)插值方式很好地處理了起伏變加密網(wǎng)格存在的一些常規(guī)差分無(wú)法處理的問(wèn)題，進(jìn)而在保證算法穩(wěn)定性的同時(shí)，增加了算法的靈活性。

3）起伏多重變加密網(wǎng)格方法通過(guò)多次改變網(wǎng)格的加密程度，使其計(jì)算效率得到顯著提升的同時(shí)保證模擬精度，并有效壓制了近地表的角點(diǎn)散射和數(shù)值頻散問(wèn)題。計(jì)算實(shí)例表明算法能很好地適應(yīng)實(shí)際復(fù)雜黃土塬介質(zhì)模型，且可通過(guò)該模擬結(jié)果分析一部分復(fù)雜黃土塬區(qū)的地震波傳播規(guī)律和波場(chǎng)特征。

綜上所述，本文提出的算法具有兼顧計(jì)算精度和效率且能很好適應(yīng)實(shí)際復(fù)雜近地表介質(zhì)的特性，在復(fù)雜區(qū)地震數(shù)據(jù)采集、觀測(cè)系統(tǒng)設(shè)計(jì)、近地表偏移成像、速度反演等方面具有一定的應(yīng)用前景。

參考文獻(xiàn)（References）：

［1］ 胡自多，劉威，雍學(xué)善，等. 三維波動(dòng)方程時(shí)空域混合網(wǎng)格有限差分?jǐn)?shù)值模擬方法［J］.地球物理學(xué)報(bào)，2021，64（8）：2809-2828.

Hu Ziduo， Liu Wei， Yong Xueshan， et al. Mixed-Grid Finite-Difference Method for Numerical Simulation of 3D Wave Equation in the Time-Space Domain［J］. Chinese Journal of Geophysics， 2021， 64（8）： 2809-2828.

［2］ 韓復(fù)興，王若雯，孫章慶，等. 地震聲波數(shù)值模擬中人工邊界條件的差別與組合［J］. 吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2022，52（1）：261-269.

Han Fuxing， Wang Ruowen， Sun Zhangqing， et al. Difference and Combination of Artificial Boundary Conditions in Seismic Acoustic Numerical Simulation［J］. Journal of Jilin University （Earth Science Edition）， 2022， 52 （1）： 261-269.

［3］ Jastram C， Behle A. Acoustic Modeling on a Grid of Vertically Varying Spacing［J］. Geophysics Prospecting， 1992， 40（2）： 157-169.

［4］ Jastram C， Tessmer E. Elastic Modelling on a Grid with Vertically Varying Spacing［J］. Geophysical Prospecting， 2010， 42（4）： 357-370.

［5］ Moczo P. Finite Difference Technique for SH Waves in 2-D Media Using Irregular Grid： Application to the Seismic Response Problem［J］. Geophysical Journal International， 1989， 99（2）： 321-329.

［6］ Wang Y， Xu J， Schuster G T. Viscoelastic Wave Simulation in Basins by a Variable-Grid Finite-Difference Method［J］. Bulletin of the Seismological Society of America， 2001， 91（6）： 1741-1749.

［7］ 張劍鋒. 彈性波數(shù)值模擬的非規(guī)則網(wǎng)格差分法［J］.地球物理學(xué)報(bào)，1998，41（增刊1）：357-366.

Zhang Jianfeng. Non-Orthogonal Grid Finite-Difference Method for Numerical Simulation of Elastic Wave Propagation［J］. Chinese Journal of Geophysics， 1998， 41（Sup. 1）： 357-366.

［8］ Pitarka A. 3D Elastic Finite-Difference Modeling of Seismic Motion Using Staggered Grids with Nonuniform Spacing［J］. Bulletin of Seismological Society of America， 1999， 89（1）： 54-68.

［9］ 李勝軍，孫成禹，倪長(zhǎng)寬，等. 聲波方程有限差分?jǐn)?shù)值模擬的變網(wǎng)格步長(zhǎng)算法［J］.工程地球物理學(xué)報(bào)，2007，4（3）：207-212.

Li Shengjun， Sun Chengyu， Ni Changkuan， et al. Acoustic Equation Numerical Modeling on a Grid of Varying Spacing［J］. Chinese Journal of Engineering Geophysics， 2007， 4（3）： 207-212.

［10］ 黃超，董良國(guó). 可變網(wǎng)格與局部時(shí)間步長(zhǎng)的高階差分地震波數(shù)值模擬［J］.地球物理學(xué)報(bào)，2009，52（1）：176-186.

Huang Chao， Dong Liangguo. High-Order Finite-Difference Method in Seismic Wave Simulation with Variable Grids and Local Time-Steps［J］. Chinese Journal of Geophysics， 2009， 52（1）： 176-186.

［11］ 黃超，董良國(guó). 可變網(wǎng)格與局部時(shí)間步長(zhǎng)的交錯(cuò)網(wǎng)格高階差分彈性波模擬［J］.地球物理學(xué)報(bào)，2009，52（11）：2870-2878.

Huang Chao， Dong Liangguo. Staggered-Grid High-Order Finite-Difference Method in Elastic Wave Simulation with Variable Grids and Local Time-Steps［J］. Chinese Journal of Geophysics， 2009， 52（11）： 2870-2878.

［12］ 孫成禹，丁玉才. 波動(dòng)方程有限差分雙變網(wǎng)格算法的精度分析［J］.石油地球物理勘探，2012，47（4）：545-551.

Sun Chengyu， Ding Yucai. Accuracy Analysis of Wave Equation Finite Difference with Dual-Variable Grid Algorithm［J］. Oil Geophysical Prospecting， 2012， 47（4）： 545-551.

［13］ Li H， Zhang W， Zhang Z， et al. Elastic Wave Finite-Difference Simulation Using Discontinuous Curvilinear Grid with Non-Uniform Time Step： Two-Dimensional Case［J］. Geophysics Journal International， 2015， 202（1）： 102-118.

［14］ Fan N， Zhao L F， Gao Y J， et al. A Discontinuous Collocated-Grid Implementation for High-Order Finite-Difference Modeling［J］. Geophysics， 2015， 80（4）： T175-T181.

［15］ 李世中，孫成禹，彭鵬鵬. 可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法地震波正演模擬［J］.石油物探，2018，57（3）：378-388.

Li Shizhong， Sun Chengyu， Peng Pengpeng. Seismic Wave Field Forward Modeling of Variable Staggered Grid Optimized Difference Coefficient Method［J］. Geophysical Prospecting for Petroleum， 2018， 57（3）： 378-388.

［16］ 姜占東，范彩偉，黎孝璋，等. 基于變網(wǎng)格高階有限差分的鬼波數(shù)值模擬研究［J］.地球物理學(xué)進(jìn)展，2021，36（1）：365-373.

Jiang Zhandong， Fan Caiwei， Li Xiaozhang， et al. Numerical Simulation of Marine Ghost Wave Based on High Order Finite Difference Method with Variable Grids［J］. Progress in Geophysics， 2021， 36（1）： 365-373.

［17］ Zhang Z G， Zhang W. Stable Discontinuous Grid Implementation for Collocated-Grid Finite-Difference Seismic Wave Modelling［J］. Geophysics Journal International， 2013， 192（3）： 1179-1188.

［18］ Pasalic D， Mcgarry R， Corp A. Convolutional Perfectly Matched Layer for Isotropic and Anisotropic Acoustic Wave Equations［C］//SEG Technical Program Expanded 2010. Denver： SEG， 2010： 2925-2929.

［19］ Komatitsch D， Martin R. An Unsplit Convolutional Perfectly Matched Layer Improved at Grazing Incidence for the Seismic Wave Equation［J］. Geophysics， 2007， 72（5）： SM155-SM167.

［20］ Roden J A， Gedney S D. Convolution PML （CPML）： An Efficient FDTD Implementation of the CFS–PML for Arbitrary Media［J］. John Wiley amp; Sons， Inc， 2000， 27（5）： 334-339.

［21］ 孫林潔，印興耀. 基于PML邊界條件的高倍可變網(wǎng)格有限差分?jǐn)?shù)值模擬方法［J］.地球物理學(xué)報(bào)，2011，54（6）：1614-1623.

Sun Linjie， Yin Xingyao. A Finite-Difference Scheme Based on PML Boundary Condition with High Power Grid Step Variation［J］. Chinese Journal of Geophysics， 2011， 54（6）： 1614-1623.

［22］ Zhang D L， Schuster G T， Ge Z. Multi-Source Least-Squares Reverse Time Migration with Topography［C］//SEG Technical Program Expanded Abstracts 2013. Houston： SEG， 2013： 3736-3740.

［23］ 解闖，宋鵬，譚軍，等. 聲波方程變網(wǎng)格有限差分正演模擬的虛假反射分析［J］.地球物理學(xué)進(jìn)展，2019， 34（2）： 639-648.

Xie Chuang， Song Peng， Tan Jun， et al. Analysis of Spurious Reflections of Variable Grid Finite Difference Forward Modeling Based on Acoustic Wave Equation［J］. Progress in Geophysics， 2019， 34（2）： 639-648.

［24］ 張志立，韓復(fù)興，孫文艷，等. 利用正演模擬實(shí)現(xiàn)面波衰減［J］. 吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2021，51（6）：1890-1896.

Zhang Zhili， Han Fuxing， Sun Wenyan， et al. Surface Wave Attenuation Method Using Forward Modeling［J］. Journal of Jilin University （Earth Science Edition）， 2021， 51 （6）： 1890-1896.

［25］ 邵廣周， 李遠(yuǎn)林， 岳亮. 主動(dòng)源與被動(dòng)源面波聯(lián)合勘探在黃土覆蓋區(qū)三維成像中的應(yīng)用［J］. 物探與化探， 2022， 46（4）： 897-903.Shao Guangzhou， Li Yuanlin， Yue Liang. Joint Application of Active and Passive Surface Wave in 3D Imaging of Loess Covered Area［J］. Geophysical and Geochemical Exploration， 2022， 46（4）： 897-903.