基于雙對角通量校正的多尺度全波形反演

楊瑞冬 黃建平* 楊振杰 慕鑫茹 王兆忠 楊珊珊

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②山東省第一地質(zhì)礦產(chǎn)勘查院，山東濟南 250100； ③中國石油塔里木油田分公司，新疆庫爾勒 841000)

0 引言

隨著地震勘探復(fù)雜程度的日益增加，對高精度速度建模提出了更高要求。20世紀80年代，Taran-tola[1]提出時間域全波形反演(Full Waveform Inversion，F(xiàn)WI)方法，即通過多次迭代更新，不斷降低觀測記錄與預(yù)測記錄的殘差，實現(xiàn)對地下速度結(jié)構(gòu)的精確反演。FWI充分利用豐富的波形信息，因此它相較于其他方法具有更高的速度建模精度。但FWI存在高度病態(tài)的非線性特性[2-3]，當預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)對應(yīng)的反射軸時差大于半個周期時，就會出現(xiàn)“周波跳躍”現(xiàn)象，此時直接求取殘差會使反演陷入“局部極值”。精確的初始速度模型能有效緩解“周波跳躍”現(xiàn)象。因此，初始模型精度是FWI的關(guān)鍵。而實際應(yīng)用中，初始模型往往無法達到較高精度。為了降低FWI對初始模型的依賴性，許多學(xué)者做了大量的研究工作。低頻信號可反演宏觀速度信息，高頻信號刻畫速度細節(jié)，從低頻到高頻的反演策略能降低FWI對初始速度信息的依賴性。Mora[4]研究表明，大炮檢距地震記錄中含有豐富的低頻信息，有效緩解了“周波跳躍”，從而獲得較高精度的反演結(jié)果。Bunks等[5]在研究中發(fā)現(xiàn)，當?shù)皖l信息豐富時FWI目標函數(shù)的局部極值較少，由此提出多尺度FWI策略，即采用從低頻逐漸反演到高頻，或從大網(wǎng)格逐漸反演到小網(wǎng)格的多尺度反演方法，可降低FWI對初始速度模型的依賴性。Pratt等[6]提出頻率域FWI方法，更易實現(xiàn)從低頻到高頻的反演策略，但該方法在低頻成分不足時，無法恢復(fù)模型的長波長結(jié)構(gòu)。Sirgue等[7]利用低通濾波器獲得炮記錄低頻信息，進而構(gòu)建初始速度場。Chen等[8]提出時間衰減的FWI方法，可降低每一步反演過程中的參數(shù)要求，從淺層向深層做反演，較好地緩解“周波跳躍”。黃建平等[9]將頻率—波數(shù)濾波引入FWI，獲得了較好反演效果。

由于地震有效頻帶寬度和去噪處理的影響，地震數(shù)據(jù)往往缺失低頻信息。Shin等[10-11]提出Laplace域FWI方法，通過加入阻尼波場獲得波場的零頻分量，該方法能反演出模型的長波長特征，但反演的深度受阻尼系數(shù)制約。Wu等[12]認為地震數(shù)據(jù)本身雖低頻成分不足，但其Hilbert包絡(luò)面中含有豐富低頻成分。據(jù)此采用包絡(luò)殘差的思想構(gòu)建目標函數(shù)，并推導(dǎo)出對應(yīng)的梯度公式，其反演結(jié)果能恢復(fù)地下速度的背景信息，再利用該速度作為常規(guī)FWI的初始模型，最終得到較理想反演結(jié)果。胡勇等[13]提出解調(diào)包絡(luò)方法重構(gòu)地震記錄中缺失的低頻信息，有效降低了FWI的非線性程度。包乾宗等[14]對炮記錄做低通濾波，并對濾波結(jié)果利用對數(shù)包絡(luò)提取地震數(shù)據(jù)中的低頻成分，反演得到較好的初始模型。Zhang等[15]利用褶積和反褶積方法在重構(gòu)炮記錄低頻數(shù)據(jù)的同時消除了震源函數(shù)的影響，一定程度上緩解了低頻數(shù)據(jù)缺失和震源函數(shù)不準確導(dǎo)致的“周波跳躍”。Guo等[16]將參考點震源信息引入平面波多尺度反演，控制射線參數(shù)，逐步恢復(fù)低頻數(shù)據(jù)，反演效果顯著提升。

