優化交錯網格差分算子的井間地震數值模擬及其頻散分析

章 陽，孟凡順

（中國海洋大學 海底科學與探測技術教育部重點實驗室，青島 266100）

優化交錯網格差分算子的井間地震數值模擬及其頻散分析

章 陽，孟凡順

（中國海洋大學 海底科學與探測技術教育部重點實驗室，青島 266100）

在井間地震有限差分數值模擬中，用離散化的高階差分方程近似連續導數的波動方程時，不可避免地會產生數值頻散，而數值頻散程度則直接影響到地震波數值模擬精度，因此為了得到清晰準確的地震波場記錄，必須盡可能地壓制數值頻散。這里在一階速度應力彈性波方程的基礎上，利用兩個約束條件構造拉格朗日函數獲取優化差分系數，與泰勒展開差分系數下的交錯網格高階差分模擬結果比較，發現改進的優化交錯網格差分算子的高階差分數值模擬能更有效地壓制數值頻散，進一步提高交錯網格高階差分數值模擬的精度，為高精度井間地震數據的波場成像、縱橫波聯合解釋等提供可靠依據。

有限差分；數值頻散；交錯網格；約束條件；優化差分算子

0 引言

地震波數值模擬技術在復雜油氣藏地震勘探開發中發揮著舉足輕重的作用，可以通過地震波數值模擬方法研究地震波在復雜介質中的傳播規律和響應特征，進而對實際地震資料的解釋以及地震波的偏移反演等處理提供可靠依據［1］。井間地震與常規的地震勘探不同，它是將震源和檢波器都放在井中，在一口井激發，相鄰井接收。井間地震正演模擬是根據給定的數學模型、震源項、地下幾何界面和物性參數研究波動方程的運動學和動力學特征。

有限差分法是目前井間地震數值模擬技術中應用最廣泛的方法，其原理是將描述介質的微分方程，用差商近似微商，進而轉化為有限差分方程進行近似求解。這種近似導致在有限差分數值模擬中不可避免地會產生數值頻散問題［2］。用離散化的高階差分方程去近似連續導數的波動方程，使相速度變成空間離散間隔和時間步長的函數，當每一波長內空間采樣太少（如粗網格情況）時，便會產生數值頻散，這是差分方程所固有的本質特征［3］。數值頻散現象即計算精度問題，對于波動方程及其差分格式已經確定的情況下，頻散現象隨著空間采樣間隔的減小而減小。Madariaga［4］首次提出交錯網格有限差分法并用于一階速度－應力彈性波方程中，數值模擬時不需要對彈性參數空間微分，減少計算量的同時有效地提高數值模擬的精度，數值頻散得到有效地改善。Dablain［5］利用Taylor級數給出聲波方程高階差分格式，成功消除空間差分頻散對數值模擬的影響，極大提高有限差分數值模擬的精度，同時還分析了時間和空間離散的頻散特征。Virieux［6－7］利用交錯網格有限差分法成功實現了各向同性介質中的SH波和P－SV波的模擬，計算精度為O（Δt2＋Δx2），同時對交錯網格與常規網格進行比較分析發現交錯網格在不增加計算量的前提下精度提高了4倍。Alford和Dablain對有限差分數值模擬的影響因素（網格大小和地震波傳播方向）進行分析，網格間距控制地震波數值模擬的精度。董良國［8］從理論上證明并分析了影響數值頻散的主要因素，同時比較常規網格與交錯網格有限差分算法產生的數值頻散。吳國枕［9］從空間差分算子時間差分算子兩方面詳細分析了有限差分數值頻散，并給出壓制數值頻散的方法即通量傳輸校正方法（FCT），取得較好的模擬效果。

本研究在前人研究成果的基礎上，對于交錯網格高階有限差分的差分系數，采取最優化原理，通過構造拉格朗日函數求取條件極值獲得優化的差分系數，比較原始差分系數和優化差分系數下的交錯網格高階有限差分模擬結果，表明優化的有限差分算法在計算效率相同的情況下能在一定程度上壓制數值頻散，提高模擬精度。

