融入先驗信息的深度學(xué)習(xí)地震子波振幅譜估計

黃顏茹 高靜懷* 陳紅靈 姜秀娣

(①西安交通大學(xué)電子與信息學(xué)部信息與通信工程學(xué)院，陜西西安 710049；②中國海洋石油總公司，北京 116503)

0 引言

反射地震學(xué)利用地面接收到的反射數(shù)據(jù)反演地下介質(zhì)的特性[1]。在反射地震學(xué)中，地震子波的估計對于處理地震數(shù)據(jù)具有重要意義[2]。頻域地震子波等于振幅譜與相位譜的乘積，因此地震子波估計包括兩部分內(nèi)容，即振幅譜和相位譜。

現(xiàn)有的地震子波估計方法可以分為兩大類，即確定型方法和統(tǒng)計型方法。Walden等[3]提出一種譜分析方法，利用測井數(shù)據(jù)與地震記錄的功率譜和交互功率譜在一定時窗內(nèi)估計子波；Tan[4]提出一種基于波動方程估計子波振幅譜的方法，該方法假設(shè)子波為最小相位；Edgar等[5]提出了一種匹配子波估計法，利用地震數(shù)據(jù)擬合出子波振幅譜，將相位置零得到初始零相位子波，結(jié)合測井反射系數(shù)可得到合成地震記錄，不斷調(diào)整子波的相位和振幅，使真實值與合成地震記錄達到最佳匹配，此時對應(yīng)的子波即為估計結(jié)果。這些確定型方法利用精確的地震子波振幅譜可加速算法收斂。另一類是統(tǒng)計型方法，這類方法通常假設(shè)反射系數(shù)序列具有白噪非高斯性。Wiggins[6]提出峰態(tài)最大化準則，用于衡量序列非高斯性的參數(shù)；Van der Baan等[7-9]對Wiggins的方法進行改進，提出了非穩(wěn)態(tài)子波估計方法。Ro-binson等[10-11]提出地震記錄的時變褶積模型，即地震子波隨地震波的傳播而變化，不同的旅行時具有不同的子波，Margrave等[12]稱之為非平穩(wěn)褶積模型。作為時變褶積模型的一種近似方法，高靜懷等[13]、Gao等[14-15]和Wang等[16]提出一種非平穩(wěn)地震記錄的自適應(yīng)劃分方法，將非平穩(wěn)地震記錄劃分為近似平穩(wěn)的多個時間間隔。對于這些方法，找到每一個時間區(qū)間的等效子波振幅譜是關(guān)鍵。

在假設(shè)反射系數(shù)序列服從白噪分布的前提下，學(xué)者們提出了一些子波振幅譜估計方法，比如對數(shù)譜平均(LSA)法[17]及譜模擬(SS)法[18]。其中，LSA法不適用于地下介質(zhì)在橫向上不連續(xù)的情況；SS法假定子波振幅譜滿足一種函數(shù)多項式，其中的參數(shù)可通過地震記錄振幅譜與函數(shù)自變量之間的最小二乘擬合確定。在實際應(yīng)用中，選擇合適的函數(shù)形式比較困難，且這些方法對子波振幅譜的形態(tài)有明確限制，得到的是類似于Ricker子波的單峰光滑曲線。此外，大量實際測井數(shù)據(jù)表明，反射系數(shù)序列不滿足隨機分布，其振幅譜符合藍色特征，即反射系數(shù)振幅譜隨頻率增大而逐漸增強，且沿頻率方向強度會發(fā)生變化，呈現(xiàn)出非白色分布的趨勢[19]，這使上述方法在實際中難以應(yīng)用。為了解決該問題，Gao等[20]提出一種壓縮映射算子方法，通過地震數(shù)據(jù)的統(tǒng)計關(guān)系估計子波振幅譜。該方法對非白反射系數(shù)序列十分有效，且不需要預(yù)估函數(shù)形式，但需要選擇合適的擬合頻率范圍。

