齊少華，劉啟元，陳九輝，郭飚

?

接收函數(shù)的曲波變換去噪

齊少華，劉啟元*，陳九輝，郭飚

中國地震局地質(zhì)研究所地震動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京100029

摘要壓制橫向非均勻地殼介質(zhì)引起的散射波場對于基于水平分層介質(zhì)模型的接收函數(shù)地殼結(jié)構(gòu)成像及其地震各向異性研究至關(guān)重要.雖然通過數(shù)據(jù)疊加和低通濾波，在一定程度上能夠壓制散射波場，但也有可能導(dǎo)致不希望的波形畸變、信息丟失或數(shù)據(jù)分辨率降低.為了避免上述問題，本文將近年來快速發(fā)展的曲波變換理論用于遠(yuǎn)震接收函數(shù)的散射噪聲壓制.與勘探地震學(xué)不同，我們面臨的主要問題在于地震臺站和震源的空間缺失導(dǎo)致的接收函數(shù)空間不均勻采樣.為此，我們將壓縮感知理論與曲波變換去噪相結(jié)合，在對缺失數(shù)據(jù)進(jìn)行波場重建的同時，實(shí)現(xiàn)散射噪聲的壓制.為了論證方法的可行性，本文進(jìn)行了噪聲壓制和波場重建的理論檢驗(yàn).并將本文方法用于處理IRIS全球臺網(wǎng)固定臺站和川西臺陣遠(yuǎn)震接收函數(shù).結(jié)果表明： 1)地殼介質(zhì)橫向不均勻引起的散射噪聲可以得到有效壓制，接收函數(shù)的信噪比得到提高，震相的可追蹤性得到改善，從而利于進(jìn)一步的接收函數(shù)反演和地震各向異性研究； 2)缺失數(shù)據(jù)可以正確重建； 3)本文的方法既可用于單臺-多事件的數(shù)據(jù)集，也可用于單個事件-陣列觀測的數(shù)據(jù)去噪，但單臺-多事件數(shù)據(jù)集的結(jié)果優(yōu)于陣列觀測的情況.

關(guān)鍵詞接收函數(shù)； 曲波變換； 壓縮感知； 橫向非均勻散射； 地殼各向異性

1引言

接收函數(shù)方法(Langston，1979)已成為地殼上地幔地震成像的一個重要方法.但是，迄今為止，接收函數(shù)的解釋仍主要基于水平分層模型(Owens et al.，1984；Duecker and Sheehan，1997；Zhu and Kanamori，2000).例如，對于接收函數(shù)反演來說，為了適應(yīng)一維水平分層模型的限制，人們不得不通過不同方位接收函數(shù)的疊加來壓制數(shù)據(jù)中橫向非均勻地殼介質(zhì)引發(fā)的散射(Liu et al.，2014).另一方面，近年來，各向異性接收函數(shù)的研究日益受到人們的重視(田寶峰等，2008；房立華和吳建平，2009；Schulte-Pelkum and Mahan，2014)，但地殼橫向非均勻介質(zhì)引起的散射給接收函數(shù)各向異性參數(shù)提取帶來了困難，而低通濾波將會降低各向異性參數(shù)的分辨能力.這并非是人們所希望的.

因此，對基于地殼水平分層模型的接收函數(shù)研究來說，在提取接收函數(shù)的過程中，環(huán)境噪聲已得到很大程度的壓制，而影響信息提取的主要噪聲來自于介質(zhì)橫向非均勻引起的散射.與環(huán)境噪聲不同，地殼橫向非均勻介質(zhì)引起的散射屬于自激式的次生波場，它通常與試圖分析研究的信號具有相近的頻帶，且在空間取向上具有較大的隨機(jī)性.但是，理論合成數(shù)據(jù)的檢驗(yàn)表明，各向異性介質(zhì)接收函數(shù)波場在空間上具有很好的連續(xù)性，或者稱為可追蹤性(Nagaya et al.，2008；房立華和吳建平，2009).本文將基于這一認(rèn)識，研究壓制接收函數(shù)中的橫向非均勻散射方法，這對推進(jìn)接收函數(shù)研究將具有重要的實(shí)際意義.

曲波變換(curvelet transform)是近年來迅速發(fā)展的多尺度、多方向高維信號分析工具(Candès and Donoho，2000，2004).相較于經(jīng)典的二維數(shù)據(jù)去噪方法，基于曲波變換的去噪，對信號的損傷最小(張軍華等，2005).在方向性特征識別上，曲波變換優(yōu)于小波變換(wavelet transform).它尤其適用于地震同相軸檢測與去噪.迄今為止，曲波變換已經(jīng)被廣泛地應(yīng)用于地震勘探數(shù)據(jù)的多次波去除(Herrmann and Verschuur，2004；張素芳等，2006)，噪聲壓制(Neelamani et al.，2008；彭才等，2008)，偏移成像(Douma and De Hoop，2007)以及地震環(huán)境噪聲的格林函數(shù)提取(Stehly et al.，2011，2015).曲波變換所具有的特性為我們壓制接收函數(shù)中橫向非均勻地殼介質(zhì)引起的散射提供了新思路.

