基于迭代收縮閾值網絡的地震數(shù)據(jù)重構研究

范 帥，邢 磊,2，李倩倩,2

(1.中國海洋大學 海底科學與探測技術教育部重點實驗室，山東 青島 266100；2.海洋國家實驗室 海洋礦產資源評價與探測技術功能實驗室，山東 青島 266071)

1 引 言

地震數(shù)據(jù)重建是勘探地震學中一個重要的研究方向，近幾十年來學者們已經提出了多種地震數(shù)據(jù)重構的方法[1,2]，如插值法。地震數(shù)據(jù)插值，也可以稱為地震數(shù)據(jù)恢復，是地震資料處理中極其重要的一環(huán)。傳統(tǒng)的地震數(shù)據(jù)插值方法可大致包含三大類：基于預測濾波器的插值算法，基于波動方程的插值算法以及基于信號分析處理的插值算法。Donoho等人(2006年)[3]提出壓縮感知理論，即信號的采樣頻率可以低于Nyquist頻率，并在采樣后通過合適的恢復算法能較為完整的恢復原信號。壓縮感知理論可以對數(shù)據(jù)同時進行采樣和壓縮，這樣可以去除大量冗余數(shù)據(jù)并且釋放有限的儲存資源，壓縮感知理論的提出推動了地震數(shù)據(jù)重構方法的發(fā)展。

Wang等人(2011年)[4]，將地震波場重構問題看作一個壓縮感知問題，提出一種基于小波變換的分段隨機子采樣方法，通過求解L1范數(shù)來解決壓縮感知問題，提高了波場恢復的質量。曹靜杰等人(2012年)[5]在Curvelet域利用梯度投影的算法求解L0范數(shù)約束下的壓縮感知問題。Gou等人(2014年)[6]通過Seislet變換對復雜地震數(shù)據(jù)進行稀疏后，應用Bregman迭代算法求解壓縮感知問題中的混合范數(shù)優(yōu)化問題。作為一種多尺度稀疏變換，Seislet變換可以將地震信號進行有效的壓縮，因此對變換域的選擇也尤為重要。隨著字典學習的火熱，周亞同等人(2014年)[7]將K-SVD字典學習方法應用于地震數(shù)據(jù)重建上，這個方法的應用是在壓縮感知框架下的。該方法通過用K-SVD字典訓練大量的地震樣本數(shù)據(jù)，并在訓練完成后得到超完備字典，并采用正則化正交匹配追蹤(ROMP,Regularization Orthogonal Matching Pursuit)來對地震數(shù)據(jù)實現(xiàn)重建，可以實現(xiàn)地震數(shù)據(jù)的自適應稀疏擴展。但基于壓縮感知的地震數(shù)據(jù)重構方法在不同變換域中的重構效果有很大區(qū)別，重構效果不穩(wěn)定。

近幾年深度學習在圖像處理和圖像修復領域體現(xiàn)出了非常好的效果， Cui等人(2014年)[8]通過對圖像使用降采樣后，將高分辨率圖像降采樣成低分辨率圖像，將降采樣后的低分辨率圖像作為訓練數(shù)據(jù)集輸入進深度卷積神經網絡中進行訓練，之后使用訓練過的網絡對圖像進行重建。Wang等人(2016年)[9]基于深度學習在圖像處理方面的優(yōu)勢，使用高質量的核磁共振圖像對神經網絡進行訓練來加速磁共振成象(MRI,Magnetic Resonance Imaging)的方法。同年，Sun等人[10]在壓縮感知可以快速磁共振成象的基礎上，通過將ADMM(ADMM,Alternating Direction Method of Multipliers)算法與深度神經網絡結合，將網絡的結構變?yōu)锳DMM算法的展開形式來進行磁共振成象的壓縮感知重建。Zhang等人(2018年)[11]將壓縮感知凸優(yōu)化算法中的迭代軟閾值算法的每一步迭代過程用深度神經網絡的方法展開，在自然圖像重建和磁共振成像(MRI)中獲得了良好的效果。Mandelli等人(2018)[12]構建了一種卷積自動編碼器的特定類型的卷積神經網絡來對損壞的地震數(shù)據(jù)進行插值和去噪。

受到圖像處理和修復的啟發(fā)，以及傳統(tǒng)的地震數(shù)據(jù)重構方法存在的運算量大、重構成本大的問題，本文通過將迭代收縮軟閾值算法和深度神經網絡結合起來，將其應用在地震數(shù)據(jù)的重構方法中，通過使用大數(shù)據(jù)集對網絡進行訓練，實現(xiàn)端到端的映射，網絡中的參數(shù)都可以通過數(shù)據(jù)集自動學習并更新，無需人為手動設置每步參數(shù)。實驗結果證明，該方法可以高效、精確地恢復地震信號。

2 迭代收縮閾值算法

