國運(yùn)東,孟凡冰,秦廣勝,李傳強(qiáng),李慶洋

(中國石油化工股份有限公司中原油田分公司物探研究院,河南濮陽 457001)

全波形反演方法(FWI)綜合利用地震數(shù)據(jù)的走時(shí)、相位、振幅等全波場信息,通過地震數(shù)據(jù)與模擬數(shù)據(jù)的差值建立目標(biāo)泛函,不斷更新地下參數(shù)模型,使得模擬地震數(shù)據(jù)越來接近觀測數(shù)據(jù),從而達(dá)到建立不同參數(shù)場的目的,FWI是一種高精度、高分辨率的物性參數(shù)建模方法。與常規(guī)的射線層析以及偏移速度分析(MVA)速度建模方法相比,FWI具有更高的分辨率以及速度建模精度,但是其在應(yīng)用中仍然存在許多問題,針對(duì)FWI不同問題的解決方法也相繼被提出[1-5]。

20世紀(jì)80年代,TARANTOLA[6]首先建立了全波形反演的理論框架,但是在缺少低頻和背景場不準(zhǔn)確的情況下,容易引起周波跳躍現(xiàn)象。BUNKS等[7]提出了分頻多尺度反演方法與網(wǎng)格多尺度反演方法,指出在大網(wǎng)格、低頻尺度下反演的背景速度模型,可以有效解決全波形反演的周波跳躍問題;SIRGUE等[8]給出了頻域多尺度FWI的頻率選取策略[9];BOONYASIRIWAT等[10]構(gòu)建了地震數(shù)據(jù)的維納濾波器,實(shí)現(xiàn)了時(shí)域FWI的多尺度方法[11]。除利用地震數(shù)據(jù)中的低頻組分反演背景速度模型外,MORA[12]和PRATT等[13]通過采用大偏移距以及透射波數(shù)據(jù)構(gòu)建背景速度模型,避免了全波形反演過程中的周波跳躍問題[14]。SHIPP等[15]通過對(duì)地震數(shù)據(jù)進(jìn)行分偏移距的反演,指出地震數(shù)據(jù)在大偏移距時(shí),由于傳播時(shí)間較長以及波形復(fù)雜性增加會(huì)加劇全波形反演的非線性。BAETEN等[16]和PLESSIX等[17]通過實(shí)際資料數(shù)據(jù)中的低頻、大偏移距組分反演試算,證實(shí)了低頻、大偏移距組分在實(shí)際地震數(shù)據(jù)反演中的重要性。但在實(shí)際地震數(shù)據(jù)采集中,想要獲取質(zhì)量非常高的低頻、大偏移距數(shù)據(jù)信息比較困難,因此許多學(xué)者發(fā)展了不同的低頻信息缺失情況下的背景速度模型構(gòu)建方法。BOZDAG等[18]利用由Hilbert變換得到的瞬時(shí)相位與瞬時(shí)振幅信息構(gòu)建新的波形反演目標(biāo)泛函,降低了全波形反演的非線性;隨后CHI等[19]提出基于地震包絡(luò)的目標(biāo)泛函,即利用波形的包絡(luò)殘差代替波形殘差,并推導(dǎo)相應(yīng)的梯度公式;WU等[20-21]通過對(duì)地震信號(hào)包絡(luò)的研究,揭示了地震信號(hào)的調(diào)制-解調(diào)機(jī)制,指出地震信號(hào)中反映地下大尺度特征的長波長分量是通過非線性的方式與短波長分量的信號(hào)調(diào)制在一起的,應(yīng)用包絡(luò)反演可以構(gòu)建這些低頻分量,從而降低全波形反演的非線性。

