引言

在地震电离层现象研究方面，越来越多的科学家发现电、 磁场效应所出现的前兆现象非常显著，而且地震孕育动力学过程并不只是局限于地球岩石圈，大量信息会通过电磁场效应等在大气层、 电离层乃至磁层中有所反映(Hayakawa，1999)．

现今，随着空间卫星观测资料的增多，国内外学者开展了大量的震例研究(Pulinets et al，2005； 吴云等，2005； Parrot et al，2006； Hao et al，2006； 余涛等，2009； 曾中超等，2009； 刘静等，2009； Liu et al，2010； 丁宗华等，2010； 张学民等，2010； Zhang et al，2010)，总结了个别强震前后电离层的异常变化特征． 同时一些专家对震前电离层扰动也进行了统计分析，Pulinets等(1998)对50次M≥5.0地震进行统计分析，发现f_oF₂(F₂层临界频率)在震前具有明显的降低现象； Liu等(2004)分析1999—2002 年台湾M≥6.0 以上地震期间电离层电子含量的变化，发现了地震前1—5 天总电子含量(total electron content，简写为TEC)的异常变化； Nemec等(2008)利用法国DEMETER(detection of electro-magnetic emissions transmitted from earthquake regions)卫星观测的VLF电场数据，对全球9000次M≥4.8地震进行分析后发现，震前0—4小时VLF电磁波会有4—6 dB的降幅； Liu等(2009)对中国1998年5月1日—2008年4月30日35次M≥6.0地震进行统计分析，发现17次M≥6.3地震前3—5天GPS TEC出现明显的下降； Le等(2011)通过分析2002—2010年全球700多次6级以上地震发现，震级越大、 离震中越近、 震源深度越浅，TEC异常越明显．

这些研究成果表明，无论是统计分析，还是个别震例研究，电离层确实在地震前存在明显的扰动现象． 本文基于DEMETER卫星高度上Langmuir探针(instrument sonde de Langmuir，简写为ISL)和等离子体分析仪(instrument d’analyse du plasma，简写为IAP)观测的原位等离子体参量，研究震前异常提取方法，并对2006年1月-2010年4月期间国内6.0级以上、 国外部分7.0级以上的29次地震进行分析，试图寻找660—710 km 高度上电离层对地震的响应，并对其异常参量、 异常特征、 异常范围及异常时间等开展研究．

1.   分析参量

本文主要分析DEMETER卫星等离子体分析仪观测的H⁺数密度(N_H⁺)、 He⁺数密度(N_He⁺)、 O⁺数密度(N_O⁺)、 离子温度(T_i)，以及Langmuir探针观测的电子数密度(N_e)、 离子数密度(N_i)和电子温度(T_e)7种参量． 依据研究结果，现将各种参量特征叙述如下．

1.1   电子数密度

现今研究表明，震前电离层异常主要表现在TEC和f_oF₂参量上，因此对于DEMETER卫星观测的电子数密度是主要的关注参量． 分析发现，对于赤道附近的地震，电子数密度反映比较灵敏，并且会出现大幅度的扰动现象．

1.2   离子数密度

DEMETER卫星的Langmuir探针虽可以观测离子数密度，但由于此载荷转换函数的一些限制，所得的离子数密度数值上并不真实． 电离层本身应近似为等离子体状态，离子数密度与电子数密度应基本相等，但DEMETER卫星实际上观测到的离子数密度要比电子数密度大10倍左右． 不过从分析中发现，离子数密度与电子数密度变化趋势基本一致，因此本参量仅作为参考．

1.3   电子温度

每条轨道上本参量在中低纬度变化幅度较小，最大不会超过1000 K； 但在中高纬度地区变化较强烈，为3000—4000 K，尤其是在磁暴期间表现更为明显． 因此对于中高纬度的震例，本参量异常不作为震前异常，或仅作为参考．

1.4   H⁺和He⁺数密度

这2个参量与太阳活动呈负相关变化，高值区一般出现在20°—50°N及10°—40°S． 在低值期，参量仅有100—300/cm³的密度量，因此小的计数变化都会造成大的相对变化量． 鉴于此，对于低值期的H⁺和He⁺数密度异常是否判断为震前异常，还要参考其它参量的变化．

1.5   O⁺数密度

在DEMETER 卫星高度上，O⁺是最主要的离子成分，其总体变化趋势与电子数密度相同，但其季节特征更为明显，夏季峰值位于北半球，冬季峰值位于南半球．

1.6   离子温度

等离子体分析仪观测的离子温度与Langmuir探针观测的电子温度变化趋势基本一致，中低纬度变化幅度小，中高纬度变化幅度大． 此外，不同于电子温度，离子温度经常出现突跳值，一般情况下主要出现在赤道以及冬季40°N左右和夏季20°S左右的地方．

2.   异常分析方法

