虛譜法無反射速度-應(yīng)力方程地震數(shù)值模擬

杜增利, 李春紅, 傅雨濛, 李永章

(1.西南石油大學(xué) 地球科學(xué)與技術(shù)學(xué)院，成都 610500；2.中國(guó)石油大慶鉆探工程公司 物探一公司，黑龍江 大慶 163357； 3.長(zhǎng)安大學(xué) 建筑學(xué)院，西安 710018)

?

虛譜法無反射速度-應(yīng)力方程地震數(shù)值模擬

杜增利1, 李春紅2, 傅雨濛3, 李永章2

(1.西南石油大學(xué) 地球科學(xué)與技術(shù)學(xué)院，成都 610500；2.中國(guó)石油大慶鉆探工程公司 物探一公司，黑龍江 大慶 163357； 3.長(zhǎng)安大學(xué) 建筑學(xué)院，西安 710018)

基于單程波波場(chǎng)延拓算子的零偏移距地震數(shù)值模擬方法難以適應(yīng)復(fù)雜地質(zhì)條件下地震波速度的空間劇烈變化，且其最大成像角度有不同程度的限制。本文依據(jù)地震反射的形成機(jī)理，重置模型空間的密度分布，使得模型空間內(nèi)縱波阻抗為常量以壓制層間多次波；采用虛譜法可得到高精度的地震波場(chǎng)空間導(dǎo)數(shù)，因此可有效抑制空間數(shù)值頻散，同時(shí)提高了數(shù)值模擬的垂向分辨率。數(shù)值模擬結(jié)果表明，得到的零偏移距地震記錄繞射波場(chǎng)完整，偏移歸位效果優(yōu)于頻率-波數(shù)域或頻率-空間域的單程波模擬結(jié)果。無反射一階應(yīng)力-速度方程零偏移距數(shù)值模擬可有效壓制多次反射波，在較大網(wǎng)格間距時(shí)，虛譜法仍可有效抑制數(shù)值頻散。

虛譜法；有限差分；數(shù)值模擬；無反射；零偏移距

地震波場(chǎng)的復(fù)雜性使得復(fù)雜構(gòu)造區(qū)地震、地質(zhì)解釋異常困難。隨著計(jì)算機(jī)技術(shù)的迅猛發(fā)展，雖然方位各向異性的黏彈性等炮域地震數(shù)值模擬技術(shù)日趨成熟，但零偏移距地震數(shù)值模擬仍是驗(yàn)證解釋方案正確性的有效手段之一[1-3]。目前，零偏移距地震記錄數(shù)值模擬多采用頻率-波數(shù)域或頻率-空間域的單程波波場(chǎng)延拓算子。

當(dāng)介質(zhì)速度在整個(gè)模型空間為常數(shù)時(shí)，相移法(phase shift, 簡(jiǎn)稱PS)能得到方程的精確解[4]，且算法無條件穩(wěn)定、背景噪聲低。Gazdag解決了速度的縱向變化問題，但仍不能適應(yīng)地層速度的橫向變化[5]；為了能同時(shí)解決速度橫向變化和大傾角問題，在相移法基礎(chǔ)上，Stoffa等提出了裂步傅里葉法(split step Fourier，簡(jiǎn)稱SSF)，在頻率-空間域增加了誤差校正項(xiàng)，使之能處理橫向速度變化較小的情況[6]。為使裂步傅里葉法能適應(yīng)更強(qiáng)的速度橫向變化，Kessinger使用多個(gè)參考速度，通過地震波場(chǎng)的Sinc插值得到任意速度的近似波場(chǎng)；但增加了計(jì)算量，導(dǎo)致計(jì)算效率降低[7]。為補(bǔ)償舍棄速度攝動(dòng)所引起的高階誤差，在裂步傅里葉法的基礎(chǔ)上，Ristow 等增加了有限差分補(bǔ)償項(xiàng)，即傅里葉有限差分法(finite Fourier difference，簡(jiǎn)稱FFD)[8]。劉文革等進(jìn)一步分析了傅里葉有限差分各階近似的精度問題[9]。

盡管諸多學(xué)者對(duì)單程波波場(chǎng)延拓算子進(jìn)行了各種近似處理，使之能夠處理速度的橫向變化，但對(duì)最大成像角度都有一定限制。而雙程波波動(dòng)方程既能處理速度的橫向突變，又無地層傾角限制，但在零偏移距地震波場(chǎng)模擬中會(huì)產(chǎn)生多次反射。為壓制層間多次波，Baysal等構(gòu)建了無反射雙程聲波方程[10]；單延明等利用無反射雙程聲波方程制作了自激自收時(shí)間剖面[11]。

