陳建友,王曉凱,楊長(zhǎng)春

(1.西安交通大學(xué),陜西西安710049;2.酒泉衛(wèi)星發(fā)射中心,甘肅酒泉732750;3.中國(guó)科學(xué)院地質(zhì)與地球物理研究所,北京100029)

近年來(lái),壓縮感知(compressive sensing,CS)技術(shù)[1]已在地震勘探領(lǐng)域引起極大關(guān)注并取得了諸多成功應(yīng)用[2-6]。多維數(shù)據(jù)規(guī)則化、多維數(shù)據(jù)去噪及壓縮是地震數(shù)據(jù)處理中的重要步驟[7-8],其中最為重要的是尋找合適的稀疏變換,這同時(shí)也是壓縮感知的重要步驟。XU等[9]于2005年提出了反泄露傅里葉變換(anti-leakage Fourier transform,ALFT)并將其應(yīng)用于地震數(shù)據(jù)規(guī)則化,2010年又將其推廣到高維數(shù)據(jù)規(guī)則化[10]。隨后幾年,ALFT及其各種發(fā)展形式在地震數(shù)據(jù)處理中得到了廣泛應(yīng)用[11-13]。梁東輝[14]將非均勻快速傅里葉變換引入反泄露傅里葉變換,以提高反泄露傅里葉變換在實(shí)際地震數(shù)據(jù)處理中的運(yùn)算速度。

在實(shí)際應(yīng)用ALFT中,尤其是在處理高維地震數(shù)據(jù)時(shí),通常進(jìn)行分窗處理以減少計(jì)算量。雖然采用分窗處理后,窗內(nèi)的數(shù)據(jù)可以近似認(rèn)為是少量平面波疊加,但是,即使在均勻網(wǎng)格的前提下,較小的窗仍然會(huì)明顯降低傅里葉變換域的分辨率,同時(shí)造成頻譜泄露,不能利用傅里葉基函數(shù)稀疏表示信號(hào)。另外,采用分窗處理雖然能降低計(jì)算量,但仍大量采用離散傅里葉變換(discrete Fourier transform,DFT),對(duì)于TB規(guī)模的地震數(shù)據(jù),計(jì)算量仍然巨大。本文針對(duì)均勻網(wǎng)格的ALFT,給出一種頻率域快速實(shí)現(xiàn)方法,大幅減少ALFT的運(yùn)算量,提高ALFT應(yīng)用于數(shù)據(jù)壓

縮或多維去噪的計(jì)算效率。

1 方法原理

1.1 均勻網(wǎng)格下傅里葉變換的頻譜泄露現(xiàn)象

假設(shè)g[n]為一個(gè)64點(diǎn)離散合成信號(hào),采樣步長(zhǎng)為1,采樣頻率為1Hz。由于實(shí)信號(hào)可通過(guò)Hilbert變換轉(zhuǎn)換為解析信號(hào),因此這里僅采用復(fù)解析信號(hào),其表達(dá)式為:

(1)

式中:f0和f1為兩個(gè)復(fù)正弦信號(hào)的頻率。當(dāng)f0和f1均為1/64Hz的整數(shù)倍時(shí)(例如為1/16Hz和5/32Hz),其實(shí)部波形及利用快速傅里葉變換(fast Fourier transform,FFT)獲得的傅里葉系數(shù)模如圖1a 和圖1b所示,圖1b能夠真實(shí)地反映信號(hào)中的兩個(gè)頻率成分。當(dāng)f0和f1都不是1/64Hz的整數(shù)倍時(shí),例如為9/128Hz和21/128Hz時(shí),其實(shí)部波形及利用FFT獲得的傅里葉系數(shù)模如圖2a和圖2b所示。

圖1 合成信號(hào)①的實(shí)部(a)及傅里葉系數(shù)模(b)

圖2 合成信號(hào)②的實(shí)部(a)及傅里葉系數(shù)模(b)

