段中鈺 李婷婷 肖 勇 王云雷 鄭桂娟

(①北京信息科技大學(xué)信息與通信工程學(xué)院，北京 100101；②東方地球物理公司國(guó)際勘探事業(yè)部，河北涿州 072751；③東方地球物理公司研究院，河北涿州，072751)

0 引言

地震數(shù)據(jù)重建是地震信號(hào)分析領(lǐng)域的一個(gè)重要研究方向。地震數(shù)據(jù)采集時(shí)由于采集環(huán)境以及采集成本等因素[1]，采集到的數(shù)據(jù)通常呈現(xiàn)不規(guī)則或稀疏分布，不完整的地震數(shù)據(jù)缺失部分地球物理信息，據(jù)其不能準(zhǔn)確地預(yù)測(cè)地下巖層的性質(zhì)。因此，為保證滿足后續(xù)地震資料處理對(duì)數(shù)據(jù)完整性的要求，必須恢復(fù)缺失的地震數(shù)據(jù)，以提高儲(chǔ)層預(yù)測(cè)的準(zhǔn)確性，促進(jìn)油氣的勘探與開發(fā)。顯然，重建缺失的地震數(shù)據(jù)具有非常重要的現(xiàn)實(shí)意義。

傳統(tǒng)的地震數(shù)據(jù)重建方法主要有基于濾波的重建方法[2-4]、基于波動(dòng)方程的重建方法[5-7]以及基于變換函數(shù)的重建方法[8-11]等三種類型。濾波重建法是通過(guò)褶積插值濾波器實(shí)現(xiàn)地震數(shù)據(jù)重建，主要用于重建規(guī)則采樣地震數(shù)據(jù)，且計(jì)算量大。波動(dòng)方程法是通過(guò)Kirchhoff積分算子實(shí)現(xiàn)地震數(shù)據(jù)重建[12]，該方法需地下信息作為先驗(yàn)條件，在地下信息未知或精度低時(shí)對(duì)重建結(jié)果影響很大。變換函數(shù)法是利用最小二乘法將信號(hào)變換到某一新的變換域，再通過(guò)反變換重構(gòu)地震數(shù)據(jù)。常用的變換函數(shù)有傅里葉變換、曲波變換、離散余弦變換等。變換函數(shù)法結(jié)構(gòu)簡(jiǎn)單，計(jì)算復(fù)雜度低，可實(shí)現(xiàn)規(guī)則和非規(guī)則采樣地震數(shù)據(jù)的重建，且不需將地下地層信息作為先驗(yàn)條件，但該方法對(duì)采樣率要求較高，造成采集成本高昂，且要求地震數(shù)據(jù)的采集規(guī)模較大。

壓縮感知理論[13-16](Compressive sensing，CS)的應(yīng)用，克服了傳統(tǒng)地震數(shù)據(jù)重建方法的缺陷，解決了對(duì)采樣率要求較高的問(wèn)題[17]，使稀疏信號(hào)可在低于奈奎斯特采樣頻率的條件下用稀疏優(yōu)化算法重建完整信號(hào)。壓縮感知地震數(shù)據(jù)重建主要分為稀疏表示階段和稀疏重建階段，重建階段相當(dāng)于稀疏優(yōu)化問(wèn)題。白蘭淑等[18]結(jié)合稀疏促進(jìn)的曲波恢復(fù)迭代算法和Bregman迭代閾值算法提出改進(jìn)的聯(lián)合迭代算法，在壓縮感知理論框架下，用Curvelet變換稀疏表示地震數(shù)據(jù)，該改進(jìn)算法能稀疏優(yōu)化重建地震數(shù)據(jù)，具有收斂快、重建信噪比高的優(yōu)點(diǎn)，但其計(jì)算量大、計(jì)算復(fù)雜度高。周亞同等[19]提出壓縮感知下K-SVD字典學(xué)習(xí)自適應(yīng)稀疏表示地震數(shù)據(jù)、正則化正交匹配追蹤(Regularized orthogonal matching pursuit，ROMP)稀疏優(yōu)化重建地震數(shù)據(jù)的方法，能加快復(fù)雜實(shí)際地震數(shù)據(jù)重建速度。但該稀疏優(yōu)化算法得到的不是全局最優(yōu)解，故重建精度較低。

