速度—應力聲波方程混合交錯網格有限差分數(shù)值模擬及逆時偏移

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

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

0 引言

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

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

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

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

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

1 混合交錯網格有限差分法的基本原理

1.1 差分離散聲波方程的導出

二維速度—應力聲波方程可表示為

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

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

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

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

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

圖2 交錯網格有限差分法的空間差分算子示意圖

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

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

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

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

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

1.2 差分系數(shù)計算

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

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

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

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

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

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

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

其中cj、βj和γj的表達式為

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

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

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

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

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

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

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

式（16）是一個范德蒙矩陣方程。

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

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

1.3 差分精度分析

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

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

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

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

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

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

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

2.1 頻散分析

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

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

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

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

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

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

2.2 穩(wěn)定性分析

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

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

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

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

3 數(shù)值模擬實例

3.1 層狀介質模型

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

圖6 層狀介質速度模型

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

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

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

表2 C-SFD 和M-SFD 法層狀介質聲波模擬的耗時和加速比對比

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

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

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

仔細對比圖9c 中的單道波形會發(fā)現(xiàn)，采用相同的時間采樣間隔Δ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)勢，需要在更長傳播距離和傳播時間的情況下，這種優(yōu)勢才能更好地體現(xiàn)出來。

3.2 復雜構造模型

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

圖10 中國塔里木盆地典型復雜模型及聲波交錯網格有限差分數(shù)值模擬單炮記錄（壓力場P）

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

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

4 混合交錯網格有限差分法在逆時偏移中的應用

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

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

圖11 中國塔里木盆地典型復雜構造模型聲波混合交錯網格有限差分RTM 剖面

5 結論

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

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

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

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

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

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

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

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