趙 強(qiáng),姜大建,朱成宏,魏哲楓

(1.頁巖油氣富集機(jī)理與有效開發(fā)國家重點(diǎn)實(shí)驗(yàn)室,北京100083;2.中國石化彈性波理論與探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京100083;3.中國石化石油勘探開發(fā)研究院,北京100083)

全波形反演(FWI)充分利用地震記錄中的動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)全波場(chǎng)信息,具有比常規(guī)射線層析類反演方法更高的建模精度[1-3]。目前,FWI已成為提高地震勘探精度的重要手段之一?？紤]到反演問題的欠定性以及全波形反演的非線性特征,FWI誤差泛函存在許多局部極值點(diǎn),所以在實(shí)際應(yīng)用時(shí)往往需要低中波數(shù)成分來降低反演的非線性和多解性[2-4]。

FWI可以綜合利用不同的波場(chǎng)信息(相位、振幅、包絡(luò)和波形等)構(gòu)建目標(biāo)反演,從而可實(shí)現(xiàn)多尺度、多特征反演,降低對(duì)于初始模型的依賴性[5-8]。其中,利用地震記錄中的反射波進(jìn)行反射波形反演(RWI)是挖掘中深部長波長背景信息的理想方式之一[9-11]。與FWI不同,為了增加平滑速度場(chǎng)情況下的反射信息,RWI利用反偏移獲得模擬地震數(shù)據(jù),伴隨源反傳播進(jìn)行梯度更新,可有效恢復(fù)長波長背景信息。按照利用信息的不同,RWI的目標(biāo)泛函一般包括波形和走時(shí)[12]。通常意義下,背景速度對(duì)地震波走時(shí)的影響大于振幅,因此地震資料中的走時(shí)誤差對(duì)背景速度的敏感程度大于振幅信息[13]。另一方面,考慮到走時(shí)信息的大尺度背景,利用走時(shí)信息構(gòu)建的反射波走時(shí)反演(RTI)的非線性程度更低[14-18],更有利于避免因初始模型不準(zhǔn)確造成的周波跳躍問題。

對(duì)于RTI而言,走時(shí)提取的準(zhǔn)確與否直接關(guān)系到反演梯度的更新方向。目前,提取走時(shí)信息的方法主要有圖像處理領(lǐng)域的動(dòng)態(tài)圖像估計(jì)(DIW)[19-20]或時(shí)延互相關(guān)[14]。時(shí)延互相關(guān)對(duì)相同空間位置單道進(jìn)行互相關(guān)計(jì)算,從而實(shí)現(xiàn)整道估計(jì)。單道全局估計(jì)難以反映走時(shí)延遲時(shí)正時(shí)負(fù)的非穩(wěn)態(tài)特征,反演結(jié)果容易出現(xiàn)周波跳動(dòng)。DIW被認(rèn)為是一種非剛性的圖像標(biāo)記方法,可提取信號(hào)的連續(xù)性變化。時(shí)差求取過程中,DIW不需要指定窗口大小,同時(shí)可以限制水平和垂直方向的變化幅度。在存在隨機(jī)噪聲的情況下,DIW通常也可以獲得較好的估計(jì)結(jié)果[19]。然而,單道估計(jì)容易造成橫向連續(xù)性較差的問題,DIW雖然可以在橫向上添加軟約束,但仍容易出現(xiàn)局部時(shí)差估計(jì)“斑狀”異常。利用多道間互相約束是解決該問題的有效方案之一,低秩近似在沿空間方向的信號(hào)增強(qiáng)中顯示出其潛力[16-17]。該類方法假設(shè)地震信號(hào)在局部區(qū)域具有相似性,存在信息冗余?；跇?gòu)造特殊的矩陣或者張量形式,可得到有效信號(hào)的低秩近似結(jié)果,并同時(shí)提高橫向連續(xù)性[21-25]。低秩近似良好的橫向約束為增加時(shí)差估計(jì)結(jié)果的道間連續(xù)性提供了更好的選擇。