迭代收縮閾值算法(ISTA,Iterative Shrinkage-Thresholding Algorithm)最早可以追溯到2004年，Daubechies等人(2004年)[13]提出了一種迭代算法，迭代收縮軟閾值算法由于其簡單性、迭代速度快、迭代所需的算法復雜度小，因此被廣泛的研究。

首先，地震數(shù)據(jù)重構問題可以同構構建一個線性方程組解決

y=Φx

(1)

其中，x為地震信號；Φ為觀測矩陣；y為觀測信號。通過使用最小二乘法可以求解上述問題：

(2)

在實際情況下，x一般都是不稀疏的，無論x是自然信號或是地震信號，它需要通過某種變換實現(xiàn)稀疏，例如稀疏基，x可以通過稀疏基來實現(xiàn)稀疏，稀疏基矩陣用Ψ來表示，s為稀疏變換后的稀疏系數(shù)。因此，公式(1)也就變成了：

y=ΦΨs

(3)

在式(3)中，y和Φ以及Ψ都是已知的值，因此只要求解出s就可以實現(xiàn)對地震數(shù)據(jù)的重建。一般而言，地震信號都是在時空域中進行表現(xiàn)，但由于Φ和Ψ都是已知的，因此引入感知矩陣的概念來簡化公式，感知矩陣A就等于A=ΦΨ。式中觀測矩陣Φ∈RM×N，Ψ∈RN×N，其中M?N，即將稀疏信號(K-Sparse)從N維空間通過線性投影到M維空間當中。

y=Φx=ΦΨs=As

(4)

求解上述問題的方法是通過對f(x)=‖Ax-b‖2求導，得到求導之后的函數(shù)：

f′(x)=2AT(Ax-b)

(5)

要獲得x的值只要讓f′(x)=0就可以求得f(x)的最小值，前提是當f(x)為凸函數(shù)，且它的局部最小值就是它的最小值時。如果矩陣A為非奇異矩陣，即矩陣的行列式不為0，有逆矩陣時，該問題有精確解，x=A-1b。若當矩陣A為奇異矩陣時，由于A沒有逆矩陣，導致方程等于0時無法求解，方程沒有精確解。因此這時需要引入懲罰項來協(xié)助進行求解。

由于無法得到精確解，因此引入懲罰項的概念，使得x可以盡量稀疏：

min‖x‖1subject to‖Ax-b‖2≤ε

(6)

minx{F(x)≡‖Ax-b‖2+λ‖x‖1}

(7)

由于式(7)中函數(shù)并不是標準的凸函數(shù)，因此需要對其進行轉換，將式(7)轉換為

=f(α)+λg(α)

(8)

這個L1范數(shù)的正則化問題通常非常難解決，通過使用迭代收縮閾值算法可以解決上述問題，具體流程如圖1所示。

圖1 迭代收縮閾值算法流程圖Fig.1 Iterative shrinking threshold algorithm flowchart

3 迭代收縮閾值網絡

擬解決以下正則化問題：

(9)

式中,ψx表示對x進行某種稀疏變換后的變換系數(shù);λ為正則化參數(shù)，一般來說是預定義的。根據(jù)迭代收縮閾值算法，式(9)可以轉換為式(10)和式(11)。解決式(10)問題的方法可以是在對式(11)進行迭代：

r(k)=x(k-1)-ρΦT(Φx(k-1)-y)

(10)

(11)

在式(10)和式(11)中，k表示ISTA迭代指數(shù);ρ是步長。式(10)是瑣碎的，得到的結果是輸入數(shù)據(jù)，而式(11)實際上是所謂近端映射的特例，即proxλφ(r(k))，此時φ(x)=‖Ψx‖1。形式上，由proxλφ(r)表示的正則化φ的近映射定義為：

(12)

由于地震數(shù)據(jù)往往是非線性的，無法通過簡單的線性變換對其進行處理。因此通過設置卷積層以及應用ReLU(整流線性單元)函數(shù)來實現(xiàn)非線性功能，將兩個卷積層與ReLU函數(shù)連接起來，實現(xiàn)對地震數(shù)據(jù)的非線性處理。

每層卷積層的卷積核不相同，如第一個卷積層(無偏置項)的卷積核的設為符合正態(tài)分布的(1×3×3)的三階張量，共有32個卷積核；第二個卷積核設為符合正態(tài)分布的(32×3×3)的三階張量，同樣具有32個，這一步驟對應著對輸入數(shù)據(jù)的稀疏變換。每一層步驟都對應著ISTA算法的一次迭代，可以看到F(x)=BReLU(Ax)。因此，將F(x)帶入到式(12)中，用F(x)替換Ψx，可以得到以下式子：

(13)

利用ISTA算法很容易就能得到式(14)的解：

(14)

圖2 F(·)函數(shù)結構Fig.2F(·) function structure