1 一階速度－應力彈性波方程有限差分數值模擬

在本研究中，主要以一階速度－應力彈性波方程為基礎。二維TI介質中，一階速度－應力彈性波方程（假設體力為零）［10］為：

其中：vi為速度分量；σxx、σzz分別為x和z方向的正應力；σxz為剪應力；ρ為介質密度；Cij為介質彈性參數。在各向同性的TI介質中，

在彈性波高階有限差分數值模擬中，實際介質是無限大的，而計算時卻是在一個有限的區域，導致在數值模擬時會引入人工邊界反射問題。本研究采取完全匹配層吸收邊界條件（PML）實現對彈性波方程數值模擬的邊界處理［13］。PML的核心思想是在計算區域外面構造有限厚度的吸收層，吸收層吸收或衰減向外傳播的波，同時在計算模型區域與吸收層的鏈接處是透明的，使得產生最小可能的虛假反射來解決模擬時的人工邊界問題。以σxx為例，根據完全匹配層的分裂思想，則時間2階，空間10階精度的交錯網格高階差分格式為：

在高階有限差分算法中，用離散化的高階差分方程近似連續導數的波動方程，必定不可避免會產生一定的誤差。隨著傳播時間的推移，下一時間的離散點處必定積累上一時間的網格離散導致的計算誤差，從而造成數值模擬結果的不穩定性。誤差大小又與所選擇的模型參數有關，參數選擇不合理，必然導致數值結果精度低，網格頻散嚴重等。因此對有限差分穩定性的分析是衡量數值模擬時參數設置是否合理的問題，本研究采用董良國［14］提出的一階速度－應力彈性波方程交錯網格高階差分穩定性分析條件進行數值模擬。

2 優化有限差分算子

在地震波數值模擬中，目前使用最普遍的是交錯網格高階有限差分法，可以有效地降低數值頻散，提高模擬精度。交錯網格高階差分算子是通過泰勒展開求取的，本研究換另外一種思路，應用拉格朗日乘數法，構造拉格朗日函數，運用條件極值求得新的交錯網格差分算子。

有限差分法數值模擬是用離散的差分算子近似替代連續的微分算子，因而會產生誤差，造成數值頻散。當離散差分算子無限接近連續微分算子時，產生的誤差就會很小，提高模擬精度，降低數值頻散?；谶@一原理，令表示在xi處的一階偏導數差分算子表示在xi處的一階偏導數微分算子，當接近“1”時，差分算子代替微分算子誤差小，近似精度高。

定義函數

根據前人的研究，應用拉格朗日乘數法，引入兩個約束條件：

條件極值為：

根據數值實驗選取l0＝0．4、β＝－12便可求得一組優化差分系數：a1＝1．250561、a2＝－0．119 811、a3＝0．031 364、a4＝－0．009 179、a5＝0．001 812。

3 數值模擬算例

3．1 各向同性單層均勻介質模型

設計一個各向同性的單層均勻介質模型，模型大小為3 000 m×3 000 m，空間網格采樣步長Δx＝Δz＝10 m，時間采樣步長Δt＝1 ms，縱波速度vp＝2 000 m／s，泊松比為0．25，密度ρ＝2 000 kg／m3，震源子波選用主頻為30 Hz的雷克子波，位于模型中央，激發井和接收井位置分別位于水平1 500 m和2 000 m處，等能量激發震源。模型參數滿足交錯網格高階差分數值模擬的穩定性條件，采用完全匹配層的吸收邊界（PML），完全匹配層的層數為δ＝50，理論反射系數R＝1e－6。數值模擬的精度為O（Δt2，Δx10），比較兩種不同差分系數的波場快照，其中（0，1 500 m）是泰勒展開差分系數的波場快照，（1 500 m，3 000 m）是優化差分系數的波長快照。圖1是t＝850 ms的混合波場的增益波場快照。

圖1 t＝850 ms時刻混合波場的增益波場快照Fig．1 The gain wave field snapshot of mixed wave field，t＝850 ms

