利用高頻GPS數據反演地震要素

張環宇,林吉航,陳庭,2

(1.武漢大學 測繪學院,武漢 430072;2.地球空間環境與大地測量教育部重點實驗室,武漢 430072)

0 引 言

在地震發生后,快速、準確地反演出地震的相關要素是地震預警工作開展的前提. 由于高頻GPS技術相較傳統技術具有諸多優點,例如能夠直接解算得到瞬時的地表三維位移、在反演震級時不存在飽和現象、誤差不隨時間積累等[1],其在地震學中的應用越來越廣泛.

目前利用高頻GPS觀測數據反演地震要素的簡要流程是:首先通過精密單點定位(PPP)或差分相對定位對高頻GPS觀測數據進行動態解算;然后使用S變換或恒星日濾波等方法對解算結果進行去噪處理;接著使用長短時窗法或閾值探測法等探測地震波到時,再利用人工識別法識別地震波振幅與周期;最后使用網格搜索法或最小二乘法等求解震中位置以及地震波傳播速度,使用經驗公式或有限斷層反演法等反演震級,再使用平均法反演地震發震時刻. 根據測站位置分布以及震動幅度,還可初步判斷地震斷層類型以及斷裂帶的走向等[2-6].

本文系統地闡述了高頻GPS數據處理以及地震要素反演的流程和相關理論,完整地實現了這一流程,提出了拾取地震中體波和面波到達時刻的雙重3-sigma法則,并設計了自動提取地震波振幅與周期的算法. 最后以2019年7月美國加州里奇克萊斯特群震中震級較大的兩場地震作為對象進行研究.

將本文的反演結果與美國地質勘探局(USGS)發布的結果進行比對,發現二者較為接近,可考慮將本文描述的方法應用于地震要素反演與預警工作之中.

1 高頻GPS觀測數據處理

1.1 數據獲取

本文使用的高頻GPS觀測數據來源于美國衛星導航系統與地殼形變觀測研究大學聯合體(UNAVCO)布設的GPS監測網絡,選用數據的觀測頻率為1 Hz,下載地址為ftp://data-out.unavco.org/pub/.

1.2 數據解算策略

PPP與差分相對定位是目前高頻GPS觀測動態解算的兩種常用方法. 由于差分相對定位需要選取1個或多個參考站. 當離震中位置較近時,參考站亦會受到地震影響;但若離震中位置過遠,基線的解算精度又會降低. 相對而言,PPP無需參考站即可直接定位,較為靈活. 因此,本文研究決定使用PPP方法解算高頻GPS觀測數據,具體的解算工作通過武漢大學衛星導航定位技術研究中心研制的Pride PPP-AR軟件[7]實現.

1.3 解算結果

最初選取了分布在震中附近的30個測站. 首先通過軟件解算得到相應的位移時間序列,觀察發現部分測站的位移時間序列在地震前后沒有明顯的起伏,推測為這些測站受到的地震影響較小,在后續的處理中沒有考慮這些震動不明顯的測站.

最終在研究2019年7月4日的地震時采用cccc、p556、p560、p568等一共8個測站;而在研究2019年7月6日的地震時采用bsry、cccc、p056、p463等共21個測站,如圖1所示.

注:三角形測站在兩次地震中均有明顯位移,菱形測站僅在7月4日地震中有明顯位移,圓形測站僅在7月6日地震中有明顯位移. 虛線代表7月4日地震中有明顯位移的測站的大致分布走向.圖1 測站分布圖

繪制選取測站對應的位移時間序列:

圖2 2019年7月4日地震前后各測站三方向位移

圖3 2019年7月6日地震前后各測站三方向位移

從圖2和圖3可以看出,解算結果的平面精度明顯優于高程精度,在后續解算過程中僅利用了平面上的位移時間序列.

1.4 位移時間序列去噪

上述解算后得到的位移時間序列中包含一定程度的噪聲,消除或減弱噪聲的影響,能提高后續工作的準確性和效率.

1.4.1S正變換

S變換是一種結合了小波變換和短時傅里葉變換的可逆時頻域定位技術[8],具有較高的時頻分辨率,被廣泛應用于信號處理領域. 對于一個時間序列函數h(t),其S變換的數學公式為[9-11]

(1)

式中：τ為控制高斯窗口在時間軸上移動的參數；f為頻率.

將S變換譜隨時間進行平均,可將其轉換為傅里葉變換譜

(2)

式中，H(f)是時間序列h(t)的傅里葉變換.

在實際處理時,需要采用離散形式進行時頻域變換.時間序列函數h(t)的離散形式表示為h[kt],k=0,1,…,N-1. 其中T為采樣時間間隔.

傅里葉變換的離散形式為[9-10]

(3)

式中，n=0,1,…,N-1.

