張聯(lián)海， 王 璐，2??， 鄭志超， 孟凡順，2

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 山東 青島 266100； 2.海底科學(xué)與探測技術(shù)教育部重點實驗室， 山東 青島 266100)

反射地震勘探是油氣勘探的主要方法之一，高信噪比和高分辨率的地震資料，是準(zhǔn)確解釋地層深部地質(zhì)構(gòu)造的必要條件之一。反褶積是提高地震資料分辨率最常用的方法，在高分辨率地震資料處理中起著至關(guān)重要的作用。

根據(jù)褶積模型，反射地震信號可以表示為地層反射系數(shù)和震源子波的褶積。由觀測的反射地震信號推測地層反射系數(shù)，理想條件下，真實的地層反射系數(shù)可以得到恢復(fù)，然而實際勘探中，震源子波的帶限的，這導(dǎo)致反褶積是一個不適定的反問題，直接求解難以得到穩(wěn)定的反演結(jié)果。合適的正則化約束可以加強(qiáng)反問題求解的穩(wěn)定性[1-2]，基于地層的分層性質(zhì)，只有較大的地層反射系數(shù)才影響反射地震信號，因此，反射系數(shù)在時間上是稀疏的，使用反射系數(shù)的稀疏性質(zhì)作為正則化約束的反褶積方法即為稀疏反褶積[3-4]?；诜瓷湎禂?shù)的稀疏性質(zhì)，Wang[5]建立了L1正則化反褶積模型，通過非單調(diào)梯度下降算法求解反射系數(shù)，然后采用遞推方法求解地層的波阻抗剖面。Gholami和Sacchi[6]將反褶積問題轉(zhuǎn)化為無約束優(yōu)化問題，此優(yōu)化問題采用Lp范數(shù)作為擬合項，L1范數(shù)作為正則化項，然后利用分裂Bregman法[7]求解，此方法在準(zhǔn)確反演稀疏反射系數(shù)的同時，可以有效壓制高斯和非高斯噪音。稀疏反褶積通常會導(dǎo)出一個Topelitz矩陣的線性系統(tǒng)[8]，基于Topelitz矩陣的稀疏分解，Wang等[9]提出了一個多道稀疏脈沖反卷積方法，將震源子波和地層反射系數(shù)都作為反演對象，二者交替求解，從而達(dá)到同時反演震源子波和反射系數(shù)的目的，此方法緩解了傳統(tǒng)方法中出現(xiàn)的空間不連續(xù)、信號受噪音影響、子波估計不準(zhǔn)確等問題。反褶積問題的求解還需使用合適的最優(yōu)化算法，如共軛梯度法[10]，迭代閾值收縮法(ISTA)[11]等，最近，國內(nèi)學(xué)者潘樹林等[12]提出自適應(yīng)步長的快速迭代閾值收縮法(FISTA)用于求解稀疏脈沖反褶積問題，F(xiàn)ISTA算法[13]是在ISTA算法的基礎(chǔ)上使用Nesterov[14-15]加速策略，使得算法具有超線性收斂速率。上述稀疏反褶積方法均屬于正則化迭代方法，對反褶積問題施加稀疏正則化約束，降低反演結(jié)果的不穩(wěn)定性，然后通過合適的優(yōu)化算法迭代求解。求解過程通常需要反復(fù)迭代，每次迭代均需進(jìn)行多次正演，因此，其計算代價較高。另外，正則化迭代方法還存在正則化參數(shù)選取困難的問題，正則化參數(shù)的作用是平衡擬合項和正則化項的權(quán)重，此參數(shù)的選取需進(jìn)行多次試驗，才能得到較合適的值。最后，正則化迭代方法的反演結(jié)果不精確，反演結(jié)果中較大的反射系數(shù)振幅有所損失，小的反射系數(shù)則難以得到恢復(fù)。

