頻率-波數域方法的發展及其在臺陣數據分析中的應用

王芳 王寶善

中國地震局地球物理研究所地震觀測與地球物理成像重點實驗室，

北京市海淀區民族大學南路5號 100081

0 引言

自1960年起，地震臺陣作為一種新型的地震學研究工具，已經廣泛應用于不同類型波場信號的檢測及地球內部不同尺度結構的成像（Birtill et al，1965；Frosch et al，1966；；Wright，1972；Doornboos et al，1972；Kvaerna，1989；Weber et al，1996；Ritter et al，2001；Krüger et al，2001）。一個地震臺陣由若干個位置接近、類型相同的地震計組成，能夠保證輸出波形的相干性，其記錄到的信號包含了由不同慢度和方位角的波疊加而成的波場信息，相對于單臺或大尺度臺網觀測有明顯的優勢（Rost et al，2002）。利用地震臺陣的特點，用頻率-波數分析（F-K）方法能夠同時獲取入射波場的慢度和方位角，因此這已發展成為一種重要的臺陣數據處理手段（Brooks et al，2009；Koper et al，2010b；Picozzi et al，2010；Gal et al，2015）。

傳統的F-K方法（lacoss et al，1969；Capon，1969a）最初主要用于單個頻率點的窄帶信號估計，Kvaerna等（1986a）提出非相干信號平均法（the incoherent averaged signal method，即IAS）并將其應用于寬頻帶信號的F-K分析中；Gal等（2014）通過加入對角加載因子（Diagonal Loading，即DL）對上述方法進行了改進，發現基于非相干平均法的F-K方法具有減小慢度估計誤差，提高F-K譜穩定性的優勢。此外，在基于臺陣進行微震噪聲的研究中，不同頻率的背景噪聲在傳播過程中會產生來自不同方位角的多次波，從而增加了噪聲源的檢測難度（Brooks et al，2009；Koper et al，2010a；Traer et al，2012；Gal et al，2015；Gal et al，2016）。傳統F-K方法的分辨率和精度有限，難以滿足檢測到所有弱源的需求。研究表明，臺陣布局包括臺陣孔徑、幾何形狀及地震臺數等是影響分辨率的重要因素（Okada，2003；Nishida et al，2008；Picozzi et al，2010）。臺陣響應函數（Array Response Function，即ARF）逐漸被引入作為衡量臺陣布局的定量化因子（Horike，1985；Halldorsson et al，2009）。為了減弱臺陣響應函數對功率譜的模糊效應，Nishida等（2008）利用Richardson-Lucy（即RL）反卷積方法（Richardson，1972；Lucy，1974）對臺陣響應函數與由最大似然法獲取的F-K圖像進行反卷積，提高了成像分辨率；Picozzi等（2010）利用RL和Tikhonov Regularization（即TR）（Bertero et al，1998）兩種反卷積方法從F-K譜中去除了臺陣響應函數，通過比較發現RL反卷積優于TR反卷積，能夠有效提高分辨率；Gal等（2016）利用原本用于天文學的CLEAN-PSF方法（Sijtsma，2007），通過不斷迭代壓制強干擾源，從而提高了對弱噪聲源的檢測能力。

各種不同類型臺陣的陸續建設，為F-K分析提供了大量的數據支持，使得該方法在理論和實踐中都得到了極大的發展。人們對F-K方法的不斷改進（Kvaerna et al，1986b；Gal et al，2014；Gal et al，2016；Picozzi et al，2010），使其應用領域也從分析地震事件拓展到研究包括噪聲在內的其他多種不同類型的波場（Rost et al，2002；Rost et al，2009；Koper et al，2010b）。本文首先介紹了F-K分析方法的發展，然后就其在檢測微弱信號源、分析噪聲特征、提取面波頻散曲線及評估臺陣設計優劣等4方面的應用進行綜述。

1 F-K分析方法

從原理上來說，F-K方法是基于頻率波數域功率譜密度函數（F-K PSDF）的計算，將臺陣觀測信號從時間-空間域轉換到頻率-波數域，從而獲取入射波場能量在不同慢度和方位角的分布。根據窗函數的不同，可分為聚束法（Lacoss et al，1969）、最大似然法（Capon，1969a）和多信號分類法（Schmidt，1986）等。由于地震臺陣記錄到的都是多源波場信號，多信號分類法很難確定信號子空間的振源數量（師黎靜，2007），在此不再介紹。聚束法和最大似然法是為了進行核爆檢測而提出的，是目前最常用的兩種F-K分析方法。

1.1 F-K譜

Lacoss等（1969）給出了基于聚束法的F-K譜估計：

