應(yīng)用自適應(yīng)矩估計的快速最小二乘逆時偏移

吳 丹 吳海莉 李 群 張向陽 劉樹仁

(①中國石油勘探開發(fā)研究院西北分院，甘肅蘭州 730030； ②中國石油天然氣集團有限公司物聯(lián)網(wǎng)重點實驗室，甘肅蘭州 730030)

0 引言

逆時偏移(RTM)互相關(guān)成像條件可闡述為：對正傳得到的震源波場和反傳得到的接收波場進行互相關(guān)，得到地下反射界面的構(gòu)造圖像[1]。實際上，該成像條件是Born正演算子的伴隨算子[2]。當炮檢距有限、數(shù)據(jù)不完美(帶限、含噪、缺失)時，常規(guī)RTM很難得到滿意的成像結(jié)果，無法滿足復(fù)雜構(gòu)造和巖性油氣藏的地震成像需求。

解決上述問題的一種嚴格做法是將地震成像看作最小二乘意義下的最優(yōu)化問題，通過最優(yōu)化算法最小化Born正演數(shù)據(jù)和觀測反射地震數(shù)據(jù)之間的差異，得到最佳匹配觀測數(shù)據(jù)的地下模型估計，即最小二乘偏移[3-6]。在最小二乘偏移中，可以選擇不同的Born正演/偏移算子，其中最精確的是雙程波算子，即最小二乘逆時偏移(LSRTM)[7-9]。LSRTM每次迭代大概需要兩次所有共炮點道集的RTM，通常需要迭代近十次，因此計算成本極高。目前主要有兩類方法可以降低計算成本：一類方法是利用炮編碼技術(shù)對疊前炮集數(shù)據(jù)進行壓縮[10-22]，降低Born正演和RTM的計算成本，如隨機相位編碼、平面波相位編碼、振幅編碼、確定性炮編碼等； 另一類方法是利用預(yù)條件算子加快迭代的收斂速度，減少迭代次數(shù)，例如近似Hessian對角矩陣[23-24]、去模糊化濾波器[25-26]、Hessian矩陣近似逆[27]、構(gòu)造導(dǎo)向濾波[28-29]等。

本文從人工智能領(lǐng)域引入自適應(yīng)矩估計(A-dam)方法[30]解決LSRTM的效率問題。Adam方法結(jié)合了深度學(xué)習(xí)中兩種常用方法的優(yōu)勢：一種是自適應(yīng)梯度(AdaGrad)方法[31]，該方法能夠自適應(yīng)地對梯度進行歸一化，改善稀疏特征的更新效果； 另一種是均方根傳播(RMSProp)方法[32]，該方法利用之前迭代的計算結(jié)果，對當前迭代的更新量和歸一化向量進行平滑處理，提高迭代收斂的穩(wěn)定性，在梯度非穩(wěn)態(tài)時效果顯著。相應(yīng)地，在LSRTM中，Adam方法有潛力壓制隨機選取部分炮集數(shù)據(jù)進行偏移而引入的成像噪聲，并補償觀測系統(tǒng)有限和復(fù)雜速度模型帶來的照明不均勻問題。另外，Adam方法也有潛力壓制炮編碼LSRTM中存在的串擾噪聲。最后通過SEG/EAGE鹽丘模型驗證了Adam方法對LSRTM的有效性。

1 方法和理論

1.1 逆時偏移

在常密度聲波介質(zhì)中，地震波的傳播過程可用頻率域波動方程描述

[?2+ω2m(x)]u(x,ω)=f(ω)δ(x-xs)

(1)

式中：x表示空間位置坐標；xs為炮點坐標；m(x)為模型慢度場的平方；u(x,ω)是頻率為ω的地震波場；f(ω)為震源函數(shù)譜。

將模型m(x)分解成兩部分

m(x)=b(x)+r(x)

(2)

