速度—應(yīng)力聲波方程混合交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬及逆時(shí)偏移

胡自多，劉威，宋家文，曾慶才，田彥燦，韓令賀

（1. 中國(guó)石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020；2. 中國(guó)石油集團(tuán)油藏描述重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州 730020；3. 東方地球物理公司研究院，河北涿州 072750；4. 中國(guó)石油勘探開發(fā)研究院，北京 100083）

0 引言

波動(dòng)方程數(shù)值模擬是研究復(fù)雜介質(zhì)中波場(chǎng)傳播規(guī)律的重要手段［1-2］，是逆時(shí)偏移（RTM）［3-4］和全波形反演（FWI）［5-6］的重要基礎(chǔ)。相比偽譜法［7-8］和有限元法［9-10］，有限差分法［11-12］具有計(jì)算效率高、占用內(nèi)存小、算法實(shí)現(xiàn)簡(jiǎn)單靈活等優(yōu)點(diǎn)，已發(fā)展成為應(yīng)用最普遍的一種波動(dòng)方程數(shù)值模擬方法［13-15］。然而，有限差分法利用差分算子近似波動(dòng)方程中的微分算子，導(dǎo)致波場(chǎng)的傳播速度與真實(shí)速度不相等，且不同頻率的波場(chǎng)傳播速度不同，即出現(xiàn)數(shù)值頻散現(xiàn)象［16］。固有的數(shù)值頻散嚴(yán)重影響有限差分法的模擬精度［17-18］，進(jìn)而影響了RTM 和FWI 的精度［19］，因此，壓制數(shù)值頻散、提高模擬精度是有限差分法的一項(xiàng)重要研究?jī)?nèi)容。

構(gòu)建更合理的差分算子和改進(jìn)差分系數(shù)計(jì)算方法是提高有限差分法模擬精度的兩條重要途徑。Dablain［16］指出，高階差分算子能夠減小數(shù)值頻散，提高模擬精度，但時(shí)間高階差分算子會(huì)使內(nèi)存占用和計(jì)算量顯著增加。此后，學(xué)者普遍保持時(shí)間2 階差分算子不變，通過設(shè)計(jì)更合理的空間差分算子以提高模擬精度。Fornberg［20］給出了基于泰勒級(jí)數(shù)展開的規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格2M階（M表示空間差分算子中與差分中心點(diǎn)等距的坐標(biāo)軸網(wǎng)格點(diǎn)的組數(shù)）空間差分算子的差分系數(shù)解析表達(dá)式，據(jù)此可以構(gòu)建采用時(shí)間2 階差分算子、空間2M階差分算子的規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格有限差分法，本文稱之為常規(guī)規(guī)則網(wǎng)格有限差分法（C-FD）和常規(guī)交錯(cuò)網(wǎng)格有限差分法（C-SFD）。CFD 和C-SFD 利用空間域頻散關(guān)系和泰勒級(jí)數(shù)展開計(jì)算差分系數(shù)，盡管空間差分算子具有2M階差分精度，但差分離散波動(dòng)方程仍然僅具有2 階差分精度。波動(dòng)方程有限差分?jǐn)?shù)值模擬通過迭代求解差分離散波動(dòng)方程實(shí)現(xiàn)，不應(yīng)通過分開逼近時(shí)間差分算子或空間差分算子達(dá)到高階差分精度，而應(yīng)該使差分離散波動(dòng)方程達(dá)到高階差分精度，才能有效提高有限差分法的模擬精度。針對(duì)二階標(biāo)量波動(dòng)方程和一階應(yīng)力—速度聲波方程有限差分模擬，一些學(xué)者提出基于時(shí)空間域頻散關(guān)系和泰勒級(jí)數(shù)展開計(jì)算差分系數(shù)，構(gòu)建了時(shí)空域規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格有限差分法，這種時(shí)空域有限差分法使得相應(yīng)的二維和三維差分離散波動(dòng)方程分別沿8 個(gè)和48 個(gè)傳播方向達(dá)到2M階差分精度，但沿其他傳播方向仍然僅具有2 階差分精度，呈現(xiàn)明顯的數(shù)值各向異性［17，21］。除了上述基于泰勒展開的差分系數(shù)算法，基于最小二乘優(yōu)化算法最小化頻散關(guān)系誤差，相速度誤差和群速度誤差求解差分系數(shù)也被廣泛采用［22-24］，這種基于最小二乘的差分系數(shù)算法能夠有效提高中高頻成分的模擬精度，但會(huì)在一定程度上損失低頻成分的模擬精度，并且最小二乘算法通過迭代優(yōu)化差分系數(shù)，計(jì)算量較大。Liu［25-26］采用最小二乘優(yōu)化相對(duì)空間域和時(shí)空域頻散關(guān)系求解差分系數(shù)，將非線性最優(yōu)化問題轉(zhuǎn)化為線性最優(yōu)化問題，求解過程不需要迭代，有效地提高了差分系數(shù)的計(jì)算效率。

