基于諧波分解的余干地震臺(tái)接收函數(shù)研究

查小惠,趙 影,曾文敬,羅 麗,郭雨帆

（江西省地震局,南昌 330026）

關(guān)鍵字:接收函數(shù)；諧波分解；各向異性

0 引言

接收函數(shù)方法[1]研究地殼結(jié)構(gòu)已被廣泛應(yīng)用。隨著觀測(cè)資料的增多,研究方法的深入,基于不同方位角接收函數(shù)的變化研究臺(tái)站下方地殼不均勻性的論文越來越多[2-5]。目前在研究橫向不均勻性中研究最多的是界面傾斜和各向異性,而各向異性中又以HTI(水平對(duì)稱軸的橫向各向同性)介質(zhì)為主。在眾多使用接收函數(shù)研究地殼不均勻性的方法中,諧波分解[6]是一種較為常用的方法。諧波分解方法的主要原理是根據(jù)傾斜界面和各向異性導(dǎo)致的Pms震相隨方位角時(shí)間變化周期存在差別[7-8]。界面傾斜導(dǎo)致Pms震相隨方位角呈2π周期變化,HTI介質(zhì)中Pms震相隨方位角呈π周期變化,諧波分解方法可以在Pms震相的相對(duì)時(shí)差中分離出這兩種疊加在一起的時(shí)間變化,從而分析臺(tái)站下方的橫向不均勻性。本文首先基于水平對(duì)稱軸各向異性介質(zhì)接收函數(shù)Pms到時(shí)理論公式和理論地層模型,分析并驗(yàn)證諧波分解方法的可行性和可靠性。然后將該方法應(yīng)用到江西余干地震臺(tái)(圖1),基于江西省地震局臺(tái)網(wǎng)中心多年來積累的連續(xù)波形資料,進(jìn)一步驗(yàn)證諧波分解方法在實(shí)際臺(tái)站中的應(yīng)用效果,探討余干臺(tái)站下方的地殼不均勻性。

1 方法及理論測(cè)試

1.1 方法原理

根據(jù)前人研究[7-8]顯示,單層水平對(duì)稱軸各向異性介質(zhì)接收函數(shù)Pms到時(shí)可以近似表達(dá)為

其中,t0為各向同性介質(zhì)Pms震相到時(shí),δt是各向異性介質(zhì)所引起的時(shí)差,θ為接收函數(shù)的后方位角,φ是各向異性快波方向。根據(jù)接收函數(shù)不同后方位角Pms震相的到時(shí)變化,可以基于諧波分解的方法,反演得到公式(1)中的模型參數(shù),理解臺(tái)站下方的介質(zhì)結(jié)構(gòu)參數(shù)。

本文首先使用理論公式和理論模型得到一系列的理論P(yáng)ms震相到時(shí),然后利用harmfit函數(shù)對(duì)該到時(shí)進(jìn)行諧波擬合,通過理論試驗(yàn),驗(yàn)證分析程序的可靠性,同時(shí)理解各反演參數(shù)間的對(duì)應(yīng)關(guān)系。最后將該方法應(yīng)用到實(shí)際臺(tái)站當(dāng)中,獲取實(shí)際臺(tái)站下方可靠的地殼結(jié)構(gòu)參數(shù)。

1.2 理論接收函數(shù)測(cè)試

1.2.1基于公式的理論測(cè)試

基于公式(1),假設(shè)t0=5,δt=0.3,φ=45°。反方位角從0°到360°,以10°為間隔,總共37道,計(jì)算出Pms震相的走時(shí)(圖2)。然后對(duì)計(jì)算出的Pms震相走時(shí)使用harmfit函數(shù)進(jìn)行諧波分析,驗(yàn)證諧波分解程序的可靠性。研究發(fā)現(xiàn),使用harmfit進(jìn)行諧波分析時(shí),需要去除Pms震相走時(shí)的直流分量,否則擬合不準(zhǔn)確。去除掉直流分量后,使用harmfit進(jìn)行4階諧波分析,得到結(jié)果參數(shù)為表1,擬合圖形為圖2,圖2中實(shí)線為公式計(jì)算的理論P(yáng)ms震相去直流后的走時(shí),虛線為使用harmfit函數(shù)擬合得到的Pms震相走時(shí)。

表1 harmfit函數(shù)4階諧波分析結(jié)果Table 1 The fourth order harmonic analysis results of harmfit function

圖1 YUG地震臺(tái)位置圖Fig.1 Location map of YUG station

圖2 理論和諧波擬合的Pms震相走時(shí)圖Fig.2 Pms phase travel time diagram fitted with harmonics

