一种基于小波变换的卫星磁场数据地震前兆异常提取方法与流程

本发明地震监测领域，具体地来讲为一种基于小波变换的卫星磁场数据地震前兆异常提取方法。

背景技术：

地震前兆指地震发生前出现的异常现象，其中包括震源及附近物质发生的地磁、地电、重力等地球物理异常。其观测方法可以包括地面观测和卫星观测两种类型，从全球范围来看地面观测台站分布较为稀疏，而卫星观测弥补了地面观测观测范围有限，且需要铺设大量台站的不足。同时岩石圈-大气层-电离层耦合机制的提出，表明地震前兆产生的异常可以通过岩石圈、大气层传播至电离层，并对其产生影响。因此，利用卫星数据提取地震前兆异常是可行的。磁场在地震孕育的过程中会发生异常变化，卫星测量的磁场数据包括标量总场和三分量矢量数据，矢量数据反映了不同方向磁场的变化情况，包含了更多可以反映地震前的异常情况的信息。

cn106021710a公开了一种基于大气电离层参数的震前卫星轨道异常识别方法，包括原始电离参数的野值剔除及标准化处理；在不同地震震级下，使用不同阈值划分概率密度函数，建立震前卫星异常轨道识别模型；使用综合评价方法得到面向多参数的轨道识别模型等内容。针对大气层电离参数数据量巨大、信息维度高等特点，给出了一种地震前兆卫星轨道异常的识别方法,可用于卫星轨道异常检测和地震前兆识别。

cn107356969a公开了一种基于卫星热红外数据及gis的地震前兆分析方法，采用热红外数据加gis的方式进行地震预报分析，能提高准确率，同时用于分析的数据来源较广，具有较好的研究效果。拓展了地震前兆的研究方向，利用gis方式可将所有对地震发生可能产生影响的量化因素全都纳入空间运算当中，不仅仅局限于地质系数以及地质构造的分析。

mehdiakhoondzadeh等利用中值与四分位法分析了2016年4月16日的厄瓜多尔地震，对2015年11月1日至2016年4月30日的swarm卫星群alpha星的磁场三分量数据分别进行处理，提取出了相应的前兆异常。angelodesantis等分析了2015年4月25日的尼泊尔地震，对地震前后一个月内的东向磁场分量异常按天进行累计，并通过s曲线进行拟合发现存在异常个数快速增加的现象，并将加速度最大处作为前兆异常天。上述技术分别对磁场三分量数据单独进行处理或仅使用其中一个分量进行处理研究，使用的数据中存在较多冗余信息或者使用的数据不全面，未能充分高效的利用三分量数据中含有的有用信息进行地震前兆异常的提取，其异常结果仅由某一分量的数据得到，未能充分全面的反映地震前兆对磁场的影响。

技术实现要素：

本发明所要解决的技术问题在于提供一种基于小波变换的卫星磁场数据地震前兆异常提取方法，通过主成分分析方法对磁场三分量数据进行降维处理，减少冗余数据并保留三分量数据的主要特征，将更多有用信息转换到降维数据中，并通过连续小波变换得到其时频域能量分布，进一步定义并求取轨道能量强度，有效提取降维数据中的异常轨道，排除非震影响因素后，得到地震前兆异常。本发明可以充分利用磁场三分量数据中的有用信息，全面的反映地震前兆对磁场数据的影响，弥补了上述技术单独使用其中某一分量不能充分使用数据有用信息和其异常提取结果不能全面反映地震前兆对磁场数据影响的不足。

本发明是这样实现的，

一种基于小波变换的卫星磁场数据地震前兆异常提取方法，该方法包括：

a、读取卫星的磁场数据，并根据标志位选取有效数据；

b、根据地磁指数选取地磁活动较平静的静磁轨道；

c、通过差分和离散小波变换对静磁轨道的磁场三分量数据进行预处理，去除幅值大和变化缓慢的静磁场部分，得到磁场数据的变化情况；

d、利用主成分分析对预处理后的磁场三分量数据进行降维处理，去除冗余信息的同时保留信号中主要的特征；