僅通過改進(jìn)差分系數(shù)算法提高模擬精度的效果有限，構(gòu)建更合理的空間差分算子是有效提高模擬精度的另一重要途徑。針對(duì)2階標(biāo)量波動(dòng)方程，Liu等［27］提出了一種菱形網(wǎng)格有限差分法，并基于時(shí)空域頻散關(guān)系和泰勒展開計(jì)算差分系數(shù)，使得差分離散波動(dòng)方程沿任意傳播方向達(dá)到2M階差分精度，有效地提高了二維標(biāo)量波動(dòng)方程的模擬精度和穩(wěn)定性，但是其空間差分算子長(zhǎng)度隨M2急劇增大，導(dǎo)致計(jì)算量巨大，計(jì)算效率低。后來(lái)，Wang 等［28］又提出組合常規(guī)十字交叉型網(wǎng)格和菱形網(wǎng)格以構(gòu)建空間差分算子，有效兼顧了計(jì)算效率和模擬精度。胡自多等［29］借鑒頻率域混合網(wǎng)格有限差分法［30-31］的構(gòu)建思路，將波動(dòng)方程中的Laplace 微分算子近似表示為直角坐標(biāo)系中坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的Laplace差分算子和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的Laplace差分算子的加權(quán)平均，建立了一種混合網(wǎng)格有限差分法，它與Wang等［28］提出的有限差分法相似，但構(gòu)建思路不同。胡自多等［32］導(dǎo)出了三維直角坐標(biāo)系中非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建Laplace差分算子的方法，進(jìn)一步構(gòu)建了三維混合網(wǎng)格有限差分法，有效提高了三維標(biāo)量波動(dòng)方程的模擬精度和穩(wěn)定性。針對(duì)速度—應(yīng)力聲波方程，Tan 等［33］提出聯(lián)合利用坐標(biāo)軸網(wǎng)格點(diǎn)和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建空間差分算子近似波動(dòng)方程中的1階空間偏導(dǎo)數(shù)，建立了一種混合交錯(cuò)網(wǎng)格有限差分法，相應(yīng)的差分離散聲波方程可以達(dá)到4階或6階差分精度，有效提高了聲波模擬精度。Ren等［34-35］進(jìn)一步發(fā)展了這種混合交錯(cuò)網(wǎng)格有限差分法，使得差分離散聲波和彈性波方程最高可以達(dá)到8階差分精度。Zhou等［36-37］對(duì)混合交錯(cuò)網(wǎng)格有限差分法做進(jìn)一步改進(jìn)，使得差分離散聲波和彈性波方程可以達(dá)到任意偶數(shù)階差分精度，并采用兩步線性優(yōu)化方法簡(jiǎn)化了基于最小二乘的優(yōu)化差分系數(shù)計(jì)算，進(jìn)一步提高了聲波和彈性波的模擬精度。然而，這類混合交錯(cuò)網(wǎng)格有限差分法不恰當(dāng)?shù)厥褂昧朔亲鴺?biāo)軸網(wǎng)格點(diǎn)的對(duì)稱性，通常將與差分中心點(diǎn)距離不相等的2組非坐標(biāo)軸網(wǎng)格點(diǎn)賦予了相同的差分系數(shù)，導(dǎo)致難以推導(dǎo)差分系數(shù)的通解。

針對(duì)速度—應(yīng)力聲波方程數(shù)值模擬，本文發(fā)展了一種改進(jìn)型的混合交錯(cuò)網(wǎng)格有限差分法（M-SFD），利用坐標(biāo)軸網(wǎng)格點(diǎn)和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建空間差分算子時(shí)，確保與差分中心點(diǎn)距離相等的1 組非坐標(biāo)軸網(wǎng)格點(diǎn)賦予相同的差分系數(shù)，與差分中心點(diǎn)距離不等的任意2 組非坐標(biāo)軸網(wǎng)格點(diǎn)賦予不同的差分系數(shù)，使得所構(gòu)建的M-SFD 更為合理，并利用時(shí)空域頻散關(guān)系和泰勒級(jí)數(shù)展開建立差分系數(shù)求解方程組，導(dǎo)出差分系數(shù)通解的解析表達(dá)式。

本文首先闡述了M-SFD 中空間差分算子的構(gòu)建方法，給出相應(yīng)的差分離散聲波方程，并導(dǎo)出差分系數(shù)通解；其次，開展頻散分析和穩(wěn)定性分析；最后，利用層狀介質(zhì)模型和中國(guó)塔里木盆地典型復(fù)雜構(gòu)造模型開展數(shù)值模擬，并將M-SFD 推廣應(yīng)用于RTM，利用塔里木復(fù)雜構(gòu)造模型模擬數(shù)據(jù)開展RTM 測(cè)試。測(cè)試表明，該方法能夠有效消除由于數(shù)值頻散造成的成像假象，從而提高深層成像精度和分辨率。

1 混合交錯(cuò)網(wǎng)格有限差分法的基本原理

1.1 差分離散聲波方程的導(dǎo)出