當n≠0時,S正變換的離散形式為

(4)

當n=0時

(5)

利用S正變換可以將原始的時域信號進行二維時頻域平面分解,一些原本在時域難以區分的信號(例如體波和背景噪聲)由于其頻率不同,通過S正變換后可以較好地區分開來.

1.4.2 閾值函數選取

頻率域常用的濾波函數有硬閾值濾波和軟閾值濾波等. 硬閾值濾波的表達式為

(6)

軟閾值濾波的表達式為

(7)

式中:x表示頻域信號,取絕對值表示信號幅值;τ即為閾值,可由下式[12]求得.

(8)

式中：Med(·)表示求取數據對應的中位數,Med(|x|)/0.6745是對信號噪聲標準差的估計；N為信號長度.

比較兩個閾值函數的表達式:硬閾值濾波不會削弱信號的幅值,但會使信號在頻率域發生突變,導致去噪后的結果出現局部抖動;軟閾值函數采用特殊的處理方式避免了“一刀切”的現象,使頻率域內的信號更平滑,但會削弱信號的幅值,造成一定的高頻信息損失,導致信號在一定程度上失真[3]. 綜合考慮,本文采用折中閾值濾波函數在頻域對信號進行濾波處理,即

(9)

式中：α為折中系數,取值在[0,1]之間. 當α=0時,該函數退化為硬閾值濾波;當α=1時,該函數退化為軟閾值濾波.

折中閾值濾波可以在保證信號平滑的同時,盡可能地減小信號高頻信息的損失. 本文選取不同的α進行實驗后,發現當α接近于1時,信號幅值會受到較大的削弱;當α接近于0時,濾波后反變換回時間域的結果存在一些局部抖動現象；當α=0時局部抖動尤其明顯.結合地震要素反演工作的內容,我們更加關心的是地震波到達時測站的位移幅值,較為緩和的局部抖動現象可以通過一定的手段避免其影響.經過試驗發現當α=0.1時效果最好,本文在計算過程中也選用了該值.

1.4.3S逆變換

連續形式的S逆變換公式如下[9-10]:

(10)

類似地,在實際處理時需要使用離散形式的逆變換,公式如下[9-10]:

(11)

通過S逆變換可將濾波后的信號重新轉換到時域.

1.4.4 去噪結果

限于篇幅本文以2019年7月4日地震前后cccc測站北方向的位移時間序列為例進行說明.

繪制cccc測站濾波前后的位移時間序列與S變換譜,如圖4所示. 觀察圖4(a)和圖4(c)可知當地震波到達測站后,在二維時頻域平面上的S變換譜的能量值有明顯的突變,同時地震波信號與背景噪聲被很好地分離開來.

S變換譜經過折中閾值函數濾波處理之后,再反變換回時間域的效果如圖4(b)和圖4(d)所示. 可以看出圖4(c)變換譜中的背景噪聲被有效去除,同時能量幅值基本保持不變;經過S逆變換返回時間域后,位移時間序列變得更加平滑.

圖4 cccc測站北方向位移時間序列與S變換譜濾波前后對比(2019年7月4日，部分)

2 地震要素反演

通過分析去噪后的高頻GPS位移時間序列,可以拾取地震波到達測站時刻,并反演出地震三要素(震中、震級、發震時刻)等相關信息.

2.1 拾取地震波到時

地震波是從震源產生的向四周輻射的彈性波. 按傳播方式可分為體波和面波兩種類型. 體波具體包括縱波(P波)和橫波(S波). 縱波的傳播速度大致為5.5～7 km/s,橫波的傳播速度大致為3～4 km/s. 地震發生后,縱波、橫波與面波先后到達測站.

由于體波振幅較小,往往無法通過GPS觀測獲得;而面波能量強、振幅較大,所以通過GPS觀測到的地表起伏往往是由面波引起的[3].

首先分別將原始的位移時間序列和經過濾波后的位移時間序列進行一階差分獲得對應的速度時間序列,同樣以2019年7月4日cccc測站北方向位移為例.

圖5 cccc測站北方向位移與速度時間序列濾波前后對比(2019年7月4日，部分)

原始的位移時間序列在一階差分后,基本消除了長周期低頻噪聲的影響,明顯更加穩定. 而經過S變換濾波后的位移時間序列在一階差分后效果更佳.

由于體波的振幅較小,往往被背景噪聲淹沒,所以從原始的位移時間序列中僅能提取出面波到時. 利用S變換和折中閾值濾波對原始位移時間序列進行處理之后,可以較好地濾除背景噪聲,提取出體波. 同時因為在1.4.2節中為了盡可能保證振幅不受干擾,選取的閾值函數折中系數較小,帶來了一些局部抖動. 所以本文研究設計了雙重3-sigma法則以先后拾取面波到時和體波到時,同時避免了由于局部抖動而對地震波到達測站的時刻產生誤判.

