陳思雅，謝凱，阮寧君，張龍,龔康奕，沈政春，李紀成

(長江大學電子信息學院，湖北荊州434023)

?

一種曲波域蒙特卡羅閾值去噪算法

陳思雅，謝凱，阮寧君，張龍,龔康奕，沈政春，李紀成

(長江大學電子信息學院，湖北荊州434023)

[摘要]為促進油氣勘探發(fā)展，有效識別地震有用信號，提出了一種曲波域蒙特卡羅閾值去噪算法。該算法利用曲波變換的多尺度與多方向性的特點，結(jié)合蒙特卡羅閾值濾波來去除隨機噪聲，再利用循環(huán)平移法抑制吉布斯現(xiàn)象，消除人為干擾，增強去噪效果。用該算法對合成地震數(shù)據(jù)與實際地震剖面分別進行去噪試驗，并與小波變換和曲波循環(huán)平移算法的試驗結(jié)果進行對比，證明了該算法能有效地保持有效波信息，增加了弱反射信號能量，很好識別有效信號，優(yōu)于小波變換去噪算法與原始曲波變換去噪算法。

[關(guān)鍵詞]循環(huán)平移；曲波變換；隨機噪聲；蒙特卡羅閾值

隨著油氣勘探開發(fā)的不斷深入，勘探開發(fā)難度不斷增加，地震勘探的任務(wù)也已由過去單純的構(gòu)造勘探發(fā)展到尋找復(fù)雜油氣藏和隱蔽油氣藏。隱蔽油氣藏[1]具有面積小、構(gòu)造復(fù)雜，在地震剖面上特征顯示不明顯的特點，因此在地震勘探中深入研究對地震弱信號的檢測與識別，提高分辨率等問題是必然的趨勢。在地震勘探中，隨機噪聲是一種頻帶較寬的干擾波，常規(guī)的去噪方法不理想。小波變換由于其良好的空域和頻域的局部特性，且對信號的奇異值敏感[2]，能反映奇異點的位置和特征。目前傳統(tǒng)的去噪方法有中值濾波法、頻域濾波法、F-K濾波法、F-X反濾波法、多項式擬合法、傅里葉變換法、Radon變換法、奇異值分解法、波動方程法及最為常用的小波變換法[3～6]。小波變換在地震信號處理實驗中雖然能夠去除大部分噪聲，但是，由于小波變換本身的缺陷，使得去噪后的地震信號同相軸的邊緣變得模糊不清、分辨率下降[7]；并且，小波變換表達二維信號時會損失邊緣細節(jié)信息，造成成像模糊[8]。針對小波變換的局限性，國外學者提出了一種具有多尺度、多方向的分析方法—曲波變換Curvelet[9，10]。1999年Candès和Donoho提出了曲波變換理論，它在對線和超平面的奇異性表示上具有優(yōu)于小波變換的特性。曲波變換是各向異性，具有很強的方向性，能同時獲得圖像平滑區(qū)域和邊緣部分的稀疏表達，且以邊緣為基本表達元素，具有完備性，能用更少的非零系數(shù)精確地表達圖像的邊緣，能更好地適合圖像的特點。將曲波分析應(yīng)用在圖像等數(shù)據(jù)的去噪上，利用其多尺度多分辨率特點，再結(jié)合閾值去噪算法，能獲得很好的效果。曲波閾值去噪法[11]的實現(xiàn)方案中存在著一定的問題，由于變換中使用的是在頻域中具有緊支集的小波函數(shù)，在頻域中計算小波變換的系數(shù)，并進行衰減，使用的也是硬閾值函數(shù)[12]，因此在去噪后的圖像中呈現(xiàn)出明顯的振鈴效應(yīng)[13]。為了減弱振鈴效應(yīng)，采用對去噪后圖像進行維納濾波[14]，但這是以模糊圖像邊緣為代價的。循環(huán)平移法[15]根據(jù)曲波自身時間周期平移不變性[16]，可以較好地消除這種效應(yīng)，保護圖像邊緣信息。為了有效去除深部隨機噪聲，筆者提出了一種基于曲波域的蒙特卡羅閾值去噪方法[17]。

1基于曲波域的蒙特卡羅閾值去噪算法

該算法分為5步，如圖1所示。

圖1　算法流程圖

1.1曲波變換