Capon（1969a）提出了基于最大似然法的F-K譜估計：

此處，f表示頻率，k為波數，N是臺陣中所包含的地震臺站數目，φnm是第n和第m個臺站記錄之間的互功率譜，Xn和Xm分別為第n和第m個臺站的坐標。

上述兩種方法最主要的不同是窗函數。聚束法在波數k0的窗函數可表示為

最大似然法的窗函數為

其中

由公式（2）～（4）可以看出，聚束法的窗函數僅與臺陣形狀有關，相當于使用了固定的空間濾波器，最大似然法的窗函數不僅依賴于臺陣形狀，而且取決于數據的質量，相當于一種優化的空間濾波器，是一種非線性方法。眾多研究表明，聚束法提供了頻率-波數功率譜漸進無偏的一致估計，而最大似然法給出了頻率-波數功率譜的最大似然估計，具有更高的分辨率（Capon，1969a；Horike，1985；師黎靜，2007；吳建明，2007）。圖1給出了利用兩種方法所獲取的三維F-K譜示例，從中可以發現，聚束法的窗函數通常有較高的旁瓣，可能會得到虛假的F-K譜，而最大似然法的分辨率較高，所獲取的F-K譜比較平滑。

圖1 兩種F-K方法估計的三維頻率-波數譜。（a）聚束法；（b）最大似然法（吳建明，2007）

在頻率f一定時，假定F-K圖像的峰值出現在波數坐標（kx，ky），則此處對應的視速度為

入射信號的后方位角為

他的胸膛劇烈地起伏。久聞云浮山一帶有妖，自己從未深信，今天老天賞臉，算是給自己撞見了活的。這能夠修煉出人腦袋的妖怪，該是有幾百年的道行？

1.2 基于非相干平均法的F-K方法

上面提到的兩種F-K方法主要應用于窄帶信號的計算。從窄帶拓展到寬頻帶有兩種方法，包括相干信號平均法（the coherent averaged signal method；Wang et al，1985；Westwood，1992）和非相干信號平均法（the incoherent averaged signal method，即IAS；Wax et al，1984；Kvaerna et al，1986a；Kvaerna et al，1986b；Baggeroer et al，1993）。前者容易混淆空間位置相近的信號源，而后者對于分析具有較強頻率依賴性的噪聲源效果很好，能夠有效提高F-K譜的穩定性。Gal等（2014）利用模擬及實際觀測的噪聲數據，驗證了非相干平均法對于減小慢度估計誤差的優勢，并通過加入對角加載因子進一步提高了F-K譜的穩定性。

基于寬頻帶非相干平均法的最大似然F-K譜（加入了對角加載因子）是若干個離散頻率點的譜總和，可表示為

其中，L是離散頻率點的個數，α表示DL參數。

圖2給出了最大似然法和非相干平均法所獲取的F-K譜，從中可以看出，盡管兩種方法都能夠檢測到v＝3.33km／s和v＝4.16km／s附近的Rg和Lg震相，但也有一定的差異，相對于前者，后者能夠檢測到更多不同類型的源，如慢度較小的次級體波噪聲源。

圖2 （a）用最大似然法法獲取的F-K譜；（b）基于非相干平均法獲取的F-K譜（Gal et al，2014）

1.3 去除臺陣響應函數的F-K方法

在過去的數十年間，全球范圍內發展了多種不同類型的臺陣用于地震學的研究，如加拿大的十字形黃刀臺陣（the Yellow knife Array，即YKA）、美國的大孔徑臺陣（the Large Aperture Seismic Array，即LASA）、德國的寬頻帶大孔徑臺陣（（the Large-aperture Grafenburg Array，即GRF）、中國的上海佘山臺陣等。臺陣的布局是影響臺陣對不同頻率和慢度信號分辨能力的重要因素，可以用臺陣響應函數（ARF）定量表示。對于一個給定的臺陣，一些通用的準則用于衡量孔徑、地震臺數、臺間距、臺陣幾何形狀等因素對臺陣響應函數的影響（Harjes et al，1973；Johnson et al，1992）?？讖绞侵概_陣的最大臺間距，孔徑控制著F-K圖像主瓣的尖銳程度，孔徑越大，能夠測量到的波數越??；地震臺站數目決定臺陣作為一個波數濾波器的分辨能力；臺間距控制F-K圖像的旁瓣位置，平均臺間距越小，在給定波速的情況下，可識別的震相波長越?。慌_陣的幾何形狀控制入射波場的方位角分布。

臺陣響應函數在波數k0處可以表示為（Horike，1985）

