波動方程初至雙差走時層析反演

趙 磊,馮 波,王華忠

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.同濟大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與智能反演成像研究組,上海200092;3.同濟大學(xué)海洋高等研究院,上海200092)

地震波走時多用于中、大尺度的速度結(jié)構(gòu)特征反演。對于有限頻帶的地震信號,地震波走時(絕對走時)受震源子波的初相位、走時測量準(zhǔn)則(如起跳時刻、包絡(luò)極值時刻,最大峰值時刻)等多種因素的影響。地震波走時的測量誤差會降低反演結(jié)果的精度及可靠性。為盡可能降低絕對走時測量誤差對反演的影響,BRUNE等[1]提出采用“相對測量值”代替“絕對測量值”以減少測量誤差?；诖朔椒y量的地震走時亦稱為雙差走時(double-difference traveltime),可以用于提高震源定位精度[2-4],反演高精度速度模型[5]。

傳統(tǒng)的雙差走時層析反演方法多基于射線理論,存在焦散及陰影區(qū)等問題[6],且反演精度不高。對于小尺度異常體(速度非均勻體的尺度小于菲涅爾體的寬度),有限頻理論[7-14]可以更好地處理地震波的一階繞射效應(yīng)[15-16]。YUAN等[17]將雙差走時測量方法引入伴隨層析(adjoint tomography),得到了更高精度的反演結(jié)果,但雙差走時敏感度核函數(shù)仍采用Born近似計算方法計算得到[11]。相較于Born近似,Rytov近似可以更好地描述由于速度擾動引起的前向散射波相位擾動[18],因而更適用于走時反演。本文利用Rytov近似構(gòu)造走時敏感度核函數(shù)[19-26],即將雙差走時測量方法與Rytov走時敏感度核函數(shù)相結(jié)合,提出了一種新的波動方程初至層析反演方法。該方法在保留Rytov近似優(yōu)點的同時,可以降低甚至消除走時測量誤差對反演結(jié)果的影響。

1 理論

1.1 雙差相位延遲敏感度核函數(shù)

頻率域標(biāo)量聲波方程可以表示為:

(1)

令擾動后的速度場為v(x),總場為u(x;ω),總場與背景場的關(guān)系表示為:

(2)

式中:ψ(x;ω)為復(fù)相位。

利用Rytov近似,將單頻諧波的相位延遲表示為:

=〈Kp(x,ω;xr,xs),Δm(x)〉M

(3)

式中:Im表示取復(fù)數(shù)的虛部;Δm(x)為慢度平方擾動(屬于模型空間M),Δm(x)=v(x)-2-v0(x)-2;xs,xr分別代表震源和檢波器橫坐標(biāo);Kp(x,ω;xr,xs)為相位延遲敏感度核函數(shù)[25]。Kp(x,ω;xr,xs)可表示為:

Kp(x,ω;xr,xs)=

(4)

(3)式所預(yù)測的相位延遲假定了u和u0由相同的震源產(chǎn)生。若震源子波未知,復(fù)相位中包含子波剩余相位差。在全波形反演中,通常用褶積[27]或反褶積[28-31]類方法消除震源子波對反演的影響。本文采用反褶積方法,引入如下反褶積波場:

(5)

式中:xi和xj為同一炮集中任意不重合的兩個檢波器橫坐標(biāo);uw(xi,ω;xs)和uw(xj,ω;xs)分別為初至波形加窗前、后的地震信號(頻率域)。uw(xi,ω;xs)可表示為:

(6)

式中:h(t)為初至波形窗函數(shù),用于提取初至震相;*為卷積運算符號;h(ω)為窗函數(shù)的頻率響應(yīng)。

令uw(xi,ω;xs)=f(ω)A(xi,ω;xs)eiφ(xi,ω;xs),uw(xj,ω;xs)=f(ω)A(xj,ω;xs)eiφ(xj,ω;xs),其中,f(ω)為震源子波頻譜,A(ω),φ(ω)分別為格林函數(shù)的振幅譜和相位譜,初至波形的相位差Δφi,j(ω;xs)可以表示為:

