方 剛，F(xiàn)OMEL Sergey，杜啟振

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東青島266580;2.Bureau of Economic Geology，John A.and Katherine G.Jackson School of Geosciences，The University of Texas at Austin，Austin，Texas 78731，USA)

地震波場延拓方法對(duì)于地震正演、偏移成像以及全波形反演具有重要作用。地震波場延拓可在時(shí)間域和頻率域進(jìn)行。對(duì)于時(shí)間域延拓而言，常用的兩類延拓方法是譜方法［1-2］和有限差分方法［3］。相對(duì)于有限差分方法，譜方法具有很高的計(jì)算精度，能夠有效壓制數(shù)值頻散［4-5］，但是較大的計(jì)算量限制了譜方法在實(shí)際生產(chǎn)中的應(yīng)用。有限差分法因其具有較高的計(jì)算效率且易于并行實(shí)現(xiàn)，被廣泛用于地震偏移成像。常規(guī)有限差分通過泰勒級(jí)數(shù)展開得到差分系數(shù)［6-7］，在數(shù)值計(jì)算中需要注意數(shù)值頻散和穩(wěn)定性。為了抑制數(shù)值頻散，空間網(wǎng)格的大小通常比地震波長小8～10倍;而為了滿足穩(wěn)定性條件，時(shí)間采樣間隔就需取得較小。因此，不管空間和時(shí)間，一般都是過采樣的，計(jì)算量大。針對(duì)地震波場延拓中有限差分方法存在的問題，許多學(xué)者提出了一系列改進(jìn)的策略，其核心思想是采用優(yōu)化的有限差分算子控制計(jì)算過程中的誤差［8-9］。Liu 和 Sen［10］基于時(shí)間－空間域的頻散關(guān)系和平面波理論提出了求解雙程波動(dòng)方程的有限差分格式，隨后又提出自適應(yīng)的變長度空間算子［11］，在保證精度的前提下節(jié)省計(jì)算量。杜啟振等［12］通過引入強(qiáng)約束條件和弱約束條件，構(gòu)造不同的Lagrange函數(shù)得到優(yōu)化的差分系數(shù)。上述優(yōu)化方法都是通過擬合精確的頻散關(guān)系得到優(yōu)化的差分系數(shù)，能夠提高計(jì)算精度和穩(wěn)定性，但本質(zhì)上這些方法還是在零頻率處的展開。Song等［13］基于Lowrank波場延拓算子［14］提出了一種新的優(yōu)化差分系數(shù)方法，其思想是通過設(shè)計(jì)有限差分系數(shù)，使之在較寬的波數(shù)范圍上匹配空間－波數(shù)域算子，因其匹配的目標(biāo)算子是基于Lowrank分解得到的，故該方法稱為Lowrank有限差分法。Fang等［15］將這種思想發(fā)展到交錯(cuò)網(wǎng)格，提出了一種在交錯(cuò)網(wǎng)格上求取優(yōu)化差分系數(shù)的交錯(cuò)網(wǎng)格Lowrank有限差分方法。筆者基于交錯(cuò)網(wǎng)格Lowrank有限差分理論，通過頻散分析比較交錯(cuò)網(wǎng)格Lowrank有限差分方法和常規(guī)有限差分方法的頻散曲線，得到交錯(cuò)網(wǎng)格Lowrank有限差分方法能夠壓制數(shù)值頻散的理論依據(jù)，并將交錯(cuò)網(wǎng)格Lowrank有限差分方法應(yīng)用于逆時(shí)偏移成像。

1 交錯(cuò)網(wǎng)格Lowrank有限差分

一階速度－應(yīng)力方程考慮速度和密度的變化，常用于地震波正演和復(fù)雜介質(zhì)的偏移成像。對(duì)于二維無耗散聲波介質(zhì)有

式中，u為以u(píng)x和uz為分量的質(zhì)點(diǎn)速度;p為聲波壓力;ρ(x)為介質(zhì)密度;c(x)為介質(zhì)的速度;x=(x，z)表示空間位置坐標(biāo)。