对于电离层数据分析现已积累了很多方法，但与研究方法较成熟的GPS TEC和f_oF₂等参量不同，DEMETER卫星数据既随时间变化也随空间变化. 并且由于1颗卫星的限制，其时空分辨率也较低，相邻轨道相距约24°，重访周期13—16天，因此对于其数据分析的方法均在探索阶段． 本研究试图尝试一些分析方法，并寻找到部分震前异常，但每种方法的科学性及可行性还有待于后期继续研究．

2.1   参量随纬度变化

将参量变化按纬度分布绘制曲线图，并将分析数据前一个月的数据按0.1°(ISL观测参量)或1°(IAP 观测参量)重新采样，计算每个纬度上的中值及标准方差，将正负2倍标准差作为上下阈值限定分析数据，寻找异常值．

2005年9月26日秘鲁北部7.6级地震前，9月22—24日在震中东南及西南侧电子数密度升高(图 1)，超出上限值70%—100%，同时发生强烈扰动． 这段时间地磁及太阳活动都很平静，其扰动应与地震有关．

图  1  2005年9月20—25日电子数密度(N_e)随纬度变化图

Figure  1.  Variation of electron density with latitude from 20 to 25 Sep. 2005

下载: 全尺寸图片 幻灯片

2.2   空间插值法

将震中周围的电子数密度数据，以2—3天的时间窗按克里金(Kriging)方法进行插值分析． 张学民等(2008)研究结果表明，DEMETER卫星电子数密度(N_e)具有较好的空间相似性，与选定参考轨道相隔1—2条轨道的数据相关系数大于0.7． 因此，我们可以将相邻的轨道做插值分析，得到电子数密度的二维空间变化图． 此外，震中±20°范围内每天仅有2条升轨(夜间)数据，为减小插值误差，本研究按2—3天形成一幅图像．

采用克里金插值的原因是，此插值方法充分吸收了空间统计的思想，认为任何空间连续性变化的属性是非常不规则的，不能用简单的平滑数学函数进行模拟，但是可以用随机表面给予较恰当的描述； 因此通过确定空间搜索半径，计算这一空间范围所有样本点的自相关和协方差，在此基础上进一步进行插值预测和标准差分析，因而达到比较客观地进行空间插值的效果． 总而言之，克里金方法考虑了数据的空间结构，能够更好地反映数据自身的特性，尤其在数据稀少时其优势更加明显．

2007年5月5日西藏改则6.1级地震前，5月1—2日震中东南侧电子数密度升高，较周围数据数值高20000/cm³； 5月3—4日高值区范围再度扩大，而震后高值异常区消失(图 2)．

图  2  2007年5月5日西藏改则6.1级地震前电子数密度(N_e)空间插值图

红色实心圆表示震中位置

Figure  2.  Space interpolation of electron density before the GaizeM6.1 earthquake in Xizang region on May 5, 2007

Red circle represents epicenter

下载: 全尺寸图片 幻灯片

2.3   空间差值法

本方法是以震中为中心，将经纬度按4°(经度)×2°(纬度)的格子进行划分，计算分析数据前一个月每个格子中的中值作为背景场． 每日观测数据也按划分好的格子求中值，与背景值相减并除以背景场中每个格子的标准差，从而得到每日相对背景场的空间演化图像．

用此方法分析2009年9月29日萨摩亚群岛8.0级地震，发现震前10小时震中东北方向电子数密度增加，涨幅超过10倍标准差，同时电子温度在最北部也呈现下降的趋势，下降幅度为6-7倍标准差，但南部降温现象不明显(图 3)． 9月29日经过本区域的仅有2条升轨数据，离震中最近的轨道未发现扰动现象，而离震中较远的东侧轨道出现异常. 这与Pulinets等(1994)的研究结果相似，异常区域一般不在震中上空而是有一定偏移．

图  3  2009年9月29日电子数密度(a)和电子温度(b)空间差值图

红色实心圆表示震中位置

Figure  3.  Space difference of electron density (a) and electron temperature (b) on Sep. 29, 200

Red circle represents epicenter

下载: 全尺寸图片 幻灯片

2.4   滑动中值法

将震中周围轨道的实测数据按0.1°(ISL观测参量)或1°(IAP观测参量)进行重采样，然后每一纬度点的数据按窗长20轨(大约10天数据)、 步长2轨(每天震中周围夜间数据约为2轨)做滑动中值处理，并取80%和20%的分位点作为上、 下阈值． 为更好地显示异常信息，未超过阈值的结果以“0”值处理，超过上限的用观测值减去上限值，超过下限的用观测值减去下限值． 为反映数据的变化量，用得到的差值除以相应的上、 下限，得到相对变化百分数，计算公式如下：

式中，dev为deviation(差值)的缩写，加“％”表示计算其变化百分数，括号内加“increase”和“decrease”分别表示上升和下降，d_o为实际观测数据，d_u为上限值，d_l为下限值. 如果d_o≤d_u，％dev(increase)=0； d_o≥d_l，％dev(decrease)=0．