有限差分法采用離散化的差分方程來逼近波動(dòng)方程，但受空間和時(shí)間采樣密度的限制，該方法不可避免地會(huì)產(chǎn)生數(shù)值頻散，為避免或減少此效應(yīng)，網(wǎng)格半徑一般都取很小并采用高階差分。地震波場(chǎng)空間導(dǎo)數(shù)的計(jì)算方法主要有高階有限差分法和虛譜法2種[12-14]。采用高階有限差分法得到的地震波場(chǎng)空間導(dǎo)數(shù)在邊界附近會(huì)出現(xiàn)振蕩；而虛譜法是通過傅里葉變換求取地震波場(chǎng)的空間導(dǎo)數(shù)，無需模型空間邊界外地震波場(chǎng)的參與。應(yīng)用傅里葉變換求取空間導(dǎo)數(shù)時(shí)，由于傅里葉變換的周期性，非光滑函數(shù)的相位譜在奈奎斯特頻率處的突變?cè)斐蓵r(shí)、空域信號(hào)的周期性震蕩，而在半網(wǎng)格點(diǎn)求導(dǎo)可有效克服Gibbs效應(yīng)[12,16]。本文應(yīng)用非均勻各向同性聲學(xué)介質(zhì)中的一階壓力-速度方程，采用虛譜法計(jì)算地震波場(chǎng)的空間導(dǎo)數(shù)實(shí)現(xiàn)零偏移距地震波場(chǎng)模擬。

1 聲波方程

1.1 聲波方程有限差分解

非均勻各向同性聲學(xué)介質(zhì)中的一階應(yīng)力-速度方程為[15-17]

(1)

其中:地震波場(chǎng)向量;P=(p,ux,uy,uz)T且p為應(yīng)力;ux,uy,uz分別為x,y,z方向的速度分量; A1,A2,A3為模量矩陣，且有

相應(yīng)的差分解為[14]

(2)

其中：n表示當(dāng)前時(shí)刻；n-表示前一時(shí)刻；Δt為時(shí)間采樣率；f={fx,fy,fz}T為外力向量；L表示求導(dǎo)運(yùn)算；i,j,k為空間坐標(biāo)。

1.2 無反射條件

采用雙程波動(dòng)方程進(jìn)行波場(chǎng)模擬時(shí)，地震波在通過介質(zhì)分界面時(shí)既有透射也有反射，而在零偏移距地震記錄正演過程中不希望有多次波出現(xiàn)，即希望地震波在通過介質(zhì)分界面時(shí)只有透射而無反射。

眾所周知，地震波在介質(zhì)分界面上產(chǎn)生反射的前提條件是聲阻抗存在差異，因此，通過人為設(shè)定介質(zhì)密度使模型空間內(nèi)聲阻抗為常量就可以有效壓制層間反射，即

ρi,j,k=v0/vPi,j,k

(3)

其中：v0為任意給定的速度常量。為了使得到的介質(zhì)密度符合地質(zhì)規(guī)律，v0一般取模型空間中的最大速度。

1.3 初始條件

根據(jù)爆炸反射面成像原理，在整個(gè)模型空間內(nèi)的所有網(wǎng)格節(jié)點(diǎn)上均放置震源，震源強(qiáng)度為該點(diǎn)的反射系數(shù)，并在t=0時(shí)刻同時(shí)引爆。采用零相位Ricker子波作為震源，其函數(shù)為

w(t)={1-2[πf(t-t0)]2}e[πf(t-t0)]2

(4)

1.4 虛譜法空間差分計(jì)算

利用傅里葉變換計(jì)算半網(wǎng)格點(diǎn)處(交錯(cuò)網(wǎng)格)地震波場(chǎng)向量的空間一階導(dǎo)數(shù)公式見文獻(xiàn)[12]。

虛譜法與高階(16階)差分法在常速介質(zhì)中的脈沖響應(yīng)對(duì)比如圖1所示，其中聲波速度設(shè)為3 km/s，x和z方向的網(wǎng)格邊長(zhǎng)均為10 m。結(jié)果表明，采用16階差分計(jì)算地震波場(chǎng)的空間導(dǎo)數(shù)，在Ricker子波主頻為40 Hz時(shí)在地震波場(chǎng)中數(shù)值頻散較明顯，子波主頻為35 Hz時(shí)地震波場(chǎng)中無數(shù)值頻散發(fā)生；而采用虛譜法計(jì)算地震波場(chǎng)的空間導(dǎo)數(shù)，在Ricker子波主頻為55 Hz時(shí)數(shù)值頻散才較為明顯，地震波場(chǎng)不產(chǎn)生數(shù)值頻散的子波最高主頻為50 Hz。由此可以說明，由于虛譜法可以得到高精度的地震波場(chǎng)空間導(dǎo)數(shù)，數(shù)值模擬時(shí)可以采用較高主頻的地震子波(本例可提高15 Hz)，因而可以得到更高垂向分辨率的模擬結(jié)果。