式中：b(x)為背景模型，是對m(x)平滑后的結(jié)果，只影響地震波場的走時，而不影響地震波場的振幅，且不會產(chǎn)生反射波；r(x)為擾動模型，不影響地震波場的走時，但入射波場遇到反射界面會產(chǎn)生反射波，且反射波的振幅與反射系數(shù)的大小有關(guān)，因此擾動模型反映的是地下模型的層位構(gòu)造信息。為描述簡化起見，本文將擾動模型稱為反射系數(shù)模型。與模型分解相對應(yīng)，地震波場u(x,ω)也可以分解成兩部分：透射波場u0(x,ω)和反射波場Δu(x,ω)。經(jīng)線性近似，式(1)可分解為

[?2+ω2b(x)]u0(x,ω)=f(ω)δ(x-xs)

(3)

[?2+ω2b(x)]Δu(x,ω)=-ω2r(x)u0(x,ω)

(4)

式(3)的解可以用Green函數(shù)表示為

u0(x,ω)=f(ω)G(x,xs,ω)

(5)

式中G(x,xs,ω)為炮點位于xs的Green函數(shù)。反射地震數(shù)據(jù)d(xr,xs,ω)可以通過求解式(4)的Δu(x,ω)并在接收點xr采樣得到。數(shù)學(xué)上，該過程可表示為

d(xr,xs,ω)=-ω2f(ω)×

(6)

該式即為線性Born正演公式。對應(yīng)地，線性Born正演的共軛算子可表示為

G*(xs,x,ω)G*(x,xr,ω)d(xr,xs,ω)

(7)

式中上標“*”表示共軛。式(7)即為互相關(guān)成像條件，也證明了互相關(guān)成像條件為線性Born正演的共軛。

1.2 炮編碼逆時偏移

從式(7)可以看出，RTM的計算成本和共炮點道集個數(shù)成正比。為了減少計算量，炮編碼技術(shù)可以對疊前共炮點道集進行壓縮。目前常見的炮編碼技術(shù)有隨機相位編碼、平面波相位編碼、振幅編碼等。Krebs等[33]的數(shù)值結(jié)果表明單點隨機時間相位編碼可以取得較好的收斂效果，因此本文將采用這種編碼方式對數(shù)據(jù)進行壓縮。該方法的編碼函數(shù)定義為

(8)

式中：ps為炮編碼實現(xiàn)次數(shù)的索引；N為炮編碼實現(xiàn)的總次數(shù)；γ是隨機符號序列，取值為+1或-1。經(jīng)過炮編碼后的地震數(shù)據(jù)為

(9)

式中dobs為觀測數(shù)據(jù)。對應(yīng)地，相同的編碼函數(shù)也用于編碼震源。由于波傳播算子和震源函數(shù)為線性關(guān)系，因此炮編碼后的震源波場可表示為

δ(x-xs)G(x,xs,ω)

(10)

同樣地，線性Born正演和輸入震源也是線性關(guān)系，因此只需要修改震源函數(shù)即可完成炮編碼下的線性Born正演，即

G(xr,x,ω)r(x)S(x,ps,ω)

(11)

對應(yīng)的炮編碼RTM公式為

(12)

從上式可以看出，炮編碼RTM的計算成本與炮編碼的次數(shù)成正比，而炮編碼次數(shù)往往遠小于炮數(shù)，因此經(jīng)過炮編碼后線性Born正演和RTM的計算成本都得到大幅降低。

1.3 基于Adam方法的最小二乘逆時偏移

由于一次反射波地震數(shù)據(jù)與反射系數(shù)呈線性關(guān)系，因此線性Born正演公式(式(6))可以寫成矩陣形式

d=Lm

(13)

式中：d為Born正演數(shù)據(jù)；m為模型參數(shù)；L為線性Born正演算子。LSRTM目標函數(shù)可表示為

(14)

式中：dobs為離散化后的觀測反射地震數(shù)據(jù)； ‖?‖2表示L2范數(shù)。

