孟閣閣,王德利,陳 鑫

(吉林大學地球探測科學與技術學院,吉林長春130026)

常規(guī)的地震資料去噪方法基于傅里葉變換[1],如頻率域濾波法、F-K濾波法等。但作為全局變換的傅里葉變換對整個時域進行平均,丟失了時間信息,對實際信號中一些非平穩(wěn)的奇異點或奇異曲線的時變局部特征描述不夠準確。作為傅里葉變換的延伸與發(fā)展,學者們提出了F-X域反褶積濾波[2]和小波變換[3]等算法。但是,F-X域反褶積濾波算法會加強地震記錄中的所有相干信息;而小波變換作為一種多尺度變換,雖然能夠更稀疏地表示一維分段函數(shù)或有界變差函數(shù),但對于高維信號的邊緣信息表達能力存在局限性。因此,需要有一種更為稀疏的變換來表示地震反射層。

Candès等[4]于1999年首次提出第一代曲波變換,它是由脊波變換衍生而來的一種多尺度幾何分析變換。曲波變換具有尺度、角度和位置分量,這使它不再僅僅局限于空間域和頻率域,也能夠在定向上實施,因此能夠對于含有奇異點的二維信號達到幾乎最優(yōu)的稀疏表示。至2004年,Candès等[5-6]又提出了第二代曲波變換,第二代的二維曲波變換無需用到脊波變換,提高了曲波變換的實現(xiàn)速度。在2004年的SEG年會上,開始出現(xiàn)曲波變換被用于地震數(shù)據(jù)多次波衰減[7]及曲波域的疊前偏移[8]等應用研究論文。由于單獨的二維曲波變換只能應用于單炮地震數(shù)據(jù),為了使曲波變換也可以應用于多炮地震數(shù)據(jù),Ying等[9]于2005年提出了三維離散曲波變換算法,為曲波變換用于分析和處理多炮數(shù)據(jù)和三維數(shù)據(jù)開啟了新的途徑。為了提高運算速度,減少計算過程中產生的冗余,Candès等[10]于2006年提出了兩種快速離散曲波變換算法。在2007年的SEG年會上,Herrmann等[11]介紹了基于曲波變換的地震數(shù)據(jù)缺失道插值方法,并將其應用于理論模型和實際三維數(shù)據(jù)體。2008年,Ramesh等[12]論述了一種基于曲波變換的隨機噪聲壓制方法。

在國內,隨著人們對曲波變換認識的加深,曲波變換在地震數(shù)據(jù)處理中的應用也逐漸引起重視。張素芳等[13]于2006年將曲波變換應用于地震數(shù)據(jù)多次波消除;2008年,彭才等[14]給出了基于曲波變換的地震數(shù)據(jù)去噪方法;同年,王德利等[15]提出了一種基于曲波變換的稀疏控制反演方法;張恒磊等[16]于2010年將曲波變換應用于地震波場分離中。至2013年,曲波變換已先后被應用于稀疏反褶積[17]和地層吸收補償[18]等問題中;馮飛等[19]還提出了基于三維曲波變換的稀疏反演一次波估計方法。總體來看,基于三維曲波變換的方法研究和應用研究在國內還處于起步階段。

基于三維曲波變換良好的稀疏性,可以用更少的系數(shù)立體全面地重構地震反射信號,我們將三維曲波變換與閾值迭代法結合,將其應用于多炮地震數(shù)據(jù)的隨機噪聲壓制中。在曲波域中用曲波系數(shù)描述地層反射系數(shù),在稀疏條件下地震數(shù)據(jù)中隨機噪聲的衰減問題可以描述為L1范數(shù)最優(yōu)化問題,對此用閾值迭代法進行求解,以達到壓制隨機噪聲的目的。采用多炮理論模型數(shù)據(jù)和實際地震資料對基于三維曲波變換的閾值迭代去噪算法進行驗證,并選用信噪比、保真度等評價參數(shù)與幾種傳統(tǒng)去噪方法的效果進行量化對比分析。

1 三維曲波變換的基本理論

曲波變換用“楔形基”代替小波變換的“塊基”來逼近二次連續(xù)可微的奇異點,所以能夠更好地描述沿邊緣的信息。三維曲波變換[9]是二維曲波變換的延伸與擴展,在尺度、角度和位置分量上具有與二維曲波變換相同的性質。