具體的拾取流程如下:

1) 針對原始的速度時間序列,使用3-sigma法則判斷拾取面波到時:若某一時刻的速度值大于之前300 s內速度值標準差的3倍,則認為在此時刻面波到達,記為T0.

2) 對于經過S變換和折中閾值濾波平滑后的速度時間序列建立搜索區間[T0-30,T0](單位:s. 當震中距在200 km以內時,體波與面波到時的理想差值大致在30 s之內),在該區間內再次使用3-sigma法則搜索體波到時,記為T1.

最終得到面波和體波到達各個測站的時刻如圖6所示.

圖6 7月4日地震中各測站體波與面波到時

注:圖6、圖7中的震中距通過USGS網站公布的地震震中位置和測站位置計算.圖7 7月6日地震中各測站體波與面波到時

從圖6和圖7可以看出地震波到達時刻與震中距之間大致呈線性正相關的關系,與實際情況相符.對地震波到時和震中距進行線性擬合,所得直線斜率的倒數即為地震波傳播速度. 7月4日地震中體波波速為3.974 km/s,面波波速為2.761 km/s;而7月6日地震中體波和面波的傳播速度則分別為3.346 km/s和2.801 km/s.

根據擬合直線斜率的倒數估算的體波與面波波速與理論上的傳播速度存在一定的差距,原因可能是地震波實際的傳播路徑較為復雜,而本文采取的震中距為地表距離,沒有考慮深度信息等.

2.2 震中位置反演

在拾取體波或者面波到達各個測站的時刻之后,便可以根據各測站的實際位置坐標來反演震中位置和地震波傳播速度.

本文采用格網搜索法[2]反演震中位置. 首先以測站網重心為中點,建立4°×4°的大范圍搜索區域,同時將格網分辨率設置為0.05°×0.05°.

然后按照行列順序依次假定各個小格網的中心為震中,計算各個測站到該格網的距離. 為了簡化計算模型,同時由于測站幾乎分布在同一平面內導致三維幾何結構強度較低[3],本文僅在地球表面上進行搜索而忽略深度的影響. 已知N個測站的經緯度坐標分別為(B1,L1),(B2,L2),…,(BN,LN),當前格網坐標為(B,L),則第i個測站到格網的距離di的表達式如下:

di=arccos(sinBisinB+

cosBicosBcos(Li-L))R.

(12)

其中,R為地球半徑.

接著計算地震波到時殘差絕對值之和為

(13)

其中地震波到時t為體波或面波到時. 體波到時等效于P波到時,所以若t為體波到時,速度v在5.5～7 km/s進行搜索;面波到時等于或者略晚于S波到時,所以若t為面波到時,速度v在3～4 km/s進行搜索. 當殘差絕對值最小時,對應的格網點坐標即為反演的震中位置.

通過搜索,得到7月4日和6日地震中體波傳播速度均為5.5 km/s,面波傳播速度均為3 km/s. 反演的震中位置與USGS網站公布的震中位置如表1所示:

表1 震中的反演位置與實際位置

反演的震中位置與參考的震中位置較為接近,但仍有千米級的偏離. 這一偏差可能是在反演震中位置時假設地震波沿各個方向的傳播速度相等時引入的. 在實際情況中由于地殼內部介質復雜且分布不均,地震波在各個方向上的傳播速度并不相等[3].

2.3 震級反演

本文采用了四個震級經驗公式反演震級并進行對比.

第一個經驗公式是在1967年國際地震學與地球內部物理學協會(IASPEI)大會上推薦的莫斯科-布拉格經驗公式(IASPEI公式)[4],該公式亦使用在國家標準《地震震級的規定》(GB 17740-2017)中:

M=lg(A/T)+1.66lg(Δ)+3.3.

(14)

式中：A為地面峰值位移，μm；T為地震波周期，s；Δ為震中距，(°).

采用的第二個震級經驗公式是古登堡面波震級公式[13]:

M=lg(A/T)+1.656lg(Δ)+1.818.

(15)

第三個震級經驗公式是Crowell等在2013年根據發生在加州和日本的數次地震的統計結果得到的經驗公式[2,14]:

(16)

第四個震級經驗公式是由Melgar等在2015年根據往期地震數據回歸分析得到的[2,15]:

(17)

值得注意的是不同于式(14)和(15),式(16)和(17)中地面峰值位移A的單位為cm,震中距Δ的單位為km,并且公式中沒有考慮地震波的周期.

地震波在某一方向上的振幅及其對應周期如圖8所示.

圖8 地震波振幅與周期示意圖

本文設計了一套自動提取地震波在東西/南北方向上振幅與周期的算法,流程如下:

1) 在地震波到達測站后,在一段時間內(如60 s)拾取一系列的東西/南北方向上位移極大值點與極小值點,記錄下各個極值點對應的時刻與位移值.

2) 對相鄰的極值點的位移值作差,最大的差值絕對值即代表振幅的2倍,對應的時刻差即為地震波周期的一半.

(18)

根據各測站拾取的地面峰值位移、地震波周期和震中距,利用經驗公式便可以計算出震級.

由于根據體波到時和面波到時搜索得到的震中位置存在略微的差別導致震中距不同,二者反演的震級存在細微的差別,本文決定取二者的平均值進行研究.

最終各個測站反演的震級如圖9所示.

圖9 7月4日各測站反演震級

圖10 7月6日各測站反演震級

取各測站反演震級的平均值作為最終的震級,7月4日使用各經驗公式反演的震級如表2所示.

表2 2019年7月4日地震反演震級對比

USGS網站上公布的震級為6.40級,數據來源于加州綜合地震網絡(California Integrated Seismic Network)布設的強震儀加速度記錄和測站臺速度記錄[16].

本文根據高頻GPS定位結果,利用IASPEI經驗公式反演的震級為6.48級,與之十分接近.

7月6日反演的震級如表3所示.

表3 2019年7月6日地震預測震級對比

7月6日地震中根據高頻GPS定位結果利用四個經驗公式反演得到的平均震級相較于USGS公布的7.10級而言都偏小,但利用IASPEI經驗公式反演得到的6.87級依然與加州綜合地震網絡測定的震級最為接近.

從表2和表3可以看出對于加州里奇克萊斯特地區,IASPEI經驗公式反演的震級效果最好.

2.4 發震時刻反演

根據每一個測站的體波或面波到時ti和震中距di以及地震波傳播速度v,可以計算出地震發震時刻Ti=ti-di/v.取各個測站推算出的發震時刻的平均值,作為最終的地震發震時刻:

(19)

最終反演的發震時刻如表4所示.

表4 地震發震時刻

反演的發震時刻與USGS測定時刻的偏差在20 s以內. 發震時刻反演偏差是由震中位置反演偏差、假設地震波傳播速度恒定等原因共同引起的.

2.5 斷層類型及斷裂帶走向

根據解算的地表位移時間序列可以看出,地震發生后測站主要在平面上發生明顯位移,在垂直方向上的位移不明顯,所以可據此大致推斷本次地震斷層類型為走滑斷層,這也與USGS的結果相同.

根據選取的測站位置分布情況可以看出,在2019年7月4日的地震中有較明顯位移的測站基本都分布在西南—東北走向上,因此推測該方向可能為斷裂帶走向,而在7月6日的地震中位移較明顯測站的分布位置較為均勻. USGS網站上公布的有關信息為斷裂可能是右旋滑動沿西北—東南走向的結果,也可能是左旋滑動沿西南—東北走向的結果. 所以可根據位移明顯的測站分布情況對地震破裂的走向進行一個初步判斷.

3 結論與展望

本文系統闡述了利用高頻GPS數據進行地震要素反演工作的原理,并且完整實現了整套流程,最終得到的結論如下:

1) 僅利用高頻GPS數據可以較為準確地反演地震的有關要素,并且結果具有很好的穩定性；

2) 本文在反演震中位置時假設地震波沿各個方向的傳播速度相等. 然而由于地殼內部介質復雜且分布不均,地震波在各個方向上的實際傳播速度是不相等的,這也是目前根據地震波反演震中和發震時刻的主要誤差來源；

3) 根據高頻GPS定位結果可以對地震的斷層方式和斷裂帶方向進行一個大致判斷,實際應用中應結合地震波數據等共同進行推算；

4) 高頻GPS定位精度受到鐘差、測量噪聲、多路徑等多種誤差的影響,高程方向的精度約為3～5 cm,水平方向的精度約為1～2 cm,這些誤差會對地震要素的反演產生一定的影響.

對于本文的后續工作,我們還進行了如下展望:

1) 震級經驗公式往往具有地域性. 在實際應用中,可以收集某一地區往期的地震數據,通過最小二乘、回歸分析等方法修改震級經驗公式中的有關常數,建立最適用于該地區的震級公式；

2) 高頻GPS的觀測精度有待進一步提高,同時解算模型和誤差模型也有待進一步優化[17-18]；

3) 考慮多源數據如高頻GPS觀測數據、地震臺速度數據、強震儀加速度數據等的融合互補,揚長避短,提升地震要素反演的準確度.

致謝:感謝為本文研究工作提供高頻GPS觀測數據的UNAVCO機構;感謝提供精密單點定位解算軟件服務的武漢大學衛星導航定位技術研究中心Pride Lab團隊.