基于深度卷積神經(jīng)網(wǎng)絡(luò)的稀疏反褶積方法?

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

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

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

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

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

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

1 方法原理

1.1 稀疏反褶積及ISTA算法

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

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

(1)

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

d=Wr+n。

(2)

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

(3)

(4)

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

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

(5)

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

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

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

(6)

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

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

使用深度學(xué)習(xí)方法求解(4)式的基本思想是建立反射地震信號(hào)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ò)通過(guò)極小化下式學(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)呈對(duì)稱(chēng)形式，分為編碼和解碼兩個(gè)部分。編碼部分使用多個(gè)卷積核對(duì)輸入地震信號(hào)提取不同的特征，并逐層增加卷積核的數(shù)量，這種多尺度分解方法是DCNN常用的方法。然后應(yīng)用最大池化操作，降低輸入數(shù)據(jù)的維度，提取主要特征，通過(guò)多次池化操作，原始數(shù)據(jù)將被映射到較低的維度，同時(shí)由于多個(gè)卷積核的作用，數(shù)據(jù)的主要特征即得到提取。解碼部分與編碼部分對(duì)稱(chēng)，使用上采樣依次增加數(shù)據(jù)的維度，同時(shí)利用多個(gè)卷積核增加網(wǎng)絡(luò)的表達(dá)能力。為避免網(wǎng)絡(luò)訓(xùn)練時(shí)出現(xiàn)梯度彌散的問(wèn)題，引入了殘差學(xué)習(xí)，即將數(shù)據(jù)越層連接，網(wǎng)絡(luò)實(shí)際學(xué)習(xí)的是輸入與輸出之差。所提網(wǎng)絡(luò)模型的具體參數(shù)為：卷積核長(zhǎng)度為5，步長(zhǎng)為1，各層卷積核數(shù)量分別為32、64、128、256、128、64、32、1，最大池化長(zhǎng)度為2，步長(zhǎng)為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í)驗(yàn)

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

2.1 模擬實(shí)驗(yàn)

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

(10)

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

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

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

圖4 楔形模型Fig.4 Wedge model

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

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

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

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

(相對(duì)誤差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í)際數(shù)據(jù)

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

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

圖11 估計(jì)的子波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é)果的單道對(duì)比結(jié)果，兩種方法均可以反演出較大振幅的反射系數(shù)，所得結(jié)果的主要差異由三個(gè)紅色方框標(biāo)記，三個(gè)紅色方框中，DCNN方法可以較好的反演出距離較近的弱反射系數(shù)，而ISTA方法未能反演出反射系數(shù)。此外，相比DCNN方法所得結(jié)果，ISTA方法得到反射系數(shù)振幅損失較大。

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

3 結(jié)語(yǔ)

本文構(gòu)建了一個(gè)DCNN模型，使用合適的訓(xùn)練集訓(xùn)練后，即可用于地震反射信號(hào)的反褶積。DCNN屬于端到端的數(shù)據(jù)驅(qū)動(dòng)方法，進(jìn)行反演時(shí)，直接輸入地震信號(hào)，輸出即為反褶積結(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í)驗(yàn)從模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)兩個(gè)方面驗(yàn)證了本文所提方法的有效性，與ISTA算法對(duì)比，DCNN方法所需時(shí)間少，反褶積結(jié)果更加準(zhǔn)確，對(duì)于距離較近的弱振幅亦能較準(zhǔn)確反演，反演得到的反射系數(shù)剖面分辨率更高。