針對(duì)傳統(tǒng)全局時(shí)延互相關(guān)存在時(shí)間方向非穩(wěn)態(tài)性差,以及DIW求取結(jié)果在空間方向連續(xù)性不足的問題,提出低秩動(dòng)態(tài)匹配約束下的反射波走時(shí)反演策略。具體地,RTI伴隨源的非穩(wěn)態(tài)時(shí)移量仍利用DIW估計(jì),同時(shí)構(gòu)造沿頻率切片數(shù)據(jù)中的低秩整形算子,對(duì)DIW道間時(shí)移量進(jìn)行橫向硬約束,緩解空間連續(xù)性差的問題;數(shù)值算例部分,利用Marmousi模型對(duì)方法的時(shí)移估計(jì)、伴隨源、反演及成像效果進(jìn)行了細(xì)致的分析討論。

1 技術(shù)方法

1.1 基于反偏移技術(shù)的反射波形反演

與全波形反演類似,反射波形反演也是在預(yù)定義的目標(biāo)函數(shù)情況下求解最優(yōu)化模型,以期實(shí)現(xiàn)模擬數(shù)據(jù)dsyn(x,t)與實(shí)測(cè)數(shù)據(jù)dobs(x,t)之間信息差異的極小化目標(biāo)泛函[26]：

(1)

式中：v和r分別是背景速度(低波數(shù)分量)和反射界面信息(高波數(shù)分量);uS和uR分別是入射波場(chǎng)及其對(duì)應(yīng)反偏移波場(chǎng);dsyn(x,t)=PuR(x,t),P是反偏移波場(chǎng)uR在檢波點(diǎn)位置處的波場(chǎng)提取函數(shù);A是波場(chǎng)傳播算子;f是預(yù)定義的目標(biāo)泛函,用以衡量模擬數(shù)據(jù)dsyn(x,t)與實(shí)測(cè)數(shù)據(jù)dobs(x,t)間的誤差。

按照(1)式所示的目標(biāo)函數(shù),反射界面及背景速度的更新梯度分別為[12]：

式中：λR和λS分別是反傳入射波場(chǎng)及其對(duì)應(yīng)反偏移波場(chǎng)。

圖1為(2)式和(3)式互相關(guān)分量示意。由(2)式及圖1a可以看出,高波數(shù)反射界面信息可通過入射波場(chǎng)的一階時(shí)間偏導(dǎo)與檢波反傳波場(chǎng)的互相關(guān)計(jì)算得到,其計(jì)算過程與傳統(tǒng)全波形反演的計(jì)算過程類似,梯度更新量主要反映界面信息;背景速度場(chǎng)梯度由兩項(xiàng)求和得到,其中,(3)式右側(cè)第1項(xiàng),即如圖1b所示結(jié)果,由震源入射波場(chǎng)與檢波側(cè)伴隨源反偏移波場(chǎng)互相關(guān)得到。圖1c所示第2項(xiàng)由震源側(cè)反偏移波場(chǎng)與檢波側(cè)伴隨源反傳波場(chǎng)互相關(guān)得到。對(duì)比圖1a、圖1b和圖1c可以看出,(2)式主要利用小角度波場(chǎng),(3)式主要利用大角度波場(chǎng)。根據(jù)MILLER等[27]和JIN等[28]給出的局部波數(shù)公式,即

圖1 公式(2)和公式(3)互相關(guān)分量示意a (2)式右側(cè)波場(chǎng)互相關(guān)計(jì)算示意; b (3)式右側(cè)第1項(xiàng)震源側(cè)反偏移互相關(guān)示意; c (3)式右側(cè)第2項(xiàng)檢波側(cè)反偏移互相關(guān)示意

式中：θ為反射角,n為單位向量,可知低頻和大反射角信息是低波數(shù)成分的直接來源,也間接證明反射波波形反演對(duì)于低波數(shù)背景速度反演的有效性。

1.2 動(dòng)態(tài)匹配約束下的反射波走時(shí)反演

