空間局部加權(quán)回歸自適應(yīng)TFPF

鄧心歡， 馬海濤*， 李月， 楊寶俊

1 吉林大學(xué)通信工程學(xué)院， 長春　130012 2 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春　130026

空間局部加權(quán)回歸自適應(yīng)TFPF

鄧心歡1， 馬海濤1*， 李月1， 楊寶俊2

1 吉林大學(xué)通信工程學(xué)院， 長春130012 2 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春130026

摘要時頻峰值濾波(TFPF)算法是一種非常有效的去噪方法.但是傳統(tǒng)的TFPF采用的單一窗長，并且僅沿時間方向進行濾波，忽略了信號的空間信息，并且TFPF近似等效成一個時不變的低通濾波器，不能追蹤快速變化的信號.針對這些問題，引入空間局部加權(quán)回歸自適應(yīng)TFPF (SLWR-ATFPF).鑒于隨機噪聲在各個位置的方向隨機性，以及有效信號在各個位置的方向確定性，首先利用空間局部加權(quán)回歸(SLWR)，對含噪信號進行空間加權(quán)，從而使加權(quán)之后的信號包含空間信息.然后，再引入凸集和Viterbi的思想，對空間加權(quán)之后的信號進行自適應(yīng)濾波.從而，完成時空域二維自適應(yīng)濾波.將SLWR-ATFPF應(yīng)用于合成記錄和實際的共炮點記錄，實驗結(jié)果表明，改進的方法與原算法相比，能夠在壓制低信噪比(SNR)隨機噪聲的同時更好地保留有效信號.

關(guān)鍵詞時頻峰值濾波(TFPF)； 自適應(yīng)； 空間局部加權(quán)回歸(SLWR)； 凸集； Viterbi算法

To solve these problems, we introduce the spatial locally weighted regression adaptive TFPF (SLWR-ATFPF). Due to the randomness of the random noise at each location and the definiteness of the valid signal at each location, we first apply the spatial locally weighted regression (SLWR) to achieve the space weighting of the noisy signal, which would make the signal contain spatial information. Then, we introduce the idea of convex sets and the Viterbi′s algorithm to finish the adaptive filtering of the spatially weighted signal. Finally, we complete the spatiotemporal adaptive filtering.

To verify the feasibility and effectiveness of the proposed approach, we use synthetic and real seismic records to test it. Through the comparison of waveforms and spectrums, it is easy to see that SLWR-ATFPF could preserve valid signals, including peaks and troughs more effectively. Moreover, the results of SLWR-ATFPF are closer to the noise-free signal both on high- and low-frequency components. In the real seismic record, the results of SLWR-ATFPF recover the reflection events more continuously and smoothly. In addition, some originally buried events are revealed clearly.

SLWR-ATFPF as an improvement method of TFPF could achieve better random noise attenuation and higher continuity and clarity of reflection events than conventional TFPF and adaptive TFPF (ATFPF). Therefore, this method has large applied space.

1引言

在地震勘探資料中，有效信號往往被湮沒在強隨機噪聲中，阻礙有效信息的恢復(fù).因此，壓制低SNR地震資料中隨機噪聲并最大限度保留有效信號，是地震勘探學(xué)界研究的難點和熱點.

