朱躍飛，曹靜杰

1.河北省戰(zhàn)略性關(guān)鍵礦產(chǎn)資源重點實驗室，河北 石家莊 050031；2.河北地質(zhì)大學 地球科學學院，河北 石家莊 050031

0 引言

地震勘探的野外數(shù)據(jù)采集經(jīng)常會由于成本、地表環(huán)境、壞道等因素造成一些地震道數(shù)據(jù)的缺失，為了避免缺失的地震道對后續(xù)數(shù)據(jù)處理的影響，需要對缺失的地震數(shù)據(jù)進行重建。地震數(shù)據(jù)重建能夠由不滿足采樣定理的數(shù)據(jù)重建出規(guī)則且滿足采樣定理的數(shù)據(jù)，避免成本昂貴的重新野外采集，消除不完整數(shù)據(jù)對疊前時間偏移[1]，三維表面多次波消除[2]和時移地震[3]的影響。地震數(shù)據(jù)重建方法有很多種，主要分為基于信號處理的方法[4-8]，基于波動方程的方法[9-11]和矩陣或張量完備的方法[12-14]，最近基于機器學習和人工智能的方法正在興起[15-16]?；谛盘柼幚淼姆椒ň哂杏嬎阈矢摺?shù)值效果穩(wěn)定和抗噪性強的優(yōu)點，是目前學術(shù)界和工業(yè)界的主流研究方向，該類方法又可以分為基于預測誤差濾波的方法[4]和基于稀疏變換的方法[5-8]。基于矩陣和張量完備的方法是最近幾年新興的研究方向，數(shù)值效果同樣可靠。在地震數(shù)據(jù)采集過程中，由于環(huán)境噪聲、儀器設(shè)備等原因，采集到的地震數(shù)據(jù)經(jīng)常含有各種噪聲，因此含噪聲地震數(shù)據(jù)的重建是地震數(shù)據(jù)重建問題中的一個新的研究分支。該方法避免了含噪聲數(shù)據(jù)的先去噪后重建，或者先重建后去噪，將重建和去噪合二為一減少了工作量并且和先重建后去噪具有相當?shù)男Ч鸞17]，最近對含噪聲數(shù)據(jù)重建的研究逐漸增多。Hennenfent和Herrmann(2006)提出了基于一種非規(guī)則采樣Curvelet變換[18]和迭代閾值法的地震數(shù)據(jù)重建和去噪；Oropeza和Sacchi(2011)基于矩陣完備理論提出了一種同時重建和去噪的加權(quán)重插入方法來消除隨機噪聲[12]，Gao等人(2013)基于傅里葉變換提出一種加權(quán)凸集投影方法來同時重建和去噪[19]，然而上述兩種算法都沒有給出求解的數(shù)學模型。曹靜杰和王本鋒(2015)證明了加權(quán)凸集投影法是求解無約束反演模型的算法，其權(quán)重因子能夠放大閾值從而實現(xiàn)去噪[17]。凸集投影法是地震數(shù)據(jù)重建中使用較多的方法，具有很好的數(shù)值效果和計算速度，但是Abma和Kabir(2006)的文章以Fourier變換為稀疏變換[20]，只是給出了算法的迭代過程，學者們對其求解的數(shù)學模型仍然模糊不清。

文章首先闡述了基于稀疏變換和矩陣完備理論的地震數(shù)據(jù)重建的模型和算法，總結(jié)了一些可以作為稀疏約束的正則算子；然后給出了凸集投影法的推導過程，證明了其是求解等式約束模型的算法，該方法只對不含噪聲數(shù)據(jù)的重建有效。對于反問題來說反演模型的建立是最重要的，當模型選擇正確后正則參數(shù)是否合適決定了反演的數(shù)值效果。在這兩者的基礎(chǔ)上選擇合適的解法才能提高反演的精度和效果。正則參數(shù)的選擇決定基于稀疏反演的含噪聲數(shù)據(jù)的重建問題的成敗，目前關(guān)于含噪聲數(shù)據(jù)重建的文章對正則參數(shù)的選擇幾乎沒有討論，因此文章列出幾種正則參數(shù)的選擇準則和方法。在數(shù)值模擬部分用一種簡單的閾值算法來求解反演模型，二維模擬數(shù)據(jù)、二維真實數(shù)據(jù)和三維真實數(shù)據(jù)的數(shù)值試驗證明了只有選擇合適的正則參數(shù)含隨機噪聲數(shù)據(jù)的重建才能達到好的效果。

1 地震數(shù)據(jù)重建模型與解法

地震數(shù)據(jù)采集過程可以表示為如下的線性過程：

Φd+n=dobs

(1)