二維速度—應(yīng)力聲波方程可表示為

式中：P=P(x，z，t) 為壓力場(chǎng)；υx=υx(x，z，t) 和υz=υz(x，z，t)分別為質(zhì)點(diǎn)振動(dòng)速度場(chǎng)的x和z分量；κ=κ(x，z)為體積模量；ρ=ρ(x，z)為介質(zhì)密度。

交錯(cuò)網(wǎng)格有限差分法求解式（1）時(shí)，波場(chǎng)變量和彈性參數(shù)定義在交錯(cuò)的網(wǎng)格位置上，如圖1 所示。波場(chǎng)變量P(x，z，t)差分離散點(diǎn)定義在空間x和z方向的半網(wǎng)格點(diǎn)、時(shí)間t方向的整網(wǎng)格點(diǎn)上（圖1紅色圓圈），其離散表達(dá)式為為整數(shù)，h為空間采樣間隔，Δt為時(shí)間采樣間隔。

圖1 二維聲波交錯(cuò)網(wǎng)格差分格式中波場(chǎng)變量及彈性參數(shù)相對(duì)位置示意圖

與C-SFD 類似，M-SFD 采用時(shí)間2 階差分算子近似式（1）中波場(chǎng)變量關(guān)于時(shí)間變量t的1 階偏導(dǎo)數(shù)，可近似表示為

C-SFD 僅利用坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建空間2M階差分算子（圖2a）近似波場(chǎng)變量關(guān)于空間的一階偏導(dǎo)數(shù)，隨著M取值的增大，新增網(wǎng)格點(diǎn)與差分中心點(diǎn)的距離也逐漸增大，對(duì)提高模擬精度的貢獻(xiàn)逐漸減小。

圖2 交錯(cuò)網(wǎng)格有限差分法的空間差分算子示意圖

本文所提M-SFD的基本思路是：聯(lián)合利用坐標(biāo)軸網(wǎng)格點(diǎn)和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建一種混合型的空間差分算子近似波場(chǎng)變量關(guān)于空間的1 階偏導(dǎo)數(shù)。圖2b～圖2e為M-SFD（N=1、2、3、4）的空間差分算子示意圖，N表示空間差分算子中與差分中心點(diǎn)等距的非坐標(biāo)軸網(wǎng)格點(diǎn)的組數(shù)。相比C-SFD，M-SFD能有效利用距離差分中心點(diǎn)更近的非坐標(biāo)軸網(wǎng)格點(diǎn)，理論上更合理。

本文提出的M-SFD 與Tan 等［33］提出的時(shí)間和空間高階交錯(cuò)網(wǎng)格有限差分法具有一定的相似性，但是他們不恰當(dāng)?shù)厥褂昧朔亲鴺?biāo)軸網(wǎng)格點(diǎn)的對(duì)稱性，將與差分中心點(diǎn)距離不相等的兩組非坐標(biāo)軸網(wǎng)格點(diǎn)賦予了相同的差分系數(shù)，例如將圖2d 中標(biāo)記為綠色②和③的兩組非坐標(biāo)軸網(wǎng)格點(diǎn)賦予了相同的差分系數(shù)，標(biāo)記為綠色②的網(wǎng)格點(diǎn)與差分中心點(diǎn)的距離為而標(biāo)記為綠色③的網(wǎng)格點(diǎn)與差分中心點(diǎn)的距離為這種不合理的賦值導(dǎo)致差分系數(shù)的解析解求解困難。本文提出的M-SFD 給任意兩組與差分中心點(diǎn)距離不等的網(wǎng)格點(diǎn)賦予不同的差分系數(shù)，理論上更合理，并且求解差分系數(shù)的解析解更容易。

M-SFD（N=1）與Tan 等［33］提出的時(shí)間4 階、空間2M階交錯(cuò)網(wǎng)格有限差分法本質(zhì)上完全相同。下面以M-SFD（N=2）為例闡述M-SFD 的基本原理。利用圖2c 中M-SFD（N=2）的空間差分算子近似方程（1）中波場(chǎng)變量關(guān)于空間變量x和z的一階導(dǎo)數(shù)，可以表示為

式中am（m=1，2，…，M）、b1、b2均為差分系數(shù)。將式（2）和式（3）代入式（1）得到

式（4）為M-SFD（N=2）對(duì)式（1）的差分離散聲波方程，同樣還可以導(dǎo)出M-SFD（N=1、3、4）對(duì)式（1）的差分離散聲波方程。

1.2 差分系數(shù)計(jì)算

合理的差分系數(shù)算法能有效提高交錯(cuò)網(wǎng)格有限差分法的模擬精度。C-SFD 利用空間域頻散關(guān)系和泰勒級(jí)數(shù)展開計(jì)算差分系數(shù)，使得空間差分算子達(dá)到2M階差分精度，但是相應(yīng)的差分離散聲波方程僅具有2階差分精度［21］。交錯(cuò)網(wǎng)格有限差分法通過迭代求解差分離散聲波方程實(shí)現(xiàn)聲波數(shù)值模擬，為了提高模擬精度，應(yīng)該設(shè)法提高差分離散聲波方程的差分精度，而不是僅僅提高空間差分算子的差分精度。