Δφi,j(ω;xs)=φ(xi,ω;xs)-φ(xj,ω;xs)≡Im[lnvi,j(ω;xs)]

(7)

顯然,子波頻率對地震信號相位的影響被抵消了。采用雙差測量方法[5],得到的觀測數(shù)據(jù)和模擬數(shù)據(jù)的雙差相位延遲可以表示為:

(8)

綜合考慮(2)式至(7)式,初至波形的相位延遲和模型擾動的線性關(guān)系可表示為:

(9)

將(9)式代入(8)式,雙差相位延遲與模型擾動的線性關(guān)系可表示為:

(10)

(11)

1.2 雙差走時敏感度核函數(shù)

有限頻帶地震信號的走時擾動可以表示為單頻諧波的相位延遲的加權(quán)疊加[22],即:

(12)

同理,雙差走時擾動與雙差相位延遲的關(guān)系可表示為:

(13)

將(10)式、(11)式與(13)式相結(jié)合,可以建立雙差走時擾動與模型擾動的線性關(guān)系:

(14)

式中:Kdd(x;xi,xj,xs)為(帶限信號或有限頻)雙差走時敏感度核函數(shù)。Kdd(x;xi,xj,xs)的表達式如下:

(15)

式中:K(x;xr,xs)為基于Rytov近似的有限頻走時敏感度核函數(shù)[24-25]。

馮波等[25]給出了Rytov近似走時敏感度核函數(shù)K(x;xr,xs)在時間-空間域的顯式計算公式為:

(16)

1.3 走時層析反問題求解

基于走時殘差L2范數(shù)的誤差泛函可以表示為:

(17)

假定每個震源有Nr道地震記錄,隨機采用第i個地震道作為參考道,其它地震道與參考道計算得到的雙差走時為ΔΔti,j,滿足1≤j≤Nr且j≠i。

采用Gauss-Newton反演算法,將模型參數(shù)更新過程表示為:

mk+1=mk-αkP[H-1(mk)g(mk)]

(18)

由(15)式可知,泛函梯度g(mk)可以用核函數(shù)與走時殘差表示為:

(19)

(20)

(19)式可以表示為:

(21)

式中:λdd(x,T-t;xs)為雙差走時伴隨震源逆時傳播產(chǎn)生的伴隨波場。

為避免直接Hessian矩陣求逆,可以通過求解如下方程的近似解獲得模型更新方向:

H(mk)Δmk=g(mk)

(22)

本文采用共軛梯度(conjugate gradient,CG)法求解方程,采用馮波等[25]提出的隱式矩陣向量乘法得到高效的Hessian向量積,無需顯式計算及存儲Hessian矩陣。具體計算公式見附錄A。

2 數(shù)值試驗

2.1 高斯擾動模型測試

我們設(shè)計了一個含有高斯異常體的速度模型v(x,z)=v0+δv(x,z),其中背景模型v0為勻速模型,高斯異常體δv描述如下:

(23)

式中:ε為速度擾動百分比;v0=2500m/s為均勻背景速度;a為高斯異常體的尺度參數(shù);(x0,z0)為高斯異常體的中心坐標(biāo),(x0,z0)=(2500m,2500m)。

采用10m×10m的網(wǎng)格離散化含有高斯異常體的速度模型,x和z方向網(wǎng)格數(shù)目均為501。高斯異常體的尺度參數(shù)a=500m,速度擾動百分比ε=10%(圖1a)。為消除采集孔徑對反演結(jié)果的影響,本文設(shè)計了一個四周激發(fā)-接收的觀測系統(tǒng),在模型四周放置100個均勻分布的震源,每邊25個震源。每炮由100個檢波器接收,檢波器均勻分布在模型四周,并與震源重合。正演采用10Hz主頻的Ricker子波(對應(yīng)波長λ0=250m,與異常體尺度參數(shù)滿足a=2λ0)。