但是，直接將曲波變換應(yīng)用于接收函數(shù)是困難的，因?yàn)榍ㄗ儞Q要求數(shù)據(jù)在各個維度上必須滿足等間距采樣.一方面受制于地震全球分布的條帶性，另一方面受具體環(huán)境條件對地震臺站(特別是流動臺站)布設(shè)的限制，接收函數(shù)在空間方向上的采樣通常是非等間距的，甚至可能存在大范圍(例如，在方位空間)數(shù)據(jù)缺失.因此，對于曲波變換去噪處理，我們需要對接收函數(shù)波場進(jìn)行正則化(均勻空間采樣).利用距離加權(quán)疊加的方法，Chai等(2015)實(shí)現(xiàn)了對空間分布相對密集且接近均勻的接收函數(shù)波場的正則化.Zhang和Zheng(2015)提出用三次樣條插值提高稀疏臺站接收函數(shù)地震成像的分辨率.作為一種探索，本文將引入由Candès等(2006b)和Donoho(2006)提出的壓縮感知(compressive sensing)理論，對接收函數(shù)波場進(jìn)行重建，進(jìn)而形成可用于曲波變換處理的數(shù)據(jù)，實(shí)現(xiàn)接收函數(shù)數(shù)據(jù)集的二維去噪.為便于理解，本文的討論將從簡要介紹曲波變換和壓縮感知的基本概念開始.

2曲波閾值法去噪

曲波可以認(rèn)為是一類具有方向性的特殊小波，它提供對曲線型邊緣的最優(yōu)稀疏表示(Candès and Donoho，2004).曲波在空間域和頻率域內(nèi)均滿足拋物尺度關(guān)系，即寬度約等于長度的平方.在空間域內(nèi)，曲波的能量在其垂向上呈現(xiàn)振蕩衰減，形態(tài)接近平面波.在頻率域內(nèi)，曲波的支撐區(qū)域是一個圍繞原點(diǎn)的楔形區(qū)域，它屬于緊支撐.為了便于直觀理解，我們利用Curvelab工具包繪制了曲波在空間域和頻率域內(nèi)的對應(yīng)關(guān)系，如圖1所示.需要說明的是，1)為了便于觀察，我們對空間域內(nèi)的曲波進(jìn)行了平移，同時，因?yàn)槎S信號的頻譜關(guān)于原點(diǎn)對稱，我們僅繪制了一半的支撐域；2)曲波在空間域和頻率域的相互正交特征取決于Fourier變換的時頻縮放特性.

圖1　曲波在空間域和頻率域的對應(yīng)關(guān)系(a) 空間域的4個尺度、方向和位置不同的曲波； (b) 頻率域中相應(yīng)曲波的支撐區(qū)域.(a)和(b)中曲波的數(shù)字編號相互對應(yīng).

根據(jù)Candès等(2006a)，曲波函數(shù)定義為

(1)

如前所述，曲波變換基為緊框架，其重構(gòu)公式具有以下形式：

(2)

其中，〈f,φj,l,k〉為信號f與曲波函數(shù)φj,l,k的內(nèi)積，c(j,k,l)為曲波系數(shù).這意味著，與信號相關(guān)的曲波系數(shù)較少且數(shù)值較大，而隨機(jī)噪聲則對應(yīng)較多的曲波系數(shù)，但數(shù)值較小.在曲波域中，信號與噪聲具有最小限度的重疊.這樣，通過閾值設(shè)定就可以對曲波系數(shù)進(jìn)行收縮和重構(gòu)，從而實(shí)現(xiàn)信號與噪聲的分離.

為此，假定由觀測數(shù)據(jù)得到的二維(空間)接收函數(shù)d符合線性矩陣方程:

(3)

這里，接收函數(shù)矩陣的行為時間序列，列為空間方位，m為不含噪聲的接收函數(shù)矩陣，n為均值為零且標(biāo)準(zhǔn)差為σ的高斯白噪聲.曲波閾值法去噪過程可以表示為

(4)

其中，C代表曲波正變換(分解矩陣)，C?代表曲波廣義逆變換(合成矩陣)，Sw為閾值函數(shù).其緊框架特性保證C?=CH，H表示共軛轉(zhuǎn)置，即曲波的廣義逆是其正變換的共軛轉(zhuǎn)置.為達(dá)到去噪的目的，我們采用所謂“硬閾值”，即把小于閾值的曲波系數(shù)置零，而閾值的選取依據(jù)對噪聲水平的估計.

3壓縮感知

傳統(tǒng)的信號采樣需要遵循Nyquist-Shannon采樣定理，即有限帶寬的信號采樣頻率必須達(dá)到其帶寬的2倍以上才能夠被精確重建.但是，近年來出現(xiàn)的壓縮感知理論突破了經(jīng)典采樣定律的限制.壓縮感知理論提出后迅速受到地震勘探領(lǐng)域的關(guān)注(HerrmannandHennenfent，2008；NaghizadehandSacchi，2010).在天然地震領(lǐng)域，壓縮感知理論也已被用于大震的震源破裂模式研究(Yaoetal.，2011，2013)和天然地震波場的逆時偏移(Shang，2014).

壓縮感知理論認(rèn)為，當(dāng)信號是稀疏的或在某個變換域內(nèi)是稀疏的，利用一個與變換基不相關(guān)的觀測矩陣(或稱為感知矩陣)可將高維信號投影到低維空間，并通過求解一個優(yōu)化問題，以高概率實(shí)現(xiàn)由少量投影重建原始信號(CandèsandWakin，2008).這意味著從遠(yuǎn)低于Nyquist-Shannon采樣頻率的樣本中也可以重建完整信號.實(shí)際上，稀疏性條件意味著信號所攜帶的信息存在某種冗余度，即信號的自由度遠(yuǎn)小于信號的長度.如果能夠通過選擇合適的變換基，使信號變換系數(shù)大部分為零或小到足以忽略，則可對該信號進(jìn)行數(shù)據(jù)壓縮.