分析圖1，當波長快照增益倍數增大時，常規的交錯網格高階有限差分數值模擬仍會產生較為嚴重的數值頻散，影響計算精度。比較混合波場快照的上下兩部分，顯然優化差分系數的交錯網格高階差分法相比泰勒展開的常規差分系數的交錯網格高階差分法，在一定程度上能更好地壓制數值頻散。對于橫波的頻散壓制效果較為明顯，對于縱波也能壓制，但效果不明顯，通過提高增益倍數能觀察到頻散的壓制。由于改進的優化差分系數交錯網格高階差分法主要是針對差分系數，通過不同的數學方法求取新的差分系數，運算效率沒有改變。因而采用優化差分系數交錯網格高階差分法進行數值模擬能達到常規交錯網格高階差分法的運算效率，同時能夠在一定程度上壓制數值頻散，提高計算精度，使波場清晰，有利于分析地震波的波場特征。

3．2 corner－edge模型

設計一個二維corner－edge模型（圖2），模型大小為5 000 m×5 000 m，空間網格采樣步長Δx＝Δz＝10 m，時間采樣步長Δt＝1 ms，第一層縱波速度為vP＝1 500 m／s，泊松比為0．25，第二層縱波速度vP＝3 500 m／s，泊松比為0．4，密度ρ＝2 000 kg／m3，震源子波選用主頻為30 Hz的雷克子波，位于模型中央，激發井和接收井位置分別位于水平2 500 m和4 000 m處，縱波震源激發。模型參數均滿足交錯網格高階差分數值模擬的穩定性條件，完全匹配層的層數為δ＝50，理論反射系數R＝1e－6。數值模擬精度為O（Δt2，Δx10），比較兩種不同差分系數的增益波長快照。圖3（a）為常規泰勒展開

圖2 二維corner－edge模型Fig．2 Two－dimensional corner－edge model

差分系數的交錯網格高階差分法記錄的混合波場的增益波場快照，圖3（b）為優化差分系數的交錯網格高階差分法記錄的混合波場的增益波場快照。觀察圖3（a）發現，常規差分系數的交錯網格高階差分數值模擬的增益波場快照出現較為嚴重的數值頻散，影響模擬的精度。當用優化差分系數的交錯網格高階差分數值模擬時，觀察圖3（b）能一定程度地壓制頻散，但仍然不能完全消除頻散。圖中標出的位置即為壓制數值頻散較為明顯的位置，效果比較好。

3．3 Marmousi模型

圖4為Marmousi模型，模型網格大小為3 401 ×701，空間網格采樣步長Δx＝Δz＝10 m，時間采樣步長Δt＝0．5 ms，模型上層是深度為900 m的海水層，下層為模型的縱波速度場。震源子波選用主頻為30 Hz的雷克子波，炮點位于（17 000 m，10 m），檢波點放在920 m深的位置，道間距10 m，縱波震源激發，記錄長度為6．5 s。完全匹配層的層數為50，理論反射系數R＝1e－6。數值模擬的精度為O（Δt2，Δx10），圖5為兩種不同差分系數混合波場的單道增益反射波地震記錄圖。

圖3 t＝1 s時刻混合波場的增益波場快照Fig．3 The gain wave field snapshot of mixed wave field，t＝1 s

圖4 Marmousi模型（橫縱坐標均為網格點數）Fig．4 Marmousi model

交錯網格高階差分法研究復雜Marmousi模型，對比常規差分系數交錯網格高階差分法和優化差分系數交錯網格高階差分法的單道增益反射波地震記錄，數值頻散有一定地改善。

圖5 Marmousi模型混合波場單道增益地震記錄圖（去直達波后反射波地震記錄）Fig．5 The single channel gain seismograms of mixed wave field with Marmousi model

4 結論

井間地震數值模擬中，用離散化的高階差分方程近似連續導數的波動方程，不可避免地產生了數值頻散問題。數值頻散會導致波形畸變或重影，嚴重影響波場特征。交錯網格高階差分法可以提高計算精度，壓制數值頻散。本研究在交錯網格高階差分法的基礎上進行改進，利用構造拉格朗日函數的方法，根據約束條件求取條件極值得到新的交錯網格差分系數。比較優化差分系數的交錯網格高階差分法和常規差分系數的交錯網格高階差分法的數值模擬結果。通過結果表明，優化差分系數法能更有效地壓制數值頻散，提高計算精度。因此在今后的地震波數值模擬時，可以采用本文改進的優化差分系數法提高計算精度，更好地壓制頻散，分析地震波場的特征。

