章 陽，孟凡順

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

優(yōu)化交錯網格差分算子的井間地震數(shù)值模擬及其頻散分析

章 陽，孟凡順

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

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

有限差分；數(shù)值頻散；交錯網格；約束條件；優(yōu)化差分算子

0 引言

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

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

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

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

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

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

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

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

2 優(yōu)化有限差分算子

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

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

定義函數(shù)

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

條件極值為：

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

3 數(shù)值模擬算例

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的雷克子波，位于模型中央，激發(fā)井和接收井位置分別位于水平1 500 m和2 000 m處，等能量激發(fā)震源。模型參數(shù)滿足交錯網格高階差分數(shù)值模擬的穩(wěn)定性條件，采用完全匹配層的吸收邊界（PML），完全匹配層的層數(shù)為δ＝50，理論反射系數(shù)R＝1e－6。數(shù)值模擬的精度為O（Δt2，Δx10），比較兩種不同差分系數(shù)的波場快照，其中（0，1 500 m）是泰勒展開差分系數(shù)的波場快照，（1 500 m，3 000 m）是優(yōu)化差分系數(shù)的波長快照。圖1是t＝850 ms的混合波場的增益波場快照。

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

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

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的雷克子波，位于模型中央，激發(fā)井和接收井位置分別位于水平2 500 m和4 000 m處，縱波震源激發(fā)。模型參數(shù)均滿足交錯網格高階差分數(shù)值模擬的穩(wěn)定性條件，完全匹配層的層數(shù)為δ＝50，理論反射系數(shù)R＝1e－6。數(shù)值模擬精度為O（Δt2，Δx10），比較兩種不同差分系數(shù)的增益波長快照。圖3（a）為常規(guī)泰勒展開

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

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

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，縱波震源激發(fā)，記錄長度為6．5 s。完全匹配層的層數(shù)為50，理論反射系數(shù)R＝1e－6。數(shù)值模擬的精度為O（Δt2，Δx10），圖5為兩種不同差分系數(shù)混合波場的單道增益反射波地震記錄圖。

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

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

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

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

4 結論

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

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

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

［3］ 李國平，姚逢昌，石玉梅，等．有限差分法地震波數(shù)值模擬的幾個關鍵問題［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］ 董良國，李培明．地震波傳播數(shù)值模擬中的頻散問題［J］．天然氣工業(yè)，2004，24（6）：53－56．

［9］ 吳國枕，王華忠．波場模擬中的數(shù)值頻散分析與校正策略［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］董良國，馬在田，曹景忠．一階彈性波方程交錯網格高階差分解法穩(wěn)定性研究［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-），男，碩士，主要從事地震波傳播理論的正反演方法研究、有限差分數(shù)值模擬和射線追蹤正演模擬，E-mail：zyzhang－08＠163．com。