原則上，低于Nyquist頻率的采樣操作將導(dǎo)致信號失真，即等間距的欠采樣會引發(fā)頻譜的重疊效應(yīng)，而隨機(jī)欠采樣則等效于頻率域內(nèi)的隨機(jī)噪聲.當(dāng)觀測矩陣與變換基不相關(guān)時，采樣不完整所引發(fā)的噪聲與原始信號在變換域內(nèi)可被分離.一般地，隨機(jī)欠采樣與常用的變換基不相關(guān).這樣，信號重建問題轉(zhuǎn)變?yōu)轭l率域內(nèi)的去噪問題.

需要關(guān)注的是，觀測矩陣相應(yīng)于對連續(xù)信號的采樣操作.接收函數(shù)在時間上的采樣滿足采樣定理，但在空間上采樣是不均勻的.若我們用d表示實(shí)際觀測的接收函數(shù)矩陣，且矩陣的行代表時間，列代表空間位置，則有

(5)

這里，m是插值后的接收函數(shù)矩陣，G為觀測矩陣，n為噪聲.

根據(jù)Naghizadeh和Sacchi(2010)，我們選擇曲波作為變換基，則有

m=CHWc,

(6)

這里，CH為曲波變換的共軛轉(zhuǎn)置，c為待求曲波系數(shù)，W為mask函數(shù).W的取值依據(jù)曲波閾值去噪中的閾值來設(shè)定.依據(jù)(6)式，曲波閾值去噪可以與接收函數(shù)波場的重建結(jié)合起來，在進(jìn)行波場重建過程中即可同時完成噪聲壓制，從而獲得去噪后的接收函數(shù).

若令A(yù)=GCHW，待求曲波系數(shù)c可以通過最小化目標(biāo)函數(shù)

(7)

求得.這里，λ是拉格朗日乘子，它反映了數(shù)據(jù)殘差與曲波系數(shù)稀疏解之間的折中關(guān)系.在求得系數(shù)c后，利用(6)式即可重建接收函數(shù)矩陣m.本文將利用SPGL1工具包(Van den Berg and Friedlander，2008)所提供的譜投影梯度法計算(7)式的L1范數(shù)最優(yōu)解.

4數(shù)值檢驗(yàn)

為檢驗(yàn)方法的可行性，我們首先考慮單個臺站記錄的不同方位接收函數(shù).為此，利用作者獨(dú)立開發(fā)的程序，我們計算了合成的各向異性接收函數(shù)，并將它作為數(shù)值檢驗(yàn)的依據(jù).需要說明的是，我們依據(jù)的算法與Levin和Park(1997)的方法沒有原則上區(qū)別.

計算各向異性理論接收函數(shù)所依據(jù)的地殼模型參數(shù)如表1所示.由表1可知，為簡單起見，我們僅考慮平面波入射單層各向異性地殼的情況.計算中采用的水平慢度為0.0618 s·km-1.大體上，這相當(dāng)于震中距為60°的情況.對于合成的接收函數(shù)，我們加入高斯白噪聲，用來模擬包含噪聲的觀測數(shù)據(jù).加入噪聲后的接收函數(shù)信噪比為0 dB.

表1　數(shù)值檢驗(yàn)所用的各向異性模型

注：α,β,ρ分別為P波速度、S波速度和密度.B和E分別為P波和S波各向異性強(qiáng)度，θ和φ分別為各向異性的傾角和方位角.

4.1曲波閾值去噪

我們考慮接收函數(shù)在空間方位上均勻采樣的情況.圖2給出了相應(yīng)的理論接收函數(shù)切向分量及其隨反方位角(back azimuth，0°～360°，間隔為1°)的變化.圖2表明，在時間―空間域內(nèi)，接收函數(shù)切向分量的同相軸有很好的可追蹤性(圖2a)，而在相應(yīng)的頻率域內(nèi)，這表現(xiàn)為主要能量僅集中在狹窄范圍內(nèi)(圖2b)，這意味著，若觀測數(shù)據(jù)中存在橫向非均勻散射，經(jīng)曲波變換就可以簡單地通過閾值去噪.

圖3給出了對接收函數(shù)(圖2a)加入隨機(jī)白噪聲后(圖3a)利用本文閾值去噪方法得到的結(jié)果(圖3b).在進(jìn)行曲波閾值去噪時，我們僅保留了大于閾值(閾值函數(shù)在最高尺度取值為4σ，而在其他尺度取值為3σ，其中σ為加入噪聲的標(biāo)準(zhǔn)差)且接近垂直方向的曲波系數(shù)(Stehly et al.，2011).比較圖2a和圖3b可知，經(jīng)過曲波閾值去噪后，原有信號可以正確恢復(fù).

4.2波場重建

如前所述，與通常的人工地震勘探不同，接收函數(shù)沿不同方位的采樣受控于震源的分布.這意味著 我們很難將曲波閾值去噪技術(shù)直接用于觀測數(shù)據(jù).為此，我們通過壓縮感知對進(jìn)行接收函數(shù)波場進(jìn)行重建，以便得到在空間上均勻采樣的接收函數(shù)波場.為了模擬這種情況，如圖4a所示，通過對圖2a的隨機(jī)抽樣，我們構(gòu)筑了方位非均勻分布的接收函數(shù)波場.圖4b給出了利用本文方法重建的接收函數(shù)波場(切向分量).重建波場的信噪比為36.3 dB，表明接收函數(shù)的振幅、極性和相位均可以正確恢復(fù)(Herrmann and Hennenfent，2008).