另一方面,波數(shù)連續(xù)性對(duì)于全波形反演非常重要[22],從理論上講,初始模型提供的低波數(shù)組分需要確保模型擾動(dòng)能夠反演在準(zhǔn)確位置,否則會(huì)發(fā)生周期跳躍[23]。對(duì)于特定的初始模型,FWI的梯度必須保證波數(shù)域中覆蓋范圍的連續(xù)性。從模型的角度來看,可以通過梯度來直接訪問特定波數(shù)分量[24],反射波全波形反演(RFWI)方法包含兩步工作流程,其中反復(fù)交替更新背景速度和高波數(shù)模型,從而能夠?yàn)镕WI和偏移成像恢復(fù)準(zhǔn)確的背景模型[23,25-27]?？赏ㄟ^對(duì)反射波進(jìn)行成像方法,適當(dāng)提取反射波場所攜帶的長波長信息,該成像可以計(jì)算對(duì)應(yīng)于高散射角的反射層位到地面的透射波路徑[28]。RFWI需要巨大的計(jì)算資源,包含F(xiàn)WI與最小二乘逆時(shí)偏移(LSRTM),一些學(xué)者使用其它方法來分離波場,如使用頻域局部散射角濾波器從平滑背景模型中恢復(fù)長波數(shù)速度擾動(dòng)[29-30],或者使用易于實(shí)現(xiàn)且有效的角度相關(guān)濾波技術(shù),實(shí)現(xiàn)在平面波域中分離低波數(shù)成分,用于反演背景速度模型[31]。GUO等[32]將平面波反演方法應(yīng)用于全波形反演的背景速度模型反演,為后續(xù)反演提供準(zhǔn)確的背景速度模型。此外,計(jì)算效率也是限制FWI應(yīng)用的關(guān)鍵因素之一,平面波編碼的反演方法是提高計(jì)算效率的有效方法之一,其計(jì)算量不再與炮數(shù)量呈線性關(guān)系,而與射線參數(shù)的數(shù)量直接相關(guān),通過正確選擇射線參數(shù),可以合理地抑制串?dāng)_噪聲,由于射線參數(shù)的數(shù)量遠(yuǎn)小于炮的數(shù)量,因此大大地提高了計(jì)算效率[33-35]。

在實(shí)際野外采集數(shù)據(jù)中缺失有效的低頻信號(hào),因此低頻數(shù)據(jù)重構(gòu)反演或不依賴低頻的低波數(shù)反演方法(例如反射波反演、包絡(luò)反演)成為近年來反演的研究熱點(diǎn),但是反射波反演也存在計(jì)算量大的問題,快速有效地利用低波數(shù)成分進(jìn)行初始速度模型反演具有重要意義。本文設(shè)計(jì)了一種適用于三維的動(dòng)態(tài)錐面波編碼的三維背景速度模型構(gòu)建方法,在較差的初始模型和不含低頻數(shù)據(jù)情況下,通過從低到高的序列中恢復(fù)波數(shù)分量,并且采用動(dòng)態(tài)參考點(diǎn)平面波策略僅需要少量的平面波,達(dá)到高效構(gòu)建背景速度模型的目的。

1 方法原理

1.1 基本原理

基于L2范數(shù)的全波形反演方法的目標(biāo)函數(shù)[6]表示為:

(1)

式中:uobs和ucal分別表示在t時(shí)刻的實(shí)際觀測地震記錄與正演模擬記錄;震源點(diǎn)為xs;接收點(diǎn)為xr;m為速度參數(shù)模型。在波形反演過程中,需要多次迭代進(jìn)行參數(shù)更新,因此計(jì)算量巨大,采用平面波編碼反演是一種提高效率的方法。

在三維速度模型反演中,平面波需要2個(gè)射線參數(shù)px與py,然后通過時(shí)間延遲疊加獲得平面波[33,35]:

Δtjk(p)=px(xsj-x0)+py(ysk-y0)

(2)