在頻率域中,二維曲波變換的支撐區(qū)域是一個“楔形”,三維曲波變換的頻率支撐建立在具有拋物型尺度屬性的楔形附近。連續(xù)三維曲波變換的基函數(shù)在頻率域中以一個球體的形式逐漸分割覆蓋(圖1),徑向上使用與二維曲波變換相同的窗口W(r)：

(1)

圖1 連續(xù)三維曲波變換分割示意圖解[9]

離散三維曲波變換的基函數(shù)在頻率域中以方塊體的形式逐漸分割覆蓋(圖2)。類似于離散二維曲波變換,離散三維曲波變換輸入的是三維笛卡爾網(wǎng)格形式f(n1,n2,n3),0≤n1,n2,n3≤n,輸出曲波系數(shù)集cD(j,l,k)：

(2)

其中,j,l∈Z,k=(k1,k2,k3)。

用笛卡爾半徑的形式表示徑向窗函數(shù)：

并且,φj(ω1,ω2,ω3)=φ(2-jω1)φ(2-jω2)φ(2-j·ω3)。其中,φ是低通圓滑函數(shù);φ是曲波變換的“母波形”。

離散三維曲波變換的基函數(shù)分割體有6個面,在精細尺度下,以第1部分(ω1>0)為例:假設(1,αl,βl)是楔形中心線的方向,角窗口函數(shù)定義為

(5)

對于其它的5個部分,交換ω1與ω2或ω1與ω3來進行類似定義。

圖2 離散三維曲波變換分割示意圖解[9]

(6)

位移參數(shù)為k的離散曲波在尺度j和角度l上由它們傅里葉變換的方式給出

(7)

其中,0≤k1

離散三維曲波變換算法實現(xiàn)步驟如下。

2 三維曲波閾值迭代法

2.1 三維曲波閾值迭代法實現(xiàn)基礎

一般的閾值去噪算法都分為信號分解、系數(shù)閾值處理、信號重構3個步驟,基于三維曲波變換的閾值迭代算法也不例外。曲波變換良好的稀疏性決定了它的優(yōu)越性,如圖3所示,將曲波變換的“楔形”基用長方形窗口近似代替,這種長方形窗口也體現(xiàn)了曲波基特有的方向性。圖3中曲線表示地震反射層。小波變換(圖3a)用類似正方形的窗口來描述地層曲線;而曲波變換(圖3b)的窗口函數(shù)是類似于長方形的。對比圖3a和圖3b可以看出,小波變換需要較多的窗口才能描述地層,并且這種“正方形”窗口函數(shù)不具有方向性,而曲波變換只用少數(shù)的曲波系數(shù)就能很好地描述地層反射信號,“長方形”窗口函數(shù)的方向性使有效信號主要集中在沿曲線方向的曲波系數(shù)上,即小傾角的曲波系數(shù)上,而隨機噪聲在任何正交基上的變換仍然是隨機地分布在整個空間上。一般來說噪聲的曲波系數(shù)值由于能量分散且幅值較小而往往小于有效信號的曲波系數(shù)值,因此只要選擇一個適當?shù)拈撝?對曲波系數(shù)進行閾值處理,就可以達到有效信號與噪聲分離的目的。

圖3 小波逼近(a)與曲波曲線逼近(b)示意圖解

2.2 三維曲波閾值迭代法實現(xiàn)步驟

Daubechies等[20]給出了一種稀疏約束下的正則化方法的求解迭代序列：

(8)

其中,

(9)

稀疏約束情況下,Elad等[21]提出了地震數(shù)據(jù)中隨機噪聲的衰減可以簡化表述為如下的約束優(yōu)化公式:

(10)

式中：y是待處理的含噪地震數(shù)據(jù);λ為拉格朗日算子;A表示三維曲波逆變換;x表示曲波系數(shù);f為最終輸出結果。Pλ問題實際上是求解L1范數(shù)最優(yōu)化的過程,我們采用軟閾值迭代法對其進行求解。

根據(jù)(8)式可知其迭代序列為

(11)

式中：AT表示三維曲波正變換;Tλ為軟閾值函數(shù),其表達式為

(12)

獲得良好去噪效果的關鍵是閾值的選取,我們根據(jù)以下兩個公式來選取閾值：