圖3　接收函數(shù)(切向分量)的曲波閾值去噪(a) 加入隨機(jī)噪聲后的接收函數(shù)； (b) 去噪后的接收函數(shù).

5觀測數(shù)據(jù)的檢驗(yàn)

對于單臺多事件觀測來說，數(shù)據(jù)中缺失部分的重建精度主要取決于地震事件的分布.通常，固定臺站由于記錄時間長，較容易獲得方位分布密集的地震事件.因此，本文首先選取IRIS全球地震臺網(wǎng)在澳大利亞東北部地區(qū)的固定臺CTAO(-20.0882°N，146.2545°E,臺網(wǎng)編號IU).它所處的地理位置有 利于獲取方位較完整的地震事件覆蓋.所用的地震為該臺站記錄的151個遠(yuǎn)震事件(自1999年2月到2015年7月，震中距30°≤Δ≤90°，Mb>5.5)，接收函數(shù)的提取采用了基于單臺的時間域迭代反褶積方法(Ligorría and Ammon，1999).

圖4　方位非均勻分布的合成接收函數(shù)(切向分量)及其波場重建(a) 由隨機(jī)抽樣得到的方位非均勻分布的接收函數(shù)波場； (b) 重建的接收函數(shù)波場.其反方位角的間隔為1°.

圖5　臺站CTAO不同方位接收函數(shù)的二維Fourier譜(a) 徑向分量(R)； (b) 切向分量(T).

圖5給出了臺站CTAO的151個遠(yuǎn)震接收函數(shù)的二維時間-方位Fourier譜.圖5顯示了分布在整個空間方位的較弱噪聲.比較圖5b和圖2b可知，與橫向均勻分層模型不同，由于實(shí)際接收區(qū)地殼介質(zhì)的橫向非均勻性，接收函數(shù)的切向分量存在著不可忽視的橫向非均勻散射.圖5a表明，這種橫向非均勻散射同樣存在于接收函數(shù)的徑向分量，盡管其散射強(qiáng)度相對較弱，它應(yīng)是非均勻空間采樣導(dǎo)致的結(jié)果.

為了客觀評價波場重建在實(shí)際數(shù)據(jù)應(yīng)用中的效果，我們首先將提取的151組接收函數(shù)(徑向與切向分量)按照反方位角1°的間隔疊加獲得123組接收函數(shù)，疊加后的事件分布如圖6c所示.我們以類似jitter欠采樣的方式(Hennenfent and Herrmann，2008)從123組接收函數(shù)中隨機(jī)抽取出35組來模擬方位缺失的數(shù)據(jù)(如圖6d中藍(lán)色和紅色圓圈所示).利用剩余的88組接收函數(shù)，以反方位角1°為間隔，進(jìn)行波場重建和去噪.最后，我們從重建后的接收函數(shù)中選擇出20組數(shù)據(jù)(如圖6d中藍(lán)色圓圈所示)，并與實(shí)際觀測接收函數(shù)進(jìn)行波形對比(如圖6a和6b所示).

圖6　臺站CTAO的觀測數(shù)據(jù)與重建數(shù)據(jù)的波形對比(a) 徑向分量的波形對比； (b) 切向分量的波形對比，其中，黑色實(shí)線代表實(shí)際觀測數(shù)據(jù)，紅色實(shí)線代表波場重建數(shù)據(jù)；(c) 按反方位角1°間隔疊加后的臺站CTAO的遠(yuǎn)震事件分布，紅色圓圈代表遠(yuǎn)震事件，綠色五角星代表臺站位置； (d) 從(c)中隨機(jī)抽取出的遠(yuǎn)震事件分布，藍(lán)色圓圈代表(a)和(b)中進(jìn)行對比的遠(yuǎn)震事件.

圖7　固定臺站CTAO實(shí)際接收函數(shù)的曲波去噪(a) 去噪前的徑向分量； (b) 去噪前的切向分量； (c) 去噪后的徑向分量； (d) 去噪后的切向分量.各方位的接收函數(shù)為以10°為間隔進(jìn)行方位分組疊加的結(jié)果.‘sum’表示全方位接收函數(shù)的疊加.

圖8　按實(shí)際地震方位抽樣的理論接收函數(shù)波場重建(a) 臺站S05記錄的地震事件分布； (b) 按實(shí)際地震方位采樣的接收函數(shù)波場(切向分量)； (c) 重建的接收函數(shù)波場(切向分量)； (d)重建合成接收函數(shù)的波場誤差.SB：四川盆地； SG：松潘地塊； CD：川滇地塊； LMS：龍門山斷裂帶； XSH：鮮水河斷裂帶； 黃色五角星：臺陣中心位置； 紅色三角：觀測臺站； 紅色圓圈：地震事件.

圖9　臺站S05實(shí)際接收函數(shù)的曲波去噪(a) 去噪前的徑向分量；(b) 去噪前的切向分量；(c) 去噪后的徑向分量；(d) 去噪后的切向分量.各方位的接收函數(shù)為以10°為間隔進(jìn)行方位分組疊加的結(jié)果.‘sum’表示全方位接收函數(shù)的疊加.