通量校正運輸(Flux Corrected Transport，F(xiàn)CT)技術(shù)最早由Boris等[17]和Book等[18]提出，用于消除不連續(xù)點附近的梯度劇烈變化，后也被應(yīng)用于地球物理領(lǐng)域。Fei等[19]將FCT思想引入逆時偏移，獲得無數(shù)值頻散的高分辨率逆時偏移結(jié)果。吳國忱等[20]用FCT技術(shù)做聲波方程正演模擬，一定程度上壓制了數(shù)值頻散。陳可洋[21]提出一種優(yōu)化的FCT正演技術(shù)，將對角元素引入波場校正，通過聲波方程正演模擬，證明該方法能更好地壓制數(shù)值頻散，提高正演模擬精度。趙威等[22]推導(dǎo)了彈性波方程的FCT公式，將FCT引入彈性波方程正演模擬。Kalita等[23]在FWI中利用FCT對波場進行校正，增加波場的低頻成分，降低了FWI對初始模型的依賴性，在模型測試和海上實際資料應(yīng)用中都取得良好效果。

本文基于FCT理論，推導(dǎo)出基于雙對角通量校正(Double Diagonal Flux Corrected Transport，DF)的二維聲波方程格式，并進行二維聲波方程正演模擬，討論其重構(gòu)低頻信息的有效性； 然后推導(dǎo)基于雙對角FCT格式的FWI梯度更新公式，通過動態(tài)控制通量校正參數(shù)實現(xiàn)多尺度反演。為驗證該方法的有效性，采用改進Marmousi模型分別進行不含噪和含噪炮記錄的FWI測試。結(jié)果表明，該方法能在反演初期構(gòu)建宏觀速度場信息，最終獲得高精度反演結(jié)果。對于含噪炮記錄，也能得到較好反演結(jié)果，證實該方法具有一定抗噪性。與其他方法[24]不同，該方法不需修改目標函數(shù)，省去了繁冗運算，只需通過簡單的波場校正和多尺度反演策略即可獲得較好的反演效果。

1 基于雙對角通量校正的多尺度全波形反演理論

時間域FWI的強非線性，要求在反演初期取得較準確的初始速度信息[25]。該初始速度信息可視為FWI的先驗信息。初始速度不準確，易導(dǎo)致“周波跳躍”現(xiàn)象，進而使反演陷入局部極值。為此，提出基于雙對角通量校正的多尺度FWI方法，忽略波場的劇烈震蕩，拓展波場的低頻信息，使目標函數(shù)收斂，同時壓制噪聲和反演假象，恢復(fù)模型的大尺度成分。

1.1 FCT基礎(chǔ)理論

Boris等[17]和Book等[18]提出的FCT理論，在流體力學(xué)研究中用于減少數(shù)值模擬時產(chǎn)生的高頻激蕩現(xiàn)象； 在地球物理領(lǐng)域主要用于解決波動方程正演模擬中產(chǎn)生的頻散。其主要步驟有三： ①計算x、z方向的擴散通量； ②對所求擴散通量得到的波場進行平滑校正； ③進行反擴散通量校正。

二維聲波方程可表示為

(1)

式中：u為時間域波場；v為地震波速度；f為震源項；t為時間。

Su=f

(2)

對式(1)采用時間2階、空間2N階差分離散，可得

(3)

式中：ai為差分系數(shù)；k、m、n分別為t、x、z方向的序號； Δx、Δz為網(wǎng)格間距。

第一步，計算k時刻x、z方向波場的擴散通量

(4)

式中ξ為擴散通量的校正參數(shù)，根據(jù)流體力學(xué)理論可知，其取值范圍是0～1。

第二步，進行波場平滑校正

(5)

第三步，進行反擴散通量計算，求出k+1時刻的x、z方向波場擴散通量

(6)

式中ψ為反擴散通量校正參數(shù)。其取值范圍也是0～1。之后進行x、z方向的反擴散通量計算

(7)

得到最終的校正波場

(8)

式中

(9)

