李海山 楊午陽(yáng) 雍學(xué)善

(①中國(guó)石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020； ②CNPC油藏描述重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州 730020)

1 問題的提出

理想地球介質(zhì)為彈性介質(zhì)，而實(shí)際地球介質(zhì)具有黏滯性，導(dǎo)致地震波在地層中傳播出現(xiàn)能量衰減、頻帶變窄、主頻降低、相位延遲等現(xiàn)象[1-3]，尤其對(duì)于近地表強(qiáng)衰減層或者是地下含氣層等強(qiáng)衰減區(qū)域，如果不考慮這些介質(zhì)的黏滯性，采用聲波方程正演模擬得到的地震波場(chǎng)與實(shí)際觀測(cè)波場(chǎng)之間的差異較大，造成全波形反演得到的介質(zhì)參數(shù)與實(shí)際介質(zhì)參數(shù)之間產(chǎn)生較大誤差[4,5]。因此，研究地震波在黏彈性介質(zhì)中傳播規(guī)律并探討?zhàn)椥越橘|(zhì)全波形反演方法具有重要意義。

通過對(duì)黏彈性介質(zhì)衰減機(jī)理及衰減規(guī)律的研究，人們提出了多種黏彈性模型，如廣義Maxwell體模型[6]、廣義Zener體模型[7]、Futterman模型[8]、廣義Kelvin-Voigt模型[9]、Kjartansson模型[10]、廣義準(zhǔn)線性固體(generalized standard linear solid，GSLS)模型[11]等。不同的黏彈性模型，其相應(yīng)的波動(dòng)方程也不同，如位移、位移—應(yīng)力、位移—速度—應(yīng)力和速度—應(yīng)力等不同形式的方程[12]； 波動(dòng)方程類型不同，則其相應(yīng)數(shù)值模擬方法也不同，實(shí)現(xiàn)的難易程度和計(jì)算量也不同。同時(shí)，研究表明地球介質(zhì)在地震頻帶范圍內(nèi)具有近似常Q特征[2]，即地下介質(zhì)品質(zhì)因子在地震頻帶范圍內(nèi)基本不隨頻率發(fā)生變化，而GSLS模型(圖1)可很好地近似這種常Q特征[13]；此外，GSLS模型相應(yīng)的波動(dòng)方程具有在時(shí)間域易于數(shù)值求解的優(yōu)點(diǎn)。因此，基于GSLS模型的黏彈性或黏聲波方程正被越來越多地應(yīng)用于正演模擬、逆時(shí)偏移及全波形反演等領(lǐng)域[14-19]。

Bai等[4，18]基于單個(gè)Maxwell體構(gòu)成的GSLS模型，采用二階黏聲波方程實(shí)現(xiàn)了時(shí)間域黏聲波全波形反演方法并進(jìn)行實(shí)際資料的全波形反演。研究表明由單個(gè)Maxwell體構(gòu)成的GSLS模型不足以近似地球介質(zhì)在地震頻帶范圍內(nèi)的常Q特征，為了更好地近似地球介質(zhì)在地震頻帶范圍內(nèi)的常Q特征，應(yīng)該采用由多個(gè)Maxwell體構(gòu)成的GSLS模型[13]，同時(shí)二階黏聲波方程不便于采用高階交錯(cuò)網(wǎng)格差分法求解。本文從Bohlen[20]基于GSLS模型得到的一階速度—應(yīng)力黏彈性波方程出發(fā)，得到一階速度—應(yīng)力黏聲波方程，并推導(dǎo)出相應(yīng)的速度梯度計(jì)算公式，同時(shí)采用共軛梯度法和高階交錯(cuò)網(wǎng)格有限差分法，實(shí)現(xiàn)了二維時(shí)間域黏聲波方程全波形反演方法，并通過模型測(cè)試驗(yàn)證了該方法的有效性。

2 一階速度—應(yīng)力黏聲波方程