時頻峰值濾波(TFPF)是消減地震勘探隨機噪聲的有效方法之一(金雷等, 2005; Lin et al.,2007; 林紅波等, 2008)，具有約束條件少，對低SNR地震資料適應(yīng)性強等優(yōu)勢.為了有效壓制噪聲且恢復(fù)有效信號，則要求信號是線性的(Boashash and Mesbah, 2004).因此，在應(yīng)用TFPF處理實際地震資料時，采用加窗的Wigner-Ville分布(WVD), 即PWVD來實現(xiàn)地震信號的局部線性化(李月等, 2009).現(xiàn)有研究表明，有效的壓制隨機噪聲和信號幅度保持對窗長的要求是矛盾的，為了得到理想的去噪效果需要大窗長，但在大窗長內(nèi)信號是非線性的，從而導(dǎo)致信號幅度衰減(Tian and Li, 2014).此外，由于地震勘探資料具有非平穩(wěn)的特征，單一窗長無法同時滿足復(fù)雜地震信號的局部線性化條件，從而導(dǎo)致該方法處理含噪信號時存在有效信號損失及噪聲壓制不足的問題(Lin et al.,2014; Zhuang et al., 2015).為此，Zeng等(2015)提出了基于凸集和Viterbi算法的自適應(yīng)TFPF(ATFPF).該方法將不同窗長的TFPF濾波結(jié)果組成一個凸集，然后在這個凸集上通過Viterbi算法尋找最小化給定的全局二次泛函所在位置的濾波結(jié)果作為有效信號的最優(yōu)估計.該方法在去噪的基礎(chǔ)上，能夠快速地跟蹤信號變化，無畸變地恢復(fù)出有效信號.但此方法也存在一定的問題：首先，該方法是時域一維濾波，并沒有考慮到信號的空間信息；其次，該方法使用的全局二次泛函是能量函數(shù)，所以當(dāng)處理較低SNR的含噪信號時，不能很好地估計出有效信號.

為了解決上述問題，本文提出了空間局部加權(quán)回歸自適應(yīng)TFPF (SLWR-ATFPF)的二維濾波算法.首先，對含噪信號進行空間局部加權(quán)回歸 (SLWR)，使處理之后的信號包含空間信息(Deng et al., 2010, 2011).然后，再對包含空間信息的含噪信號進行基于凸集和Viterbi算法的自適應(yīng)TFPF，從而完成時空二維自適應(yīng)濾波.

2空間局部加權(quán)回歸自適應(yīng)TFPF

2.1空間局部加權(quán)回歸 (SLWR)

本文所考慮的地震勘探資料可以用如下方程建模：

(1)

(2)

式中，Crk為道間回歸系數(shù)，Ctk為時間回歸系數(shù).

由于地震勘探信號具有非平穩(wěn)的性質(zhì)，我們引入局部的思想.對含噪信號s進行SLWR(玄海燕等, 2013).選用的局部空間的大小為m×n(m，n均為奇數(shù)).則局部道間回歸系數(shù)crk可以通過局部加權(quán)最小二乘得到.具體的過程如下：

設(shè)對應(yīng)于采樣點t(n-1)/2上r1～rm地震道的觀測值為s(n-1)/2,j=[s(n-1)/2,1,…,s(n-1)/2,m]，進行加權(quán)最小二乘擬合，即使得

令R為Vandermonde矩陣，βr為擬合系數(shù)矩陣，Wr為加權(quán)矩陣，且需要滿足文獻(Cleveland, 1979)中的4個要求.利用加權(quán)最小二乘理論，得到的道間加權(quán)回歸值為(虞樂和肖基毅, 2012)：

(3)

則在r(n-1)/2,j處的擬合值為：

(4)

上述求解crk的過程可以通過相應(yīng)的變換進行簡化.如下所示：

由于加權(quán)矩陣Wr為對角陣，因此可以對其進行Cholesky分解，得到：

(5)

把(5)式代入(3)式中，可以得到：

S(n-1)/2,j=(SrR)[(SrR)T(SrR)]-1(SrR)Ts(n-1)/2,j,

(6)

對公式(6)的證明見附錄A.令B=SrR,Pr=B(BTB)-1BT,則得到S(n-1)/2,j=Prs(n-1)/2,j.不難發(fā)現(xiàn)，Pr為正交投影矩陣(龐善起, 2005).對矩陣B進行正交三角分解，得到：

(7)

則

Pr=B(BTB)-1BT

(8)

綜上可知，局部道間回歸系數(shù)：crk=Pr(:,(m-1)/2).

從而，完成對局部道間回歸系數(shù)求解的簡化.同理可以得到沿采樣點方向的局部時間回歸系數(shù)ctk=Pt((n-1)/2,:)，其中Pt為時間方向上的正交投影矩陣.則