式(14)的反問題可以通過梯度下降(GD)法求解。在深度學(xué)習(xí)領(lǐng)域，一種更加有效的最優(yōu)化算法是小批量梯度下降法。與GD法相比，該算法隨機地打亂樣本集，然后只使用樣本集的一部分樣本(即“小批量”的來歷)計算梯度，從而快速地更新網(wǎng)絡(luò)參數(shù)。在深度神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)過程中，往往需要大量的訓(xùn)練樣本以避免過擬合，小批量梯度下降法只需要使用小部分樣本就可以對模型參數(shù)進行更新，相比全樣本集GD法，效率往往高很多。

但是小批量GD往往容易劇烈震蕩。為解決這一問題，動量法應(yīng)運而生。該方法考慮了之前迭代的梯度信息，對當前迭代次數(shù)中的梯度進行平滑。更準確地說，動量法采用指數(shù)加權(quán)函數(shù)對前面所有的梯度進行平均，得到新的梯度，即

(15)

式中vt-1、vt分別表示第t-1和第t次迭代的梯度；β1是控制梯度平滑程度的超級參數(shù)，β1∈[0,1)。β1越大，更新方向越平滑。如果β1=0，上式就退化為常規(guī)的梯度下降法； 如果β1太大，就會平滑掉一些有用的更新信息。大多數(shù)情況下，β1=0.9是一個比較合理的選擇。

式(15)通常會產(chǎn)生邊界效應(yīng)，需對平滑后的梯度進行校正。校正后的梯度為

(16)

式(16)的校正后更新方向梯度沒有考慮非穩(wěn)態(tài)性。具體到LSRTM方法，它沒有考慮計算的梯度方向存在照明不均勻的問題。RMSProp方法應(yīng)用梯度平方的指數(shù)加權(quán)平均近似照明矩陣

(17)

式中：st-1、st分別表示第t-1和第t次迭代的照明矩陣；β2為控制照明矩陣平滑程度的超級參數(shù)，與β1類似，一般取0.9。同樣地，也需要校正邊界效應(yīng)

(18)

然后，利用RMSProp估計的照明矩陣歸一化校正后的梯度(式(16))，得到模型的更新量

(19)

式中ε為穩(wěn)定因子，以保證除法的穩(wěn)定性。最終，反射系數(shù)模型更新可表示為

mt=mt-1+λΔm

(20)

式中λ在深度學(xué)習(xí)中稱為學(xué)習(xí)率，在數(shù)值優(yōu)化中稱為步長。步長可以隨著迭代次數(shù)的增加逐步減小

(21)

式中：λt-1、λt分別為第t-1和第t次迭代的步長；η為衰減率，本文將其設(shè)置為1。

圖1為Adam法LSRTM計算流程。該算法共

圖1 Adam法LSRTM計算流程

有3層循環(huán)：①最內(nèi)層循環(huán)是使用每個小批量共炮點道集數(shù)據(jù)計算梯度； ②中層循環(huán)則是利用Adam方法更新反射系數(shù)模型； ③最外層循環(huán)執(zhí)行的是打亂并劃分的共炮點道集數(shù)據(jù)，以及自適應(yīng)地減小步長。

1.4 炮編碼Adam法LSRTM

在動態(tài)炮編碼LSRTM的每次迭代過程中，雖然更新了炮編碼，但是經(jīng)過偏移后提供的有效梯度信息具有相干性，而動態(tài)編碼引入的串擾噪聲不相干。由于Adam法采用指數(shù)加權(quán)平均對前面所有迭代步中計算得到的梯度和近似照明矩陣進行平滑，因此Adam法具有壓制串擾噪聲的能力。為此，本文修改Adam法以適用于炮編碼LSRTM。修改后的計算流程和梯度下降法動態(tài)炮編碼LSRTM類似，唯一的不同之處在于Adam方法采用自適應(yīng)矩估計法計算梯度。完整的算法流程如圖2所示。

圖2 炮編碼Adam法LSRTM計算流程

2 數(shù)值試驗