M-SFD 基于時(shí)空域頻散關(guān)系和泰勒級(jí)數(shù)展開計(jì)算差分系數(shù)，將有助于提高差分離散聲波方程的差分精度。下面以M-SFD（N=2）為例，闡述M-SFD 的差分系數(shù)計(jì)算方法。均勻介質(zhì)中，速度—應(yīng)力聲波方程（1）具有如下形式的離散平面波解

式中：AP、Aυx和Aυz為平面波的振幅；k為波數(shù)；θ為平面波傳播方向與x軸正向的夾角；ω為角頻率。

將式（5）代入式（4）得到

消去式（6）中的AP，Aυx和Aυz，并且考慮到ω=vk和κ=ρv2，得到

式中：v為介質(zhì)中聲波的傳播速度rant條件數(shù)。

式（7）為M-SFD（N=2）給出的差分離散聲波方程的頻散關(guān)系，也稱為時(shí)空域頻散關(guān)系。泰勒展開其中的三角函數(shù)，得到

其中cj、βj和γj的表達(dá)式為

使式（8）左右兩邊k2xk2zh2的系數(shù)對(duì)應(yīng)相等，得到

使式（8）左右兩邊k2xk4zh4（或k4xk2zh4）的系數(shù)對(duì)應(yīng)相等，得到

使式（8）左右兩邊k2j+2xh2j（或k2j+2zh2j）(j=0，1，2，…，M-1) 的系數(shù)對(duì)應(yīng)相等，得到

根據(jù)式（12）可得出c0=±1，當(dāng)c0從1 變?yōu)?1，相應(yīng)的差分系數(shù)為a1，a2，…，aM；b1變?yōu)槠湎喾磾?shù)，對(duì)最終結(jié)果沒有影響。這里取c0=1，并進(jìn)行計(jì)算和推導(dǎo)可以得到

由式（13）可知c0=1，c1=r2，將它們代入式（10）和式（11）并聯(lián)立求解，可得到

將式（13）代入式（9）得到

式（15）可改寫為如下矩陣方程

式（16）是一個(gè)范德蒙矩陣方程。

聯(lián)合求解方程（14）和（16）得到

式（17）為M-SFD（N=2）的差分系數(shù)通解，利用同樣的方法，可以導(dǎo)出M-SFD（N=1、3、4）的差分系數(shù)通解，見附錄A。

1.3 差分精度分析

可以將差分離散聲波方程的頻散關(guān)系誤差函數(shù)關(guān)于空間采樣間隔h或時(shí)間采樣間隔Δt的最小冪指數(shù)定義為差分離散聲波方程的差分階數(shù)。根據(jù)M-SFD（N=2）給出的差分離散聲波方程的頻散關(guān)系式（7），定義誤差函數(shù)EM-SFD(N=2)，表達(dá)式為

結(jié)合式（8）～式（12），EM-SFD(N=2)可化簡(jiǎn)為

式（19）表明，EM-SFD(N=2)關(guān)于空間采樣間隔h的最小冪指數(shù)為6，因此M-SFD（N=2）給出的差分離散聲波方程具有6 階差分精度。式（19）給出的EM-SFD(N=2)的表達(dá)式較為復(fù)雜，還可以從差分系數(shù)求解過程直接分析M-SFD（N=2）給出的差分離散聲波方程的差分精度。根據(jù)M-SFD（N=2）的差分系數(shù)求解過程可以看出，EM-SFD(N=2)可以看成關(guān)于h的多項(xiàng)式，差分系數(shù)會(huì)使EM-SFD(N=2)中的系數(shù)均為零，EM-SFD(N=2)中系數(shù)不為零的最低次項(xiàng)為因此，M-SFD（N=2）給出的差分離散聲波方程具有6階差分精度。

同樣地， M-SFD（N=1）給出的差分離散聲波方程的頻散關(guān)系誤差函數(shù)EM-SFD(N=1)可以看成關(guān)于h的多項(xiàng)式，差分系數(shù)能夠使得EM-SFD(N=1)中的系數(shù)為零，EM-SFD(N=1)中系數(shù)不為零的最低次項(xiàng)為和給出的差分離散聲波方程具有4 階差分精度；M-SFD（N=3）給出的差分離散聲波方程的頻散關(guān)系誤差函數(shù)EM-SFD(N=3)可以看成關(guān)于h的多項(xiàng)式，差分系數(shù)會(huì)使得EM-SFD(N=3)中的系數(shù)均為零，EM-SFD(N=3)中系數(shù)不為零的最低次項(xiàng)為，因此M-SFD（N=3）給出的差分離散聲波方程具有6 階差分精度；M-SFD（N=4）給出的差分離散聲波方程的頻散關(guān)系誤差函數(shù)EM-SFD(N=4)可以看成關(guān)于h的多項(xiàng)式，差分系數(shù)會(huì)使得EM-SFD(N=4)中0，1，2，…，M-1)，的系數(shù)均為零，EM-SFD(N=4)中系數(shù)不為零的關(guān)于h的最低次項(xiàng)為和因此M-SFD（N=4）給出的差分離散聲波方程具有8 階差分精度。理論上，繼續(xù)增大N的取值，M-SFD 給出的差分離散聲波方程可以達(dá)到任意偶數(shù)階差分精度。C-SFD 給出的差分離散聲波方程僅具有2 階差分精度［21］，因此，M-SFD 能更有效地提高數(shù)值模擬精度。