機器學(xué)習(xí)作為一種數(shù)據(jù)驅(qū)動方法，可在不需要物理理論知識的情況下處理非線性問題。機器學(xué)習(xí)最初是在視覺領(lǐng)域發(fā)展起來的，并逐漸在地球物理領(lǐng)域取得了較好的應(yīng)用效果。目前，人工神經(jīng)網(wǎng)絡(luò)(ANN)及其最新的變體網(wǎng)絡(luò)(即深度學(xué)習(xí))可通過大量訓(xùn)練數(shù)據(jù)學(xué)習(xí)、擬合復(fù)雜非線性函數(shù)，效果非常好。Wang等[21]利用Hopfield神經(jīng)網(wǎng)絡(luò)實現(xiàn)了最小預(yù)測誤差反褶積和地震子波的估計。Chen等[22]提出一種模型驅(qū)動交替耦合的深度網(wǎng)絡(luò)結(jié)構(gòu)，實現(xiàn)了地震子波與反射系數(shù)的反演。唐杰等[23]基于多分辨率的U-Net神經(jīng)網(wǎng)絡(luò)實現(xiàn)了地震數(shù)據(jù)斷層檢測，準確率較高。

本文提出一種融入先驗信息的卷積神經(jīng)網(wǎng)絡(luò)估計地震子波振幅譜的新方法，將子波振幅譜的光滑性作為先驗信息，對地震數(shù)據(jù)進行預(yù)處理，擬合得到的結(jié)果更準確。將本文方法獲得的子波譜與廣泛使用的譜模擬方法[18]獲得的子波譜進行比較，發(fā)現(xiàn)前者對于白色和非白反射系數(shù)均適用，且當子波譜為非單峰(如可控震源子波)時也可準確估計。最后，將本文方法應(yīng)用于一個疊后地震剖面，結(jié)果表明該方法可有效提高地震數(shù)據(jù)的分辨率。

1 方法原理

1.1 褶積模型

常用的地震數(shù)據(jù)處理方法基于傳統(tǒng)的褶積模型，有賴于六個基本假設(shè)[24]。若假設(shè)子波是緊支撐的，則地震記錄在時域中可表示為

(1)

式中：s(t)表示地震記錄，t為時間；w(t)表示地震子波；r(t)為反射系數(shù)序列；n(t)是隨機噪聲;τ表示積分運算中的位移量。其中w(t)、r(t)和n(t)均屬于L2(R)空間。對式(1)兩端同時進行傅里葉變換，可得

(2)

1.2 融入先驗信息的深度學(xué)習(xí)地震子波振幅譜估計算法

通常地，地震子波振幅譜由地震記錄振幅譜經(jīng)非線性運算獲得，因此對式(2)等號兩端同時求振幅譜，即可得到地震記錄的振幅譜

(3)

為了使深度神經(jīng)網(wǎng)絡(luò)向正確的梯度更新方向進行訓(xùn)練，并充分學(xué)習(xí)子波振幅譜信息，可利用子波振幅譜的一些先驗信息。Neidell[25]指出，地震子波振幅譜基本為光滑、有限支集的函數(shù)(正頻率或負頻率)?；谶@一先驗信息，在已知子波振幅譜滿足光滑分布的前提下，對式(3)得到的地震數(shù)據(jù)進行預(yù)處理，并將處理結(jié)果作為輸入數(shù)據(jù)融入神經(jīng)網(wǎng)絡(luò)。

由于地震記錄譜是不規(guī)則的，因此根據(jù)Gao等[20]提出的壓縮映射算子法，通過累積分布函數(shù)中的積分運算，可實現(xiàn)地震數(shù)據(jù)的預(yù)處理，從而得到一個相對光滑的地震子波譜。預(yù)處理的具體過程可表示為

(4)

(5)

