與地震相關(guān)的電離層三頻信標(biāo)CT診斷——實驗及初步結(jié)果*

廖 力 趙庶凡

1) 中國地震局地球物理研究所， 北京100081 2) 中國地震局地震預(yù)測研究所， 北京100036

與地震相關(guān)的電離層三頻信標(biāo)CT診斷
——實驗及初步結(jié)果*

廖 力1)※趙庶凡2)

1) 中國地震局地球物理研究所， 北京100081 2) 中國地震局地震預(yù)測研究所， 北京100036

已有研究表明地震發(fā)生前會出現(xiàn)電離層電子密度異常。 三頻信標(biāo)探測通過對電離層進(jìn)行快速掃描， 反演得到電子密度的空間分布， 能夠更清晰地反映地震電離層異常特征， 目前用于地震觀測的三頻信標(biāo)觀測剛剛開展。 本研究介紹一種基于三頻信標(biāo)的電離層反演算法， 采用120°E觀測鏈的數(shù)據(jù)進(jìn)行電離層層析成像， 并與電離層垂測觀測結(jié)果進(jìn)行對比， 研究觀測鏈路上低緯和高緯電離層的時變特征。 結(jié)果表明， 本研究采用的電離層層析方法獲得的F2層電子密度相對于初始輸入的IRI模型更接近垂測觀測， 達(dá)到了一定的反演精度， 可反映電離層的空間及時間變化特征。

三頻信標(biāo); 電離層層析成像; 地震

引言

近年研究發(fā)現(xiàn)， 地震前和地震期間會出現(xiàn)電離層參數(shù)異常、 電離層電子總含量TEC發(fā)生變化、 電離層F層擾動等現(xiàn)象。 通過衛(wèi)星監(jiān)測地震危險區(qū)的電離層變化， 可獲得與地震有關(guān)的電離層變化特征， 可能發(fā)展成為地震預(yù)報的一種有效手段。

已有研究表明， 地震發(fā)生前會引起電離層TEC(電離層總電子含量)異常， 并能被GNSS接收機(jī)、 電離層雙頻接收機(jī)以及三頻信標(biāo)接收機(jī)等多種觀測手段探測到。 Calais和Minster[1]最早采用GPS探測北嶺(North bridge)地震后電離層的TEC情況。 Liu等[2-3]分析了1999—2001年臺灣M≥6.0以上地震期間電離層TEC變化， 發(fā)現(xiàn)地震前1～5天TEC有顯著的減少。 祝芙英等[4]利用GPS數(shù)據(jù)對汶川地震前的TEC異常進(jìn)行了分析， 也發(fā)現(xiàn)震前1～6天(除5月9日電離層TEC出現(xiàn)異常增加外)時間內(nèi)多次出現(xiàn)顯著的TEC異常減少。 陳必焰等[5]利用GPS數(shù)據(jù)對2011年日本9.0級進(jìn)行了電離層反演， 認(rèn)為2月28日UT14：00～16：00出現(xiàn)的電離層電子密度值異常減小、 3月2日UT08：00～14：00、 3日UT00：00～06：00和4日UT12：00～20：00出現(xiàn)的電離層電子密度值異常增大極有可能與日本地震有關(guān)。 趙必強(qiáng)等[6]認(rèn)為， TEC減小基本可歸結(jié)為地震引起的電離層效應(yīng)。 因此， 利用觀測資料進(jìn)行電離層反演分析TEC變化， 對于識別地震前兆異常具有重要意義。

相對于其他觀測手段， 三頻信標(biāo)在VHF、 UHF和L頻段上輸出頻率穩(wěn)定相位相關(guān)的載頻信號， 并經(jīng)全向天線向預(yù)定覆蓋區(qū)域輻射， 通過地面接收設(shè)備在單站或多站進(jìn)行信號幅度和相位測量， 電離層層析掃描技術(shù)可以更高精度地獲得電離層電子總含量、 電離層閃爍和電子密度及電離層不規(guī)則體信息， 因此， 反演得到電子密度分布更為精確， 可以探測到更小尺度的電離層擾動[7]。 本文介紹了一種三頻信標(biāo)電離層反演算法及其實現(xiàn)方法， 利用120°E鏈路上三頻信標(biāo)數(shù)據(jù)進(jìn)行了電離層層析成像， 并與電離層垂測實驗的結(jié)果進(jìn)行了對比， 同時對算法進(jìn)行了電離層電子密度時變特征的應(yīng)用研究。

1 三頻信標(biāo)反演算法

