基于IGS觀測數據的智利Mw8.2級地震分析

程喆,余晨(1.中國地質大學信息工程學院,湖北武漢 430074; .中國人民解放軍61175部隊,湖北武漢 430074)

基于IGS觀測數據的智利Mw8.2級地震分析

程喆1,2?,余晨2
(1.中國地質大學信息工程學院,湖北武漢 430074; 2.中國人民解放軍61175部隊,湖北武漢 430074)

摘 要:利用2014年智利Mw8.2級地震前后各3d的震中周邊40個IGS站的連續觀測數據和IGS精密星歷,計算獲得遠場同震地表位移,結果顯示,智利地震同震水平位移影響范圍達1 500 km～2 000 km;經過對其中2個IGS站的高頻觀測數據進行解算,獲得地震發生瞬時測站動態運動軌跡,在此基礎上,使用S變換,過濾背景噪聲,有效拾取P波初至時刻,反演出P波在地殼中的傳播速度,并與IRIS提供的參考值進行對比和分析,說明結果的可靠性。

關鍵詞:智利地震;同震位移;高頻數據;S變換;P波初至

1　引 言

智利時間2014年4月1日20時46分于智利伊基克西北約99 km處海域(19.610°S,70.769°W)發生Mw8.2級地震,震源深度25 km。該地震發生在納斯卡板塊和南美板塊交匯處,在該區域納斯卡板塊向東俯沖到南美板塊下部,兩者相對匯聚速率達到7.2 cm/ a(http:/ / earthquake.usgs.gov)。

由于震源遠離陸地,板塊破裂情況無法直接通過野外地質考察或航攝影像得出,只能通過余震定位、地震波震相識別和GPS分析進行。王敏等利用全球GPS觀測資料,計算出2004年印尼蘇門答臘島Mw9.3級地震的同震永久位移影響范圍為6 000 km ～7 000 km[1]。“中國地殼運動觀測網絡”項目組利用GPS獲取的地表形變資料,測定了2008年汶川Ms8.0級地震的同震位移場[2]。楊少敏等利用IGS和陸態網數據計算了2011年日本Mw9.0級地震的遠場同震位移[3,4].繆淼等根據該地區的地質構造、地球物理資料分析了2012年印尼Mw8.6級地震的孕震機理[5],然而目前分析2012年印尼Mw8.6級地震同震位移的資料仍比較少。

近年來,隨著GPS接收機存儲容量的不斷增大,高頻(1 Hz)和超高頻(2 Hz～50 Hz)GPS技術開始出現,GPS作為地震儀的補充,彌補了地震儀速度或加速度積分容易出錯和振幅容易飽和的缺點,Irwan、方榮新、蘇小寧等人利用高頻GPS對2003年Tokachi-oki M8.1級地震,2004年Sumatra-Andaman M9.3級巨震,2008年汶川M8.0級地震,2011年日本宮城M9.0級巨震進行了研究,得出了高精度、高分辨率的地表位移時間序列[6～10],吳繼忠、彭方喜、徐韶光等人基于高頻GPS得出日本宮城地震和汶川地震中S波的平均傳播速度[11～13],彭懋磊等人得出日本宮城地震中Rayleigh波的傳播速度[14],然而S波和Rayleigh波在地震波中破壞力最大,事后拾取其初至時刻,社會效益不大,只有破壞力較小,第一個達到測站的P波的初至時刻,快速分析其頻率、振幅等信息,確定地震影響范圍和隨后而來的S波的達到時刻,才能實現地震預警。張小紅[15]等人利用高頻GPS基于PPP進行了P波初至拾取,取得了比較好的結果,然而目前利用高頻GPS基于TRACK進行P波初至拾取的資料比較少。

2　同震位移

2.1利用GAMIT/ GLOBK處理GPS數據

本文利用GAMIT/ GLOBK軟件對智利Mw8.2級地震前后各3 d(3月30日～4月4日)震中周邊40個IGS站的觀測數據進行統一處理。由于地震發生于UTC時間2014年4月1日23:46:47,因此4月1日當天的數據只采用了UTC時間23:00以前的數據。GAMIT解算設置如表1所示。

GAMIT解算設置 表1