E(S)=E(crksctk)

=E(crk(x+g)ctk)

=E(crkxctk)+E(crkgctk).

(9)

由于隨機噪聲在各個位置的方向隨機性，以及有效信號在各個位置的方向確定性，可以得到：

(10)

(11)

因此，經(jīng)過SLWR處理之后得到的空間加權(quán)信號中，有效信號部分得到了增強，隨機噪聲得到了衰減.2.2TFPF

TFPF是一種基于時頻分析理論的噪聲消減算法(Boashash and Mesbah, 2004).以(1)式含噪信號的模型為例，簡述一下TFPF的實現(xiàn)過程.

首先，對含噪信號s進行調(diào)頻，使其成為調(diào)頻信號的瞬時頻率(IF)，

(12)

其中，μ>0為調(diào)頻系數(shù).

然后，計算關(guān)于z的PWVD：

(13)

最后，從(W)i,f中得到峰值處的瞬時頻率作為x的估計.

(14)

簡化上述三個步驟，TFPF最終可以等效成一個低通濾波(Zeng et al., 2015)，表示成如下形式：

(15)

(16)

2.3SLWR-ATFPF

conv(H)={hc∶hc=∑αvhLv,αv≥0,∑αv=1},

(17)

顯然，conv(H)為凸集，hc為不同窗長TFPF沖激響應(yīng)的凸組合(Zeng et al., 2015).根據(jù)卷積運算的線性特點，空間加權(quán)信號S的估計集合Yc為：

(18)

也為凸集.且Yc=conv(Y), Y={yLv=hLv*S,hLv∈H}.從Yc中選擇y等效于從conv(H)中選擇hc.再根據(jù)Hilbert投影定理，可知某定義在n的距離泛函J(y)一定可以在Yc上取得全局最小值(Boyd and Vandenberghe, 2004).

令定義在Yc上的距離泛函J(y)為：

J(y)=‖A(y-S)‖2+‖D(By)‖2, l≤y≤u

(19)

J(y)的表達式中，第一項表示對觀測信號的跟蹤，第二項表示對觀測信號的平滑.A，B控制著跟蹤與平滑程度之間的平衡.為了使y能夠自適應(yīng)于信號的變化，一個合理的想法是選取參數(shù)為(Zeng et al., 2015)：

(20)

其中，A體現(xiàn)了對Yc內(nèi)信號最大衰減幅度的跟蹤，B體現(xiàn)了對噪聲最大衰減的跟蹤.hui和hli分別為上下邊界處等效成的濾波器，σ為噪聲的標(biāo)準(zhǔn)差的估計值.由于該最優(yōu)化問題往往是病態(tài)的，故選擇Viterbi算法作為求解工具(Forney, 1973).SLWR-ATFPF的具體步驟如下：

(1) 對含噪信號s進行SLWR處理，得到空間加權(quán)信號S.

(2) 給定濾波器序列H，求解空間加權(quán)信號S的估計集合Yc，以及Yc的上下邊界u，l.

(3) 在每個時間采樣點處取區(qū)間[li,ui]的分割點集Pi，使得Pi={yi,k:yi,k=li+k(ui-li)/N,k=0,1,2,…,N}，其中N為層數(shù)，以yi,k為頂點，可以組成網(wǎng)格圖1.

圖1　SLWR-ATFPF的網(wǎng)格圖Fig.1　Trellis diagram of SLWR-ATFPF

(4) 根據(jù)泛函J(y)定義相鄰頂點的距離如下表達式所示：

(21)

其中，yi,a，yi-1,b分別為Pi和Pi-1中的頂點，且其方向為Pi-1指向Pi.Si∈S，為對應(yīng)空間加權(quán)信號中，第i個采樣點處的值.

(5) 根據(jù)Viterbi算法，在網(wǎng)格圖中尋找距離最短路徑y(tǒng)，也就是使得J(y)最小的最優(yōu)解，則y即為有效信號x的最優(yōu)估計.