交替乘子方向算法(Alternating direction method of multipliers，ADMM)[20-22]具有簡(jiǎn)化復(fù)雜性、高重建速度和良好重建質(zhì)量的優(yōu)點(diǎn)，在地震數(shù)據(jù)領(lǐng)域引起了廣泛關(guān)注。Sternfels等[23]將ADMM算法與壓縮感知理論結(jié)合實(shí)現(xiàn)含噪地震數(shù)據(jù)的降噪。Zhang等[24]將ADMM算法應(yīng)用于缺失地震數(shù)據(jù)的恢復(fù)，實(shí)現(xiàn)壓縮感知地震數(shù)據(jù)的精確重構(gòu)。李慧等[25]在此基礎(chǔ)上提出壓縮感知框架下K-SVD字典學(xué)習(xí)稀疏表示地震數(shù)據(jù)，利用ADMM算法優(yōu)化重建地震數(shù)據(jù)，達(dá)到提高地震數(shù)據(jù)重建精度的目的。

盡管ADMM算法在地震數(shù)據(jù)重建中效果顯著，但由于該算法中平衡因子無(wú)法根據(jù)地震數(shù)據(jù)自適應(yīng)地進(jìn)行調(diào)整，易損害重建精度。另外，在基于ADMM算法的壓縮感知地震數(shù)據(jù)重建模型中，由于輸入訓(xùn)練樣本較少，算法參數(shù)較多以及模型復(fù)雜度較高等問(wèn)題，容易出現(xiàn)過(guò)擬合現(xiàn)象。針對(duì)上述問(wèn)題，本文對(duì)ADMM算法進(jìn)行改進(jìn)，提出平方正則交替乘子方向算法(Square regular-alternating direction method of multipliers，SR-ADMM)。一方面在ADMM算法迭代中加入平方正則項(xiàng)；另一方面自適應(yīng)選取平衡因子。最終將SR-ADMM算法與壓縮感知理論相結(jié)合，實(shí)現(xiàn)基于壓縮感知的SR-ADMM算法地震數(shù)據(jù)重建。通過(guò)模擬地震數(shù)據(jù)、實(shí)際單炮數(shù)據(jù)以及大慶油田實(shí)際地震數(shù)據(jù)對(duì)本文所提算法進(jìn)行了驗(yàn)證。

1 壓縮感知基本理論

壓縮感知理論具有三個(gè)重要的前提條件：地震數(shù)據(jù)的稀疏性、測(cè)量矩陣的不相干性以及合適的優(yōu)化算法。壓縮感知的數(shù)據(jù)重建模型為

y=Φx

(1)

式中：x∈RN，為原始信號(hào)；Φ∈RM×N，為測(cè)量矩陣；y∈RM，為測(cè)量矩陣下原始信號(hào)的觀測(cè)值。該式表示原始信號(hào)從N維信號(hào)降維到M維獲得觀測(cè)值，信號(hào)降維過(guò)程可用圖1所示圖像表示。

圖1 原始信號(hào)降維

由于式(1)為欠定方程，有無(wú)窮多的解，很難直接通過(guò)觀測(cè)值y恢復(fù)出原始信號(hào)x。因此可用一個(gè)與測(cè)量矩陣不相關(guān)的變換矩陣對(duì)原始信號(hào)做稀疏表示，得到稀疏矩陣。其數(shù)學(xué)表達(dá)式為

x=DΘ

(2)

式中：D∈RN×N，為稀疏變換矩陣；Θ為稀疏系數(shù)，是長(zhǎng)度為N的一維向量。式(2)的原始信號(hào)稀疏表示過(guò)程可用圖2所示圖像表示。

圖2 信號(hào)的稀疏表示

因此觀測(cè)值y可以表示為

y=Φx=ΦDΘ=AΘ

(3)

式中A=ΦD，為M×N階感知矩陣。A滿足約束等距性質(zhì)，即可通過(guò)求解L0范數(shù)優(yōu)化問(wèn)題很好地從觀測(cè)值y中恢復(fù)出原始信號(hào)x

(4)