e、对降维后的磁场数据进行连续小波变换，通过小波系数定义轨道能量强度，并利用其进行异常轨道提取；

f、排除非地震因素影响导致的异常；

g、输出地震前兆异常的提取结果。

进一步地，步骤a包括：读取卫星观测结果中的磁场数据的磁场三分量bx、by和bz，并根据数据的标志位去除无效数据。

进一步地，步骤b根据地磁指数选取静磁轨道包括：通过地磁指数ap和dst去除受地磁活动干扰较大的轨道，静磁轨道满足的条件为：当前时刻ap＜12，|dst|≤20；前23个小时ap≤32，|dst|≤30。

进一步地，步骤c包括先对选取的磁场三分量数据按轨道分别进行一阶差分处理得到磁场变化量，再分别对差分处理后数据进行离散小波变换，去除其中的低频成分。

进一步地，离散小波变换获取每层高频成分和低频成分的公式为：

其中，j为分解层数，l表示低频成分，h表示高频成分，g为低通滤波器，h为高通滤波器，db为磁场差分数据，每条轨道磁场三分量差分数据与低频成分的残差为预处理结果。

进一步地，步骤d，所述的利用主成分分析卫星磁场数据降维处理包括：

将预处理后每条轨道的磁场三分量数据按照时间序列表示为：

xi＝[xi(1),xi(2)...,xi(m)]，i＝1,2,...,n

其中m为轨道长度，n为数据维度，得到矩阵y的表达式为：

计算矩阵y的协方差矩阵cy(n×n)的元素γuv，计算公式为：

其中，xiu和xiv分别为矩阵y中的第i行u列元素和第i行v列元素；和分别是第u列和第v列元素的均值；

计算协方差矩阵的特征值和特征向量：

cy＝rλr^t其中，λ(n×n)为从大到小排列的特征值对角矩阵，r(1×n)为对应特征值的特征向量，将特征值表示为λ1,λ2,...,λn(λ1＞λ2＞...＞λn)；

利用矩阵r将矩阵y的原始数据进行线性投影得到主成分φ，

φ＝r·y＝[φ1,φ2,...,φn]^t

其中，φ1,φ2,...,φn为第1至第个n主成分(1×m)，前k个主成分的贡献率由其对应特征值计算：

选取使贡献率达到60％以上的前k个主成分代表磁场数据的主要特征，并作为降维处理的结果。

进一步地，步骤e对降维后的磁场数据进行连续小波变换，通过小波系数定义轨道能量强度，并利用其进行异常轨道提取包括：

通过连续小波变换将信号由时域变换到时频域；

在时频域计算信号每一时刻的能量强度；

使用滑动矩形窗对数据进行平均化，以确保异常持续一定的时间，排除突变性异常；

利用该条轨道上的最大平均能量强度反映其异常变化情况。

进一步地，还包括：计算所有轨道能量强度的均值μz和标准差σz，并设阈值为thz＝μz+kk×σz，当轨道的能量强度大于设置的阈值时，该轨道被认为是异常轨道。

本发明与现有技术相比，有益效果在于：

本发明基于小波变换的卫星磁场数据地震前兆异常提取方法，通过差分得到磁场变化量，利用离散小波变换去除周期性的幅值大的静磁场变化趋势得到磁场变化情况；利用主成分分析对三分量磁场数据进行降维处理，提取三分量数据的主要特征，将更多有用信息保留在降维数据中；最后对数据进行连续小波变换得到其时频域能量分布，进一步定义并求取轨道能量强度，突出信号的异常变化，能够得到很好的异常提取效果，并保证提取异常的持续时间，排除突变性异常。本发明通过该过程能够充分利用磁场三分量中包含的信息，有效的提取出具有持续时间的地震前兆异常变化。

附图说明

图1为基于小波变换的卫星磁场数据地震前兆异常提取方法流程图；

图2为厄瓜多尔地震震中、影响区域和研究区域示意图；

图3为卫星磁场三分量数据原始曲线，其中a为bx分量，b为by分量，c为bz分量；