式(4)是對式(3)得到的地震記錄振幅譜|S(ω)|進行歸一化處理，然后利用式(5)對歸一化結(jié)果S0(ω)進行積分運算，得到平滑后的地震記錄振幅譜F0(ω)。該預(yù)處理過程具有低通濾波的作用。

預(yù)處理之后，將F0(ω)作為神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)，經(jīng)卷積神經(jīng)網(wǎng)絡(luò)的非線性擬合，最終可估計出較準確的地震子波振幅譜。

本文搭建的深度網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示。網(wǎng)絡(luò)的輸入層為551×1的平滑地震記錄譜，得到的網(wǎng)絡(luò)輸出為551×1的子波振幅譜。需特別注意的是，該網(wǎng)絡(luò)的輸入、輸出可隨輸入長度的大小進行調(diào)整。這里的卷積步長為1，卷積核大小為5×1，采用較大的卷積核可增強對輸入信號的平滑效果。由于線性整流函數(shù)(Rectified Linear Unit, ReLU)相比于Sigmoid和雙曲正切(Tanh)激活函數(shù)的收斂更快，不易產(chǎn)生梯度消失，更適用于深度神經(jīng)網(wǎng)絡(luò)。因此，在除最后一層以外的每層卷積操作后，都使用ReLU函數(shù)進行激活。

圖1 DNN估計子波振幅譜結(jié)構(gòu)圖

模型參數(shù)優(yōu)化的目標函數(shù)為

(6)

損失函數(shù)選用均方誤差(MSE)函數(shù)，該函數(shù)易訓(xùn)練且具有穩(wěn)定解，而且與所要解決問題的目標函數(shù)形式一致。

梯度下降選用自適應(yīng)矩估計(Adaptive Moment Estimation, Adam)梯度最速下降法。

通過深度學(xué)習(xí)網(wǎng)絡(luò)，由地震數(shù)據(jù)得到地震子波振幅譜的過程可描述為

(7)

式中GΘ是非線性擬合算子。該算子通過一個12層深度神經(jīng)網(wǎng)絡(luò)(DNN)在訓(xùn)練過程中不斷優(yōu)化參數(shù)集Θ。增加神經(jīng)網(wǎng)絡(luò)的深度可以增加網(wǎng)絡(luò)的參數(shù)種類，以便更好地學(xué)習(xí)目標特征。但是，過多地增加網(wǎng)絡(luò)層數(shù)會出現(xiàn)梯度消失、過擬合等問題。因此，在同時考慮網(wǎng)絡(luò)訓(xùn)練和測試精度的情況下，將輸入層、12個卷積層(Conv)、歸一化層(BN)、輸出層依次連接，構(gòu)成一個12層的深度卷積網(wǎng)絡(luò)(圖1)。估計子波振幅譜的算法實現(xiàn)流程如表1所示。

表1 融入先驗信息的基于DNN的子波振幅譜估計流程

1.3 訓(xùn)練樣本的構(gòu)建

深度學(xué)習(xí)是一種數(shù)據(jù)驅(qū)動的方法，因此需要生成數(shù)據(jù)集對網(wǎng)絡(luò)進行訓(xùn)練和測試。為增強算法的泛化性與可靠性，采用四種不同類型的地震子波及服從高斯白噪分布和非白噪分布的兩種不同類型的反射系數(shù)序列生成數(shù)據(jù)集。

本文模型數(shù)據(jù)采用的樣本地震子波分別為Ricker子波、廣義地震子波[26]、可控震源子波和帶通子波，參數(shù)設(shè)置見表2。其中，根據(jù)Wang[26]提出的廣義地震子波公式，對參考頻率ω0、時域子波中心點τ0以及分數(shù)值u進行設(shè)置。表中的前兩個子波屬于單峰子波，后兩個屬于非單峰子波，且子波振幅譜在形態(tài)上差異較大。本文暫且不考慮時變子波的情況。