三維曲波閾值迭代法去噪處理的步驟見圖4 所示,圖中i表示當前外循環(huán)次數(shù);j表示當前內循環(huán)次數(shù);L表示總的外循環(huán)次數(shù);l表示總的內循環(huán)次數(shù);x表示曲波系數(shù)。

圖4 三維曲波閾值迭代法去噪處理流程

3 計算實例

3.1 模擬數(shù)據(jù)測試算例

為了說明三維曲波變換較其它常規(guī)方法擁有更好的去噪效果和保真能力,先對模擬多炮地震數(shù)據(jù)進行測試。所用模擬數(shù)據(jù)有181炮,每炮有181道記錄,共1125個采樣點,采樣間隔為0.0016s。在模擬原始數(shù)據(jù)中加入不同程度的隨機噪聲,針對高信噪比模擬數(shù)據(jù)和低信噪比模擬數(shù)據(jù)兩種情況,分別采用傳統(tǒng)中值濾波方法、F-X反褶積方法、二維曲波閾值迭代法和三維曲波閾值迭代法進行去噪處理。

圖5a給出了模擬原始數(shù)據(jù)的第83炮原始記錄;圖5b是加噪后高信噪比模擬記錄;圖5c到圖5f 分別給出了上述4種方法的去噪處理結果。

為了計算去噪后數(shù)據(jù)的信噪比和保真度,更好地對比不同方法的去噪效果,定義如下評價參數(shù)。

1) 信噪比:

(15)

2) 峰值信噪比：

圖5 模擬原始數(shù)據(jù)(a)和加噪后高信噪比數(shù)據(jù)(SNR=5.5691dB)(b)的中值濾波法(c),F-X反褶積法(d),二維曲波閾值迭代法(e)和三維曲波閾值迭代法(f)去噪處理結果

(16)

式中：s0表示原始無噪聲數(shù)據(jù);s1表示加噪聲數(shù)據(jù)或去噪后數(shù)據(jù);M代表數(shù)據(jù)的行數(shù);N代表數(shù)據(jù)的列數(shù)。

3) 均方差：

(17)

4) 平均絕對誤差:

(18)

式中：v0i為原始無噪數(shù)據(jù);vi為去噪后數(shù)據(jù)。

5) 去噪能力：

(19)

式中:SNR0為原始含噪數(shù)據(jù)的信噪比;Enoise為去噪后數(shù)據(jù)中的噪聲能量;Enoise0為原始含噪數(shù)據(jù)中的噪聲能量。

表1給出了對含噪高信噪比模擬數(shù)據(jù)利用4種方法去噪處理結果計算的信噪比SNR,峰值信噪比PSNR,均方差MSE,平均絕對誤差MAE和去噪能力DA的對比。

由圖5和表1可以看出,對于圖5b所示信噪比較高的含噪模擬數(shù)據(jù)(SNR=5.5691dB;PSNR=33.3821dB),中值濾波法去噪結果(圖5c)效果差強人意(SNR=5.9236dB;PSNR=33.73631dB),

有效信息損失嚴重(MSE=0.455×10-5)。F-X反褶積法去噪結果(圖5d)由于有效信號與相干噪聲同時加強,所以信噪比提高量不大(SNR=6.7077dB;PSNR=34.5207dB),保真效果也不太理想(MSE=0.3790×10-5)。二維曲波閾值迭代法取得了不錯的去噪效果(SNR=13.5072dB;PSNR=41.3202dB;MSE=0.0790×10-5)(圖5e),但是同相軸曲率較大處信息模糊。經(jīng)三維曲波閾值迭代法去噪后(圖5f),信噪比得到了顯著提升(SNR=15.9672dB;PSNR=43.7801dB),同相軸連續(xù)、清晰,且有效信息損失也是最小的(MSE=0.0450×10-5),獲得了滿意的效果。

圖6a和圖6b分別為采用二維曲波閾值迭代法和三維曲波閾值迭代法去噪前、后的差值剖面,在圖中箭頭所示的位置可以明顯地看出,經(jīng)二維曲波閾值迭代去噪后的差值剖面上同相軸殘余較多,而經(jīng)三維曲波閾值迭代去噪后的差值剖面上有效信息明顯減少?？梢?三維曲波閾值迭代法不僅去噪效果好,而且有效信息的損失相對較小。

