曹靜杰， 王本鋒

1 石家莊經(jīng)濟(jì)學(xué)院， 石家莊 050031 2 中國(guó)石油大學(xué)(北京)，油氣資源與勘探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249

?

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

曹靜杰1， 王本鋒2

1 石家莊經(jīng)濟(jì)學(xué)院， 石家莊 050031 2 中國(guó)石油大學(xué)(北京)，油氣資源與勘探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249

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

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

1 引言

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

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

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

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

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

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

Φd+n=dobs，

(1)

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

ΦΨx+n=dobs.

(2)

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

(3)

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

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

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

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

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

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

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

(5)

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

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

(6)

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

(7)

其中τk為第k個(gè)閾值參數(shù). 對(duì)公式(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)的推導(dǎo)可以參考文獻(xiàn)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ù)上述推導(dǎo)過(guò)程，權(quán)重因子α實(shí)際上為問(wèn)題(6)中擬合誤差項(xiàng)的系數(shù). 若問(wèn)題(6)兩邊除以α，則問(wèn)題 (6)可以變?yōu)?/p>

(11)

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

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

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

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

4 數(shù)值模擬

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

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

圖1 (a)原始模擬數(shù)據(jù)；(b)隨機(jī)采樣的模擬數(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) 隨機(jī)采樣的炮數(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 無(wú)噪聲模擬數(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é)果稍微好一點(diǎn)，這是因?yàn)樵诓蓸游恢?，POCS(uN)的插值結(jié)果被采樣dobs代替了，但是這僅僅對(duì)無(wú)噪聲數(shù)據(jù)是可行的.對(duì)于含噪聲數(shù)據(jù)這樣得到的解是不可接受的，下面的含噪聲數(shù)據(jù)插值的例子證明了這一點(diǎn).4.3 含噪聲模擬數(shù)據(jù)的插值圖7a是含有零均值隨機(jī)噪聲的模擬數(shù)據(jù)，該數(shù)據(jù)的信噪比為4.1617 dB，圖7b是圖7a的采樣結(jié)果，信噪比為1.6023 dB，采樣道數(shù)為128.兩種算法的最小閾值為0.17，迭代步數(shù)N為 60. 圖8a是POCS (uN)的插值結(jié)果，由于仍然含有較嚴(yán)重噪聲，因此圖8a的結(jié)果不能令人滿意，圖8b是圖8a與原始數(shù)據(jù)(圖 1a)的誤差.圖9a給出了POCS(uN)插值結(jié)果的第3道，該道為采樣道，圖9b是POCS(uN)插值的第5道，這是一個(gè)未被采樣的道，可以看到采樣道的插值結(jié)果不是很好. 圖10a是POCS(dN) 的插值結(jié)果，圖10b是圖10a與原始數(shù)據(jù)的誤差. 兩種算法的CPU時(shí)間，信噪比和相對(duì)誤差見表3.

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

為了驗(yàn)證同時(shí)插值和去噪的可靠性，我們首先用POCS(uN)方法對(duì)圖7b進(jìn)行插值，然后對(duì)插值結(jié)果去噪. POCS(uN)方法的最小閾值為0.01，最大迭代步數(shù)N為60. 插值結(jié)果見圖11a，信噪比為5.0561 dB.對(duì)圖11a采用閾值去噪的結(jié)果見圖11b，信噪比為12.6132 dB.很明顯，圖10a和圖11b的信噪比相當(dāng)，因此當(dāng)參數(shù)選擇合適時(shí)，同時(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í)地震數(shù)據(jù)的插值

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

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

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

5 結(jié)論

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

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

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.

附中文參考文獻(xiàn)

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

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

陸基孟. 1993. 地震勘探原理. 東營(yíng): 石油大學(xué)出版社.

王彥飛. 2007. 反演問(wèn)題的計(jì)算方法及其應(yīng)用. 北京: 高等教育出版社.

張華, 陳小宏. 2013. 基于jitter采樣和曲波變換的三維地震數(shù)據(jù)重建. 地球物理學(xué)報(bào), 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

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

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

10.6038/cjg20150826.

10.6038/cjg20150826

P631

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

曹靜杰，王本鋒. 2015. 基于一種改進(jìn)凸集投影方法的地震數(shù)據(jù)同時(shí)插值和去噪.地球物理學(xué)報(bào)，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.