表1 給出了當(dāng)差分離散聲波方程達(dá)到相同階的差分精度時(shí)，本文的M-SFD 與Tan 等［33］構(gòu)建的時(shí)間和空間高階交錯(cuò)網(wǎng)格差分格式需要的非坐標(biāo)網(wǎng)格點(diǎn)數(shù)。對(duì)比可以看出，達(dá)到相同階的差分精度，本文的M-SFD 需要的非坐標(biāo)軸網(wǎng)格點(diǎn)數(shù)更少，因而計(jì)算量更小，計(jì)算效率更高。

表1 相同階差分精度本文方法與Tan 等（2014）方法所需的非坐標(biāo)軸網(wǎng)格點(diǎn)數(shù)

2 頻散分析和穩(wěn)定性分析

2.1 頻散分析

數(shù)值頻散能直接反映交錯(cuò)網(wǎng)格有限差分法的數(shù)值模擬精度。本文引入歸一化相速度誤差函數(shù)εph(kh，θ)定量描述數(shù)值頻散的大小。根據(jù)M-SFD（N=1）給出的差分離散聲波方程的頻散關(guān)系式（7）和相速度的定義的表達(dá)式為

εph(kh，θ)=0 時(shí)，相速度與真實(shí)速度相等，無(wú)數(shù)值頻散；εph(kh，θ)>0時(shí)，相速度偏大，呈現(xiàn)時(shí)間頻散；εph(kh，θ)<0時(shí)，相速度偏小，呈現(xiàn)空間頻散。同理，可導(dǎo)出C-SFD 和M-SFD（N=2、3、4）的歸一化相速度誤差εph(kh，θ)的表達(dá)式。

圖3 給出了不同M值時(shí)C-SFD 和M-SFD 的歸一化相速度誤差εph(kh，θ)隨kh和θ變化的特征?？梢钥闯觯孩費(fèi)=2 時(shí)，C-SFD 呈現(xiàn)明顯的空間頻散，相速度誤差達(dá)到-4.0%（圖3a 左）；M增至5 時(shí)，空間頻散明顯減小，相速度誤差最大為-1.5%，但出現(xiàn)較明顯的時(shí)間頻散，相速度誤差最大為2.0%（圖3b左）；M增至8，空間頻散基本消失，時(shí)間頻散進(jìn)一步增大，相速度誤差最大為2.0%（圖3c 左）。②M=2時(shí)，M-SFD 呈現(xiàn)明顯的空間頻散，相速度誤差達(dá)到-4.0%（圖3a 右）；M增至5 時(shí)，空間頻散明顯減小，相速度誤差最大為-3.0%（圖3b 右）；M增至8 時(shí)，空間頻散進(jìn)一步減小，相速度誤差最大為-1.0%（圖3c 右）。③M取值較小時(shí)（如M=2），M-SFD 和C-SFD 均具有較嚴(yán)重的數(shù)值頻散，模擬精度均較低；M取值較大時(shí)（如M=8），M-SFD 的數(shù)值頻散明顯小于C-SFD，因此M-SFD 的模擬精度優(yōu)于C-SFD。④如果將數(shù)值頻散誤差的絕對(duì)值小于1‰定義為高精度。M的取值從2 增至8 時(shí)，C-SFD 的高精度模擬區(qū)域基本沒有增大，而M-SFD 的高精度模擬區(qū)域逐步增大。M取值大于5 時(shí)，M-SFD 的高精度模擬區(qū)域明顯大于C-SFD。⑤M的取值從2增至8時(shí)，C-SFD和M-SFD 的空間頻散均逐步減小，說(shuō)明增加空間差分算子中坐標(biāo)軸網(wǎng)格點(diǎn)數(shù)有助于減小空間頻散；M取值相同時(shí)，M-SFD 的時(shí)間頻散均小于C-SFD，說(shuō)明增加空間差分算子中非坐標(biāo)軸網(wǎng)格點(diǎn)數(shù)有助于減小時(shí)間頻散。

圖3 不同M 值時(shí)C-SFD（左）和M-SFD（N=1）（右）相速度頻散等值線圖（r=0.3）