式中稀疏系數(shù)Θ的L0范數(shù)表示該向量中所有非零元素的個(gè)數(shù)。但L0范數(shù)的最優(yōu)化問(wèn)題是一個(gè)非確定性多項(xiàng)式問(wèn)題，在稀疏變換矩陣D與測(cè)量矩陣Φ不相干時(shí)，求解Θ的L0范數(shù)問(wèn)題等價(jià)于求解Θ的L1范數(shù)，故可將式(4)轉(zhuǎn)化為求解L1范數(shù)優(yōu)化問(wèn)題

(5)

(6)

2 SR-ADMM算法

適用的重構(gòu)算法是壓縮感知精確重構(gòu)的重要前提。ADMM算法的優(yōu)點(diǎn)在于將對(duì)偶上升法的可分解性與乘子法的收斂性質(zhì)結(jié)合起來(lái)，通過(guò)對(duì)偶上升法使線性稀疏項(xiàng)不斷接近最優(yōu)解，保證在懲罰參數(shù)較小的情況下滿足精度要求；克服乘子法中變量二次項(xiàng)無(wú)法拆分的問(wèn)題，將原始問(wèn)題拆分求解，可實(shí)現(xiàn)并行運(yùn)算且提高運(yùn)算效率。使用ADMM算法求解壓縮感知地震數(shù)據(jù)重建問(wèn)題的數(shù)學(xué)模型為

(7)

式中λ表示平衡因子，用于控制前、后兩項(xiàng)之間的權(quán)重。

引入輔助變量Z

(8)

該式的增廣拉格朗日函數(shù)為

(9)

根據(jù)上式可知，Θ的更新就是在固定參數(shù)Z和λ的情況下，求解增廣拉格朗日函數(shù)的最小值問(wèn)題，即

L′Θ=AT(Aθk+1-y)+λk+ρ(θk+1-zk)

(10)

求解可知ADMM算法第k+1次的迭代表達(dá)式為

(11)

式中：u=λ/ρ，其中ρ表示懲罰參數(shù)；Sλ/ρ為軟閾值函數(shù)，定義為

(12)

以K-SVD字典學(xué)習(xí)對(duì)地震數(shù)據(jù)進(jìn)行稀疏表示，并用ADMM算法解決稀疏優(yōu)化問(wèn)題，構(gòu)成基于ADMM算法和K-SVD字典學(xué)習(xí)的壓縮感知地震數(shù)據(jù)重建模型。在該模型中，由于存在輸入訓(xùn)練樣本較少、ADMM算法迭代參數(shù)較多及模型復(fù)雜度較高等問(wèn)題，容易出現(xiàn)過(guò)擬合現(xiàn)象，因此增加平方鄰近項(xiàng)正則化式中的子問(wèn)題，防止過(guò)擬合。正則化交替乘子方向算法公式如下

(13)

另外，由于ADMM算法中平衡因子的取值只能根據(jù)經(jīng)驗(yàn)選取，若平衡因子λ較大，則說(shuō)明稀疏項(xiàng)權(quán)重更大，所求值更稀疏而不會(huì)過(guò)多考慮約束；反之，重建數(shù)據(jù)更接近于原始信號(hào)。為了使λ更精確，從而獲得更高重建精度，本文將λ的選取改進(jìn)為可自適應(yīng)方式。對(duì)平衡因子的改進(jìn)方案如下

λ1，a=λ×αl，λ2，a=λ×αl-1,…

(14)

λ1，b=λ/αl，λ2，b=λ/αl-1,…

(15)

首先根據(jù)經(jīng)驗(yàn)設(shè)定初始平衡因子λ。兩式中α的作用是選取λ的鄰近值，可選擇略大于1或略小于1的數(shù)，如α=0.95。式(14)和式(15)分別表示對(duì)初始平衡因子做乘法和除法，獲得多個(gè)近似值。接著分別用所獲近似值以式(13)進(jìn)行迭代實(shí)現(xiàn)重建，選取使地震數(shù)據(jù)重建誤差最小的平衡因子為最優(yōu)平衡因子，通過(guò)一輪迭代就可得到平衡因子最優(yōu)值；之后，繼續(xù)完成后續(xù)迭代過(guò)程。兩式中常數(shù)l表征近似解的個(gè)數(shù)，若l取值越大，就能獲得數(shù)量更多的平衡因子近似值，理論上重建效果更好。但l過(guò)大，會(huì)造成時(shí)間空間負(fù)荷，因此需根據(jù)地震數(shù)據(jù)實(shí)情選取。

