曹靜杰， 王本鋒

1 石家莊經(jīng)濟學院， 石家莊 050031 2 中國石油大學(北京)，油氣資源與勘探測國家重點實驗室， 北京 102249

?

基于一種改進凸集投影方法的地震數(shù)據(jù)同時插值和去噪

曹靜杰1， 王本鋒2

1 石家莊經(jīng)濟學院， 石家莊 050031 2 中國石油大學(北京)，油氣資源與勘探測國家重點實驗室， 北京 102249

基于稀疏反演的地震插值方法是一種重要的插值方法，然而大多數(shù)這類方法只針對無噪聲數(shù)據(jù)或者高信噪比數(shù)據(jù)插值.實際上，地震數(shù)據(jù)含有各種噪聲，使得插值問題變得更加困難.凸集投影方法是一種高效的插值算法，但是對于含噪聲數(shù)據(jù)的插值效果不理想，針對含噪聲數(shù)據(jù)提出的加權(quán)凸集投影方法能夠?qū)崿F(xiàn)同時插值和去噪，但是除了最小閾值需要認真選取外，增加一個權(quán)重因子來實現(xiàn)去噪功能.本文由迭代閾值算法推導出加權(quán)凸集投影方法，證明其是解無約束優(yōu)化問題的一種方法，加權(quán)因子可以看作擬合誤差項的系數(shù).本文還提出了一種改進的凸集投影方法，與原始凸集投影方法相比該方法不需要增加任何計算量，只要通過閾值的選擇來進行插值和去噪.數(shù)值模擬證明了該算法的計算效率，并且對含噪聲數(shù)據(jù)能夠?qū)崿F(xiàn)較好的插值效果；先插值后去噪的結(jié)果證明了同時去噪和插值算法的可靠性和穩(wěn)定性.

稀疏反演； 去噪； 插值； 反演

1 引言

由于采集成本、地形、壞道等因素的限制，野外采集的地震數(shù)據(jù)經(jīng)常不滿足香農(nóng)采樣定理(陸基孟， 1993). 不完整的數(shù)據(jù)會產(chǎn)生采集腳印，嚴重地影響疊前時間偏移(Symes，2007)、三維表面多次波消除(Berkhout and Verschuur，1997)和全波形反演(Bunks et al.，1995)的效果，因此需要采用插值方法來得到完備的波場信息. 地震數(shù)據(jù)插值不僅能夠提供可靠的波場數(shù)據(jù)，還可以指導采集設計，節(jié)約采集成本(Yilmaz，1987).目前，地震插值方法主要包括預測濾波方法(Spitz， 1991)，基于稀疏變換的方法(Liu and Sacchi，2004；Zwartjes，2005；Herrmann and Hennenfent，2008；Xu et al.，2010)，基于波動方程的方法(Ronen，1987；Stolt，2002)和矩陣降秩方法(Oropeza and Sacchi，2011；Kreimer and Sacchi，2012；Ma，2013). 預測濾波方法專門針對規(guī)則不完備采樣數(shù)據(jù)的插值(Spitz，1991)；基于波動方程的方法需要地下的參數(shù)信息并且計算量較大. 隨著壓縮感知理論引入到圖像和三維數(shù)據(jù)的填充，矩陣/張量完備方法也被用于地震數(shù)據(jù)插值(Oropeza and Sacchi，2011；Kreimer and Sacchi，2012；Ma， 2013)，該方法將插值問題看作圖像填充問題，基于完整的矩陣或張量的秩最小化來實現(xiàn)地震數(shù)據(jù)插值. 基于稀疏變換的方法利用各種變換，比如 Fourier 變換(Liu and Sacchi，2004；Zwartjes，2005；Xu et al.， 2010)，Radon 變換(Trad et al.，2002)，seislet 變換(Fomel and Liu，2010)，Dreamlet變換(耿瑜等，2012)，curvelet 變換(Cao et al.，2011；Wang et al., 2011；曹靜杰等，2012；Mansour et al.，2013)等作為稀疏變換來壓縮地震數(shù)據(jù)，由于其高效性和穩(wěn)定性，該方法得到了很大的發(fā)展.

野外地震數(shù)據(jù)經(jīng)常包含不同類型的噪聲，這些噪聲可以分為兩類：隨機噪聲和相干噪聲.目前的插值方法很少直接研究含噪聲數(shù)據(jù)的插值.對于含隨機噪聲的采樣不足數(shù)據(jù)，一般是先插值后去噪. Hennenfent 和 Herrmann (2006)提出了一種基于非規(guī)則采樣curvelet 變換和閾值法的地震數(shù)據(jù)插值和去噪；Oropeza 和Sacchi (2011)提出了一種同時插值和去噪的重加權(quán)替換方法，通過加權(quán)來消除隨機噪聲，Gao 等(2013)提出一種加權(quán)凸集投影(POCS)方法來同時去噪和插值.