經過解算后,得到各自對應于不同坐標框架的6個單日松弛解,利用GLOBK將這些單日解與斯克里普斯軌道與永久陣列中心(Scripps Orbit and Permanent Array Center,SOPAC)計算出的全球IGS站的6個單日松弛解(igs1～igs6)合并,綜合考慮40個站點的分布情況和數據質量,選擇其中21個作為框架點,采用七參數將單日解轉換到統一的ITRF2008框架下,形成具有統一基準的單日解。單日解絕大多數站點水平精度優于5 mm,垂直精度優于10 mm。

2.2計算與結果分析

在地殼運動分析中,GPS點在ti時刻的位置分量可以描述為[2]:

其中,y(t0)為t0時刻的位置,v表示測站長期運動速度,aj、fj、φj分別為年周期和半年周期變化的振幅、頻率和初相,c為同震位移,在以對數形式的函數代表的震后弛豫模型中,p、τ分別為位移系數和弛豫時間常數,teq為地震發生時刻,ri為觀測噪聲。這里由于參與解算的數據時間跨度比較短(6 d),期間震后的蠕滑形變和周期性變化等非構造形變可以忽略不計,求解同震位移c。解算結果能夠反映地震瞬間的地表位移情況,計算精度取決于單日解的精度。通過計算,得出27個IGS站水平方向的同震位移場,如圖1所示。

圖1　智利Mw8.2級地震引起的遠場水平同震位移場

同震位移計算結果表明,地震的影響范圍只局限在震中區域附近1 500 km～2 000 km范圍內,且隨著震中距的增大同震位移逐漸減小。距離震中最近的IQQE站向北移動1.2 cm,向西移動42.7 cm,AREQ站向北移動0.3 cm,向西移動0.2 cm,UNSA站向北移動0.2 cm,向西移動0.4 cm。

選取2014年2-10月在中山大學中山眼科中心海南省眼科醫院屈光中心接受全飛秒激光手術的近視散光患者。按患者選擇的手術方式將他們分為FLEx組28例（56眼）和SMILE組33例（66眼）。本研究獲得中山眼科中心海南省眼科醫院倫理委員會的批準（批號：2014-005），所有患者術前均簽署手術知情同意書。

3　動態運動軌跡

TRACK是GAMIT/ GLOBK的運動學分析模塊,采用雙差模式進行高精度動態相對定位,逐個歷元進行差分計算,獲得站點的三維坐標和相應誤差,從而得出站點動態運動軌跡。本文選取震中距由小到大的IQQE站(68 km)、AREQ站(460 km)、UNSA站(689 km)、BOGT 站(2 724 km)和ISPA站(4 011 km)1-Hz高頻GPS數據,使用TRACK分析,通過選取合適的固定站[16,17],處理從UTC時間4月1日23:00～4月2日01:00數據,在所得結果中取UTC時間4月1日23:46～4月2日00: 01之間15 min的結果,得到各個測站北、東、垂直方向動態軌跡如圖2所示,測站從上到下按照震中距由小到大排列。可以看出,地震發生后,距離震中最近的IQQE站最先感知,隨著地震波的傳播,其能量逐漸衰弱,各站各分量均有不同程度的變化,變化大小主要取決于震中距及地震帶破裂方向。

圖2已經詳細的反映了測站在N、E、U方向的運動過程,為了更直觀地刻畫測站運動軌跡,圖3中給出了UTC時間4月1日23:46～4月2日00:01之間15 minIQQE和AREQ水平向的運動軌跡,兩站初始位置都為(0,0)。可以看出,地震波到達前,測站坐標非常集中,地震波到達后,能量逐漸聚集,測站發生劇烈抖動,偏離初始位置,地震波衰退后,測站逐漸回彈到初始位置附近,坐標又趨于集中,初始位置和最終位置之間的差異即為該站高頻同震位移,列于表2中。

圖2　同震位移時間序列(15 min)

圖3　IQQE和AREQ水平方向運動軌跡

高頻GPS測站同震位移 表2

從表2可以看出,高頻同震位移方向與基于年積日平均坐標的測站永久變形相同,但大小有差異,說明主要位移是由主震引起,差異由主震誘發的多次余震及震后緩慢蠕動造成。

4　地震波震相識別

4.1S變換原理及初至識別機理