3 基于壓縮感知的SR-ADMM地震 數(shù)據(jù)重建

本文基于壓縮感知的SR-ADMM算法地震數(shù)據(jù)重建的流程如圖3所示。

圖3 重建流程圖

在壓縮感知框架下，以K-SVD字典學(xué)習(xí)對(duì)地震數(shù)據(jù)做稀疏表示，用SR-ADMM算法稀疏優(yōu)化實(shí)現(xiàn)缺失地震數(shù)據(jù)重建。缺失數(shù)據(jù)重建具體步驟為：

(1)輸入訓(xùn)練樣本，對(duì)重建模型進(jìn)行訓(xùn)練，得到初始字典D0、初始稀疏系數(shù)Θ0；

(2)輸入原始缺失地震數(shù)據(jù)x，用初始字典D0對(duì)該地震數(shù)據(jù)進(jìn)行稀疏表示；

(4)以K-SVD字典學(xué)習(xí)算法更新稀疏變換字典D，用更新后的字典D對(duì)步驟(3)初步重建地震數(shù)據(jù)做稀疏表示；

4 地震數(shù)據(jù)仿真驗(yàn)證

將信噪比作為地震數(shù)據(jù)重建性能的評(píng)判指標(biāo)

(16)

4.1 模擬地震數(shù)據(jù)重建仿真驗(yàn)證

驗(yàn)證模型的模擬地震數(shù)據(jù)(圖4)共76道，522個(gè)采樣點(diǎn)，采樣間隔為1ms。圖5為該炮點(diǎn)集隨機(jī)采樣數(shù)據(jù)。ADMM算法中平衡因子根據(jù)多次實(shí)驗(yàn)后得出，選取λ=0.15；SR-ADMM算法平衡因子則根據(jù)模擬地震數(shù)據(jù)自適應(yīng)地選取。

圖4 原始模擬地震數(shù)據(jù)

OMP算法是壓縮感知框架下常用于解決優(yōu)化問(wèn)題的一種稀疏優(yōu)化算法。分別采用OMP(圖6a)、ADMM(圖6b)和SR-ADMM(圖6c)三種算法對(duì)圖5所示隨機(jī)采樣模擬地震數(shù)據(jù)進(jìn)行重建。

圖5 隨機(jī)采樣模擬地震數(shù)據(jù)

上述三種算法的重建信噪比和均方誤差如表1所示。將這三種算法重建結(jié)果(圖6，深藍(lán)色波形)與原始地震數(shù)據(jù)(黑色波形)進(jìn)行對(duì)比，并將缺失較嚴(yán)重的第19～第33道數(shù)據(jù)進(jìn)行局部放大?？梢妼?duì)于該隨機(jī)缺失模擬地震數(shù)據(jù)，基于壓縮感知理論的三種算法都能實(shí)現(xiàn)數(shù)據(jù)重建，OMP算法(圖6a)重建效果一般，從隨機(jī)缺失的第21道～第26道數(shù)據(jù)，易見缺失部分重建振幅與原始振幅存在誤差，且該算法對(duì)于未缺失部分(如第45道)數(shù)據(jù)產(chǎn)生影響，重建的均方誤差相對(duì)其他兩種方法較大，重建信噪比相對(duì)其他兩種方法最低；ADMM算法(圖6b)重建效果較好，不會(huì)對(duì)未缺失部分?jǐn)?shù)據(jù)產(chǎn)生影響，但缺失數(shù)據(jù)的振幅尚未得到完全恢復(fù)；SR-ADMM算法(圖6c)重建的效果最好，對(duì)數(shù)據(jù)細(xì)節(jié)的恢復(fù)更完備，同相軸清晰可見，重建數(shù)據(jù)振幅與原始地震數(shù)據(jù)振幅一致，且該算法重建信噪比最高，重建均方誤差最小。

表1 模擬數(shù)據(jù)重建信噪比和均方誤差

圖6 三種算法重建結(jié)果與原始數(shù)據(jù)對(duì)比(a)OMP；(b)ADMM；(c)SR-ADMM