二維空間中，設(shè)空間域變量為x,頻率域變量為w,頻率域極坐標為r和θ,l為方向參數(shù)。假設(shè)存在光滑、非負且實值的半徑窗與角度窗，且滿足容許性條件：

(1)

(2)

則對每個j，頻域窗Uj定義為：

(3)

令母Curvelet函數(shù)為：

(4)

(5)

式中，旋轉(zhuǎn)角度θl=2π×2-j/2,l=0,1，…，且0≤θl≤2π；平移參數(shù)k=(k1，k2)∈Z2,Rθ表示角度為θ時的旋轉(zhuǎn)矩陣。則函數(shù)f(x)∈L2(R2)的Curvelet變換可表示為其與其他Curvelet函數(shù)φj,l,k(x)的內(nèi)積：

(6)

其在頻率域中的定義為：

(7)

曲波重構(gòu)公式定義為：

(8)

式中，〈f(x),cμ〉是函數(shù)f(x)的Curvelet系數(shù)；cμ表示由多指標μ=(j,l,k1,k2)確定的Curvelet族；M是指標積。

1.2循環(huán)平移增強去噪效果

不論在小波變換還是曲波變換去噪時，如果變換缺乏平移不變性，則會在不連續(xù)點領(lǐng)域產(chǎn)生明顯的人工效應(yīng)，導致信號失真。這種信號失真與不連續(xù)點的位置密切相關(guān)，即對某些不連續(xù)點處理得好，在其領(lǐng)域內(nèi)就不會產(chǎn)生吉布現(xiàn)象，但在其他相鄰近的不連續(xù)點的領(lǐng)域內(nèi)卻會產(chǎn)生明顯的吉布斯現(xiàn)象。為消除這種效應(yīng)，采用了循環(huán)平移方法，借助曲波變換里小波的時間周期平移不變性，對原信號數(shù)據(jù)組成的二維矩陣沿行和列方向進行平移，每次移動一行一列，最大平移次數(shù)為行數(shù)乘以列數(shù)。對平移后的信號作去噪處理，再反平移回去，在進行多次循環(huán)平移后，對所有得到的去噪后的信號作平均。由于每次平移后的信號作閾值去噪處理后使得吉布斯現(xiàn)象出現(xiàn)在不同地方，因此通過對行和列方向上的每組平移量得到的去噪結(jié)果作線性平均，可以在一定程度上消除這種振鈴效應(yīng)。

1.3蒙特卡羅閾值

小波去噪對噪聲方差的估計用得較多的是中值估計，在曲波去噪用得較多的是蒙特卡羅閾值。蒙特卡羅是一種以概率統(tǒng)計理論為指導的一類非常重要的數(shù)值計算方法，是使用隨機數(shù)(或更常見的偽隨機數(shù))來解決很多計算問題的方法，其公式為T=k×σ×σj[18]。式中k是一個依賴于尺度和方向的常數(shù)，不同的尺度和方向下的曲波可取不同的值。利用蒙特卡羅方法自適應(yīng)估計各個尺度方向上的閾值，可一定程度上避免常規(guī)統(tǒng)一閾值對有效信號的影響。在最粗尺度上人為選取閾值比例系數(shù)k為2.2，在最細尺度上k取2.5。σ為加入圖像中的噪聲空間域標準差的估計值，在這里取為20。σj為對均值為0，方差為1的高斯白噪聲進行曲波變換后，進行蒙特卡羅測試得到的每一尺度對應(yīng)的系數(shù)標準差的近似值。由于曲波分析不像傅里葉變換或者小波變換是范數(shù)不變的，噪聲的離散曲波系數(shù)與曲波的參數(shù)有關(guān)，如用F表示曲波矩陣，白噪聲服從均值為0，方差為1的正態(tài)分布，那么通過蒙特卡羅仿真法來計算系數(shù)個體的方差近似值σj，通過簡單計算一些標準白噪聲圖像的曲波變換來獲得。

在這里取的是中心元素為1，其他元素全為0的規(guī)范矩陣，來代表均值為0，方差為1的標準白噪圖像。因為加入的隨機噪聲服從正態(tài)分布，所以在進行曲波變換后，噪聲的曲波系數(shù)也服從正態(tài)分布且均值為0。可以通過對規(guī)范矩陣求標準差來估計含噪圖像里的噪聲標準差，從而確定閾值。具體計算方法為：求取s尺度w方向的系數(shù)矩陣的范數(shù)，再除以該矩陣包含的元素個數(shù)，就得到了噪聲的標準差的估值。根據(jù)上述公式求出的T就是在不同尺度和方向上的蒙特卡羅閾值。

