0 引言

关于地震—电离层耦合关系的研究，一般认为始于1965年美国阿拉斯加大地震，此次地震首次证实电离层扰动与地震的发生存在某种关联（Leonard et al，1965）。此后，国内外学者基于地震—电离层耦合效应的前兆异常信息判断进行了大量研究，结果表明，由地震引起的电离层异常确实存在，一般在发震前一周内出现，具有明显的短临特性（Liu et al，2000；Pulinets et al，2009；Le et al，2011；刘静等，2014；张学民等，2014）。在众多电离层监测手段中，基于GPS观测的垂直电子总含量（简称VTEC，下文统称VTEC）被广泛应用于临震电离层扰动现象的监测和研究。例如：姚璐等（2014）利用全球电离层VTEC数据和COSMIC掩星探测电子密度剖面，分析2010年4月13日青海玉树M_S 7.0地震前电离层异常变化，研究表明4月13日震中南部存在VTEC异常增强现象，且具有磁共轭效应；王新志等（2014）利用CODE全球电离层VTEC数据和国际参考电离层模型（IRI）理论数据，对比分析2013年四川芦山M_S 7.0地震前后电离层异常变化，表明芦山地震前后震中电离层VTEC出现异常扰动，全球异常扰动存在自西向东移动趋势及赤道对称现象；王泽民等（2016）利用震中区域IGS站点数据和CODE全球电离层VTEC数据，分析2015年4月25日尼泊尔M_S 8.1地震前后电离层VTEC异常变化，研究表明，排除太阳、地磁活动及地球物理现象短时扰动影响，尼泊尔地震前电离层VTEC值在4月23日表现出显著的正异常。可见，利用全球电离层VTEC数据能够有效监测大地震发生前电离层异常扰动变化。然而，全球电离层VTEC数据结果延迟约3天，且全球GPS站点空间分辨率较低，导致解算结果不准确，对于地震短临预测、预报效力存在明显不足。因此，利用GPS数据，从不同空间尺度分析地震发生前后电离层异常扰动变化，快速、有效地提取与地震相关的电离层地球物理异常信息显得更为重要。

本文针对2019年6月17日15:09（UT）（北京时间23:09）四川长宁M_S 6.0地震发生前后电离层的演化特征，选取CODE全球电离层格网数据、中国地壳运动观测网络GPS数据、震中区域站点GPS数据，分别计算全球、中国区域VTEC异常分布及震中区域站点VTEC时序变化，尝试基于不同空间尺度，系统分析地震前后电离层异常扰动特征变化，为进一步加强川滇区域临震电离层异常扰动监测及地震分析预报提供震例依据。

1 数据

选取CODE全球电离层格网数据、中国地壳运动观测网络GPS数据、震中区域站点GPS数据，尝试以3种空间尺度分析2019年长宁地震前后电离层异常扰动变化。

（1）采用CODE全球电离层格网数据，选取2019年5月17日至6月20日四川长宁地震前后共35天全球VTEC格网数据（经纬度分辨率5°×2.5°），分析地震前后全球电离层异常特征变化。

（2）采用中国地壳运动观测网络60个均匀分布基准站相同时间段GPS数据，选择球谐函数作为电离层VTEC的拟合模型，建立大规模高精度格网电离层模型（林剑等，2009；蔡华等，2014），获取中国区域VTEC格网数据（经纬度分辨率1°×1°），分析地震前后震中及中国区域上空电离层异常特征变化。

（3）采用震中区域相同时间段GPS站点（宜宾、马边、冕宁、宁南）数据，通过单站点VTEC解算（熊波，2012），获取地震前后站点上空VTEC值，分析该区域站点地震前后VTEC时间序列变化。震中区域GPS站点分布见图 1，图中：五角星为长宁震中，三角形为震中附近4个GPS站点。

2 异常提取方法 2.1 中国区域电离层格网数据解算

对于中国区域电离层格网数据解算，可采用以下方法。电离层VTEC模型如下：利用GPS双频接收机2个频率载波相位观测值及其加载的伪距码观测值，通过载波相位平滑伪距观测量，形成电离层残差组合（Dach et al，2015），公式如下

$ P_{4}=40.28\left(\frac{1}{f_{1}^{2}}-\frac{1}{f_{2}^{2}}\right) \mathrm{TEC}+c \cdot \Delta b^{\mathrm{s}}+c \cdot \Delta b_{\mathrm{r}} $

(1)

式中：P₄为载波相位平滑伪距，c为光速，f₁、f₂为卫星载波频率，Δb^s是卫星系统硬件延迟，Δb_r是接收机硬件延迟，TEC为倾斜路径总电子含量。

采用球谐函数对电离层VTEC进行建模，数学表达式（Dach et al，2015）为

$ \operatorname{VTEC}(\beta, s)=\sum\limits_{n=0}^{n_{\max }} \sum\limits_{m=0}^{n} \widetilde{P}_{n m}(\sin \beta)\left(a_{n m} \cos (m \cdot s)+b_{n m} \sin (m \cdot s)\right) $

(2)