式中:(xsj,ysk)表示沿x測線第j炮與沿y測線第k炮的炮點(diǎn)位置;(x0,y0)表示合成平面波在地表的原始參考位置;x0與y0分別表示三維空間上的x軸(主側(cè)線方向)位置與y軸(聯(lián)絡(luò)側(cè)線方向)位置。應(yīng)用平面波編碼,通過改變射線參數(shù)p的動(dòng)態(tài)平面波編碼的全波形反演目標(biāo)函數(shù)可定義為[33]:

(3)

本文提出的動(dòng)態(tài)平面波編碼由射線參數(shù)p控制波場的角度,波形反演的散射角信息常常與射線參數(shù)p相關(guān),而在三維地震數(shù)據(jù)的波形反演中,地震波數(shù)的高低與散射角的大小密切相關(guān),散射角的大小與射線參數(shù)的絕對(duì)值大小密切相關(guān),常規(guī)平面波FWI含有兩個(gè)方向的射線參數(shù),因此很難控制與特定散射角大小的射線參數(shù),為了在相同的射線參數(shù)p下(控制某一散射角)得到多個(gè)梯度,本文選取合適的射線參數(shù)和參考點(diǎn)構(gòu)造特殊的平面波,即錐面波,并用錐面波對(duì)三維地震炮記錄編碼,得到錐面波超級(jí)炮記錄,其中錐面波的生成函數(shù)為:

Δt[(xk,yk):p,(x0,y0)]=

(4)

式中:Δt表示時(shí)間延遲量;p為射線參數(shù),其值越大表示錐面波的傳播方向與地面夾角越大,與傳統(tǒng)的平面波編碼相比,其不再分兩個(gè)方向,并且與地震波場的入射角相關(guān);第k炮位置為需要編碼的炮點(diǎn)位置,原始參考位置(x0,y0)為每次迭代隨機(jī)選取的位置。

合成的錐面波超級(jí)炮記錄可以表示為:

(5)

利用(5)式中的錐面波函數(shù)對(duì)震源子波編碼,進(jìn)行三維有限差分正演得到錐面波炮記錄,錐面波對(duì)震源子波編碼具體如下:

(6)

三維錐面波速度模型反演的目標(biāo)函數(shù)表示為:

(7)

根據(jù)數(shù)據(jù)殘差利用最速下降法或者共軛梯度法確定梯度,由于本文采用動(dòng)態(tài)編碼的方式,兩次反演迭代之間的梯度場相差較大,因此采用最速下降法進(jìn)行梯度更新。使用伴隨狀態(tài)方法[36],目標(biāo)函數(shù)相對(duì)于模型參數(shù)的梯度可以通過以下方法獲得:

(8)

1.2 反演流程

基于動(dòng)態(tài)錐面波編碼的三維背景速度模型構(gòu)建方法的流程如圖1所示。相比于傳統(tǒng)的全波形反演,本文方法每次反演都需要選取不同的射線參數(shù),在同一射線參數(shù)下,利用不同參考點(diǎn),對(duì)炮記錄和震源進(jìn)行相同的動(dòng)態(tài)錐面波編碼,滿足對(duì)這一射線參數(shù)的反演精度,再逐步減小射線參數(shù),達(dá)到背景速度模型的反演目標(biāo),利用得到的背景速度模型進(jìn)行常規(guī)全波形反演得到最終的速度模型。

圖1 動(dòng)態(tài)錐面波編碼的全波形反演流程

2 模型測試

2.1 簡單模型測試

簡單模型如圖2a所示。模型主要有3層介質(zhì),其中第二層中心部分存在一個(gè)洼陷,其速度大小與第一層相同,模型尺寸(x,y和z方向)為1500m×1000m×1300m,網(wǎng)格間距為10m,模型速度范圍為3000～4000m/s。

計(jì)算參數(shù):采用主頻15Hz的雷克子波,采樣時(shí)間長度為1.2s,采樣間隔0.5ms,采用時(shí)間2階、空間8階的有限差分交錯(cuò)網(wǎng)格進(jìn)行正演模擬,施加完全匹配吸收(PML)邊界條件。觀測系統(tǒng):共有176炮(16×11),炮間隔10m,第一炮位置為(10m,0,0),地表全接收。采用線性梯度模型(圖2b)作為反演的初始模型,反演過程中假定震源子波已知。

