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

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

0 引言

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

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

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

機(jī)器學(xué)習(xí)作為一種數(shù)據(jù)驅(qū)動(dòng)方法，可在不需要物理理論知識(shí)的情況下處理非線性問題。機(jī)器學(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ò)實(shí)現(xiàn)了最小預(yù)測誤差反褶積和地震子波的估計(jì)。Chen等[22]提出一種模型驅(qū)動(dòng)交替耦合的深度網(wǎng)絡(luò)結(jié)構(gòu)，實(shí)現(xiàn)了地震子波與反射系數(shù)的反演。唐杰等[23]基于多分辨率的U-Net神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)了地震數(shù)據(jù)斷層檢測，準(zhǔn)確率較高。

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

1 方法原理

1.1 褶積模型

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

(1)

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

(2)

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

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

(3)

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

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

(4)

(5)

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

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

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

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

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

(6)

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

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

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

(7)

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

表1 融入先驗(yàn)信息的基于DNN的子波振幅譜估計(jì)流程

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

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

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

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

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

(8)

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

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

2 應(yīng)用效果

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

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

(9)

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

(10)

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

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

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

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

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

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

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

2.1.3 非單峰子波振幅譜的估計(jì)

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

圖5 可控震源子波振幅譜估計(jì)

圖6 帶通子波振幅譜估計(jì)

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

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

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

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

圖7 實(shí)際資料DNN法估計(jì)的地震子波振幅譜

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

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

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

3 結(jié)束語

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

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

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

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

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

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

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

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