1.5 邊界條件

采用Collino等給出的一階應(yīng)力-速度彈性波數(shù)值模擬的完美匹配層吸收邊界條件(PML)[18]，其基本思想就是人為地將地震波場(chǎng)分解為垂直和平行于傳播方向的2組波，并使其中的沿界面法向方向傳播的地震波場(chǎng)在衰減帶內(nèi)快速衰減，以達(dá)到減少邊界反射的目的[17-18]。

2 數(shù)值模擬實(shí)例

2.1 SEG/EAGE模型

SEG/EAGE模型(圖2-A)主要針對(duì)大傾角地層和高速膏巖層邊界及其下伏地層成像問題而設(shè)計(jì)，其參數(shù)為：長(zhǎng)51 560 ft (1 ft=30.48 cm)，深12 000 ft，Δx=Δz=40 ft，x方向網(wǎng)格數(shù)1 290，z方向網(wǎng)格數(shù)301。

圖1 常速介質(zhì)中的脈沖響應(yīng)對(duì)比Fig.1 Comparison of snap shots in constant velocity medium(A)子波主頻為35 Hz的16階差分法結(jié)果; (B)子波主頻為40 Hz的16階差分法結(jié)果; (C)子波主頻為50 Hz的虛譜法結(jié)果; (D)子波主頻為55 Hz的虛譜法結(jié)果

單程波波動(dòng)方程零偏移距地震數(shù)值模擬采用相移加插值算法，為保證插值的準(zhǔn)確度，參考速度最多取7個(gè)，其合成記錄只有一次反射波和繞射波(圖2-B)。由于無最大傾角限制，常規(guī)的雙程波一階應(yīng)力-速度方程零偏移距數(shù)值模擬結(jié)果中繞射波場(chǎng)更加豐富(圖2-C中藍(lán)色箭頭所示)；但多次反射波較強(qiáng)(圖2-C中紅色箭頭所示)，多種波場(chǎng)的疊加效應(yīng)導(dǎo)致中深層背景噪聲較強(qiáng)(圖2-C中紅色矩形框所示)。與常規(guī)的雙程波方法結(jié)果相比，無反射條件的雙程波零偏移距數(shù)值模擬結(jié)果中多次反射得到了較好的壓制(圖2-D中紅色箭頭所示)，進(jìn)而減弱了中深層的背景噪聲(圖2-D中紅色矩形框所示)。

為驗(yàn)證地震數(shù)值模擬結(jié)果的準(zhǔn)確性，采用裂步傅里葉加插值法進(jìn)行深度偏移處理。結(jié)果表明，受單程波波場(chǎng)延拓算子最大成像傾角的限制，相移加插值法單程波正演、裂步傅里葉加插值法深度偏移結(jié)果在大傾角鹽丘邊界處成像較差(圖3-A，紅色箭頭所示)，且有偏移過量造成偏移劃弧現(xiàn)象(圖3-A，藍(lán)色箭頭、藍(lán)色方框所示)。由于常規(guī)的雙程波正演結(jié)果中多次波發(fā)育，裂步傅里葉加插值法深度偏移后多次波同樣成像(圖3-B，藍(lán)色箭頭所示)，在剖面上造成地質(zhì)假象，影響地質(zhì)人員正確認(rèn)識(shí)地下地質(zhì)特征；同時(shí)，多次波和繞射波場(chǎng)疊加導(dǎo)致深度偏移不能準(zhǔn)確歸位，中深層信噪比較低(圖3-B，藍(lán)色方框所示)。由于無傾角限制，無反射雙程波正演、裂步傅里葉加插值法深度偏移后，圖3-C中紅色箭頭所示位置的大傾角鹽丘邊界成像質(zhì)量明顯提高，且無偏移過量的劃弧現(xiàn)象(圖3-C，藍(lán)色箭頭、藍(lán)色方框所示)。

圖2 SEG/EAGE模型及零偏移距正演結(jié)果對(duì)比Fig.2 SEG/EAGE velocity model and comparison of zero offset simulation results(A)SEG/EAGE速度模型; (B)單程波數(shù)值模擬結(jié)果; (C)常規(guī)雙程波動(dòng)方程數(shù)值模擬結(jié)果; (D)無反射雙程波動(dòng)方程數(shù)值模擬結(jié)果。1 ft=30.48 cm