由于地震數(shù)據(jù)插值需要處理大量數(shù)據(jù)， 因此屬于大規(guī)模計算問題， 需要快速穩(wěn)定的插值算法. 目前凸集投影(POCS)算法是工業(yè)界一種重要的插值方法(Abma and Kabir，2006)，該方法計算效率很高，容易實現(xiàn)，并且具有很好的插值效果.但是對含噪聲數(shù)據(jù)的插值，由于采樣的噪聲被帶入到插值結(jié)果，因此張華和陳小宏(2013)證明該方法的結(jié)果不太理想. 加權(quán)POCS算法作為其改進能夠?qū)崿F(xiàn)去噪功能，但是除了需要調(diào)節(jié)最小的閾值外，增加的權(quán)重因子同樣也需要認真選取.

本文針對含隨機噪聲數(shù)據(jù)的插值問題進行研究.我們從一般的0范數(shù)正則化和優(yōu)化問題出發(fā)，由迭代閾值方法推導出了加權(quán)POCS方法，證明了加權(quán)POCS方法是解無約束優(yōu)化問題的一種方法，權(quán)重因子可以看作是擬合誤差項的系數(shù)(Oropeza and Sacchi，2011；Gao et al.，2013).同時，我們還提出了一種改進的POCS方法來實現(xiàn)同時插值和去噪，該方法不需要增加任何的計算量，能夠適合含噪聲數(shù)據(jù)的插值.數(shù)值模擬表明，改進的POCS方法在保持POCS方法計算效率高的同時，能夠很好地重建無噪聲的數(shù)據(jù)和含噪聲的數(shù)據(jù).通過和先插值后去噪的結(jié)果比較，證明了同時插值和去噪的穩(wěn)定性.

2 地震數(shù)據(jù)的同時插值和去噪

假設地震數(shù)據(jù)的采集過程可以表示為以下的數(shù)學模型：

Φd+n=dobs，

(1)

其中Φ是采樣算子矩陣，d是完整的波場數(shù)據(jù)，n為可加性噪聲，此處假設為隨機白噪聲，dobs為采集的地震數(shù)據(jù). 如果完整的波場數(shù)據(jù)d在某個變換域是稀疏的，即x=ΨTd是某個變換域的系數(shù)， 是稀疏的或可壓縮的，其中ΨT是一個正交變換或者緊框架，則公式(1)可以表示為

ΦΨx+n=dobs.

(2)

定義A=ΦΨ，則公式 (2)的解可由以下優(yōu)化問題得到

(3)

其中λ是用來平衡擬合誤差項和x的稀疏性約束項的正則參數(shù)，R(x)是關(guān)于x的算子，是對解的稀疏性的描述. 地震數(shù)據(jù)插值屬于大規(guī)模計算問題，因此需要快速穩(wěn)定的算法. 在眾多稀疏反演方法中，迭代閾值類算法是一種有效的解法. 如果R(x)=‖x‖1，則可以用迭代軟閾值法解問題(3)(Daubechiesetal.，2004)； 如果R(x)=‖x‖0，則迭代硬閾值法是解問題(3)的一個很好的策略(BlumensathandDavies，2008).理論上，存在唯一的λ使得問題 (3) 能得到最優(yōu)解(王彥飛，2007). 但是，迭代閾值類方法屬于局部優(yōu)化方法，為了防止得到局部解，在初始迭代時參數(shù)λ應當選擇較大的值，然后逐漸減小到合適的值.

現(xiàn)有的基于問題(3)的大多數(shù)插值方法只是研究無噪聲數(shù)據(jù)或高信噪比數(shù)據(jù)的插值，實際上問題(3) 同樣適合于含有隨機噪聲數(shù)據(jù)的插值，如果參數(shù)λ選擇合適就能夠?qū)崿F(xiàn)同時插值和去噪.參數(shù)λ和迭代閾值類算法中的閾值關(guān)系緊密，直接影響閾值算法的結(jié)果.迭代閾值類算法一般都是在變換域中進行閾值運算，本質(zhì)上，變換域中的閾值處理是一種去噪過程，因此基于閾值類算法的插值可以看作變換域的去噪. 參數(shù)λ與噪聲的能量有關(guān)，較大的噪聲對應于較大的λ. 如果λ選擇合適，則能夠達到插值效果和去噪的最佳平衡. 然而不同數(shù)據(jù)的噪聲能量是不同的，因此一般很難選擇精確的λ.