地震波中體波分為P波(縱波)和S波(橫波),它們在巖石圈中的實際傳播速度取決于巖石的密度和內在彈性。總的說來,P波在地殼中運動速度比S波快, P波第一個到達目標測站,使地面發生上下振動,破壞性較弱。S波第二個到達,使地面發生前后、左右抖動,破壞性較強。地震預警的基本思想是利用P波與S波傳播過程中的時間差,當P波到來時,快速準確判斷波形,分析其初至時刻,預估隨后而來的S波初至時刻和影響范圍,確定預警級別和范圍。通常預警時間都只有九秒至數十秒,但這對于減輕人員和財產損失起著十分重要的作用[15]。

S變換(S -Transform)是由美國地球物理學家Stockwell[18]等在1996年提出的一種無損可逆的線性時頻分析方法,它克服了短時傅里葉變換時窗函數固定、對信號時間定位能力差的缺點,同時其基本小波不必滿足容許性條件,是一種優于短時傅里葉變換和小波變換的時間-頻率分析技術,在信號處理[19]、地震分析[20]領域得到了廣泛應用。利用S變換進行時間序列的頻譜分析,一方面能夠提高對背景噪聲成分的時間定位能力,便于實現噪聲的濾波處理;另一方面可以增加頻譜系數的個數,從而改善地震波信號的統計特性[21]。

Stockwell等將時域信號h(t)的S變換定義為

其中,S(τ,f)為變換后得到的時頻譜,ωf(t)為高斯時窗函數,τ表示時窗函數中心對應的時間,f表示頻率(單位Hz)。對于一維函數h(t),變量τ和其中一個頻率fi可以確定一個S(τ,fi),稱其為一個時窗的信號。因此,當頻率取不同值時,S變換可以獲得不同頻率信號頻譜隨時間的變化特征。

對單道地震記錄經S變換后,得到時間、頻率的聯合表示,在地震波到來前,觀測到的一般都是背景噪聲,而地震波到來后,地震信號與背景噪聲混雜在一起,這樣振幅能量的突變點就是初至起跳處。如果信噪比比較高,用肉眼可以在時間域上識別振幅能量在初至起跳處前后的差異;如果信噪比低,初至起跳處前后的能量差異小,在時間域無法識別,可以利用地震波與背景噪聲頻率的差異,使得初至前后的能量在時頻圖上具有不同的空間分布,同時地震信號的時頻譜有明顯的能量聚集的特性加以區分,有效拾取P波初至時間。

4.2提取背景噪聲

首先取IQQE和AREQ兩個測站震前15 min(900歷元)的高頻時間序列進行一階求差,消除長周期噪聲,獲得測站速度波形,對速度波形進行S變換,將時間域信息轉換到頻率域,判斷背景噪聲能量大小,從而確定提取P波初至時刻的能量閾值,得到如圖4所示頻譜圖,第一行為IQQE站,第二行為AREQ站,從左至右依次為北、東、垂直方向,橫軸為歷元數,縱軸為頻率。

圖4　IQQE站和AREQ站震前速度波形S變換

可以看出,震前速度波形背景噪聲頻率主要分布在0.2 Hz～0.5 Hz之間,十分雜亂,沒有規律,垂直向背景噪聲能量明顯大于水平向,這是由于GPS垂直方向比水平方向觀測誤差大造成的,北、東、垂向背景噪聲能量分別為1.5×10-4、1×10-4和2×10-4,這些數值將作為P波初至拾取的能量閾值,低于或等于該值的能量信號將在后續分析中予以過濾去除。

4.3拾取P波初至時刻

分別取IQQE和AREQ站UTC時間23:46～23: 51各5 min的高頻時間序列,經過S變換后過濾背景噪聲,獲得頻譜圖和速度序列圖如圖5、圖6所示,其中從上至下是頻譜圖和速度圖,從左至右依次是北、東、垂向,橫軸是UTC時間,縱軸是頻率或速度,頻譜圖中能量低于或等于背景噪聲的部分用白色表示。

可以看出,當P波剛剛到達時,能量比較弱,速度序列圖變化不明顯,信噪比較低,肉眼很難根據速度變化識別出初至時刻,隨著破壞力較大的S波和面波到達,能量逐漸增強到最強(對應于速度序列圖變化最快的部分),維持一段時間后,隨著地震波的減弱而逐漸衰退,最后恢復到背景噪聲。拾取的P波初至時刻列于表3中,同一測站由N、E、U三個方向拾取的時刻并不相同,差異在10 s之內。