RWI一般假設(shè)(2)式所求取的成像剖面僅存在由于背景速度誤差所導(dǎo)致的深度差異,振幅、頻率和相位信息需保持相對(duì)正確。因此,反射波形反演求解過程中往往引入最小二乘法求解(2)式[12],這進(jìn)一步增加了反射波形反演的計(jì)算量。另一方面,背景速度對(duì)于地震波走時(shí)的影響程度強(qiáng)于振幅信息,因此,走時(shí)信息對(duì)于背景速度誤差更加敏感。鑒于此,反射波走時(shí)反演成為獲取背景速度的有效手段。與波形反演不同,走時(shí)反演采用走時(shí)相關(guān)目標(biāo)泛函度量反演誤差[12],則有：

(5)

LUO等[14]對(duì)(5)式中的走時(shí)目標(biāo)泛函敏感核進(jìn)行了詳細(xì)推導(dǎo),本文不作贅述。根據(jù)LUO等[14]給出的推導(dǎo)結(jié)果,檢波波場(chǎng)伴隨源可表示為：

從(6)式可以看出,與RWI相比,RTI的伴隨源主要取決于時(shí)移量τ的求取精度,是在迭代過程中產(chǎn)生正確梯度方向的關(guān)鍵點(diǎn)。而震源、檢波波場(chǎng)構(gòu)造形式以及梯度量的更新則與RWI基本無差別。然而,傳統(tǒng)意義上單道互相關(guān)只能測(cè)量全局相似性或選定窗口內(nèi)的局部相似性。當(dāng)時(shí)窗內(nèi)時(shí)移量快速變化時(shí),可能會(huì)出現(xiàn)周波跳動(dòng)的問題。動(dòng)態(tài)圖像估計(jì)(DIW)是對(duì)語音識(shí)別的經(jīng)典動(dòng)態(tài)匹配算法的有效擴(kuò)展。按照MA等[20]給出的定義,基于DIW的時(shí)移量求取可表示為如下的最優(yōu)化問題,即

式中：σt,σR,σS是時(shí)延量求解在時(shí)間方向、檢波點(diǎn)方向以及炮點(diǎn)方向的軟約束條件,目的是避免異常值的出現(xiàn)。

1.3 低秩動(dòng)態(tài)匹配約束時(shí)移求解

當(dāng)初始模型準(zhǔn)確度較低時(shí),成像界面同相軸的連續(xù)性及信噪比勢(shì)必受到影響,隨之造成反偏移地震記錄中存在明顯的散射波。此時(shí)即便DIW在求取過程中存在3個(gè)維度的約束條件,其時(shí)移估計(jì)結(jié)果的橫向連續(xù)性同樣難以保證,容易造成走時(shí)剖面的斑狀特征[20],時(shí)差求取結(jié)果穩(wěn)定性變差,容易在梯度更新時(shí)陷入局部極值。

一般情況下,地震記錄中有效反射波存在局部線性平面波特征,在數(shù)學(xué)意義上是低秩的。為了克服水平不連續(xù)性,本文對(duì)DIW求取的時(shí)移剖面引入低秩整形算子約束水平方向的時(shí)移量。假設(shè)時(shí)移τ(x,t)由一系列斜率為p1,p2,…,pM的線性同相軸求取得到,其頻率域形式S(x,ω)為：

式中：A(ω)是圓頻率,為ω信號(hào)的振幅信息。

令θm=e-jωpmΔx,則圓頻率為ω的信號(hào)在位置Δx,2Δx,…,NΔx可表示為：

以(9)式單頻信號(hào)為基礎(chǔ),選擇參數(shù)L=floor(N/2)+1和K=N-L+1為行數(shù)和列數(shù),構(gòu)造如下形式的Hankel或Toeplitz矩陣。

HL[S(x,ω)]=A(ω)·

對(duì)(10)式進(jìn)行奇異值分解(SVD),則有：

HL[S(x,ω)]=A(ω)·

由(11)式可以看出,當(dāng)?shù)卣鹦盘?hào)由一系列線性同相軸組成時(shí),地震信號(hào)等價(jià)于一系列單斜率地震同相軸的疊加,且該矩陣的秩由斜率pm的個(gè)數(shù)決定,而與采集道數(shù)和采樣時(shí)間無關(guān)。低秩硬約束可表示為核函數(shù)的求解形式‖H[S(x,ω)]‖*,可通過帶限奇異值分解求解。此時(shí)低秩約束小的動(dòng)態(tài)時(shí)移求解可表示為：

