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

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

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

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

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

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

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

0 前言

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

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

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

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

1 基本原理

1.1 一階準(zhǔn)P波方程

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

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

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

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

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

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

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

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

向前差分格式：

z方向空間一階差分8點(diǎn)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 緊致交錯(cuò)網(wǎng)格差分格式的構(gòu)造

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

式中 Δx、Δz為空間離散間隔；Δt為時(shí)間離散間隔；i、j為空間離散點(diǎn)序號(hào)；n為時(shí)間離散點(diǎn)序號(hào)。

1.5 PML吸收邊界條件

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

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

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

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

2 數(shù)值模擬

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

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

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

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

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

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

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

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

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

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

2.2 層狀介質(zhì)模型

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

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

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

［3］ 李萬(wàn)萬(wàn)，裴正林.井間地震彈性波傳播特征數(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ì)中地震波場(chǎng)的高階同位網(wǎng)格有限差分模擬［J］.地球物理學(xué)報(bào)，2009，52（6）：1536.

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

［7］ 韓令賀，何兵壽.VTI介質(zhì)一階準(zhǔn)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］王書(shū)強(qiáng)，楊頂輝，楊寬德.彈性波方程的緊致差分方法［J］.清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2002，42（8）：1128.

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

［14］葛永斌，劉國(guó)濤.三維Helmholtz方程的高階隱式緊致差分方法［J］.工程數(shù)學(xué)學(xué)報(bào)，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］董良國(guó)，馬在田，曹景忠，等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法［J］.地球物理學(xué)報(bào)，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］何兵壽，張會(huì)星.VTI介質(zhì)中準(zhǔn)P波方程的數(shù)值解法［J］.煤炭學(xué)報(bào)，2006，31（4）：446.

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

［22］吳律.層析基礎(chǔ)及其在井間地震中的應(yīng)用［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ù)值模擬與層析成像研究。