林清瑩, 于忠海, 張永明, 陳 鵬, 張 碩

(1. 濟南市勘察測繪研究院, 濟南 250101；2. 西安科技大學測繪科學與技術學院, 西安 710054；3. 山東農業(yè)大學資源與環(huán)境學院, 山東 泰安 271018)

0 引言

地震是由于地殼快速釋放能量所造成的震動，對地表具有很強的破壞性，給人類造成嚴重的危害，因此探測地震引起的異?，F(xiàn)象具有重要的意義。地震發(fā)生時釋放的能量不僅會破壞地球表面，還會傳播至電離層并引起電子密度(Ne)和總電子含量(Total Electron Content，TEC)的變化等異?，F(xiàn)象。自1964年美國阿拉斯加大地震中首次探測到地震電離層異常擾動以來，電離層異常與地震之間的耦合關系已經引起了廣泛關注[1]。隨著大量全球和區(qū)域衛(wèi)星導航系統(tǒng)(Global Navigation Satellite System，GNSS)網的建設和使用，GNSS技術已廣泛地應用于地球動力學研究中，突破了傳統(tǒng)地震監(jiān)測技術的時空限制。

通過全球定位系統(tǒng)(Global Positioning System，GPS)觀測數(shù)據(jù)探測發(fā)現(xiàn)，地震發(fā)生時激起的地震波向上傳播至電離層高度，會引起電離層總電子含量的變化，產生同震電離層擾動(Coseismic Ionospheric Disturbances,CID)現(xiàn)象。E.Astafyeva等研究了1994年Kurile Mw8.1級地震同震電離層異?，F(xiàn)象，發(fā)現(xiàn)CID出現(xiàn)的時間與距震中的距離存在相關性[2]。通過CID的傳播速度，可以研究電離層擾動的激發(fā)條件，如陳鵬和陳家君發(fā)現(xiàn)2004年蘇門答臘Mw8.6級和Mw8.2級地震CID傳播速度分別為0.73km/s與0.69km/s，與聲波在0～450km內的平均傳播速度一致[3]。Liu J.Y.等發(fā)現(xiàn)1999年Chi-Chi Mw7.6級地震引起的兩種CID傳播速度分別為2.5km/s和370m/s，認為是由瑞利波和重力波引發(fā)的擾動[4]。CID的傳播方向也是一個重要特征，K.Heki等發(fā)現(xiàn)因受地磁場的影響,北半球的CID表現(xiàn)出向南傳播的趨勢，擾動幅度一般不超過2個 TECU，擾動峰值點出現(xiàn)在震后10～20min，持續(xù)幾分鐘后逐漸消失[5-6]，M.N.Cahyadi等發(fā)現(xiàn)破裂帶走向也會對CID傳播方向產生影響[7-8]。同震地表垂直位移被認為是引發(fā)CID現(xiàn)象的主要原因，如Jin S.等發(fā)現(xiàn)瑞利波引起的CID主要受到垂直地面運動的影響，并會直接影響CID波形[9-10]。E. Astafyeva等發(fā)現(xiàn)CID的幅度和持續(xù)時間與地震的震級存在明顯的正相關關系[11]。此外，低軌衛(wèi)星等空基數(shù)據(jù)也被用于CID研究。Ma X.等利用cosmic數(shù)據(jù)，發(fā)現(xiàn)汶川地震期間TEC和NmF2時空分布變化存在同震響應[12]。

GPS數(shù)據(jù)為探測電離層活動提供了高時空分辨率數(shù)據(jù)支持，這有助于深入探測地震電離層異常效應，促進地震-電離層耦合機制的研究。本文對2018年1月23日在阿拉斯加灣科迪亞克島東南部地區(qū)發(fā)生的Mw7.9級地震進行了研究，據(jù)美國地質調查局(United States Geological Survey，USGS)公布，此次地震是由于太平洋板塊淺層巖石圈斷層走滑造成的，沿西南偏西或西北偏北方向發(fā)生破裂。本文采用美國連續(xù)運行參考站(Continuously Operating Reference Station，CORS)發(fā)布的GPS觀測數(shù)據(jù)，解算孕震區(qū)周圍傾斜總電子含量(Slant Total Electron Content，STEC)數(shù)據(jù),利用奇異譜分析(Singular Spectrum Analysis，SSA)的方法提取電離層擾動信號，對CID現(xiàn)象和較高緯度地區(qū)衛(wèi)星選取進行進一步分析，并利用SWARM衛(wèi)星垂直總電子含量(Vertical Total Electron Content，VTEC)數(shù)據(jù)對震中上空的總電子含量變化進行研究。