式中：VTEC(β，s)为垂直总电子含量，n_max为球谐展开的最大阶数，β为穿刺点地理纬度，s为日固系中穿刺点经度，$\widetilde{P}_{n m}$为n阶m次归一化缔合勒让德函数，a_nm和b_nm为待估计电离层模型参数。TEC与VTEC数学表达式（Dach et al，2015）为

$ \mathrm{TEC}=\mathrm{VTEC} \cdot F(Z) $

(3)

$ F(Z)=\frac{1}{\cos Z^{\prime}} \quad \sin Z^{\prime}=\frac{R}{R+H} \sin (\alpha Z) $

(4)

式中：F(Z)为电离层投影函数，单层模型高度H取350 km，α为0.978，R为6 387 km，Z为卫星天顶距。

将式（2）、式（3）、式（4）代入式（1），形成最终解算方程（谢益炳等，2014），则

$ \begin{aligned} P_{4}=& 40.28\left(\frac{1}{f_{1}^{2}}-\frac{1}{f_{2}^{2}}\right) F(Z) \sum\limits_{n=0}^{n_{\max }} \sum\limits_{m=0}^{n} \widetilde{P}_{n m}(\sin \beta) \cdot \\ &\left(a_{n m} \cos (m \cdot s)+b_{n m} \sin (m \cdot s)\right)+c \cdot \Delta b^\text{s}+c \cdot \Delta b_{\mathrm{r}} \end{aligned} $

(5)

式中：a_nm、b_nm、Δb^s、Δb_r为待估计未知量，采用每2 h多项式模型与单日固定值的接收机硬件延迟和卫星硬件延迟共同建立方程，利用最小二乘求解，解算每2 h模型参数和单日硬件延迟，即可解算中国区域电离层格网数据（王泽民等，2016）。选择球谐函数作为电离层VTEC拟合模型，算法成熟，解算结果可信度较高。

2.2 单站点VTEC解算

中国区域电离层格网数据能够反映区域电离层异常变化趋势，但时间分辨率较低，为了获取高分辨率VTEC数据结果进而分析GPS站点上空VTEC时序变化，采用电离层GPS-TEC反演算法获取斜向TEC，通过映射函数获取站点上空真实VTEC数据，时间分辨率为15 min。

相对TEC的测量精度为10¹³数量级，绝对TEC的测量精度为10¹⁶数量级，利用伪距观测的高可靠性和载波相位观测的高精度，采取载波相位观测量平滑伪距观测量等相关措施，可获得传播路径上高精度斜向TEC。对于给定的接收机和卫星，设每段载波相位连续观测的数据弧段上的测量历元数为N，对于数据弧段上的第i个历元，由差分伪距获得的绝对TEC数据为TECSa_i，同一时刻由差分载波相位获得的相对TEC数据为TECSr_i，偏差B利用TECSa_i、TECSr_i数据所在弧段上所有数据点之差的平均值进行估算（熊波，2012），即

$ B=\frac{1}{N} \sum\limits_{i=1}^{N}\left(\mathrm{TECS} a_{i}-\mathrm{TECS} r_{i}\right) $

(6)

进一步得到传播路径上的斜向TECS_i，表达式如下

$ \mathrm{TECS}_{i}=B+\mathrm{TECS} r_{i} $

(7)

斜向TEC采用电离层薄层模型，通过电离层投影函数最终获取单站点高精度VTEC，电离层GPS-TEC反演算法流程见图 2。

2.3 异常检验方法

选用滑动四分位距法提取VTEC异常信息，选取前15天相同时刻的观测值求取相应四分位距（IQR）及中位数（M），建立上下边界阈值M ±1.5×IQR，当同一时刻观测值超出上或下边界时，该观测值视为异常，且高于上边界为正异常，低于下边界为负异常。滑动四分位距法M ±1.5×IQR的判断阈值，约为标准差的2倍，检测置信度为95%。

从二维空间分布上了解长宁地震前后全球及中国区域电离层VTEC异常分布，分析VTEC异常的ΔTEC空间分布，计算方法如下