三頻信標(biāo)反演電子密度的計算流程是通過三頻信標(biāo)測量獲得從發(fā)射站到接收站射線路徑上的總電子密度含量， 然后通過迭代算法反演電子密度。 三頻信標(biāo)機(jī)的測量值僅有多普勒相位值， 只能獲得相對電子含量， 要得到絕對電子含量需要求解相位積分常數(shù)。 本文將從雙頻信標(biāo)的雙站法出發(fā)介紹相位積分常數(shù)的求解方法。

1.1 相位積分常數(shù)的求法

衛(wèi)星(S)至接收機(jī)(R)路徑的幾何關(guān)系圖如圖1所示。S1表示從接收站R1到衛(wèi)星的射線路徑，S2表示從接收站R2到衛(wèi)星的射線路徑，hS表示從地面到衛(wèi)星的垂直高度。

圖1 將斜向全電子含量轉(zhuǎn)換成垂直全電子含量示意圖

對于接收站R1， 斜向全電子含量IS1定義為沿傳播路徑上電子密度的線積分， 與投影的垂直全電子含量IV1的關(guān)系式為：

(1)

其中，D1稱為路徑S1上斜向全電子含量與垂直全電子含量的轉(zhuǎn)換函數(shù)，χ1為路徑S1與接收站R1垂線方向的夾角。 衛(wèi)星過境期間斜向全電子含量IS1與差分多普勒相位關(guān)系式表示為：

(2)

(3)

其中，D1k和IV1k分別為第k個測點與R1對應(yīng)射線的轉(zhuǎn)換函數(shù)和垂直TEC。 同理， 對于接收點R2有：

(4)

其中，φd2(tk)是接收機(jī)R2觀測給出的雙頻差分多普勒相位；φ2(0)為待求的R2接收站的未知積分常數(shù)；D2k和IV2k為第k個測點與R2對應(yīng)射線的轉(zhuǎn)換函數(shù)和垂直TEC。

式(3)和式(4)可以變?yōu)椋?/p>

(5)

近似地認(rèn)為， 對于地面不同位置的兩個接收站， 兩個接收站的垂向TEC相等。 因此， 對于雙站法IV1k≈IV2k， 利用最小平方誤差限定IV1k和IV2k兩者間的誤差， 兩站法的目標(biāo)函數(shù)可表示為：

(6)

從而代入式(2)可由觀測得到的差分相對TEC， 得到各個站的雙頻絕對TEC， 如下式所示：

(7)

推廣到m個接收站， 也可按相同的概念求解多個站的初始相位， 即所謂的多站法(multi-station method)。 多站法的目標(biāo)函數(shù)為公式(8)。

(8)

通過求解E的極小值也可得到各個站的位置相位積分常數(shù)φi(0)(i=1, …,m)， 從而求得各個站的雙頻絕對TEC。

(9)

1.2 代數(shù)重建算法

我們采用代數(shù)重建算法(ART)來進(jìn)行重建影像的微調(diào)。 代數(shù)重建算法需要給定一個模擬的初始值， 進(jìn)行多次迭代獲得最終的結(jié)果。 反演電離層電子密度的ART算法為：

(10)

其中， ai表示系數(shù)矩陣A的第i行行向量， di表示觀測值向量D的第i個元素， X(k)表示對未知數(shù)的初值X(0)的第k輪迭代結(jié)果。 λk為松弛因子， 取值為0～2， 對于含測量誤差的數(shù)據(jù)， 適當(dāng)選取松弛因子能改善重建圖像的質(zhì)量和迭代的效率。 一般情況下， 在迭代過程中松弛因子可以取一個定值。 由于ART并不保證迭代結(jié)果為正， 考慮實際的物理情形， 需要在ART算法中加上電子密度非負(fù)的約束。

2 層析結(jié)果與垂測資料的對比實驗

利用臺灣低緯電離層層析成像網(wǎng)(Low-latitude Ionosphere Tomography Network， LITN)的極軌衛(wèi)星信標(biāo)觀測數(shù)據(jù)進(jìn)行電離層CT反演， 為提高反演精度需盡可能選取接收低軌衛(wèi)星信號重復(fù)較多的數(shù)據(jù)， 因此， 本研究選取了2012年5月9日位于臺灣的4個位置相鄰的三頻信標(biāo)接受臺站數(shù)據(jù)進(jìn)行電離層層析成像， 臺站位置如圖2所示。

本研究以2012年5月9日地方時18點22分的國際參考電離層模型(International Reference Ionosphere， IRI)作為反演的初始數(shù)據(jù)， 利用ART算法進(jìn)行了電離層層析成像。 反演區(qū)域為10°N～30°N， 121.4°E， 水平分辨率為0.5°， 垂直方向為100～800 km， 垂直分辨率為20 km， 反演過程中取至少通過10條射線的網(wǎng)格參與計算。 為避免選取到未完全穿透反演區(qū)域的射線， 去除掉接收站中傾角小于15°的射線。 反演結(jié)果如圖3所示， 反演的收斂過程如圖4所示， 從圖4可以看出， 經(jīng)過約40步迭代計算之后， 觀測值與反演值之差的二范數(shù)不再減小， 計算結(jié)果得到收斂。