進一步對低信噪比模擬數(shù)據(jù)進行測試。為了便于比較,這里仍然針對其中的第83炮模擬原始記錄(圖7a)。圖7b是加噪后的低信噪比模擬數(shù)據(jù)(SNR=-8.4118dB),分別采用中值濾波法(圖7c),F-X反褶積法(圖7d),二維曲波閾值迭代法(圖7e)和三維曲波閾值迭代法(圖7f)進行去噪處理。表2是對含噪低信噪比模擬數(shù)據(jù)利用4種方法去噪處理結果計算的信噪比SNR,峰值信噪比PSNR,均方差MSE,平均絕對誤差MAE和去噪能力DA。

表1 高信噪比模擬數(shù)據(jù)4種方法去噪結果的評價參數(shù)

圖6 高信噪比模擬數(shù)據(jù)二維曲波閾值迭代法(a)和三維曲波閾值迭代法(b)去噪前、后的差值剖面

圖7 模擬原始數(shù)據(jù)(a)和加噪后低信噪比數(shù)據(jù)(SNR=-8.4118dB)(b)的中值濾波法(c),F-X反褶積法(d),二維曲波閾值迭代法(e)和三維曲波閾值迭代法(f)去噪處理結果

中值濾波 F-X反褶積二維曲波三維曲波SNR/dB-7.9800 0.4805 4.59218.3245PSNR/dB 19.8324 28.293532.4052 36.1375MSE 1.1170×10-4 0.1590×10-40.0620×10-4 0.0260×10-4MAE 0.7530×10-20.2780×10-20.1660×10-2 0.1100×10-2DA0.1560 1.61302.8080 3.6380

由圖7和表2可以看出,對于含有較多隨機噪聲的地震數(shù)據(jù)(SNR=-8.4118dB;PSNR=19.4021dB),中值濾波法去噪的效果(圖7c)不太理想,剖面中仍有大量的噪聲(SNR=-7.980dB;PSNR=19.8324dB),并且有效信號損失較多(MSE=1.1170×10-4)。經(jīng)F-X反褶積法去噪后(圖7d),信噪比雖然有所提高(SNR=0.4805dB;PSNR=28.1935dB),但是有效信號信息明顯減弱(MSE=0.1590×10-4)。二維曲波閾值迭代法結果(圖7e)雖然去除了大量噪聲(SNR=4.5921dB;PSNR=32.4052dB),剖面視覺效果優(yōu)于圖7c和圖7d,但是從圖7e中0.5～1.0s兩個箭頭處可以看出,二維曲波閾值迭代法在同相軸曲率較大處損失較大(MSE=0.0620×10-4),出現(xiàn)了假頻;從圖7e中1.5s上、下的兩個箭頭處可以看出,二維曲波閾值迭代法很難保證弱有效信號的真實信息,不利于后續(xù)數(shù)據(jù)處理工作的進行;從圖7e 整體可以看出,隨著信噪比的降低,二維曲波閾值迭代法去噪后的剖面背景中存在一些明顯的小曲波分量,這將干擾后續(xù)的剖面解釋。而經(jīng)三維曲波閾值迭代法去噪后(圖7f),信噪比得到很大的提高(SNR=8.3245dB;PSNR=36.1375dB),有效信號的損失還是最小(MSE=0.1590×10-4)。

圖8a和圖8b分別為采用二維曲波閾值迭代法和三維曲波閾值迭代法去噪前、后的差值剖面?？梢悦黠@看出,經(jīng)二維曲波閾值迭代法去噪后的差值剖面上存有較明顯的同相軸殘余,而經(jīng)三維曲波閾值迭代法去噪后的差值剖面上看不出明顯的有效信息?？梢娂词乖诘托旁氡鹊那闆r下,三維曲波閾值迭代法的去噪效果仍然要好于其它3種方法,更符合高信噪比、高保真度的地震資料去噪處理要求。

通過以上模擬數(shù)據(jù)測試結果的對比分析可以看出,相比于傳統(tǒng)的去噪算法,無論原始數(shù)據(jù)的信噪比高低與否,三維曲波閾值迭代法都能獲得令人滿意的去噪效果。

3.2 實際地震資料處理試驗