將式(1)中涉及到的物理量定義在交錯(cuò)網(wǎng)格上，聲壓定義在主網(wǎng)格上，速度定義在半網(wǎng)格點(diǎn)上，可以得到式(1)的離散形式:

其中

式中的偏微分算子可以表示為kx-space算子［15-16］:

其中，F(xiàn)表示空間的Fourier變換;kx和kz表示空間的頻率分量，算子中的指數(shù)函數(shù)項(xiàng)使得計(jì)算點(diǎn)沿著空間網(wǎng)格坐標(biāo)軸的正向或負(fù)向交錯(cuò)半個(gè)網(wǎng)格，這里用“+”表示沿著坐標(biāo)軸正向交錯(cuò)半個(gè)網(wǎng)格，用“-”表示沿著坐標(biāo)軸負(fù)方向交錯(cuò)半個(gè)網(wǎng)格。直接求解式(3)需要的計(jì)算量約為O(N2x)，其中Nx為計(jì)算區(qū)域的網(wǎng)格點(diǎn)數(shù)，對(duì)于實(shí)際的地震勘探數(shù)據(jù)，這樣的計(jì)算量是難以承受的。筆者根據(jù)Lowrank有限差分的思想［13］，提出計(jì)算式(3)的有限差分格式，可以降低計(jì)算量。

根據(jù)Lowrank分解的思想，可將式(4)定義的矩陣分解為矩陣的乘積［14］:

其中，W1=Wx(x，km)，由Wx中的某m列組成;W2=Wx(xn，k)，由Wx中的某n行組成，m和n是依賴于矩陣秩的整數(shù)，通常m和n遠(yuǎn)小于計(jì)算點(diǎn)數(shù)Nx。A是一個(gè)對(duì)角陣。Fomel等［14］給出了Lowrank分解式(5)的理論基礎(chǔ)和詳細(xì)的實(shí)施步驟。在Lowrank分解的基礎(chǔ)上，Song等［13］提出W2可以進(jìn)一步分解為兩個(gè)矩陣的乘積:

式中，B為一個(gè)L×Nx的矩陣，L為有限差分算子的長度。本文中選取B為如下形式:

其中 ξ =(ξ1(l)，ξ2(l))是一個(gè)二維向量，l=0，1，2，…L;Λ =(k1Δ1，k2Δ2)T，這里 Δj是 j方向空間網(wǎng)格間距;kj是波數(shù)向量在j方向的分量;對(duì)于二維情況，j=1，2 分別表示空間 x，z方向。

式(6)中C是N×L的矩陣，由下式計(jì)算得到:

其中“?”表示矩陣的廣義逆。事實(shí)上，式(8)給出了如下優(yōu)化問題的最小二乘解:

將式(6)代入式(5)得

考慮到式(7)，則式(10)的分量形式為

定義上式中關(guān)于m和n的求和為新的Nx×L矩陣:

其矩陣表示形式為

則

即

根據(jù)歐拉公式，式(16)可進(jìn)一步表示為

可得到求解x方向偏導(dǎo)數(shù)算子的交錯(cuò)網(wǎng)格Lowrank有限差分格式

其中 xR=(x1+l1Δ1，x2+l2Δ2)，xL=(x1-(l1-1)Δ1，x2+l2Δ2)，G(x，l)為有限差分的差分系數(shù)，該差分系數(shù)隨著模型自適應(yīng)地變化，使式(19)具有較高的計(jì)算精度。

類似地可以推導(dǎo)出式(3)中其他偏導(dǎo)數(shù)對(duì)應(yīng)的交錯(cuò)網(wǎng)格Lowrank有限差分算子。對(duì)于沿著x軸負(fù)方向交錯(cuò)的空間導(dǎo)數(shù)算子，其交錯(cuò)網(wǎng)格Lowrank有限差分算子為