其中Φ是采樣矩陣，d是未知的完整波場數(shù)據(jù)，n是可加噪聲，假設(shè)n為隨機白噪聲，dobs為采樣的不完整數(shù)據(jù)。如果完整的波場數(shù)據(jù)d在某個變換域是稀疏或者可壓縮的，即x=Ψd是稀疏或者可壓縮的，則公式(1)可以表示為：

ΦΨ*x+n=dobs

(2)

其中Ψ*為Ψ的共軛轉(zhuǎn)置(一般選擇Ψ是正交變換或者緊框架)。由于采樣數(shù)據(jù)的不完整性，公式(2)是欠定方程組，存在無窮多解?；趚的稀疏性假設(shè)，公式(2)的稀疏解可以通過解如下的優(yōu)化問題找到：

(3)

其中σ是噪聲的能量，0表示非零元素的個數(shù)。0不是連續(xù)可微的，因此問題(3)不存在多項式解法，常用1范數(shù)代替0范數(shù)作為目標函數(shù)，建立如下的基追蹤問題[21]：

(4)

當x求出后，則可以由Ψ*x求出d。

基于矩陣完備理論的重建方法是另一種常見的重建方法，該方法需要假設(shè)同相軸是線性的，假設(shè)完整的地震數(shù)據(jù)可以重排為一個低秩矩陣，根據(jù)重排矩陣的低秩性約束來建立矩陣優(yōu)化模型。地震數(shù)據(jù)可以重排為塊Toeplitz或塊Hankel矩陣[11]，除此之外還有一些其他的重排方式[21-22]。假設(shè)完整的地震數(shù)據(jù)d重排為Md，采樣矩陣Φ重排為MΦ，噪聲n重排為Mn，采樣數(shù)據(jù)dobs重排為Mdobs，公式(1)可以變?yōu)椋?/p>

MΦ·Md+Mn=Mdobs

(5)

根據(jù)矩陣的低秩約束，可以建立如下的矩陣優(yōu)化問題：

minΛ(Md)0s.t.MΦ·Md-MdobsF≤ρ

(6)

其中Λ(Md)0表示Md的奇異值的個數(shù)，Λ(Md)為由Md的奇異值組成的向量，該矩陣優(yōu)化問題同樣不存在多項式解法，一般采用矩陣Md的核范數(shù)來代替Λ(Md)0，建立如下模型：

minMd*s.t.MΦ·Md-MdobsF≤ρ

(7)

其中Md*=Λ(Md)1表示Md的核范數(shù)，·F表示矩陣的F范數(shù)。

當?shù)卣饠?shù)據(jù)噪聲很弱時，可以假設(shè)σ≈0，則問題(4)變?yōu)椋?/p>

minx1s.t.ΦΨ*x=dobs

(8)

如果以時間-空間域的變量為未知數(shù)，公式(8)又可以表示成：

minΨd1s.t.Φd=dobs

(9)

同理當ρ≈0時，問題(7)可以退化為：

minMd*s.t.MΦ·Md=Mdobs

(10)

或者

minF(d)*s.t.Φd=dobs

(11)

其中F(d)=Md表示對數(shù)據(jù)d重新排成Md的過程。

當?shù)卣饠?shù)據(jù)中的噪聲不可忽略時，σ>0或者ρ>0，此時直接解問題(4)和問題(7)相對困難，一般求解它們的無約束形式：

(12)

或

(13)

其中λ為正則參數(shù)，用來平衡數(shù)據(jù)擬合誤差和解的稀疏性(或者低秩性)。反演理論表明對于每個σ和ρ，存在唯一的λ使得約束優(yōu)化問題(4)，(7)和無約束優(yōu)化問題(12)，(13)存在相同的解[28-23]。對于含噪聲數(shù)據(jù)重建，一般求解的是(12)，(13)或者它們的變形：

(14)

和

(15)