采用SEG/EAGE 鹽丘模型測試本文提出的算法。該模型的橫向和縱向空間采樣間隔均為30m。圖3a為該模型的速度場，由式(2)可計算其反射系數(shù)模型(圖3b)。采用主頻為15Hz、時間采樣間隔為1ms的Ricker子波進行有限差分線性Born正演?？偣舱萘?63炮，起始炮點位于(0，0)，間隔為120m； 采用645個接收點的固定排列接收，起始點位于(0，0)，接收點間隔為30m。

圖4a為RTM經(jīng)拉普拉斯濾波后的結(jié)果，可以看出成像結(jié)果分辨率較低，且鹽下成像振幅較弱。圖4b為迭代10次后共軛梯度(CG)法的LSRTM結(jié)果，與圖4a相比，提高了分辨率，而且一定程度上恢復(fù)了鹽下的反射振幅，但仍然與真實反射系數(shù)(圖3b)相差較遠。通過迭代更多次數(shù)，可以進一步改善LSRTM的成像結(jié)果，但10次迭代的計算成本已相當于20次全部共炮點道集的RTM。圖4c和圖4d均為迭代10次后的Adam法LSRTM成像結(jié)果，其中前者的批量參數(shù)為16，后者為4，后者的成像結(jié)果遠優(yōu)于前者，說明批量參數(shù)對成像效果影響很大。對不同批量參數(shù)進行測試，最終選擇成像效果最好的批量參數(shù)4。對比圖4d與圖4b可以看出，在相同的迭代次數(shù)下，Adam法的效果遠優(yōu)于CG法，從而驗證了Adam法比共軛梯度法的收斂速度快得多。

圖5是CG法、批量參數(shù)為16的Adam法、批量參數(shù)為4的Adam法的LSRTM數(shù)據(jù)擬合誤差曲線和模型誤差曲線，可見批量參數(shù)為4時Adam算法具有最快的收斂速度。進一步對比三種方法LSRTM結(jié)果在橫向10km處的反射系數(shù)波形(圖6)，批量參數(shù)為4時的Adam法偏移結(jié)果最接近于真實反射系數(shù)模型。

圖3 鹽丘速度(a)和反射系數(shù)(b)模型

圖4 RTM和LSRTM結(jié)果對比(a)RTM+拉普拉斯濾波； (b)CG法LSRTM； (c)Adam法LSRTM(批量參數(shù)為16)； (d)Adam法LSRTM(批量參數(shù)為4)

圖5 三種LSRTM法的誤差曲線對比(a)數(shù)據(jù)擬合誤差； (b)反射系數(shù)平方誤差

圖6 三種方法LSRTM結(jié)果與真實反射系數(shù)的單道波形對比

結(jié)合Adam方法和炮編碼可以進一步提高LSRTM的計算效率。圖7a為應(yīng)用炮編碼GD法進行LSRTM迭代163次后的成像結(jié)果，其計算成本相當于2次全部共炮點道集的RTM，但已經(jīng)優(yōu)于迭代10次的LSRTM結(jié)果(圖4b)。圖7b為應(yīng)用炮編碼Adam法進行LSRTM迭代163次的成像結(jié)果。兩種方法的計算成本大致相同，但炮編碼Adam法LSRTM的結(jié)果要優(yōu)于GD法，計算結(jié)果十分接近真實的反射系數(shù)模型。

圖7 炮編碼LSRTM反射系數(shù)成像結(jié)果(a)GD法； (b)Adam法

圖8為炮編碼LSRTM的數(shù)據(jù)擬合誤差曲線和模型誤差曲線，可以看出，相比于GD法，炮編碼Adam法LSRTM具有更快的收斂速度。進一步對比橫向10km處的反射系數(shù)波形，炮編碼Adam法LSRTM的成像結(jié)果更接近于真實反射系數(shù)模型(圖9)。

