2. 中南大学地球科学与信息物理学院，湖南 长沙 410083；
3. 武汉大学测绘学院，湖北 武汉 430079
,
,
一、引 言
近年来，震前电离层TEC(总电子含量)的扰动及地震电离层耦合机制的研究成为热点之一[1]。 越来越多的研究证实了由地震引起的电离层扰动确实存在。 Calais和Minster首次利用GPS探测1994年Northridge地震前后的电离层TEC异常扰动，发现其震中区上空的电离层TEC震后几分钟出现了异常变化[2, 3]。刘正彦等发现靠近震中区域上空的电离层TEC在集集地震期间呈下降趋势[4]。刘正彦在分析了1999—2002 年台湾地区Ms5.0级以上的地震后，发现震前1—5 d TEC有明显的负异常现象，检出异常的震例占总震例的80%[5]。此外，许多学者通过对汶川地震电离层特征的分析研究，均发现了震前电离层的异常扰动[6, 7, 8, 9]。
本文利用日本及周边IGS(International GNSS Service)测站提供的观测数据，计算获得了日本区域内某震例的震前VTEC(vertical total electron content)时间序列，并结合地震前后太阳和地磁场活动，采用滑动平均法和非滑动平均法两种统计学方法对电离层异常现象进行了详细研究，均发现了很可能与地震相关的震前电离层异常现象。此外，本文将两种方法得出的相关结果进行对比分析，验证了本文结果的可靠性，并探讨了两种方法的适用性。
二、电离层VTEC异常检验方法 1. VTEC求解方法基于电离层薄壳模型，即假设电离层中所有的自由电子都集中在一个高度为 hion的与地球同心的无限薄的球面上，通过高度映射函数，求解VTEC值。hion取值接近电子密度最大高度处，文中取350 km。信号在卫星到接收机之间路径上的总电子含量TEC计算公式如下[10]
式中，A=40.3 m2/s2；
1、2为相位观测值；c为光速；N为整周模糊度；tIFB为频率间的偏差；tTGD为卫星码偏差(可由GPS导航电文得到)；LG称为有偏差估计量。TEC也可以类似式(1)用双频伪距观测值ρ1和ρ2进行计算，将相位观测值和伪距观测值分别得出的求解TEC的公式进行组合，可估算出整周模糊度N。设无偏差估计量LGU=LG+N。高度映射函数Eθ求解如下[10]
式中，RE为地球的平均半径，取值为6371 km；θ为卫星高度角。重新排列式(1)，将VTEC=TEC×Eθ代入式(1)，并进一步推得表达式为 其中 式(3)给出了每站的观测量LGU与未知的VTEC、tIFB之间的线性关系。对于每一个测站，在任何的观测历元，只需要接收至少两颗卫星的有效观测数据，依据式(3)利用最小二乘就可以计算出VTEC和tIFB。 2. VTEC异常检验方法
在统计学分析中，平均数是量的平衡，当所描述的数据资料中出现异常值时，利用平均数进行统计分析数据集中趋势就会受到影响产生偏差。中位数是个数的平衡，如果数据资料中异常值个数不超过总个数的一半，中位数就不受影响。由上可知，中位数与平均数这两种用来描述数据资料中心的统计量在分析包含异常值的数据资料时，中位数更加可靠。为了有效地对数据资料进行统计学分析，除了要有描述集中趋势的统计量，还要有描述数据离散程度的统计量。标准差是各数据偏离平均数距离的平均数，同样受到异常值的影响。四分位距(IQR)是一种稳健统计技术中用于表示数据离散程度的一个量。按照标准化正态分布进行统计计算四分位距与标准差的比值为1.34，文中取IQR=1.34 σ，σ为标准差。因此，本文采用中位数及四分位距值进行数据资料的统计分析，以减少异常值对分析结果的影响。
本文构建背景参考值的算法有滑动平均法和非滑动平均法。首先选取时间窗口长度，考虑到VTEC变化具有季节效应，窗口时间长度不宜过长，而且本文主要研究震前的电离层VTEC异常现象，因此选取的用于检验的数据不应受到包含震中及震后等因素的影响，因此本文选择15 d作为时间窗口的长度，重点对震前15 d的VTEC时间序列进行分析。滑动平均法在分析每天的数据时，取其前15 d及当天的同一时刻VTEC值的中位数x 作为背景参考值(如分析16日12∶00时刻的数据时，选取从1—15日及当天16个12∶00时刻的VTEC值)，相关的k倍的四分位距值IQR作为异常检验边界。非滑动平均法只需计算震前15 d的VTEC值，在分析每天的数据时，将这15 d数据的同一时刻的VTEC中位数x作为背景参考值，相关的k倍的四分位距值IQR作为异常检验边界。由此可见，选取相同时间段数据进行异常分析时，非滑动平均法的计算量较小。本文定义x±1.5IQR作为VTEC异常检验上下边界的判断阈值，约为标准差的2倍。为避免检测到由如雷暴、移动性电离层扰动、电离泡等其他地球物理现象引起的短时间扰动，本文认为如果观测到的VTEC值至少在2 h时间内连续超过上限或低于下限时，这样的一天才可以视为有异常现象发生，该检验的置信度为95.4%[5]。下文具体震例分析时，会对两种统计学方法得出的结果作出对比，并探讨两种方法的适用性。
三、震例分析基于上述的理论方法，本文选取了发生在日本的某震例进行分析验证，地震的参数信息参见表 1。GPS数据选取了日本及周边区域的8个IGS测站(aira、chan、ksmv、stk2、suwn、tskb、usud、yssk)，使用震前30 d的GPS观测数据，通过计算得到相应测站上空的电离层VTEC时间序列，单位为TECU。
太阳活动水平的改变、地磁活动异常、天气变化等因素都能够引起电离层的扰动。为了研究地震电离层扰动，检验由地震引起的电离层VTEC异常变化，就必须排除上述背景因素的影响。本文主要查看了2011年3月1日至3月13日(年积日：60—72 d)期间Kp指数和F10.7指数这两个能够反映地磁活动和太阳活动水平的因子(http://www.swpc.noaa.gov)。在年积日60 d、61 d、62 d、69 d、70 d、71 d，地磁场比较活跃；整个3月份太阳活动都达到了中等活跃水平，但年积日66 d、67 d有强烈的太阳活动，并且67 d有太阳质子事件发生。除了上述的几天外，背景因素对电离层VTEC影响较小。
2. VTEC时间序列分析试验绘制了aira、chan、ksmv、stk2、suwn、tskb、usud、yssk 8个GPS观测站在地震前后几天的VTEC时间序列图，限于篇幅，本文选择aira测站震前10 d、地震当天和震后两天的数据为例，部分结果如图 1所示(黑色实线代表VTEC观测值，灰色实线代表上、下限，竖直点线代表发震时刻，填充部分分别表示VTEC值超出上边界及下边界的值，下文同)。从图中可以明显看出，两个测站上空的VTEC值在年积日64 d、67 d都有明显的上升趋势，均升至或超出上边界线，图 1(a)中年积日64 d、67 d及图 1(b)中年积日67 d出现了正异常现象(VTEC值至少在2 h时间内连续超过上限)。由图 1可知，滑动平均法进行异常分析得出的结果(如图 1(a)所示)与由非滑动平均法进行异常分析得出的结果(如图 1(b)所示)总体趋势一致，验证了所得结果的可靠性，但图 1(a)的结果更加明显。由本例可知，采用非滑动平均法能够进行异常整体趋势的考查，但辨别异常情况时，采用滑动平均法分析能够得出更加清晰有效的结果。
 |
| 图 1 aira站上空VTEC时间序列 |
综合对背景因素的分析，年积日67 d的正异常现象很可能与太阳质子流事件有关，为了进一步判断该正异常现象与地震的相关性，本文给出了年积日67 d 06∶00—12∶00UT电离层异常在全球的空间分布，如图 2所示。其中，电离层异常为VTEC值与异常检验上下边界值的差值，VTEC值为IGS提供的全球电离层格网数据(ftp://cddisa.gsfc.nasa.gov/pub/gps/products/ionex)。如果VTEC观测值与上边界值的差值小于0或与下边界值的差值大于0时，则VTEC异常值为0；当VTEC观测值与上边界值的差值大于0时，该差值为电离层正异常；当VTEC观测值与下边界值的差值小于0时，该差值为电离层负异常。由图 2可清晰看出，VTEC异常的分布及变化情况，显著的震区电离层正异常现象主要集中在15°N—35°N，100°E—140°E范围内，并在磁赤道共轭区出现了相似的异常现象，异常位置靠近震中偏赤道的方向，并呈现向磁赤道漂移及向西移动的趋势，异常峰值和幅度在08∶00UT达到最大是太阳活动引起的电离层异常一般表现为更大地理范围的异常，而本文分析的异常现象主要集中在震中附近，且持续时间较长，证实了观测到的电离层异常现象与本次地震密切相关。
 |
| 图 2 3月8日06∶00—12∶00UT时2 h间隔全球VTEC异常分布(单位：TECU) |
此外，年积日64 d的正异常现象基本可以排除太阳活动及地磁活动等背景因素的影响，认为最有可能与本次地震有关。结合上文对年积日67 d电离层正异常的分析，本文认为3月5日及3月8日的电离层正异常现象最有可能与2011年3月11日日本地震相关。
四、结束语基于IGS站提供的GPS数据，采用滑动平均法和非滑动平均法对选取的某地震进行震前VTEC异常分析，发现了很可能与地震相关的较为明显的异常现象，既证明了利用GPS技术可以进行震前VTEC异常的探测，也为地震电离层的研究提供了新的参考资料。此外，本文对背景参考值的选取有较为详细的说明，认为采用中位数及四分位距值进行数据资料的统计分析，可以更加有效地减少异常值对分析结果的影响。通过对两种统计学方法的对比分析，在验证了计算结果可靠性的同时，对两种方法的适用性进行了探讨，本文结果表明，滑动平均法适用性强，能够更有效地进行异常分析，非滑动平均法适用于作整体异常趋势的分析，计算量小，方法简便。
目前，震前电离层异常扰动机理仍然处于探索阶段。至今为止，仍然没有提出很好的物理解释模型，太阳活动及地磁活动对电离层异常的影响还不能准确判断，本文通过对全球电离层异常空间分布进行分析，初步讨论了太阳活动剧烈环境下的电离层异常情况，为确定电离层异常与地震之间的相关性提供了新的证据。
| [1] | 陈鹏. GNSS 电离层层析及震前电离层异常研究[J]. 测绘学报, 2013, 42(3): 474-474. |
| [2] | CALAIS E, MINSTER J B. GPS Detection of Ionospheric TEC Perturbations Following the January 17, 1994, Northirdge Earthquake[J]. Geophysical Research Letters, 1995, 22(9):1045-1048. |
| [3] | LIU J Y, CHEN Y I, CHUO Y J, et al. Variations of Ionospheric Total Electron Content during the Chi-Chi Earthquake[J]. Geophysical Research Letters, 2001, 28(7): 1383-1386. |
| [4] | LIU J Y, CHUO Y J,SHAN S J, et al. Pre-earthquake Ionospheric Anomalies Registered by Continuous GPS TEC Measurements[J]. Annales Geophysicae,EGU,2004,22(5):1585-1593. |
| [5] | YU T, MAO T, WANG Y G, et al. Study of the Ionospheric Anomaly before the Wenchuan Earthquake[J]. Chinese Science Bulletin, 2009, 54(6): 1080-1086. |
| [6] | LIU J Y, CHEN Y I, CHEN C H, et al. Seismoionospheric GPS Total Electron Content Anomalies Observed before the 12 May 2008 Mw7. 9 Wenchuan Earthquake[J]. Journal of Geophysical Research: Space Physics (1978—2012), 2009, 114(A4). DOI:10.1029/2008JA013698. |
| [7] | PULINETS S A, BONDUR V G, TSIDILINA M N, et al. Verification of the Concept of Seismoionospheric Coupling under Quiet Heliogeomagnetic Conditions: Using the Wenchuan (China) Earthquake of May 12, 2008, as an Example[J]. Geomagnetism and Aeronomy, 2010, 50(2): 231-242. |
| [8] | 马新欣. 利用GPS数据研究地震与电离层TEC异常[D]. 成都: 成都理工大学, 2009. |
| [9] | DAUTERMANN T,CALAIS E,HAASE J, et al. Investigation of Ionospheric Electron Content Variations before Earthquakes in Southern California, 2003—2004[J]. Journal of Geophysical Research,2007, 112, B02106.DOI:10.1029/2006JB004447. |
| [10] | 李建勇. 利用GPS数据研究地震电离层TEC异常[D]. 北京: 中国地震局地震预测研究所, 2011. |