近幾年，深度學(xué)習(xí)[16]方法因其強(qiáng)大的學(xué)習(xí)能力受到廣泛關(guān)注，其中回歸型深度神經(jīng)網(wǎng)絡(luò)在反問題的求解方面也取得了大量進(jìn)展，如信號去噪[17]，圖像反卷積[18]，圖像重建[19]等。國內(nèi)外學(xué)者已將深度學(xué)習(xí)推廣到地球物理反演中。Kim和Nakata[20]以地震反射系數(shù)反演為例探討了地球物理反演和機(jī)器學(xué)習(xí)之間的聯(lián)系和區(qū)別。從數(shù)學(xué)的角度，二者的目標(biāo)均為解決不適定的非線性反問題，且都需應(yīng)用合適的優(yōu)化算法。地球物理反演首先需要根據(jù)一定的物理規(guī)律建立反演模型，然后加入合理的先驗假設(shè)降低反演問題的不適定性，最后利用優(yōu)化算法求解。機(jī)器學(xué)習(xí)則屬于數(shù)據(jù)驅(qū)動類方法，通過大量的數(shù)據(jù)訓(xùn)練，從數(shù)據(jù)中學(xué)習(xí)輸入到輸出的非線性關(guān)系，然后利用學(xué)習(xí)到的非線性關(guān)系對測試數(shù)據(jù)進(jìn)行預(yù)測，從而達(dá)到求解反問題的目的。Yu等[21]利用DCNN對地震信號進(jìn)行去噪，傳統(tǒng)的信號去噪方法多是基于模型的，而DCNN則是數(shù)據(jù)驅(qū)動型方法，無需任何先驗假設(shè)，經(jīng)過訓(xùn)練的網(wǎng)絡(luò)可以自適應(yīng)的去除隨機(jī)噪聲、線性噪聲、多次波等多種干擾噪音。然而上述工作使用的仍然是“黑盒式”的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，網(wǎng)絡(luò)的可解釋性較差。近幾年，對于反問題的求解，眾多研究者開始關(guān)注深度學(xué)習(xí)網(wǎng)絡(luò)與傳統(tǒng)正則化迭代算法之間的關(guān)系，并由此構(gòu)建高效的、可解釋的深度神經(jīng)網(wǎng)絡(luò)。Gregor和LeCun[22]研究了ISTA算法和共享分層神經(jīng)網(wǎng)絡(luò)之間的相似性，并證明了多層神經(jīng)網(wǎng)絡(luò)可以作為一種快速的近似稀疏編碼器。類似的，Jin等[23]提出一個深度卷積神經(jīng)網(wǎng)絡(luò)用于稀疏圖像重建，并解釋了此網(wǎng)絡(luò)與傳統(tǒng)迭代算法的關(guān)系。

受到上述研究工作的啟發(fā)，針對傳統(tǒng)正則化迭代算法存在的計算代價高，正則化參數(shù)難以選取，反演結(jié)果不精確等問題，本文構(gòu)建一個深度卷積神經(jīng)網(wǎng)絡(luò)，用于地震反射信號的反褶積。首先使用訓(xùn)練數(shù)據(jù)集對此網(wǎng)絡(luò)進(jìn)行訓(xùn)練，使得其可以從數(shù)據(jù)中自動學(xué)得輸入的地震反射信號與輸出的地震反射系數(shù)之間的映射關(guān)系，然后利用訓(xùn)練好的網(wǎng)絡(luò)對反射系數(shù)進(jìn)行預(yù)測。所建網(wǎng)絡(luò)融合了多尺度分解和殘差學(xué)習(xí)[24]等技巧，提高了網(wǎng)絡(luò)的表達(dá)能力。

1 方法原理

1.1 稀疏反褶積及ISTA算法

由地震信號的褶積模型知，一道地震信號可表示為震源子波與反射系數(shù)的褶積，即

d(t)=ω(t)×r(t)+n(t)。

(1)

式中：d(t)為觀測到的疊后地震信號；ω(t)為震源子波；n(t)表示噪聲；×表示褶積運(yùn)算。(1)式還可寫成矩陣方程的形式，即

d=Wr+n。

(2)

