孟凡順，張 亮，李景巖，李洋森

（1.中國海洋大學 海洋地球科學學院，山東 青島 266100；2.中國海洋大學 “海底科學與探測技術(shù)”教育部重點實驗室，山東 青島 266100；3.中海石油（中國）有限公司 湛江分公司，廣東 湛江 524057）

準P波方程緊致交錯網(wǎng)格井間地震波場模擬及邊界條件

孟凡順1，2，張 亮3，李景巖1，2，李洋森1，2

（1.中國海洋大學 海洋地球科學學院，山東 青島 266100；2.中國海洋大學 “海底科學與探測技術(shù)”教育部重點實驗室，山東 青島 266100；3.中海石油（中國）有限公司 湛江分公司，廣東 湛江 524057）

研究井間地震波場的形成過程以及波場的傳播機理、規(guī)律，對于指導實際井間地震勘探有著重要的意義。基于具有垂直對稱軸的橫向各向同性（VTI）介質(zhì)中的一階準P波方程，應用具有無條件穩(wěn)定性質(zhì)的緊致交錯網(wǎng)格隱式差分格式求解該方程。重點研究了緊致交錯網(wǎng)格求解該方程的完全匹配層（PML）吸收邊界條件，在此基礎(chǔ)上實現(xiàn)了VTI介質(zhì)中一階準P波方程的井間地震波場的正演模擬。數(shù)值算例表明：緊致交錯網(wǎng)格能精準模擬VTI介質(zhì)中準P波的傳播過程，得到高精度的正演結(jié)果。一階準P波方程能以足夠的精度描述VTI介質(zhì)中準P波特征。完全匹配層吸收邊界能有效地解決人工邊界問題，是一種高效的邊界吸收算法。

井間地震，緊致交錯網(wǎng)格；VTI介質(zhì)；準P波方程；完全匹配層

0 前言

井間地震是將震源與檢波器都置于井中進行地震波觀測的一種新的物探方法。該方法可以獲得較高分辨率的地震信號，與地面地震互補可以大幅度地提高復雜陸相儲層的描述精度［1、2］。但是當?shù)叵碌刭|(zhì)比較復雜時，尤其是地質(zhì)呈現(xiàn)各向異性時，井間地震波場類型多，而且各種類型波之間相互干涉或疊加，形成了非常復雜的井間地震波場，這對于井間地震波場的識別和分離以及成像帶來了很大的麻煩［3］。因此應用正演模擬技術(shù)研究復雜地層中井間地震波場特征和傳播規(guī)律，具有非常重要的理論和實際意義。

目前，各向異性介質(zhì)中地震波的正演模擬，主要通過求解彈性波方程實現(xiàn)［4～6］。該方法用三分量矢量場描述地震波場，每個分量都包含縱波和橫波，因此能模擬比較豐富的波場信息。但也具有十分明顯的缺陷［7］，主要表現(xiàn)在：各向異性介質(zhì)彈性波方程包含多個分量，這就必然導致計算速度慢、效率低，并且這些彈性參數(shù)物理意義不明確，在實際生產(chǎn)中這些參數(shù)也很難得到。目前在實際生產(chǎn)中仍以縱波勘探為主，即使在各向異性介質(zhì)中，也仍然只記錄縱波信息。在目前技術(shù)的情況下，多分量理論數(shù)據(jù)的應用受到了限制。為克服以上缺陷，Alkhalifan［8］假設(shè)橫波速度為零，推導出各向異性介質(zhì)中的準P波方程。該方程對彈性波方程能有效的近似，更加有利于研究各向異性介質(zhì)中準P波的傳播規(guī)律。隨后，Hestholm［9］推導了該方程的一階速度～應力方程形式，并采用高階差分實現(xiàn)了該方程的數(shù)值模擬。

在數(shù)值模擬領(lǐng)域中，由于有限差分法計算效率高，容易解決邊界問題，編程簡單等優(yōu)點被廣泛地應用。Madariage［10］首次提出了交錯網(wǎng)格有限差分算法，極大地提高了有限差分數(shù)值模擬的精度。Lele［11］提出了高階緊致網(wǎng)格有限差分算法，緊致網(wǎng)格是一種隱式差分格式，與顯式差分格式不同的是其具有無條件穩(wěn)定的性質(zhì)，其主要優(yōu)勢在于應用較少的網(wǎng)格節(jié)點就可以得到高精度的模擬記錄，而且較容易解決邊界問題［12～14］。Nagarajan［15］應用緊致交錯網(wǎng)格有限差分成功解決了大渦問題。Bendiks［16］應用緊致交錯網(wǎng)格有限差求解了Navier－Stokes方程，結(jié)果表明緊致交錯網(wǎng)格有限差分能夠很好地壓制數(shù)值頻散。