图4a为卫星磁场三分量差分数据与低频成分曲线，其中，a为bx分量，b为by分量，c为bz分量；

图4b为卫星磁场三分量预处理结果曲线，其中，a为bx分量，b为by分量，c为bz分量；

图5为卫星磁场三分量数据第一主成分曲线；

图6a为卫星磁场降维数据小波变换系数结果图；

图6b为卫星磁场降维数据平均能量强度变化曲线；

图7为卫星磁场降维数据轨道能量强度变化曲线；

图8a为2016年2月29日轨道3地震前兆异常结果图；

图8b为2016年3月14日轨道1地震前兆异常结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明基于小波变换对卫星磁场数据进行地震前兆异常的提取，首先读取卫星的磁场数据，并根据标志位选取有效数据；再根据地磁指数选取地磁活动较平静的轨道；通过差分和离散小波变换对静磁轨道的磁场数据进行预处理，去除幅值大和变化缓慢的静磁场部分，得到磁场数据的变化情况；利用主成分分析对预处理后的磁场三分量数据进行降维处理，去除冗余信息的同时保留信号中最主要的特征，使信号中包含更多有用信息；对降维后的磁场数据进行连续小波变换，通过小波系数定义轨道能量强度，突出异常变化，保证提取异常的持续时间，排除突变性异常，并利用其进行异常轨道提取；排除非地震因素影响导致的异常；输出地震前兆异常的提取结果。本发明能够有效的利用小波变换对降维后的卫星磁场数据进行地震前兆异常进行提取。

参见图1所示，一种基于小波变换的卫星磁场数据地震前兆异常提取方法，包括以下步骤：

a、读取卫星的磁场数据，并根据标志位选取有效数据；

b、根据地磁指数选取地磁活动较平静的静磁轨道；

c、通过差分和离散小波变换对静磁轨道的磁场三分量数据进行预处理，去除幅值大和变化缓慢的静磁场部分，得到磁场数据的变化情况；

d、利用主成分分析对预处理后的磁场三分量数据进行降维处理，去除冗余信息的同时保留信号中主要的特征；

e、对降维后的磁场数据进行连续小波变换，通过小波系数定义轨道能量强度，并利用其进行异常轨道提取；

f、排除非地震因素影响导致的异常

g、输出地震前兆异常的提取结果；

步骤a，读取卫星的磁场数据，并根据标志位选取有效数据，是读取卫星观测结果中的磁场三分量bx、by和bz，分别指的是：北向、东向和垂直分量，并根据数据的标志位去除无效数据，仅使用有效的磁场数据，标志位是表征磁场数据质量情况的参数，可以通过标志位判断仪器测量时是否出现问题，去除由于仪器测量有误导致的无效数据，为现有的方法，此处不再赘述。由于测量数据的时间和空间均会随卫星飞行发生变化，故按照轨道对数据进行处理和异常提取。

步骤b，根据地磁指数选取静磁轨道，是通过地磁指数ap和dst去除受地磁活动干扰较大的轨道，避免提取的异常由地磁活动引起。地磁指数反映了一段时间内的地磁活动强度，静磁轨道需要满足的条件为：当前时刻ap＜12，|dst|≤20；前23个小时ap≤32，|dst|≤30。

步骤c，利用差分和离散小波变换卫星磁场数据预处理，是先对选取的磁场三分量数据按轨道分别进行一阶差分处理得到磁场变化量，再分别对差分数据进行离散小波变换，去除其中的低频成分，突出数据的变化情况。卫星测量磁场为地磁场，其中包含主磁场、岩石圈磁场、变化磁场和感应磁场四部分。地磁场幅值为10⁴nt量级，而地震前兆引起的变化主要表现在岩石圈磁场中，只有几nt甚至更小，相对于总体幅值来说非常微弱，直接对原始测量磁场进行操作则无法提取异常。故首先通过一阶差分得到磁场数据的变化量，由于测量的地磁场本身存在周期性变化，差分之后的结果依旧存在相对于异常变化较大的低频成分。利用离散小波变换对差分后的数据进行分解，通过去除数据中的低频成分消除周期性变化的影响，最后得到磁场的变化情况，其结果可用于主成分分析等后续操作。

