黃韻博, 黃建平, 李振春, 劉博文

(中國石油大學(華東)地球科學與技術學院,山東青島 266580)

地震偏移成像的目標是通過地震數(shù)據(jù)對地下構造對應的反射系數(shù)模型進行精確成像[1]。提高成像的分辨率對于非常規(guī)油氣勘探及小尺度地質(zhì)目標的刻畫十分重要。但是,常規(guī)的偏移算子只是正演算子的轉(zhuǎn)置,而不是它的逆,無法完全逆轉(zhuǎn)地震波傳播效應[2],因此常規(guī)偏移方法存在分辨率低、照明不均、偏移假象以及低頻噪音嚴重等問題[3-4]。為了得到更加精確的成像結果,可以將成像問題放入最小二乘反演框架,利用迭代最小二乘法求解,即最小二乘偏移[5-8]。近年來,最小二乘偏移已被證明能夠提高不完整數(shù)據(jù)[9]成像質(zhì)量和均衡不均勻地震波照明。最小二乘偏移包括模型域最小二乘偏移[10-12]和數(shù)據(jù)域最小二乘偏移[13-20]。模型域的最小二乘偏移是通過少量的迭代計算出Hessian矩陣的逆矩陣,將其直接作用于常規(guī)偏移結果,得到高精度成像結果。但是,Hessian矩陣十分龐大,對存儲量要求極高[21]。數(shù)據(jù)域的最小二乘偏移則是通過迭代最小二乘法,逐步更新反射系數(shù)模型,這種方式不需要存儲Hessian矩陣,但是需要大量迭代計算量。近年來,深度學習被廣泛應用于地球物理領域以解決斷層識別[22]、地震數(shù)據(jù)去噪[23]以及初至拾取[24]等問題,并取得了良好的效果。筆者采用循環(huán)對抗神經(jīng)網(wǎng)絡[25]來近似替代Hessian矩陣的逆矩陣,將訓練后的網(wǎng)絡直接作用于逆時偏移[26-27]的結果,快速提高成像質(zhì)量。選取3個速度模型,通過有限差分正演和逆時偏移成像得到訓練數(shù)據(jù)中的輸入數(shù)據(jù),通過雷克子波與反射系數(shù)褶積得到訓練集中的輸出數(shù)據(jù)。同時,再選取另外兩個模型作為驗證集,測試網(wǎng)絡訓練后的預測效果,驗證集不參與模型訓練。

1 方法原理

1.1 偏移理論

對于給定的反射系數(shù)m和線性正演算子L,觀測數(shù)據(jù)d可以表示為

d=Lm.

(1)

地震成像的目標是尋找準確的m來擬合觀測數(shù)據(jù)d,常規(guī)逆時偏移成像可以表示為

mmig=LTd.

(2)

式中,mmig為逆時偏移成像結果的矩陣表示;LT為逆時偏移算子的矩陣表示。LT只是正演算子的轉(zhuǎn)置,并不是它的逆,故逆時偏移的成像結果并不能精確表示真實反射系數(shù)模型。

對真實反射系數(shù)m的求解可以放在最小二乘框架下進行,本文中定義最小二范數(shù)的目標泛函f(m)為

(3)

由此可以得到反射系數(shù)m的最小二乘解為

m=(LTL)-1LTd.

(4)

這里(LTL)-1即Hessian矩陣的逆矩陣,式中LTd為逆時偏移的結果,因此逆時偏移成像結果和真實反射系數(shù)模型之間的關系可以表示為

m=H-1mmig.

(5)

常規(guī)二維聲波方程表達式為

(6)

式中,σ(x)=1/ν(x)為介質(zhì)的慢度場;ν(x)為對應的介質(zhì)速度場;fs(ω)為震源項,表示脈沖震源點在xs處,圓頻率為ω?；诜囱輪栴}近似線性的假設,對于這個二維聲波方程,其Hessian矩陣的表達式可以寫成如下形式:

GT(xs,y,ω)×G(x,xr,ω)GT(y,xr,ω)}dω.

(7)

式中,G(xs,x,ω)為激發(fā)點的格林函數(shù);G(x,xr,ω)為接收點xr的格林函數(shù);fs(ω)為震源子波在頻率域的表示。格林函數(shù)項反映了成像過程中地震波的照明問題,而震源子波項則對成像分辨率產(chǎn)生影響。