圖10　川西臺陣31°N測線19個臺站接收函數(shù)曲波閾值去噪(a) 原始的徑向分量； (b) 去噪后的徑向分量； (c) 原始的切向分量； (d) 去噪后的切向分量.圖的左側(cè)標(biāo)出了臺站代碼.

觀察可知，波場重建很好地恢復(fù)了缺失部分的數(shù)據(jù).與此同時，波場重建校正了觀測數(shù)據(jù)中的一些振幅和相位的異常，重建后的接收函數(shù)隨反方位角振幅和相位的變化更為連續(xù)和合理，這是因?yàn)榍ㄩ撝等ピ肱c波場重建是同時進(jìn)行的.因此，對于單臺多事件觀測來說，本文提出的曲波變換去噪方法不僅能夠有效壓制接收函數(shù)中的隨機(jī)散射噪聲，而且在地震事件分布較好的情況下或地震事件相對密集的一些區(qū)域內(nèi)，可以對方位缺失的數(shù)據(jù)進(jìn)行重構(gòu).

圖7給出了CTAO臺站123組接收函數(shù)的觀測數(shù)據(jù)和曲波去噪后數(shù)據(jù)按反方位角10°的間隔分組疊加的結(jié)果.對圖7我們可有以下觀察：1)利用本文給出的曲波閾值去噪方法，我們不但得到了完整(方位)的接收函數(shù)波場，而且由于橫向散射波場得到壓制，經(jīng)過曲波去噪的接收函數(shù)(徑向和切向)相較于原始數(shù)據(jù)顯示了更為清晰、合理的震相追蹤效果，表明去噪并未削弱有效信息的能量；2)對比去噪處理前后全方位接收函數(shù)疊加的結(jié)果(sum)，它們的波形沒有明顯的不同；3)對方位缺失部分預(yù)測的接收函數(shù)振幅和相位的追蹤與其他部分是協(xié)調(diào)的，這為后續(xù)研究提供更多可用數(shù)據(jù).

相對于固定臺站，流動地震臺站由于觀測時間較短(除地理位置的因素之外)，通常很難獲得方位覆蓋密集的事件分布.同時，對于中國大陸布設(shè)的臺站來說，往往存在如圖8a所示的地震事件在西北方向的嚴(yán)重缺失(原因是地中海地震帶的地震活動性相對較弱，Mb>5.0地震事件較少).為此，需要探討方位嚴(yán)重缺失對波場重建帶來的影響.

實(shí)際上，實(shí)際數(shù)據(jù)方位的缺失導(dǎo)致對實(shí)際數(shù)據(jù)按反方位角的采樣難以滿足隨機(jī)性.為此，作為方法檢驗(yàn)，我們按照實(shí)際數(shù)據(jù)的方位分布對理論接收函數(shù)波場進(jìn)行抽樣并重建.我們選用了川西臺陣中的一個臺站(S05)，圖8a給出了它的地理位置及相應(yīng)的震中分布.有關(guān)川西臺陣和臺站詳細(xì)的信息可參考相關(guān)文章(劉啟元等，2008).

圖8a表明，臺站S05的接收函數(shù)方位除了空間上的不均勻分布，在西北方向反方位角200°～300°的區(qū)間還存在明顯的方位缺失.由于其方位缺失超過了曲波長度，它導(dǎo)致缺失的數(shù)據(jù)難以精確重建(Naghizadeh and Sacchi，2010).作為對壓縮感知理論在方位缺失情況下波場重建能力的檢驗(yàn)，圖8d給出了重建的理論接收函數(shù)與原有的合成接收函數(shù)波場(圖2a)的差異.圖8d表明，利用本文方法可以正確重建方位密集區(qū)的波場，但對于方位嚴(yán)重缺失的區(qū)域，重建的接收函數(shù)存在大約6%～8%的振幅誤差，并導(dǎo)致重建波場的信噪比下降(10.4 dB).

為了進(jìn)一步檢驗(yàn)本文曲波閾值去噪方法能否用于流動觀測數(shù)據(jù)，并對結(jié)果進(jìn)行評價.我們對臺站S05記錄的280個(Mb>5.0)接收函數(shù)進(jìn)行曲波去噪，接收函數(shù)的提取采用了基于臺陣觀測的三分量接收函數(shù)的多道最大或然性反褶積方法(劉啟元和Kind，2004)，其反方位角(back azimuth)范圍為30°～311°，震中距范圍為30° ≤Δ≤ 90°.

圖9給出了臺站S05接收函數(shù)的觀測數(shù)據(jù)和重建數(shù)據(jù)按反方位角10°的間隔分組疊加的結(jié)果.圖9表明：1)經(jīng)過曲波去噪的接收函數(shù)(徑向和切向)相較于原始數(shù)據(jù)顯示了更為清晰、合理的震相追蹤效果；2)對于反方位角200°～300°的嚴(yán)重缺失部分(已超過曲波長度)，接收函數(shù)的振幅不能精確恢復(fù)，但這并不破壞該區(qū)域之外的接收函數(shù)；3)對比去噪處理前后全方位接收函數(shù)疊加的結(jié)果(sum)，它們的波形沒有明顯不同；4)對方位嚴(yán)重缺失部分預(yù)測的接收函數(shù)相位追蹤與其他部分是協(xié)調(diào)的，表明其相位信息是有效的.