离散小波变换获取每层高频成分和低频成分的公式为：

其中，j为分解层数，l表示低频成分，h表示高频成分，g为低通滤波器，h为高通滤波器，db为磁场差分数据。每条轨道磁场三分量差分数据与低频成分的残差为预处理结果。

步骤d，利用主成分分析卫星磁场数据降维处理，是通过求取三分量数据的主成分对数据进行降维处理，并尽可能保留多维数据的主要特征。地震前兆异常产生的磁场变化应该在磁场的不同分量上都有反映，利用主成分分析对磁场三分量数据进行降维处理，可以去除冗余信息使后续处理的计算过程得到简化，同时保留三分量磁场中最主要的特征。降维结果同时包含了三个分量数据的主要特征，对其进行分析处理比单独对某个分量进行处理结果更全面，能够更好地反应地震前兆信息。具体包括：

将预处理后每条轨道的磁场三分量数据按照时间序列表示为：

xi＝[xi(1),xi(2)...,xi(m)]，i＝1,2,...,n(2)

其中m为轨道长度，n为数据维度。则可得到矩阵y的表达式为：

计算矩阵y的协方差矩阵cy(n×n)的元素γuv，计算公式为：

其中，xiu和xiv分别为矩阵y中的第i行u列元素和第i行v列元素；和分别是第u列和第v列元素的均值。

计算协方差矩阵的特征值和特征向量：

cy＝rλr^t(5)

其中，λ(n×n)为从大到小排列的特征值对角矩阵，r(1×n)为对应特征值的特征向量，将特征值表示为λ1,λ2,...,λn(λ1＞λ2＞...＞λn)。

利用矩阵r将原始数据进行线性投影可得到主成分φ。

φ＝r·y＝[φ1,φ2,...,φn]^t(6)

其中，φ1,φ2,...,φn为第1至第个n主成分(1×m)，前k个主成分的贡献率由其对应特征值计算：

选取使贡献率达到60％以上的前k个主成分代表磁场数据的主要特征，并作为降维处理的结果。

步骤e，对卫星磁场降维数据进行连续小波变换，并利用得到的小波系数求取轨道能量强度对异常轨道进行提取。降维处理得到的结果包含了三分量的数据信息，但是其中的异常并不明显，直接通过时域进行异常提取存在困难。故通过连续小波变换将信号由时域变换到时频域，在时频域计算信号每一时刻的能量强度，突出其异常变化，并使用滑动矩形窗对数据进行的平均化，以确保异常持续一定的时间，排除突变性异常，最后利用该条轨道上的最大平均能量强度反映其异常变化情况，此时进行异常提取会得到较为理想的效果。通过连续小波变换并计算轨道能量强度，可以非常有效的提取出降维数据中的异常轨道。

连续小波变换的公式为：

其中a为尺度因子，b为平移因子，ψ为母小波函数，f(t)为待处理数据。每条轨道磁场降维数据的连续小波变换结果为一矩阵：

其中m为轨道的长度，n由选取的尺度因子a决定。

根据能量守恒定理可以从频域计算信号的能量，故可以利用时频域的小波系数计算能量强度表征磁场降维数据每一时刻的能量，以突出磁场的异常变化。由小波系数计算能量强度的公式为：

则得到每条轨道的能量强度数据序列为：ze＝[ze(1),ze(2),...,ze(m)]。

为确保异常持续一定的时间，排除突变性异常，使用窗长为l，步长为s的矩形窗在轨道上滑动，求取每个窗内能量强度的均值，并取整个轨道所有均值的最大值为该条轨道的能量强度zemax，代表该条轨道具有的能量强度，得到所有轨道能量强度的序列为：

zemax＝[zemax(1),zemax(2),...,zemax(p)](10)

其中p为静磁轨道个数。

计算所有轨道能量强度zemax的均值μz和标准差σz，并设阈值为thz＝μz+kk×σz，其中kk为经验系数，当轨道的能量强度大于设置的阈值时，该轨道被认为是异常轨道。