如式(9)所示，在求取反擴散通量的過程中，需要求取最大值、最小值、絕對值、sign函數(shù)，這些函數(shù)均為不可微函數(shù)，而地震信號卻為連續(xù)信號。因此，為了在FWI中獲得有效的反演梯度，當將FCT引入FWI時，僅對波場進行擴散通量校正，不需考慮上述第三步。

1.2 雙對角FCT理論

在地震波場數(shù)值模擬中，相鄰點波場間有很強的相互關(guān)系。因此，考慮對角線方向的波場對校正目標點波場的影響，獲得最佳波場校正效果，在FCT過程中加入雙對角方向波場值，計算各方向擴散通量，雙對角通量校正的擴散通量可表示為

(10)

利用式(10)求出不同方向擴散通量，對地震波場進行校正，雙對角FCT對應(yīng)的校正波場表示為

(11)

根據(jù)文獻[23]，為確定雙對角通量校正參數(shù)ξ的取值范圍，采用改進Marmousi模型進行二維聲波正演模擬。如圖1所示，為便于區(qū)分直達波與反射波，在Marmousi模型上層添加速度為2400m/s的地層。改進Marmousi模型的網(wǎng)格橫縱數(shù)目為596×180，網(wǎng)格間距Δx=Δz=8m，整個模型橫向距離為4768m，縱向深度為1440m。震源加載20Hz雷克子波，震源坐標為(2320m，8m)，沿地表均勻布設(shè)596個檢波器，檢波器間距為8m，采樣間隔為0.5ms，采樣時間為2.25s，采用30層PML邊界條件對邊界反射進行處理。

圖1 改進Marmousi模型

圖2分別表示常規(guī)正演與不同校正參數(shù)所對應(yīng)的單炮記錄。從圖中可看出，隨著ξ的逐漸增大單炮記錄逐漸平滑。

對圖2中各單炮記錄抽取2000m處數(shù)據(jù)，進行頻譜分析，頻譜分析結(jié)果如圖3所示。從圖3可見，采用雙對角通量校正技術(shù)獲得的地震記錄信號頻率隨ξ的增大而降低，而當ξ>0.08時，出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象，無法獲得有效的炮記錄，因此，雙對角通量校正取ξ=0.08能獲得最佳信號拓展低頻效果。對FWI而言，低頻信息的增多能緩解“周波跳躍”現(xiàn)象，降低反演問題的非線性程度，從而有效地降低FWI對初始模型的依賴性[26]。

圖2 常規(guī)正演與雙對角FCT的單炮記錄

圖3 常規(guī)正演與雙對角FCT(ξ=0.02、0.04、0.06、0.08)對應(yīng)的頻譜分析

1.3 基于雙對角通量校正的多尺度FWI

FWI通過優(yōu)化算法不斷迭代降低觀測數(shù)據(jù)與預(yù)測數(shù)據(jù)之間的殘差，使反演結(jié)果不斷逼近真實速度信息，其中求取反演梯度和步長是FWI的關(guān)鍵。

本文采用L2范數(shù)作為目標函數(shù)

(12)

式中：dobs為觀測數(shù)據(jù)；dcal為預(yù)測數(shù)據(jù)；ns為震源數(shù)量；nr為檢波點數(shù)量。利用Plessix[27]給出的伴隨狀態(tài)法可知L2范數(shù)目標函數(shù)的反演梯度為

(13)

將式(2)代入式(13)，可得

(14)

上兩式中：S是波場延拓算子；S-1是波場反傳算子；R是將波場限制到檢波點位置的算子； δd=dcal-dobs為預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)之間的殘差，即反演的伴隨源。因此，反演梯度為正傳波場關(guān)于時間上的二階導(dǎo)數(shù)和伴隨源的反傳波場在時間上的內(nèi)積。

(15)

本文利用三點拋物插值法求取反演步長，該方法需已知三個步長值α1、α2、α3和三個對應(yīng)的目標函數(shù)值E(α1)、E(α2)、E(α3)，一般設(shè)α1=0，且滿足

α1<α2<α3,E(α1)>E(α2)，E(α2)

(16)

假設(shè)目標函數(shù)E(α)與步長滿足以下二次關(guān)系

(17)

求取E(α)的極值：E′(α)=2aα+b=0