式中：d∈RN；r∈RN和n∈RN分別表示d(t)，r(t)和n(t)的離散采樣；W∈RN×N為Toeplitz矩陣，且矩陣元素Wi,j=ωi-j+1；ωi-j+1是震源子波ω(t)的第i-j+1個采樣點。反褶積即從觀測到的反射地震信號d中推測反射系數(shù)r。由于子波是帶限的，由(2)式反演反射系數(shù)r是一個不適定的反問題，直接求解難以得到穩(wěn)定的反演結(jié)果，因此引入合適的正則化約束是必要的?？紤]到的反射系數(shù)的稀疏性質(zhì)，引入L1正則化，可得反演反射系數(shù)的L1正則化模型如下

(3)

(4)

上式中第一項通常稱為擬合項，本文采用L2范數(shù)衡量擬合度，第二項為正則化項。(3)式或(4)式的求解可采用ISTA算法，ISTA為迭代優(yōu)化算法，第k+1步的迭代格式為

rk+1=Τλtk(rk-tkWΤ(Wrk-d))。

(5)

Τα(r)i=(|ri|-α)+sgn(ri)。

sgn為符號函數(shù)。進(jìn)一步的，(5)式可寫為

rk+1=Τα(Ark+Bd)。

(6)

式中：α=λtk；A=I-tkWΤW；B=tkWΤ。注意到，上式中正則化參數(shù)λ，步長tk均需人工選擇，實際實驗時，需反復(fù)嘗試不同的數(shù)值；此外ISTA為迭代算法，需多次迭代才能最終得到最優(yōu)解，計算代價較高。

1.2 深度卷積神經(jīng)網(wǎng)絡(luò)方法

使用深度學(xué)習(xí)方法求解(4)式的基本思想是建立反射地震信號d和(4)式最優(yōu)解r*之間的映射關(guān)系，其數(shù)學(xué)模型如下：

r*=Net(d,θ)。

(7)

式中：r*=argminrE(r)；Net表示神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)；θ為待學(xué)習(xí)的網(wǎng)絡(luò)參數(shù)集合。網(wǎng)絡(luò)通過極小化下式學(xué)習(xí)參數(shù)θ：

(8)

(9)

本文所建網(wǎng)絡(luò)模型參考了U型網(wǎng)絡(luò)[26]的結(jié)構(gòu)(U型網(wǎng)絡(luò)最初主要用于醫(yī)學(xué)圖像分割)，網(wǎng)絡(luò)結(jié)構(gòu)呈對稱形式，分為編碼和解碼兩個部分。編碼部分使用多個卷積核對輸入地震信號提取不同的特征，并逐層增加卷積核的數(shù)量，這種多尺度分解方法是DCNN常用的方法。然后應(yīng)用最大池化操作，降低輸入數(shù)據(jù)的維度，提取主要特征，通過多次池化操作，原始數(shù)據(jù)將被映射到較低的維度，同時由于多個卷積核的作用，數(shù)據(jù)的主要特征即得到提取。解碼部分與編碼部分對稱，使用上采樣依次增加數(shù)據(jù)的維度，同時利用多個卷積核增加網(wǎng)絡(luò)的表達(dá)能力。為避免網(wǎng)絡(luò)訓(xùn)練時出現(xiàn)梯度彌散的問題，引入了殘差學(xué)習(xí)，即將數(shù)據(jù)越層連接，網(wǎng)絡(luò)實際學(xué)習(xí)的是輸入與輸出之差。所提網(wǎng)絡(luò)模型的具體參數(shù)為：卷積核長度為5，步長為1，各層卷積核數(shù)量分別為32、64、128、256、128、64、32、1，最大池化長度為2，步長為1，網(wǎng)絡(luò)結(jié)構(gòu)的具體細(xì)節(jié)如圖1所示。

圖1 U型網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.1 U-net architecture

2 數(shù)值實驗

數(shù)值實驗包括兩個方面，其分別使用模擬數(shù)據(jù)和實際數(shù)據(jù)，兩部分實驗均與ISTA方法對比。所有實驗均在個人電腦上進(jìn)行，電腦裝配Intel i7-4790 CPU, 主頻為3.60 GHz，內(nèi)存為32 GB，GPU為NVIDIA GeForce GTX 1070。

2.1 模擬實驗