表2 地震子波參數(shù)設(shè)置

樣本中的反射系數(shù)序列分別采用高斯白噪聲隨機數(shù)序列、藍色噪聲序列[27-28]及服從α-stable分布的隨機數(shù)序列[29-31]。第一種反射系數(shù)序列服從高斯白噪分布，后兩者屬于服從非白噪分布的反射系數(shù)序列。其中，藍色噪聲可用白噪聲經(jīng)過一個ARMA(1,1)系統(tǒng)[28]模擬，本文取ARMA(1,1)系統(tǒng)為

(8)

此外，根據(jù)Pavel等[30]提出的產(chǎn)生α-stable分布隨機數(shù)的方法，生成反射系數(shù)序列。

根據(jù)褶積模型將地震子波與反射系數(shù)序列進行褶積計算，可得到合成地震記錄。

2 應(yīng)用效果

2.1 模型數(shù)據(jù)

采用合成地震記錄對DNN進行訓(xùn)練，將估計的子波振幅譜與Rosa等[18]提出的譜模擬方法(SS法)得到的結(jié)果進行比較。定義誤差函數(shù)

(9)

衡量估計得到的子波譜與真實子波譜的相似程度。式中：ASSW表示子波譜；下標“estimated”和“true”分別表示估計值和真實值。誤差函數(shù)的最大值定義為

(10)

下面分析不同類型反射系數(shù)序列情況下得到的估計子波振幅譜。

2.1.1 滿足白噪分布的反射系數(shù)序列

圖2 反射系數(shù)序列服從白噪分布的子波振幅譜估計

上述分析表明，當反射系數(shù)序列服從高斯白噪分布時，利用DNN法和SS法均可從地震記錄振幅譜中提取較準確的子波振幅譜。然而，更多的實驗表明，SS法對地震記錄振幅譜的有效頻率及位置十分敏感，且在真實子波未知的情況下，該方法中多項式里的擬合參數(shù)較難確定，為子波振幅譜的估計帶來很大限制。

2.1.2 滿足非白噪分布的反射系數(shù)序列

圖3 反射系數(shù)序列服從藍噪分布時的子波振幅譜估計

圖4 反射系數(shù)序列服從α-stable分布時的子波振幅譜估計結(jié)果

2.1.3 非單峰子波振幅譜的估計

上述測試主要針對單峰子波振幅譜的估計，本節(jié)分析本文算法對于非單峰子波振幅譜的估計效果。圖5和圖6分別是可控震源子波和帶通子波兩種非單峰子波振幅譜的估計結(jié)果。

圖5 可控震源子波振幅譜估計

圖6 帶通子波振幅譜估計

綜上所述，無論是高斯白噪，還是藍色噪聲或α-stable分布的反射系數(shù)序列，本文提出的DNN法都能較準確地估計地震子波振幅譜，尤其能夠?qū)崿F(xiàn)非單峰子波振幅譜的準確估計，而現(xiàn)有的其他地震子波振幅譜估計方法則無法實現(xiàn)。

2.2 實際數(shù)據(jù)

為了檢驗本文提出的DNN方法對實際數(shù)據(jù)的處理效果，對中國海洋石油總公司的一個實際疊后地震剖面進行高分辨處理。該疊后地震剖面的時間采樣點為551，采樣間隔為2ms，共1001個地震道。對于每個地震道取其中150個采樣點近似為平穩(wěn)的地震數(shù)據(jù)。將實際資料作為測試數(shù)據(jù)，并利用基于模型數(shù)據(jù)訓(xùn)練好的網(wǎng)絡(luò)進行子波提取效果測試。隨機抽取了實際資料中的4道數(shù)據(jù)，進行子波振幅譜估計。