下面，我們進(jìn)一步考慮多臺陣列觀測的情況.對于二維線性觀測剖面來說，接收函數(shù)的曲波閾值去噪在一定程度上類似于地震勘探的情況(Herrmann and Hennenfent，2008；彭才等，2008).兩者的不同在于，用于天然地震觀測的線性臺陣難以保證空間上的等間隔采樣.這對于本文考慮的曲波閾值去噪來說并沒有原則上的困難.與單個臺站的情況不同，這種觀測方式吸引人的地方在于可以避免方位缺失帶來的問題.

作為一個例子，我們考慮沿川西臺陣31°N測線上19個臺站的接收函數(shù).臺站分布如圖8a所示.其剖面總長度約為420 km.臺站間距約為5～40 km(劉啟元等，2009).圖10給出了川西臺陣31°N測線接收函數(shù)去噪結(jié)果及其與原始數(shù)據(jù)的對比.所用的地震為2006年12月7日的千島群島MS6.4地震(震中位置46.26°N，153.80°E，深度16 km，發(fā)震時刻19h 10m 23.2s UTC).提取接收函數(shù)時，仍采用了三分量多道最大或然性反褶積方法(劉啟元和Kind，2004).由于該事件在臺站S04和S14的記錄質(zhì)量問題，我們用另外兩個地震事件的相應(yīng)數(shù)據(jù)加以替換.對臺站S04，我們選用了2007年4月9日的千島群島MS5.6地震(震中位置48.38°N，154.84°E，深度44 km，發(fā)震時刻10h 18m 2.3s UTC).對臺站S14，我們選用了2007年10月25日的千島群島以東MS5.9地震(震中位置45.96°N，153.63°E，深度10 km，發(fā)震時刻為13h 50min 4.0s UTC).它們的震中距和方位角與被替換數(shù)據(jù)大體相同.

對于圖10可以有如下觀察：1)接收函數(shù)的信噪比得到明顯提高，同時震相沿測線的可追蹤性增強(qiáng)；2)與松潘—甘孜地塊內(nèi)的臺站不同(S05—S13)，四川盆地內(nèi)各臺站(S14—S19)的接收函數(shù)去噪后顯示了簡單的波形特征，這意味著四川盆地具有較為均勻且各向異性較弱的地殼，而松潘—甘孜地塊內(nèi)各臺站接收函數(shù)較強(qiáng)的能量則與該地殼內(nèi)的低速層及下地殼的變形方式有關(guān)(Liu et al.，2014)；3)觀察零時刻前的背景噪聲，并與圖9相比，我們發(fā)現(xiàn)對于陣列觀測方式來說，其去噪效果不如單臺觀測的方式(例如，臺站S05)的好.

6結(jié)論與討論

通過把曲波變換和壓縮感知引入接收函數(shù)的數(shù)據(jù)處理，本文發(fā)展了一種二維多道接收函數(shù)的噪聲壓制方法.數(shù)值檢驗(yàn)和實(shí)際觀測數(shù)據(jù)的測試結(jié)果表明，本文發(fā)展的接收函數(shù)曲波變換去噪方法有很好的可行性，它既可以用于單個臺站觀測的情況，也可以用于多臺陣列觀測的情況.上述兩種情況，其壓制橫向非均勻散射的效果是明顯的.這有利于基于水平分層地殼模型的接收函數(shù)研究，這意味著以往所需要的疊加數(shù)據(jù)量可以減少，從而減小由于疊加造成的多次波能量削弱，并有助于縮短流動觀測所需的時間.

對以地殼各向異性研究為目標(biāo)的接收函數(shù)偏振分析來說(齊少華等，2009)，這意味著大幅度減少了散射噪聲的干擾，有助于提高接收函數(shù)偏振分析的穩(wěn)定性.由于本文討論的去噪方法可以同時彌補(bǔ)受地震方位控制的數(shù)據(jù)缺失，這有助于依據(jù)接收函數(shù)隨方位變化規(guī)律的地殼各向異性研究(田寶峰等，2008；Schulte-Pelkum and Mahan，2014).

但是，本文結(jié)果表明，單臺觀測數(shù)據(jù)的去噪效果優(yōu)于多臺陣列觀測的情況.其原因在于，單個臺站不同方位接收函數(shù)波場中，橫向散射波在空間取向上的隨機(jī)性強(qiáng)于線性剖面測線的情況.單事件的陣列觀測對接收函數(shù)波場的空間采樣依賴于臺站布設(shè)方式.相比之下，單臺多方位接收函數(shù)波場的空間采樣更加接近隨機(jī)分布.理論上，曲波去噪趨向于對相干性較強(qiáng)信息的保護(hù)，而對在空間取向上隨機(jī)性強(qiáng)的橫向散射波趨向于壓制.

對于我國大部分地區(qū)來說，全球地震分布導(dǎo)致的接收函數(shù)方位缺失不利于利用本文方法正確預(yù)測缺失方位的數(shù)據(jù)，但這并不妨礙單純的去噪目的.如何進(jìn)一步克服這種困難有賴于今后進(jìn)一步的深入研究.

致謝作者感謝評閱人對本文提出的寶貴建議，這使本文得以進(jìn)一步完善.本文所用數(shù)據(jù)來源于中國地震局地質(zhì)研究所國家重點(diǎn)實(shí)驗(yàn)室和IRIS全球地震臺網(wǎng).本項(xiàng)研究使用了Curvelab，Spot和SPGL1工具包.Curvelab：http:∥www.curvelet.org/.Spot：http:∥www.cs.ubc.ca/labs/scl/spot/.SPGL1：https:∥www.math.ucdavis.edu/～mpf/spgl1/.