這里 xR=(x1+(l1-1)Δ1，x2+l2Δ2)，xL=(x1-l1Δ1，x2+l2Δ2)。對(duì)于沿著 z軸正方向交錯(cuò)的空間導(dǎo)數(shù)算子，其交錯(cuò)網(wǎng)格Lowrank有限差分算子為

其中 xU=(x1+l1Δ1，x2+l2Δ2)，xD=(x1+l1Δ1，x2-(l2-1)Δ2)。對(duì)于沿著z軸負(fù)方向交錯(cuò)的空間導(dǎo)數(shù)算子，其交錯(cuò)網(wǎng)格Lowrank有限差分算子為

其中 xU=(x1+l1Δ1，x2+(l2-1)Δ2)，xD=(x1+l1Δ1，x2-l2Δ2)。

交錯(cuò)網(wǎng)格Lowrank有限差分算子式(19)～(22)是對(duì)交錯(cuò)網(wǎng)格kx-space算子(3)的近似。理論上有限差分算子的長度L越長，有限差分算子就越逼近于kx-space算子，計(jì)算精度越高。不難發(fā)現(xiàn)，交錯(cuò)網(wǎng)格Lowrank有限差分式(19)～(22)的計(jì)算量為O(LNx)，遠(yuǎn)小于直接計(jì)算kx-space算子式(3)的計(jì)算量(O())。相對(duì)于常規(guī)有限差分，Lowrank有限差分具有更高的精度，能夠有效地壓制頻散。因此，雖然Lowrank有限差分和常規(guī)有限差分具有相同的計(jì)算量，但是在實(shí)際的地震波場數(shù)值延拓過程中，Lowrank有限差分可以取較大的時(shí)間步長(通?？梢匀〉匠Ｒ?guī)有限差分時(shí)間步長的1.5～2倍)，在達(dá)到相同的計(jì)算精度時(shí)，節(jié)省了計(jì)算量。

2 交錯(cuò)網(wǎng)格Lowrank有限差分精度分析

為了分析交錯(cuò)網(wǎng)格Lowrank有限差分的計(jì)算精度，用平面波理論研究不同速度下式(19)的頻散關(guān)系，將一階速度 應(yīng)力方程的平面波解

代入式(2)，并且考慮交錯(cuò)網(wǎng)格Lowrank有限差分的計(jì)算格式(19)，根據(jù)頻散關(guān)系ω|k|v，定義交錯(cuò)網(wǎng)格Lowrank有限差分方法對(duì)應(yīng)的相速度

根據(jù)式(24)，可以計(jì)算不同介質(zhì)速度對(duì)應(yīng)的相速度。采用相速度和介質(zhì)速度的比值r=vLFD/v刻畫有限差分格式的數(shù)值頻散。如果r=1，表明有限差分格式不存在數(shù)值頻散;如果r偏離1越大，表明數(shù)值頻散越嚴(yán)重，計(jì)算精度越低。

圖1比較了不同介質(zhì)速度(200，2500，3000和3500 m/s)下，常規(guī)8階交錯(cuò)網(wǎng)格有限差分和交錯(cuò)網(wǎng)格Lowrank有限差分的頻散曲線，其中時(shí)間采樣間隔Δt=1 ms，空間采樣間隔Δx=10 m。對(duì)于一維的8階交錯(cuò)網(wǎng)格Lowrank有限差分，L=4。

從圖1可以看出，常規(guī)的交錯(cuò)網(wǎng)格有限差分只在零頻率(波數(shù))附近是精確的，隨著頻率增大，常規(guī)有限差分的頻散越來越嚴(yán)重。交錯(cuò)網(wǎng)格Lowrank有限差分能夠在較寬的波數(shù)(小于70% 的Nyquist波數(shù))范圍內(nèi)精確地刻畫頻散關(guān)系，這是該方法能夠較好地壓制數(shù)值頻散，具有較高精度的原因。

圖1 不同速度的頻散曲線Fig.1 Dispersion curves with different velocity