為了驗證本文雙差走時反演方法的有效性,我們首先對每個震源子波引入隨機延遲,并采用有限差分方法模擬觀測地震記錄。本文采用SEISCOPE數(shù)值優(yōu)化軟件包(SEISCOPE optimization toolbox[31])中的截斷牛頓算法(對于走時層析,截斷牛頓算法退化為高斯牛頓算法)求解雙差走時目標(biāo)函數(shù)。根據(jù)附錄A中的隱式矩陣向量積計算公式計算梯度與Hessian向量積。初始模型采用勻速背景(v0=2500m/s),震源子波采用0.5s延遲的Ricker子波(10Hz主頻)。目標(biāo)函數(shù)的終止準(zhǔn)則為檢測目標(biāo)函數(shù)的相對誤差(σ=J(mk)/J(m0))小于預(yù)先給定的門檻值σ=1.0×10-4。步長采用線性搜索估計,每個高斯牛頓方向采用2次CG內(nèi)迭代求解(16)式。經(jīng)過10次高斯牛頓方向迭代求解,得到的速度擾動(圖1b)最大值為249.7m/s,與真實速度擾動的最大值250.0m/s基本一致。

接著我們根據(jù)相同的反演參數(shù),利用絕對走時[25]層析反演方法進行反演,將正演觀測數(shù)據(jù)得到的隨機延遲Ricker子波作為震源,得到的速度擾動如圖1c所示?？梢钥闯?采用雙差走時與絕對走時層析反演方法得到的速度擾動均與真實速度模型吻合較好。因此對于完備的觀測系統(tǒng),雙差走時層析反演與常規(guī)反演結(jié)果基本一致,證明了雙差走時層析反演方法的有效性。

圖1 含有高斯異常體的速度模型(a)、采用雙差走時(b)與絕對走時(c)層析反演方法得到的速度擾動

2.2 Overthrust模型測試

為進一步測試本文方法在近地表速度建模中的效果,我們用SEG/EAGE Overthrust速度模型[32]正演觀測地震數(shù)據(jù)。原始速度模型在x和z方向網(wǎng)格數(shù)目分別為801和187,網(wǎng)格間距為25m。為更好地反演速度模型左側(cè)的推覆構(gòu)造,我們將速度模型左右側(cè)分別拓展,得到一個橫向長度為25000m的SEG/EAGE Overthrust近地表速度模型(圖2a)。我們設(shè)計了一個陸上單邊偏移距觀測系統(tǒng),在地表放置91個均勻分布的震源,炮間距為200m,每炮由117道檢波器接收,道間距為25m(最小偏移距100m,最大偏移距3000m)。觀測地震記錄由有限差分算法計算得到,地震子波主頻為8Hz且存在隨機延遲的Ricker子波。

由于觀測地震記錄較為復(fù)雜,本文采用自動初至拾取算法得到觀測數(shù)據(jù)和模擬數(shù)據(jù)的初至,并計算其雙差走時(參考道為最小偏移距地震道)。初始模型采用線性遞增的速度模型,當(dāng)z=0時,v=2400m/s,當(dāng)z=3000m時,v=4000m/s,震源子波采用8Hz主頻的Ricker子波(子波延遲時間為0.2s),采用截斷牛頓算法計算模型更新方向,每個方向1次CG內(nèi)迭代求解(16)式,并用線性搜索方法估計迭代步長。在反演過程中,采用了多偏移距反演策略(從最長偏移距(3000m)逐漸減少到最小偏移距(500m)),最終采用雙差走時層析反演方法得到的速度模型如圖2b所示。可以看出,大尺度的近地表速度結(jié)構(gòu)特征恢復(fù)得較好。圖2c為利用該反演方法得到的速度擾動,主要表現(xiàn)為近地表中-大尺度的速度更新,最大有效反演深度約為150m。

圖3為不同深度處速度橫向抽線對比結(jié)果,顯然雙差走時層析反演結(jié)果在淺層的分辨率較高,甚至能分辨中、小尺度的速度異常。隨著深度增加,其反演精度逐漸降低(有效反演深度由最大偏移距及速度結(jié)構(gòu)共同決定)。