選取某海區(qū)的實際含噪地震資料進行三維曲波閾值迭代法去噪處理試驗。該實際地震資料共236炮,每炮有236道記錄,采樣點數(shù)為1250個。圖9a為其中的第130炮原始記錄,由圖9a可見,聲的痕跡,而三維曲波閾值迭代法處理結果由于隨機噪聲的存在,原始地震記錄信噪比較低,2.0s以下同相軸較難分辨,而且記錄中還存在大量涌浪噪聲。圖9b到圖9e分別為經(jīng)中值濾波法、F-X反褶積法、二維曲波閾值迭代法和三維曲波閾值迭代法去噪處理后的結果。圖10a和圖10b分別為該實際炮記錄用二維曲波閾值迭代法和三維曲波閾值迭代法去噪前、后的差值剖面。

圖8 低信噪比模擬數(shù)據(jù)二維曲波閾值迭代法(a)和三維曲波閾值迭代法(b)去噪前、后的差值剖面

在采用中值濾波法去噪后的剖面上(圖9b)幾乎很難分辨出同相軸,其處理結果基本上不能用于實際生產。F-X反褶積法處理結果(圖9c)比中值濾波算法的去噪效果要好一些,但是中部弱反射信息被嚴重剝蝕。對比二維曲波閾值迭代法(圖9d)和三維曲波閾值迭代法(圖9e)去噪處理結果可以看出：①在0.5～1.0s的白色箭頭所指處,即同相軸的邊緣部分,經(jīng)二維曲波閾值迭代法處理后出現(xiàn)了假頻,產生了較大的模糊,從圖10的差值剖面上也可以看出這一點;②對比最底部的兩處藍色箭頭所指處可見,經(jīng)二維曲波閾值迭代法去噪后的剖面中有明顯的大傾角曲波分量存在;③原始記錄在160道附近有一個豎直的涌浪噪聲(圖中紅色箭頭所示),在圖9d中仍然可以看出涌浪噪聲的痕跡,而三維曲波閾值迭代法處理結果(圖9e)在去除隨機噪聲的同時也一并去除了涌浪噪聲,同相軸清晰連續(xù);④原始記錄中與橫坐標成銳角方向有一規(guī)則弱反射信息(圖中黃色箭頭所示),這在原始含噪聲剖面中很難分辨,而在圖9e中這一規(guī)則干擾被清晰地反映了出來,這也體現(xiàn)了三維曲波閾值迭代去噪算法優(yōu)于其它3種算法,這一點對于實際生產中提取弱反射信息至關重要。

圖9 某海區(qū)實際地震資料原始炮記錄(a)及其中值濾波法(b),F-X反褶積法(c),二維曲波閾值迭代法(d)和三維曲波閾值迭代法(e)去噪處理結果

綜上所述,對于多炮地震數(shù)據(jù)的隨機噪聲壓制,三維曲波閾值迭代算法是一種有效的去噪方法,去噪效果明顯優(yōu)于其它傳統(tǒng)算法。

圖10 實際地震炮記錄二維曲波閾值迭代法(a)和三維曲波閾值迭代法(b)去噪前、后的差值剖面

5 結束語

三維曲波變換是一種多尺度分析方法,它增加的尺度、角度和位置分量給三維曲波變換帶來了其它變換不可及的稀疏性,使閾值迭代法的求解精度更高。經(jīng)模擬數(shù)據(jù)和實際地震資料去噪處理結果的對比分析可以看出,三維曲波閾值迭代法較常規(guī)去噪方法不僅去噪效果突出,而且去噪后有用信息損失較少,是一種有效的多炮地震數(shù)據(jù)隨機噪聲去除方法。

同時也應該指出,由于三維曲波變換會產生較高的冗余度,計算三維數(shù)據(jù)所需的計算時間和存儲空間都是巨大的,這一問題限制了其在實際生產中的應用。如何進一步簡化其結構,管理數(shù)據(jù),減少冗余度,還有待進一步研究。此外,我們采用了三維曲波變換與軟閾值迭代法結合的算法,雖然去噪效果比較滿意,但是由于軟閾值算法的先天缺陷可能會造成信號邊緣模糊,不利于具有尖銳邊緣的地震資料的處理,因此還需要尋找一種更好的閾值選取方法。

參 考 文 獻

[1] 張軍華,呂寧,田連玉,等.地震資料去噪方法綜合評述[J].石油地球物理勘探,2005,40(增刊)：121-127