Bohlen[20]基于GSLS模型(圖1)得到各向同性介質(zhì)中的一階速度—應(yīng)力黏彈性波動(dòng)方程，并進(jìn)行三維黏彈性波正演模擬方法研究。考慮到理想聲介質(zhì)中不存在剪切應(yīng)力，則由一階速度—應(yīng)力黏彈性波動(dòng)方程可得到如下的一階速度—應(yīng)力黏聲波方程

(1)

式中：p為壓力場(chǎng)；vi為波場(chǎng)速度分量；ρ為介質(zhì)密度；fp為震源項(xiàng)；L為GSLS模型中Maxwell體的個(gè)數(shù)；rl和τσl分別為第l個(gè)Maxwell體對(duì)應(yīng)的記憶變量和應(yīng)力松弛時(shí)間；τp為縱波松弛時(shí)間；M為松弛模量，且有

(2)

其中：符號(hào)R()表示取實(shí)部；ω0為參考頻率。若品質(zhì)因子Q已知，則可利用Blanch等[13]或Bohlen[20]給出的方法計(jì)算出最優(yōu)應(yīng)力松弛時(shí)間τσl和縱波松弛時(shí)間τp。式(1)相應(yīng)的應(yīng)力—位移黏聲波波動(dòng)方程為

(3)

式中ui為波場(chǎng)位移分量。若設(shè)置Maxwell體的個(gè)數(shù)為1，則該式就轉(zhuǎn)化為Bai等[4]采用的二階黏聲波方程。

3 目標(biāo)函數(shù)及模型參數(shù)梯度

全波形反演是利用全波場(chǎng)信息，通過一定的優(yōu)化方法估計(jì)地下介質(zhì)參數(shù)模型，使正演地震記錄與觀測(cè)地震記錄達(dá)到最佳擬合[21]。設(shè)目標(biāo)函數(shù)為波場(chǎng)殘差能量，即

(4)

式中δu=u-uobs，u=(vx,vz,p,rl)T為地震波場(chǎng)向量。目標(biāo)函數(shù)對(duì)模型參數(shù)求偏導(dǎo)，有

(5)

(6)

同理可由數(shù)據(jù)空間小的擾動(dòng)得到模型空間總的變化[22]

(7)

(8)

(9)

可見目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度等于由數(shù)據(jù)空間小的擾動(dòng)引起的模型空間的變化量。如果在式(3)中的各變量和參數(shù)都施加一階擾動(dòng)，同時(shí)把模型參數(shù)變化引起的數(shù)據(jù)空間的擾動(dòng)看作是新的震源，結(jié)合式(6)～式(9)推導(dǎo)可得松弛模量的梯度計(jì)算公式(詳細(xì)推導(dǎo)見附錄A)。松弛模量的梯度為

(10)

式中pford和pback分別為正演模擬壓力場(chǎng)與殘差逆時(shí)反傳壓力場(chǎng)。據(jù)式(2)縱波速度與松弛模量之間的關(guān)系，可得到縱波速度梯度

(11)

4 全波形反演策略

要使目標(biāo)函數(shù)達(dá)到最小，目前基本上都采用局部?jī)?yōu)化方法。在局部最優(yōu)化方法中，雖然共軛梯度法只利用一階導(dǎo)數(shù)信息，卻克服了最速下降法收斂慢的缺點(diǎn)，避免了牛頓法需要存儲(chǔ)和計(jì)算Hessian矩陣并求逆的缺點(diǎn)，因此得到廣泛的應(yīng)用[23，24]。本文采用共軛梯度法更新縱波速度模型，即

vP，n+1=vP，n+αφn

(12)

式中：n為迭代次數(shù)；α為更新步長(zhǎng)；φn為共軛方向

(13)

通過FR方法[22]，可計(jì)算得到

(14)

