与地震相关的电离层三频信标CT诊断——实验及初步结果*

廖 力 赵庶凡

1) 中国地震局地球物理研究所， 北京100081 2) 中国地震局地震预测研究所， 北京100036

与地震相关的电离层三频信标CT诊断
——实验及初步结果*

廖 力1)※赵庶凡2)

1) 中国地震局地球物理研究所， 北京100081 2) 中国地震局地震预测研究所， 北京100036

已有研究表明地震发生前会出现电离层电子密度异常。 三频信标探测通过对电离层进行快速扫描， 反演得到电子密度的空间分布， 能够更清晰地反映地震电离层异常特征， 目前用于地震观测的三频信标观测刚刚开展。 本研究介绍一种基于三频信标的电离层反演算法， 采用120°E观测链的数据进行电离层层析成像， 并与电离层垂测观测结果进行对比， 研究观测链路上低纬和高纬电离层的时变特征。 结果表明， 本研究采用的电离层层析方法获得的F2层电子密度相对于初始输入的IRI模型更接近垂测观测， 达到了一定的反演精度， 可反映电离层的空间及时间变化特征。

三频信标; 电离层层析成像; 地震

引言

近年研究发现， 地震前和地震期间会出现电离层参数异常、 电离层电子总含量TEC发生变化、 电离层F层扰动等现象。 通过卫星监测地震危险区的电离层变化， 可获得与地震有关的电离层变化特征， 可能发展成为地震预报的一种有效手段。

已有研究表明， 地震发生前会引起电离层TEC(电离层总电子含量)异常， 并能被GNSS接收机、 电离层双频接收机以及三频信标接收机等多种观测手段探测到。 Calais和Minster[1]最早采用GPS探测北岭(North bridge)地震后电离层的TEC情况。 Liu等[2-3]分析了1999—2001年台湾M≥6.0以上地震期间电离层TEC变化， 发现地震前1～5天TEC有显著的减少。 祝芙英等[4]利用GPS数据对汶川地震前的TEC异常进行了分析， 也发现震前1～6天(除5月9日电离层TEC出现异常增加外)时间内多次出现显著的TEC异常减少。 陈必焰等[5]利用GPS数据对2011年日本9.0级进行了电离层反演， 认为2月28日UT14：00～16：00出现的电离层电子密度值异常减小、 3月2日UT08：00～14：00、 3日UT00：00～06：00和4日UT12：00～20：00出现的电离层电子密度值异常增大极有可能与日本地震有关。 赵必强等[6]认为， TEC减小基本可归结为地震引起的电离层效应。 因此， 利用观测资料进行电离层反演分析TEC变化， 对于识别地震前兆异常具有重要意义。

相对于其他观测手段， 三频信标在VHF、 UHF和L频段上输出频率稳定相位相关的载频信号， 并经全向天线向预定覆盖区域辐射， 通过地面接收设备在单站或多站进行信号幅度和相位测量， 电离层层析扫描技术可以更高精度地获得电离层电子总含量、 电离层闪烁和电子密度及电离层不规则体信息， 因此， 反演得到电子密度分布更为精确， 可以探测到更小尺度的电离层扰动[7]。 本文介绍了一种三频信标电离层反演算法及其实现方法， 利用120°E链路上三频信标数据进行了电离层层析成像， 并与电离层垂测实验的结果进行了对比， 同时对算法进行了电离层电子密度时变特征的应用研究。

1 三频信标反演算法

三频信标反演电子密度的计算流程是通过三频信标测量获得从发射站到接收站射线路径上的总电子密度含量， 然后通过迭代算法反演电子密度。 三频信标机的测量值仅有多普勒相位值， 只能获得相对电子含量， 要得到绝对电子含量需要求解相位积分常数。 本文将从双频信标的双站法出发介绍相位积分常数的求解方法。

1.1 相位积分常数的求法

卫星(S)至接收机(R)路径的几何关系图如图1所示。S1表示从接收站R1到卫星的射线路径，S2表示从接收站R2到卫星的射线路径，hS表示从地面到卫星的垂直高度。

图1 将斜向全电子含量转换成垂直全电子含量示意图

对于接收站R1， 斜向全电子含量IS1定义为沿传播路径上电子密度的线积分， 与投影的垂直全电子含量IV1的关系式为：

(1)