2試驗結(jié)果分析

2.1模型數(shù)據(jù)試驗

圖2(a)是合成地震剖面，圖2(b)在圖2(a)的基礎(chǔ)上進行加噪處理已形成含隨機噪聲的合成地震剖面，用3種不同的去噪算法分別對合成地震剖面進行處理，通過視覺比較以及峰值信噪比對比處理效果。其中，峰值信噪比PNSR計算公式如下：

(9)

(10)

式中，MSE是均方差，I和K分別代表含噪模型數(shù)據(jù)和去噪后的模型數(shù)據(jù)，兩者都是大小均為500×50的二維信號。

圖2　模型數(shù)據(jù)試驗對比圖

圖3　峰值信噪比對比圖

圖2以原始模型數(shù)據(jù)圖2(a)作為參照，在其進行加噪處理后得到圖2(b)。其中圖2(c)、(d)、(e)分別是進行小波去噪、原始曲波去噪及曲波循環(huán)平移去噪得的結(jié)果。與原始模型數(shù)據(jù)圖2(a)相比，小波去噪效果圖2(c)中噪聲依然十分明顯，而原始曲波去噪效果圖2(d)與曲波循環(huán)平移去噪效果圖2(e)與原始模型數(shù)據(jù)十分接近，噪聲很少，這說明曲波去噪效果比小波效果好很多。將曲波循環(huán)平移效果圖2(e)和原始曲波變換去噪效果圖2(d)對比，前者去噪效果更好，從圖中畫圈標記看出振鈴效應(yīng)也得到了改善。圖3是峰值信噪比對比圖，從圖3中可以得知，在相同噪聲系數(shù)下，曲波循環(huán)平移去噪后的峰值信噪比最大，原始曲波變換其次，小波去噪后的最小。通過比較圖2和圖3，曲波變換循環(huán)平移去噪法效果最好，原始曲波變換其次，小波變換最差，且曲波變換循環(huán)平移去噪法改善了振鈴效應(yīng)，增強去噪效果。

2.2實際地震數(shù)據(jù)

基于理論分析結(jié)果，進一步將該算法應(yīng)用于實際地震剖面，可以看出圖4(a) 實際地震原剖面中存在很強的隨機干擾，針對剖面中隨機干擾，用3種算法進行處理，分析處理效果。

圖4　實際地震剖面試驗對比圖

在原始實際地震剖面圖4(a)的基礎(chǔ)上分別進行了小波去噪、原始曲波去噪和曲波循環(huán)平移去噪的處理，對應(yīng)的試驗處理效果如圖4中的(b)、(c)、(d)。對比圖4中的(b)、(c)、(d)可以發(fā)現(xiàn)，曲波變換處理較小波變換處理后剖面噪聲有大幅度的減少，且將曲波循環(huán)平移效果圖4(d)與原始曲波去噪效果圖4(c)對比，可看出圖4(d)更清晰，且圖中畫圈標記的部分說明圖4(d)的同相軸更連續(xù)，改善了振鈴效應(yīng)。以上圖像對比說明曲波循環(huán)平移去噪法效果最好，有效波保真度好，且改善了振鈴效應(yīng)，增強去噪效果。

3結(jié)論

通過對仿真地震數(shù)據(jù)與實際地震剖面處理效果表明，筆者所采用的曲波變換，作為一種新興的多尺度分析方法，與蒙特卡羅閾值結(jié)合去噪，可以很好的衰減隨機噪聲，實現(xiàn)了最大程度保留有效反射系數(shù)的前提下分離隨機噪聲的目的，再利用循環(huán)平移處理，消除人為干擾，增強了去噪效果。在實際工區(qū)地震數(shù)據(jù)處理過程中，其數(shù)據(jù)量大而復(fù)雜，并且算法多次利用曲波變換相關(guān)算法還涉及到多個循環(huán)嵌套，使整個算法的處理速度相對比較慢。在此基礎(chǔ)上，采用 GPU/CPU協(xié)同并行處理這一技術(shù)可以加速和優(yōu)化整個算法。

[參考文獻]

[1]Chen Longlong,Liu Yifei. Research on Main Constraints in Sustainable Development of China Oil-Gas Upstream Industry[J]. Energy Procedia,2012,14: 325～330.