圖7展示不同地震道的估計子波振幅譜?？梢园l(fā)現(xiàn)，DNN法從實際地震資料得到的擬合包絡(luò)合理、準確，能夠較好地描述地震子波形態(tài)。將該子波振幅譜用于維納反褶積[33]，可得到高分辨率結(jié)果(圖8)。本文只采用子波振幅進行反褶積，即零相位反褶積。

圖7 實際資料DNN法估計的地震子波振幅譜

圖8是該實際地震剖面的高分辨率處理結(jié)果，其中第270道為測井數(shù)據(jù)，用以驗證處理結(jié)果的準確性。由圖8b可以看出，采用DNN法估計的子波振幅譜經(jīng)反褶積處理后，資料的縱向分辨率得到明顯提高，與測井數(shù)據(jù)有較好的一致性，尤其是黃色方框區(qū)域，展現(xiàn)出較好的橫向連續(xù)性。從圖8c所示的傳統(tǒng)方法處理結(jié)果可見，剖面上存在明顯的橫向不連續(xù)現(xiàn)象及噪聲，對實際資料的縱向分辨率沒有明顯改善。

圖8 實際資料地震剖面反褶積結(jié)果

因此，經(jīng)本文方法處理的資料分辨率得到明顯提高，橫向連續(xù)性也好，實現(xiàn)了對地震資料的保真高分辨率處理。

3 結(jié)束語

本文提出一種融入先驗信息的深度學(xué)習(xí)地震子波振幅譜估計方法。首先通過已知地震子波振幅譜光滑這一先驗信息，對地震數(shù)據(jù)進行預(yù)處理；然后，經(jīng)過一個12層的卷積神經(jīng)網(wǎng)絡(luò)，對輸入地震數(shù)據(jù)進行非線性擬合；最后，提取地震子波振幅譜。得到以下結(jié)論。

(1)在本文DNN子波振幅譜估計方法中，需要對地震數(shù)據(jù)進行必要的預(yù)處理，因為實際的測井數(shù)據(jù)既包含噪聲，也包含測量誤差[34]。

(2)本文方法不需要任何假設(shè)，不僅能夠?qū)畏搴头菃畏宓卣鹱硬ㄕ穹V進行很好的估計，還避免了各種多項式的參數(shù)估計以及有效頻段估計所帶來的計算誤差，具有較好的抗噪性能。模型算例和實際數(shù)據(jù)的測試都證明了這一點。

(3)由于傳統(tǒng)的地震子波振幅譜估計方法是基于統(tǒng)計學(xué)提出的，因此對于不同分布特征的反射系數(shù)以及子波的形態(tài)有一定限制；深度學(xué)習(xí)作為一種數(shù)據(jù)驅(qū)動方法，不依賴于地震數(shù)據(jù)的統(tǒng)計特性，因此不需要對反射系數(shù)類型作假設(shè)。

尚需指出，本文方法雖能準確地估計地震子波振幅譜，但仍具有以下幾方面發(fā)展空間：

(1)本文提出的融合先驗信息的深度神經(jīng)網(wǎng)絡(luò)是一種有監(jiān)督的機器學(xué)習(xí)方法，其估計精度依賴于訓(xùn)練樣本中地震子波種類的個數(shù)，若訓(xùn)練樣本中地震子波的類型不夠充足，則訓(xùn)練的模型泛化性不能得到保證。筆者單位擬與油田合作，根據(jù)多塊實際數(shù)據(jù)工區(qū)重新建立子波樣本標簽，針對實際中復(fù)雜的地震數(shù)據(jù)，以使該方法具有更好的泛化性；

(2)本文方法無需考慮時變子波振幅譜的提取，因此對于非平穩(wěn)地震數(shù)據(jù)，可采用非平穩(wěn)地震數(shù)據(jù)劃分[16]等方法對實際地震數(shù)據(jù)先作分段處理，在后續(xù)的研究中則需考慮時變性；

(3)反射系數(shù)序列服從非白分布情況下的子波振幅譜估計對于Q值的提取也至關(guān)重要，需做進一步的研究。