結合式(5)和式(7)分析可得,逆時偏移結果受到Hessian矩陣的影響,表現(xiàn)出成像分辨率低、照明不均及偏移噪音等問題。求解得到Hessian矩陣的逆矩陣,將其作用于逆時偏移的成像結果,便可以得到真實的反射系數(shù)模型,這種方式即模型域最小二乘逆時偏移。但是,直接求解Hessian矩陣的逆需要很大的計算量和存儲量,對于實際應用來說適用性不好。另一種方式是數(shù)據(jù)域最小二乘逆時偏移,即通過迭代算法不斷更新反射系數(shù),使得模擬出來的數(shù)據(jù)和觀測數(shù)據(jù)之間的誤差達到最小或期望值。數(shù)據(jù)域最小二乘逆時偏移同樣需要大量的計算成本。本文思想與模型域最小二乘逆時偏移類似,采用循環(huán)對抗神經(jīng)網(wǎng)絡來近似表征Hessian矩陣的逆矩陣,將訓練后的網(wǎng)絡直接作用于逆時偏移的成像結果,得到高質(zhì)量的成像剖面。

1.2 循環(huán)對抗神經(jīng)網(wǎng)絡

對抗神經(jīng)網(wǎng)絡(GAN)在源域向目標域的轉(zhuǎn)換問題上有良好的適用性,常被用于風格遷移。常規(guī)的神經(jīng)網(wǎng)絡是通過反傳算法尋找一個輸入到輸出的最佳映射關系,網(wǎng)絡的本質(zhì)可以看作一個復雜的高階多項式。對抗神經(jīng)網(wǎng)絡一般有兩部分網(wǎng)絡組成,一部分作為生成器,另一部分作為鑒別器。在網(wǎng)絡的訓練過程中,生成器充當常規(guī)神經(jīng)網(wǎng)絡的角色,將輸入轉(zhuǎn)換為和輸出類似的數(shù)據(jù),而鑒別器會以生成器的輸出數(shù)據(jù)作為輸入,嘗試區(qū)分哪些是生成的數(shù)據(jù),哪些是訓練集里的真實數(shù)據(jù)。通過這種對抗式的訓練,生成器和鑒別器都在不斷優(yōu)化,最終得到期望效果的生成器作為預測網(wǎng)絡。

循環(huán)對抗神經(jīng)網(wǎng)絡由2個鏡像對稱的對抗神經(jīng)網(wǎng)絡組成,共同構成一個環(huán)形網(wǎng)絡。如圖1所示,cycleGAN有2個生成器和2個鑒別器。對于源域數(shù)據(jù)X,由紅色虛線代表的單向GAN生成Y1和X1,Y1進入鑒別器鑒別真?zhèn)巍７催^來,作為目標域的數(shù)據(jù)Y也可以作為輸入,由綠色虛線代表的另一個單向GAN生成X2和Y2,X2進入鑒別器做真?zhèn)舞b別。生成器由編碼器、轉(zhuǎn)換器和解碼器構成,編碼-解碼結構與大多數(shù)卷積神經(jīng)網(wǎng)絡類似,用卷積層和反卷積層來實現(xiàn)特征提取和還原,轉(zhuǎn)換器的結構由多層Resnet來實現(xiàn)。鑒別器由基本的卷積神經(jīng)網(wǎng)絡來實現(xiàn),其判別圖像需要從輸入圖像上提取特征,最后添加一維輸出的卷積層來確定所提取的特征是否屬于特定類別。

圖1 循環(huán)對抗神經(jīng)網(wǎng)絡結構Fig.1 Structure of cycle-consistent adversarial neural network

1.3 損失函數(shù)

損失函數(shù)是衡量網(wǎng)絡輸出結果好壞的標尺,同時也是反傳過程中修正張量梯度的依據(jù)。對于一個給定大小的訓練批次,假設源域的數(shù)據(jù)分布為px,目標域的數(shù)據(jù)分布為py,則每一個訓練數(shù)據(jù)x和y的數(shù)據(jù)分布可以表示為x～px和y～py。因此常規(guī)GAN的目標函數(shù)可以表示為

LGAN(GXY,DY,X,Y)=Ey～py[logDY(y)]+Ex～px[log(1-DY(G(x)))].

(8)

