朱生旺1，2 曲壽利1 魏修成1 劉春園3
（1.中國石化石油勘探開發研究院，北京100038；2.中國地質大學（北京），北京100083；3.中國石油大學（北京），北京102249）
摘要 研究復雜介質中地震波傳播規律及地震響應特征需要在細小網格剖分下進行彈性波方程數值模擬計算，而小網格下的數值模擬計算將帶來巨大的計算量問題，采用變網格計算是減少計算量的有效途徑。本文給出壹種變網格差分計算的實現方法，在局部復雜介質區域采用細網格計算，其余區域采用粗網格計算，在兩種網格的過渡區通過改變差分算子和波場插值實現波傳播的過渡銜接。該實現方法的特點是：①在所有計算節點均采用交錯網格計算方式；②估算空間壹階導數時，保持差分計算節點的對稱性；③可選擇任意階計算精度。該方法的計算結果表明，粗細網格的過渡不會對地震波傳播模擬帶來影響，從而達到了既減少計算量又保證計算精度的目的。
關鍵詞 彈性波方程 數值模擬 變網格差分 頻散
Elastic Modeling on A Staggered Grid with Varying Spacing
ZHU Sheng-wang1，2，QU Shou-li1，WEI Xiu-cheng1，LIU Chun-yuan3
（1.Exploration ＆ Production Research lnstitute，SlNOPEC，Beijing100083；2.China University of Geoscience，Beijing100083；3.China University of Petroleum，Beijing102249）
Abstract It is necessary to use small grid in elastic wave equation modeling for researching the law of seismic wave propagation and feature of seismic response in complex medium.It is an effective scheme to calculate on a staggered grid with varying spacing which can reduce the enormous calculation amount caused by using the small grid calculation means.This article proposes the finite difference method in which small grid is used for local complex medium，and big grid is used for others.In transition region between two kinds of grid，transient conjugation of wave propagation is implemented by changing the difference operator and applying wave field interpolation.The trait of this way includes：（1）Keeping staggered-grid calculation mode in all operation node；（2）Maintaining The symmetry of operation node in estimating spatial first derivate；（3）Any precision order available in finite difference calculation.The result of the way indicates that the simulation of seismic wave propagation will not be affected by using two kinds of grid，and as a result，not only the calculation amount decreases，but also the high calculation accuracy remains.
Key words elastic wave equation modeling staggered-grid dispersion
在復雜油氣藏的勘探開發中，地震數值模擬（正演）發揮著越來越重要的作用。利用地震方法進行儲層預測，首先要研究儲層的地震響應特征，找出與儲層最為相關的地震屬性，進而優選儲層預測的技術方法，地震數值模擬是研究儲層地震響應特征不可缺少的手段。
對復雜儲層或地質體進行地震數值模擬需在足夠小的空間網格上進行計算，以保證地質模型的局部變化細節能夠準確地反映到數值模擬結果中來，即確保數值模擬結果能夠較為精確地反映介質的小尺度非均質性引起的波場變化。在正演計算中，為保證計算收斂，空間和時間方向的采樣率必須滿足收斂條件，由彈性波方程有限差分正演的收斂條件［1，2］知，減小空間計算網格尺度的同時，時間方向的遞推計算步長也必須大體上做相應比例的減小，即空間計算網格點數增加 n倍，意味著計算量將增加 n2倍。因此，若對整個模型采用細小網格計算所面臨的問題是計算量太大，即使是二維模型的正演計算，譬如使用1m×1m 或更小尺度的計算網格，當采用能夠反映實際地震反射深度的地質模型時，模擬多次覆蓋采集，其計算量亦十分驚人。實際中的數值模擬壹般是針對某壹特定的地質目標（或儲層），且地質目標僅占據整個模型很小的壹部分，對計算網格僅要求在該地質目標處細化，於是，采用變網格計算，即在目標地質體區域采用較小的計算網格，而在該區域以外采用較大的計算網格就可以大大地減少計算網格節點的數量，從而可在很大程度上降低計算量。
為減少計算量，提高數值模擬的計算效率，Jastram C和Tessmer E提出了單方向上變網格計算方法［3］，文獻［2］亦強調采用變網格計算技術，但沒有給出具體實現方法。采用變網格進行差分計算的關鍵是解決好網格尺度變化處的波場過渡銜接問題，保證在過渡區不會因網格尺度變化產生明顯的計算噪聲。本文基於二維彈性波方程數值模擬的交錯網格差分算法，通過引入不等間距節點情況下壹階導數的高階精度差分計算式，結合波場的高階精度插值和粗細網格過渡區的計算節點選取技巧，給出彈性波方程變網格正演計算的實現方法，應用結果表明，該方法能夠得到滿意的計算精度，在過渡區不會產生明顯的因網格尺度變化帶來的附加計算噪聲。
1 利用對稱任意間距節點估算壹階導數的高階精度差分近似式
在沿時間方向的遞推計算過程中，提高波場函數空間壹階導數的估計精度是彈性波方程數值模擬研究的關鍵內容之壹。對於固定差分網格，通過應用等距節點高階精度的差分算子可實現空間導數的高精度估算［1，2］。對於局部采用細網格的非固定差分網格情況，在粗、細網格區域內部可分別應用上述等距節點的高階精度的差分算子，而在粗、細網格過渡區，本文采用下面導出的具有任意節點間距的差分算子估算壹階導數值。
設在x方向上以點x0為中心對稱分布2N個網格節點，其x坐標分別為x0-qNΔx/2，…，x0-q1Δx/2，x0+q1Δx/2，…，x0+qNΔx/2，如圖1 所示，這裏Δx為節點間的最小間距。由2N個節點處的函數值f（x0-qNΔx/2），…，f（x0-q1Δx/2），f（x0+q1Δx/2），…，f（x0+qNΔx/2）可給出函數f（x）在點x0處的壹階導數估值。由泰勒展開有：
圖1 估算壹階導數的對稱計算節點示意圖
油氣成藏理論與勘探開發技術
油氣成藏理論與勘探開發技術
兩式相減，得
油氣成藏理論與勘探開發技術
略去（1）式中誤差項O（（qiΔx/2）2N+1），將精確的f（i）（x0）換成f（i）（x0）的估計 （x0），並寫成矩陣形式，則有
油氣成藏理論與勘探開發技術
記
油氣成藏理論與勘探開發技術
油氣成藏理論與勘探開發技術
設A-1為A 的逆，即 AA-1=I，則（A-1）TAT=I，即（A-1）T為 AT的逆，從而有 AT（A-1）T=I，右乘向量e1得
AT（A-1）Te1=e1 （3）
（3）由方程（2）知
油氣成藏理論與勘探開發技術
記向量（A-1）Te1=（c1，c2，…，cN）T，代入式（4）則有
油氣成藏理論與勘探開發技術
且由式（3）及（A-1）Te1=（c1，c2，…，cN）T知cn（n=1，2，…，N）滿足方程
油氣成藏理論與勘探開發技術
式（5）即為用對稱任意間距節點計算壹階導數的計算式，cn（n=1，2，…N）通過求解方程（6）得到。若取qn=2n-1（n=1，2，…，N），則cn（n=1，2，…，N）即為等距節點情況的壹階導數具有2N階精度差分近似的差分系數。
將f（x）做傅氏變換變換到波數域，即f（x） F（k），對f（x）壹階導數進行傅氏變換則有f（1）（x） .i2πkF（k），即f（x）可看作是對f（x）進行的空間線性濾（1）波結果，其濾波響應為i2πk。用f（x）離散點的值估算?f/?x得 （x），亦看作是對f（x）進行的空間線性濾波結果，記其濾波響應為 （k），由式（5）知 （k）= cnsin（πkqnΔx）。記
油氣成藏理論與勘探開發技術
顯然，用 （x）作為f（1）（x）的估計，其精度可由 （k）逼近i2πk的程度，即e（k）來反映。
2 壹階導數的變網格有限差分估算
采用變網格計算，以x方向的壹階導數估算為例來說明兩種網格過渡處的處理方法。設細網格節點間距為Δx，粗網格節點間距為mΔx，即粗網格節點間距是細網格節點間距的m倍。為表述方便，不妨取m=3，采用10階精度，即N=5，如圖2所示。圖中 表示粗網格節點； 表示粗網格中的加密節點，其與粗網格節點壹起組成細網格節點；“+”為求導點位置。①～⑥為粗細網格過渡區的6個壹階導數計算點，計算該6個點處的壹階導數所用的網格節點分別在圖2中標示，其中①及其左邊大網格區的各計算點按粗網格計算；⑥及其右邊的細網格區各計算點按細網格計算；②～⑤4個點需按式（5）計算。
圖2 兩種網格過渡區計算示意圖
壹般取2N階差分計算精度時，在粗、細網格分界的細網格壹側***有N-1個計算點需要采用如式（5）給出的對稱不等間距節點差分算子進行壹階導數估算，每壹個點的差分算子不同，加上粗網格和細網格內部計算點的差分算子，這樣***有N+1個不同的差分算子。N+1個不同的差分算子可事先算好存於內存數組中供每次計算調用。
圖3 10階精度（N=5）不同計算點差分算子對應的濾波響應
①～⑥分別為圖2中對應計算點的差分算子的濾波響應；⑦為理想的濾波響應
圖3 給出的是圖2中的①～⑥求導計算點的估算壹階導數差分算子的濾波響應。對於計算點⑥，由於參加計算的節點為細網格中的相鄰2N個網格節點，顯然其差分算子的濾波響應與理想的濾波響應最接近，精度最高；對於計算點①，參加計算的節點為粗網格中的相鄰2N個網格節點，其差分算子的濾波響應僅在低波數與理想的濾波響應接近，精度最低；而對於計算點②～⑤，精度介於計算點①與⑥之間，且隨著參加計算的近距離點增加，其精度越來越高。
在粗細網格過渡區壹階導數的估計精度不壹致可能會引起反射，附加計算噪聲，現對此做出分析。
用有限差分方法估算空間壹階導數的精度由差分頻散決定，有限差分算子的濾波響應只能在低波數段逼近壹階導數算子的濾波響應，在高波數段將出現嚴重的頻散。設kα為具有可接受差分頻散水平的最高波數，即在k≤kα的低波數區， （k）與i2πk足夠逼近，頻散很小，對最終正演結果的影響可忽略不計。記kα=αkN，這裏kN=1/（2Δx），0＜α＜1，α值由差分算子決定，如：對均勻網格時的10 階精度的差分算子，可認為α約為0.6。又設模擬地震波場的最大波數為kmax，它由子波最高頻率fmax與最低速度vmin決定，即kmax=fmax/vmin。為了使有限差分數值模擬計算結果不受差分頻散的顯著影響，必須選取Δx滿足α/（2Δx）=kα＞kmax=fmax/vmin，即
油氣成藏理論與勘探開發技術
式（7）是保證數值模擬計算精度的空間離散步長Δx應滿足的條件，只要滿足條件（7），則在波數範圍［0，kmax］內差分頻散誤差可忽略。
采用變網格計算，粗網格尺度也應該滿足式（7），即
油氣成藏理論與勘探開發技術
否則波場在粗網格區傳播就會產生強的頻散噪聲。式（8）成立，則條件（7）必然滿足，由前面對兩種網格過渡區的差分算子精度的分析結論知，此時過渡區的差分算子均滿足精度要求，即kα＞kmax。既然所有差分算子都能夠滿足不產生明顯的差分頻散噪聲的要求，差分算子變化也將不會引起明顯的計算誤差。
因此，對於兩種網格過渡區波場傳播的銜接問題，上述分析的結論是只要粗網格計算滿足精度要求，即不產生明顯的差分頻散，則按本文的計算方法，在兩種網格的過渡區不會出現明顯的反射噪聲，波場傳播的模擬精度基本不受計算網格變化的影響。
3 彈性波方程變網格有限差分正演模擬的實現
3.1 實現方法
油氣成藏理論與勘探開發技術
以二維彈性波方程數值模擬計算來說明變網格有限差分法的實現。利用二維彈性波方程進行有限差分正演計算，在交錯網格中的遞推算式為
油氣成藏理論與勘探開發技術
上式中， 和 分別表示 （x，z，t） 和 ?z（x，z，t） 的有限差分估計，余類推。在式（10）中，ρ，c11，c13，c33及c55都是空變的，為書寫簡潔，略去了其空間下標。計算中對邊界的處理是，在頂邊界采用自由邊界條件，在其他數值邊界則采用PML吸收邊界條件［4，6］。
圖4 粗、細兩種網格中計算節點的布排
在壹定的網格尺度下，交錯網格計算具有精度高的優點。當采用變網格，即在模型局部用細網格剖分進行計算時，關鍵是要合理布置計算節點，以保證在所有計算節點均能采用交錯網格計算方式。圖4 以粗細網格尺度比等於3 為例，給出本文采用的計算節點布排方式，不難看出，要保持交錯網格的計算方式，x和z方向粗細網格尺度比 mx和 mz都必須為奇數。
對波場遞推計算中空間壹階導數的有限差分估算，根據節點所處位置，分為 3種情況分別處理（為簡化圖示，不妨設差分算子長度N=2，即4階精度）：
（1）對細網格區域之外（圖4實線框之外）的計算節點采用粗網格進行計算；
（2）對細網格區域內部（圖4虛線框內）的計算節點采用細網格進行計算；
（3）對細網格邊界區域（圖4虛線框與實線框之間）的計算節點按圖2描述的方法進行計算，但在計算前對參與計算的某些點（圖4 中空心點所在的位置）的波場值需要通過內插得到，本文采用拉格朗日插值方法，內插階數與差分算子的階數相同。
3.2 算例
首先，用壹個均勻介質模型檢驗網格變化是否對模擬波場傳播產生明顯的影響。
如圖 5 所示，模型大小為 1000m×1000m，vP=3000m/s，vS=1800m/s，ρ=2.4g/cm3，vz震源位於模型正中位置，采用兩種網格劃分進行正演計算。第壹種網格劃分如圖5（a）所示，在225～275m深度段（圖中陰影部分）用1m×1m的細網格，其余部分用5m×5m的粗網格，采用該變網格進行正演計算，圖5（b）和（c）分別是160ms時刻的x分量和z分量波場切片；第二種網格劃分如圖5（d）所示，在725～775m深度段（圖中陰影部分）用1m×1m的細網格，其余部分用5m×5m的粗網格，采用該變網格進行正演計算，圖5（e）和（f）分別是160ms時刻的x分量和z分量波場切片。比較兩組計算結果可以看到，圖（b）與（e），（c）與（f）幾乎沒有差別，即局部采用細網格對波場傳播沒有帶來明顯的影響。為了看得更清楚，將沿過震源的垂直線上（圖5（d）中虛線所示）的波傳播情況顯示於圖6，從向上、向下兩個方向傳播的波的壹致性可知，波通過細網格區域時沒有產生反射噪聲。
圖5 變網格計算與波場切片
其次，通過壹個含不同尺度圓形洞的介質模型正演結果觀察變網格差分計算的效果。
如圖7所示，模型大小為1500m×600m，vP=5000m/s，vS=3000m/s，ρ=2.6g/cm3，在模型中部500m深度上分布3個半徑分別為5m，10m，20m的3個圓形洞，洞間相距200m，洞內vP=1800m/s，ρ=1.2g/cm3。模擬野外地震觀測，炮間距10m，道間距5m，排列長度為1500m，中間激發，激發震源處於模型中部地表，激發子波主頻為40Hz，采用變網格計算，粗網格尺度為5m×5m，細網格（圖7上的白色虛線框內）尺度為1m×1m，***得150個炮記錄，抽出零偏移距剖面，其結果如圖8所示。圖9（a）是激發點位於模型中部的壹炮記錄；圖9（b）是與圖9（a）相對應的對整個模型用1m×1m細網格計算的炮記錄。將這兩個炮記錄進行對比可看到二者差別甚微，證實采用變網格計算基本沒有降低計算精度。但采用變網格計算卻極大地減少了計算量，就該模型來說，如果對整個模型采用1m×1m的網格進行計算，其計算量約為采用上述變網格計算的計算量的15倍。不難看出，如果對更大尺度的模型進行正演計算，采用變網格計算對減少計算量就顯得更有必要。
圖6 變網格對波場傳播影響測試
圖7 含圓洞的介質模型
4 結論
在對復雜介質進行數值模擬計算時，采用變網格計算是解決計算精度與計算量之間矛盾的有效途徑。理論分析和數值模擬結果表明，本文給出的適應變網格計算的彈性波方程正演模擬實現方法保持了交錯網格高階差分的精度高特點，在粗細網格過渡區不產生人為反射噪聲，是壹種在保證計算精度前提下減少計算量的有效實用方法。
圖8 圓洞模型的零偏移距剖面
圖9 變網格與固定小網格計算的炮記錄對比
參考文獻
［1］董良國.壹階彈性波方程交錯網格高階差分解.地球物理學報，2000，43（3）：411～419.
［2］牟永光，裴正林.三維復雜介質地震數值模擬.北京：石油工業出版社，2005.
［3］Jastram C，Tessmer E.Elastic modeling on a grid with vertically varying spacing.Geophys.Prosp.，1994，42：357～370.
［4］Francis Collino，Chrysoula Tsogka.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in an isotopic heterogeneous media.Geophysics，2001，66（1）：294～307
［5］Cunha C A.Elastic modeling in discontinuous media.Geophysics，1993，58（12）：1840～1851.
［6］Rune Mittet.Free-surface boundary conditions for elastic staggered-grid modeling schemes.Geophysics，2002，67（5）：1616～1623.