Oropeza 和 Sacchi (2011)提出一種基于多道奇異譜分析的重加權(quán)替換方法，其中地震數(shù)據(jù)被重新排列為塊Hankel/Toeplitz矩陣，通過降低Hankel/Toeplitz矩陣的秩來進行插值，秩的作用和λ的作用類似. 該算法同樣通過逐漸降低秩來實現(xiàn)插值，因此多道奇異譜分析方法的思想和迭代閾值法的思想是一致的.

3 加權(quán)POCS方法和改進的POCS方法

凸集投影法是地震插值領(lǐng)域最常用的方法之一，該方法計算效率高，結(jié)果穩(wěn)定，非常適合大規(guī)模計算，但該方法對含噪聲數(shù)據(jù)的插值效果不理想，為此Gao等(2013) 提出一種加權(quán)POCS方法來插值含噪聲數(shù)據(jù).下面我們對加權(quán)POCS方法進行推導，證明其是解無約束優(yōu)化問題的一種方法，權(quán)重因子可以看作是擬合誤差項的系數(shù).Yang等(2012)證明了對于插值問題，迭代硬閾值法可以推導出POCS方法.我們借鑒其中的思想推導出加權(quán)POCS方法.POCS方法可以用硬閾值迭代，也可以用軟閾值迭代，取決于選用的正則算子是0范數(shù)還是1范數(shù). 本文采用硬閾值進行推導，軟閾值的推導過程和硬閾值的推導過程是一樣的.迭代硬閾值法的迭代公式為

(4)其中N是迭代次數(shù)，xk為第k個迭代解，

是硬閾值運算，τ是閾值.τk隨著迭代次數(shù)的增大而減小.τk的下降方式對插值結(jié)果有一定的影響，但不是關(guān)鍵的，關(guān)鍵的參數(shù)是最小閾值(正則參數(shù))的大小，對于含噪聲數(shù)據(jù)，應該仔細地選擇最小閾值，對于無噪聲數(shù)據(jù)，最小閾值可以很小.τk的下降方式可以采用線性下降或者指數(shù)下降，本文采用指數(shù)下降方式. 根據(jù)公式(4)，迭代硬閾值方法逐漸更新變換域的系數(shù)xk，當?shù)玫阶罱K的迭代解xN時，只要通過反變換就能得到時間-空間域的插值結(jié)果.POCS方法的迭代公式為

(5)

其中I是元素全部為1的矩陣，該方法通過更新時間-空間域的解uk來得到最優(yōu)解. 如果POCS法的初始解為0，則POCS方法第一次迭代的解為dobs，因此POCS方法的第一次迭代就得到一個好的初始解. 雖然POCS方法插值的信噪比在前幾次迭代提高的很快，但是多次迭代后其最終解和迭代硬閾值法的解是一樣的(Yangetal.，2012).

加權(quán)POCS方法利用加權(quán)的采樣數(shù)據(jù)和加權(quán)的迭代解來構(gòu)造新的迭代解(Gaoetal.，2013)，其中的權(quán)重因子是一個關(guān)鍵的參數(shù). 下面給出其推導過程，定義

(6)

=Hτk(xk-αAT(Axk-dobs))，

(7)

其中τk為第k個閾值參數(shù). 對公式(7)兩邊乘以Ψ，并且定義dk=Ψxk，可以得到

dk+1=Ψxk+1

=ΨHτk(xk+α(ΦΨ)T(dobs-Φdk))

=ΨHτk(ΨTdk+α(ΦΨ)T(dobs-Φdk))

=ΨHτk(ΨTdk+αΨT(ΦTdobs-ΦTΦdk))

=ΨHτk(ΨT(dk+αdobs-αΦdk))

=ΨHτk(ΨT(αdobs+(I-αΦ)dk)).

(8)

關(guān)于(8)的推導可以參考文獻Yang等(2012). 定義uk=αdobs+(I-αΦ)dk，則公式 (8)可以變?yōu)?/p>

dk+1=ΨHτk(ΨT(uk)).

(9)

因此，

uk+1=αdobs+(I-αΦ)dk+1

=αdobs+(I-αΦ)dk+1

=αdobs+(I-αΦ)ΨHτk(ΨT(uk))