Zhang J H,Lv N,Tian L Y,et al.Comprehensive review of seismic data denoising[J].Oil Geophysical Prospecting,2005,40(S1)：121-127

[2] 康冶,于承業(yè),賈臥,等.f-x域去噪方法研究[J].石油地球物理勘探,2003,38(2):136-138

Kang Y,Yu C Y,Jia W,et al.f-x domain denoising research[J].Oil Geophysical Prospecting,2003,38(2):136-138

[3] 王勇,酈軍,張奎鳳.基于小波變換的地震信號降噪處理[J].石油物探,1998,37(3)：72-76

Wang Y,Li J,Zhang K F.The noise elimination processing of seismic signal based on wavelet transform[J].Geophysical Prospecting for Petroleum,1998,37(3)：72-76

[4] Candès E,Donoho D.Curvelets:a surprisingly effective nonadaptive representation of objects with edges[M].Tennessee Nashville:Vanderbilt University Press,1999：1-10

[5] Candès E J,Donoho D L.New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities[J].Communications on Pure and Applied Mathematics,2004,57(2):219-266

[6] Candès E,Guo F.New multiscale transforms,minimum total variation synthesis:applications to edge-preserving image reconstruction[J].Signal Processing,2002,82:1519-1543

[7] Herrmann F J,Verschuur E.Curvelet-domain multiple elimination with sparseness constraints[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1333-1336

[8] Huub D,Maarten V.Wave-character preserving pre-stack map migration using curvelets[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,961-964

[9] Ying L,Demanet L,Candès E.3-D discrete curvelet transform[J].International Society for Optics and Photonics,2005,5914(13):344-354

[10] Candès E,Demanet L,Donoho D,et al.Fast discrete curvelet transforms[J].SI-AM Multi-scale Modeling and Simulation,2006,5(3):861-899

[11] Herrmann F J,Wang D L,Hennenfent G.Multiple prediction from incomplete data with the focoused curvelet transform[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007,2505-2509

[12] Ramesh N,Domimique G,Mohamed T,et al.Coherent and random noise attenu-ation using the curvelet transform[J].The Leading Edge,2008,27(2):240-248

[13] 張素芳,徐義賢,雷棟.基于Curvelet變換的多次波去除技術[J].石油地球物理勘探,2006,41(3)：262-265

Zhang S F,Xu Y X,Lei D.Multiple-eliminated technique based on Curvelet transform[J].Oil Geophysical Prospecting,2006,41(3):262-265

[14] 彭才,常智,朱仕軍.基于曲波變換的地震數(shù)據(jù)去噪方法[J].石油物探,2008,47(5):461-464

Peng C,Chang Z,Zhu S Y.Noise elimination method based on curvelet transform[J].Geophysical Prospecting for Petroleum,2008,47(5):461-464

[15] Wang D L,Rayan S,Yilmaz?,et al,Bayesian wavefield separation by transform-domain sparsity promotion[J].Geophysics,2008,73(5):A33-A38

[16] 張恒磊,劉天佑.Curvelet變換及其在地震波場分離中的應用[J].煤田地質與勘探,2010,38(1)：76-80

Zhang H L,Liu T Y.Curvelet transform and its application in seismic wave field separation[J].Coal Geology & Exploration,2010,38(1)：76-80

[17] 孟大江,王德利,馮飛,等.基于Curvelet變換的稀疏反褶積[J].石油學報,2013,34(1):107-114

Meng D J,Wang D Li,Feng F,et al.Spare deconvolution based on the Curvelet transform[J].Acta Petrolei Sinica,2013,34(1):107-114

[18] Wang D,Sun J,Meng D,et al.Absorption compensation based on curvelet transform[J].Journal of Seismic Exploration,2013,22(1)：19-23

[19] Feng F,Wang D L,Zhu H,et al.Estimating primaries by sparse inversion of the 3D Curvelet transform and the L1-norm constraint[J].Applied Geophysics,2013,10(2)：201-209

[20] Daubechies I,Defrise M,de Mol C.An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J].Communications on Pure and Applied Mathematics,2004,57(11):1413-1457

[21] Elad M,Starck J L,Querre P,et al.Simultaneous cartoon and texture image inpainting using morphological component analysis(MCA)[J].Applied and Computational Harmonic Analysis,2005,19(3):340-358