作者在本文根據(jù)VTI介質(zhì)中的一階準P波方程，研究了該方程在緊致交錯網(wǎng)格中求解的高階有限差分算法。重點研究了緊致交錯網(wǎng)格一階準P波方程的PML邊界方程，在此基礎(chǔ)上實現(xiàn)了VTI介質(zhì)中準P波方程的正演模擬，得到了高精度的合成記錄。

1 基本原理

1.1 一階準P波方程

二維VTI介質(zhì)中準P波方程的一階速度～應力方程為［9］：

其中 Vx、Vz分別為x、z方向的質(zhì)點速度；ρ為密度；ψ、κ、ζ為中間變量。

1.2 緊致交錯網(wǎng)格時間差分近似

1.3 緊致交錯網(wǎng)格空間差分近似

設(shè)f（x）是連續(xù)函數(shù)，并且存在n階偏導數(shù)，令

2 N點2 N＋2階緊致交錯網(wǎng)格有限差分格式為［18］式（4）。

其中 Δx是空間采樣間隔，待定系數(shù)a、bn可以通過將處Taylor公式展開后解方程組求得。8點10階緊致交錯網(wǎng)格的待定系數(shù)為：a＝0.257 89、b1＝0.889 87、b2＝ 0.216 12、b3＝－0.004 7、b4＝0.000 15。例如，在x方向應力σxx的空間一階差分8點10階近似為：

其中 A和B是n×n階系數(shù)矩陣，F(xiàn)j和fj是n×1階列向量。

向前差分格式：

z方向空間一階差分8點10階近似為：

其中 Fi＝ ［Fi，1Fi，2…Fi，n－1Fi，n］；fi＝ ［fi，1fi，2…fi，n－1fi，n］是1×n維行向量，AT是A 的轉(zhuǎn)置，向前差分B＝B1，向后差分B＝B2。

1.4 緊致交錯網(wǎng)格差分格式的構(gòu)造

緊致交錯網(wǎng)格單元的空間交錯策略見下頁圖1，利用公式（5）、公式（6）及網(wǎng)格的空間分布即可得VTI介質(zhì)中一階準P波方程的緊致交錯網(wǎng)格時間二階，空間8點十階有限差分格式為：

式中 Δx、Δz為空間離散間隔；Δt為時間離散間隔；i、j為空間離散點序號；n為時間離散點序號。

1.5 PML吸收邊界條件

根據(jù)PML分裂原理［19］針對二維VTI介質(zhì)中的一階準P波方程，應用緊致交錯網(wǎng)格有限差分導出了適應該方程的吸收邊界方程，見式（8）～式（16）。

其中 式（12）和式（13）分別對應x方向和z方向，其余各式在x方向和z方向均可使用；Px、Pz分別為P在x方向和z方向上的分量；d（x）、d（z）為x與z方向的吸收因子。

作者以式（14）～式（16）為例給出其時間二階精度，空間8點十階精度的緊致交錯網(wǎng)格有限差分格式，即：

其中 di、dj為x與z方向吸收因子的離散值。

2 數(shù)值模擬

2.1 均勻VTI介質(zhì)模型

首先設(shè)計一個400×400的均勻VTI介質(zhì)模型，來檢驗緊致交錯網(wǎng)格有限差分數(shù)值模擬算法的可行性。合成計算的空間步長和時間步長分別為Δx＝ Δz＝10m，Δt＝1ms，介質(zhì)參數(shù)v ＝2 401m／s、vv＝ 2 450m／s、η ＝ 0.518、ρ ＝2 150kg／m3。震源為主頻30Hz的Ricker子波位于模型的中間，如下頁圖2所示。

采用時間二階精度，空間八點十階精度對上述模型進行模擬，得到的波場快照如圖3所示（見下頁）。波場快照中外層傳播較快的即為準P波，內(nèi)層呈菱角狀傳播較慢的波即為由計算誤差產(chǎn)生的P2波［20］。

由分析波場快照可以得出：緊致交錯網(wǎng)格隱式格式有限差分模擬的地震記錄數(shù)值頻散較小，模擬的效果比較有利于波場的識別、速度分析、偏移成像等。緊致交錯網(wǎng)格隱式格式準P波方程的PML吸收邊界，能有效地吸收人工邊界產(chǎn)生的虛假反射波，而且與內(nèi)層網(wǎng)格的高階格式相匹配，避免了吸收邊界區(qū)與內(nèi)層網(wǎng)格交界面由于差分階數(shù)的突變而產(chǎn)生的偽反射，是一種高效的邊界吸收算法。