=αdobs+(I-αΦ)ΨHτk(ΨT(uk)).

(10)這即為Gao等(2013)文中的加權(quán)POCS的迭代公式.

根據(jù)上述推導過程，權(quán)重因子α實際上為問題(6)中擬合誤差項的系數(shù). 若問題(6)兩邊除以α，則問題 (6)可以變?yōu)?/p>

(11)

很明顯，正則算子的系數(shù)λ變?yōu)棣?α.因此，加權(quán)POCS方法是解無約束優(yōu)化問題的方法，權(quán)重因子可以看作為擬合誤差項的系數(shù). 與之不同的是，原始POCS方法實際上是解等式約束優(yōu)化問題的方法， 這是兩者最大的不同之處，關(guān)于原始POCS方法的理論推導將在其他文章中闡述.

由公式(5)可知，POCS方法得到的解uN是dobs+(I-Φ)dN，其中dN=Ψ(xN)是時間-空間域的迭代解. 這里的(I-Φ)dN的作用是將采樣點位置的元素置為0，而未被采樣的點的位置的元素置為dN，則uN的采樣點位置的解為觀測數(shù)據(jù)dobs，未被采樣的位置的解為dN.如果觀測數(shù)據(jù)dobs含有噪聲， 這樣就會將噪聲帶入插值結(jié)果，因此POCS方法對于含噪聲數(shù)據(jù)是不合適的.加權(quán)POCS算法的解在采樣位置為αdobs+(1-α)dN，在未被采樣的位置為dN，采樣位置的解仍然受到含噪聲的觀測數(shù)據(jù)dobs的影響，此外該方法中除了需要認真選擇最小閾值外，還需要權(quán)重因子α來得到好的結(jié)果.由公式(6)可知，實際上權(quán)重因子可以消去，只要選擇合適的λ或最小閾值就可以實現(xiàn)同時插值和去噪，如果最小閾值選擇合適，并且返回dN為最終解，則可以實現(xiàn)同時插值和去噪，因此我們得到一種改進的POCS算法.該算法的迭代公式為

dk+1=ΨHτk(ΨT(dobs+(I-Φ)dk)),k=1,2,…,N，

(12)其中的迭代解為dk，屬于時間-空間域的解，當多次迭代后，返回dN為最終解.下面我們稱原始POCS為POCS(uN)，改進的POCS方法為POCS(dN)，對這兩個算法進行比較，驗證改進算法對含噪聲數(shù)據(jù)和不含噪聲數(shù)據(jù)插值的結(jié)果.

4 數(shù)值模擬

本節(jié)前兩個例子用來證明改進POCS 方法對無噪聲數(shù)據(jù)的插值效果；后兩個例子用來證明該方法對含噪聲數(shù)據(jù)插值的效果.我們用先插值后去噪的結(jié)果來驗證同時插值和去噪的可靠性.本文采用 curvelet 變換為稀疏變換(Candès et al.，2006).

4.1 無噪聲模擬數(shù)據(jù)的插值

圖1 (a)原始模擬數(shù)據(jù)；(b)隨機采樣的模擬數(shù)據(jù)Fig.1 The original synthetic data (a) and randomly sampled synthetic data (b)

圖2 (a) POCS (uN) 方法插值結(jié)果；(b) 圖2a與原始數(shù)據(jù)的誤差Fig.2 Interpolation of POCS (uN) method (a) and difference between Fig.2a and the original data (b)

圖3 (a) POCS(dN) 方法的插值結(jié)果；(b) 圖3a與原始數(shù)據(jù)的誤差Fig.3 Interpolation of POCS(dN) method (a) and difference between Fig.3a and the original data (b)

圖4 (a) 原始炮數(shù)據(jù)； (b) 隨機采樣的炮數(shù)據(jù)Fig.4 The original field data (a) and randomly sampled field data (b)

圖5 (a) POCS(uN) 插值結(jié)果； (b) 圖5a與原始數(shù)據(jù)的誤差Fig.5 Interpolation of POCS (uN) method (a) and difference between Fig.5a and the original data (b)

圖6 (a) POCS(dN) 方法的插值結(jié)果；(b) 圖6a與原始數(shù)據(jù)的誤差Fig.6 Interpolation of POCS (dN) method (a) and difference between Fig 6a and the original data (b)

圖7 (a) 噪聲模擬數(shù)據(jù)； (b) 含噪聲采樣數(shù)據(jù)Fig.7 The noisy synthetic data (a) and randomly sampled noisy synthetic data (b)