(15)

4 數(shù)值實驗

4.1 模擬數(shù)據(jù)

使用規(guī)則采樣以及隨機采樣對模擬數(shù)據(jù)進行采樣。模擬數(shù)據(jù)大小總共為1 101×20 000，也就是將100炮數(shù)據(jù)沿x軸排列，每炮200道接收，總共的道數(shù)為20 000道，時間為1 101 ms。在這將100炮數(shù)據(jù)體重新抽樣排列，每炮抽取100道以及100道的前800 ms傳播時間作為輸入數(shù)據(jù)，并將這100炮按照放炮次序進行排列，得到100×800×100的模擬數(shù)據(jù)體。

本次實驗使用的地質模型為12層地層形成的地質體，震源為50 Hz的Ricker子波，模擬計算地震波在12層地質模型中得到的地震波波場，計算方法為聲波方程。由于深度神經網絡需要大量數(shù)據(jù)集，因此100炮地震數(shù)據(jù)顯然不足以得到一個優(yōu)秀的訓練模型，數(shù)據(jù)量過小容易造成過擬合。因此，在這里將100×800×100的數(shù)據(jù)體中隨機選取95炮，數(shù)據(jù)體大小為95×800×100，使用隨機測量矩陣對這95炮數(shù)據(jù)進行數(shù)據(jù)集的擴充，即對這95炮數(shù)據(jù)體中每一炮都進行10次隨機的采樣，使得95炮數(shù)據(jù)體變?yōu)?50×800×100的數(shù)據(jù)體，將該數(shù)據(jù)體作為ISTA-Net中的訓練集；同理，在抽取95炮數(shù)據(jù)體后，遺留的5炮可以作為測試集，在950×800×100的地震數(shù)據(jù)體訓練完畢后，將剩余5炮的數(shù)據(jù)導入，作為測試集進行測試，并將最終重構的結果進行保存。

將數(shù)據(jù)導入ISTA-Net中，預設總共的對整個ISTA-Net進行200次的迭代運算，在每次ISTA-Net中設置6個階段層，每一批次導入20個訓練塊，總共950炮數(shù)據(jù)，因此設置950個訓練塊。在每一個訓練周期完成后，使用學習率為0.000 1 的Adam算法對ISTA-Net的參數(shù)進行優(yōu)化，訓練200次后終止。圖3表示的是模擬地震記錄第5炮，其色標柱是振幅值。圖3～圖14的色標柱相同。圖4表示的為對其進行隨機采樣50 % 后得到的采樣矩陣。

圖3 模擬地震記錄第4炮數(shù)據(jù)Fig.3 The 4th shot data of simulated seismic records

圖4 隨機采樣50 %的第4炮數(shù)據(jù)Fig.4 Randomly sample 50 % of the 4th shot data

此次模擬數(shù)據(jù)的實驗的結果使用信噪比(SNR，Signal to Noise Ratio)進行衡量，衡量標準SNR的公式為

(16)

在式(15)中，參數(shù)p代表最終訓練得到的重構結果，r為輸入的原始地震數(shù)據(jù)，通過計算SNR可以得到重構效果的優(yōu)劣，SNR越大，則表示重構后的地震數(shù)據(jù)的質量越高，效果越好。

圖5是在使用迭代收縮閾值算法在使用Curvelet變換后的圖像，圖6表示的是Curvelet重構數(shù)據(jù)誤差，Curvelet重構數(shù)據(jù)的SNR=12.691 2。從圖5中可以看出，基于Curvelet變換后的重構數(shù)據(jù)同相軸較為連續(xù)， SNR較高，但重構誤差相對較大，且在多道有連續(xù)缺失的位置，重構的效果并不好。

圖5 Curvelet重構數(shù)據(jù)Fig.5 Curvelet reconstruction data

圖6 Curvelet重構誤差Fig.6 Curvelet reconstruction error

圖7中表示的是利用迭代收縮閾值網絡方法重構的地震數(shù)據(jù)，圖8中是利用該方法的重構誤差，使用迭代收縮閾值網絡方法重構的數(shù)據(jù)的SNR=20.898 1。從圖中可以明顯看出，重構后的地震數(shù)據(jù)地震同相軸連續(xù)性好，無論是直達波還是反射波，其重構效果較好且絕對誤差較小，連續(xù)缺失地震道的位置重構效果依然很好。

圖7 迭代收縮閾值網絡重構數(shù)據(jù)Fig.7 Reconstruction data of iterative shrinkage threshold network

圖8 迭代收縮閾值網絡重構誤差Fig.8 Reconstruction error of iterative shrinkage threshold network

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