反演過程中采用固定步長(zhǎng)，保持Q模型和密度參數(shù)不變；在每次迭代過程中，正演波場(chǎng)和殘差逆時(shí)反傳波場(chǎng)用式(1)計(jì)算得到，用式(11)計(jì)算縱波速度梯度，波場(chǎng)計(jì)算采用8階交錯(cuò)網(wǎng)格有限差分法；為更好地近似常Q特征，設(shè)置GSLS模型中Maxwell體的個(gè)數(shù)為3[1]； 給定GSLS模型參數(shù)τσl，可利用Blanch等[13]給出的方法得到給定頻帶范圍內(nèi)的最優(yōu)縱波松弛時(shí)間τp。此外，若在全波形反演過程中令式(1)中的記憶變量及與黏彈性相關(guān)的參數(shù)為零，即可進(jìn)行二維聲波全波形反演。

5 方法模型測(cè)試

圖2a是從3D SEG/EAGE逆掩推覆模型中抽取出來的一個(gè)二維速度模型，網(wǎng)格尺寸為400m×187m，空間步長(zhǎng)為12m； 圖2b是對(duì)圖2a進(jìn)行平滑后得到的，作為全波形反演的初始速度模型； 圖2c是由圖2a所示速度模型映射后生成的Q模型；圖2d是對(duì)圖2c進(jìn)行平滑而生成的平滑Q模型。Q模型在近地表衰減最強(qiáng)，從淺到深衰減逐漸減弱。

利用圖2a所示速度模型和圖2c所示Q模型生成二維黏聲波觀測(cè)記錄，震源子波選取主頻為15Hz的雷克子波，炮點(diǎn)與接收點(diǎn)均置于距地面24m位置處，共采集80炮，每炮設(shè)置400個(gè)檢波點(diǎn)，炮間距為60m，檢波點(diǎn)距為12m，記錄長(zhǎng)度為1.5s，時(shí)間采樣間隔為1ms，圖3b是其第40個(gè)單炮記錄。采用相同參數(shù)，利用圖2a所示速度模型生成二維聲波觀測(cè)記錄，圖3a是對(duì)應(yīng)的第40個(gè)單炮記錄。與聲波炮記錄比較可見，黏聲波炮記錄隨旅行時(shí)增加能量快速衰減，到深層已看不到弱反射信息。

圖4給出了圖3中兩個(gè)矩形區(qū)域正演炮記錄相應(yīng)的振幅譜，其中藍(lán)線對(duì)應(yīng)聲波正演炮記錄，紅線對(duì)應(yīng)黏聲波正演炮記錄，明顯可見地震波的衰減隨著旅行時(shí)的增大而增強(qiáng)，同時(shí)可見地震波的主頻隨著旅行時(shí)的增大而降低。

首先，以聲波正演炮記錄為觀測(cè)記錄，從圖2b所示初始速度模型出發(fā)進(jìn)行聲波全波形反演，迭代50次后得到圖5a所示的縱波速度反演結(jié)果。然后，以黏聲波正演炮記錄為觀測(cè)記錄，分別進(jìn)行聲波及黏聲波全波形反演，迭代50次后得到圖5b和圖5c所示的縱波速度反演結(jié)果，黏聲波全波形反演過程中采用圖2c所示的真實(shí)Q模型。由圖5a可見，對(duì)于聲波觀測(cè)記錄，采用聲波全波形反演方法會(huì)得到較好反演結(jié)果，而由圖5b和圖5c可見，對(duì)于黏聲波觀測(cè)記錄，采用聲波全波形反演方法會(huì)得到錯(cuò)誤結(jié)果，反演得到縱波速度模型偏離真實(shí)縱波速度模型較遠(yuǎn)，只有采用黏聲波全波形反演方法才會(huì)得到較準(zhǔn)確反演結(jié)果。

在進(jìn)行實(shí)際資料全波形反演時(shí)，準(zhǔn)確的Q模型是得不到的，只能得到相對(duì)準(zhǔn)確的平滑Q模型，為研究本文反演方法對(duì)Q模型的敏感性，以黏聲波正演炮記錄為觀測(cè)記錄，采用圖2d所示的平滑Q模型進(jìn)行黏聲波全波形反演，反演結(jié)果見圖5d。比較圖2c和圖2d可見雖然平滑Q模型與真實(shí)Q模型有一定的偏差，但比較圖5c與圖5d可見采用平滑Q模型能夠得到與采用真實(shí)Q模型接近相同的反演結(jié)果，這表明在進(jìn)行實(shí)際資料全波形反演時(shí)，只要提供足夠準(zhǔn)確的平滑Q模型，就能顯著提高反演結(jié)果的準(zhǔn)確性。