圖8 (a) POCS(uN)方法插值結(jié)果； (b) 圖8a與圖1a的誤差Fig.8 Interpolation of POCS(uN) method (a) and difference between Fig.8a and Fig.1a (b)

圖9 (a) POCS(uN)插值結(jié)果的第3道； (b) POCS(uN)插值結(jié)果的第5道Fig.9 The 3rd trace of POCS(uN) interpolation (a) and the 5th trace of POCS (uN) interpolation (b)

圖10 (a) POCS (dN) 方法的插值結(jié)果； (b) 圖10a與原始數(shù)據(jù)的誤差Fig.10 Interpolation of POCS (dN) method (a) and difference between Fig.10a and Fig.1a (b)

圖11 (a) POCS (uN)插值結(jié)果，最小閾值為 0.01； (b) 圖11a去噪結(jié)果Fig.11 Interpolation of POCS (uN) with least threshold 0.01 (a) and denoising result of Fig.11a (b)

圖12 (a) 含噪聲炮數(shù)據(jù)， 信噪比為5.1492 dB； (b) 采樣的數(shù)據(jù)， 信噪比為2.34 dBFig.12 The noisy field data， SNR=5.1492 dB (a) and randomly sampled noisy field data， SNR=2.34 dB (b)

圖13 (a) POCS(uN)方法插值結(jié)果； (b) 圖13a的殘差Fig.13 Interpolation of POCS (uN) method (a) and residual of Fig.13a (b)

圖14 (a) POCS(dN)方法的插值結(jié)果； (b) 圖14a的殘差Fig.14 Interpolation of POCS(dN) method (a) and residual of Fig.14a (b)

圖15 (a) POCS(uN)插值結(jié)果，最小閾值1； (b) 圖15a的去噪結(jié)果Fig.15 Interpolation of POCS (uN) method with least threshold 1 (a) and denoising result of Fig.15a (b)

表1 無噪聲模擬數(shù)據(jù)POCS(uN)方法和POCS(dN)方法的插值結(jié)果比較Table 1 Interpolation of POCS (uN) and POCS(dN) for noiseless synthetic data

這里POCS(uN)方法比POCS(dN)方法的結(jié)果稍微好一點，這是因為在采樣位置，POCS(uN)的插值結(jié)果被采樣dobs代替了，但是這僅僅對無噪聲數(shù)據(jù)是可行的.對于含噪聲數(shù)據(jù)這樣得到的解是不可接受的，下面的含噪聲數(shù)據(jù)插值的例子證明了這一點.4.3 含噪聲模擬數(shù)據(jù)的插值圖7a是含有零均值隨機噪聲的模擬數(shù)據(jù)，該數(shù)據(jù)的信噪比為4.1617 dB，圖7b是圖7a的采樣結(jié)果，信噪比為1.6023 dB，采樣道數(shù)為128.兩種算法的最小閾值為0.17，迭代步數(shù)N為 60. 圖8a是POCS (uN)的插值結(jié)果，由于仍然含有較嚴重噪聲，因此圖8a的結(jié)果不能令人滿意，圖8b是圖8a與原始數(shù)據(jù)(圖 1a)的誤差.圖9a給出了POCS(uN)插值結(jié)果的第3道，該道為采樣道，圖9b是POCS(uN)插值的第5道，這是一個未被采樣的道，可以看到采樣道的插值結(jié)果不是很好. 圖10a是POCS(dN) 的插值結(jié)果，圖10b是圖10a與原始數(shù)據(jù)的誤差. 兩種算法的CPU時間，信噪比和相對誤差見表3.

表2 無噪聲真實數(shù)據(jù)POCS(uN)方法和POCS(dN)方法插值結(jié)果比較Table 2 Interpolation of POCS(uN) and POCS(dN) for noiseless field data

為了驗證同時插值和去噪的可靠性，我們首先用POCS(uN)方法對圖7b進行插值，然后對插值結(jié)果去噪. POCS(uN)方法的最小閾值為0.01，最大迭代步數(shù)N為60. 插值結(jié)果見圖11a，信噪比為5.0561 dB.對圖11a采用閾值去噪的結(jié)果見圖11b，信噪比為12.6132 dB.很明顯，圖10a和圖11b的信噪比相當，因此當參數(shù)選擇合適時，同時插值和去噪是穩(wěn)定的.

表3 含噪聲合成數(shù)據(jù)POCS (uN)方法和POCS(dN)方法插值結(jié)果比較Table 3 Interpolation of POCS(uN) and POCS(dN) for noisy synthetic data