雖然圖1a和圖2a均為正弦信號(hào),但是其振幅譜(圖1b 和圖2b)差別非常大。究其原因,圖1a所示信號(hào)的兩個(gè)分量均能夠用64個(gè)相互正交的傅里葉基函數(shù)中的一個(gè)表示(當(dāng)f0和f1為1/64Hz的整數(shù)倍時(shí)),因此可以準(zhǔn)確稀疏地表示該信號(hào);而圖2a所包含的兩個(gè)分量均不能用64個(gè)相互正交的傅里葉基函數(shù)中的一個(gè)表示,必須采用這64個(gè)相互正交的傅里葉基函數(shù)來(lái)表示圖2a所示信號(hào),造成單頻信號(hào)的能量泄露到其它頻率,即頻譜泄露的假象。

1.2 均勻網(wǎng)格下的反泄露傅里葉變換

盡管XU等[9-10]是針對(duì)數(shù)據(jù)規(guī)則化中的非均勻網(wǎng)格提出的反泄露傅里葉變換,但是該算法對(duì)于均勻網(wǎng)格同樣適用。

假設(shè)g[n]為待分析信號(hào),G[f]為利用DFT計(jì)算的g[n]傅里葉系數(shù):

(2)

2) 利用DFT計(jì)算r[n]的傅里葉系數(shù)R[f];

3) 找出R[f]中的最大模值及其對(duì)應(yīng)的fk,并記錄fk和對(duì)應(yīng)的R[fk],

(3)

(4)

5) 更新r[n],

(5)

6) 如果r[n]的能量小于設(shè)定的能量閾值ε,則退出迭代過(guò)程,否則返回步驟2)。

通過(guò)上述迭代,由各個(gè)單頻信號(hào)相加的信號(hào)可獲得稀疏表示(高維情況下通常進(jìn)行分窗處理,在分析窗內(nèi)可以近似認(rèn)為僅有少數(shù)幾個(gè)平面波混合在一起),并可恢復(fù)出任意其它網(wǎng)格點(diǎn)的信號(hào)。

如果一維信號(hào)包含有N個(gè)采樣點(diǎn),頻率域采用4N個(gè)采樣點(diǎn),則迭代一次的運(yùn)算復(fù)雜度為O[4N·lg(4N)]。假設(shè)共迭代M次,則上述過(guò)程的運(yùn)算復(fù)雜度大約為O[4MNlg(4N)]。

1.3 頻率域均勻網(wǎng)格ALFT的快速實(shí)現(xiàn)方法

觀察上述ALFT實(shí)現(xiàn)過(guò)程可知,步驟2)為運(yùn)算量最大的步驟。如果能將該步驟的運(yùn)算量減少,則整個(gè)算法的運(yùn)算量會(huì)大幅降低。上述過(guò)程中之所以需要步驟2),是因?yàn)槊看蔚?信號(hào)r[n]都有所變化,導(dǎo)致需要利用DFT重新計(jì)算R[f]。如果R[f]在每一次迭代之前是已知的,則該步驟可以省去,從而大幅減少計(jì)算量。

對(duì)公式(5)做DFT,得:

(6)

其中,

(7)

為SINC函數(shù)。因此,我們可以在頻率域快速實(shí)現(xiàn)均勻網(wǎng)格的ALFT,其步驟簡(jiǎn)述如下:

2) 利用DFT計(jì)算g[n]的傅里葉系數(shù)G[f],并將R[f]初始化為G[f];

3) 尋找R[f]中模值最大的一個(gè)及其對(duì)應(yīng)的fk,并記錄fk和對(duì)應(yīng)的R[fk],

(8)

(9)

5) 通過(guò)減去SINC函數(shù)更新R[f],

(10)

6) 如果R[f]的能量小于設(shè)定閾值ε,則退出迭代過(guò)程,否則返回步驟3)。