圖3 第40炮正演記錄 (a)聲波炮記錄； (b)黏聲波炮記錄

圖4 對(duì)應(yīng)圖3中兩個(gè)矩形區(qū)域炮記錄的振幅譜 (a)藍(lán)色矩形框； (b)紅色矩形區(qū)域

圖5 全波形反演迭代50次得到的縱波速度模型 (a)聲波觀測(cè)—聲波全波形反演； (b)黏聲波觀測(cè)—聲波全波形反演； (c)黏聲波觀測(cè)— 黏聲波全波形反演—真實(shí)Q模型； (d)黏聲波觀測(cè)—黏聲波全波形反演—平滑Q模型

圖6a和圖6b分別為模型水平方向1.8km和3.0km位置處以黏聲波正演炮記錄為觀測(cè)記錄，反演前、后的單道速度曲線比較，可更明顯地看出黏聲波反演結(jié)果較接近真實(shí)模型，而聲波全波形反演結(jié)果誤差較大，且深度越大偏離真實(shí)模型越遠(yuǎn)。

圖6 真實(shí)模型及反演縱波速度單道曲線比較 (a)1.8km位置速度曲線； (b)3.0km位置速度曲線

圖7是以黏聲波正演炮記錄為觀測(cè)記錄，進(jìn)行全波形反演時(shí)歸一化數(shù)據(jù)殘差隨迭代次數(shù)的變化曲線，由圖可見聲波全波形反演方法和黏聲波反演方法都收斂了，但聲波全波形反演方法得到了錯(cuò)誤的反演結(jié)果，而黏聲波全波形反演方法得到的反演結(jié)果是正確的，迭代50次后觀測(cè)記錄與擬合記錄之間的誤差較小。因此，全波形反演中，反演誤差和收斂曲線只可作為一種參考，不能作為衡量反演結(jié)果準(zhǔn)確性的標(biāo)準(zhǔn)。

圖7 歸一化誤差隨迭代次數(shù)變化曲線

考慮到實(shí)際資料全波形反演中噪聲干擾是不可避免的，對(duì)于陸上資料該問題尤為嚴(yán)重，為分析隨機(jī)噪聲對(duì)本文反演方法的影響，在黏聲波觀測(cè)記錄中加入信噪比為10的隨機(jī)噪聲(圖8a)，采用圖2d所示的平滑Q模型進(jìn)行黏聲波全波形反演(圖9)。比較圖3b與圖8b可見，觀測(cè)炮記錄與反演炮記錄中反射同相軸的旅行時(shí)和波形匹配較好；同時(shí)，比較圖2a與圖9可見，當(dāng)觀測(cè)記錄中含有一定信噪比的隨機(jī)噪聲時(shí)，也能得到相對(duì)較好的縱波速度反演結(jié)果。

需要注意的是，本文在反演過程中采用的子波是準(zhǔn)確子波，反演方法的噪聲抑制能力也有待進(jìn)一步提高，為進(jìn)一步提高反演方法的性能，下一步可采用褶積型目標(biāo)函數(shù)[25]，并引入模型參數(shù)的正則化項(xiàng)[26]。

圖8 含噪黏聲波炮記錄(a)及其全波形反演炮記錄(b)

圖9 含噪黏聲波炮記錄黏聲波反演結(jié)果

6 結(jié)論