[2] Xuezhi Zhao,Bangyan Ye. Similarity of signal processing effect between Hankel matrix-based SVD and wavelet transform and its mechanism analysis[J]. Mechanical Systems and Signal Processing,2008，234: 1062～1075.

[3] 莫延鋼，朱振今,朱偉強，等. F-K濾波適用公式研究[J]. 科協(xié)論壇(下半月)，2013(4)：125～126.

[4] 胡小健，張超龍. 利用Randon變換進行圖像分割與特征提取[J]. 仲愷農(nóng)業(yè)技術(shù)學院學報，2006，19(1)：44～46.

[5] Galiana-Merino J J，Rosa-Herranz J L，Rosa-Cintas S，et al. SeismicWaveTool: Continuous and discrete wavelet analysis and filtering for multichannel seismic data[J]. Computer Physics Communications，2013，184(1)：162～171.

[6] Xinhuan Deng，Haitao Ma，Yue Li，et al. Seismic random noise attenuation based on adaptive time-frequency peak filtering[J]. Journal of Applied Geophysics，2015，113(10)：31～37.

[7] 張博. 基于稀疏變換的地震數(shù)據(jù)去噪方法研究[D]. 杭州：浙江大學，2013.

[8] 繆彥舒. 基于Curvelet變換的高密度地震弱信號檢測與去噪方法研究[D]. 北京：中國石油大學，2010.

[9] Chunxue Huang，Wenxing Bao. A Remote Sensing Image Fusion Algorithm Based on the Second Generation Curvelet Transform and DS Evidence Theory[J]. Journal of the Indian Society of Remote Sensing，2014，423：645～650.

[10] Jianping Zhao，Huimin Pang，Jindong Xu. Digital Image Watermarking Algorithm Based on Fast Curvelet Transform[J]. Journal of Software Engineering and Applications，2010， 3(10)：939～943.

[11] Ziwei Chen，Shun Wang，Guangyou Fang，et al. Ionograms denoising via curvelet transform[J]. Advances in Space Research，2013， 52(7)：1289～1296.

[12] Cui Huimin，Zhao Ruimei，Hou Yanli. Improved Threshold Denoising Method Based on Wavelet Transform[J]. Physics Procedia，2015，33(1)：1354～1359.

[13] Behzad S，Rohan D，Troy L，et al. Towards elimination of the dark-rim artifact in first-pass myocardial perfusion MRI: Removing gibbs ringing effects using optimized radial imaging.[J]. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine/Society of Magnetic Resonance in Medicine， 2014，72(1)：124～136.

[14] Larissa C S R，Marcelo A C V，Homero S，et al. Mammographic Image Denoising and Enhancement Using the Anscombe Transformation， Adaptive Wiener Filtering， and the Modulation Transfer Function[J]. Journal of Digital Imaging，2013，26(2)：183～197.

[15] Bhadauria H S，Dewal M L. Analysis of effect of cycle spinning on wavelet- and curvelet-based denoising methods on brain CT images[J]. Journal of the Chinese Institute of Engineers，2014，37(7)：939～945.

[16] Cavusoglu B. Multiscale texture retrieval based on low-dimensional and rotation-invariant features of curvelet transform[J]. EURASIP Journal on Image and Video Processing， 2014，2014(22)：1～19.

[17] Moon B，Jun J Y，Lee J H，et al. Robust Image Denoising Using a Virtual Flash Image for Monte Carlo Ray Tracing[J]. Computer Graphics Forum，2013， 32(1)：139～151.

[18]張恒磊，劉天佑.Curvelet域蒙特卡羅估計的隨機噪聲衰減[J].西南石油大學學報(自然科學版)，2011，33(4)：64～68.

[編輯]計飛翔

[文獻標志碼]A

[文章編號]1673-1409(2016)13-0034-05

[中圖分類號]P631.4

[作者簡介]陳思雅(1994-)，女，碩士生，現(xiàn)主要從事地震信號處理方面的研究工作；通信作者：謝凱，pami2009@163.com。

[基金項目]中國石油創(chuàng)新基金項目(2010D-5006-0304)；湖北省大學生創(chuàng)新創(chuàng)業(yè)訓練項目(2014012)。

[收稿日期]2016-01-10

[引著格式]陳思雅，謝凱，阮寧君，等.一種曲波域蒙特卡羅閾值去噪算法[J].長江大學學報(自科版)，2016,13(13)：34～38.