可知最佳步長為：αopt=-b/2a。根據(jù)之前給定的三個試探步長和對應(yīng)的目標函數(shù)值可知最佳步長表達式為

(18)

在求得反演梯度和最佳步長后，利用下式更新速度模型

vn=vn-1+αopt?vE(v)DF

(19)

式中：vn為更新后的速度；vn-1為更新前的速度。

基于雙對角通量校正的多尺度FWI過程中針對正傳波場和殘差的反傳波場進行校正，求取梯度，修改過程較為簡單，因此能夠推廣至其他目標函數(shù)的反演方法。

1.4 基于雙對角通量校正的多尺度FWI流程

本文提出基于雙對角通量校正的多尺度FWI方法對現(xiàn)有FWI程序有很強的適應(yīng)性。此方法不需改變現(xiàn)有FWI框架。在反演算法中僅加入校正后的波場，并且加入校正波場所產(chǎn)生的計算成本對于整個FWI程序而言可忽略不計。較大的ξ能增加波場中的低頻成分，減少“周波跳躍”現(xiàn)象，為反演提供一個較平滑的梯度。通過控制通量校正參數(shù)ξ的值逐漸減小直至0，實現(xiàn)多尺度反演，上一個ξ值得到的反演結(jié)果作為下一個ξ值反演的初始模型，當ξ=0時，進行常規(guī)的FWI。在初始模型速度信息不準確的情況下，獲得較好的反演結(jié)果。

圖4 基于雙對角通量校正的多尺度FWI流程

2 模型測試

2.1 不精確初始速度模型條件下的FWI

為了測試基于雙對角通量校正的多尺度FWI方法的有效性，利用修改Marmousi模型(圖1)進行模型測試。在深度8m處均勻布設(shè)60個炮點，炮間距為80m。沿地表均勻布設(shè)596個檢波器，檢波器間距為8m。震源加載20Hz雷克子波，采樣間隔為0.5ms，采樣時長為2.25s。

圖5為反演的初始速度模型。從地表到深度240m范圍內(nèi)速度恒為2400m/s，深度240m以下速度開始線性增加，范圍是2400～5800m/s，初始模型中不含真實模型的構(gòu)造信息。

圖5 初始速度模型

分別用常規(guī)FWI(圖6)和基于雙對角通量校正的多尺度FWI(圖7)進行90次迭代，得到反演結(jié)果。基于雙對角通量校正的多尺度FWI是通過逐漸降低ξ值實現(xiàn)多尺度反演。用該方法分別在ξ={0.08、0.06、0.04、0.03、0.02、0.01、0}時迭代{6、8、8、8、10、15、35}次，得到最終結(jié)果，其中ξ=0.08、0.04、0.02、0的反演結(jié)果如圖7所示。從圖7(基于本文方法反演結(jié)果)可看出，隨著ξ值的逐漸降低，反演結(jié)果逐步向真實模型收斂，反演較好地恢復(fù)了模型的長波長速度信息，為ξ=0的常規(guī)FWI提供了較準確的初始速度模型。最終ξ=0的反演結(jié)果，恢復(fù)了真實模型的主要構(gòu)造，包括三個斷層和下部的背斜等細小的構(gòu)造。而從圖6的常規(guī)反演結(jié)果可知，常規(guī)FWI能反演得到淺層信息，但因初始速度模型不準確，導(dǎo)致深層和細小構(gòu)造信息的反演不精確。

圖6 常規(guī)FWI結(jié)果

圖7 不同通量校正參數(shù)的反演結(jié)果

圖8和圖9分別為ξ=0.08和利用常規(guī)FWI進行第一次迭代所獲得的梯度場，對比可看出本文方法能在反演初期提供一個較平滑的反演梯度，從而恢復(fù)了模型長波長信息。

圖8 ξ=0.08的第一次迭代梯度

圖9 常規(guī)FWI的第一次迭代梯度剖面

為了定量分析本文方法較常規(guī)FWI方法的準確程度，抽取初始模型、真實模型和兩種反演方法在1200m和3200m兩處最終結(jié)果的速度曲線(圖10、圖11)并進行對比。