為了驗證準P波方程對VTI介質(zhì)中縱波的近似程度，作者在本文做了如下試驗，震源點位于模型中間，接收點位于震源右下角水平和垂直距離均為1 000m處，其中彈性波方程正演模擬采用的參數(shù)和一階準P波方程相一致，可通過計算得出［21］。并且彈性波方程模擬時采用縱波源激發(fā)，彈性參數(shù)同上。

后面圖4為在相同彈性參數(shù)的情況下，分別應用一階準P波方程和彈性波方程得到的理論地震圖，分析地震圖可以得到以下認識：

（1）緊致交錯網(wǎng)格有限差分算法產(chǎn)生的數(shù)值頻散很小，主要表現(xiàn)在接收的振動曲線平滑很少出現(xiàn)劇烈抖動。

（2）一階準P波方程在動力學上對VTI介質(zhì)中的準P波具有很高的近似精度。主要表現(xiàn)在一階準P波方程和彈性波方程在準P波的起跳時間上一致，而且兩種方程得到振動曲線準P波部份重合度很高。

圖1 緊致交錯網(wǎng)格空間單元交錯示意圖Fig.1 Elementary cells of compact staggered grids

（3）彈性波方程模擬時，即使采用縱波源激發(fā)VTI介質(zhì)中也并非只產(chǎn)生qP波，還將產(chǎn)生與qP波耦合的qSV波；準P方程產(chǎn)生與彈性波方程同樣的qP波，還將產(chǎn)生由計算誤差引起的P2波。

2.2 層狀介質(zhì)模型

模型如圖5所示（見下頁），模型計算區(qū)域為400×400，空間采樣間隔Δx＝Δz＝10m，時間采樣間隔Δt＝1ms，激發(fā)井和接收井分別位于水平1 500m和2 500m處，井間距1 000m，道間距10m，共400道。記錄時間長度t＝1 500ms，炮點深度為1 200m和2 000m，介質(zhì)的分界面位于垂直深度1 500m、2 500m處，介質(zhì)“1”為各向同性介質(zhì)，介質(zhì)“2、3”為VTI介質(zhì)。介質(zhì)的相關(guān)參數(shù)見后面表1。其它模擬參數(shù)同均勻VTI介質(zhì)模型。

表1 層狀介質(zhì)模型參數(shù)Tab.1 The physical parameters of layer media

通過分析圖6、圖7井間地震層狀介質(zhì)模型的共炮點道集記錄特征可以得出：

（1）緊致交錯網(wǎng)格有限差分方法能高精度的模擬井間地震波場。各向同性介質(zhì)中井間地震直達波同向軸為雙曲線，而VTI介質(zhì)中井間地震直達波的同向軸不再是雙曲線形態(tài)。因為當?shù)貙颖憩F(xiàn)為各向同性時，地震波傳播的速度在不同方向上為恒定不變的，時距曲線方程為雙曲線方程。當?shù)貙颖憩F(xiàn)為VTI介質(zhì)時，地震波傳播的速度在不同方向上是變化的，此時方程就不是雙曲線方程了，表現(xiàn)在共炮點道集記錄上直達波就不是雙曲線形態(tài)。

（2）震源位于各向同性介質(zhì)中時，只產(chǎn)生縱波向四周傳播，波傳播到模型交界面處將產(chǎn)生透射縱波和反射縱波（見圖6標注），界面處不產(chǎn)生轉(zhuǎn)換橫波。這是一階準P波方程模擬和彈性波方程模擬井間波場時的主要區(qū)別，彈性波方程模擬時將產(chǎn)生各種轉(zhuǎn)換橫波［22］。對比結(jié)果表明，一階準P波方程能夠?qū)TI介質(zhì)中準P波進行很好的近似。

（3）震源位于VTI介質(zhì)中時，不僅產(chǎn)生縱波還產(chǎn)生實際介質(zhì)中并不存在的P2波，此波向上傳播到各向同性介質(zhì)交界面處將產(chǎn)生反射波（見圖7中標注I波）；當向下傳播到VTI介質(zhì)的交界面處將產(chǎn)生反射波（見圖7中標注II波）和透射波（見圖7中標注III波）。得知P2波只產(chǎn)生于VTI介質(zhì)中，在各向同性介質(zhì)中不會產(chǎn)生P2波，因此實際地質(zhì)條件下，我們只需將震源放在沉積相對比較穩(wěn)定可以近似認為各向同性的地層中進行激發(fā)，就可以有效地消除P2波的影響。

3 結(jié)論

應用緊致交錯網(wǎng)格有限差分對井間地質(zhì)模型進行數(shù)值模擬，結(jié)果表明：

（1）緊致交錯網(wǎng)格能準確模擬VTI介質(zhì)中準P波的傳播過程，得到的井間地震波場齊全、清晰?？朔藦椥圆▓稣菽M波場復雜、計算效率低、彈性參數(shù)不明確等缺點。