模擬實驗采用人工模擬的數(shù)據(jù)進(jìn)行對比試驗。使用DCNN方法進(jìn)行反褶積之前，需要訓(xùn)練集對網(wǎng)絡(luò)進(jìn)行訓(xùn)練。訓(xùn)練數(shù)據(jù)由兩部分組成，即地震信號和反射系數(shù)。離散反射系數(shù)長度為256，采樣率為1 ms，于是可使用稀疏隨機(jī)向量代替，向量中每個反射系數(shù)幅值范圍為[-1, 1]，位置是隨機(jī)的。反射地震信號由反射系數(shù)和主頻為35 Hz的Ricker子波褶積得到。為使網(wǎng)絡(luò)得到充分訓(xùn)練，本文使用10 000組訓(xùn)練數(shù)據(jù)，1 000組測試數(shù)據(jù)用于測試DCNN方法的有效性，其中一組如圖2所示。網(wǎng)絡(luò)結(jié)構(gòu)已由上節(jié)給出，訓(xùn)練時，為防止網(wǎng)絡(luò)出現(xiàn)過擬合現(xiàn)象，每層均舍棄20%的參數(shù)。為加快訓(xùn)練速度，每次輸入64組數(shù)據(jù)，采用步長為1×e-4的Adam算法[27]進(jìn)行訓(xùn)練，訓(xùn)練200次即終止。經(jīng)訓(xùn)練后的網(wǎng)絡(luò)即可用于地震信號的反褶積，反褶積結(jié)果的相對誤差由下式給出

(10)

圖2 訓(xùn)練數(shù)據(jù)Fig.2 Training data

圖3 ISTA及DCNN所得結(jié)果Fig.3 Results of ISTA and DCNN

進(jìn)一步的，為測試本文所提方法的有效性，設(shè)計一個四層的楔形模型，模型包含50道數(shù)據(jù)，采樣率為1 ms，采樣點數(shù)為256，如圖4所示。采用主頻為35 Hz的Ricker子波與反射系數(shù)褶積得到地震記錄，如圖5所示。反演時，使用訓(xùn)練好的深度卷積神經(jīng)網(wǎng)絡(luò)，將正演得到的地震記錄作為輸入，網(wǎng)絡(luò)的輸出即為反演得到的反射系數(shù)，如圖6所示，其相對誤差Err=0.046 3。圖7對比了圖6中第一道反射系數(shù)與真實反射系數(shù)的差別，其相對誤差Err=0.085 6。由圖6和圖7可知，DCNN方法反演得到的結(jié)果準(zhǔn)確度較高，反射系數(shù)振幅真實反射系數(shù)振幅大小一致。作為對比，采用ISTA算法迭代求解(4)式，迭代格式由(5)式給出，取正則化參數(shù)λ=0.9。ISTA算法的反演結(jié)果分別由圖8、圖9給出，由圖可知，ISTA算法可以壓縮子波，反演反射系數(shù)的目的，然而由圖8可知ISTA算法對子波的壓縮并不充分，圖8中楔形左上角處(即紅色方框處)的反射系數(shù)未能得到正確的反演，其相對誤差為Err=0.734 8，圖9亦顯示了ISTA算法未能反演出相距較近的兩個反射系數(shù)，且反射系數(shù)振幅損失較大，其相對誤差為Err=0.798 0。綜合對比圖4至圖9可知，ISTA方法反演的反射系數(shù)振幅損失較大，對于子波的壓縮并不充分，并且不能有效分辨距離較近的兩個反射系數(shù)，這使得楔形模型的左上角處出現(xiàn)了誤差。而DCNN方法可以較準(zhǔn)確的反演反射系數(shù)，即使相距較近的兩個反射系數(shù)也能準(zhǔn)確反演。

圖4 楔形模型Fig.4 Wedge model

圖5 楔形模型的正演記錄Fig.5 Forward modeling record of wedge model

(相對誤差Err=0.046 3。Relative errorErr=0.046 3.)圖6 DCNN所得反射系數(shù)Fig.6 Reflection coefficient obtained by DCNN

(相對誤差Err=0.085 6。Relative errorErr=0.085 6.)圖7 DCNN所得的單道反射系數(shù)Fig.7 Single channel reflection coefficient obtained by DCNN