圖3 SEG/EAGE模型疊后深度偏移結(jié)果對(duì)比Fig.3 Comparison of post-stack depth migration in SEG/EAGE model(A)相移加插值法單程波正演、裂步傅里葉加插值法深度偏移結(jié)果; (B)常規(guī)雙程波正演、裂步傅里葉加插值法深度偏移結(jié)果; (C)無反射雙程波正演、裂步傅里葉加插值法深度偏移結(jié)果。1 ft=30.48 cm

2.2 Sigsbee模型

采用Sigsbee模型(圖4-A)測(cè)試數(shù)值模擬方法對(duì)地質(zhì)界面劇烈變化和高速膏巖層邊界的適應(yīng)性，模型參數(shù)為：長(zhǎng)16 000 ft (1 ft=30.48 cm)，深6 000 ft，Δx=Δz=10 ft，x方向網(wǎng)格數(shù)1 600，z方向網(wǎng)格數(shù)600。

數(shù)值模擬結(jié)果表明，常規(guī)的雙程波一階應(yīng)力-速度方程零偏移距數(shù)值模擬結(jié)果(圖4-C)較單程波正演結(jié)果(圖4-B)在深層繞射波更完整(圖4-C，藍(lán)色箭頭所示)，但多次波發(fā)育(圖4-C，藍(lán)色方框所示)；無反射條件的雙程波零偏移距數(shù)值模擬結(jié)果中多次反射得到了較好的壓制(圖4-D，藍(lán)色方框所示)，中深層信噪比顯著提高。采用裂步傅里葉加插值法進(jìn)行的深度偏移處理結(jié)果表明，無反射雙程波正演、裂步傅里葉加插值法深度偏移后陡傾角鹽丘邊界成像質(zhì)量明顯提高(圖5，藍(lán)色方框所示)。

圖4 Sigsbee速度模型及零偏移距正演結(jié)果對(duì)比Fig.4 Sigsbee velocity model and comparison of zero offset simulation results(A)速度模型; (B)單程波數(shù)值模擬結(jié)果; (C)常規(guī)雙程波動(dòng)方程數(shù)值模擬結(jié)果;(D)無反射雙程波動(dòng)方程數(shù)值模擬結(jié)果。1 ft=30.48 cm

圖5 Sigsbee模型疊后深度偏移結(jié)果對(duì)比Fig.5 Comparison of post-stack depth migration in Sigsbee model(A)相移加插值法單程波正演、裂步傅里葉加插值法深度偏移結(jié)果;(B)無反射雙程波正演、裂步傅里葉加插值法深度偏移結(jié)果。1 ft=30.48 cm

3 結(jié) 論

a.使模型空間內(nèi)聲阻抗為常量進(jìn)行全波一階應(yīng)力-速度方程零偏移距數(shù)值模擬可有效壓制多次反射波，在較大網(wǎng)格間距時(shí)，虛譜法仍可有效抑制數(shù)值頻散。

b. 無反射雙程波一階應(yīng)力-速度方程零偏移距數(shù)值模擬可以適應(yīng)速度的空間突變和陡傾角地層，是復(fù)雜構(gòu)造區(qū)零偏移距地震波場(chǎng)模擬的有效手段。

[1] Claerbout J F. Imaging the Earth’s Interior [M]. London: Black Well Scientific Publications Inc, 1985.

[2] Sherif R E, Geldart L P. Exploration Seismology [M]. New York: Cambridge University Press, 1995.

[3] Aki K, Richards P. Quantitative Seismology [M]. San Francisco: W H Freeman Co, 1980.

[4] Stolt R H. Migration by Fourier transform[J]. Geophysics, 1978, 43(1): 23-48.

[5] Gazdag J. Wave equation migration with the phase-shift method[J]. Geophysics, 1978, 43(7): 1342-1351.

[6] Stoffa P L, Fokkema J T, de Luna Freire R M,etal. Split-step Fourier migration[J]. Geophysics, 1990, 55(4): 410-421.

[7] Kessinger W. Extended split-step Fourier migration[C]//1992 SEG Annual Meeting. Society of Exploration Geophysicists, 1992: 917-920.

[8] Ristow D, Rühl T. Fourier finite-difference migration[J]. Geophysics, 1994, 59(12): 1882-1893.

[9] 劉文革，賀振華，黃德濟(jì)，等. 高精度傅立葉有限差分法模型正演[J]. 石油地球物理勘探，2007，42(6)：629-633.

