TSAnalysis軟件在遼寧地區面波頻散提取中的應用

雷 晨,趙龍梅,朱葉琳,王 喬

（遼寧省地震局，遼寧 沈陽 110034）

0 引言

進行基于地震面波的成像，首先需要從觀測臺站記錄到的大量觀測數據中，可靠地提取出對于研究有價值的基礎資料，最重要、最繁瑣也最基礎的工作就是面波相速度頻散曲線的測量。20世紀50年代，傅里葉變換技術應用于面波頻散測量，60年代初，數值濾波技術和時間變量濾波技術在面波頻散測量中得到應用，此后逐漸發展建立了面波頻時分析的基礎。Landisman等[1]提出，可以利用雙臺數據互相關提取雙臺間的面波群速度、相速度。時至今日，用雙臺波形互相關法來提取面波的頻散曲線已經是常用方法之一，面波頻散測量工作可以借助多種較為成熟的軟件開展，姚華建[2]提出了一種基于圖像分析的雙臺面波相速度頻散曲線提取方法，實現了快速、直觀、準確地頻散測量。近年來，基于遼寧及周邊地區固定地震臺站記錄的資料，開展了一系列地震面波層析成像相關研究工作，在基礎數據處理過程中，使用了姚華建編寫的TSAnalysis（雙臺相速度頻散曲線快速提取軟件）進行資料處理，得到了較好的應用效果。

1 雙臺互相關方法原理

按照提取一條頻散曲線所用到的資料來劃分，常用的頻散測量方法有單臺法、雙臺法和雙事件法等。所謂雙臺法，顧名思義，是基于成對臺站觀測數據開展相關研究的方法，且兩個臺站幾乎位于從事件到臺站的同一大圓路徑上；而互相關，是兩個數據之間關于時間的一個函數，反映兩個數據的相匹配程度。傳統的單臺法在測量中需要用到關于震源的大量參數，比如震源機制、發震時刻、震中位置、初始相位等等，而測定這些參數的同時會引入大量誤差項，拉低頻散的測量精度。相比之下，使用雙臺互相關方法處理觀測數據，就不用考慮震源的這些參數，可以在保留頻散信息的同時去掉震源的影響，從而可以提高很多測量精度[3-6]。雙臺互相關法的原理如下。

雙臺數據互相關譜函數為：

對于n階振型，上式(3)則可以重寫為：

由式（3）、（4）可知，雙臺數據經過互相關可以消除震源的相位，只保留下頻散信息。式（4）中第一項可以反映出面波在雙臺之間介質傳播的各階振型的頻散信息，而第二項為交叉振型項。如果假設在統計意義下，不同振型相互獨立，不同階振型互相關后疊加會相互抵消，那么式中的第二項可以忽略不計；走時信息屬于已知量，可以消去，式（4）可以變為：

當然，雙臺法也有一些應用限制，首先面波路徑要為大圓路徑，即所選雙臺與所選地震震中幾乎位于同一大圓路徑上，并且雙臺在震中的同側；其次所選雙臺之間的距離要遠小于震中距，即要使用遠震事件面波觀測資料；最后，所研究區域臺站分布不能太稀疏，還要篩選數量足夠且方位角分布較好的地震事件。

2 基于圖像分析的雙臺相速度頻散曲線快速提取軟件

TSAnalysis是使用Matlab編寫的一款圖形界面軟件[2,8]，界面如圖1所示。在遼寧地區面波成像研究中使用到該軟件的主要功能有：地震資料的預處理—包括地震波形讀取，數據降采樣，去除零漂、儀器響應，振幅歸一化；窄帶濾波器生成及濾波；基于多重濾波法計算雙臺群速度頻散；基于圖像分析提取雙臺相速度頻散曲線等。

圖1 TSAnalysis軟件圖形界面Fig.1 TSAnalysis graphic interface

在地震資料的預處理過程中，通過程序參數設置控制每一步計算，首先對大量采樣率為100Hz的地震事件數據降采樣，重采樣頻率為1Hz；然后對數據資料進行零漂校正、去儀器響應、波形歸一化等一系列處理；在窄帶濾波器的生成和濾波環節，選取了加Kaiser窗的有限沖擊響應濾波器，為獲得最好的窄帶濾波效果，即盡可能抑制高階面波和其他波的影響，很好的分離出基階面波，Kaiser窗的特征參數β值取為9。

TSAnalysis軟件的核心部分為計算能夠直接反映相速度和周期之間關系的振幅矩陣并用圖像的方法顯示，從而基于圖像清晰直觀的確定實際頻散曲線。為了更好的控制整個過程的正確性和精度，在使用多重濾波法處理雙臺記錄得到每個臺站群速度到時后，只選擇到時連續清晰的雙臺資料進行窄帶濾波；在采用二次樣條插值方法計算雙臺互相關振幅矩陣時，只選取互相關最明顯且相速度值符合天然地震頻散特征范圍的區域進行計算和插值變換。Yao et al[8,9]為提取高質量的相速度頻散曲線采用了多種技巧，使這種基于圖像分析技術提取瑞雷面波相速度頻散的方法在保證正確性和精度的同時非常直觀易于操作。

3 資料處理實例