本文采用基于GSLS黏彈性模型的二維一階速度—應(yīng)力黏聲波方程，推導(dǎo)出相應(yīng)的縱波速度梯度計(jì)算公式，并采用高階交錯(cuò)網(wǎng)格有限差分法和共軛梯度法，實(shí)現(xiàn)了二維時(shí)間黏聲波全波形反演方法。該方法由于采用由多個(gè)Maxwell體構(gòu)成的GSLS模型，可以更好地近似地球介質(zhì)的常Q特征；由于在時(shí)間域?qū)崿F(xiàn)，因此波場(chǎng)正演和殘差逆時(shí)反傳比較直接快速。

模型測(cè)試結(jié)果驗(yàn)證了該方法的可行性和有效性。以黏聲波正演炮記錄為觀測(cè)記錄時(shí)，與聲波全波形反演結(jié)果相比，采用黏聲波全波形反演方法能夠得到較為準(zhǔn)確的縱波速度反演結(jié)果，可見黏彈性對(duì)全波形反演結(jié)果的準(zhǔn)確性和分辨率有較大的影響，因此建議在強(qiáng)衰減區(qū)域采用黏彈性或者黏聲波全波形反演方法。

附錄A 松弛模量梯度的推導(dǎo)

根據(jù)攝動(dòng)理論，正文式(3)中的各波場(chǎng)變量和模型參數(shù)都施加一階擾動(dòng)，即

(A-1)

則可得到新的應(yīng)力—位移黏聲波方程

(A-2)

式中新的震源項(xiàng)

(A-3)

如果把模型參數(shù)變化引起的擾動(dòng)看作是新的震源，則式(A-2)中位移變量δui的解為

δui(x,t)

(A-4)

式中Gi(x,t;x′,t′)為格林函數(shù)。如果令密度為常數(shù)，即只進(jìn)行縱波速度的全波形反演，則ΔF=0，將式(A-3)代入式(A-4)，有

(A-5)

對(duì)比式(6)可得

(A-6)

由式(7)可得

(A-7)

進(jìn)一步整理，有

(A-8)

定義新的波場(chǎng)

(A-9)

則式(A-8)可寫為

(A-10)

式中波場(chǎng)ψi(x,t)是由波場(chǎng)殘差δu′在接收點(diǎn)位置開始逆時(shí)反傳產(chǎn)生的。將式(A-10)中的隱式和展開，有

(A-11)

對(duì)于二維情形有

(A-12)

考慮到黏聲波波場(chǎng)正演模擬和殘差波場(chǎng)逆時(shí)反傳都采用一階速度—應(yīng)力黏聲波方程，所以式(A-12)中的位移分量需要用應(yīng)力來替換。利用應(yīng)力與位移關(guān)系式可得

(A-13)

式中pford和pback分別為正演模擬壓力場(chǎng)與殘差逆時(shí)反傳壓力場(chǎng)。

[1] Robertsson J O A，Blanch J O and Symes W W.Viscoelastic finite-difference modeling.Geophysics，1994，59(9)：1444-1456.

[2] Liu H P，Anderson D L and Kanamori H.Velocity dispersion due to an elasticity；implications for seismology and mantle composition.Geophysical Journal International，1976，47(1)：41-58.

[3] 陳志德，趙忠華，王成.黏滯聲學(xué)介質(zhì)地震波吸收補(bǔ)償疊前時(shí)間偏移方法.石油地球物理勘探，2016，51(2)：325-333. Chen Zhide，Zhao Zhonghua，Wang Cheng. Prestack time migration with compensation for absorption of viscous-elastic media.OGP，2016，51(2)：325-333.

[4] Bai J，Yingst D，Bloor R et al.Waveform inversion with attenuation.SEG Technical Program Expanded Abstracts，2012，31：1-5.

[5] 楊午陽(yáng)，王西文，雍學(xué)善等.地震全波形反演方法研究綜述.地球物理學(xué)進(jìn)展，2013，28(2)：766-776. Yang Wuyang，Wang Xiwen，Yong Xueshan et al.The review of seismic full waveform inversion method.Progress in Geophysics，2013,28(2)：766-776.

[6] Emmerich H and Korn M.Incorporation of attenuation into time domain computations of seismic wave fields.Geophysics，1987，52(9)：1252-1264.