Liu W G, He Z H, Huang D J,etal. Forward modeling by Fourier finite difference approach[J]. Oil Geophysical Prospecting, 2007, 42(6): 629-633. (In Chinese)

[10] Baysal E, Kosloff D D, Sherwood J W C. A two-way nonreflecting wave equation[J]. Geophysics, 1984, 49(2): 132-141.

[11] 單延明，吳清嶺，于承業(yè)，等. 復(fù)雜地質(zhì)構(gòu)造波動(dòng)方程數(shù)值模擬[J].物探化探計(jì)算技術(shù)，2002，24(1)：22-26.

Shan Y M, Wu Q L, Yu C Y,etal. Wave equation modeling for complex media[J]. Techniques for Geophysical and Geochemical Exploration, 2002, 24(1): 22-26. (In Chinese)

[12] Corrêa G J P, Spiegelman M, Carbotte S,etal. Centered and staggered Fourier derivatives and Hilbert transforms[J]. Geophysics, 2002, 67(5): 1558-1563.

[13] Chung E T, Lam C Y, Qian J. A staggered discontinuous Galerkin method for the simulation of seismic waves with surface topography[J]. Geophysics, 2015, 80(4): T119-T135.

[14] Amundsen L, Robertsson J O A. Wave equation processing using finite-difference propagators, Part 1: Wavefield dissection and imaging of marine multicomponent seismic data[J]. Geophysics, 2014, 79(6): T287-T300.

[15] 杜增利，徐峰，高宏亮.虛譜法交錯(cuò)網(wǎng)格地震數(shù)值模擬[J].石油物探，2010，49(5)：430-437.

Du Z L, Xu F, Gao H L. Pseudo spectral seismic wavefield simulation with staggered grid [J]. Geophysical Prospecting for Petroleum, 2010, 49(5): 430-437. (In Chinese)

[16] 牟永光，裴正林.三維復(fù)雜介質(zhì)地震數(shù)值模擬[M].北京：石油工業(yè)出版社，2005.

Mu Y G, Pei Z L. Seismic Numerical Modeling for 3-D Complex Media[M]. Beijing: Petroleum Industry Press, 2005. (In Chinese)

[17] 杜增利，李亞林，尹成，等. 虛譜法速度應(yīng)力方程地震數(shù)值模擬[J].石油地球物理勘探，2009，44(5)：637-641.

Du Z L, Li Y L, Yin C,etal. Pseudo spectral seismic simulation with one order stress-velocity equation[J]. Oil Geophysical Prospecting, 2007, 42(6): 637-641. (In Chinese)

[18] Collino F, Tsogka C. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J]. Geophysics, 2001, 66(1): 294-307.

Non-reflection numerical simulation of zero offset seismic wavefield with velocity-stress equation based on pseudo spectral method

DU Zeng-li1, LI Chun-hong2, FU Yu-meng3, LI Yong-zhang2

1.SchoolofGeosciencesandTechnology,SouthwestPetroleumUniversity,Chengdu610500,China; 2.No.1GeophysicalExplorationCompany,DaqingDrilling&ExplorationEngineeringCorporation,CNPC,Daqing163357,China; 3.SchoolofArchitecture,Chang’anUniversity,Xi’an710018,China

Since the wavenumber inzaxis may not be calculated directly and correctly when P-wave velocity varies dramatically inx-yplane, zero offset seismic numerical simulation based on one-way wave equation is not suitable in the complex geological areas, and the maximum imaging angle of wavefield extrapolation operators based on one-way wave equation is limited. According to the formation mechanism of seismic reflection, the density distribution is calculated manually to make the P-wave impedance is constant in the model space, thus the multiple reflections generated from impedance interfaces are decreased. The precision of spatial derivative of the seismic wave field calculated by pseudo spectral method is higher than that of finite difference, therefore, the numerical dispersion is suppressed effectively, and the vertical resolution is enhanced. Numerical simulation results show that the diffraction wave field in zero offset seismic data is abundant, and the effect of seismic migration is better than that of frequency wave number domain or frequency space domain based on one-way wave equation.

pseudo spectral method; finite difference; numerical simulation; non-reflection; zero offset seismic data

10.3969/j.issn.1671-9727.2016.05.12

1671-9727(2016)05-0617-06

2015-10-14。

國(guó)家科技重大專項(xiàng)(2011ZX05049)。

杜增利(1969-)，男，博士，副研究員，研究方向：地球物理勘探, E-mail:duzl0896@126.com。

P631.4

A