4.4 含噪聲真實地震數(shù)據(jù)的插值

圖12a是圖4a添加隨機白噪聲的結(jié)果，信噪比為5.1492 dB，圖12b是圖12a隨機采樣的結(jié)果，信噪比為2.34 dB. 選擇的最小閾值為16，最大迭代步數(shù)N為 60. 圖13a是POCS(uN)的插值結(jié)果，圖13b是圖13a的殘差.圖14a是 POCS(dN)的插值結(jié)果，圖14b是圖14a的殘差. 兩種算法的CPU 時間，信噪比SNR和相對誤差見表4. 由圖13a可知，POCS(uN)方法的插值結(jié)果是不可接受的，但是POCS(dN)方法可以得到較滿意的結(jié)果. 為了驗證同時插值和去噪的可靠性，圖15a是用POCS(uN)方法的插值結(jié)果，所選的最小閾值為1，信噪比為4.9458 dB，圖15b是圖15a閾值法的去噪結(jié)果，信噪比為8.9228 dB.上述結(jié)果存在一些噪聲，這是因為有效信號的部分curvelet系數(shù)和噪聲的curvelet系數(shù)耦合在一起了，僅僅通過閾值類去噪算法無法消除這種現(xiàn)象.

根據(jù)以上的數(shù)值模擬，POCS(dN)方法的計算效率和結(jié)果得到了驗證，其不僅僅適合去無噪聲數(shù)據(jù)的插值，而且適合含噪聲數(shù)據(jù)的插值.

表4 含噪聲真實數(shù)據(jù)POCS (uN)方法和POCS(dN) 方法插值結(jié)果比較Table 4 Interpolation of POCS (uN) and POCS(dN) for noisy field data

5 結(jié)論

本文由迭代閾值類方法推導出了加權(quán)POCS方法，證明其是解無約束優(yōu)化問題的一種方法，權(quán)重因子可以看作擬合誤差項的系數(shù)，同時也放大了最小閾值. 本文還提出了一種改進的 POCS方法來實現(xiàn)同時插值和去噪，該方法只含有一個超參數(shù)(閾值)，通過選擇合適的閾值就能夠?qū)崿F(xiàn)插值和去噪. 數(shù)值模擬表明，新的方法比原始POCS方法不增加任何計算量，在保持高效計算的同時能夠很好地插值無噪聲數(shù)據(jù)和含噪聲數(shù)據(jù). 通過和先插值后去噪結(jié)果的比較，證明了同時插值和去噪的可靠性.

基于迭代閾值類算法進行同時插值和去噪的關(guān)鍵是選擇合適的正則參數(shù)λ.λ和噪聲的能量有直接關(guān)系，雖然可以針對每個數(shù)據(jù)的噪聲能量來估計λ，但是這樣做效率較低，另外本文采用最大迭代步數(shù)來控制收斂，這個迭代步數(shù)是根據(jù)我們的經(jīng)驗選擇的，更實用的方法是選擇一種自適應的停機準則來控制迭代停止.為了便于我們的推導，我們采用二維數(shù)據(jù)進行驗證，本文的算法和結(jié)果對高維數(shù)據(jù)同樣適用.本算法不僅僅適用于curvelet變換，其他的稀疏變換也適用于該算法框架.未來的研究方向之一是實現(xiàn)含噪聲高維地震數(shù)據(jù)的插值和基于其他小波類變換的波場插值.關(guān)于POCS方法是解等式約束優(yōu)化問題的方法的證明將在其他文章中闡述.

Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm.Geophysics, 71(6): E91-E97.Berkhout A J, Verschuur D J. 1997. Estimation of multiple scattering by iterative inversion, Part I: Theoretical considerations.Geophysics, 62(5): 1586-1595.Blumensath T, Davies M E. 2008. Iterative thresholding for sparse approximations.JournalofFourierAnalysisandApplications, 14(5-6): 629-654. Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473.

Candès E, Demanet L, Donoho D, et al. 2006. Fast discrete curvelet transforms.SIAMJournalonMultiscaleModelingandSimulation, 5(3): 861-899.

Cao J J, Wang Y F, Zhao J T, et al. 2011. A review on restoration of seismic wavefields based on regularization and compressive sensing.InverseProblemsinScienceandEngineering, 19(5): 679-704.

Cao J J, Wang Y F, Yang C C. 2012. Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization.ChineseJournalofGeophysics(in Chinese), 55(2): 596-607, doi: 10.6038/j.issn.0001 5733.2012.02.22.Daubechies I, Defrise M, de Mol C. 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.CommunicationsonPureandAppliedMathematics, 57(11): 1413-1457.