References

Candès E J,Donoho D L.2000.Curvelets—a surprisingly effective nonadaptive representation for objects with edges.∥ Cohen A,Rabut C,Schumaker L L,eds.Curve and Surface Fitting: Saint-Malo.Nashville: Vanderbilt University Press,105-120.

Candès E J,Donoho D L.2004.New tight frames of curvelets and optimal representations of objects with piecewise-C2singularities.Commun.Pure Appl.Math.,57(2): 219-266.

Candès E J,Demanet L,Donoho D L,et al.2006a.Fast discrete curvelet transforms.Multiscale Model.Simul.,5(3): 861-899.

Candès E J,Romberg J,Tao T.2006b.Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information.IEEE Trans.Inform.Theory,52(2): 489-509.

Candès E J,Wakin M B.2008.An introduction to compressive sampling.IEEE Signal Proc.Mag.,25(2): 21-30.

Chai C P,Ammon C J,Maceira M,et al.2015.Inverting interpolated receiver functions with surface wave dispersion and gravity: Application to the western U.S.and adjacent Canada and Mexico.Geophys.Res.Lett.,42(11): 4359-4366.

Donoho D L.2006.Compressed sensing.IEEE Trans.Inform.Theory,52(4): 1289-1306.

Douma H,de Hoop M V.2007.Leading-order seismic imaging using curvelets.Geophysics,72(6): S231-S248.

Duecker K G,Sheehan A F.1997.Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track.J.Geophys.Res.,102(B4): 8313-8327.

Fang L H,Wu J P.2009.Effects of dipping boundaries and anisotropic media on receiver functions.Progress in Geophys.(in Chinese),24(1): 42-50.

Hennenfent G,Herrmann F J.2008.Simply denoise: Wavefield reconstruction via jittered undersampling.Geophysics,73(3): V19-V28.

Herrmann F J,Verschuur E.2004.Curvelet-domain multiple elimination with sparseness constraints.∥ Expanded abstracts of 74th SEG Annual Internat.Mtg.,1333-1336.

Herrmann F J,Hennenfent G.2008.Non-parametric seismic data recovery with curvelet frames.Geophys.J.Int.,173(1): 233-248.

Langston C A.1979.Structure under Mount Rainier,Washington,inferred from teleseismic body waves.J.Geophys.Res.,84(B9): 4797-4762.

Levin V,Park J.1997.P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation.Geophys.J.Int.,131(2): 253-266.

Ligorría J P,Ammon C J.1999.Iterative deconvolution and receiver-function estimation.Bull.Seism.Soc.Am.,89(5): 1395-1400.

Liu Q Y,Kind R.2004.Multi-channel maximal likelihood deconvolution method for isolating 3-component teleseismic receiver function.Seismology and Geology (in Chinese),26(3): 416-425.

Liu Q Y,Chen J H,Li S C,et al.2009.The MS8.0 Wenchuan earthquake: preliminary results from the western Sichuan mobile seismic array observations.Seismology and Geology (in Chinese),30(3): 584-596.

Liu Q Y,Li Y,Chen J H,et al.2009.Wenchuan MS8.0 earthquake: preliminary study of the S-wave velocity structure of the crust and upper mantle.Chinese J.Geophys.(in Chinese),52(2): 309-319.

Liu Q Y,van der Hilst R D,Li Y,et al.2014.Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults.Nature Geoscience,7: 361-365.

Nagaya M,Oda H,Akazawa H,et al.2008.Receiver functions of seismic waves in layered anisotropy media: application to the estimate of seismic anisotropy.Bull.Seism.Soc.Am.,98(6): 2990-3006.

Naghizadeh M,Sacchi M D.2010.Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data.Geophysics,75(6): WB189-WB202.

Neelamani R,Baumstein A I,Gillard D G,et al.2008.Coherent and random noise attenuation using the curvelet transform.The Leading Edge,27(2): 240-248.

Owens T J,Zandt G,Taylor S R.1984.Seismic evidence for an ancient rift beneath the Cumberland Plateau,Tennessee: A detailed analysis of broadband teleseismic P waveforms.J.Geophys.Res.,89(B9): 7783-7795.

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

Qi S H,Liu Q Y,Chen J H,et al.2009.Wenchuan earthquake MS8.0: Preliminary study of crustal anisotropy on both sides of the Longmenshan faults.Seismology and Geology (in Chinese),31(3): 377-388.

Schulte-Pelkum V,Mahan K H.2014.A method for mapping crustal deformation and anisotropy with receiver functions and first results from USArray.Earth Planet.Sci.Lett.,402: 221-233.

Shang X F.2014.Inverse scattering: theory and application to the imaging of the Earth′s seismic discontinuities [Ph.D.thesis].Massachusetts: Massachusetts Institute of Technology.

Stehly L,Cupillard P,Romanowicz B.2011.Towards improving ambient noise tomography using simultaneously curvelet denoising filters and SEM simulations of seismic ambient noise.Comptes Rendus Geoscience,343(8-9): 591-599.

Stehly L,Fromnet B,Campillo M,et al.2015.Monitoring seismic wave velocity changes associated with the Mw7.9 Wenchuan earthquake: increasing the temporal resolution using curvelet filters.Geophys.J.Int.,201(3): 1939-1949.

