1 引言

与地震相关的电离层异常的发现源于Leonard和Barnes^[1]对美国阿拉斯加大地震的研究.近几年在利用电磁卫星观测记录提取中强地震前兆异常已取得了重大进展^[2-10]，但在时空背景趋势研究上较少^[6]，其主要原因之一可能是电磁卫星资料存在明显的时空背景变化特征^[7].在原始图像上，地震前兆短临异常变化常常淹没在卫星轨道运行相关的强背景场中，有可能存在的中强地震相关的电离层微弱异常变化信息则无法显现.我们对电磁卫星的时空背景进行了系统研究，获得了同样的看法^{[11, 12]}.由于卫星运行轨道差异带来的电磁卫星时空背景场对地震前兆信息提取造成的强干扰，目前大多采用卫星轨道对比技术来识别可能存在的地震短临前兆异常变化，但是这样的处理并不能够很好地消除强背景干扰^[2-10].

为了从电磁卫星观测资料的强时空背景变化趋势中有效地挖掘和识别出一些中强地震前电离层可能存在的微弱的地震前兆异常变化信息，需要研究有效的高分辨率的定量识别技术.近年来，Rundle教授在地震活动研究中，提出了具有统计物理意义的图像信息方法(PatternInformatics Method,PI方法)^[13-15]，并应用于地震活动的中长期统计预测研究中，获得了较高的时空分辨能力并取得了较好的空间预测效果^[16-18].此方法所拥有的时空融合能力满足电磁卫星资料时空动态变化特征信息提取的需求.由于PI方法用于地震前兆电磁异常的分析和异常提取的研究很少，因此，通过实际电磁卫星资料处理实践，对PI方法进行必要改进，以满足既有效消除电磁卫星资料的时空背景变化趋势，又能较好地突出短期细微变化信息，从而才实现与地震相关的电离层参量时空动态演化图像异常的定量识别.本文基于改进型PI方法(简称MPI方法)，结合汶川大地震，借助法国DEMETER 电磁卫星的观测资料，通过合理的时空网格划分和数据统计变换，获得电离层参量时空演化图像，进而给出汶川大地震前电磁卫星资料的MPI时空细微变化图像.通过方法测试和特征信号分析，提取到汶川大地震可能存在的前兆孕震信息，证实MPI方法的有效性.

2 PI方法基本原理及改进

Rundle(2000)^[13]提出的PI方法是通过把研究区域进行时空网格划分，对每个网格的地震活动构建其频度时间序列，形成地震活动的时空演化图像，然后计算空间上每个网格里显著地震事件发生的概率，扣除背景概率后可获得发震概率高的网格，即“地震热点".下面我们具体介绍该方法的基本原理^[18]并给出其改进型.

2.1 地震活动研究中PI方法的基本原理

我们对所关注的地震活动区域进行网格化，将其划分为Nb个边长为Δx的网格(图 1a).对每个网格创建一个时间序列Ni(t)，即第i个网格单位时间内发生的震级大于震级下限M_c 的地震数量.在任意一个空间子区(图 1a)的时间进程上，定义背景参考t_b-t₁､短临变化t₁ -t₂､未来预测t₂ -t₃ 三个时间段(图 1b).其中，t₀ 为研究资料的起始时刻，t_b以Δt为步长从t₀ 时刻向t₁ 时刻滑动.一般情况下，取t₂ -t₁ =t₃ -t₂，而预测的结果在未来预测时间段t₂ -t₃ 内有效.在两个时刻t_b和t之间(其中t_b<t)，地震强度函数I_i(t_b，t)定义为从t_b到t的单位时间发生在网格i中震级大于M_c 的平均地震数目

(1)

为了对不同时段的地震强度函数I_i(t_b，t)进行比较，这些地震强度函数应具有相同的统计属性，因此，将地震活动强度函数I_i(t_b，t)进行归一化处理

(2)

其中，〈I_i(t_b，t)〉为所有网格在t时刻地震活动强度函数的平均值，σ(t_b，t)是所有网格在t时刻地震活动强度函数的标准偏差.进而，利用标准化强度函数之差来判断网格i内的地震活动是否异常