根據(jù)表1可以看到,1,3,4階的振幅都為零,根據(jù)harmfit函數(shù)說明,得到的二階諧波表達(dá)式應(yīng)該為AMP*cos(2t+ALP),其中AMP為振幅,ALP為相位,單位是弧度。二階的振幅為0.15,振幅的2倍為0.3,和假設(shè)的δt=0.3一致。相位信息為1.57,為弧度結(jié)果,接近pi/2,二階表達(dá)式可以寫成0.15*cos(2*t+1.57),根據(jù)cos函數(shù)的公式cos(θ-pi)=-cos(θ)?？梢詫⒍A表達(dá)式寫成-0.15*cos[2(θ-pi/4)],和(1)式對(duì)照認(rèn)為,harmfit諧波分解結(jié)果和假設(shè)的參數(shù)是一致的。通過該理論公式和harmfit的函數(shù)擬合,驗(yàn)證了harmfit函數(shù)諧波分解的正確性,并知道如何理解harmfit的諧波分解參數(shù)。

1.2.2 基于模型的理論測(cè)試

為了更深入的理解各向異性地殼模型產(chǎn)生的接收函數(shù)Pms震相走時(shí)變化,諧波分解方法對(duì)地殼各向異性的分析能力,本文設(shè)定了兩個(gè)地殼模型(表2)。使用諧波分解方法分析理論地層模型的Pms震相走時(shí)數(shù)據(jù),驗(yàn)證了程序的可用性并更深入的理解了相關(guān)參數(shù)的對(duì)應(yīng)關(guān)系。

表2 計(jì)算接收函數(shù)的地層模型Table 2 Strata models used for receiver function modeling

對(duì)A1和A2兩個(gè)地層模型,使用RAYSUM程序[9]合成理論的RTZ三分量地震記錄,射線參數(shù)設(shè)為0.05 s/km,高斯脈沖寬度設(shè)為0.5 s,采樣間隔設(shè)為0.01 s,總共5024個(gè)采樣點(diǎn),360°的反向方位角以10°為間隔觀測(cè),總共37道接收。再將理論合成的RTZ三分量地震記錄利用迭代反褶積[10-11]得到R分量的接收函數(shù),接收函數(shù)高斯低通濾波參數(shù)設(shè)為3,大約對(duì)應(yīng)1 Hz的拐角頻率。求取得到理論接收函數(shù)以后使用腳本將反方位角、射線參數(shù)等信息寫入接收函數(shù)頭文件,方便后續(xù)處理。

圖3 A1模型Ps震相Fig.3 Ps phase of A1 model

圖4 A1模型Ps震相對(duì)齊時(shí)差曲線Fig.4 Ps phase alignment time difference curve of A1 model

根據(jù)Ps震相的到時(shí),選取接收函數(shù)時(shí)間窗3.5 s～6 s進(jìn)行分析(圖3)。以反方位角為0°的接收函數(shù)為參考道,其它道接收函數(shù)向該道接收函數(shù)對(duì)齊,計(jì)算出每道相對(duì)第一道的時(shí)間偏移量,可以看出各道接收函數(shù)的時(shí)間偏移連線顯示明顯的余弦曲線特征(圖4)。需要說明的是,時(shí)間偏移曲線的余弦特征和參考道的選取沒有關(guān)系,選取任何一個(gè)反方位角的接收函數(shù)作為參考道,得到的時(shí)間偏移曲線余弦特征不變,但是相對(duì)到時(shí)在數(shù)值上存在一個(gè)常量的整體的平移。利用harmfit函數(shù)對(duì)相對(duì)時(shí)差進(jìn)行諧波分析,可以得到不同階余弦函數(shù)的振幅和相位(表3),下面我們對(duì)諧波分解結(jié)果進(jìn)行分析。

表3 模型A1和A2諧波分解結(jié)果Table 3 Harmonic decomposition results of model A1 and A2

Raysum軟件中各向異性的定義是基于百分比,假設(shè)各向異性強(qiáng)度K可以寫為其中Vmax為快波速度,Vmin為慢波速度；快慢波走時(shí)差dt可以表達(dá)為其中L為各向異性層厚度。將(1)式代入(2)式可得(3)式L*K=dt*Vmin(3)。在A1和A2模型中L為20 km,模型中給出的是平均速度為3.5 km/s,各向異性強(qiáng)度K為5%時(shí),計(jì)算可得Vmin約為3.41 km/s,可以算出理論的時(shí)間延遲應(yīng)該為0.29 s左右。