圖3為不同三維地震炮記錄以及數(shù)據(jù)殘差。利用GPU進(jìn)行3D正演模擬,使用10Hz主頻的雷克子波作為目標(biāo)子波,進(jìn)行維納濾波后的單炮記錄如圖3a 所示,炮點(diǎn)位于圖2a的(500m,500m,10m)處,圖上能清晰識(shí)別直達(dá)波與反射波等有效信息;圖3b 為基于公式(6)合成的震源進(jìn)行正演模擬的炮記錄;圖3c為根據(jù)公式(5)合成的錐面波超級(jí)炮記錄,可以發(fā)現(xiàn),由于不同單炮在時(shí)間和空間上的重疊,模擬波場較為復(fù)雜,表層能觀測到能量較強(qiáng)的直達(dá)波場;圖3d為基于公式(5)合成的錐面波超級(jí)炮記錄與基于公式(6)合成的震源進(jìn)行正演模擬的炮記錄的殘差,可以看出,直達(dá)波能量削弱,模擬波場更為復(fù)雜。將兩種超級(jí)炮記錄分別作為輸入,為后續(xù)全波形反演提供基礎(chǔ)數(shù)據(jù)。

圖2 反演采用的速度模型(洼陷模型)a 真實(shí)模型; b 梯度場的初始模型

圖3 不同三維地震炮記錄a 單炮記錄; b 利用合成震源正演模擬的炮記錄; c 根據(jù)公式(5)合成的錐面波超級(jí)炮記錄; d 觀測數(shù)據(jù)與模擬數(shù)據(jù)的殘差

為了測試本文方法的有效性,本節(jié)采用維納濾波多尺度的FWI方法進(jìn)行對(duì)比,由于三維的計(jì)算量巨大,采用隨機(jī)極性編碼的方法將所有合成一個(gè)超級(jí)炮記錄進(jìn)行反演,不同頻率子波維納濾波后,炮記錄的頻率發(fā)生對(duì)應(yīng)變化,濾波后的地震記錄與對(duì)應(yīng)主頻子波產(chǎn)生的正演波場進(jìn)行匹配達(dá)到頻率多尺度反演的目的,當(dāng)然極性編碼也存在串?dāng)_噪聲的影響,但是本節(jié)主要測試背景場的構(gòu)建能力,因此暫不考慮串?dāng)_的影響,并且將極性編碼的多震源多尺度FWI記為MFWI(Multi-source encoding multi-scale FWI)方法。

常規(guī)多震源編碼FWI經(jīng)過40次迭代后反演的速度模型如圖4所示(截面:x=750m,y=750m,z=600m)。由于采用的是梯度場作為輸入的初始速度模型,與真實(shí)速度模型存在較大的差異,且地震數(shù)據(jù)的主頻比較高,波形反演過程中陷入局部極小值,使得洼陷構(gòu)造的區(qū)域反演不清晰,無法有效地更新速度模型,洼陷結(jié)構(gòu)很難恢復(fù),只能看到洼陷的大體形態(tài)。由于上層反演不準(zhǔn)確,洼陷下面層位的反演也不準(zhǔn)確,變得彎曲。

圖4 常規(guī)多震源編碼FWI三維洼陷模型反演結(jié)果(10～15Hz)