Fomel S, Liu Y. 2010. Seislet transform and seislet frame.Geophysics, 75(3): V25-V38.

Gao J J, Stanton A, Naghizadeh M, et al. 2013. Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction.GeophysicalProspecting, 61(S1): 138-151.

Geng Y, Wu R S, Gao J H. 2012. Dreamlet compression of seismic data.ChineseJournalofGeophysics(in Chinese), 55(8): 2705-2715, doi: 10.6038/issn.0001-5733.2012.08.023.

Hennenfent G, Herrmann F J. 2006. Seismic denoising with nonuniformly sampled curvelets. Computing in Science and Engineering, 8: 16-25. Herrmann F J, Hennenfent G. 2008. Non-parametric seismic data recovery with curvelet frames.GeophysicalJournalInternational, 173(1): 233-248.Kreimer N, Sacchi M D. 2012. A tensor higher-order singular value decomposition for prestack seismic data noise reduction and interpolation.Geophysics, 77(3): V113-V122.

Liu B, Sacchi M D. 2004. Minimum weighted norm interpolation of seismic records.Geophysics, 69(6): 1560-1568.

Lu J M. 1993. Principle of Seismic Exploration (in Chinese). Dongying: Petroleum University Press. Ma J W. 2013. Three-dimensional irregular seismic data reconstruction via low-rank matrix completion.Geophysics, 78(5): V181-V192.Mansour H, Herrmann F J, Yilmaz O. 2013. Improved wavefield reconstruction from randomized sampling via weighted one-norm minimization.Geophysics, 78(5): V193-V206.

Oropeza V, Sacchi M D. 2011. Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis.Geophysics, 76(3): V25-V32.

Ronen J. 1987. Wave-equation trace interpolation.Geophysics, 52(7): 973-984.

Spitz S. 1991. Seismic trace interpolation in the f-x domain.Geophysics, 56(6): 785-794. Stolt R H. 2002. Seismic data mapping and reconstruction.Geophysics, 67(3): 890-908. Symes W W. 2007. Reverse time migration with optimal checkpointing.Geophysics, 72(5): SM213-SM221. Trad D O, Ulrych T J, Sacchi M D. 2002. Accurate interpolation with high-resolution time-variant Radon transforms.Geophysics, 67(2): 644-656. Wang Y F. 2007. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press. Wang Y F, Cao J J, Yang C C. 2011. Recovery of seismic wavefields based on compressive sensing by an l1-norm constrained trust region method and the piecewise random subsampling.GeophysicalJournalInternational, 187(1): 199-213.

Xu S, Zhang Y, Lambaré G. 2010. Antileakage Fourier transform for seismic data regularization in higher dimensions.Geophysics, 75(6): WB113-WB120.Yang P L, Gao J H, Chen W C. 2012. Curvelet-based POCS interpolation of nonuniformly sampled seismic records.JournalofAppliedGeophysics, 79: 90-99.

Yilmaz O. 1987. Seismic Data Processing, Investigations in Geophysics, Vol. 2. Tulsa, Okla: Society of Exploration Geophysicists. Zhang H, Chen X H. 2013. Seismic data reconstruction based on jittered sampling and curvelet transform.ChineseJournalofGeophysics(in Chinese), 56(5): 1637-1649, doi: 10.6038/cjg20130521.

Zwartjes P M. 2005. Fourier reconstruction with sparse inversion [Ph. D. Thesis]. Delft: Delft University of Technology.

附中文參考文獻

曹靜杰, 王彥飛, 楊長春. 2012. 地震數(shù)據(jù)壓縮重構(gòu)的正則化與零范數(shù)稀疏最優(yōu)化方法. 地球物理學報, 55(2): 596-607, doi: 10.6038/j.issn.0001 5733.2012.02.22.

耿瑜, 吳如山, 高靜懷. 2012. 基于Dreamlet變換的地震數(shù)據(jù)壓縮理論與方法. 地球物理學報, 55(8): 2705-2715, doi: 10.6038/issn.0001-5733.2012.08.023.

陸基孟. 1993. 地震勘探原理. 東營: 石油大學出版社.

王彥飛. 2007. 反演問題的計算方法及其應用. 北京: 高等教育出版社.

張華, 陳小宏. 2013. 基于jitter采樣和曲波變換的三維地震數(shù)據(jù)重建. 地球物理學報, 56(5): 1637-1649, doi: 10.6038/cjg20130521.

