融入先驗信息的深度學習地震子波振幅譜估計

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

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

0 引言

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

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

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

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

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

1 方法原理

1.1 褶積模型

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

(1)

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

(2)

1.2 融入先驗信息的深度學習地震子波振幅譜估計算法

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

(3)

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

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

(4)

(5)

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

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

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

圖1 DNN估計子波振幅譜結構圖

模型參數優化的目標函數為

(6)

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

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

通過深度學習網絡，由地震數據得到地震子波振幅譜的過程可描述為

(7)

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

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

1.3 訓練樣本的構建

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

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

表2 地震子波參數設置

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

(8)

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

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

2 應用效果

2.1 模型數據

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

(9)

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

(10)

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

2.1.1 滿足白噪分布的反射系數序列

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

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

2.1.2 滿足非白噪分布的反射系數序列

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

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

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

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

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

圖6 帶通子波振幅譜估計

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

2.2 實際數據

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

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

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

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

圖8 實際資料地震剖面反褶積結果

因此，經本文方法處理的資料分辨率得到明顯提高，橫向連續性也好，實現了對地震資料的保真高分辨率處理。

3 結束語

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

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

(2)本文方法不需要任何假設，不僅能夠對單峰和非單峰地震子波振幅譜進行很好的估計，還避免了各種多項式的參數估計以及有效頻段估計所帶來的計算誤差，具有較好的抗噪性能。模型算例和實際數據的測試都證明了這一點。

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

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

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

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

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