其中，D1称为路径S1上斜向全电子含量与垂直全电子含量的转换函数，χ1为路径S1与接收站R1垂线方向的夹角。 卫星过境期间斜向全电子含量IS1与差分多普勒相位关系式表示为：

(2)

(3)

其中，D1k和IV1k分别为第k个测点与R1对应射线的转换函数和垂直TEC。 同理， 对于接收点R2有：

(4)

其中，φd2(tk)是接收机R2观测给出的双频差分多普勒相位；φ2(0)为待求的R2接收站的未知积分常数；D2k和IV2k为第k个测点与R2对应射线的转换函数和垂直TEC。

式(3)和式(4)可以变为：

(5)

近似地认为， 对于地面不同位置的两个接收站， 两个接收站的垂向TEC相等。 因此， 对于双站法IV1k≈IV2k， 利用最小平方误差限定IV1k和IV2k两者间的误差， 两站法的目标函数可表示为：

(6)

从而代入式(2)可由观测得到的差分相对TEC， 得到各个站的双频绝对TEC， 如下式所示：

(7)

推广到m个接收站， 也可按相同的概念求解多个站的初始相位， 即所谓的多站法(multi-station method)。 多站法的目标函数为公式(8)。

(8)

通过求解E的极小值也可得到各个站的位置相位积分常数φi(0)(i=1, …,m)， 从而求得各个站的双频绝对TEC。

(9)

1.2 代数重建算法

我们采用代数重建算法(ART)来进行重建影像的微调。 代数重建算法需要给定一个模拟的初始值， 进行多次迭代获得最终的结果。 反演电离层电子密度的ART算法为：

(10)

其中， ai表示系数矩阵A的第i行行向量， di表示观测值向量D的第i个元素， X(k)表示对未知数的初值X(0)的第k轮迭代结果。 λk为松弛因子， 取值为0～2， 对于含测量误差的数据， 适当选取松弛因子能改善重建图像的质量和迭代的效率。 一般情况下， 在迭代过程中松弛因子可以取一个定值。 由于ART并不保证迭代结果为正， 考虑实际的物理情形， 需要在ART算法中加上电子密度非负的约束。

2 层析结果与垂测资料的对比实验

利用台湾低纬电离层层析成像网(Low-latitude Ionosphere Tomography Network， LITN)的极轨卫星信标观测数据进行电离层CT反演， 为提高反演精度需尽可能选取接收低轨卫星信号重复较多的数据， 因此， 本研究选取了2012年5月9日位于台湾的4个位置相邻的三频信标接受台站数据进行电离层层析成像， 台站位置如图2所示。

本研究以2012年5月9日地方时18点22分的国际参考电离层模型(International Reference Ionosphere， IRI)作为反演的初始数据， 利用ART算法进行了电离层层析成像。 反演区域为10°N～30°N， 121.4°E， 水平分辨率为0.5°， 垂直方向为100～800 km， 垂直分辨率为20 km， 反演过程中取至少通过10条射线的网格参与计算。 为避免选取到未完全穿透反演区域的射线， 去除掉接收站中倾角小于15°的射线。 反演结果如图3所示， 反演的收敛过程如图4所示， 从图4可以看出， 经过约40步迭代计算之后， 观测值与反演值之差的二范数不再减小， 计算结果得到收敛。

图2 参与反演的三频信标台站

表1 IRI模型、 反演结果与垂测数据对比

为检验反演结果， 我们将反演结果与台湾花莲地区(23.6°N， 121.4°E)的垂测数据的电子浓度峰值高度hmF2及电子浓度峰值频率foF2 进行比较， 其中foF2由公式(11)可得， 由表1可知， IRI模型与反演结果的电子浓度峰值与垂测数据基本一致， 电子浓度的峰值频率， 反演结果相对初始输入的IRI模型更接近于垂测数据的测量值， 但可以看到， 迭代结果非常依赖于初始值。 上述采用像素基的反演方法非常依赖迭代初值， 因此， 为获得更好的电离层层析成像结果， 应尽量提高初始值的准确度， 比如目前已有一些研究采用的先用函数基模型获得电子密度初始分布， 再利用像素基模型对初始分布进行二次迭代[9]。

图3 ART算法反演结果

图4 ART算法收敛过程横坐标为迭代次数， 纵坐标为迭代误差