利用RL反卷積方法可以去掉臺陣響應函數的影響，從而獲取更高分辨率的F-K譜。RL反卷積方法可以表示為

其中，A表示臺陣響應函數，m是迭代次數，g為原始的F-K譜，Eδm是經過m次迭代后的正則化解。

去除PSF旁瓣后的互功率譜矩陣可以表示為

其中，Ф是約束因子，控制去除能量所占的比例（例如，Ф＝0.05表示去除最強信號源能量的5%），wmax為每對互功率譜矩陣的權重因子矢量，w*max是權重因子矢量的復共軛，Pmax代表功率譜中的最大值。

去除最強信號源后，重構功率譜為

其中，M表示迭代次數，ФPimax是第i次迭代時清除掉的能量。

最后得到的F-K譜為

2 F-K方法在臺陣數據分析中的應用

本文主要就F-K方法在檢測微弱信號源、分析噪聲特征、提取面波頻散曲線、評估臺陣設計優劣等4方面的應用進行綜述。

2.1 檢測不同類型的微弱信號源

地震臺陣能夠在一定條件下記錄到一致性較好的波形，采用F-K方法處理臺陣數據，可以區分不同方向入射波場的能量分布，聚束具有相同頻率、相同慢度但不同入射方向的信號，因而能夠大幅度降低全球地震、核爆等不同類型事件的檢測下限，這對于識別被強震掩蓋的微弱事件具有良好的效果。

Li等（2016）利用ChinArray臺陣數據波形識別2013年的俄羅斯隕石墜落事件時，發現由于同期湯加地震事件的干擾，很難從原始波形中分辨出相對微弱的隕石墜落信號。利用Vespagram及F-K方法對湯加地震事件的面波信號和俄羅斯隕石墜落事件進行分離，發現最終得到的震相方位角和慢度與隕石墜落事件產生的面波一致。進而應用同樣的方法分析了F-net的觀測記錄，粗略定位了該事件的位置。說明F-K臺陣處理方法對于識別較遠處的微弱事件有較好的效果。

Gupta（1990）應用寬頻帶F-K方法，通過滑動時窗法計算了各個窗口與包括核爆P波初至在內的窗口的F-K譜殘差，識別出2個不同臺陣附近的近地表散射源，說明F-K方法有助于定位臺陣附近的散射源。佟玉霞（2002）將同樣的方法應用于上海地震臺陣記錄到的8個遠震事件，識別出了來自相似方位和距離的短周期面波，推斷是由遠震P波在上海臺陣南面的浙北平原與浙南山區交界區及北面的蘇北盆地邊緣發生散射引起的。

2.2 分析微震噪聲特征

現代地震儀的連續記錄中，僅有一小部分是地震事件，而大多是由海洋波和自然界生物活動引起的微震噪聲。眾多研究表明，這些往往被視為干擾信號的噪聲記錄是進行地下精細結構成像，監測斷裂帶、火山等重要的四維地質構造，研究與海洋風暴相關的氣候變化等的一種有效工具（Bromirski，1999；Grevemeyer et al，2000；Bromirski et al，2002；Shapiro，2004；Roux et al，2005；Gerstoft et al，2006；Chevrot et al，2007；Brenguier et al，2008a，2008b；Bensen et al，2008；Stutzmann et al，2009；Aster et al，2010；Zhan et al，2010；K?hler et al，2011；Tkal?i′cc et al，2012）。而噪聲信號的特征對于利用其開展這些科學研究具有重要意義，F-K方法可以根據頻率、慢度及傳播方向分離不同的信號，是一種研究噪聲源特征的有效方法。

美國LASA于1965年完成建設并投入使用，使得利用F-K方法在進行噪聲特征分析上取得了一系列突破（Toks?z et al，1968；Lacoss et al，1969；Capon，1969b；Haubrich et al，1969；Cessaro et al，1989）。研究表明，1～3s只有遠震體波噪聲，起源于遠洋風暴，與區域的活動構造關聯很?。?～7s除了遠震體波之外，還存在瑞雷波占主導地位的高階面波信號；長周期7～33s只有面波信號，可分別從垂直和水平分量中識別出基階瑞雷面波和勒夫面波，Toks?z等（1968）推斷勒夫面波是由于加諸在固體介質上的剪切力或者傳播路徑上的地殼不規則性引起瑞雷波發生轉換而產生的。同時，如果是在近海岸區域或大型湖泊附近，基階瑞雷面波的頻率會提高。其中，Cessaro等（1989）利用最大似然F-K方法，聯合ALPA臺陣（the Alaskan Long Period Array），不僅獲取了20s的長周期噪聲特征，而且采取三角測量的方法定位到了大西洋和太平洋風暴兩個噪聲源，利用滑動窗F-K分析方法定位到的噪聲源方位角比較穩定，并不隨海洋風暴位置的變化而變化，因此認為微震噪聲源與遠洋風暴的位置并不存在直接的關聯，而是由于海洋波及其在海洋基底反射波的相互干涉引起的。