根據(jù)harmfit函數(shù)說明,得到的2階諧波表達(dá)式應(yīng)該為AMP*cos(2t+ALP),ALP為弧度,可以認(rèn)為使用harmfit求解兩個(gè)理論模型的時(shí)間延遲分別為0.30 s和0.28 s,和理論時(shí)間0.29 s的時(shí)間延遲相符,快軸方向分別約為-6°和39°,和理論的0°和45°相近但有誤差?；诶碚撃Ｐ偷挠?jì)算,我們認(rèn)為該方法可以較好的提取單層水平對(duì)稱軸各向異性地殼模型的地殼結(jié)構(gòu)信息,后續(xù)將在實(shí)際資料處理中,進(jìn)一步驗(yàn)證。

2 實(shí)際臺(tái)站處理

基于前文理論公式和理論模型的諧波分解結(jié)果可以知道,只要精確獲取不同方位角的Ps震相相對(duì)時(shí)間偏差,通過諧波分解,可以得到各向異性的快波方向和時(shí)間延遲。但在實(shí)際處理過程中,接收函數(shù)存在噪聲,且接收函數(shù)方位角覆蓋不完整,這些都會(huì)影響到不同方位角相對(duì)時(shí)差的提取,進(jìn)而影響諧波分解的可靠性。本文基于江西余干地震臺(tái),細(xì)致分析了方位角覆蓋的完整性問題,基于波形互相關(guān)方法提取Pms的相對(duì)時(shí)差。然后利用諧波分解方法對(duì)臺(tái)站接收函數(shù)進(jìn)行了分析,對(duì)結(jié)果進(jìn)行了解釋。

2.1 方位角完整性分析

為了分析地震事件方位角覆蓋情況,本次以余干地震臺(tái)為中心,挑選出2011年到2019年5.5級(jí)以上地震,震中距在30°到95°之間記錄的遠(yuǎn)震總共1837個(gè)。以10為間隔,將反方位角分為36個(gè)區(qū)間,統(tǒng)計(jì)每個(gè)區(qū)間接收函數(shù)的數(shù)量(表4)。同時(shí)繪制地震震中分布圖(圖5)。根據(jù)表4和圖5可以明顯看出,地震事件的方位角分布是非常不均勻的。反方位角從50°到110°只有5個(gè)地震事件,地震事件為個(gè)位數(shù)的反方位角區(qū)間也還有9個(gè),如果這些地震事件求取的接收函數(shù)形態(tài)不好,會(huì)導(dǎo)致這些方位角的接收函數(shù)約束缺失,這對(duì)我們使用不同方位角接收函數(shù)反演地殼結(jié)構(gòu)是非常不利的。對(duì)于某些方位角,如120°～130°,地震事件數(shù)量達(dá)到665條,存在很大的冗余。雖然地震事件分布很不均勻,但地殼橫向不均勻性造成的Pms震相周期較大,根據(jù)采樣定理,當(dāng)前的方位角覆蓋完全可以對(duì)對(duì)信號(hào)進(jìn)行恢復(fù)。

表4 YUG臺(tái)2011—2019不同方位角接收函數(shù)數(shù)量Table 4 Number of receiver functions in different azimuths of YUG station from 2011 to 2019

圖5 2011—2019年遠(yuǎn)震震中分布圖Fig.5 Epicenter distribution of teleseismic from 2011 to 2019

圖6 使用的地震事件震中分布圖Fig.6 Epicenter distribution map of earthquakes used

在實(shí)際處理過程中,基于地震事件震中方位角分布的不均勻性,本文采取如下接收函數(shù)挑選方法。首先基于接收函數(shù)的形態(tài),在所有接收函數(shù)中挑選出Ps震相較為清晰,波形質(zhì)量較好的接收函數(shù),然后對(duì)于某一個(gè)方位角范圍,在挑選好的多個(gè)接收函數(shù)中挑選出一條高質(zhì)量的可代表該方位角接收函數(shù)特征的事件參與不同方位角的接收函數(shù)反演。而對(duì)于接收函數(shù)很少的方位角區(qū)間,逐條挑選,確保方位角的覆蓋能達(dá)到最大。

2.2 實(shí)際臺(tái)站處理

本次使用的地震資料時(shí)間跨度較長(zhǎng),使用了2007—2017年的地震資料。選擇震中距在30°～95°之間的5.5級(jí)以上遠(yuǎn)震,截取P波初至前20 s和初至后80 s事件波形進(jìn)行接收函數(shù)計(jì)算,帶通濾波頻率為0.1～2.0 Hz,迭代次數(shù)設(shè)置為100、接收函數(shù)的高斯低通濾波系數(shù)設(shè)為3.0(該系數(shù)對(duì)應(yīng)的拐角頻率約為1 Hz)。挑選到267條接收函數(shù),然后每個(gè)方位角選擇1條接收函數(shù)參與反方位角計(jì)算,總共有20條接收函數(shù)參與計(jì)算,參與計(jì)算的地震震中分布圖見圖6。