對原始Born正演數(shù)據(jù)加入隨機噪聲，其信噪比為-7.7dB。圖10為加噪數(shù)據(jù)的炮編碼GD法和Adam法LSRTM的結(jié)果，可見Adam算法的鹽下成像結(jié)果更清晰，整體分辨率更高。圖11為兩種方法的數(shù)據(jù)擬合誤差曲線和模型誤差曲線，可以看出Adam算法的收斂速度更快。兩種方法成像結(jié)果在橫向10km處的反射系數(shù)波形與真實反射系數(shù)的對比如圖12所示，可見在數(shù)據(jù)存在噪聲的情況下，炮編碼Adam法LSRTM成像結(jié)果更接近于真實反射系數(shù)模型。需要注意的是，反演類成像方法對數(shù)據(jù)的信噪比一般要求較高，因此建議在反演成像之前對數(shù)據(jù)做去噪處理，以盡可能減小噪聲對反演結(jié)果的影響。

圖8 炮編碼兩種方法LSRTM法的誤差曲線對比(a)數(shù)據(jù)擬合誤差； (b)反射系數(shù)平方誤差

圖9 炮編碼兩種方法LSRTM結(jié)果與真實反射系數(shù)的單道波形對比

圖10 含噪數(shù)據(jù)炮編碼兩種方法LSRTM結(jié)果(a)GD法； (b)Adam法

應(yīng)用聲波方程有限差分法正演數(shù)據(jù)作為觀測數(shù)據(jù)分析炮編碼GD法和Adam法LSRTM對正演算子誤差的敏感性。圖13為兩種方法的成像結(jié)果，可見Adam算法的對鹽下成像結(jié)果更清晰，整體的分辨率更高。圖14為這兩種方法的數(shù)據(jù)擬合誤差曲線和模型誤差曲線，可以看出Adam算法的收斂速度更快。對比橫向10km處的兩種方法的反射系數(shù)波形(圖15)，可見炮編碼Adam法LSRTM反演結(jié)果更接近真實反射系數(shù)模型。

圖12 含噪數(shù)據(jù)炮編碼兩種方法LSRTM結(jié)果與真實反射系數(shù)的單道波形對比

圖13 有限差分法正演數(shù)據(jù)炮編碼兩種方法LSRTM結(jié)果(a)GD法； (b)Adam法

圖14 有限差分法正演數(shù)據(jù)炮編碼兩種方法LSRTM的誤差曲線對比(a)數(shù)據(jù)擬合； (b)反射系數(shù)

使用偏移速度為準確速度的95%測試炮編碼GD法和Adam法LSRTM的對速度誤差的敏感性。圖16為兩種方法對Born正演數(shù)據(jù)的成像結(jié)果。可以看出，不管是GD法還是Adam法，都無法對鹽丘模型準確成像，表明LSRTM對速度精度要求很高。

圖15 有限差分法正演數(shù)據(jù)兩種方法炮編碼LSRTM結(jié)果與真實反射系數(shù)的單道波形對比

圖16 偏移速度有誤差時炮編碼兩種方法LSRTM結(jié)果(a)GD法； (b)Adam法

3 結(jié)論

本文首次采用深度學(xué)習(xí)領(lǐng)域中的優(yōu)化算法——Adam算法解決LSRTM的計算成本問題。SEG/EAGE 鹽丘模型數(shù)值實驗表明，通過結(jié)合Adam算法和炮編碼技術(shù)，只需要消耗大約兩次完整炮集數(shù)據(jù)RTM的計算成本，就能夠得到比RTM分辨率更高、更精確的反射系數(shù)結(jié)果，從而使LSRTM應(yīng)用于實際地震數(shù)據(jù)的成像成為可能，為復(fù)雜構(gòu)造成像提供更完善的信息。

另外，本文進一步分析了Adam法LSRTM對隨機噪聲、正演算子誤差、速度誤差的敏感性，結(jié)果表明：在存在隨機噪聲和正演誤差的情況下，Adam方法的收斂速度都要優(yōu)于GD法； 但是在速度存在誤差的情況下，Adam方法與其他方法一樣，都無法得到準確的成像結(jié)果。

值得注意的是，直接從實際數(shù)據(jù)中的波形信息估計反射系數(shù)模型十分困難。深度學(xué)習(xí)也許能夠從地震數(shù)據(jù)中提取更可靠的特征，提供對地下巖性參數(shù)更加有效的估計，這也是下一步的研究方向。