步骤f，电离层测量参数受到多种因素的影响，例如季节、地理位置、热层风、电离层扰动以及地磁扰动。为进一步确认提取异常由地震引起，须排除其他因素影响导致的异常。步骤b已经通过地磁指数去除受地磁扰动影响的轨道，由于非震因素产生的影响属于全球性干扰，故通过去除同一轨道上既在地震影响区域内出现又在地震影响区域外的异常，以及出现在地震影响区域内的但延伸到影响区域外的异常或仅出现在地震影响区域外的异常，排除非震因素影响导致的异常。

步骤g，通过非震因素影响的排除，将最终得到的异常轨道作为地震前兆异常轨道进行输出。

下面结合实施实例对本发明进行详细说明。

针对2016年4月16日发生的7.8级厄瓜多尔地震，以swarm卫星群alpha星的磁场三分量数据(1hz)为例。根据dobrovolsky's公式r＝10^0.43m(m为地震的震级)选取厄瓜多尔地震的研究区域，以震中为中心以4518.8km为边长的正方形区域为研究区域，选取震前90天至震后30天，即2016年1月17日至5月16日为研究的时间范围。厄瓜多尔地震震中、地震影响区域以及研究区域如图2所示，其中五角星为震中，圆形区域为地震影响区域，正方形区域为研究区域。

a、录入swarm卫星群alpha星2016年1月17日至5月16日磁场三分量数据(1hz)，数据按北向、东向和垂直分量分别记为bx、by和bz。挑选出经过研究区域的轨道，并根据数据中的标志位去除无效数据，仅使用有效的磁场三分量数据。由于测量数据的时间和空间均会随卫星飞行发生变化，故按照轨道对数据进行处理和异常提取。以2016年4月1日轨道1为例，磁场三分量的原始曲线如图3所示，可看出磁场原始数据幅值非常大，无法观察到其微弱的变化，直接对其进行操作无法提取异常。

b、利用研究时间范围内的地磁指数ap和dst，按照静磁条件：当前时刻ap＜12，|dst|≤20；前23个小时ap≤32，|dst|≤30。去除受地磁活动干扰较大的轨道，避免提取的异常由地磁活动造成，得到静磁条件下226条轨道的磁场三分量数据。

c、对选取的磁场三分量数据按轨道分别进行一阶差分处理得到磁场变化量，再分别对差分数据进行离散小波变换，去除其中的低频成分，突出数据的变化情况。离散小波变换获取每层高频成分和低频成分的公式为：

针对磁场数据本身的特点，取分解层数j为1至3，并选择db4小波基的低通滤波器和高通滤波器作为g和h，l表示低频成分，h表示高频成分，db是磁场差分数据。将计算出的第3层的低频成分db3,l作为去除项，每条轨道磁场三分量的差分数据与第3层低频成分的残差为预处理结果。

以2016年4月1日轨道1为例，其差分数据和低频成分分别如图4a中的实线和虚线所示，可看出差分数据中依旧存在幅值较大的低频周期成分，需要进行进一步去除，离散小波变换计算的低频成分与差分数据对应良好，相减后不会产生伪异常。利用离散小波变换去除低频成分可以有效去除静磁场背景，突出磁场本身的变化情况。差分数据减去低频成分的残差曲线如图4b所示，由结果图可知该步处理可以很好的得到磁场的变化情况。

d、将预处理后每条轨道的磁场数据按时间序列表示为

xi＝[xi(1),xi(2)...,xi(m)]，i＝1,2,...,n(2)

其中m为轨道长度，n为数据纬度取3。则可得到矩阵y的表达式为：

计算矩阵y的协方差矩阵cy(3×3)的元素γuv，计算公式为：

其中，xiu和xiv分别为矩阵y中第i行的第u列和第v列数据；和分别是第u列和第v列数据的均值。

计算协方差矩阵的特征值和特征向量：

cy＝rλr^t(5)

式中，λ(3×3)为从大到小排列的特征值对角矩阵，r(1×3)是对应特征值的特征向量，将特征值表示为λ1,λ2,λ3(λ1＞λ2＞λ3)。