首先對(duì)20條接收函數(shù)進(jìn)行動(dòng)校正,校正到震中距為67°的標(biāo)準(zhǔn)。然后按反方位角排列,突出顯示2.3 s～4.3 s的波形,可以明顯看到Pms震相存在到時(shí)差(圖7)。

圖7 余干臺(tái)Ps震相Fig.7 Ps phase of YUG station

圖8 余干臺(tái)Ps震相對(duì)齊時(shí)差曲線Fig.8 Ps phase alignment time difference curve of YUG station

利用波形互相關(guān)方法將余干臺(tái)20道接收函數(shù)Ps震相對(duì)齊(圖8),得到不同方位角的Ps震相相對(duì)走時(shí)差,然后進(jìn)行諧波分解,4階諧波分解結(jié)果見表5?；诒?的二階諧波分析結(jié)果,可以認(rèn)為地殼各向異性的時(shí)間延遲約為0.28 s,快波方向?yàn)?3°。但是諧波分解結(jié)果顯示,1階和4階的振幅系數(shù)也較大,這表明臺(tái)站下方的結(jié)構(gòu)較為復(fù)雜,僅僅使用水平軸的各向異性地層模型很難解釋,可能要考慮傾斜界面甚至其它更復(fù)雜的地殼結(jié)構(gòu)模型才可以完全解釋不同方位角的Ps震相時(shí)間偏差。

表5 余干臺(tái)4階諧波分解結(jié)果Table 5 The fourth order harmonic decomposition results of Yugan station

3 結(jié)論和討論

本文基于各向異性介質(zhì)下Pms震相到時(shí)的理論公式和理論地層模型,驗(yàn)證了諧波分解方法求取地殼各向異性參數(shù)的可行性和可靠性。對(duì)江西余干地震臺(tái)實(shí)際資料處理,認(rèn)為余干臺(tái)地殼各向異性的快波方向約為83°,時(shí)間延遲約為0.28 s。在處理中發(fā)現(xiàn),Pms相對(duì)到時(shí)的拾取精度對(duì)諧波分解結(jié)果的影響是很大的。本文通過波形互相關(guān)方法對(duì)齊Pms震相,求取相對(duì)時(shí)差,發(fā)現(xiàn)并不是每個(gè)反方位角Pms震相的波峰都可以對(duì)齊。因?yàn)椴ㄐ位ハ嚓P(guān)計(jì)算的是給定時(shí)間區(qū)間內(nèi)波形的整體相似度。韓明等[5]提出可以通過提取不同方位角Pms震相波峰最大值的時(shí)間,然后進(jìn)行諧波分解,本文進(jìn)行了嘗試,但結(jié)果相差不大。余干臺(tái)接收函數(shù)諧波分解的結(jié)果顯示,2階諧波并不能完全擬合Pms震相的相對(duì)時(shí)差,表明臺(tái)站下方還存在傾斜等其它橫向不均勻性。王瓊等[12]研究認(rèn)為傾斜界面的存在不會(huì)對(duì)快波方向產(chǎn)生影響,但是會(huì)影響延遲時(shí)間的計(jì)算的。所以本文得到快波方向參數(shù)比延遲時(shí)間具有更高的準(zhǔn)確性。

楊曉瑜[13]使用SKS分裂方法研究了長(zhǎng)江中下游地區(qū)地幔的各向異性,結(jié)論為快波方向總體呈現(xiàn)E-W向,平均快波方向?yàn)?6.19°,和本文得到的余干臺(tái)地殼各向異性快波方向較為一致,支持華南地區(qū)殼幔耦合變形的觀點(diǎn)。石玉燕等[14]對(duì)華東地區(qū)部分臺(tái)站進(jìn)行了SKS分裂研究,包括余干地震臺(tái)。文中的得到的余干地震臺(tái)快波方向?yàn)?02°,時(shí)間延遲為0.9 s?？觳ǚ较蚝捅疚牡慕Y(jié)果存在一定的差別,需進(jìn)一步分析。本文得到的時(shí)間延遲小于0.3 s,結(jié)合前人得到的SKS時(shí)間延遲,支持地幔的各向異性強(qiáng)于地殼的觀點(diǎn)。與前人結(jié)果的對(duì)比進(jìn)一步證明了本文方法和結(jié)論的可靠性,未來將該方法進(jìn)一步用于江西地區(qū)的其它臺(tái)站研究,為理解該區(qū)的形變機(jī)制提供更多的約束。

致謝：本文使用的諧波分解基礎(chǔ)函數(shù)為harmfit,來 自GITHUB,num-matlab-master軟 件包,作者為Fran?ois Beauducel。繪圖過程中使用了EikeRietsch的Seislab工具箱和GMT軟件,審稿專家對(duì)文章的修改提出了建設(shè)性意見,這里一并表示感謝。