(相對誤差Err=0.734 8。Relative errorErr=0.734 8.)圖8 ISTA所得反射系數(shù)Fig.8 Reflection coefficient obtained by ISTA

(相對誤差Err=0.798 0。Relative error Err=0.798 0.)圖9 ISTA所得的單道反射系數(shù)Fig.9 Single channel reflection coefficient obtained by ISTA

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

本部分實驗采用實際數(shù)據(jù)檢驗所提方法的有效性。截取的某工區(qū)部分實際數(shù)據(jù)如圖10所示，此數(shù)據(jù)體包括53道地震反射信號，時間長度為256 ms，采樣率為1 ms。利用統(tǒng)計方法估計的震源子波如圖11所示。得到震源子波后，使用隨機(jī)產(chǎn)生反射系數(shù)序列與震源子波褶積得到模擬地震信號，由反射系數(shù)序列和模擬地震信號組成DCNN方法的訓(xùn)練數(shù)據(jù)集，訓(xùn)練輸入為模擬地震信號，輸出為反射系數(shù)序列。訓(xùn)練好的網(wǎng)絡(luò)模型即可直接用于反褶積，結(jié)果如圖12所示，由圖可知，相比地震數(shù)據(jù)剖面，反褶積得到的反射系數(shù)剖面中子波寬度得到壓縮，分辨率較高，不僅大的反射系數(shù)得到恢復(fù)，較小的反射系數(shù)也被反演得到。作為對比，仍采用ISTA算法迭代求解(4)式，正則化參數(shù)λ=0.9時的計算結(jié)果如圖13所示，由圖可知，相比地震數(shù)據(jù)剖面，反射系數(shù)剖面的分辨率有所提高，然而子波寬度未能得到充分壓縮，且僅有較大的反射系數(shù)得到恢復(fù)，較小的反射系數(shù)并未出現(xiàn)。

圖10 實際地震數(shù)據(jù)Fig.10 Actual seismic data

圖11 估計的子波Fig.11 Estimated wavelet

圖12 DCNN方法的反褶積結(jié)果Fig.12 Deconvolution results of DCNN method

圖13 ISTA算法的反褶積結(jié)果Fig.13 Deconvolution results of ISTA algorithm

圖14給出了DCNN方法和ISTA方法反褶積結(jié)果的單道對比結(jié)果，兩種方法均可以反演出較大振幅的反射系數(shù)，所得結(jié)果的主要差異由三個紅色方框標(biāo)記，三個紅色方框中，DCNN方法可以較好的反演出距離較近的弱反射系數(shù)，而ISTA方法未能反演出反射系數(shù)。此外，相比DCNN方法所得結(jié)果，ISTA方法得到反射系數(shù)振幅損失較大。

圖14 DCNN和ISTA的反褶積結(jié)果對比Fig.14 Comparison of deconvolution results of DCNN and ISTA

3 結(jié)語

本文構(gòu)建了一個DCNN模型，使用合適的訓(xùn)練集訓(xùn)練后，即可用于地震反射信號的反褶積。DCNN屬于端到端的數(shù)據(jù)驅(qū)動方法，進(jìn)行反演時，直接輸入地震信號，輸出即為反褶積結(jié)果。傳統(tǒng)神經(jīng)網(wǎng)絡(luò)僅僅是建立了輸入端到輸出端的映射關(guān)系，網(wǎng)絡(luò)結(jié)構(gòu)的可解釋性較差，本文闡述了所建的深度卷積神經(jīng)網(wǎng)絡(luò)模型與ISTA算法結(jié)構(gòu)上的相似性，網(wǎng)絡(luò)具有一定的可解釋性。數(shù)值實驗從模擬數(shù)據(jù)和實際數(shù)據(jù)兩個方面驗證了本文所提方法的有效性，與ISTA算法對比，DCNN方法所需時間少，反褶積結(jié)果更加準(zhǔn)確，對于距離較近的弱振幅亦能較準(zhǔn)確反演，反演得到的反射系數(shù)剖面分辨率更高。