選取2013年11月25日5:56:53在151.2oE，45.57oN發生的矩震級6級地震事件作為計算實例，地震震源深度35km。

首先加載全部臺站記錄的地震波形文件和所有設備的儀器響應文件，設置重采樣頻率為1Hz。設置加Kaiser窗的窄帶濾波器的測試中心周期為40s，Kaiser窗的特征參數β值取為9，濾波器帶寬設置為1s，于是生成的窄帶濾波器如圖2所示。

圖2 窄帶濾波器實例(a.時域；b.頻域)Fig.2 Narrowband filter example(a.Time domain, b.Frequency domain)

在程序篩選符合位于同一大圓路徑上的兩個臺站時，限定兩個最大偏移角α和β分別為3o和5o（α表示震中到兩個臺站的方位角差，β表示震中到較近臺站和兩臺站連線之間的方位角差）。此實例中所選取地震記錄是位于雙臺大圓路徑上的DHT(敦化臺)和ZUH(遵化臺)兩個臺站的波形記錄，經過程序降采樣和去儀器響應并窄帶濾波后提取的波形如圖3所示，DHT臺站距離震中較近，可見圖中上部為DHT臺站記錄波形，下部為ZUH臺站記錄波形。

圖3 經程序初步處理的大圓路徑上雙臺波形記錄Fig.3 Waveform of dual stations on great circle path

圖4 群到時圖像Fig.4 Image of group arrival time

根據此前試算遼寧地區數據時的經驗，利用現有數據提取周期15s至120s的相速度頻散是可行的，所以使用軟件提取頻散過程中對周期范圍略放大，起始周期10s，終止周期130s，周期間隔設置為1s。首先由提取的雙臺記錄得到雙臺群到時圖像，如圖4所示，顏色標識歸一化之后每個周期波包能量情況，紅色代表高能量，藍色相反，綠色曲線即為程序追蹤的群到時曲線。圖中可見紅色區域比較清晰且連續，說明這一組數據的質量合格，可以繼續進行下一步計算。加移動窗過濾可以進一步濾除高階振型和噪聲的影響，于是得到更加清晰的群速度到時圖像，如圖5所示。

圖5 加移動窗處理后的群到時圖像Fig.5 Image of group arrival time by add moving time window

至此，可以通過程序計算得到互相關振幅矩陣并直接用圖像法顯示，如圖6所示，可以很清晰的看到互相關最為明顯的區域(即黃顏色最深的區域)，并操作程序提取相速度頻散曲線。如圖中黑色實線為程序自動搜索得到的頻散曲線，紅色點為依照周期范圍設置所提取的滿足信噪比條件的各周期頻散點，這些數據可進入下一步分析。

圖6 基于圖像分析的雙臺相速度頻散曲線提取實例Fig.6 Example for the determination of inter-station Rayleigh wave phase velocity dispersion curves based on the image analysis technique

4 資料及處理結果

依照本次研究區域范圍，向遼寧省地震局監測中心和國家測震臺網數據備份中心[10-11]申請了遼寧省內及周邊鄰省共77個固定臺站寬頻帶地震計記錄的2012年1月至2014年6月間，矩震級滿足大于5.5級且小于7.5級，震源深度小于100km，震中距滿足大于10o小于100o的全部700多個地震事件垂直向波形數據。經過計算和挑選，最終從數百個地震事件中篩選出瑞利面波發育較好，具有比較完備的震中距和方位角覆蓋，且頻散特征明顯的105個地震事件用于本次研究。

圖7 雙臺路徑分布圖Fig.7 Distribution of paths

使用上述基于圖像分析的雙臺相速度頻散曲線提取方法，借助TSAnalysis軟件，共提取超過3000條周期范圍在15s到120s內的相速度頻散曲線(包括重復路徑)。經過仔細篩選和計算對比，仍保留了超過2600條的高質量重復路徑相速度頻散曲線，可見這種方法的質量控制效果還是較為高效及可信的。通過對重復路徑頻散資料取算數平均數，最后獲得超過1100條獨立雙臺路徑相速度頻散資料，對于整個遼寧地區及周邊的覆蓋情況如圖7所示，圖中亮黃色標識遼寧省邊界。圖8展示了從15s到120s之間不同中心周期的獨立路徑分布情況，整體覆蓋情況均滿足下一步工作開展的需求。

圖8 各周期獨立路徑數Fig.8 Number of paths at each interest period

5 結論

本文利用遼寧及鄰省固定地震臺觀測資料，使用TSAnalysis軟件，應用基于圖像分析的雙臺相速度頻散曲線提取方法處理了105個地震事件的波形記錄，成功提取了1100條獨立雙臺路徑相速度頻散資料，對于整個遼寧地區及周邊區域的覆蓋情況良好，高質量的相速度頻散資料為后續開展的地震面波層析成像工作奠定了良好基礎。相速度頻散的識別和分析是很多項地球物理研究的基礎，基于圖像分析的雙臺相速度頻散曲線提取方法除了本文展示的天然地震面波數據處理，還可以應用于噪聲成像及其他研究的基礎數據處理。相速度頻散的有效提取對于更好的研究淺部構造及深部地球動力學過程有著積極的意義。