含隨機(jī)噪聲地震數(shù)據(jù)重建的稀疏反演模型與算法

朱躍飛，曹靜杰

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

0 引言

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

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

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

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

Φd+n=dobs

(1)

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

ΦΨ*x+n=dobs

(2)

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

(3)

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

(4)

當(dāng)x求出后，則可以由Ψ*x求出d。

基于矩陣完備理論的重建方法是另一種常見的重建方法，該方法需要假設(shè)同相軸是線性的，假設(shè)完整的地震數(shù)據(jù)可以重排為一個(gè)低秩矩陣，根據(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的奇異值的個(gè)數(shù)，Λ(Md)為由Md的奇異值組成的向量，該矩陣優(yōu)化問題同樣不存在多項(xiàng)式解法，一般采用矩陣Md的核范數(shù)來代替Λ(Md)0，建立如下模型：

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

(7)

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

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

minx1s.t.ΦΨ*x=dobs

(8)

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

minΨd1s.t.Φd=dobs

(9)

同理當(dāng)ρ≈0時(shí)，問題(7)可以退化為：

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

(10)

或者

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

(11)

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

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

(12)

或

(13)

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

(14)

和

(15)

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

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

(16)

(17)

令步長(zhǎng)α等于μ，則：

(18)

(19)

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

(20)

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

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

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

2 數(shù)值模擬

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

算法1：

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

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

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

將Morozov偏差原則和迭代次數(shù)結(jié)合起來作為停機(jī)準(zhǔn)則。算法1的閾值下降方式采用指數(shù)下降[19]，τ的下降對(duì)應(yīng)于正則參數(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)含隨機(jī)噪聲的數(shù)據(jù)，信噪比為6.3963 dB；(b)含隨機(jī)噪聲數(shù)據(jù)的頻譜Fig.2 (a)Noisy synthetic data with SNR=6.3963 dB;(b)F-K spectrum of Fig.2a

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

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

圖11 (a)含噪聲的三維數(shù)據(jù)體；(b)隨機(jī)采樣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é)論與討論

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

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