圖4 給出了不同M及N時(shí)，M-SFD 的歸一化相速度誤差εph(kh，θ)隨kh和θ變化的情況。需要注意M-SFD（M=10；N=1、2、3、4）和M-SFD（M=20；N=1、2、3、4）兩組子圖中等值線的刻度不同，色標(biāo)代表的相速度誤差范圍也不同?？梢钥闯觯孩費(fèi)=10，N=1、2、3、4 時(shí)，M-SFD 的數(shù)值頻散特征相似，基本能將相速度誤差的絕對(duì)值控制在1‰以內(nèi)（圖4 左）。②M=20，N=1時(shí)，M-SFD 主要表現(xiàn)出時(shí)間頻散，相速度誤差最大為3.0‰（圖4a右）；N增至2 時(shí)，時(shí)間頻散顯著減?。▓D4b 右）；N增至4 時(shí)，時(shí)間頻散進(jìn)一步減?。▓D4d右）；M=20，N=4時(shí)，M-SFD 基本能將相速度誤差的絕對(duì)值控制在0.1‰以內(nèi)（圖4d右）。③對(duì)于一般的高精度模擬，要求將相速度誤差的絕對(duì)值控制在1‰以內(nèi)，建議采用M-SFD（N=1），M的取值為10 左右；如果對(duì)模擬精度要求極高，要求將相速度誤差的絕對(duì)值控制在0.1‰以內(nèi)，可以采用MSFD（N=4），M的取值約為20。因此，M-SFD 應(yīng)根據(jù)模擬精度要求，合理選擇M和N的取值以兼顧計(jì)算效率。④固定M的取值，N的取值為2 和3 時(shí)，MSFD 的相速度誤差分布規(guī)律基本一致，這是因?yàn)榇藭r(shí)M-SFD 給出的差分離散聲波方程均為6階差分精度，見表1.

圖4 M=10（左）及M=20（右）時(shí)不同N 值的M-SFD 相速度頻散等值線（r=0.3）

2.2 穩(wěn)定性分析

根據(jù)M-SFD（N=2）給出的差分離散聲波方程的頻散關(guān)系式（7）得到

式中S為穩(wěn)定因子。

式（22）為M-SFD（N=2）給出的差分離散聲波方程的穩(wěn)定性條件表達(dá)式。同樣地，可以導(dǎo)出C-SFD和M-SFD（N=1、3、4）的穩(wěn)定性條件表達(dá)式。圖5 給出了穩(wěn)定性條件約束下，r的最大取值隨M的變化曲線，稱為穩(wěn)定性曲線?？梢钥闯觯孩貱-SFD 和M-SFD（N=1、2、3，4）的穩(wěn)定性均隨M取值的增大而降低；②M取值相同時(shí)，M-SFD（N=1、2、3，4）的穩(wěn)定性強(qiáng)于C-SFD；③固定M取值，M-SFD 的穩(wěn)定性會(huì)隨N取值的增大而增強(qiáng)，因此，增加空間差分算子中非坐標(biāo)軸網(wǎng)格點(diǎn)數(shù)能增強(qiáng)穩(wěn)定性；④N的取值為2 和3 時(shí)，M-SFD 的穩(wěn)定性完全相同，這是因?yàn)镸-SFD（N=2、3）給出的差分離散聲波方程均為6 階差分精度，見表1。

圖5 C-SFD 和M-SFD（N=1、2、3、4）的穩(wěn)定性曲線

3 數(shù)值模擬實(shí)例

3.1 層狀介質(zhì)模型

圖6給出了一個(gè)包含5個(gè)反射界面的層狀介質(zhì)速度模型，模型尺寸為12 km×12 km，模型空間采樣間隔h=15 m，網(wǎng)格數(shù)為801×801。震源采用主頻為22 Hz 的Ricker 子波，位于（0.15 km，0.15 km）。CSFD（M=12）采用時(shí)間采樣間隔Δt=0.5、1.0 ms；M-SFD（M=10；N=1）采用Δt=1.5 、2.0 ms；MSFD（M=10；N=2，3，4）采用Δt=2.0 ms 進(jìn)行模擬。圖7 和圖8 分別給出了C-SFD 和M-SFD（M=10；N=1）模擬得到的3.6 s 時(shí)刻壓力場(chǎng)P和質(zhì)點(diǎn)振動(dòng)速度場(chǎng)x分量υx的波場(chǎng)快照。圖9a給出了M-SFD（M=10；N=2）模擬生成的單炮記錄。為對(duì)比方便，圖9b 和圖9c 分別為C-SFD（M=12）和M-SFD（M=10；N=1，2，3，4）采用不同時(shí)間采樣間隔模擬生成的炮集記錄中第781道不同時(shí)間區(qū)間范圍內(nèi)的波形對(duì)比。

圖6 層狀介質(zhì)速度模型

圖7 不同參數(shù)C-SFD（上）及M-SFD（下）層狀介質(zhì)3.6 s 時(shí)刻壓力場(chǎng)（P）模擬波場(chǎng)快照

圖9 M-SFD 層狀介質(zhì)模型模擬炮集記錄（a）和兩種方法不同時(shí)區(qū)單道波形對(duì)比（b）及其放大顯示（c）