圖2 參與反演的三頻信標(biāo)臺站

表1 IRI模型、 反演結(jié)果與垂測數(shù)據(jù)對比

為檢驗反演結(jié)果， 我們將反演結(jié)果與臺灣花蓮地區(qū)(23.6°N， 121.4°E)的垂測數(shù)據(jù)的電子濃度峰值高度hmF2及電子濃度峰值頻率foF2 進(jìn)行比較， 其中foF2由公式(11)可得， 由表1可知， IRI模型與反演結(jié)果的電子濃度峰值與垂測數(shù)據(jù)基本一致， 電子濃度的峰值頻率， 反演結(jié)果相對初始輸入的IRI模型更接近于垂測數(shù)據(jù)的測量值， 但可以看到， 迭代結(jié)果非常依賴于初始值。 上述采用像素基的反演方法非常依賴迭代初值， 因此， 為獲得更好的電離層層析成像結(jié)果， 應(yīng)盡量提高初始值的準(zhǔn)確度， 比如目前已有一些研究采用的先用函數(shù)基模型獲得電子密度初始分布， 再利用像素基模型對初始分布進(jìn)行二次迭代[9]。

圖3 ART算法反演結(jié)果

圖4 ART算法收斂過程橫坐標(biāo)為迭代次數(shù)， 縱坐標(biāo)為迭代誤差

(11)

3 層析算法在120°E鏈路電離層電子密度時變特征中的應(yīng)用研究

中國地震局地震預(yù)測研究所在溫州和上海分別架設(shè)了三頻信標(biāo)接收機(jī)ITS33s， 兩個站是對臺灣低緯電離層層析成像網(wǎng)的觀測鏈路向中高緯地區(qū)的擴(kuò)展。 本研究利用該鏈路的觀測數(shù)據(jù)， 研究鏈路上低緯及中緯的電離層電子密度背景在不同地方時、 不同季節(jié)的變化規(guī)律。 圖5為利用2012年9月17日接收的COSMOS 2414衛(wèi)星信標(biāo)數(shù)據(jù)進(jìn)行反演研究反演結(jié)果。 可見地方時下午16點電子密度最大， 在15°N附近有很大的峰值， 在所有高度都呈現(xiàn)單峰狀態(tài)； 10點電子密度也很大、 強(qiáng)度次之并在F2層高度出現(xiàn)雙峰結(jié)構(gòu)， 雙峰的中心位置約在35°N附近； 太陽落下后的地方時21時和24時電子密度很小， 但在F2層高度仍可見清晰的雙峰結(jié)構(gòu)， 中心位置較10點南移約在25°左右； 夜間3點電子密度最小， 并且雙峰結(jié)構(gòu)完全消失。

比較了2012年9月地方時16點的反演結(jié)果與12月地方時16點的反演結(jié)果， 結(jié)果如圖6所示， 發(fā)現(xiàn)相比于電子密度顯著的日變化， 電子密度的季節(jié)變化不大， 這與IRI模型的計算結(jié)果一致。

圖5 2012年9月17日日變化(LT=10， 16， 21， 24， 3)

圖6 2012年9月17日地方時16點與12月7日地方時16點對比

4 結(jié)論

近幾年法國的DEMETER衛(wèi)星、 俄羅斯的COMPASS-2衛(wèi)星、 美國的QuakeSat微小衛(wèi)星等相關(guān)科研項目中都包括了探測災(zāi)害性事件(地震、 強(qiáng)磁暴等)引發(fā)的大氣層、 電離層和磁層的異?，F(xiàn)象。 日本及中國臺灣學(xué)者也對地震和電離層異常的關(guān)系進(jìn)行過研究， 證實了兩者之間的關(guān)聯(lián)性， 以及通過電離層觀測進(jìn)行地震異常研究的可行性。 三頻信標(biāo)是一種重要的空間環(huán)境探測星載設(shè)備， 可對電離層進(jìn)行快速的層析掃描從而反演得到不同尺度的電離層不均勻結(jié)構(gòu)的參數(shù)特征， 開展基于三頻信標(biāo)的地震電離層異常研究， 可以為地震前兆監(jiān)測提供可能的新手段。 此外， 中國地震局牽頭即將發(fā)射的電磁衛(wèi)星搭載了三頻信標(biāo)發(fā)射機(jī)， 發(fā)展獨立自主的基于三頻信標(biāo)電離層層析算法符合地震電磁衛(wèi)星工程的發(fā)展需要。