由于SR-ADMM算法中平衡因子是根據(jù)地震數(shù)據(jù)自適應(yīng)地選取，所選平衡因子使數(shù)據(jù)重建誤差最小，從而達(dá)到提高重建精度的目的。而且，基于該算法的壓縮感知數(shù)據(jù)重建模型防止了數(shù)據(jù)的過(guò)擬合現(xiàn)象，一定程度上也提高了數(shù)據(jù)重建精度。因此，針對(duì)隨機(jī)缺失模擬地震數(shù)據(jù)，采用SR-ADMM算法重建的信噪比高于OMP和ADMM算法，即SR-ADMM算法的重建精度高于經(jīng)典ADMM算法。

4.2 實(shí)際地震數(shù)據(jù)仿真驗(yàn)證

選取的M油田實(shí)際單炮地震數(shù)據(jù)(圖7)共有170道，采樣點(diǎn)數(shù)為1300，時(shí)間采樣率為1ms；其缺失率為40%的單炮地震數(shù)據(jù)如圖8a所示。分別使用OMP、ADMM及SR-ADMM三種算法對(duì)缺失數(shù)據(jù)(圖8a)進(jìn)行重建(圖8b～圖8d)。ADMM算法中平衡因子據(jù)多次試驗(yàn)后設(shè)定為λ=0.17，SR-ADMM算法的平衡因子據(jù)實(shí)際地震數(shù)據(jù)自適應(yīng)地選取。三種算法對(duì)應(yīng)的重建信噪比和均方誤差如表2所示。

圖7 實(shí)際地震數(shù)據(jù)單炮記錄

表2 實(shí)際數(shù)據(jù)重建信噪比和均方誤差

結(jié)合表2和圖8b～圖8d可看出：OMP算法(圖8b)重建效果欠理想，同相軸連續(xù)性較差，重建信噪比較低，均方誤差較大；ADMM算法(圖8c)重建效果比OMP算法稍好，大致實(shí)現(xiàn)了缺失數(shù)據(jù)的重建，但重建后的數(shù)據(jù)附近出現(xiàn)抖動(dòng)，重建信噪比較高，均方誤差中等；SR-ADMM算法(圖8d)重建效果很好，同相軸更光滑，重建振幅連續(xù)性良好，在三種算法中重建信噪比最高，均方誤差最小。

圖8 實(shí)際缺失地震數(shù)據(jù)及三種算法重建后單炮記錄(a)實(shí)際缺失數(shù)據(jù)；(b)、(c)、(d)對(duì)應(yīng)OMP、ADMM、SR-ADMM算法

為了更直觀地對(duì)比三種優(yōu)化算法的重建效果，將隨機(jī)采樣原始數(shù)據(jù)中缺失較多的第59～第109道進(jìn)行局部放大(圖9a)，其中第59～第67道、第72～第74道、第76道、第78～第87道、第90～第95道、第100道及第104道為缺失數(shù)據(jù)。分別采用OMP(圖9b)、ADMM(圖9c)和SR-ADMM(圖9d)三種算法得到相應(yīng)重建結(jié)果的局部放大圖。

圖9 原始數(shù)據(jù)及OMP、ADMM、SR-ADMM三種算法重建結(jié)果局部放大(a)原始地震數(shù)據(jù)局部放大；(b)、(c)、(d)對(duì)應(yīng)OMP、ADMM、SR-ADMM三種算法重建結(jié)果的局部放大

對(duì)比圖9a與圖9b可看出部分缺失數(shù)據(jù)實(shí)現(xiàn)了重建(圖9b藍(lán)色箭頭所指)，部分缺失數(shù)據(jù)未得到恢復(fù)(紅色箭頭所指)，且重建數(shù)據(jù)振幅與原始數(shù)據(jù)有差距；對(duì)比圖9a與圖9c可見，缺失數(shù)據(jù)雖得到了恢復(fù)(圖9c藍(lán)色箭頭)，但重建后的振幅與原始數(shù)據(jù)存在誤差；對(duì)比圖9a與圖9d發(fā)現(xiàn)，所有缺失數(shù)據(jù)都得到有效恢復(fù)(圖9d藍(lán)色箭頭)，且重建數(shù)據(jù)振幅與原始數(shù)據(jù)一致。

因此，采用SR-ADMM算法可以很好地重建實(shí)際地震數(shù)據(jù)復(fù)雜的細(xì)節(jié)，重建性能較好。