如圖2a和圖2b所示,不準(zhǔn)確的初始模型會(huì)形成散射波,而在時(shí)差估計(jì)中使用該反射波會(huì)產(chǎn)生短波長擾動(dòng),影響DIW結(jié)果的準(zhǔn)確性(圖2c)。圖2d為低秩約束下的時(shí)差估計(jì)結(jié)果。對(duì)比圖2c和圖2d可以看出,由初始模型不準(zhǔn)確引起的時(shí)差剖面水平不連續(xù)得到了有效緩解,更符合走時(shí)的大尺度特征。

圖2 觀測(cè)記錄與反偏移記錄時(shí)差剖面對(duì)比a 觀測(cè)記錄; b 使用線性初始模型的反偏移記錄; c 使用DIW估計(jì)的時(shí)差; d 低秩約束下的時(shí)差估計(jì)結(jié)果

2 應(yīng)用實(shí)例

利用如圖3a所示的Marmousi重采樣模型進(jìn)行測(cè)試分析,線性初始模型如圖3b所示,模型大小為4650m×1480m。線性模型速度值與準(zhǔn)確模型速度值大小相對(duì)關(guān)系如圖3c所示。其中,數(shù)值-1表示初始模型速度偏小,數(shù)值1表示初始模型速度偏大,數(shù)值0表示初始模型速度準(zhǔn)確。速度反演震源采用主頻為10Hz的雷克子波,偏移成像采用主頻為20Hz的雷克子波,震源位于地表下10m位置,檢波點(diǎn)位于地表。炮點(diǎn)間距為100m,檢波點(diǎn)間距為10m。

圖3 Marmousi模型測(cè)試示意a 重采樣Marmousi模型; b 線性初始模型; c 初始模型與準(zhǔn)確模型速度大小相對(duì)關(guān)系

雙邊放炮,最大偏移距為1500m,與最大深度相當(dāng)。采用10階有限差分法模擬,縱橫向空間采樣間隔為10m,時(shí)間采樣間隔為1ms,模擬時(shí)長為1.6s。為了壓制邊界反射,在模型四周均設(shè)置20層吸收邊界。

為驗(yàn)證低秩動(dòng)態(tài)匹配約束的有效性,抽取走時(shí)反演過程中的一炮進(jìn)行分析,結(jié)果如圖4所示。對(duì)比圖4a 與圖4b可以看出,觀測(cè)地震記錄與反偏移得到的地震記錄相似度極高,表明(2)式得到的反偏移結(jié)果的相位信息相對(duì)準(zhǔn)確,為后續(xù)走時(shí)匹配的準(zhǔn)確性提供較好數(shù)據(jù)基礎(chǔ)。圖4c和圖4d分別為常規(guī)DIW及低秩DIW得到的時(shí)差剖面?？梢钥闯?即使常規(guī)DIW得到的時(shí)差剖面存在橫向道間約束,同一同相軸、不同道間的“掛面條”現(xiàn)象較明顯。低秩DIW得到的時(shí)差剖面的橫向連續(xù)性得到極大增強(qiáng)。圖4e所示的常規(guī)DIW方法得到的伴隨源同相軸存在明顯的錯(cuò)斷現(xiàn)象。圖4f 所示的低秩DIW方法得到的伴隨源的連續(xù)性明顯增加,更有利于后續(xù)反傳波場(chǎng)的穩(wěn)定性。與此同時(shí),從圖4a和圖4b的比較中還可以看出,在真實(shí)模型和初始模型得到的地震記錄中,同一界面的遠(yuǎn)偏移距反射波的時(shí)差已達(dá)到1/2波長,在這種情況下,本文方法仍可以較好地避免周期跳躍的問題,進(jìn)一步增加反射波走時(shí)反演的穩(wěn)定性。

圖4 不同方法計(jì)算得到的時(shí)差及伴隨源剖面a 觀測(cè)數(shù)據(jù); b 反偏移數(shù)據(jù); c 常規(guī)DIW時(shí)差剖面; d 低秩DIW時(shí)差剖面; e 常規(guī)DIW得到的伴隨源; f 低秩DIW得到的伴隨源