3仿真實驗與應(yīng)用

應(yīng)用SLWR-ATFPF處理單炮反射地震資料模型.采用Ricker子波構(gòu)建模擬地震記錄，包含兩個反射軸，主頻分別為45 Hz，30 Hz，對應(yīng)的速度為2100 m·s-1，2600 m·s-1，采樣頻率為500 Hz.加入高斯白噪聲使SNR為-5 dB，如圖2a所示.分別采用TFPF(窗長為7)，ATFPF(H={h1,h13,h53,h1023})以及本文提出的SLWR-ATFPF(H={h1,h13,h53,h1023})對含噪信號進行處理, 得到的結(jié)果如圖 2所示.對比圖2a—2d，不難發(fā)現(xiàn)圖2d軸向更連續(xù).為了使效果更加明顯，我們?nèi)〕龅?3道進行幅度和頻譜的對比，如圖3所示.通過對比不難發(fā)現(xiàn)，SLWR-ATFPF在波峰、波谷處幅度的保持，以及頻譜成分的跟蹤均優(yōu)于前兩種方法.這主要是由于TFPF僅采用單一窗長，而且沒有包含道間相關(guān)信息；ATFPF雖然在時間方向上可以得到最優(yōu)解，但是當(dāng)SNR較低時，基于能量泛函J(y)的最優(yōu)解判斷，就會把有效信號當(dāng)成噪聲去掉；而SLWR-ATFPF克服了上述問題，通過時空二維濾波，先將信號進行空間加權(quán)，再處理，從而得到較為理想的結(jié)果.表1列舉了不同方法處理之后的SNR，不難發(fā)現(xiàn)SLWR-ATFPF處理之后得到的SNR最高.

表1　信噪比定量對比(單位：dB)

將本文提出的SLWR-ATFPF應(yīng)用于168道，4000采樣點，采樣頻率為1000 Hz共炮點記錄(見圖4).其中，圖4a為原始記錄，圖4b為窗長為9的TFPF處理結(jié)果，圖4c是ATFPF(H={h1,h9,h49,h499})處理的結(jié)果，圖4d是本文提出方法(H={h1,h9,h49,h499})的處理結(jié)果.對比四幅圖片，不難看出SLWR-ATFPF處理的結(jié)果中對隨機噪聲有更好的壓制效果，為了體現(xiàn)改進算法的保幅效果，我們對方框里的結(jié)果進行了放大，對應(yīng)的放大圖如圖 5所示.通過對比不難發(fā)現(xiàn)，SLWR-ATFPF處理之后，同相軸更加連續(xù).

圖2　合成地震記錄結(jié)果比較(a) 含噪合成地震記錄; (b) TFPF處理結(jié)果；(c) ATFPF處理結(jié)果; (d) SLWR-ATFPF處理結(jié)果.Fig.2　Comparison of results of synthetic records(a) Noisy synthetic seismic records; (b) Result of TFPF; (c) Result of ATFPF; (d) Result of SLWR-ATFPF.

圖3　單道合成地震記錄結(jié)果比較(a) 23道45 Hz的波形; (b) 23道45 Hz的頻譜；(c) 23道30 Hz的波形; (d) 23道30 Hz的頻譜.Fig.3　Comparison of results of single channel synthetic records(a) Waveform of 45 Hz of 23rd channel; (b) Spectrum of 45 Hz of 23rd channel; (c) Waveform of 30 Hz of 23rd channel; (d) Spectrum of 30 Hz of 23rd channel.

圖4　共炮點記錄處理結(jié)果對比(a) 實際地震勘探記錄； (b) TFPF處理結(jié)果； (c) ATFPF處理結(jié)果； (d) SLWR-ATFPF處理結(jié)果.Fig.4　Comparison of common shot point seismic data processing results(a) Real seismic record; (b) Result of TFPF; (c) Result of ATFPF; (d) Result of SLWR-ATFPF.

