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

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

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

0 引言

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

1 方法及理論測試

1.1 方法原理

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

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

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

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

1.2.1基于公式的理論測試

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

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

圖1 YUG地震臺位置圖Fig.1 Location map of YUG station

圖2 理論和諧波擬合的Pms震相走時圖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)式對照認(rèn)為,harmfit諧波分解結(jié)果和假設(shè)的參數(shù)是一致的。通過該理論公式和harmfit的函數(shù)擬合,驗證了harmfit函數(shù)諧波分解的正確性,并知道如何理解harmfit的諧波分解參數(shù)。

1.2.2 基于模型的理論測試

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

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

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

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

圖4 A1模型Ps震相對齊時差曲線Fig.4 Ps phase alignment time difference curve of A1 model

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

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

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

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

2 實際臺站處理

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

2.1 方位角完整性分析

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

表4 YUG臺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ù)的形態(tài),在所有接收函數(shù)中挑選出Ps震相較為清晰,波形質(zhì)量較好的接收函數(shù),然后對于某一個方位角范圍,在挑選好的多個接收函數(shù)中挑選出一條高質(zhì)量的可代表該方位角接收函數(shù)特征的事件參與不同方位角的接收函數(shù)反演。而對于接收函數(shù)很少的方位角區(qū)間,逐條挑選,確保方位角的覆蓋能達(dá)到最大。

2.2 實際臺站處理

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

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

圖7 余干臺Ps震相Fig.7 Ps phase of YUG station

圖8 余干臺Ps震相對齊時差曲線Fig.8 Ps phase alignment time difference curve of YUG station

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

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

3 結(jié)論和討論

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

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

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