用矩阵r将原始数据进行线性投影可得到主成分φ。

φ＝r·y＝[φ1,φ2,φ3]^t(6)

其中，φ1,φ2,φ3为第1至第3个主成分(1×m)，前k个主成分的贡献率由其对应特征值计算：

经计算发现第一主成分的贡献率ζ1＞60％，故选取第一主成分φ1代表磁场数据的主要特征，并作为降维处理的结果。

以2016年4月1日轨道1为例，其降维结果如图5所示，与图4b对比可以看出，主成分降维后的数据保留了信号中最主要的特征，且去除了冗余的信息使后续处理的计算过程得到简化，有利于后续异常的识别和提取。

e、对磁场降维数据进行连续小波变换，将信号由时域变换到时频域。连续小波变换的公式为：

根据磁场降维数据的特性，令尺度因子a为512，ψ为db4小波基函数，f(t)为φ1。每条轨道降维数据的连续小波变换结果为一矩阵：

其中m为轨道的长度，n由尺度因子a决定为512。以2016年4月1日轨道1为例，磁场降维数据的小波变换结果如图6a所示，可以得到该轨道数据在时频域的能量分布情况，用于求取能量强度表征每一时刻的能量。

根据能量守恒定理可以从频域计算信号的能量，故利用时频域的小波系数计算能量强度表征磁场降维数据每一时刻的能量，以突出磁场的异常变化。由小波系数计算能量强度的公式为：

则得到每条轨道的能量强度数据序列为：ze＝[ze(1),ze(2),...,ze(m)]。

为确保异常持续一定的时间，排除突变性异常，使用窗长为l＝20，步长为s＝20的矩形窗在轨道上滑动，求取每个窗内能量强度的均值。以2016年4月1日轨道1为例，其平均能量强度曲线如图6b所示，对比图5的降维数据可知，求取数据的能量强度可以突出异常变化，求取均值可以反映磁场能量在该段时间内的平均强度，消除突变异常的影响。

取整个轨道所有均值的最大值为该条轨道的能量强度zemax，代表该条轨道具有的能量强度，得到所有轨道能量强度的序列为：

zemax＝[zemax(1),zemax(2),...,zemax(226)](10)

其中p为静磁轨道个数为226。

计算所有轨道能量强度zemax的均值μz和标准差σz，取kk＝3设阈值为thz＝μz+3×σz，当轨道的特征值大于设置的阈值时，该轨道被认为是异常轨道。异常提取结果如图7所示，其中竖直虚线为地震当天，水平点划线自下而上分别为均值μz和阈值thz＝μz+3×σz，星号为超出阈值的异常轨道。图中震前超出阈值的轨道分别为2016年1月25日轨道3、2月6日轨道3、2月29日轨道3和3月14日轨道1。从结果图可看出通过连续小波变换并计算轨道能量强度可以突出异常，有效提取所有静磁轨道中的异常轨道。

f，电离层测量参数受到多种因素的影响，例如季节、地理位置、热层风、电离层扰动以及地磁扰动。为进一步确认提取异常由地震引起，须排除其他因素影响导致的异常。步骤b已经通过地磁指数去除受地磁扰动影响的轨道，由于非震因素产生的影响属于全球性干扰，故通过去除同一轨道上既在地震影响区域内出现又在地震影响区域外的异常2016年1月25日轨道3，以及出现在地震影响区域内的但延伸到影响区域外的异常或仅出现在地震影响区域外的异常2016年2月6日轨道3，排除非震因素影响导致的异常。

g，通过非震因素影响的排除，最后只有2016年2月29日轨道3和2016年3月14日轨道1为异常，如图8a和8b所示，由提取结果可知通过本发明的处理过程可以去磁场数据中的幅值大且变化缓慢的静磁场得到磁场数据的变化情况，并充分利用磁场三分量数据包含的信息，有效的突出异常变化，保证提取异常的持续时间，排除突变性异常，得到良好的异常提取结果。将2月29日轨道3和3月14日轨道1作为厄瓜多尔地震前兆异常轨道进行输出。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。