(11)

3 层析算法在120°E链路电离层电子密度时变特征中的应用研究

中国地震局地震预测研究所在温州和上海分别架设了三频信标接收机ITS33s， 两个站是对台湾低纬电离层层析成像网的观测链路向中高纬地区的扩展。 本研究利用该链路的观测数据， 研究链路上低纬及中纬的电离层电子密度背景在不同地方时、 不同季节的变化规律。 图5为利用2012年9月17日接收的COSMOS 2414卫星信标数据进行反演研究反演结果。 可见地方时下午16点电子密度最大， 在15°N附近有很大的峰值， 在所有高度都呈现单峰状态； 10点电子密度也很大、 强度次之并在F2层高度出现双峰结构， 双峰的中心位置约在35°N附近； 太阳落下后的地方时21时和24时电子密度很小， 但在F2层高度仍可见清晰的双峰结构， 中心位置较10点南移约在25°左右； 夜间3点电子密度最小， 并且双峰结构完全消失。

比较了2012年9月地方时16点的反演结果与12月地方时16点的反演结果， 结果如图6所示， 发现相比于电子密度显著的日变化， 电子密度的季节变化不大， 这与IRI模型的计算结果一致。

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

图6 2012年9月17日地方时16点与12月7日地方时16点对比

4 结论

近几年法国的DEMETER卫星、 俄罗斯的COMPASS-2卫星、 美国的QuakeSat微小卫星等相关科研项目中都包括了探测灾害性事件(地震、 强磁暴等)引发的大气层、 电离层和磁层的异常现象。 日本及中国台湾学者也对地震和电离层异常的关系进行过研究， 证实了两者之间的关联性， 以及通过电离层观测进行地震异常研究的可行性。 三频信标是一种重要的空间环境探测星载设备， 可对电离层进行快速的层析扫描从而反演得到不同尺度的电离层不均匀结构的参数特征， 开展基于三频信标的地震电离层异常研究， 可以为地震前兆监测提供可能的新手段。 此外， 中国地震局牵头即将发射的电磁卫星搭载了三频信标发射机， 发展独立自主的基于三频信标电离层层析算法符合地震电磁卫星工程的发展需要。

本文首先介绍了电离层信标绝对TEC的解算方法， 以及基于ART算法的电离层层析成像方法； 接着将利用上述电离层信标层析成像方法反演的电离层电子密度结果与电离层垂测实验的数据进行了对比， 并结合在温州和上海的三频信标观测数据和台湾LITN观测网的数据对120°E链路的电子密度时变特征进行了初步研究。 结果发现， 本文采用的反演算法驱使电离层电子密度成像结果相对于初始输入的IRI模型， 更接近电离层垂测实验的观测值， 达到了一定的反演精度， 可反映电离层的空间及时间变化特征。 由于本文采用的反演算法对初始值的输入较为敏感， 要获得更高精度的反演结果， 可考虑采用更接近真实的电子密度观测值作为反演的初值， 或使用函数基方法获得初始分布作为输入。 另外， 受限于观测数据， 本文只验证了层析成像结果中某一台站位置的F2层最大电子密度， 在后续研究中需对反演区域中更多网格的结果进行检验。

[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异常分析. 大地测量与地球动力学， 2008， 28(6): 16-21

[5] 陈必焰， 戴吾蛟， 蔡昌盛， 等. 利用电离层层析技术探测日本9.0级地震前电离层异常. 大地测量与地球动力学， 2011， 31(6): 11-14

[6] 赵必强， 万卫星， 王敏， 等. 震前电离层扰动研究进展及汶川地震前电离层变化. 科技导报， 2008， 26(11)： 30-34

[7] 赵海生. 三频信标电离层扰动探测及层析成像关键技术研究. 杭州： 杭州电子科技大学， 2010

[8] 赵海生， 许正文， 吴健， 等. 三频卫星信标测量TEC方法探讨及数值模拟. 装备环境工程， 2010， 7(4): 9-13

[9] 欧明， 甄卫民， 刘裔文， 等．一种基于LEO卫星信标的电离层层析成像新算法．地球物理学报， 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项目和中国地震局地球物理研究所基本科研业务费(DQJB16B11)共同资助。

P315.72+2；

A；

10.3969/j.issn.0235-4975.2016.12.006