Tian B F,Li J,Wang W M,et al.2008.Crust anisotropy of Taihangshan mountain range in north China inferred from receiver functions.Chinese J.Geophys.(in Chinese),51(5): 1459-1467.

Van den Berg E,Friedlander M P.2008.Probing the Pareto frontier for basis pursuit solutions.SIAM J.Sci.Comput.,31(2): 890-912.

Yao H J,Gerstoft P,Shearer P M,et al.2011.Compressive sensing of the Tohoku-Oki Mw9.0 Earthquake: Frequency-dependent rupture modes.Geophys.Res.Lett.,38(20): doi: 10.1029/2011GL049223.

Yao H J,Shearer P M,Gerstoft P.2013.Compressive sensing of frequency-dependent seismic radiation from subduction zone megathrust ruptures.Proc.Natl.Acad.Sci.USA,110(12): 4512-4517.

Zhang J H,Zheng T Y.2015.Receiver function imaging with reconstructed wavefields from sparsely scattered stations.Seismol.Res.Lett.,86(1): 165-172.

Zhang J H,Lü N,Tian L Y,et al.2005.An overview of the methods and techniques for seismic data noise attenuation.Progress in Geophysics (in Chinese),20(4): 1083-1091.

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

Zhu L P,Kanamori H.2000.Moho depth variation in southern California from teleseismic receiver functions.J.Geophys.Res.,105(B2): 2969-2980.

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

房立華,吳建平.2009.傾斜界面和各向異性介質(zhì)對接收函數(shù)的影 響.地球物理學(xué)進(jìn)展,24(1): 42-50.

劉啟元,Kind R.2004.分離三分量遠(yuǎn)震接收函數(shù)的多道最大或然性反褶積方法.地震地質(zhì),26(3): 416-425.

劉啟元,陳九輝,李順成等.2008.汶川MS8.0地震: 川西流動地震臺陣觀測數(shù)據(jù)的初步分析.地震地質(zhì),30(3): 584-596.

劉啟元,李昱,陳九輝等.2009.汶川MS8.0地震: 地殼上地幔S波速度結(jié)構(gòu)的初步研究.地球物理學(xué)報,52(2): 309-319.

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

齊少華,劉啟元,陳九輝等.2009.汶川MS8.0地震: 龍門山斷裂兩側(cè)地殼各向異性的初步研究.地震地質(zhì),31(3): 377-388.

田寶峰,李娟,王衛(wèi)民等.2008.華北太行山區(qū)地殼各向異性的接收函數(shù)證據(jù).地球物理學(xué)報,51(5): 1459-1467.

張軍華,呂寧,田連玉等.2005.地震資料去噪方法、技術(shù)綜合評述.地球物理學(xué)進(jìn)展,20(4): 1083-1091.

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

(本文編輯胡素芳)

基金項(xiàng)目國家自然科學(xué)基金(41274060)資助.

作者簡介齊少華，男，1983年生，博士研究生，主要從事地震各向異性和接收函數(shù)研究.E-mail:ricown@163.com *通訊作者劉啟元，男，1945年生，研究員，主要從事寬頻帶地震學(xué)和地球深部結(jié)構(gòu)研究.E-mail:qyliu@ies.ac.cn

doi:10.6038/cjg20160311 中圖分類號P315

收稿日期2015-05-14，2016-01-10收修定稿

Attenuation of noise in receiver functions using curvelet transform

QI Shao-Hua，LIU Qi-Yuan*，CHEN Jiu-Hui，GUO Biao

StateKeyLaboratoryofEarthquakeDynamics，InstituteofGeology，ChinaEarthquakeAdministration，Beijing100029，China

AbstractSuppressing the scattering induced by the laterally heterogeneous media is an important problem for the imaging of the crustal structure and its anisotropy based on the stratified model from receiver functions.Although the scattering is able to be suppressed,to some degree,with stacking technique and low-pass filtering，these may lead to undesired waveform distortion,information loss or resolution reduction.To avoid these problems，we make use of the curvelet transform technique,which is developing rapidly in recent years,to reduce the scattering field in the receiver functions.Unlike exploration seismology，our major challenge is the uneven spatial sampling of the receiver functions,which is caused by the spatially incomplete distribution of stations and earthquake events.To overcome the difficulty，we incorporate the compressed sensing with the curvelet-based denoising to realize simultaneously the denoising and wavefield reconstruction.To verify our method,we test the denoising and wavefield reconstruction with synthetic receiver functions，and then apply our method to the observed data at one of the IRIS stations and the western Sichuan array,respectively.The results show： 1) Our method is effective in suppressing the scattering induced by the laterally heterogeneous media of the crust,so that the signal-to-noise ratio and the spatial traceability of the receiver functions are enhanced significantly; this may help the waveform imaging of the crustal structure and anisotropic parameters; 2) Missing data can be reconstructed correctly; 3) Our method can be applied to cases including single station and array observations,but it is more effective in the single station case than the array observations.

KeywordsReceiver function； Curvelet transform； Compressive sensing； Laterally heterogeneous scattering； Crustal anisotropy

齊少華，劉啟元，陳九輝等.2016.接收函數(shù)的曲波變換去噪.地球物理學(xué)報，59(3):884-896,doi:10.6038/cjg20160311.

Qi S H,Liu Q Y,Chen J H,et al.2016.Attenuation of noise in receiver functions using curvelet transform.Chinese J.Geophys.(in Chinese),59(3):884-896,doi:10.6038/cjg20160311.