然后利用錐面波多尺度全波形反演方法(CFWI)進(jìn)行測試,結(jié)果如圖5所示。測試中,平面波入射方向與z方向的夾角的正弦值從0.8減少到0.4,當(dāng)數(shù)據(jù)誤差收斂到一定數(shù)值,變化較小時(shí),對(duì)同一射線參數(shù)p停止反演。測試中僅合成了1個(gè)錐面波,并且參考震源點(diǎn)隨機(jī)地分布在表面上,每次迭代都改變。使用錐面波FWI(CFWI)進(jìn)行反演結(jié)果如圖5b所示,與利用MFWI反演(主頻5～10Hz)結(jié)果(圖5a)對(duì)比,可以觀察到:①通過本文CFWI方法反演的背景速度模型與含有低頻反演結(jié)果相似。②反演的洼陷構(gòu)造的區(qū)域大體形態(tài)比較準(zhǔn)確,并且洼陷構(gòu)造下面的水平地層歸位也相對(duì)準(zhǔn)確。本文對(duì)單次梯度場的反演時(shí)間也進(jìn)行了對(duì)比。采用相同的配置,CFWI方法需要8.6min,MFWI方法需要8.5min,傳統(tǒng)的單炮疊加方法需要1393.2min,CFWI方法與MFWI方法的計(jì)算效率基本相近。由于本文采用了動(dòng)態(tài)編碼的CFWI,僅需要合成一個(gè)超級(jí)炮記錄,就能達(dá)到反演背景速度模型的目的,因此大大提高了計(jì)算效率。

圖5 三維洼陷模型采用不同多震源編碼多尺度FWI的反演結(jié)果a 常規(guī)多震源編碼多尺度FWI(MFWI)(主頻5～10Hz); b 基于錐面波的CFWI(主頻10Hz)

利用CFWI反演的背景速度模型為后續(xù)的MFWI提供了較為準(zhǔn)確的背景速度模型。這樣,使用相同子波的常規(guī)MFWI通過增加高波數(shù)分量進(jìn)一步提高速度模型的反演精度(圖6a),基本可以達(dá)到全頻帶多尺度的MFWI的反演效果(圖6b)。

圖6 不同背景場反演的三維洼陷模型結(jié)果a CFWI+MFWI; b 常規(guī)全頻帶多尺度FWI(MFWI)(主頻5～15Hz)

2.2 復(fù)雜模型測試

3D逆沖模型(Overthrust模型)如圖7a所示,由于模型過大,本文抽稀了一部分進(jìn)行3D多震源多尺度FWI(MFWI)測試。測試中,深度方向采樣點(diǎn)為100,Inline方向(x方向)采樣點(diǎn)為201,Crossline方向(y方向)采樣點(diǎn)為101,采樣間隔為10m,模型速度范圍為2428～6000m/s。該模型存在一個(gè)Inline方向的大型逆掩斷層構(gòu)造以及一系列的斷塊,縱向界面上發(fā)育有多條河道,整個(gè)模型的地質(zhì)構(gòu)造比較復(fù)雜。

計(jì)算參數(shù):以主頻15Hz雷克子波作為震源子波,地震記錄的時(shí)長為1.2s,采樣間隔0.5ms,采用時(shí)間2階、空間8階的有限差分交錯(cuò)網(wǎng)格進(jìn)行正演模擬,施加PML邊界條件。觀測系統(tǒng):共有231炮(21×11),炮間隔10m,第一炮位置為(10m,0,0),地表全接收,共20301個(gè)檢波點(diǎn),檢波點(diǎn)間隔為10m。采用線性梯度模型作為反演的初始模型(圖7b),反演過程中給定震源子波。

圖7 三維Overthrust速度模型(a)和梯度場的初始速度模型(b)