2008年5月12日汶川8.0级地震前，对震中±20°范围内的离子温度数据做滑动中值分析，将11°—51°N范围内的数据进行重采样，其原始观测数据、 差值数据及变化百分数时序曲线分别如图 4a，b，c所示． 图中显示，5月9日2条轨道的离子温度均增加800 K左右，最大涨幅110％，异常区域分别位于震中东北2100 km和西侧550 km． Zhao等(2008)研究表明，5月9日GPS TEC也出现异常增加的现象，但位置在震中东南及其共轭区，比本研究离子温度异常偏南．

图  4  2008年5月12日汶川M_S8.0地震前离子温度滑动中值图

(a)原始观测数据; (b)差值数据; (c)变化百分数时序曲线

Figure  4.  Moving median values of ion temperature before the WenchuanM_S8.0 earthquake on May 12, 2008

(a) Original data; (b) Difference data; (c) Time-sequence curve of percentage change

下载: 全尺寸图片 幻灯片

2.5   相关性分析

震前电离层电磁场的变化会引起电子、 离子数密度的变化. 数密度增加，电子、 离子的库仑及弹性碰撞增多，需要将热能转换为动能弥补损失的能量，因此温度会下降； 电磁场变化加热电离层，驱使受热区的带电粒子运动，从而使温度升高区的数密度下降．

无论上述哪种解释，震前电离层的带电粒子数密度及温度都会呈负相关变化，因此本研究对电子数密度及电子温度进行了相关性分析．

2009年3月19日汤加7.9级地震前，3月16日在震中东北1100 km处，电子数密度与电子温度呈现负相关，如图 5中红色框所示，相关系数下降，低于-0.5．

图  5  2009年3月16日电子数密度与电子温度相关系数图

Figure  5.  Correlation coefficient between electron density and electron temperature on Mar. 16, 2009

下载: 全尺寸图片 幻灯片

3.   震例研究

3.1   震例事件

本研究对2006年1月—2010年4月期间的国内6.0级以上、 国外部分7.0级以上的29次地震进行了分析，具体震例信息见表 1．

表  1  震例信息表

Table  1.  Information of earthquakes

下载: 导出CSV 

| 显示表格

3.2   震例分析结果

本研究将上述5种方法应用于29个震例分析中，所选数据为震中周围±20°范围震前10天、 震后2天内的升轨轨道数据． 由于现在还无法定量排除磁暴等因素对电离层的影响，因此对于磁暴日(Dst＜-30，Kp≥4^-)的异常暂不作为震前异常． 由于分析参量及方法较多，因此一条轨道大于等于2个参量出现异常时，本研究才将其视为震前异常．

根据上述原则及各种分析方法的异常阈值，我们发现29个震例中有15次地震出现了震前异常. 每次地震的异常参量、 异常特征、 异常范围、 异常时间及所采用的分析方法见表 2．

表  2  震例统计结果

Table  2.  Case study results of earthquake-related anomalies

下载: 导出CSV 

| 显示表格

4.   讨论与结论

1)地震异常时间基本为震前1—7天，分析中仅有2007年9月12日苏门答腊8.5级地震及2010年1月12日海地7.3级地震前9天及8天出现异常．

2)异常并不出现在震中正上空，多偏向磁赤道地区； 异常有的在震中东侧，有的在震中西侧，范围从数百公里至2000 km左右． 向磁赤道的偏移，应与扰动带电离子沿磁力线的漂移有关． 依据静电场的地震电离层耦合机理(Takeuchi et al，2006)，如果震前地面积累正电荷，那么在向上电场的作用下，带电粒子在E×B的作用下会向西偏移，从而造成震中西侧数密度的会集； 而如果地面积累负电荷，则上述情况相反． 另外，地震多沿断层发生，而孕震过程中可能沿主破裂面出现异常也可能沿次生破裂面出现异常，并且如果断层走向沿东西方向，那么在断层的哪一端出现异常都是可能的． 由分析发现，异常最近距离也在数百公里以上. 其原因一是由于每日震中周围相邻轨道经度间距大约24°左右，空间采样率较低； 二是如上所述，异常较震中会有一定的偏移．

3)Pulinets和Boyarchuk(2004)曾提及，震前等离子体数密度异常幅度一般情况下在±30%左右，部分会超过100%． 本次震例分析发现，异常幅度变化较大，有百分之几十的，也有超过100%的． 这种变化是相对变化，如果本底较低，小的变化量也会造成大的变化幅度，可能与这个原因有关．

4)DEMETER作为一颗太阳同步卫星，观测数据既随时间变化也随空间变化，应用其数据开展震前异常提取方法的研究尚属探索阶段． 本研究将参量随纬度变化、 空间插值、 空间差值、 滑动中值及相关性分析等5种方法应用于29个具体震例研究，提取出了部分地震前的震前异常，说明这些方法的有效性． 但这些分析方法中仍存在很多不足及需要不断完善的地方，后续工作中应继续探索和改进．