除了上述的1范數(shù)作為稀疏約束，還可以采用Cauchy函數(shù)，Huber函數(shù)[24]，1范數(shù)的光滑近似函數(shù)[25]，p-范數(shù)(0

凸集投影法是一種重要的重建算法，該算法計算效率高，數(shù)值效果穩(wěn)定，是工業(yè)界常用的方法。該方法的每次迭代由兩部分組成，第一步是當前迭代解在變換域的閾值運算，第二步通過往凸集上的投影升級迭代解。Abma和Kabir(2006)將該方法引入地震數(shù)據(jù)處理領(lǐng)域時只給出了求解過程[20]，并沒有給出求解的數(shù)學模型。對于重建問題來說該方法求解的是問題(9)，即凸集投影法解的是等式約束問題。

(16)

(17)

令步長α等于μ，則：

(18)

(19)

dk+1=dobs+(I-Φ)Ψ*Tμ(Ψdk)

(20)

根據(jù)上面關(guān)于凸集投影法的推導可知軟閾值運算是專門針對1范數(shù)正則化產(chǎn)生的，變換域的閾值運算可以看作變換域的去噪，閾值類算法實現(xiàn)含噪聲數(shù)據(jù)重建的關(guān)鍵是選擇合適的閾值，閾值和上述的正則參數(shù)是等價的，該參數(shù)和噪聲的能量有關(guān)，當噪聲能量較大時，選擇較大的閾值，當噪聲能量較小時選擇較小的閾值。

經(jīng)典的迭代軟閾值法是解問題(12)的方法[27]，SPGL1方法通過解問題(4)來實現(xiàn)稀疏解的求取[23]，這兩種方法都適合含隨機噪聲數(shù)據(jù)的重建。曹靜杰和王本鋒(2015)證明了加權(quán)凸集投影法是解無約束優(yōu)化問題(14)的方法[17]，根據(jù)加權(quán)凸集投影法的推導過程可知，Oropeza和Sacchi(2011)中加權(quán)重插入算法是解矩陣優(yōu)化問題(15)的算法[12]。該算法通過對重排矩陣Md的奇異值不斷的作閾值運算來重建低秩矩陣，秩的作用如同變換域的系數(shù)個數(shù)，該方法同樣要確定最小特征值的大小和權(quán)重，而權(quán)重可以放大最小特征值，因此該方法和加權(quán)凸集投影法的思想是一致的。Wang等(2011)，Cao和Wang(2014)提出用信賴域算法解以1范數(shù)為信賴域的模型[28-29]，同樣能夠?qū)崿F(xiàn)稀疏解的求取，證明了在合適的模型下，信賴域方法同樣可以作為稀疏解法。除此之外還有很多解上述模型的算法，請參考Cao等[25]。

對于無噪聲地震數(shù)據(jù)重建，同樣可以解上面的無約束優(yōu)化形式，其中的正則參數(shù)只要取非常小的正數(shù)，因此對于不含噪聲數(shù)據(jù)的重建，正則參數(shù)的選取不存在問題。但是當噪聲能量不可忽略時，正則參數(shù)的作用舉足輕重，目前的含噪聲地震重建的文章和算法大都沒有詳細討論如何確定合適的正則參數(shù)。論文對正則參數(shù)的選擇準則進行梳理，列出幾種正則參數(shù)的選擇方法。λ是σ的隱函數(shù)，即兩者不存在顯示關(guān)系。這兩種參數(shù)都可以看作是正則參數(shù)，在反演理論中，存在兩種策略來估計正則參數(shù)：先驗的策略和后驗的策略[30]。先驗策略指的是在計算之前確定正則參數(shù)，但是先驗策略只有理論分析意義，實用性不大。后驗策略指的是在求解過程中估計正則參數(shù)，該類方法的實用性較強，比如Morozov偏差原則，廣義偏差原則，交叉檢驗原則，L曲線等方法[30]，前三種方法都需要知道噪聲的能量。即使正則參數(shù)λ給定，在求解時也往往采用同倫法求解，其思路是求解λ逐漸下降的一系列無約束模型直到給定的λ，用上一個模型的解作為下一個模型的初始解，這樣可以防止得到局部最優(yōu)解。

2 數(shù)值模擬

下面給出算法1求解1范數(shù)正則化的模型(14)來證明當反演模型給定時，正則參數(shù)正確與否決定含噪聲數(shù)據(jù)的重建和去噪的效果。

算法1：

第一步：輸入采樣數(shù)據(jù)dobs，采樣矩陣Φ，稀疏變換Ψ，最大迭代次數(shù)N，噪聲能量的估計ε，最小閾值τmin，令k=0。

第三步：通過軟閾值運算求出變換域的系數(shù)αk+1=Tτk(ΨT(dk+rk))，返回時間-空間域的解dk+1=Ψαk+1。

第五步：輸出最終結(jié)果df=dN。

將Morozov偏差原則和迭代次數(shù)結(jié)合起來作為停機準則。算法1的閾值下降方式采用指數(shù)下降[19]，τ的下降對應于正則參數(shù)λ的下降。

圖1 (a)不含噪聲模擬數(shù)據(jù)；(b)不含噪聲模擬數(shù)據(jù)的F-K譜Fig.1 (a)Noise-free synthetic data;(b)F-K spectrum of Fig.1a

圖2 (a)含隨機噪聲的數(shù)據(jù)，信噪比為6.3963 dB；(b)含隨機噪聲數(shù)據(jù)的頻譜Fig.2 (a)Noisy synthetic data with SNR=6.3963 dB;(b)F-K spectrum of Fig.2a

圖3 (a)隨機采樣的含噪聲數(shù)據(jù)，信噪比2.6498 dB；(b)隨機采樣的含噪聲數(shù)據(jù)的F-K譜Fig.3 (a)Randomly sampled noisy data with SNR=2.6498 dB; (b)F-K spectrum of of Fig.3a