圖5　IQQE站震時頻譜圖

圖6　AREQ站震時頻譜圖

各測站P波初至時刻對比 表3

4.4初至時刻檢驗

為了驗證S變換提取P波初至時刻的準確性,將IRIS公布的附近地震臺站地震儀觀測P波初至時刻與本文進行對比分析,結果列于表3中,其中臺站LVC距離IQQE站345 km,臺站NNA距離AREQ站 700 km,根據臺站到震中的距離和IRIS公布的P波初至時刻,計算出P波在地殼中傳播距離,分別用LVC、NNA監測到的P波初至時刻內插出P波傳播到IQQE和AREQ站的時刻,與本文計算結果進行對比。

從表3可知,LVC內插得出IQQE站P波初至時刻參考值為23:46:56,NNA內插得出AREQ站P波初至時刻參考值為23:47:46,測站S變換的拾取結果都要比IRIS提供的參考值晚2 s～5 s,這是因為首先,受精度限制,GPS本身測量噪聲較大,在測站較遠,信噪比較低的情況下,不容易將背景噪聲與地震信號十分準確的分離出來,而地震儀的信噪比比高頻GPS要高;其次,S變換中背景噪聲能量閾值主要是通過目視頻譜確定,難免存在一定誤差,而閾值的選取對于拾取時刻的影響可能達到數秒;第三,N、E、U三方向拾取時刻的不同是由各個方向對地震波的敏感程度不同引起的,不能簡單求平均,確定整體初至時刻;第四,地震儀兩測站到高頻GPS兩測站有一定距離,期間由于巖石的密度和內在彈性的不同,P波傳播速度會發生一定變化,簡單的內插與外推,會有誤差存在。使用S變換拾取P波初至時刻與地震儀數據不完全一致,結果大體上吻合較好。

另一方面,由S變換反演出P波傳播速度為6.80 km/ s和7.30 km/ s,這與P波在地殼中平均傳播速度為5 km/ s～8 km/ s的范圍一致。

5　結 語

本文利用智利Mw8.2級地震前后各3 d的震中周邊40個IGS站的連續觀測數據和IGS精密星歷,計算獲得遠場同震位移,結果顯示其同震水平位移影響范圍達1 500 km～2 000 km;對其中2個IGS站的高頻數據進行解算,獲得地震發生瞬時二維和三維動態運動軌跡以及高頻同震位移,說明高頻同震位移與基于年積日平均坐標的測站永久變形的差異由主震誘發的多次余震及震后緩慢蠕動造成。

在此基礎上,對測站速度波形使用S變換,過濾背景噪聲,在北、東、垂直三個方向有效拾取P波初至時刻,并與IRIS提供的參考值進行對比和分析,發現同一測站三個方向S變換拾取時刻各不相同,差異在5 s之內;兩個測站拾取時刻都比參考值要晚2 s～5 s,總體上吻合較好,反演出P波在地殼中的傳播速度在6 km/ s～8 km/ s,與合理值較符合,說明了結果的可靠性,研究成果對于今后利用GPS數據進行地震預警有一定參考意義。

參考文獻

[1] 王敏,張培震,沈正康等.全球定位系統(GPS)測定的印尼蘇門達臘巨震的遠場同震地表位移[J].科學通報, 2006,51(3):365～368.

[2] 國家重大科學工程“中國地殼運動觀測網絡”項目組.GPS測定的2008年汶川Ms8.0級地震的同震位移場[J].中國科學D輯:地球科學,2008,38(10):1195～1206.

[3] 楊少敏,聶兆生,賈志革等.GPS解算的日本Mw9.0級地震的遠場同震地表位移[J].武漢大學學報信息科學版, 2011,36(11):1336～1339.

[4] 王敏,李強,王凡等.全球定位系統測定的2011年日本宮城Mw9.0級地震遠場同震位移[J].科學通報,2011, 56(20):1593～1596.

[5] 繆淼,朱守彪.2012年北蘇門答臘西海域Mw8.6地震的孕震機理及其對地震活動性的影響[J].地震工程學報, 2013,35(2):278～288.