如果一維信號(hào)包含有N個(gè)采樣點(diǎn),頻率域采用4N個(gè)采樣點(diǎn),由于上述實(shí)現(xiàn)方法避免了每次迭代中使用DFT,因此迭代一次的運(yùn)算復(fù)雜度為O(4N),迭代M次的運(yùn)算復(fù)雜度大約為O(4MN)。與前述常規(guī)實(shí)現(xiàn)方法相比,本文給出的實(shí)現(xiàn)方法運(yùn)算量大幅下降。

2 應(yīng)用算例

利用本文ALFT頻率域高效實(shí)現(xiàn)方法計(jì)算了圖2a 所示合成信號(hào)②的傅里葉系數(shù)模,結(jié)果如圖3所示。與圖2b相比,圖3明顯改善了頻譜泄露現(xiàn)象,同時(shí)對(duì)原始信號(hào)的表示更為稀疏。采用35Hz的Ricker子波合成了一個(gè)三維數(shù)據(jù)體,該數(shù)據(jù)體包含有8個(gè)傳播方向不同的平面波,時(shí)間方向采樣間隔為0.5ms,共有512個(gè)采樣點(diǎn);x方向和y方向采樣間隔均為5m,每個(gè)方向均有64個(gè)采樣點(diǎn)。整個(gè)數(shù)據(jù)體共有2097152個(gè)采樣點(diǎn)。在能量閾值ε為2.0×10-5的限制下,利用本文ALFT高效實(shí)現(xiàn)方法對(duì)該三維數(shù)據(jù)體進(jìn)行稀疏表示,僅需要2616個(gè)系數(shù)(為原始數(shù)據(jù)采樣點(diǎn)數(shù)的1/800)就可以基本恢復(fù)出三維合成數(shù)據(jù)。在同樣的誤差限制下,采用DFT則需要

53796個(gè)系數(shù)才能基本恢復(fù)出三維合成數(shù)據(jù),所需系數(shù)的個(gè)數(shù)大約為本文ALFT方法的20倍。圖4和圖5對(duì)比了利用本文ALFT方法和DFT方法壓縮三維合成數(shù)據(jù)并進(jìn)行恢復(fù)的誤差。其中圖4a與圖5a

圖3 利用本文ALFT方法計(jì)算的合成信號(hào)②的傅里葉系數(shù)模

圖4 利用本文ALFT方法和DFT方法壓縮三維合成數(shù)據(jù)恢復(fù)的誤差比較①a 原始剖面①; b 利用2616個(gè)ALFT系數(shù)恢復(fù)的誤差剖面; c 利用53796個(gè)DFT系數(shù)恢復(fù)的誤差剖面

圖5 利用本文ALFT方法和DFT方法壓縮三維合成數(shù)據(jù)恢復(fù)的誤差比較②a 原始剖面②; b 利用2616個(gè)ALFT系數(shù)恢復(fù)的誤差剖面; c 利用53796個(gè)DFT系數(shù)恢復(fù)的誤差剖面

為三維合成數(shù)據(jù)的兩個(gè)原始剖面,圖4b和圖5b為利用本文ALFT方法保留2616個(gè)系數(shù)恢復(fù)的兩個(gè)剖面的誤差;圖4c和圖5c為利用DFT方法保留53796個(gè)系數(shù)恢復(fù)的兩個(gè)剖面的誤差。可以看出,兩種方法恢復(fù)的誤差基本在一個(gè)量級(jí)上,但是采用本文ALFT方法所需要的系數(shù)僅為DFT方法的1/20。

3 結(jié)束語(yǔ)

本文提出了一種均勻網(wǎng)格反泄露傅里葉變換的頻率域高效實(shí)現(xiàn)方法。該方法直接在頻率域?qū)崿F(xiàn)迭代搜索過(guò)程,大大提高了ALFT的計(jì)算效率。合成數(shù)據(jù)應(yīng)用算例表明,利用本文ALFT高效實(shí)現(xiàn)方法可以明顯減少頻譜泄露現(xiàn)象,提高稀疏表示性能。