圖5　局部共炮點記錄處理結(jié)果對比(a) 局部實際地震勘探記錄； (b) TFPF處理局部放大結(jié)果； (c) ATFPF處理局部放大結(jié)果； (d) SLWR-ATFPF處理局部放大結(jié)果.Fig.5　Comparison of local common shot point seismic data processing results(a) Local real seismic record; (b) Local enlarged result of TFPF; (c) Local enlarged result of ATFPF; (d) Local enlarged result of SLWR-ATFPF.

4結(jié)論

本文針對傳統(tǒng)的TFPF僅采用單一窗長且沒用利用空間相關(guān)信息的問題，研究了SLWR對TFPF的影響，并提出了SLWR-ATFPF.通過分析TFPF的特點，將SLWR處理之后的空間加權(quán)信號作為輸入，再進行ATFPF，從而完成時空二維濾波.仿真實驗和共炮點記錄處理結(jié)果表明，本文提出的算法降低了濾波后有效信號在波峰和波谷處的衰減，具有更好的保幅性，同時對噪聲有良好的壓制作用；在實際的地震資料處理中，該方法也得到了較好的處理效果.

附錄A對公式(6)的證明

References

Boashash B, Mesbah M. 2004. Signal enhancement by time-frequency peak filtering.IEEETransactionsonSignalProcessing, 52(4): 929-937. Boyd S, Vandenberghe L. 2004. Convex Optimization. Cambridge: Cambridge University Press.

Cleveland W S. 1979. Robust locally weighted regression and smoothing scatterplots.JournaloftheAmericanStatisticalAssociation, 74(368): 829-836.

Deng X Y, Yang D H, Peng J M, et al. 2010. Noise reduction and drift removal using least-squares support vector regression with the implicit bias term.Geophysics, 75(6): V119-V127.

Deng X Y, Yang D H, Liu T, et al. 2011. Study of least squares support vector regression filtering technology with a new 2D Ricker wavelet kernel.JournalofSeismicExploration, 20(2): 161-176.

Forney G D. 1973. The Viterbi algorithm.ProceedingsoftheIEEE, 61(3): 268-278.

Jin L, Li Y, Yang B J. 2005. Reduction of random noise for seismic data by time-frequency peak filtering.ProgressinGeophysics(in Chinese), 20(3): 724-728.

Lax P D. 2002. Functional Analysis. New York: John Wiley.

Li Y, Lin H B, Yang B J, et al. 2009. The influence of limited linearization of time window on TFPT under the strong noise background.ChineseJournalofGeophysics(in Chinese), 52(7): 1899-1906, doi: 10.3969/j.issn.0001-5733.2009.07.025. Lin H B, Li Y, Yang B J. 2007. Recovery of seismic events by time-frequency peak filtering. ∥ Proceedings of the IEEE International Conference on Image Processing. San Antonio, TX: IEEE, 5: V-441-V-444. Lin H B, Li Y, Ye W H, et al. 2008. Denoising technique based on time-frequency peak filtering and its application.ProgressinGeophysics(in Chinese), 23(6): 1953-1957.

Lin H B, Li Y, Yang B J, et al. 2014. Seismic random noise elimination by adaptive time-frequency peak filtering.IEEEGeoscienceandRemoteSensingLetters, 11(1): 337-341.

Pang S Q. 2005. A class of orthogonal projection matrices and related orthogonal arrays.ActaMathematicaeApplicataeSinica(in Chinese), 28(4): 668-674. Tian Y N, Li Y. 2014. Parabolic-trace time-frequency peak filtering for seismic random noise attenuation.IEEEGeoscienceandRemoteSensingLetters, 11(1): 158-162.

Xuan H Y, Li Q, Yang N N. 2013. Variable bandwidth and local estimation of geographically and temporally weighted regression model.JournalofGansuLianheUniversity(NaturalSciences) (in Chinese), 27(4): 1-4.