[6] Irwan M,Kimata F,Hirahara K,et al.Measuring ground deformations with 1-Hz GPS data:the 2003 Tokachi-oki earthquake(preliminary report).Earth,Planets and Space,2004, 56(3):389～393.

[7] Ohta Y,Meilano I,Sagiya T,et al.Large surface wave of the 2004 Sumatra - Andaman earthquake captured by the very long baseline kinematic analysis of 1-Hz GPS data.Earth, Planets and Space,2006,58(2):153～157.

[8] 方榮新,施闖,王廣興等.利用高頻GPS確定的大地震震中和震級研究:2008年汶川8.0級地震應用結果[J].中國科學D輯:地球科學,2014,44(1):90～97.

[9] 蘇小寧,孟國杰,廖華等.利用高頻GPS數據研究汶川地震引起的近場地表運動學特征[J].大地測量與地球動力學,2013,33(SI):11～15.

[10] 黃勇,楊少敏,劉剛等.2011年日本Mw9.0地震的高頻GPS動態形變研究[J].大地測量與地球動力學,2012, 32(6):33～36.

[11] 吳繼忠,吳文壇.利用高頻GPS進行地震動態變形分析及地震定位[J].大地測量與地球動力學,2012,32(2): 20～23.

[12] 彭方喜,吳云.高頻GPS動態監測地表形變的試驗與研究[J].大地測量與地球動力學,2012,32(S):109～112.

[13] 徐韶光,熊永良,廖華等.汶川地震同震形變的靜態和動態分析[J].大地測量與地球動力學,2010,30(3):27 ～31.

[14] 彭懋磊,劉剛,聶兆生等.日本Mw9.0強震的高頻GPS形變波特征初探[J].大地測量與地球動力學,2011,31 (6):6～10.

[15] 張小紅,郭斐,郭博峰等.利用高頻GPS進行地表同震位移監測及震相識別[J].地球物理學報,2012,55(6): 1912～1918.

[16] 范士杰,劉焱雄,張健.PPP與GAMIT/ TRACK在地震監測中的應用[J].測繪科學,2013,38(2):184～186.

[17] 殷海濤,肖根如,張磊等.TRACK高頻GPS定位中震時參考站的選取方法[J].大地測量與地球動力學,2012, 32(4):15～19.

[18] Stockwell R G,Mansinha L,Lowe R P.Localization of the complex spectrum:The S transform.IEEE Transactions on Signal Processing,1996,44(4):998～1001.

[19] 景建恩,魏文博,陳海燕等.基于廣義S變換的大地電磁測深數據處理[J].地球物理學報,2012,55(12):4015 ～4022.

[20] 陳昕,李新安,劉威序等.基于S變換的初至拾取方法及應用[J].中國煤炭地質,2009,21(11):59～63.

[21] 劉琦,張晶.S變換在汶川地震前后應變變化分析中的應用[J].大地測量與地球動力學,2011,31(4):6～9.

Analysis of Chile Mw8. 2 Earthquake Based on Observations of IGS Stations

Cheng Zhe1,2,Yu Chen2
(1.Faculty of Information Engineering,China University of Geosciences,Wuhan 430074,China; 2.PLA 61175,Wuhan 430074,China)

Abstract:Based on IGS precise ephemeris and 3-day IGS observations before and after the great Mw8.2 Chile earthquake,we presented the coseismic deformation field composed of 40 sites for Indonesia and surrounding area.The results show that the scope of co-seismic horizontal displacement up to 1500-2000km;we solved high-frequency observations of 2 IGS stations to obtain the dynamic trajectory,furthermore,we used the S-transform algorithm for background noise Filtering and the primary wave(P-wave) arrival picking from high-rate GPS seismic waves,inversed the P-wave propagation velocity in the earth's crust,we compared and analyzed it with reference values provided by IRIS(Incorporated Research Institutions for Seismology) to illustrate the reliability of the result.The results of research could be used for GNSS data monitoring and forecasting on future as a reference.

Key words:chile earthquake;coseismic displacement;high-frequency data;s-transform;p-wave arrival time

文章編號:1672-8262(2015)01-73-06中圖分類號:P228

文獻標識碼:A

收稿日期:?2014—09—01

作者簡介:程喆(1989—),男,碩士,助理工程師,從事GPS數據處理科研與生產工作。