[7] Schiessel H，Metzler R，Biumen A et al.Generalized viscoelastic models：their fractional equations with solutions.Journal of Physics A： Mathematical and General，1995，28(23)：6567-6584.

[8] Futterman W I.Dispersive body waves.Journal of Geophysical Research，1962，67(13)：5279-5291.

[9] Jiang F and Mehta A J.Mudbanks of the Southwest Coast of India Ⅳ： Mud viscoelastic properties.Journal of Coastal Research,1995，11(3)：918-926.

[10] Kjartansson E.Constant-Q wave propagation and attenuation.Journal of Geophysical Research，1979，84(B9)：4737-4748.

[11] Carcione J M，Kosloff D and Kosloff R.Wave propagation simulation in a linear viscoacoustic medium.Geophysical Journal International，1988，93(2)：393-407.

[12] Cao D.Viscoelastic wave equations based on different rheological model in time domain modeling.Beijing 2014 International Geophysical Conference & Exposition，2014，673-676.

[13] Blanch J O，Robertsson J O A and Symes W W.Mode-ling of a constant Q：Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique.Geophysics，1995，60(1)：176-184.

[14] Groos L，Sch?fer M，F(xiàn)orbriger T et al.On the significance of viscoelasticity in a 2D full waveform inversion of shallow seismic surface.74th EAGE Conference & Exhibition Incorporating SPE EUROPEC 2012，Copenhagen，Denmark，2012.

[15] Zhang G and Gao J.Applying the t-method to visco-elastic forward modeling.SEG Technical Program Expanded Abstracts，2013，32：3531-3536.

[16] Dutta G，Lu K，Wang X et al.Attenuation compensation in least-squares reverse time migration using the visco-acoustic wave equation.SEG Technical Program Expanded Abstracts，2013，32：3723-3725.

[17] Bai J and Yingst D.Q estimation through waveform inversion.75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013，London，UK，2013.

[18] Bai J，Yingst D，Bloor R et al.Viscoacoustic waveform inversion of velocity structures in the time domain.Geophysics，2014，79(3)：R103-R119.

[19] 李金麗，劉建勛，姜春香等.黏聲VTI介質(zhì)正演模擬與照明分析.石油地球物理勘探，2017，52(5)：906-914. Li Jinli，Liu Jianxun，Jiang Chunxiang et al.Forward modeling and illumination analysis in visco-acoustic VTI media.OGP，2017，52(5)：906-914.

[20] Bohlen T.Parallel 3-D viscoelastic finite difference seismic modeling.Computers & Geosciences，2002，28(8)：887-899.

[21] 白璐，韓立國(guó)，張盼等.最小平方濾波法時(shí)間域全波形反演.石油地球物理勘探，2016，51(4)：721-729. Bai Lu，Han Liguo，Zhang Pan et al.Time-domain full waveform inversion based on least square filter.OGP，2016，51(4)：721-729.

[22] Tarantola A.Inverse Problem Theory and Methods for Model Parameter Estimation.Society for Industrial and Applied Mathematics，2005.

[23] 張廣智，孫昌路，潘新朋等.快速共軛梯度法頻率域聲波全波形反演.石油地球物理勘探，2016，51(4)：730-737. Zhang Guangzhi，Sun Changlu，Pan Xinpeng et al.Acoustic full waveform inversion in the frequency domain based on fast conjugate gradient method.OGP，2016，51(4)：730-737.

[24] Kamei R and Pratt R G.Inversion strategies for visco-acoustic waveform inversion.Geophysical Journal International，2013,194(2)：859-884.

[25] Choi Y and Alkhalifah T.Source-independent time-domain waveform inversion using convolved wavefields： Application to the encoded multisource waveform inversion.Geophysics，2011，76(5)：R125-R134.

[26] Asnaashari A，Brossier R，Garambois S et al.Regula-rized seismic full waveform inversion with prior model information.Geophysics，2013，78(2)：R25-R36.