交錯(cuò)網(wǎng)格Lowrank有限差分法可用于實(shí)現(xiàn)高精度的時(shí)間域波場外推。比較常規(guī)交錯(cuò)網(wǎng)格有限差分和交錯(cuò)網(wǎng)格Lowrank有限差分進(jìn)行波場延拓的效果?？紤]光滑變化的二維模型，介質(zhì)速度介于800和4563 m/s之間，其分布由函數(shù)

給出;介質(zhì)的密度介于1100和2563 g/cm3之間，其分布由函數(shù)

給出。采用主頻為20 Hz的雷克子波，在模型中間位置激發(fā)，Δx= Δz=10 m，Δt=1.5 ms。

圖2為常規(guī)交錯(cuò)網(wǎng)格有限差分和交錯(cuò)網(wǎng)格Lowrank有限差分得到的波場快照。圖3為圖2所示的波場快照的中間道?？梢钥闯觯Ｒ?guī)交錯(cuò)網(wǎng)格有限差分精度較低，所構(gòu)建的波場存在嚴(yán)重的數(shù)值頻散，而交錯(cuò)網(wǎng)格Lowrank有限差分具有較高的精度，可以構(gòu)建出頻散很弱的波場。

圖2 不同方法在速度和密度變化的模型中正演得到的波場快照Fig.2 Wavefield snapshots in a variable velocity and density field

圖3 圖2所示波場快照的中間道比較Fig.3 Comparison of the middle trace from wavefield snapshot shown in fig.2

3 交錯(cuò)網(wǎng)格Lowrank有限差分法逆時(shí)偏移

逆時(shí)偏移［17-19］采用雙程波算子進(jìn)行波場外推，能夠構(gòu)建出全波場的信息，適于復(fù)雜構(gòu)造的地震成像。波場延拓是逆時(shí)偏移成像的核心，直接決定逆時(shí)偏移成像精度和計(jì)算效率。給出用交錯(cuò)網(wǎng)格Lowrank有限差分構(gòu)建地震波場，實(shí)現(xiàn)逆時(shí)偏移的方法步驟。

3.1 PML邊界條件

在地震逆時(shí)偏移成像時(shí)需要在一個(gè)有限區(qū)域內(nèi)模擬地震波在半無限大地球介質(zhì)中傳播的情形。實(shí)際計(jì)算中，基于交錯(cuò)網(wǎng)格的有限差分法波場延拓通常采用PML邊界條件［20］，對(duì)傳播至計(jì)算區(qū)域之外的波場進(jìn)行吸收衰減。本文基于交錯(cuò)網(wǎng)格Lowrank有限差分，實(shí)現(xiàn)PML邊界條件。在PML邊界內(nèi)的一階速度 應(yīng)力方程為

式中的阻尼項(xiàng)［20］可以取為

其中δj表示j方向PML層的寬度，α0為阻尼參數(shù)，可根據(jù)經(jīng)驗(yàn)選取。式(25)中的偏導(dǎo)數(shù)算子采用式(19)中定義的Lowrank有限差分算子計(jì)算。

3.2 逆時(shí)偏移實(shí)現(xiàn)步驟

逆時(shí)偏移的實(shí)現(xiàn)包括三個(gè)步驟:震源波場的正向延拓，檢波波場的逆時(shí)延拓和成像條件的應(yīng)用。利用交錯(cuò)網(wǎng)格Lowrank有限差分方法實(shí)現(xiàn)逆時(shí)偏移，具體步驟如下:

(1)根據(jù)式(5)對(duì)混合域算子(4)應(yīng)用Lowrank分解，得到對(duì)應(yīng)的分解矩陣;

(2)利用步驟(1)得到的分解矩陣，根據(jù)式(8)計(jì)算不同點(diǎn)上的有限差分系數(shù)G(x，l)，它依賴于速度模型以及時(shí)、空采樣間隔;

(3)利用G(x，l)，根據(jù)公式(2)和(10)進(jìn)行波場的正向和反向外推;

(4)在波場外推過程中，當(dāng)波場進(jìn)入PML邊界計(jì)算區(qū)域之后，采用公式(25)處理邊界內(nèi)的波場;