圖4 (a)閾值為0.2時的重建結(jié)果，信噪比13.7644 dB；(b)圖4a的F-K譜Fig.4 (a)Interpolated result with the least threshold value 0.2，SNR=13.7644 dB;(b)F-K spectrum of Fig.4a

圖5 (a)閾值為0.01的重建結(jié)果，信噪比7.1445 dB；(b)圖5a的F-K譜.Fig.5 (a)Interpolated result with the least threshold value 0.01，SNR=7.1445 dB;(b)F-K spectrum of Fig.5a

圖6 (a)圖5(a)與原始數(shù)據(jù)的誤差；(b)圖6(a)的F-K譜.Fig.6 (a)Difference between Fig.5a and the original noiseless data;(b)F-K spectrum of Fig.6a

圖7(a)是一個時間采樣點為276，含有232道的二維地震剖面，時間采樣率為4毫秒，道間距為25米，該數(shù)據(jù)含有一定程度的噪聲，圖7(b)是該數(shù)據(jù)的F-K譜；圖8(a)為圖7(a)隨機采樣50%道的結(jié)果，圖8(b)為其F-K譜；圖9(a)為閾值法的結(jié)果，由圖9(a)可以看出不僅實現(xiàn)了重建而且提高了信噪比。圖9(b)為其F-K譜。最小閾值為0.04，該閾值經(jīng)過測試得到，采用60步迭代。圖10(a)為原始數(shù)據(jù)與重建結(jié)果之差，幾乎不存在有效信號，圖10(b)為圖10(a)的頻譜。

圖7 (a)含噪聲二維剖面；(b)圖7(a)的F-K譜Fig.7 (a)Noisy 2D field data;(b)F-K spectrum of Fig.7a

圖8 (a)隨機采樣50%地震道的二維剖面；(b)圖8(a)的F-K譜Fig.8 (a)Randomly sampled data with 50% traces missing;(b)F-K spectrum of Fig.8a

圖9 (a)迭代閾值法的重建結(jié)果；(b)圖9(a)的F-K譜.Fig.9 (a) Interpolated result with the thresholding method;(b)F-K spectrum of Fig.9a

圖10 (a)重建結(jié)果與原始數(shù)據(jù)的誤差；(b)圖10(a)的F-K譜Fig.10 (a) Difference between the interpolated data and the original data;(b)F-K spectrum of Fig.10a

下面通過一個三維數(shù)據(jù)試驗來進一步驗證文章的結(jié)論。雖然對于三維數(shù)據(jù)也可以用二維變換對時間或頻率切片執(zhí)行，但是這樣做沒有利用數(shù)據(jù)在三維空間的連續(xù)性，因此論文采用三維變換為稀疏變換[32]。圖11(a)為含噪聲的三維數(shù)據(jù)，大小為，時間采樣率為4毫秒，道間距為12.5米。圖11(b)為圖11(a)的采樣，隨機的采樣50%道，圖12(a)為基于本文的重建方法的結(jié)果，閾值經(jīng)過了精心選擇，該方法不僅重建了地震道，而且消除了隨機噪聲，提高了信噪比，圖12(b)為原始數(shù)據(jù)和重建結(jié)果的誤差，幾乎不存在有效信號。

圖11 (a)含噪聲的三維數(shù)據(jù)體；(b)隨機采樣50%的結(jié)果Fig.11 (a)Noise contained 3D data;(b)50% randomly sampled result of Fig.11a

圖12 (a)迭代閾值法重建的結(jié)果；(b)重建結(jié)果與圖11(a)的誤差Fig.12 (a)Interpolation result with the thresholding method;(b)the difference between Fig.12a and Fig.11b

3 結(jié)論與討論

文章對基于稀疏變換和矩陣完備理論的地震數(shù)據(jù)重建問題的數(shù)學模型和算法進行了歸納總結(jié)，給出了凸集投影算法的一種推導過程，證明凸集投影法是求解等式約束反演模型的方法，闡明了其不適合含噪聲數(shù)據(jù)重建的原因。正則參數(shù)的選擇是基于稀疏反演的含噪聲數(shù)據(jù)重建的關(guān)鍵，文章介紹了幾種正則參數(shù)的選擇方法和原則，闡明了只有正則參數(shù)合適算法才能得到合適的解。最后用迭代閾值法求解含噪聲數(shù)據(jù)重建試驗證明正則參數(shù)選擇的重要性。

對于反演問題，正則算子和正則參數(shù)是決定反演成敗的關(guān)鍵，更多的研究應該集中于這兩方面。目前正則參數(shù)一般由經(jīng)驗和嘗試所得，這會耗費一定的時間，研究含噪聲數(shù)據(jù)重建的自適應方法，在沒有過多人為干預的情況下實現(xiàn)含噪聲數(shù)據(jù)的重建和去噪具有重要意義。