(3)

为了减少地震活动随机波动或噪声的影响，需要计算所有标准化地震强度函数ΔI_i(t_b，t₁，t₂)的平均变化

(4)

由(4)式就可以定义和计算未来地震发生在第i个网格内的概率函数P_i(t₀，t₁，t₂)，即地震强度函数平均变化的平方

(5)

最后，可以确定强地震发生在第i个网格内的概率增益为该网格的概率值减去所有网格概率的平均值

(6)

其中，〈P_i(t₀，t₁，t₂)〉为所有网格发震概率的平均值.对于ΔP_i(t₀，t₁，t₂)>0的网格值取常用对数，即

(7)

这里需要设定合理的预测阈值W，以确定热点地区，并利用ROC(Receiver Operating Characteristic,接收者操作特征)方法检验预测时段的地震发生情况^[13-15].在地震活动中，一般对于小于0 的其他空间值则取为0，形成所有的空间网格值，作为空间预测值.最终，对于具有正值的网格(地区)，即ΔP_i(t₀，t₁，t₂)>0网格(地区)命名为热点(危险区)，并假设在t₂ 到t₃ 的时段内，同时考虑摩尔原则^[13-15]，热点区将发生M_f 以上的地震，其中M_f ≥M_c+2.

2.2 用于电磁卫星观测资料处理的MPI方法

由于电磁卫星资料存在强时空背景，提取与地震前兆有关的异常信息非常困难^{[10, 11]}.地震活动中的PI方法需要进行如下改进才能有效地消除时空背景场，提取出有效异常信息:

(1) 创建时间序列Ni(t)时，取消累积时序计算方式以及空间上的摩尔原则(图 1a)^[13-18]，直接采用微分方式计算，提高和突出局部变化信息的敏感度；(2)基于复杂地震前兆异常形态^{[19, 20]}与电磁卫星电离层参量异常变化形态不同^[4-10]，我们取消公式(5)的平方形式，采用指数形式P_i(t₀，t₁，t₂)=e^{ΔI_i(t_b，t₁，t₂)}.在不改变异常形态情况下，只调整数值相对大小，从而可以起到数值的拉伸与放大作用，突出了细微变化成分和微弱异常信息；

(3) 采用动态滑动方式以满足时间进程上的可对比性，同时又能消除电磁卫星的空间不均匀性.具体而言，基于子区上所有时间的数值进行空间归一化方式，代替公式(2)所表达的时间切面的归一化方式，以消除空间背景的差异性.虽然计算公式相同，但归一化方式的改变，其意义与实质不同，这样可以在消除时空变化背景趋势的同时突出短期变化成分；

(4) 取消公式(7)中的常用对数方式和空间值取阈值变换，形成更客观的､连续的､线性变化的前兆空间演化场，每个时刻的空间图像更具有直观的可对比性；

(5) 以固定的时间窗长和空间步长的滑动方式，形成MPI时空图像整体性定量描述参量，从而在消除卫星数据背景变化趋势的同时，又能实现动态刻画时空微弱变化信息的目标，在强背景变化趋势下提取和观察可能存在的微弱异常变化信息.

2.3 基于MPI方法的电磁卫星观测资料处理一般步骤

MPI方法处理电磁卫星资料的具体步骤如下:

(1) 从法国DEMETER 卫星中心提取电磁卫星观测资料；

(2) 确定电磁卫星时空研究对象，进行四维时空网格剖分；

(3) 基于电磁卫星数据的实际演化特点，确定合适的统计量，如平均值或中值；

(4) 按照四维时空网格进行统计量的计算，形成时空演化图像；

(5) 进行合理的时空插值，尽可能恢复和逼近原始的电离层某一参量､某一时段､连续的空间演化场，主要是考虑到一颗电磁卫星运行的周期和轨道的动态分布特点和固有缺陷，有些网格内可能无值，有些网格内平均值可能过小；