(本文編輯 胡素芳)

An improved projection onto convex sets method for simultaneous interpolation and denoising

CAO Jing-Jie1， WANG Ben-Feng2

1ShijiazhuangUniversityofEconomics，Shijiazhuang050031，China2StateKeyLaboratoryofPetroleumResourcesandProspecting，ChinaUniversityofPetroleum，Beijing102249，China

Seismic interpolation can provide a complete wavefield from incomplete wavefield information and is crucial for some processings such as multiple elimination, wave equation migration and denoising. Currently, most interpolation methods are suitable for high signal-to-noise ratio data. However, filed seismic data often contain different kinds of noise, especially the random noise. The existence of noise can deteriorate the result of interpolation. The weighted projection onto convex sets(POCS) method is a method for random noise contained data interpolation. However, there is no specific explanation of this method except the iteration formula. This paper analyzed the algorithm theory of the weighted POCS method, and proved that it is a method for unconstrained optimization problem, the weighting factor can be regarded as the coefficient of the fitting error. Meanwhile, we also proposed an improved POCS method for simultaneous interpolation and denoising. This algorithm requires no additional computational effort compared with the original POCS method.Seismic interpolation especially the high dimensional interpolation belongs to large-scale computational problems, thus fast algorithms should be chosen to improve the efficiency of seismic interpolation. Sparse transform methods change the interpolation problem into a sparse optimization based on the assumption that seismic data can be expressed sparsely in some transform domain. POCS method is such a method with high efficiency and robust results. However, this method can not get acceptable results for random noise contained data interpolation. The weighted POCS method is a method for random noise contained data interpolation on the basis of the POCS method. But there is no specific explanation of this method up to now. This paper deduced the weighted POCS method from the iterative thresholding method, we proved that the weighted POCS method is a method for unconstrained optimization, the weighted factor can be seen as the coefficient of the fitting error. However, the original POCS method is aiming at an equality constrained optimization. This is the essential difference between them. Meanwhile, we also proposed an improved POCS method with the same computational complexity as the original POCS method. The improved POCS method can realize interpolation and denoising simultaneously by adjustment of the least thresholds.In order to test the interpolation and denoising capability of the improved POCS method, four data sets were used to make numerical experiments. The first experiment is a synthetic noiseless data interpolation. The dimension of the synthetic data is 256×256 with trace distance 25 m and time sampling rate 4 ms. The original POCS method is chosen as the benchmark. Both the computational efficiency and numerical result are comparable for them. The second data is field data with 115 traces. The trace distance is 25 m and the time sampling rate is 4 ms. This experiment also proved that both of them can get comparable results for noiseless data interpolation. The third example is a noise contained data interpolation. The original POCS method can not get acceptable result no matter how to choose the least thresholds, while the improved POCS method can provide good result by choosing the least thresholds carefully. The fourth example also proved that the original POCS method failed in the noise contained field data interpolation, and the improved POCS method is suitable for noise contained data interpolation. Numerical experiment on first-interpolation-then-denoising proved the rationality and robustness of simultaneous interpolation and denoising.We deduced the weighted POCS method from the iterative thresholding method, and proved that it is a method for unconstrained optimization, the weighting factor is the coefficient of the fitting error. However, the original POCS method is a method for equality constrained optimization. This is the essential difference between them. We also proposed an improved POCS method to realize simultaneous interpolation and denoising. This method just need to adjust the least thresholds to realize noiseless data and random noise contained data interpolation. The computational effort of the improved POCS method is the same as the original POCS method but with more broad applications.

Sparse inversion； Denoise； Interpolation； Inversion

國家自然科學基金項目(41204075)，河北省自然科學基金(D2014403007)和河北省普通高等學校青年拔尖人才項目(BJ2014049)共同資助．

曹靜杰，男，副教授，1982年生，2011年于中國科學院地質(zhì)與地球物理研究所獲博士學位，主要從事地震信號處理和地球物理反演算法研究. E-mail：cao18601861@163.com

10.6038/cjg20150826.

10.6038/cjg20150826

P631

2014-09-17，2015-05-26收修定稿

曹靜杰，王本鋒. 2015. 基于一種改進凸集投影方法的地震數(shù)據(jù)同時插值和去噪.地球物理學報，58(8):2935-2947,

Cao J J, Wang B F. 2015. An improved projection onto convex sets method for simultaneous interpolation and denoising.ChineseJ.Geophys. (in Chinese),58(8):2935-2947,doi:10.6038/cjg20150826.