(5)應(yīng)用成像條件得到偏移成像結(jié)果。

4 Marmousi模型試算

Marmousi模型于1988年根據(jù)Cuanza盆地的構(gòu)造特征建立，包含劇烈的速度橫向變化和不整合界面等精細(xì)的構(gòu)造，常用于測(cè)試偏移方法對(duì)復(fù)雜構(gòu)造的成像能力，其速度模型如圖4所示。為了驗(yàn)證本文中方法的有效性，用該模型進(jìn)行偏移成像測(cè)試，采用主頻為20 Hz的Ricker子波作為震源，一共激發(fā)65炮，炮間隔為100 m，最大炮檢距2400 m，記錄的時(shí)間長度為4 s，時(shí)間采樣間隔為1.5 ms，水平和垂直方向的采樣間隔均為10 m，用4階交錯(cuò)網(wǎng)格Lowrank有限差分方法進(jìn)行波場延拓，計(jì)算所需的秩為4。圖5為常規(guī)交錯(cuò)網(wǎng)格有限差分和交錯(cuò)網(wǎng)格Lowrank有限差分得到的單炮記錄，圖6為圖5局部的放大顯示。

圖4 Marmousi速度模型Fig.4 Marmousi velocity model

圖5 不同方法得到的單炮記錄Fig.5 Shot records generated by different methods

可以看出，常規(guī)的交錯(cuò)網(wǎng)格有限差分法的結(jié)果(圖5(a))存在數(shù)值頻散(白色箭頭所示)，如果將其用于逆時(shí)偏移會(huì)影響成像精度，而交錯(cuò)網(wǎng)格Lowrank有限差分法能效壓制數(shù)值頻散，得到精度較高的波場記錄(圖5(b))。

圖7為交錯(cuò)網(wǎng)格Lowrank有限差分法逆時(shí)偏移成像結(jié)果。可以看出，大的構(gòu)造起伏以及局部細(xì)小的反射界面都得到比較精細(xì)的刻畫，振幅能量均一，油藏所在位置同相軸極性正確。交錯(cuò)網(wǎng)格Lowrank有限差分法逆時(shí)偏移具有較高的計(jì)算效率，采用1.5 ms的時(shí)間間隔，單炮偏移平均耗時(shí)77.02 s。為了避免數(shù)值頻散，常規(guī)有限差分法逆時(shí)偏移，需要選擇1 ms的時(shí)間間隔，單炮偏移需要109.34 s。相對(duì)于常規(guī)逆時(shí)偏移，交錯(cuò)網(wǎng)格Lowrank有限差分法逆時(shí)偏移的計(jì)算效率提高了42%。

5 結(jié)論

(1)交錯(cuò)網(wǎng)格Lowrank有限差分方法不僅進(jìn)一步提高了Lowrank有限差分的精度，而且可以處理變密度問題。

(2)交錯(cuò)網(wǎng)格Lowrank有限差分方法能在較寬的波數(shù)范圍內(nèi)精確擬合頻散關(guān)系，從而有效地壓制數(shù)值頻散;具有更高的計(jì)算精度，能夠準(zhǔn)確構(gòu)建地震波場，是實(shí)現(xiàn)高精度成像的基礎(chǔ);具有較高的計(jì)算穩(wěn)定性，在保證計(jì)算精度的前提下，可以選用較大的時(shí)間步長，從而節(jié)省計(jì)算量和存儲(chǔ)量。

(3)交錯(cuò)網(wǎng)格Lowrank有限差分法逆時(shí)偏移成像數(shù)值算例表明，該方法適用于變化劇烈的速度場，可對(duì)復(fù)雜構(gòu)造進(jìn)行高精度成像，在保證成像精度的前提下節(jié)省計(jì)算量。

圖6 圖5中部分區(qū)域放大顯示Fig.6 Zooms in part of fig.5

圖7 交錯(cuò)網(wǎng)格Lowrank有限差分法逆時(shí)偏移成像結(jié)果Fig.7 Image of reverse time migration using staggered grid Lowrank finite difference