在訓練過程中,生成器GXY試圖生成更像反射系數(shù)的GXY(x),而鑒別器DY則試圖區(qū)分開生成的GXY(x)和目標域中真實的反射系數(shù)y。因此,生成器的目標是使目標函數(shù)最小化,而鑒別器的目標是使目標函數(shù)最大化。對于第一個GAN,其最終目標是minGmaxDYLGAN(GXY,DY,X,Y),同樣地,對于另一個GAN,其目標為minGmaxDYLGAN(GXY,DX,Y,X)。

由上述目標函數(shù)的約束,理論上GAN可以將任意的源域映射至目標域。但是,只有這個約束可能導致生成器只是生成一個具有反射系數(shù)特點的結果,而結構內(nèi)容與源域數(shù)據(jù)并不一樣,即通過訓練只實現(xiàn)了風格的遷移,并沒有保證內(nèi)容不發(fā)生變化。對于任意一個源域數(shù)據(jù)x,需要保證x→GXY(x)→GYX(GXY(x))≈x,這也稱為前向循環(huán)一致性。同樣地,對于任意一個目標域數(shù)據(jù)y,也需要其滿足反向循環(huán)一致性,即:y→GYX(y)→GXY(GYX(y))≈y。因此cycleGAN引入了循環(huán)一致性誤差函數(shù)如下:

(9)

通過上述誤差函數(shù)的約束,對于源域或者目標域的任一數(shù)據(jù),經(jīng)過兩個生成器的循環(huán),最后可以有效還原至原始數(shù)據(jù)特征。但是,當訓練完網(wǎng)絡進行應用時,風格完成遷移但內(nèi)容發(fā)生改變的情況還是會發(fā)生。通過分析可知,應用時只采用了一個生成器作為預測網(wǎng)絡,而循環(huán)一致性誤差函數(shù)沒有對中間變量進行約束。因此,本文中引入身份誤差函數(shù)對中間變量進行約束,其表達式為

(10)

因此對cycleGAN進行訓練時,其整體誤差函數(shù)如下所示:

L(GXY,GXY,DX,DY)=LGAN(GXY,DY,X,Y)+LGAN(GYX,DX,Y,X)+Lcycle(GXY,GYX)+Lidentity(GXY,GYX).

(11)

1.4 訓練集制作

如圖2所示,本文中選取了BP2.5D模型、SEAM模型以及Pluto模型作為基礎模型來制作訓練數(shù)據(jù)。通過有限差分正演和逆時偏移得到三個模型對應的成像結果。同時,基于速度模型生成了反射系數(shù)模型。考慮到反射系數(shù)模型與偏移成像結果的頻帶范圍差異性較大,采用褶積的方式,用雷克子波對反射系數(shù)模型做褶積,作為目標域數(shù)據(jù)。最后采用裁切、重采樣等數(shù)據(jù)增強方式得到了300對訓練數(shù)據(jù)。部分訓練數(shù)據(jù)如圖3所示,第一排為源域的逆時偏移成像數(shù)據(jù),第二排為目標域的褶積結果。

圖2 用于生成訓練集的速度模型和對應的逆時偏移成像結果Fig.2 Velocity model used to generate training set and corresponding reverse time migration imaging results

圖3 部分用于訓練的輸入輸出數(shù)據(jù)Fig.3 Some input and output data used for training

1.5 訓 練

循環(huán)對抗神經(jīng)網(wǎng)絡的訓練流程如圖4所示,對于源域至目標域的單向GAN,訓練時先固定生成器GXY,單獨訓練鑒別器DY,提升鑒別器對生成器輸出結果的鑒別能力,然后再固定鑒別器DY,單獨訓練生成器GXY,使其生成更具有目標域風格的結果,繞開鑒別器的識別。通過這種對抗式的訓練方式,能夠極大提升生成器GXY向目標域映射的能力。同理,對于另一個目標域反映射至源域的單向GAN,也采用同樣的對抗式訓練。當兩個GAN訓練一輪后,將所有損失函數(shù)的值累加,判別其是否滿足期望的誤差大小,若大于期望誤差,則進入下一輪訓練,相反,若小于期望誤差,則結束訓練。最后,將生成器GXY作為期望的預測網(wǎng)絡,用于近似替代Hessian矩陣的逆,對新的逆時偏移成像結果進行預測處理,輸出高質(zhì)量的成像結果。