圖5為不同時(shí)差估計(jì)方法得到的走時(shí)反演梯度及迭代20次反演結(jié)果。對(duì)比圖5a和圖5b可以看出,不同時(shí)差估計(jì)結(jié)果得到的反演梯度方向與圖3c基本一致,即模型左側(cè)速度向減小反向變化,模型右側(cè)速度向增大方向變化。然而,受圖4e伴隨源估計(jì)結(jié)果中橫向連續(xù)性差的影響,圖5a梯度剖面中的局部異常點(diǎn)增多,圖5b所示反演結(jié)果中的高陡異常構(gòu)造增多。當(dāng)在伴隨源時(shí)差估計(jì)過程中引入低秩約束后,圖5c所示梯度平滑度大大提高,黑框區(qū)域更加符合走時(shí)反演的大尺度特性,對(duì)應(yīng)圖5d反演結(jié)果的連續(xù)性明顯提高。

圖6為不同初始模型的常規(guī)全波形反演結(jié)果。圖6a所示的反演結(jié)果表明單尺度全波形反演難以有效更新背景低波數(shù)速度,圖6b為利用圖5d作為初始模型的反演結(jié)果,反演結(jié)果改善明顯,與圖3a所示準(zhǔn)確模型的深度一致性更好。圖7a和圖7b分別展示了z方向400m和x方向1000m位置處準(zhǔn)確模型、初始模型、直接FWI結(jié)果以及低秩RTI+FWI反演結(jié)果?？梢钥闯?低秩RTI+FWI反演結(jié)果的整體趨勢(shì)與準(zhǔn)確模型一致性高,尤其是淺層速度吻合度高。圖8為不同速度模型的逆時(shí)偏移成像結(jié)果。與圖8a利用準(zhǔn)確速度模型得到的成像結(jié)果相比,圖8b初始線性模型成像結(jié)果存在明顯的深度差,具體體現(xiàn)在模型左側(cè)速度偏大區(qū)域成像深度“下拉”偏深、模型右側(cè)速度偏小區(qū)域成像深度“上翹”偏淺。即使利用單尺度全波形反演模型后,圖8c所示的成像結(jié)果的改善較小,說明背景走時(shí)信息未得到有效更新。利用圖5d作為初始模型的FWI反演結(jié)果,圖8d所示的對(duì)應(yīng)逆時(shí)偏移成像結(jié)果深度誤差明顯減小,左側(cè)“下拉”、右側(cè)“上翹”現(xiàn)象得到極大修正,說明本文方法可有效更新背景低波數(shù)成分,后續(xù)全波形反演結(jié)果更加準(zhǔn)確。

圖6 不同初始模型的常規(guī)全波形反演結(jié)果 a 線性初始模型的常規(guī)FWI反演結(jié)果; b 以圖5d作為初始模型的FWI反演結(jié)果

圖8 不同速度模型的逆時(shí)偏移成像結(jié)果a 準(zhǔn)確模型的逆時(shí)偏移成像結(jié)果; b 線性初始模型的逆時(shí)偏移成像結(jié)果; c 圖6a模型的逆時(shí)偏移成像結(jié)果; d 圖6b模型的逆時(shí)偏移成像結(jié)果

3 結(jié)論

長波長低波數(shù)成分在波形反演中起著至關(guān)重要的作用,同時(shí)也影響到高波數(shù)成分反演的穩(wěn)定性。針對(duì)反射波走時(shí)反演中的時(shí)差構(gòu)建環(huán)節(jié),構(gòu)建了一種低秩動(dòng)態(tài)匹配約束的時(shí)差提取方法,進(jìn)而提高了走時(shí)反演伴隨源求取和反演的精度。數(shù)值模擬結(jié)果表明,所提方法得到的時(shí)差剖面具有更好的穩(wěn)定性,反演梯度橫向連續(xù)性更好,更符合地質(zhì)沉積的結(jié)構(gòu)特征,反演結(jié)果連續(xù)性更優(yōu)。利用本文結(jié)果作為初始模型的全波形反演及逆時(shí)偏移成像結(jié)果均證明了其有效性。在三維及彈性波情況下的拓展應(yīng)用是后續(xù)重點(diǎn)研究方向。