表2 統(tǒng)計(jì)了C-SFD（M=12）和M-SFD（M=10；N=1，2，3，4）采用不同時(shí)間采樣間隔進(jìn)行模擬的耗時(shí)和加速比，模擬總時(shí)長(zhǎng)為9.0 s，且均采用Inter Xeon Gold 5218 CPU 處理器。

表2 C-SFD 和M-SFD 法層狀介質(zhì)聲波模擬的耗時(shí)和加速比對(duì)比

圖7和圖8中兩組波場(chǎng)快照顯示：C-SFD（M=12）采用Δt=0.5 ms進(jìn)行模擬的波場(chǎng)快照中無(wú)明顯數(shù)值頻散，當(dāng)時(shí)間采樣間隔增至Δt=1.0 ms時(shí)，波場(chǎng)快照中出現(xiàn)明顯的時(shí)間數(shù)值頻散；M-SFD（M=10；N=1，2）采用Δt=2.0 ms進(jìn)行模擬得到的波場(chǎng)快照均無(wú)明顯數(shù)值頻散。

圖9b和圖9c中的單道波形顯示：C-SFD（M=12）采用Δt=1.0 ms進(jìn)行模擬，波形畸變嚴(yán)重，采用Δt=0.5 ms 進(jìn)行模擬，波形仍存在一定程度的畸變；MSFD（M=10；N=1）采用Δt=2.0 ms 進(jìn)行模擬，波形存在一定程度的畸變，采用Δt=1.5 ms進(jìn)行模擬，波形僅存在微小畸變；M-SFD（M=10；N=2，3，4）采用Δt=2.0 ms進(jìn)行模擬，波形畸變基本可以忽略，模擬波形與參考波形基本重合。

綜合圖7～圖9 給出的數(shù)值頻散特征及表2 給出的計(jì)算效率統(tǒng)計(jì)結(jié)果，可以看出：M-SFD（M=10；N=1）采用Δt=1.5 ms比C-SFD（M=12）采用Δt=0.5 ms 模擬的數(shù)值頻散更小，精度更高，且計(jì)算效率提高了1.96 倍；M-SFD（M=10；N=2）采用Δt=2.0 ms 比M-SFD（M=10；N=1）采用Δt=1.5 ms模擬，數(shù)值頻散更小，精度更高，且計(jì)算效率更高。

仔細(xì)對(duì)比圖9c 中的單道波形會(huì)發(fā)現(xiàn)，采用相同的時(shí)間采樣間隔Δt=2.0 ms，M-SFD（M=10；N=3）與M-SFD（M=10；N=2）的模擬精度基本一致，M-SFD（M=10；N=4）相比M-SFD（M=10；N=2，3）在模擬精度方面存在十分微弱的優(yōu)勢(shì)，需要在更長(zhǎng)傳播距離和傳播時(shí)間的情況下，這種優(yōu)勢(shì)才能更好地體現(xiàn)出來(lái)。

3.2 復(fù)雜構(gòu)造模型

圖10a 為中國(guó)塔里木盆地典型復(fù)雜構(gòu)造速度模型，模型尺寸為18 km×7.875 km，模型空間采樣間隔h=15 m，網(wǎng)格數(shù)為1201×526。震源采用主頻為25 Hz 的Ricker 子波，位于（9 km，0.15 km），C-SFD（M=10）和M-SFD（M=8；N=1）分別采用時(shí)間采樣間隔Δt=1.0 ms 和Δt=1.5 ms 進(jìn)行模擬。C-SFD（M=10）和M-SFD（M=8；N=1）的空間差分算子均包含20 個(gè)網(wǎng)格點(diǎn)，單次時(shí)間迭代的計(jì)算量基本相同，模擬時(shí)長(zhǎng)相同時(shí)，計(jì)算效率與時(shí)間采樣間隔成正比。

圖10 中國(guó)塔里木盆地典型復(fù)雜模型及聲波交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬單炮記錄（壓力場(chǎng)P）

圖10b 給出了M-SFD（M=8；N=1）采用采Δt=1.5 ms 模擬生成的單炮記錄。為了對(duì)比方便，圖10c～圖10f 給出了C-SFD（M=10）和MSFD（M=8；N=1）不同時(shí)間采樣間隔模擬生成單炮的局部放大顯示。對(duì)比可以看出：C-SFD（M=10）采用Δt=1.0 ms 和Δt=1.5 ms 模擬生成的單炮記錄中均存在明顯的時(shí)間頻散；M-SFD（M=8；N=1）采用兩種Δt模擬生成的單炮記錄中均無(wú)明顯數(shù)值頻散。

復(fù)雜構(gòu)造模型模擬結(jié)果表明：計(jì)算效率基本相同時(shí)，M-SFD 比C-SFD 能更有效地壓制數(shù)值頻散，模擬精度更高；M-SFD 可采用更大的時(shí)間采樣間隔以提高計(jì)算效率，并保持更高的模擬精度。