圖4 循環(huán)對抗神經(jīng)網(wǎng)絡訓練流程Fig.4 Training process of cycle-consistent adversarial neural network

網(wǎng)絡訓練的迭代次數(shù)設置為200次,學習率為0.0002,在進行到第100次迭代后,學習率逐步下降。網(wǎng)絡的循環(huán)一致性誤差權重為10,身份誤差的權重設置為5?；趐ytorch平臺,使用單個NVIDIA RTX5000 GPU訓練,訓練總耗時10 h。

2 模型試算

2.1 Marmousi模型預測

圖5(a)為Marmousi速度模型,對其進行有限差分正演和逆時偏移后,得到的偏移成像結果如圖5(c)所示。由于受到Hessian矩陣的影響,RTM成像結果在淺層有較多的低頻噪音。同時,受模型深部高速屏蔽層的影響,模型底部構造成像分辨率低。圖5(d)是網(wǎng)絡預測的結果,對比圖5(c)可得,預測出的成像結果整體分辨率有了較大提升,反射層位更加精細。淺層的低頻噪音大大減弱,在兩側(cè)有少量殘留,訓練好的網(wǎng)絡所得成像結果中的炮點痕跡不可見。深層的成像質(zhì)量也得到了提升,但是也有一些錯誤成像出現(xiàn),這是因為RTM的成像結果中存在一些雜亂的成像能量,被網(wǎng)絡誤認為層位并恢復出來。

圖5 Marmousi模型預測結果Fig.5 Prediction results of Marmousi model

2.2 Sigsbee2A模型預測

圖6是Sigsbee2A模型的預測結果對比。圖6(a)是速度模型,模型中有橫向展布較大的高速鹽丘。同樣的,對其進行有限差分正演和逆時偏移后,得到圖6(c)所示的逆時偏移成像結果。Sigsbee2A模型的成像剖面同樣存在著明顯的低頻噪音和炮點痕跡,在速度差異較大的分界面附近以及鹽丘內(nèi)部存在著一些偏移假象和噪音,同時,由于鹽丘對地震波的屏蔽作用,鹽下成像能量十分微弱。圖6(d)是網(wǎng)絡對圖6(c)的預測結果,對比圖6(c)可以看出,預測結果的成像分辨率得到了提升,成像能量照明均勻,淺層沒有低頻噪音影響和炮點痕跡消失,鹽下成像質(zhì)量提升明顯,鹽丘邊緣和內(nèi)部的偏移假象也得到了壓制。在成像質(zhì)量提升的同時,預測結果中的繞射點成像效果變差,這可能是因為訓練集的數(shù)據(jù)不包含繞射點的成像數(shù)據(jù),網(wǎng)絡未能學習到與繞射點相關的特征。針對此問題,后期可以在數(shù)據(jù)集中加入具有繞射點相關特征的數(shù)據(jù),提升網(wǎng)絡的泛化能力。

圖6 Sigsbee2A模型預測結果Fig.6 Prediction results of Sigsbee2A model

3 結 論

(1)逆時偏移成像算子可以看作是正演算子的轉(zhuǎn)置,而不是逆,其成像結果為真實反射系數(shù)與Hessian矩陣的褶積,導致逆時偏移成像結果存在低頻噪音、照明不足以及分辨率低等問題。直接求解Hessian矩陣的逆矩陣在計算量和存儲量上十分困難。

(2)循環(huán)對抗神經(jīng)網(wǎng)絡在畫風遷移方面表現(xiàn)優(yōu)秀,其循環(huán)一致性誤差函數(shù)和身份誤差函數(shù)保證了輸出的風格遷移,內(nèi)容不變。通過構建偏移成像數(shù)據(jù)集對網(wǎng)絡進行訓練,訓練后的網(wǎng)絡可以有效地表征Hessian矩陣的逆矩陣。

(3)將訓練好的網(wǎng)絡直接作用于逆時偏移成像結果,預測出的成像剖面低頻噪音得到了顯著的改善,整體成像分辨率有較大提升,高速體屏蔽下的弱成像能量得到了補償,整體成像能量分布比較均衡。雖然該方法在訓練網(wǎng)絡階段比較耗時,但其網(wǎng)絡預測階段的計算量遠小于模型域最小二乘逆時偏移,能夠快速、高效地提高成像質(zhì)量。