從圖10和圖11可看出，基于本文方法所得到的曲線較常規(guī)FWI的更接近于真實速度模型，對速度突變更敏感。將從圖10和圖11兩種方法得到的速度曲線分別與真實速度模型求取相關(guān)系數(shù)。圖10中常規(guī)FWI方法的相關(guān)系數(shù)為0.8216，基于本文方法的相關(guān)系數(shù)為0.8921。圖11中常規(guī)FWI方法的相關(guān)系數(shù)為0.8723，本文方法的相關(guān)系數(shù)為0.9365。因此，采用基于雙對角通量校正方法進行反演能有效降低FWI問題的非線性程度，降低對初始模型的依賴性，提高反演結(jié)果的準確性。

圖10 1200m處的速度曲線對比

圖11 3200m處的速度曲線對比

2.2 含噪觀測記錄的測試

實際應(yīng)用中，在采集原始地震數(shù)據(jù)時，往往會存在各種不確定干擾因素，使得采集的地震數(shù)據(jù)含有噪聲[28-30]。為了驗證基于本文方法的抗噪性，在利用真實模型得到觀測炮記錄后，再使用SU軟件在觀測記錄中加入能量為原始信號最大振幅4%的高斯噪聲，獲得加入噪聲后的觀測記錄(圖12)。

圖12 加入高斯噪聲后的觀測記錄

反演參數(shù)與無噪聲反演相同。常規(guī)FWI結(jié)果如圖13所示，淺層有效信號較強區(qū)域有一定的反演效果，而在深層有效信號被噪聲淹沒，速度信息無法恢復(fù)。基于本文方法結(jié)果如圖14所示。對比圖14a～圖14d可看出，在反演初期ξ較大時，噪聲對反演結(jié)果造成一定的干擾，淺層出現(xiàn)斑狀速度異常，但恢復(fù)了模型的長波長信息。為下一階段反演打下良好基礎(chǔ)。隨著ξ的不斷減小，反演結(jié)果逐漸收斂于真實模型。為了對比反演結(jié)果的準確性，抽取兩種方法1600m和3600m處的反演結(jié)果(圖15、圖16)進行對比。從圖中能看出基于本文方法結(jié)果更逼近于真實速度曲線，分別對圖15和圖16兩種方法結(jié)果與真實速度求取相關(guān)系數(shù)。圖15中本文方法的相關(guān)系數(shù)為0.9462，常規(guī)FWI的相關(guān)系數(shù)為0.9204。圖16中本文方法的相關(guān)系數(shù)為0.9199，常規(guī)FWI的相關(guān)系數(shù)為0.8197。從兩種反演結(jié)果的相關(guān)系數(shù)看，當存在一定噪聲干擾時，本文方法依舊能獲得較好結(jié)果，表明它具有一定抗噪優(yōu)勢。

圖14 存在高斯噪聲時不同通量校正參數(shù)的反演結(jié)果

圖15 1600m處的速度曲線對比

圖16 3600m處的速度曲線對比

3 結(jié)論

基于雙對角通量校正的多尺度FWI是對常規(guī)時間域FWI方法的一種改進。通過從大到小動態(tài)調(diào)整通量校正參數(shù)ξ的取值，實現(xiàn)多尺度的反演策略，獲得準確的反演結(jié)果。應(yīng)用該方法進行模型測試、抗噪性分析，得到如下認識和結(jié)論：

(1)采用雙對角通量校正方法對波場進行校正，拓展波場中的低頻信息，低頻信息可靠，能夠緩解初始速度模型不精確所導(dǎo)致的“周波跳躍”現(xiàn)象。

(2)基于雙對角通量校正的多尺度FWI較常規(guī)FWI而言，在反演初期，能從高頻主導(dǎo)的地震記錄中提供較平滑的梯度場，恢復(fù)長波長速度信息，在初始模型較差時，也能得到高精度FWI結(jié)果。

(3)當?shù)卣鹩涗浐幸欢ㄔ肼晻r，該方法依舊能獲得較好反演效果，說明該方法具有較好抗噪性。

(4)基于本文方法，運用到FCT技術(shù)中的擴散和校正技術(shù)，對波場進行校正，并不需對目標函數(shù)做修改，對現(xiàn)有的FWI程序有較強可移植性。