1 數(shù)據(jù)和處理方法

1.1 電離層數(shù)據(jù)

本文采用CORS系統(tǒng)發(fā)布的GPS觀測數(shù)據(jù)，反演阿拉斯加地震孕震區(qū)周圍的震發(fā)1h內STEC數(shù)據(jù)，處理后精度達到0.01個TECU[10]，可用于阿拉斯加震后的電離層擾動探測。CORS系統(tǒng)在美國有1900多個觀測站，并且在阿拉斯加地區(qū)分布非常密集，為研究電離層異常擾動提供了充足的數(shù)據(jù)，CORS系統(tǒng)在阿拉斯加地區(qū)的站點分布如圖1所示。利用采樣間隔為30s和15s的GPS數(shù)據(jù)反演的TEC中均發(fā)現(xiàn)了異?，F(xiàn)象，為了提高精度，本文選取了采樣間隔15s的GPS觀測數(shù)據(jù)(ftp://data-out.unavco.org/pub/rinex/obs/)作為實驗數(shù)據(jù)，進行CID探測。

圖1 CORS站點分布、震源球表示震中位置

利用GPS相位觀測數(shù)據(jù)計算相對STEC的關系式

STECΔ=

(1)

其中，f1、f2為載波頻率；L1、L2為載波相位觀測值；λ1、λ2為載波的波長；N1、N2為相位觀測值的整周模糊度；ε為觀測噪聲。由于本文采用高精度GPS數(shù)據(jù)，求取STEC的相對精度很高，為了減少電離層水平梯度的影響，衛(wèi)星的截止高度角選取為15°。

采用歐洲航天局(European Space Agency，ESA)發(fā)布的SWARM衛(wèi)星Level2級TEC數(shù)據(jù)(https://swarm-diss.eo.esa.int/)，VTEC數(shù)據(jù)采樣間隔為1Hz，將相同經緯度格網點的VTEC數(shù)據(jù)求均值，作為相應高度的電離層總電子含量。

1.2 太陽-地磁活動

太陽和地磁環(huán)境的變化會造成電離層活動的異常干擾[14-16]，在分析阿拉斯加地震同震電離層異常效應前，應首先對太陽-地磁活動狀況進行判斷。Dst指數(shù)可以反映地磁活動變化情況，正常情況下在-30nT～30nT范圍內，而太陽輻射指數(shù)F10.7可以反映太陽活動的變化情況，通常太陽活動平靜時F10.7指數(shù)小于90SFU。本文利用日本京都地磁數(shù)據(jù)中心提供的磁暴時變化指數(shù)Dst(http://wdc.kugi.kyoto-u.ac.jp/dstae/index.html)數(shù)據(jù)和空間環(huán)境預報中心提供的太陽輻射指數(shù)F10.7(http://www.sepc.ac.cn/F107Index_chn.php)的觀測數(shù)據(jù)對太陽-地磁活動狀況進行分析。

圖2給出了包括震發(fā)當天在內10天的地磁Dst指數(shù)和太陽射電量F10.7的變化。由圖2可知，在包括震發(fā)在內的10天時間里，|Dst|<30，地磁場比較平靜；同時太陽射電量F10.7指數(shù)<90SFU，太陽活動在這個時期很穩(wěn)定，說明震發(fā)當天太陽-地磁活動對電離層造成的影響很小。

(a) Dst指數(shù)

1.3 方法

從STEC中提取地震電離層擾動信號，需要去除由GPS衛(wèi)星的沿軌運動、電離層的空間變化和時間變化等引起的電離層TEC背景的趨勢項。因此，本文采用SSA的方法進行提取處理，該方法通過對一維時間序列的分解與重構，識別出原始序列中的不同信號，提取出趨勢分量、周期分量與噪聲分量等。

(1)構建遲滯矩陣

對反演出的震發(fā)時段UTC 9-10h STEC數(shù)據(jù)，按照15s的采樣頻率，將1h的STEC數(shù)據(jù)長度N定義為240，將其組成x1,x2,…,xN的時間序列{x}，建立時滯矩陣。

(2)奇異值分解(Singular Value Decomposition，SVD)

定義矩陣S=XXT，XT為X的轉置矩陣，設λ1,…,λM為矩陣S的特征值，U1,…,UM分別為λ1,…,λM對應的特征向量，其中λ1≥…≥λM≥0。

(3)分組

將初等矩陣{1,…,d}劃分為m個不相交的子集I1,I2,…,Im，其中I={i1,…,ip}。軌跡矩陣X的SVD可以表達為：XI=XI1+…+XIm。分組即指對I1,I2,…,Im的確定過程。

(4)對角平均化

將分組得到的矩陣轉換為一系列新的長度為N的時間序列。長度為N的時間序列定義為重建成分(Reconstruction Components，RC)，原始序列即所有重建成分之和。圖3給出了前15重構成分的w-correlation圖。由圖3可知，分解的前5項重構向量與其他重構向量獨立性較好。

圖3 STEC數(shù)據(jù)SSA前15重構成分的w-correlation

(5)獲取擾動信號

利用SVD得到l個重構分量和對應的特征向量λi(1 ≤i≤l)，將重構分量的特征值λK按照大小排列λ1≥λ2≥ …≥λl≥ 0。特征值的大小對應的特征向量代表了TEC信號的變化趨勢的大小?？紤]到前5個重構分量的特征值累積貢獻率超過95%，截取前5個較大的特征值，找到對應的xk之和，即可以重建出充分反映TEC原序列的整體特征，即

(2)

定義包含電離層短周期信號的重構STEC時間序列主成分為STECmain，為提取擾動信號提供高精度的背景參考值。

用原始STEC時間序列減去STECmain，得到電離層STEC擾動信號

ΔSTEC=STEC-STECmain

(3)

2 結果與分析

2.1 CID現(xiàn)象

選取震中附近測站和GPS PRN05號衛(wèi)星的STEC數(shù)據(jù)，基于SSA方法提取震發(fā)時刻UTC 9:00-10:00 ΔSTEC時間序列。圖4給出了各測站的ΔSTEC時間序列和衛(wèi)星軌跡分布圖。由圖4(a)可知，在震后0.5h內，測站ab07、ac42、av06和av10的ΔSTEC序列出現(xiàn)了明顯的N形波擾動，且先后出現(xiàn)2個峰值點。在測站ac42的ΔSTEC序列中，UTC 9:39開始出現(xiàn)擾動，經過了2次由負到正的變化，第1個峰值點出現(xiàn)在UTC 9:43，擾動幅度達到了0.15個TECU;經過幾分鐘后，在UTC 9:49出現(xiàn)了第2個峰值點，擾動幅度為0.06個TECU，在UTC 9:54后恢復平靜，整個擾動過程持續(xù)了15min左右。第1個峰值點的SIP(Sub-Ionospheric Point)距離震中177.98km，第2個峰值點的SIP距離震中186.7km，且SIP位于震中西南方向區(qū)域,如圖4(b)所示。

(a) ΔSTEC時間序列

2.2 觀測衛(wèi)星選取

由于地震產生的聲波會向各個方向傳播，能量在傳播過程中不斷衰減，因此在基于STEC擾動分析同震電離層異?，F(xiàn)象過程中，不僅要保證測站對于觀測衛(wèi)星有連續(xù)的觀測信號，還要選取合適的GPS觀測衛(wèi)星，才能確保阿拉斯加地震同震電離層異常探測的有效性。

本文選取測站av40的STEC擾動結果進行分析，如圖5(a)所示，在UTC 9:00-10:00經過測站上空的6顆具有連續(xù)觀測信號的衛(wèi)星中，只有PRN05號衛(wèi)星在UTC 9:40探測出了明顯的擾動現(xiàn)象，擾動幅度達到0.12個TECU，峰值點出現(xiàn)在UTC 9:44和UTC 9:49，擾動過程持續(xù)了15min，其他觀測衛(wèi)星未發(fā)現(xiàn)明顯的擾動現(xiàn)象。圖5(b)給出了6顆衛(wèi)星的SIP軌跡，可以看出只有PRN05號衛(wèi)星SIP軌跡在地震附近區(qū)域。

(a) ΔSTEC時間序列

由于聲波高度角低時，其在對流層和中氣層中接近水平方向傳播；高度角高時，聲波沿著垂直方向到達電離層，再沿水平方向在電離層中傳播[5]，只有高度角較大的聲波才能傳播到電離層高度引起電離層擾動。聲波只有通過衛(wèi)星和測站的連線路徑，才能探測到震中上空附近電離層擾動情況，而高緯度地區(qū)測站觀測到的衛(wèi)星數(shù)目較少且SIP大多距離測站超過2000km，因此選擇合適的GPS衛(wèi)星對該地區(qū)同震電離層研究特別重要。

2.3 CID傳播速度和擾動幅度

為了進一步分析CID傳播特征，利用高度為350km的電離層單層模型獲得CID的位置信息。在分析STEC擾動過程中，由衛(wèi)星接收機連線與電離層薄層的交點可以獲得電離層穿刺點(Ionosph-eric Pierce Point, IPP)位置。IPP在地面上的投影點為層下點SIP，計算CID峰值點處的SIP經緯度，根據(jù)CID出現(xiàn)的時間和對應的SIP點震中距即可計算CID的水平傳播速度。

圖6給出了2018.1.23 UTC 9:00-10:00 震中附近ΔSTEC擾動序列的時間-距離圖。根據(jù)CID峰值點出現(xiàn)的時間和相應的位置，利用最小二乘進行線性擬合，斜率近似等于震后UTC 9:00-10:00出現(xiàn)的CID傳播速度。

圖6 阿拉斯加地震ΔSTEC變化走勢(UTC 9：00-10：00),藍色虛線表示震發(fā)時刻，黑線代表根據(jù)CID峰值點擬合的趨勢線

由圖6可以看出，CID按照特定的速度進行傳播，并且發(fā)現(xiàn)CID振幅隨著距離增加而減小。CID在震后8min距離震中42.91km處最先被探測到，出現(xiàn)在震后10～20min內，和2011年Mw9.1級日本地震、2008年Mw7.8級汶川地震的CID出現(xiàn)時間一致[17]。CID的傳播速度為(2.62±0.04)km/s，符合地震產生的瑞利波向上傳播到電離層高度引起的CID水平傳播速度[5]。

2.4 CID方向性

利用GPS跟蹤站的觀測信號，可以得到GPS跟蹤站與PRN05、23、31衛(wèi)星的SIP軌跡，同時得到跟蹤站與衛(wèi)星連線的STEC擾動信息。本文基于阿拉斯加地震震中附近CORS網測站，并分析這些站點上空的電離層擾動情況，結果如圖7所示。

圖7 UTC 09:00-10:00阿拉斯加地區(qū)SIP點分布圖，震中上方和西南方向擾動情況及SIP分布，顏色深淺代表電離層擾動幅度

圖7給出了震發(fā)當天UTC 9:00-10:00在震中上方和西南方向探測到的電離層異常和SIP軌跡，可以看到明顯的CID，擾動幅度達到0.16個TECU，出現(xiàn)在距離震中42.9～488km范圍內。圖4(b)給出了震中東北方向的電離層擾動和SIP軌跡，并沒有探測到明顯的CID信號，CID傳播表現(xiàn)出了明顯的方向差異性。

由圖7可以看出，CID從震中沿著西南方向傳播，一方面由于瑞利波的傳播受到地震破裂帶的影響，引起CID傳播方向和破裂方向保持一致；另一方面，由于CID現(xiàn)象的存在主要是通過帶電粒子的運動形式反映的，地磁場施加的庫侖力使其更容易沿著磁力線運動，因此在北半球的CID表現(xiàn)出向南傳播的趨勢。

為了進一步分析瑞利波傳播至電離層高度引起的擾動情況，圖8給出了震中上方和西南方向CID峰值點分布位置，可以看出地震上空的擾動值最大，達到0.16個TECU，CID幅度沿著西南方向隨著距離增加而減小，然后逐漸消失，這進一步證明了瑞利波引起的CID沿著西南方向傳播，并且CID能量以較快速度衰減，符合瑞利波特性[17]。

圖8 震中上方和西南方向CID峰值點示意圖，震源球代表震發(fā)位置，顏色深淺代表峰值點幅度

2.5 SWARM衛(wèi)星數(shù)據(jù)

利用軌道高度為460km的SWARM-A衛(wèi)星觀測得到的VTEC數(shù)據(jù)，選取了包括震發(fā)當天在內前后3天的觀測數(shù)據(jù),對震發(fā)上空電離層異常進行研究。SWARM衛(wèi)星在固定軌道高度運行，在經過同一區(qū)域時相對于地面高度一致，且穿過電離層高度。

圖9(a)給出了衛(wèi)星經過的SIP軌跡的位置及震中位置的分布圖，可以看到,SWARM-A衛(wèi)星經過震中和西南方向區(qū)域上空，持續(xù)時間為3min左右。震發(fā)當天從UTC 20:08(黑色軌跡)開始經過震中區(qū)域，與相同區(qū)域前一天UTC 20:43(紅色軌跡)和后一天UTC 19:17(藍色軌跡)總電子含量進行對比。圖9(b)給出了前后3天在相同高度444～446km上的VTEC，可以看出震發(fā)當天與前一天和后一天相比，總電子含量明顯變小。在排除了太陽和地磁活動影響的情況下，電子含量明顯降低可能是由于地震造成的。

(a) 衛(wèi)星軌跡及震中位置

3 結論

本文基于CORS提供的GPS觀測數(shù)據(jù)，解算震中上空高精度的STEC數(shù)據(jù)，利用SSA方法提取電離層擾動信息，分析阿拉斯加Mw7.9級地震CID現(xiàn)象。

分析表明，震后8min出現(xiàn)了明顯的電離層擾動，先后出現(xiàn)2個峰值點，擾動幅度達到0.16個TECU，擾動時間持續(xù)了15min。CID在距震中600km的范圍內傳播，傳播速度為2.62km/s,與瑞利波水平方向傳播速度一致，利用頻譜分析發(fā)現(xiàn)異常時刻中心頻率為4.87mHz,符合瑞利波頻率。在阿拉斯加地震的CID現(xiàn)象中，發(fā)現(xiàn)擾動具有明顯的方向性，CID傳播方向受到地震破裂帶和地磁場的影響。CID在西南方向傳播的過程中，能量以較快速度衰減，符合瑞利波特性。結合聲波傳播高度角和衛(wèi)星SIP軌跡，發(fā)現(xiàn)緯度較高地區(qū)的CID探測需要先選取合適的觀測衛(wèi)星。最后,利用SWARM衛(wèi)星的VTEC數(shù)據(jù)，發(fā)現(xiàn)震發(fā)當天震中上空區(qū)域總電子含量明顯減小，進一步證明了阿拉斯加地震存在電離層擾動現(xiàn)象。

上述研究表明，利用SSA方法處理由GPS數(shù)據(jù)反演得到高精度的STEC數(shù)據(jù)，可以探測電離層異常效應。此外，需要充足的高時空分辨率的電離層數(shù)據(jù)進行更詳細的探測，并結合大量震例研究地震釋放能量方式對電離層的影響，進一步了解地震電離層耦合機制。