［1］ 牟永光，裴正林．三維復雜介質地震數值模擬［M］．北京：石油工業出版社，2005．

［2］ 蔡其新，何佩軍．有限差分數值模擬的最小頻散算法及其應用［J］．石油地球物理勘探，2003，38（3）：247－251．

［3］ 李國平，姚逢昌，石玉梅，等．有限差分法地震波數值模擬的幾個關鍵問題［J］．地球物理學進展，2011，26（2）：469－476．

［4］ MADARIAGA R．Dynamic of an expanding circular fault［J］．BSSA，1976，66（3）：639－666．

［5］ DABLAIN M A．The application of high-order differencing to scalar wave equation［J］．Geophysics，1986，51（1）：54－66．

［6］ VIRIEUX J．SH-wave propagation in heterogeneous media：Velocity-stress finite difference method［J］．Geophysics，1984，49（11）：1933－1941．

［7］ VIRIEUX J．P-SV wave propagation in heterogeneous media：Velocity-stress finite difference method［J］．Geophysics，1986，5（4）：1889－1893．

［8］ 董良國，李培明．地震波傳播數值模擬中的頻散問題［J］．天然氣工業，2004，24（6）：53－56．

［9］ 吳國枕，王華忠．波場模擬中的數值頻散分析與校正策略［J］．地球物理學進展，2005，20（1）：58－65．

［10］董良國，馬在田，等．一維彈性波方程交錯網格高階差分解法［J］，地球物理學報，2000，43（3），411－419．

［11］陸基孟，王永剛．地震勘探原理［M］．北京：中國石油大學出版社，2011．

［12］劉喜武．彈性波場論基礎［M］．青島：中國海洋大學出版社，2008．

［13］COLLINO F，TSOGKA C．Applications of the perfectly matched absorbing layer model to the linear elastic dynamic problem in anisotropic heterogeneous media［J］．Geophysics，1999，66（1）：294－307．

［14］董良國，馬在田，曹景忠．一階彈性波方程交錯網格高階差分解法穩定性研究［J］，地球物理學報，2000，43（6）：856－864．

Optimized staggered grid finite－difference numerical simulation of well seismic and analysis of numerical dispersion

ZHANG Yang，MENG Fan-shun
（Key Laboratory of Submarine Geosciences and Prospecting Techniques，Ministry of Education，Ocean University of China，Qingdao 266100，China）

In the finite－difference numerical simulation of well seismic，when using discretized higher－order differential equation approximate continuous derivative wave equation，it will inevitably lead to numerical dispersion．As numerical dispersion directly affects the precision of numerical simulation of seismic wave，so in order to get more clear and accurate records of seismic wave field，we must try best to suppress numerical dispersion．In this paper，we are based on an order velocity－stress elastic wave equation，introducing two constraint conditions and construct a Lagrange function，then obtain new differential coefficients by extremely conditions．Analyze and compare the numerical simulation of an order velocity－stress elastic wave equation with high－order staggered grid differential operators by Taylor expansion and optimized differential operators，we conclude that the numerical simulation using by optimized staggered grid differential operators can more effectively suppress the numerical dispersion and improve the accuracy of simulation．and provide a reliable basis for wave field imaging of high precision seismic data in well and interpretation between P wave and S wave．

finite-difference；numerical dispersion；staggered grid；constraint condition；optimized differential operator

O 175．8

A

10．3969／j．issn．1001－1749．2014．05．17

1001-1749（2014）05-0606-07

2014-02-18 改回日期：2014-06-06

國家基金課題（41174157）

章陽（1989-），男，碩士，主要從事地震波傳播理論的正反演方法研究、有限差分數值模擬和射線追蹤正演模擬，E-mail：zyzhang－08＠163．com。