本次實驗使用的數(shù)據(jù)來源于神狐海域的一條拖纜測線的數(shù)據(jù)，拖纜數(shù)據(jù)中存在著多種噪聲，在對數(shù)據(jù)進行采樣之前，通過時頻分析發(fā)現(xiàn)數(shù)據(jù)中存在低頻干擾，因此在采樣前對數(shù)據(jù)使用了一個帶通濾波將低頻干擾濾掉。整個數(shù)據(jù)體大小超過5個G，因此效仿模擬數(shù)據(jù)的截取，將拖纜數(shù)據(jù)截取部分用于此次實驗。拖纜數(shù)據(jù)單炮由360道接收，將其中100炮的數(shù)據(jù)提取出來，并且每一炮中提取51道到150道的數(shù)據(jù)，時間方向上選擇1 500 ms到2 300 ms處地震同相軸較為明顯的部分，將其截取并導出，數(shù)據(jù)體大小為100×800×100。其中第4炮的地震數(shù)據(jù)如圖9所示，圖10為隨機采樣50 %后的實際拖纜數(shù)據(jù)。數(shù)據(jù)產生流程與模擬數(shù)據(jù)產生流程相似，隨機抽取100炮中的95炮數(shù)據(jù)，對它進行重采樣10倍后，可以得到950×800×100的數(shù)據(jù)體，這個數(shù)據(jù)體作為訓練集進行訓練，而剩余的5炮作為測試集，測試集數(shù)據(jù)體大小為5×800×100。

圖9 實際拖纜數(shù)據(jù)Fig.9 Actual streamer data

圖10 隨機采樣50 %的拖纜數(shù)據(jù)Fig.10 Randomly sample 50 % of the streamer data

圖11是在迭代收縮閾值算法使用Curvelet變換后的重構數(shù)據(jù)，圖12為使用Curvelet變換后的重構誤差，重構后的數(shù)據(jù)SNR=5.286 0。從圖11中可以看出，基于Curvelet變換后的實際重構數(shù)據(jù)同相軸不連續(xù)，SNR較低，且在多道有連續(xù)缺失的位置以及噪聲較大的位置，重構的效果并不好。

圖11 Curvelet重構數(shù)據(jù)Fig.11 Curvelet reconstruction data

圖12 Curvelet重構誤差Fig.12 Curvelet reconstruction error

圖13中表示的是使用迭代收縮閾值方法重構的實際拖纜數(shù)據(jù)，圖14表示的是利用該方法重構的誤差，重構后的數(shù)據(jù)的SNR=11.073 6。從圖中可以明顯看出，重構后的地震數(shù)據(jù)的地震同相軸連續(xù)性較好，數(shù)據(jù)上半部分直達波的恢復效果較好，絕對誤差相對較小，連續(xù)缺失地震道的位置重構效果較好，重構時去除了一部分噪聲，突出了有效同相軸。

圖13 迭代收縮閾值網絡重構數(shù)據(jù)Fig.13 Reconstruction data of iterative shrinkage threshold network

圖14 迭代收縮閾值網絡重構誤差Fig.14 Reconstruction error of iterative shrinkage threshold network

5 結 論

基于復雜地質條件下的地震數(shù)據(jù)缺失的問題，結合傳統(tǒng)的凸優(yōu)化算法以及近年來在圖像處理效果較好的深度神經網絡方法得到迭代收縮閾值方法，并將其應用在模擬數(shù)據(jù)和實際拖纜數(shù)據(jù)的地震數(shù)據(jù)重構上，通過大量的實驗，主要得到以下幾點結論：

1)通過將迭代收縮閾值算法與深度神經網絡相結合，分析發(fā)現(xiàn)迭代收縮閾值算法的每一步迭代過程可以用深度神經網絡的每一層進行替換，因此可以將迭代收縮閾值算法映射到神經網絡中去，結合深度神經網絡的特點，克服了迭代閾值收縮算子需要預設步長、迭代次數(shù)以及變換域不符的問題，通過將迭代收縮閾值算法映射到深度神經網絡中，各個參數(shù)可以實現(xiàn)端到端的學習，不需要對各個參數(shù)進行預設，其可以通過神經網絡進行自我更新學習，并將稀疏變換矩陣用稀疏變換算子代替，采用非線性可學習的稀疏變換，可以對數(shù)據(jù)進行更好的稀疏操作。

2)通過將迭代收縮閾值網絡重構方法應用于模擬數(shù)據(jù)和實際拖纜數(shù)據(jù)上，發(fā)現(xiàn)該方法在各個方面的重構效果都要好于傳統(tǒng)的迭代收縮閾值算法的表現(xiàn)，尤其是該方法在隨機采樣時對連續(xù)缺失道的恢復效果遠好于其他方法。在實際地震數(shù)據(jù)中迭代收縮閾值網絡重構方法可以用于對含噪數(shù)據(jù)進行一些去噪操作，使得同相軸可以更好的凸顯出來。該方法的重構后的數(shù)據(jù)分辨率較好，SNR值較高，且同相軸連續(xù)。