圖2 橫向長度為25000m的SEG/EAGE Overthrust近地表速度模型(a)、采用雙差走時層析反演方法得到的速度模型(b)及速度擾動(c)

圖3 不同深度處速度橫向抽線對比結(jié)果

2.3 二維陸上實際地震資料測試

采用中國某地二維陸上實際地震資料(共1638炮,雙邊接收,最大偏移距7132m;炮間隔60m,道間距10m)進行測試,該地震資料采用初至自動拾取方法得到(拾取起震時刻)。初始模型采用線性遞增速度模型(起伏地表之上用空氣速度填充),震源子波為主頻為10Hz的Ricker子波,基于相同的初至拾取標(biāo)準(zhǔn)得到模擬數(shù)據(jù)的初至,并采用多偏移距反演策略(最大偏移距分別為3000m,2000m,1000m,500m,每個尺度迭代8次),得到如圖4所示的最終反演結(jié)果,可以看出最大有效反演深度約為700m。從反演結(jié)果可以看出,山間低速帶得到了較好的刻畫(近地表速度低至800m/s,低速帶厚度約為200m),這與野外實際地表露頭結(jié)果相符,速度模型左側(cè)的平原地帶,反演結(jié)果表現(xiàn)出明顯的成層性。

圖4 二維陸上實際地震資料反演結(jié)果

3 結(jié)論

為減小未知地震子波波形(或延遲時)及走時測量方法對(絕對)初至?xí)r間檢測帶來的誤差,本文對觀測數(shù)據(jù)和模擬數(shù)據(jù)采用相同的初至拾取方法并計算雙差走時,以降低甚至消除走時測量誤差對反演結(jié)果的影響。雖然根據(jù)公式推導(dǎo),嚴(yán)格的雙差走時應(yīng)該采用加權(quán)相位延遲計算得到,但在實際應(yīng)用中直接采用(自動或人工)拾取的初至計算雙差走時得到的結(jié)果仍然是穩(wěn)健的。

相較于傳統(tǒng)波動方程走時反演方法中引入的一階Born近似(要求速度異常體的尺度和擾動強度都足夠小),本文結(jié)合雙差走時測量及Rytov近似構(gòu)建相位延遲敏感度核函數(shù),可以更好地預(yù)測由于(大尺度)速度擾動引起的前向散射波的相位擾動,因此降低了對初始模型精度的要求。

在反問題數(shù)值求解方面,本文由于引入了基于隱式矩陣向量積的高斯-牛頓迭代算法,故僅需波動方程Born模擬及逆時偏移算法即可高效計算高斯-牛頓搜索方向,無需顯式計算和存儲Hessian矩陣,因此適用于求解大規(guī)模計算問題。

附錄A 隱式矩陣向量乘法

1) 矩陣向量積Kddp(p為模型空間中的向量p=p(x))。

根據(jù)(15)式,Kdd中任意一行與模型空間中的向量p的內(nèi)積可以表示為:

(Kddp)(xi,xj,xs)=〈K(x;xi,xs)-K(x;xj,xs),p〉M=〈K(x;xi,xs),p〉M-〈K(x;xj,xs),p〉M

(A1)

馮波等[25]證明,走時敏感度核函數(shù)與模型空間向量的內(nèi)積可以轉(zhuǎn)化為波場空間中的內(nèi)積:

(Kp)(xr,xs)=〈uq(xr,t;xs),up(xr,t;xs)〉T

(A2)

將(A2)代入(A1),有:

(Kddp)(xi,xj,xs)=〈uq(xi,t;xs),up(xi,t;xs)〉T-〈uq(xj,t;xs),up(xj,t;xs)〉T

(A3)

(A4)

其中λdd(x,T-t;xs)由如下雙差走時伴隨震源產(chǎn)生:

(A5)

順序計算(A3)及(A4)式,無需顯式計算和存儲敏感度核函數(shù)及Hessian矩陣,可以直接獲得Hessian向量積。