隨后，越來越多不同類型臺陣的建立更加豐富了F-K方法在噪聲分析中的應用。Bungum等（1985）利用F-K方法分析了NORSAR（the Norwegian Seismic Array）臺陣數據，提取了周期低至0.125s的信號，相速度介于3～5km／s，認為是瑞雷面波。Friedrich等（1998）分析了GRF臺陣的寬頻帶F-K譜，發現主微震（12～20s）和次級微震（6～11s）都以基階勒夫波和瑞雷波為主，且兩個頻帶內勒夫波與瑞雷波之間的能量比是不同的，前者比是6∶5，后者比是1∶4。Essen等（2003）對GRF臺陣4～10s的F-K分析結果顯示，該頻帶對應于3～3.5km／s的相速度，屬于基階瑞雷面波。Gerstoft等（2006）利用聚束F-K方法研究了颶風發生時期的150個南加州臺站記錄，在4～6s的頻帶發現了相速度為11.7km／s的來自上地幔的P波。Koper等（2008）和Koper等（2009）分別在CMAR（the Chilang Mai Array）的2s附近和YKA臺陣1～3s的F-K圖像上，觀測到了遠震P波和PKP波。

基于眾多單個臺陣的研究結果，Koper等（2010b）利用全球18個地震臺陣連續一年的垂直向噪聲記錄，分析了短周期噪聲（0.25～2.5 s）的構成。圖3給出了一個圓形臺陣的F-K分布圖。幾個不同半徑的圓分別對應于不同的慢度：4.4s／deg，13.9s／deg，27.8s／deg，37.1s／deg。圖4給出了利用18個臺陣2007～2008年的噪聲記錄獲得的相速度統計結果，占主導地位的是勒夫波，推斷來自淺層或近海岸的海洋區域；超過25%的是各種類型的P波，包括上地幔的轉換波（P，14%），近震P波（Pn／Pg，6%）和外核折射波（PKP，8%），推斷其來源之一是遠離海岸線的北太平洋地區，該區域是一個穩定的長期微震噪聲激發源，已經激發了近40年的P波能量；基階瑞雷面波Rg占據的比例最小，可能是由于其衰減較快，只有接收臺陣距離激發源較近時方能觀測到。

圖3 一個圓形臺陣的F-K慢度分布示例圖

圖4 利用18個臺陣2007～2008年的噪聲記錄統計獲得的相速度分布（Koper et al，2010b）

2.3 提取面波頻散曲線

地震面波在傳播過程中最重要的特點之一是存在頻散現象，這一特性是由地下介質的結構和非彈性決定的。因此，利用面波的頻散曲線，通過合適的反演方法，便可以獲取地下介質的速度結構。假設面波是噪聲的主要成分，利用F-K方法將空間-時間域的噪聲記錄轉換到頻率-波數域后，不同頻率的面波在二維波數空間上會出現相應的峰值；確定該峰值的位置即可得到波數，從而能夠計算出不同頻率面波的相速度，獲取面波的頻散曲線。如果研究三分量的連續噪聲記錄，便可以同時提取瑞雷波和勒夫波的頻散曲線。

Okada（1987）利用F-K方法估算場地的瑞雷波頻散曲線，并反演了場地的剪切波速度結構。Huang等（1990）利用最大似然法分析了臺灣SMART地區的噪聲記錄，得到了該地區各個頻率下的F-K譜，并獲取了剪切波速度結構。Kavand等（2006）應用最大似然法研究了伊朗東南部Bam市區的噪聲記錄，獲得了瑞雷面波頻散曲線，并反演了該地區的剪切波速度結構。Tokeshi等（2006）針對兩個場地，通過比較模擬波場和實際波場所獲取的瑞雷波頻散曲線，認為F-K方法用于提取頻散曲線足夠可靠。張碧星等（2005）采用F-K方法對瑞雷波多模式頻散曲線進行了研究和分析，通過選取實際中經常碰到的3種分層介質模型，詳細分析了接收道數目、源檢距等參數對頻散曲線的影響，發現當源檢距大于半波長時利用F-K方法可以得到可靠的頻散曲線。吳建明（2007）利用最大似然法估算了廈門2個不同場地的面波頻散曲線，分析了窗函數、功率譜估算方法對F-K圖像的影響。Picozzi等（2010）利用Richardson-Lucy（RL）和Tikhonov Regularization（TR）的2種反卷積方法從F-K圖像中去除了臺陣響應函數，獲取了更高分辨率的F-K譜，并提取到了3.5～13Hz范圍內的瑞雷波頻散曲線。