本文首先介紹了電離層信標(biāo)絕對TEC的解算方法， 以及基于ART算法的電離層層析成像方法； 接著將利用上述電離層信標(biāo)層析成像方法反演的電離層電子密度結(jié)果與電離層垂測實驗的數(shù)據(jù)進(jìn)行了對比， 并結(jié)合在溫州和上海的三頻信標(biāo)觀測數(shù)據(jù)和臺灣LITN觀測網(wǎng)的數(shù)據(jù)對120°E鏈路的電子密度時變特征進(jìn)行了初步研究。 結(jié)果發(fā)現(xiàn)， 本文采用的反演算法驅(qū)使電離層電子密度成像結(jié)果相對于初始輸入的IRI模型， 更接近電離層垂測實驗的觀測值， 達(dá)到了一定的反演精度， 可反映電離層的空間及時間變化特征。 由于本文采用的反演算法對初始值的輸入較為敏感， 要獲得更高精度的反演結(jié)果， 可考慮采用更接近真實的電子密度觀測值作為反演的初值， 或使用函數(shù)基方法獲得初始分布作為輸入。 另外， 受限于觀測數(shù)據(jù)， 本文只驗證了層析成像結(jié)果中某一臺站位置的F2層最大電子密度， 在后續(xù)研究中需對反演區(qū)域中更多網(wǎng)格的結(jié)果進(jìn)行檢驗。

[1] Calais E， Minster J B. GPS detection of ionospheric perturbations following the January 17， 1994， Northridge earthquake. Geophys. Res. Lett.， 1995， 22(9): 1045-1048

[2] Liu J Y， Chen Y I， Pulinets S A， et al. Seismo-ionosphere signatures prior toM≥6.0 Taiwan earthquake. Geophys. Res. Lett.， 2000， 27(19): 3113-3116

[3] Liu J Y， Chen Y I， Chuo Y J， et al. Variations of ionospheric total electron content during the Chi-Chi earthquake. Geophys. Res. Lett.， 2001， 28(7): 1383-1386

[4] 祝芙英， 吳云， 林劍， 等. 汶川MS8.0地震前電離層TEC異常分析. 大地測量與地球動力學(xué)， 2008， 28(6): 16-21

[5] 陳必焰， 戴吾蛟， 蔡昌盛， 等. 利用電離層層析技術(shù)探測日本9.0級地震前電離層異常. 大地測量與地球動力學(xué)， 2011， 31(6): 11-14

[6] 趙必強(qiáng)， 萬衛(wèi)星， 王敏， 等. 震前電離層擾動研究進(jìn)展及汶川地震前電離層變化. 科技導(dǎo)報， 2008， 26(11)： 30-34

[7] 趙海生. 三頻信標(biāo)電離層擾動探測及層析成像關(guān)鍵技術(shù)研究. 杭州： 杭州電子科技大學(xué)， 2010

[8] 趙海生， 許正文， 吳健， 等. 三頻衛(wèi)星信標(biāo)測量TEC方法探討及數(shù)值模擬. 裝備環(huán)境工程， 2010， 7(4): 9-13

[9] 歐明， 甄衛(wèi)民， 劉裔文， 等．一種基于LEO衛(wèi)星信標(biāo)的電離層層析成像新算法．地球物理學(xué)報， 2015， 35(10)： 3469-3480

Research and application of computer ionospheric tomography algorithm using tri-band beacon data

Liao Li1)， Zhao Shufan2)

1) Institute of Geophysics， China Earthquake Administration， Beijing 100081， China 2) Institute of Earthquake Science， China Earthquake Administration， Beijing 100036， China

The past research has shown that there is ionospheric abnormal before earthquake. Tri-band beacon (TBB) could scan the ionosphere quickly， and the spatial distribution of electron density which reflects the features of seismic ionospheric anomaly more clearly could be got through the inversion. Currently， the TBB observation is applied to seismologic research. In this research， a computer ionospheric tomography (CIT) algorithm based on TBB has been introduced. Using the chain observation data on 120°E， CIT has been done， and results have been compared with the ionosonde data， meanwhile the ionospheric characteristics at low latitude and high latitude have been also studied. Results have shown that the CIT method used in this research could get the electron density in F2 layer with more accuracy than the initial input， and can reflect the spatial and time variation of the ionosphere correctly.

tri-band beacon; computer ionospheric tomography; earthquake

2016-06-06； 采用日期： 2016-10-08。

※通訊作者： 廖力， e-mail: 24817395@qq.com。 基金項目： 國家國際科技合作對俄科技合作專項(2014DFR21280)、 ISSI-BeiJing項目和中國地震局地球物理研究所基本科研業(yè)務(wù)費(fèi)(DQJB16B11)共同資助。

P315.72+2；

A；

10.3969/j.issn.0235-4975.2016.12.006