(6) 合理设置MPI方法计算处理时需要的一些参量，如窗长和步长；

(7) 对原始时空图像的滑动处理，形成MPI时空动态演化图像；

(8) 对获得的MPI时空动态演化图像的系统分析.

3 MPI方法在汶川大地震前电离层时空异常识别中的应用 3.1 法国DEMETER电磁卫星

法国DEMETER 电磁卫星是全球第一颗专门用于地震预测研究的卫星，其于2004年6月29日发射升空.该卫星运行在太阳同步的圆形轨道上，倾角98.23°，高度660km.卫星的有效载荷包括感应式磁力仪(IMSC)､电场探测仪(ICE)､等离子体分析仪(IAP)､朗缪尔探针(ISL)和高能粒子探测仪(IDP)，其搭载的电磁观测项目非常丰富.一般常用于地震前兆异常分析的有电子浓度､离子浓度､电子温度等电离层观测参量^[21-24].

目前，除了法国的科学家之外，美国､日本､俄罗斯､墨西哥､意大利､中国和印度等国家的相关领域的科学家们，也纷纷开始了DEMETER 电磁卫星观测数据在地震科学方面的应用研究.通过大量中强震例研究表明，地震前的电离层前兆异常，通常出现在大地震前较短的时间内，在比较接近地震震中上空的区域中^[25-31].

3.2 由汶川地区上空电离层观测资料计算时空动态演化图像

北京时间2008年5月12日14时28分在中国四川省汶川县(31.00°N，103.40°E)发生了Ms8.0级地震，震源深度为14km.迄今为止，已有不少学者从各个方面对这次大地震的成因机制､深层过程和动力学响应､震源过程以及震前电离层异常进行了深入的系统研究^{[9, 10, 32-38]}.在此，我们基于MPI技术，针对电子浓度､离子浓度､电子温度3 个电磁卫星观测参量，分析汶川大地震前电离层的时空动态演化特征.

本文采用2008 年1 月1 日至2008 年6 月23日之间的观测资料.图 2是2008年1月电子温度的空间分布图像，可以看出卫星运行轨道所形成的强背景干扰，而很难从中识别出有用的地震信息.在实际资料处理之前，我们将坏的和不稳定的观测资料进行了剔除.根据卫星运行周期和轨道分布特征，对电离层观测资料进行时､空采样.通过多次对比试验，确定时间采样间隔为1 天.在空间方面，以汶川大地震震中为研究中心，采用20°×20°的空间地理范围，纬度和经度采样间隔均为1°.这样，我们可以得到1天×1°×1°的时空网格，进而统计落入每个时空四维格子中电离层参量数值的总和，求出相应的平均值.对无值格子(小于10个数的，根据经验视为无值)，基于卫星运行时空特征，进行具有17天窗长因子的克里金空间插值^[39].通过以上平均值统计和考虑时间窗因素的克里金空间插值重构技术，获得电离层参量连续的时空分布状态图像.进一步利用MPI方法对该时空分布状态图像进行处理，获得MPI时空分布图像，并用于提取与此次地震有关的电离层前兆信息.

3.3 基于不同t₂-t₁ 时间尺度的时空动态滑动MPI图像计算

为了实时跟踪汶川大地震前的电离层参量空间图像随时间动态演化的细节特征，我们分别将获得的电子浓度､离子浓度､电子温度3个电离层参量的时空动态演化图像(图 3-5中的a图)，采用1天滑动步长，通过MPI方法滑动处理，分别得到t₂-t₁ 为15天､20天和25天的MPI时空动态演化图像(图 3-5中的b､c､d图).

3.4 MPI方法处理前后的时空动态演化图像综合对比分析

通过图 3-5动态演化图像的系统对比分析，不难看出以下事实:

(1) 原始的时空动态演化图像中，只能够看到时空背景趋势性､规则性变化(图 3-5 中的a图).在全球各参量时空背景分析中，均可以明显看到时间性和空间性变化规律^{[11, 12]}，而在相对小区域的研究范围之内，电磁卫星运行轨道的周期性和空间性动态变化不明显，也看不到其他有意义的变化信息^[11].电子浓度(图 3 中的a图)与离子浓度(图 4中的a图)的时空变化特征表现基本一致，而两者与电子温度(图 5 中的a图)的时空变化特征差别较大，但都无法用于提取汶川大地震之前可能存在的异常变化信息.

(2) MPI方法在处理时空数据时，为了消除动态变化的背景趋势及随机干扰，通过相同的开始时间，不同的结束时间，分别形成同一子区的2个时间序列，并通过两个时间序列相同部分的相减和t₂-t₁区间剩余部分的求平均，消除可能存在的动态的背景趋势与随机因素.该处理过程不需要任何模型，具有一定的自适应性.因此，通过MPI方法对每幅图像(实际上包含一定的时窗)处理之后(图 3-5中的b､c､d图)，获得的是一个平均滑动结果，已基本消除了明显的时空背景趋势成分，仅仅保留并突出了最基本的细微变化.时空背景消除实际效果较好，达到了预想目标；而其他方法处理，背景趋势总是无法消除的^{[11, 12]}.

(3) 从保留的MPI空间图像随时间的动态演化过程中(图 3-5中的b､c､d图)不难看出:虽然不同参量的MPI时间滑动空间图像存在着一定的差异性，但在汶川大地震爆发前30天，都具有一个共同的特点，即基本上看不到明显的短期变化信息；但在震前10-20天之间，则明显存在持续的､显著的异常变化信息；而在震前10天左右这些异常变化信息又基本上消失了.这一现象有可能与震前地震活动中的“小震平静现象"或岩石破裂声发射试验中存在的岩石大破裂前的平静现象^{[40, 41]}具有一定的相似性，极可能在机制上存在着某种内在的孕震关系，具有发震时间的指示意义；震后则无任何异常变化信息.这样的变化特征也大致与地面震中周围电阻率快速异常变化阶段出现的时间点基本一致^[42]，这可能是大震前岩石变形甚至微破裂产生的可能电磁辐射通过大气不断传递到电离层的作用结果.

(4) 对于3 个电磁卫星观测参量，没有经过MPI处理前，时空图像分布具有一定的差别，至少电子温度与其他参量表现特征不一样；而经过MPI处理后，异常信息形态差别并不大，具有一定的相似性.由此可以说明这些卫星参量时空背景成分虽然不一样，但包含的与地震相关的异常变化信息差别具有很好的时空一致性.因而，从电子浓度､离子浓度､电子温度3个电离层参量中获得的MPI动态滑动时空图像中均可以看到，汶川大地震前存在的显著异常变化和时空分布规律，对于大地震发生的时间和地点具有一定的指示意义.

4 讨论

基于DEMETER 电磁卫星观测资料的时空演化特点，本文对PI方法改进和引入整体滑动平均处理形成MPI方法，并首次用于电磁卫星资料的时空动态演化分析中.通过实际资料处理，结果表明MPI方法能够有效地提取与地震前兆有关的电离层参量的时空特征异常.这里对研究中涉及到的关键性问题进行讨论:

(1) 不同于地震活动研究中主要考虑频度的增强与减弱双边异常方式^[18]，而电磁卫星数据基本属于单边异常方式，可以采用指数形式，因而在本研究中，将公式(5)中平方式用指数式P_i(t₀，t₁，t₂)=e^{ΔI_i(t_b，t₁，t₂)} 进行代替.当然也可以采用其他形式的数值拉伸方式，只要与异常形态一致即可.如果开展其他前兆的MPI演化分析，要根据特定的前兆异常的实际特点，如异常的向上､向下､双边或单边等特征，采用不同的表达方式，以恰当地描述异常的实际变化形态，也才有可能取得较好的背景趋势消除与短临异常提取处理效果；