2.4 評估臺陣設計的優劣

臺陣作為地震觀測網的重要組成部分，引起了地震學家的高度重視。出于不同的觀測目的，如核爆識別，繪制區域地震活動圖，監測遠震等，全球范圍內已經建設或正在籌建各種不同類型的臺陣。臺陣設計對于節約成本，確保各種監測目的實現的重要性不言而喻。作為一種重要的臺陣數據處理方法，F-K譜的分辨率是衡量臺陣布局是否合理的有效工具。

臺陣布設中需要考慮到臺陣形狀、測點數量、孔徑大小、臺間距等多種因素，可以用臺陣響應函數在頻率-波數域的值來定量描述這些因素的影響。Schweitzer等（2002）比較了十字形的YKA和小孔徑的ARCES兩個臺陣的響應函數，發現ARCES的幾何形狀比較完美，對不同方位角的分辨率很好，且旁瓣距離主瓣較遠，但由于孔徑較小，不容易識別具有微小差異的波數。而YKA由于孔徑比較大導致臺陣響應函數的主瓣很窄，在測量視速度時具有更高的分辨率，但是由于其幾何形狀的原因，在不同的方位角具有不同的分辨率。

吳建明（2007）利用模擬的平穩波場分析了5種形狀的臺陣，包括半徑變化的圓形臺陣、半徑不變的圓形臺陣、十字形臺陣、混合型臺陣及菱形臺陣，考察了相應F-K譜的穩定性，據此評估了各個臺陣的觀測能力。在測點數量相同的情況下，十字形臺陣和菱形臺陣得到的F-K譜比圓形臺陣及混合型臺陣更為穩定。

3 討論與展望

本文介紹了F-K分析的基本原理及其各種改進方法，并綜述了F-K方法在檢測微弱信號源、噪聲特征分析、提取面波頻散曲線、臺陣設計等4方面的應用。總體來說，F-K臺陣處理方法具有較高的分辨率，其方法的不斷改進也在拓寬它的應用領域。

（1）本文總結了幾種F-K分析方法，包括傳統的聚束法和最大似然法，基于非相干平均法的F-K方法，利用Richardson-Lucy（RL）反卷積和CLEAN-PSF去除臺站響應函數的F-K法。傳統的聚束法和最大似然法的分辨率有限，容易漏檢一些微弱的噪聲源；而基于非相干平均法的F-K法能夠拓寬所檢測噪聲源的頻帶，從而識別更多的微震噪聲。然而，眾多研究表明，臺陣里有限的采樣點及不合理的排列方式往往會限制F-K圖像的空間分辨率及其對不同頻率面波的分辨能力。以上幾種方法均沒有考慮到臺陣布局對F-K譜的影響，利用它們對不規則的臺陣進行F-K成像時會影響結果的分辨率。去除臺陣響應函數的F-K分析方法能夠彌補上述幾種方法的不足，可以在一定程度上提高F-K成像的分辨率，使得如圓形、L型、十字型或隨機分布的臺陣都能夠得到有效利用。

（2）除了文中所述的應用之外，考慮到F-K分析能夠降低地震等事件的識別下限，利用多個臺陣便可以實現對微小地震的精定位。此外，F-K分析能夠提高得微弱信號的信噪比，已經在定位近地表散射源中有了初步的應用，如何將其用于識別更多的反射波和散射波等微弱信號，從而拓寬地震臺陣在地球內部的精細結構成像方面的應用，需要更深入的探索。

（3）像其他大多臺陣方法一樣，F-K分析是基于平面波入射的假設，臺陣下方的各向異性有可能會改變波前并破壞信號的相干性，從而改變F-K分析的結果。研究高頻信號時采用的小孔徑臺陣的地下介質往往比較均勻，因而獲得的入射波方位角較穩定，但由于高頻微震信號源的變化，導致其視慢度值不穩定，可以通過滑動窗法獲取高頻噪聲的時空變化特性；研究低頻信號時需要臺間距相對較大的臺陣，視慢度值較穩定，但可能由于傳播路徑上的橫向不均勻性而造成方位角的偏離，長時間疊加可能有助于穩定F-K圖像的結果。

致謝審稿專家對本文提出了寶貴的修改意見和建議，謹表謝忱。