$ \Delta \mathrm{TEC}=\left\{\begin{array}{ll} K-L_{1} & \left(K>L_{1}\right) \\ K-L_{2} & \left(K<L_{2}\right) \\ 0 & \left(L_{2} \leqslant K \leqslant L_{1}\right) \end{array}\right. $

(8)

其中，ΔTEC为VTEC的扰动异常，L₁和L₂分别为VTEC背景上下边界阈值，K为待分析时刻VTEC值。

3 VTEC异常分布及时间序列变化 3.1 空间活动性

电离层参数变化与太阳活动、磁暴及空间天气变化紧密相关，有必要将空间环境变化与孕震环境变化对电离层扰动的影响区分出来。选取2019年6月1日至6月20日空间活动性指数数据，绘制F_10.7、K_p、D_st指数变化曲线，结果见图 3（其中F_10.7 =120 sfu、K_p= 4、D_st= -30 nT分别表示异常阈值），分析此次长宁地震前后空间天气变化对电离层异常的影响。

由图 3可见，反映太阳活动影响的指数F_10.7在整个时段内变化幅度较小，数值在60—80 sfu稳定变化，表明该时段太阳辐射影响相对较小；地磁活动指数K_p和D_st的变化表明，2019年6月全球地磁活动总体平稳，其中6月9日K_p指数超过5，地磁扰动较大，其余时段K_p和D_st指数变化较小，地磁活动较为平静，未发生明显磁暴现象。

3.2 全球VTEC异常分布

根据CODE全球电离层格网数据，采用滑动四分位距法提取异常信息，获取VTEC异常的ΔTEC空间分布。图 4给出2019年6月14日08:00— 14:00 UT（世界时）全球ΔTEC分布，图中红色五角星为长宁震中位置。由图 4可见，6月14日12:00 UT在震中区域附近存在明显的电离层异常扰动现象，持续时间较短，电离层异常区域的局部性、集中性和显著性特征较为明显。参考6月14日空间活动性指数数据，可以排除太阳活动、磁暴对电离层异常扰动的影响。

3.3 中国区域VTEC异常分布

根据中国地壳运动观测网络GPS数据，选择球谐函数作为电离层VTEC的拟合模型，获取中国区域VTEC格网数据，并采用滑动四分位距法提取异常信息，获取VTEC异常的ΔTEC空间分布。图 5给出2019年6月14日08:00 — 14:00（UT）中国区域ΔTEC分布，其中红色五角星为长宁震中位置。由图 5可见，6月14日震中区域附近存在明显的电离层正异常，08:00（UT）（北京时间16:00）异常出现，10:00—12:00（UT）（北京时间18:00—20:00）震中区域附近电离层异常达到最大，震中位于异常区域北侧，且电离层异常分布存在沿磁力线向磁道偏移的趋势，14:00（UT）（北京时间22:00）异常向西移动逐渐消失。该分析结果与全球ΔTEC分布基本吻合，表明长宁地震前3天，震中区域附近出现连续电离层正异常扰动，其短临特征明显。

3.4 震中临近站点上空VTEC时间序列变化

为了从时间尺度验证震中区域电离层异常变化与地震发生的相关性，根据震中区域GPS站点数据，通过单站点VTEC数据解算，获得2019年6月1日至6月20日震中区域附近4个GPS站点（宜宾、马边、冕宁、宁南）电离层VTEC时间序列曲线，见图 6，其中红色曲线为VTEC观测值，蓝色曲线为VTEC判断阈值上下限，箭头标志为电离层异常扰动。由图 6可见，宜宾（YBIN）、马边（SCMB）、冕宁（SCMN）、宁南（SCNN）4个站点6月1日、6月4日、6月14日电离层异常扰动明显，查看F_10.7、K_p、D_st指数，发现均未超过阈值。

为进一步分析震中区域GPS站点上空电离层异常扰动原因，详细查看全球及中国区域VTEC异常分布，文中给出2019年6月1日和4日的04:00 —10:00 UT中国区域ΔTEC分布，结果见图 7、图 8，图 8中红色五角星为2019年6月4日中国台湾台东海域地震震中。

由图 7可见，6月1日中国区域电离层异常扰动不明显，未表现出地震—电离层耦合的局部性、集中性和显著性特征，可能与空间天气及其他因素有关。由图 8可见，6月4日中国区域电离层正异常扰动明显，且全球ΔTEC分布并未出现显著的多区域电离层异常现象。因此，该异常与地震相关的可能性较大。根据以往震例统计，电离层异常扰动显著出现在震前一周，而该异常出现在长宁地震前13天，与长宁地震的相关性较小。根据中国地震台网测定，北京时间2019年6月4日17时46分，中国台湾台东海域发生M_S 5.8地震，电离层异常出现在震前4—6小时，是否与此次地震具有相关性需进一步分析。

然而，6月14日全球、中国区域VTEC异常分布均显示，在震中区域附近出现明显的电离层正异常扰动，同时在长宁地震发生前后，该异常区域并未发生其他强震事件。可见，6月14日震中区域附近电离层异常扰动与长宁M_S 6.0地震有一定相关性。

4 结论

通过对2019长宁M_S 6.0地震前后全球、中国区域VTEC异常分布及震中区域站点VTEC时间序列变化分析，可得出如下结论。

2019长宁M_S 6.0地震发生前3天（6月14日），全球、中国区域VTEC异常分布均显示，在震中区域上空出现明显的电离层异常扰动现象，异常空间分布呈现向磁赤道偏移趋势，该异常现象与地震—电离层耦合特征一致。单站点VTEC结果从时间尺度验证了震中区域电离层异常变化与长宁M_S 6.0地震有一定相关性，可能是地震孕育过程中地震破裂带化学物质激发累积，产生大气电场异常，进而引起电离层异常扰动，进一步佐证了地震—电离层效应的短临特征和局地性特点。

利用GPS观测站数据，快速获取全球、中国区域VTEC异常分布及川滇区域GPS站点上空VTEC时间序列变化，可实时有效地监测川滇区域电离层异常扰动，尝试将电离层地球物理异常信息用于川滇区域临震监测预报。

感谢CODE提供全球电离层格网数据。