（2）完全匹配層吸收邊界能夠有效解決人工邊界問題，并且和計算區(qū)域交界面處產(chǎn)生虛假反射很小，是一種理想的吸收邊界條件。

（3）地震各向異性的研究，使地震學理論向?qū)嶋H地球介質(zhì)波動理論研究邁出了一大步。隨著對地震精度要求的不斷提高，基于各向異性的處理解釋技術(shù)也受到人們越來越多的重視。一階準P波方程能高精度地描述VTI介質(zhì)中準P波的特征，具有較高的近似精度，在縱波勘探領(lǐng)域有著廣闊的應用前景。

（4）在井間二維VTI介質(zhì)中準P波的數(shù)值模擬方面，將震源置于各向同性介質(zhì)中激發(fā)可以消除計算誤差產(chǎn)生的P2波影響，但是還需進一步研究更好的去除P2波的方法，使震源位于VTI介質(zhì)時也能去除P2波的影響。

［1］ 周建宇，何惺華，李安夏，等.羅家地區(qū)井間地震方法與效果［J］.石油地球物理勘探，2001，36（6）：745.

［2］ 孔慶豐，王延光，左建軍，等.井間地震有限角度疊加方法研究與應用［J］.石油地球物理勘探，2007，42（3）：256.

［3］ 李萬萬，裴正林.井間地震彈性波傳播特征數(shù)值模擬［J］.物探與化探，2008，32（2）：207.

［4］ ZHANG J F.VERSCHUUR D J.Elastic wave propagation in heterogeneous anisotropic media using the lumped finite－element method［J］.Geophysics，2002，67（2）：625.

［5］ 祝賀君，張偉，陳曉非.二維各向異性介質(zhì)中地震波場的高階同位網(wǎng)格有限差分模擬［J］.地球物理學報，2009，52（6）：1536.

［6］ 汪利民，徐義賢，江貴.井中和井間地震波場正演模擬［J］.石油物探，2009，48（2）：146.

［7］ 韓令賀，何兵壽.VTI介質(zhì)一階準P波方程正演模擬及邊界條件［J］.石油地球物理勘探，2010，45（6）：819.

［8］ ALKHALIFAH T.An acoustic wave equation for anisotropic media［J］.Geophysics，2000，65（4）：1239.

［9］ HESTHOLM S.Acoustic VTI modeling using high－order finite－differences［J］.Geophysics，2009，74（5）：T67.

［10］MADARIAGA R.Dynamics of an expanding circular fault［J］.BSSA，1976，66（3）：639.

［11］LELE S K.Compact finite difference schemes with spectral－like resolution［J］.J.Comput Phys，1992，103（1）：16.

［12］王書強，楊頂輝，楊寬德.彈性波方程的緊致差分方法［J］.清華大學學報：自然科學版，2002，42（8）：1128.

［13］馬廷福，金濤，曹福軍，等.三維雙調(diào)和方程的高階緊致差分格式及其多重網(wǎng)格方法［J］.蘭州理工大學學報，2010，36（3）：142.

［14］葛永斌，劉國濤.三維Helmholtz方程的高階隱式緊致差分方法［J］.工程數(shù)學學報，2010，27（5）：853.

［15］NAGARAJAN S.LELE S K，F(xiàn)ERZIGER J H.A robust high－order compact method for large eddy simulation［J］.Comput.Phys，2003，191（2）：392.

［16］BENDIKS.A staggered compact finite difference formulation for the compressible Navier－Stokes equations［J］.Comput.Phys，2005，208（2）：675.

［17］董良國，馬在田，曹景忠，等.一階彈性波方程交錯網(wǎng)格高階差分解法［J］.地球物理學報，2000，43（3）：411.

［18］DU Q Z，LIN B，HOU B.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered－grid finite difference scheme［J］.Applied Geophysical，2009，6（1）：42.

［19］COLLINO F，TSOGKA C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in an isotropic heterogeneous media［J］.Geophysics，2001，66（1）：294.

［20］何兵壽，張會星.VTI介質(zhì)中準P波方程的數(shù)值解法［J］.煤炭學報，2006，31（4）：446.

［21］TARIQ A.An acoustic wave equation for anisotropic media［J］.Geophysics，2000，65（4）：1239.

［22］吳律.層析基礎(chǔ)及其在井間地震中的應用［M］.北京：石油工業(yè)出版社，1997.

P 631.4

A

10.3969／j.issn.1001－1749.2012.05.02

1001—1749（2012）05—0510—08

2011－12－06 改回日期：2012－05－17

孟凡順（1960－），男，博士后，教授，主要從事地震波數(shù)值模擬與層析成像研究。