4.3 三維實(shí)際地震數(shù)據(jù)應(yīng)用

上述模擬和實(shí)際地震單炮記錄驗(yàn)證了本文方法的有效性，進(jìn)一步將本文方法應(yīng)用于中國(guó)石油大慶油田實(shí)際三維數(shù)據(jù)。從該三維地震數(shù)據(jù)的一條測(cè)線中選取200道(圖10)，時(shí)間采樣率為1ms，采樣點(diǎn)為500個(gè)。圖11a為隨機(jī)缺失30%的該三維數(shù)據(jù)實(shí)例。ADMM算法中平衡因子據(jù)多次實(shí)驗(yàn)后設(shè)定為λ=0.20，SR-ADMM算法的平衡因子基于實(shí)際地震數(shù)據(jù)自適應(yīng)地選取。

圖10 三維地震數(shù)據(jù)實(shí)例

同樣分別采用OMP、ADMM及SR-ADMM三種算法對(duì)隨機(jī)采樣實(shí)際三維地震數(shù)據(jù)進(jìn)行重建，得到相應(yīng)的重建結(jié)果(圖11b～圖11d)?？梢奜MP算法的數(shù)據(jù)重建(圖11b)效果一般，未能很好地實(shí)現(xiàn)缺失數(shù)據(jù)的恢復(fù)重建，同相軸欠連續(xù)；ADMM算法重建數(shù)據(jù)(圖11c)效果相對(duì)較好，同相軸連續(xù)性提高，但存在同相軸抖動(dòng)現(xiàn)象；SR-ADMM算法很好地實(shí)現(xiàn)了三維數(shù)據(jù)重建(圖11d)，缺失位置的數(shù)據(jù)與相鄰位置同相軸的走向相符合，且很連續(xù)。

圖11 實(shí)際缺失三維地震數(shù)據(jù)及三種算法重建結(jié)果(a)實(shí)際缺失數(shù)據(jù)；(b)、(c)、(d)對(duì)應(yīng)OMP、ADMM、SR-ADMM算法

為便于比較三種算法重建效果，將三種重建結(jié)果的第50～第60道與原始數(shù)據(jù)進(jìn)行放大對(duì)比(圖12)。三種算法的重建信噪比及均方誤差如表3所示。

綜合分析表3和圖12所示局部放大圖，可以看出OMP算法恢復(fù)地震數(shù)據(jù)(圖12a)的振幅與原始地震數(shù)據(jù)(黑色)相差較大，且重建信噪比較低，均方誤差較大；ADMM算法重建(圖12b)的振幅與原始地震數(shù)據(jù)更契合，其重建信噪比高于OMP算法，均方誤差較?。欢鳶R-ADMM算法重建后(圖12c)振幅與原始地震數(shù)據(jù)幾乎一致，重建信噪比最高，均方誤差最小。顯然，SR-ADMM算法在實(shí)際三維地震數(shù)據(jù)重建中也有效，且可提高實(shí)際三維地震數(shù)據(jù)的重建精度。

表3 三維實(shí)際數(shù)據(jù)重建信噪比和均方誤差

圖12 三種算法重建結(jié)果(藍(lán)色)與原始數(shù)據(jù)(黑色)的波形對(duì)比(a)OMP；(b)ADMM；(c)SR-ADMM

5 結(jié)束語(yǔ)

本文提出基于壓縮感知的SR-ADMM算法地震數(shù)據(jù)重建方法。該算法在原始ADMM算法基礎(chǔ)上加入正則項(xiàng)正則化迭代中子問(wèn)題，并且自適應(yīng)地選取平衡因子，利用該算法解決最優(yōu)化問(wèn)題實(shí)現(xiàn)缺失地震數(shù)據(jù)的重建。仿真驗(yàn)證表明利用本文算法能有效地實(shí)現(xiàn)缺失地震數(shù)據(jù)，且其重建精度高于常規(guī)壓縮感知地震數(shù)據(jù)重建優(yōu)化算法。

需要指出的是，數(shù)據(jù)的稀疏性是壓縮感知數(shù)據(jù)重建的前提，對(duì)重建效果有很大影響，數(shù)據(jù)越稀疏，重建效果越好。因此，如何更好的對(duì)地震數(shù)據(jù)進(jìn)行稀疏表示是下一步的研究方向。