4 混合交錯(cuò)網(wǎng)格有限差分法在逆時(shí)偏移中的應(yīng)用

將M-SFD 作為逆時(shí)偏移算法中的波場(chǎng)傳播算子，以圖10a 中的復(fù)雜構(gòu)造模型進(jìn)行逆時(shí)偏移測(cè)試。震源采用主頻為25 Hz 的Ricker 子波。C-SFD（M=15）采用非常小的時(shí)間采樣間隔Δt=0.1 ms，模擬生成150 炮無(wú)數(shù)值頻散的炮集資料作為逆時(shí)偏移的輸入道集，每炮600道接收，炮間距120 m，道間距30 m。

分別以C-SFD（M=10）和M-SFD（M=8；N=1）作為逆時(shí)偏移算法中的波場(chǎng)傳播算子，時(shí)間采樣間隔Δt=1.5 ms，采用互相關(guān)成像條件，并對(duì)成像結(jié)果做Laplace 濾波去除低頻噪聲。圖11 給出了C-SFD（M=10）和M-SFD（M=8；N=1）的逆時(shí)偏移剖面。對(duì)比可以看出：C-SFD（M=10）的偏移成像結(jié)果中，深層同相軸存在大量由于數(shù)值頻散造成的成像假象；M-SFD（M=8；N=1）的偏移成像結(jié)果中，由于數(shù)值頻散造成的成像假象消失，深層同相軸能量更強(qiáng)，分辨率更高。因此，相比C-SFD，M-SFD 作為逆時(shí)偏移的波場(chǎng)傳播算子能有效提高深層構(gòu)造的成像精度和分辨率。

圖11 中國(guó)塔里木盆地典型復(fù)雜構(gòu)造模型聲波混合交錯(cuò)網(wǎng)格有限差分RTM 剖面

5 結(jié)論

本文聯(lián)合利用坐標(biāo)軸網(wǎng)格點(diǎn)和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建空間差分算子近似波場(chǎng)的一階空間偏導(dǎo)數(shù)，建立了適用于速度—應(yīng)力聲波方程模擬的M-SFD，并基于時(shí)空域頻散關(guān)系和泰勒級(jí)數(shù)展開導(dǎo)出了差分系數(shù)通解。通過差分精度分析、頻散分析、穩(wěn)定性分析、數(shù)值模擬和逆時(shí)偏移試驗(yàn)，得到以下結(jié)論：

（1）C-SFD 給出的差分離散聲波方程僅具有二階差分精度，M-SFD（N=1、2、4）給出的差分離散聲波方程可達(dá)到4、6、8階差分精度，繼續(xù)增大N的取值，理論上可以達(dá)到任意偶數(shù)階差分精度；

（2）計(jì)算效率基本相同時(shí)，M-SFD 比C-SFD 能更有效地壓制數(shù)值頻散，模擬精度更高；M-SFD 可通過采用更大的時(shí)間采樣間隔以提高計(jì)算效率，且模擬精度仍然高于C-SFD；

（3）M-SFD 作為逆時(shí)偏移的波場(chǎng)傳播算子能夠更有效地消除數(shù)值頻散造成的成像假象，進(jìn)而提高深層的成像精度和分辨率。

附錄A M-SFD（N=1，3，4）的差分離散聲波方程及差分系數(shù)通解

與文中M-SFD（N=2）的差分系數(shù)求解過程類似，利用M-SFD（N=1）對(duì)聲波方程（式（1））進(jìn)行差分離散可以導(dǎo)出相應(yīng)差分離散聲波方程，再將離散平面波解代入差分離散聲波方程導(dǎo)出相應(yīng)的時(shí)空域頻散關(guān)系，并利用泰勒公式對(duì)頻散關(guān)系中的三角函數(shù)進(jìn)行級(jí)數(shù)展開，然后，使得方程左右兩邊的系數(shù)對(duì)應(yīng)相等以構(gòu)建關(guān)于差分系數(shù)的方程組，求解此方程組可以得到M-SFD（N=1）的差分系數(shù)a1，a2，…，aM；b1的通解為

利用泰勒公式對(duì)M-SFD（N=3）給出的差分離散聲波方程的時(shí)空域頻散關(guān)系中的三角函數(shù)進(jìn)行級(jí)數(shù)展開，然后，使得方程左右兩邊0，1，2，…，M-1)的系數(shù)對(duì)應(yīng)相等以構(gòu)建關(guān)于差分系數(shù)的方程組，求解此方程組可以得到M-SFD（N=3）的差分系數(shù)a1，a2，…，aM；b1，b2，b3的通解為

利用泰勒公式對(duì)M-SFD（N=4）給出的差分離散聲波方程的時(shí)空域頻散關(guān)系中的三角函數(shù)進(jìn)行級(jí)數(shù)展開，然后使得方程左右兩邊或的系數(shù)對(duì)應(yīng)相等以構(gòu)建關(guān)于差分系數(shù)的方程組，求解此方程組可以得到差分系數(shù)a1,a2,…,aM;b1,b2,b3,b4的通解為