Yu L, Xiao J Y. 2012. Implementation of robust locally weighted regression algorithm in data mining.ComputerKnowledgeandTechnology(in Chinese), 8(7): 1493-1495.

Zeng Q, Li Y, Deng X H. 2015. Adaptive time-frequency peak filtering based on convex sets and the Viterbi algorithm.IEEEGeoscienceandRemoteSensingLetters, 12(6): 1267-1271.

Zhuang G H, Li Y, Liu Y P, et al. 2015. Varying-window-length TFPF in high-resolution radon domain for seismic random noise attenuation.IEEEGeoscienceandRemoteSensingLetters, 12(2): 404-408.

附中文參考文獻

金雷, 李月, 楊寶俊. 2005. 用時頻峰值濾波方法消減地震勘探資料中隨機噪聲的初步研究. 地球物理學(xué)進展, 20(3): 724-728.

李月, 林紅波, 楊寶俊等. 2009. 強隨機噪聲條件下時窗類型局部線性化對TFPF技術(shù)的影響. 地球物理學(xué)報, 52(7): 1899-1906, doi: 10.3969/j.issn.0001-5733.2009.07.025.

林紅波, 李月, 葉文海等. 2008. 時頻峰值濾波去噪技術(shù)及其應(yīng)用. 地球物理學(xué)進展, 23(6): 1953-1957.

龐善起. 2005. 一類正交投影矩陣及其相關(guān)正交表. 應(yīng)用數(shù)學(xué)學(xué)報, 28(4): 668-674.

玄海燕, 李琪, 楊娜娜. 2013. 時空加權(quán)回歸模型的變窗寬局部估計. 甘肅聯(lián)合大學(xué)學(xué)報: 自然科學(xué)版, 27(4): 1-4.

虞樂, 肖基毅. 2012. 數(shù)據(jù)挖掘中強局部加權(quán)回歸算法實現(xiàn). 電腦知識與技術(shù), 8(7): 1493-1495.

(本文編輯何燕)

Spatial locally weighted regression adaptive time-frequency peak filtering

DENG Xin-Huan1, MA Hai-Tao1*, LI Yue1, YANG Bao-Jun2

1CollegeofCommunicationandEngineering,JilinUniversity,Changchun130012,China2CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China

AbstractThe time-frequency peak filtering (TFPF) algorithm is an effective method for random noise attenuation. Recently, TFPF has been widely applied to seismic random noise attenuation. Whereas, TFPF still has some shortcomings. (1) The conventional TFPF uses a fixed window length for all frequency components. As a consequence, serious loss of the valid information or insufficient suppression of the noise is unavoidable. (2) TFPF carries out filtering only along the time direction, which lacks consideration of the spatial correlation. (3) TFPF is approximately equivalent to a time-invariant low-pass filer, which means that the frequency components of signal higher than some cut-off frequency would be attenuated.

KeywordsTime-frequency peak filtering (TFPF); Adaptive; Spatial locally weighted regression (SLWR); Convex sets; Viterbi algorithm

基金項目國家自然科學(xué)基金重點項目(41130421)資助.

作者簡介鄧心歡, 男, 吉林省人, 漢, 吉林大學(xué)信號與信息處理專業(yè)研究生，研究方向為地震信號處理.E-mail:dengxinhuan@yeah.net *通訊作者馬海濤，男，副教授，從事地震勘探信號處理.E-mail:mahaitao04@yeah.net

doi:10.6038/cjg20160525 中圖分類號P631

收稿日期2015-01-23，2016-03-04收修定稿

鄧心歡， 馬海濤， 李月等. 2016. 空間局部加權(quán)回歸自適應(yīng)TFPF.地球物理學(xué)報，59(5):1824-1830,doi:10.6038/cjg20160525.Deng X H, Ma H T, Li Y, et al. 2016. Spatial locally weighted regression adaptive time-frequency peak filtering.ChineseJ.Geophys. (in Chinese),59(5):1824-1830,doi:10.6038/cjg20160525.