盡管本文是以一維數(shù)據(jù)為例推導(dǎo)了ALFT的頻率域高效實(shí)現(xiàn)方法,但是要將其推廣到高維并非難事。下一步的研究工作是針對(duì)五維數(shù)據(jù)去噪及數(shù)據(jù)壓縮進(jìn)行測(cè)試,并在以后的工作中不斷發(fā)現(xiàn)問(wèn)題。

致謝:感謝西安交通大學(xué)陳文超教授提供的建議及幫助。

[1] CANDES E J,WAKINM B.An introduction to compressive sampling[J].IEEE Signal Processing Magazine,2008,25(2):21-30

[2] 李翔.基于壓縮感知技術(shù)的全波形反演[J].石油物探,2017,56(1):20-25

LI X.Full-waveform inversion from compressively recovered updates[J].Geophysical Prospecting for Petroleum,2017,56(1):20-25

[3] 王華忠,馮波,王雄文,等.壓縮感知及其在地震勘探中的應(yīng)用[J].石油物探,2016,55(4):467-474

WANG H Z,FENG B,WANG X W,et al.Compressed sensing and its application in seismic exploration[J].Geophysical Prospecting for Petroleum,2016,55(4):467-474

[4] 陳祖慶,王靜波.基于壓縮感知的稀疏脈沖反射系數(shù)譜反演方法研究[J].石油物探,2015,54(4):459-466

CHEN Z Q,WANG J B.A spectral inversion method of sparse-spike reflection coefficients based on compressed sensing[J].Geophysical Prospecting for Petroleum,2015,54(4):459-466

[5] 劉偉,曹思遠(yuǎn),崔震.基于壓縮感知和TV準(zhǔn)則約束的地震資料去噪[J].石油物探,2015,54(2):180-187

LIU W,CAO S Y,CUI Z.Random noise attenuation based on compressive sensing and TV rule[J].Geophysical Prospecting for Petroleum,2015,54(2):180-187

[6] 周松,呂堯,呂公河,等.基于壓縮感知的非規(guī)則地震勘探觀測(cè)系統(tǒng)設(shè)計(jì)與數(shù)據(jù)重建[J].石油物探,2017,56(5):617-625

ZHOU S,LV Y,LV G H,et al.Irregular seismic geometry design and data reconstruction based on compressive sensing[J].Geophysical Prospecting for Petroleum,2017,56(5):617-625

[7] 張良,韓立國(guó),劉爭(zhēng)光,等.基于壓縮感知和Contourlet變換的地震數(shù)據(jù)重建方法[J].石油物探,2017,56(6):804-811

ZHANG L,HAN L G,LIU Z G,et al.Seismic data reconstruction based on compressed sensing and Contourlet transform[J].Geophysical Prospecting for Petroleum,2017,56(6):804-811

[8] 馮飛,王征,劉成明,等.基于Shearlet變換的稀疏約束地震數(shù)據(jù)重建[J].石油物探,2016,55(5):682-691

FENG F,WANG Z,LIU C M,et al.Seismic data reconstruction based on sparse constraint in the Shearlet domain[J].Geophysical Prospecting for Petroleum,2016,55(5):682-691

[9] XU S,ZHANG Y,DON P,et al.Antileakage Fourier transform for seismic data regularization[J].Geophysics,2005,70(4):87-95

[10] XU S,ZHANG Y,GILLES L.Antileakage Fourier transform for seismic data regularization in higher dimensions[J].Geophysics,2010,75(6):113-120

[11] MICHEL S,ANDREAS K,ALAN V.Anti-alias anti-leakage Fourier transform[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:3249-3252

[12] MOSTAFA N,INNANEN K A.Two-dimensional fast generalized Fourier interpolation of seismic records[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-4

[13] MICHEL S,YAN Z M,MARTIN B,et al.Matching pursuit Fourier interpolation using priors derived from a second data set[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:3651-3654

[14] 梁東輝.基于傅里葉變換的地震數(shù)據(jù)規(guī)則化和插值[D].杭州:浙江大學(xué),2015

LIANG D H.Seismic data normalization and interpolation based on Fourier transform[D].Hangzhou:Zhejiang University,2015