(2)从本研究可以看到(图 3-5的b､c､d图)，对于不同的t₂-t₁ 时间变化尺度，不同时间滑动的时空图像结果具有细微的差别，说明MPI图像包含的变化信息量与t₂-t₁ 大小的选取具有一定的关系.针对不同的t₂-t₁ 时间尺度，当然也包含背景时段长度t₁-t₀ 及两者的比例关系，通过反复的对比试验，表明MPI图像蕴涵的变化信息差别并不大，且具有一定的稳定性.因此，只要合理选取t₂-t₁ 的大小，就能得到稳定的结果；

(3) MPI方法与PI方法在技术上和表达的意义上都存在一定本质差别.PI方法主要用于地震活动中长期危险性概率预测，其主要依据是当前的地震活动水平对未来的地震危险强度有直接的和内在的关系，这是危险性概率预测的内在基础；MPI方法对应的是地震前兆观测量分析，而前兆观测量与未来强震之间的关系很复杂，具有不确定性，其与地震活动的类比关系是明显不一样的，在电磁卫星处理中是没有概率的意义，只是作为一个时空整体描述参量.通过历史背景的变化规律，消除当前的背景趋势，突出短临变化以提取可能存在的短临异常.当然，这些电磁卫星参量在中强地震之前确实存在一定的前兆异常和预报意义^{[5-8, 10]}，也是本文利用MPI处理卫星资料的基础.不过，对其两者的直接关系还需要进行深入研究和进一步分析，特别对其异常的传递机制更要深入研究，也需要更多的震例来证明这种关系和机制，才有可能提高电磁卫星资料在地震预警中的应有作用.当然有效的地震预警应该是多地震前兆信息综合运用的结果.

(4) 目前，对电磁卫星观测参量的研究很多，但大多采用空间轨道的对比技术，缺乏有效的周期性时空背景场消除能力.本研究从时空背景角度出发，则可以弥补这方面的不足，并就空间地震异常的提取和时空动态演化提出新的分析思路.实际上，在一些不多的研究中也已发现^[6]，通过电磁卫星时空背景的初步分析，也存在类似现象和结果，表明并不是震前几天时空背景异常才能够出现，异常变化早在地震前1个月左右就已产生.因此，对于汶川大地震来说，本文的结果也并非个案.同时该结果与地面的地震前兆异常出现的时间对比，差别也并不大，这极可能与震前平静和岩石的破裂现象都具有一定的内在关系，是这些内在现象的真实外在表现.因此，汶川大地震之前异常出现的时间点特征，即在震前较早时间产生，也具有一定的合理性和内在性，尽管其可能的机制并不完全了解.当然，对产生这种差异性的原因，需要进一步深入研究.

5 结论

基于以上电磁卫星观测资料时空背景的MPI定量计算､实际分析和系统讨论，可以获得以下基本认识:

(1) MPI方法对电磁卫星数据处理的基本核心思想是通过已有的长时间尺度观测资料(t₁-t₀)与当前短时间尺度的资料(t₂-t₁)的合理搭配，基于两者的趋势相关性､同一空间子区中两个不同长度时间序列的非固定模式的自适应相减，在消除当前数据中存在的动态背景趋势的同时，又能突出当前阶段的短期变化异常信息，具有一定的自适应性；

(2) 通过实际电磁卫星资料分析证明，MPI方法既考虑了长期的背景场趋势变化，又突出了近期的短期性变化；既消除电离层因太阳活动等因素产生的空间异常随机变化，又增强了微弱的短期变化结果的显著性､稳定性与客观性.因而，MPI是一种时空综合处理的融合方法，可以用于定量描述电磁卫星时空背景消除与短期变化信息的有效挖掘；

(3) 汶川大地震实际处理结果表明，MPI方法通过全时空分析，消除了背景趋势和随机噪声，获得一定时间窗内具有平均效应的异常分布，可以客观地识别和反映汶川巨大地震之前有可能存在的背景异常变化过程，从而证明了高空电磁卫星客观上确实能够观测到地壳内部地震在孕育过程中岩石破裂产生的可能异常变化^[10].这种现象有可能与地震破裂过程中变化的物理性质有关.当然，对于造成该变化过程的传递方式与物理机制，还有待于今后进一步深入研究.