［1］ TAL-EZER H，KOSLOFF D，KOREN Z.An accurate scheme for seismic forward modeling［J］.Geophysical Prospecting，1987，35:479-490.

［2］ RESHEF M，KOSLOFF D，EDWARDS M，et al.Threedimensional acoustic modeling by the Fourier method［J］.Geophysics，1988，53(9):1175-1183.

［3］ ETGEN J.High order finite-difference reverse time migration with the two way nonreflecting wave equation［R］.SEP-48:Stanford Exploration Project，1986:133-146.

［4］ CHU C，STOFFAP L.Apseudospectral-finite difference hybrid approach for large-scale seismic modeling and rtm on parallel computers［C］.78th Annual InternationalMeeting，SEG，Expanded Abstracts，2008，2087-2091.

［5］ ETGEN J，BRANDSBERG-DAHL S.The pseudo-analytical method:application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation［C］.79th Annual InternationalMeeting，SEG，2009，2552-2556.

［6］ DABLAIN M A.The application of high-order differencing to the scalar wave equation［J］.Geophysics，1986，51(1):54-66.

［7］ KINDELAN M，KAMEL A，SGUAZZERO P.On the construction and efficiency of staggered numerical differentiators for the wave equation［J］.Geophysics，1990，55(1):107-110.

［8］ TAKEUCHI N，GELLER R J.Optimally accurate second order time-domain finite difference scheme for computing synthetic seismograms in 2-D and 3-D media［J］.Phys Earth Planet Int，2000，119:99-131.

［9］ CHU C，STOFFA P L，SEIF R.3D elastic wave modeling using modified high-order time stepping schemes with improved stability condition［C］.79th Annual International Meeting，SEG，Expanded Abstracts，c2009.

［10］ LIU Y，SEN M K.A new time-space domain high-order finite-difference method for the acoustic wave equation［J］.Journal of Computational Physics，2009，228(23):8779-8806.

［11］ LIU Y，SEN M K.Finite-difference modeling with adaptive variable-length spatial operators［J］.Geophysics，2011，76(4):T79-T89.

［12］ 杜啟振，白清云，李賓.橫向各向同性介質(zhì)優(yōu)化差分系數(shù)法地震波場數(shù)值模擬［J］.石油地球物理勘探，2010，45(2):170-176.

DU Qi-zhen，BAI Qing-yun，LI Bin.Seismic wavefield modeling with optimal finite difference method in VTI media［J］.Oil Geophysical Prospecting，2010，45(2):170-176.

［13］ SONG X，F(xiàn)OMEL S，YING L.Low-rank finite-differences and low-rank Fourier finite-differences for seismic wave extrapolation［J］.Geophysical Journal International，2013，192:1-10.

［14］ FOMEL S，YING L，SONG X.Seismic wave extrapolation using low-rank symbol approximation［J］.Geophysical Prospecting，2013，61(3):526-536.

［15］ FANG G，F(xiàn)OMEL S，DU Q，HU J.Seismic wave extrapolation on a staggered grid using Lowrank decomposition and Lowrank finite-differences［C］.83th Annual InternationalMeeting， SEG， Expanded Abstracts，c2013.

［16］ TABEI M，MAST T D，WAAG R C.A k-space method for coupled first-order acoustic propagation equations［J］.The Journal of the Acoustical Society of America，2002，111:53-63.

［17］ BAYSAL E，KOSOFF D D，SHERWOOD J W C.Reverse time migration［J］.Geophysics，1983，48(11):1514-24.

［18］ WHITMORE D.Iterative depth migration by backward time propagation［C］.53th Annual International Meeting，SEG，Expanded Abstracts，c1983.

［19］ MCMECHAN G A.Migration by extrapolation of timedependent boundaryvalues［J］.Geophysical Prospecting，1983，31:413-420.

［20］ COLLINO F，TSOGKA C.Application of the perfectly matched absorbing layer model to the linear elastic dynamic problem in anisotropic heterogeneous media［J］.Geophysics，2001，66(1):294-307.