本文以含有低頻的MFWI結(jié)果作為標(biāo)準(zhǔn)進(jìn)行對(duì)比。反演結(jié)果如圖8所示。圖8b中左側(cè)部分?jǐn)鄬右约昂拥赖葟?fù)雜構(gòu)造,復(fù)雜區(qū)域的背景場反演比較準(zhǔn)確,河道也刻畫得比較清楚,驗(yàn)證了本文CFWI反演策略的有效性。應(yīng)用不同的背景速度模型進(jìn)行多尺度FWI反演,結(jié)果如圖9所示。常規(guī)MFWI經(jīng)過40次迭代后反演的速度模型如圖9a所示。由于梯度場的初始模型與真實(shí)模型差異較大,并且數(shù)據(jù)的主頻比較高,因此反演陷入局部極小值。河道的結(jié)構(gòu)很難恢復(fù),這時(shí)逆掩斷層區(qū)域會(huì)出現(xiàn)強(qiáng)假象,深部的層位也位于不正確的位置,并且斷層被強(qiáng)烈的低噪聲污染。與初始模型的反演結(jié)果相比(圖9a),應(yīng)用CFWI+MFWI獲得的最終速度模型比較準(zhǔn)確(圖9c),圖上包括了斷層和河道在內(nèi)的大部分主要特征。主要斷層和河道結(jié)構(gòu)得到更好的反演,反射層位已經(jīng)接近真實(shí)位置,動(dòng)態(tài)編碼壓制串?dāng)_噪聲效果較好,幾乎與全頻帶多尺度MFWI結(jié)果相近(圖9b)。相對(duì)傳統(tǒng)三維單炮疊加的FWI方法,利用動(dòng)態(tài)CFWI+MFWI,需要合成的超級(jí)炮數(shù)量大大減少,顯著提高了計(jì)算效率。

圖8 三維逆沖模型背景場反演結(jié)果a 常規(guī)多震源編碼多尺度FWI(MFWI)(5～10Hz); b 基于錐面波的CFWI(主頻:10Hz)

圖9 三維逆沖模型最終反演結(jié)果a MFWI的反演結(jié)果(10～15Hz); b MFWI的反演結(jié)果(5～15Hz); c CFWI+MFWI的反演結(jié)果

3 結(jié)論及討論

本文針對(duì)波形反演中地震數(shù)據(jù)缺少低頻的情況下背景速度模型反演較為困難的問題,從低波數(shù)的反演角度出發(fā),應(yīng)用新的平面波合成方法,發(fā)展了一種基于平面波的有效低波數(shù)背景速度模型構(gòu)建方法,通過模型試算和對(duì)比分析,得到如下幾點(diǎn)認(rèn)識(shí)。

1) 動(dòng)態(tài)平面波多尺度全波形反演方法可以對(duì)單個(gè)射線參數(shù)合成更多的平面波,從而反演得到背景速度模型。在每次迭代選用少量不同的平面波也可以反演出較好的背景速度場,大大減小了計(jì)算量。

2) 將動(dòng)態(tài)平面波方法與傳統(tǒng)的FWI相結(jié)合,可以實(shí)現(xiàn)從低波數(shù)到高波數(shù)的速度模型反演,即使在缺少低頻的情況下,也可以實(shí)現(xiàn)較高精度的速度模型反演,具有較好的收斂性和穩(wěn)定性。

3) 針對(duì)三維觀測系統(tǒng),本文提出的錐面波編碼多震源多尺度全波形反演方法(CFWI+MFWI),通過控制入射角的方式可以完成三維背景速度模型的構(gòu)建。由于錐面波超級(jí)炮的個(gè)數(shù)遠(yuǎn)小于原始的普通炮集個(gè)數(shù),因此該方法大大降低了計(jì)算量。

雖然本文提出的一種利用平面波的直接反演方法,通過直接利用低波數(shù)反演背景速度模型,避免了直接進(jìn)行波數(shù)分離的過程,并且相對(duì)于反射波波形反演,避免了最小二乘的反演過程,但還存在觀測系統(tǒng)的要求嚴(yán)格問題。此外在背景場準(zhǔn)確情況下,多尺度反演策略在后期準(zhǔn)確速度模型反演仍具有重要的意義,而且多震源方法存在固有的串?dāng)_噪聲對(duì)中高波數(shù)反演的影響也比較明顯,這也是后續(xù)要研究的主要方向。