阅读说明：本技术 一种基于非负矩阵分解的卫星磁场数据地震异常检测方法 (Satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization ) 是由 朱凯光 樊蒙璇 李凯艳 贺小丹 池成全 于紫凝 孙慧慧 于 2019-08-26 设计创作，主要内容包括：本发明为一种基于非负矩阵分解的卫星磁场数据地震异常检测方法,包括：根据标志位去除无效数据减去CHAOS-6磁场模型数据,对结果求一阶差分得到差分数据；对差分数据进行短时傅里叶变换构造磁场数据对应的非负时频幅值矩阵；利用非负矩阵分解方法对非负时频幅值矩阵进行分解,将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离；根据能量比选出因地震产生的局部影响分量,通过超限率的方法对该分量进行异常轨道判断；对每日异常轨道个数进行累计,通过其偏离背景拟合直线的程度检测地震异常。本发明可以保留并利用测量得到的所有数据对地震进行研究,同时得到与地震活动更为相关的分量有效的进行地震异常检测。(The invention relates to a satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization, which comprises the following steps: subtracting CHAOS-6 magnetic field model data from the removal invalid data according to the zone bit, and solving a first-order difference of the result to obtain differential data; carrying out short-time Fourier transform on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data; decomposing a non-negative time-frequency amplitude matrix by using a non-negative matrix decomposition method, and separating a local influence component generated due to an earthquake and a global influence component generated due to solar activity and geomagnetic activity; selecting local influence components generated by earthquake according to the energy ratio, and judging abnormal orbits of the components by an overrun method; and accumulating the number of the abnormal orbits every day, and detecting earthquake abnormity according to the deviation degree of the abnormal orbits from the background fitting straight line. The invention can reserve and utilize all data obtained by measurement to study the earthquake, and simultaneously obtains components more relevant to earthquake activities to effectively detect the earthquake abnormity.)

一种基于非负矩阵分解的卫星磁场数据地震异常检测方法

技术领域

本发明属于地震监测领域，具体地来讲为一种基于非负矩阵分解的卫星磁场数据地震异常检测方法。

背景技术

地震异常指地震发生前后出现的异常现象，其中包括震源附近的地磁、地电、重力等地球物理异常。研究人员往往通过测量上述物理量，并研究其异常情况实现对进行地震异常的检测。测量方式主要包括地面台站和卫星监测两种，卫星监测方式因其测量范围广，且已被证明可以检测到地震异常现象而被广泛研究。其中卫星测量的磁场数据被大量研究，目前主要通过选用夜间数据同时去除其中太阳活动指数和地磁指数较高的数据排除太阳和地磁活动对数据的影响，并直接使用未剔除数据进行异常检测。

现提出一种基于非负矩阵分解的卫星磁场数据地震异常检测方法。非负矩阵分解方法可以将任意一高维非负矩阵分解为两个低维的非负矩阵W和H，使得 W与H的乘积近似的逼近原高维矩阵V，其中W称为基矩阵，H称为权重矩阵。基矩阵W可以看成原矩阵V的一组基底，其列向量表征了数据的特征；而权重矩阵H为矩阵V在基向量矩阵空间的投影，其行向量表示基矩阵W对应列向量随时间变化的幅值。其中原始矩阵V中的元素均为非负数，其分解矩阵W和 H中的元素也均为非负数，这使得分解结果具有真实的物理意义。基于该特性将卫星Y分量磁场数据预处理后进行短时傅里叶变换得到的非负时频矩阵作为原始矩阵进行分解，可得到数据的几个时频分布特征及其在时域上的幅值。将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离，不需要去除太阳和地磁活动较强的天的数据，同时可以得到与地震活动更为相关的分量有效的进行地震异常检测。

CN101793972A公开了一种强地震短临预测卫星热红外亮温异常技术，主要利用极轨卫星、静止气象卫星和卫星星座微波进行观测，得到卫星热红外云图，并通过大气模型进行校正获得实际温度值，根据地形和气团形态排除掉地形、天气因素导致的干扰信息，最后得到地震前兆热红外温度异常特性及其演变规律对地震三要素短临预报的关系。

泽仁志玛等利用DEMETER卫星记录的变化磁场数据统计研究了2005-2009 年北半球7级以上强震前后空间磁场的扰动特征。根据地震震级和发生时间选取的研究空间和时间范围，选用夜间数据去除太阳的影响，剔除地磁活动较强的卫星观测数据，使用剩余数据对地震进行研究。利用同期观测数据构建背景场，通过磁场在地震时段相对于背景场的扰动幅度观测地震异常。

Mehdi Akhoondzadeh等研究了2016年4月16日的厄瓜多尔地震，研究了2015 年11月1日至2016年4月30日的Swarm卫星群Alpha星的磁场三分量数据，并计算每日数据的中值，通过所有中值的四分位数判断该天是否为异常天，将结果中地磁活动影响较大的天去除，认为剩余异常天可能由地震引起。

上述技术使用传统的方式去除太阳和地磁活动对数据的影响，需要剔除部分数据并直接对剩余数据进行异常检测.

发明内容

本发明所要解决的技术问题在于提供一种基于非负矩阵分解的卫星磁场数据地震异常检测方法，可以保留并利用测量得到的所有数据对地震进行研究，同时得到与地震活动更为相关的分量有效的进行地震异常检测。

本发明是这样实现的，

一种基于非负矩阵分解的卫星磁场数据地震异常检测方法，该方法包括：

a、读取卫星的磁场数据，并根据标志位去除无效数据后作为原始数据；

b、原始数据减去CHAOS-6磁场模型数据，并对结果求一阶差分得到差分数据；

c、对步骤b的差分数据进行短时傅里叶变换构造磁场数据对应的非负时频幅值矩阵；

d、利用非负矩阵分解方法对磁场数据的非负时频幅值矩阵进行分解，将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离；

e、根据能量比选出因地震产生的局部影响分量，并通过超限率的方法对该分量进行异常轨道判断；

f、对每日异常轨道个数进行累计，并通过其偏离背景拟合直线的程度检测地震异常；

g、输出偏离程度指数曲线结果图。

进一步地，步骤a，读取卫星的磁场数据包括北向X分量B_x，东向Y分量B_y和垂直Z分量B_z三个分量的数据。进一步地，步骤b，具体为：用原始数据的Y 分量磁场数据减去CHAOS-6磁场模型计算的背景场，得到残差数据，根据研究Y分量磁场数据最容易受到岩石圈变化的影响反映地震异常。此处选择原始数据的Y分量磁场数据。

再对残差数据取一阶差分去除剩余小背景值得到磁场的数据的变化量，采用式(1)计算：

i＝1，2，3，...L-1

其中B_y为原始Y分量磁场数据，B_m为CHAOS-6磁场模型计算的背景场， t_i为第i个数据对应的时刻，L为每条轨道数据的长度，

为减模型求一阶差分后的预处理结果。

进一步地，步骤c，对差分数据进行短时傅里叶变换构造磁场数据对应的非负时频幅值矩阵包括：

通过短时傅里叶变换对预处理后的磁场数据进行时频变换得到包括幅值信息和相位信息的时频矩阵；

取幅值构成的时频幅值矩阵作为后续非负矩阵分解的输入矩阵。

进一步地，步骤d，具体包括：通过非负矩阵分解方法将任意一高维非负时频幅值矩阵分解为两个低维的非负矩阵W和H，使得W与H的乘积近似的逼近原高维矩阵V，其中W称为基矩阵，H称为权重矩阵，得到数据的时频分布特征矩阵及时频分布特征矩阵在时域上的幅值；

非负矩阵分解方法的数学模型为：

V≈WH (6)

其中V∈R^≥0为一M×N的非负矩阵，V∈R^≥0为一M×R的非负矩阵， V∈R^≥0为一R×N的非负矩阵，分解特征个数R的选择满足R＜MN/(M+N)；

基于K-L散度制定非负矩阵分解的目标函数，其中基于K-L散度的目标函数为：

i＝1，2，3...M

j＝1，2，3，...N

其中V_ij表示矩阵V的第i行、j列元素，[WH]_ij表示基矩阵W和权重矩阵 H相乘得到矩阵的第i行、j列元素；

利用最小行列式准则对权重矩阵H进行约束，通过权重矩阵H的向量张成空间的体积最小进行约束，权重矩阵H的向量张成空间的体积的求解公式为：

则添加最小行列式约束的K-L散度非负矩阵分解目标函数为：

i＝1，2，3，...M

j＝1，2，3，...N

其中det表示求矩阵的行列式，α为平衡参数；

对目标函数采用乘性迭代方法进行优化求解，通过交替更新变量求得最优解，利用梯度下降法计算更新公式：

其中令μ_w和μ_h分别为：

最终更新迭代公式为：

其中

表示矩阵元素的点乘；

标准化公式为：

设置迭代停止条件为目标函数小于设置的阈值Th，同时设置迭代次数的上限为N_iter，若迭代N_iter次后目标函数的值依旧不小于所设置的阈值则强行停止迭代；

得到每条轨道分离得到的基矩阵W为：

权重矩阵H

为：

进一步地，步骤e，根据能量比选出因地震产生的局部影响分量，通过超限率的方法对该分量进行异常轨道判断，具体为：

通过计算非负时频幅值矩阵分解得到的权重矩阵H的各个行向量中研究区域内能量与整条轨道能量之比，选择能量比最大的分量HH作为因地震产生的局部影响分量；

再选用反应数据能量平均水平的均方根作为阈值参数，根据最大的分量 HH的值是否大于整条轨道均方根的k倍进行异常点判断，认为异常点只出现在研究区域内，且对应标志位显示异常点不是由于已知电离层活动引起的轨道为异常轨道。

进一步地，步骤f，对每日异常轨道个数进行累计，并通过其偏离背景拟合直线的程度检测地震异常具体包括：

通过离差标准化对每日异常轨道与每日所有轨道个数进行归一化处理，

计算得到每日异常轨道个数进行累计归一化结果为：

MN_ano＝[MN_ano(1)，MN_ano(2)，MN_ano(3)，...MN_ano(d)] (26)

每日所有轨道累计归一化结果通过最小二乘作直线拟合，得到背景拟合直线；

计算每日所有轨道个数进行累计归一化结果为：

MN_nom＝[MN_nom(1)，MN_nom(2)，MN_nom(3)，...MN_nom(d)] (27)

对MN_nom进行最小二乘线性拟合，拟合直线为：

y＝a₀+a₁x (28)

其中参数a₀和a₁估计求解公式为：

其中x为[1，2，3...d]，d为天数，x_i为x中的第i个元素；y为MN_nom，y_i为 MN_nom中的第i个元素；x_iy_i为x中的第i个元素与MN_nom中的第i个元素相乘；经计算得到直线方程参数

和

通过每日异常轨道个数累计曲线偏离背景拟合直线的程度来显示地震变化过程，偏离线性拟合程度计算公式为：

i＝1，2，3，...d

其中Y_i为MN_ano中的第i个元素；

ρ是表示每日异常轨道累计曲线偏离背景拟合直线程度的指数，越大表示数据因受地震影响产生的异常越大。

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

本发明基于非负矩阵分解的卫星磁场数据地震异常检测方法，通过减去 CHAOS-6磁场模型数据并进行一阶差分，去除磁场数据背景得到磁场数据变化量；利用短时傅里叶变换构造磁场数据对应的非负时频幅值矩阵；通过非负矩阵分解对构造的时频幅值矩阵进行分解，根据数据的时频特性自适应分解出相应的频率分布特征矩阵和权重矩阵，将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离；并通过能量比选取因地震产生的局部影响分量进行研究，使用超限率的方法对该分量进行异常轨道判断；最后对每日异常轨道个数进行累计，并通过其偏离线性拟合的程度检测地震异常。本发明通过非负矩阵分解对数据的非负时频幅值矩阵进行分解，根据磁场数据的时频特性自适应分离出太阳和地磁活动对数据的干扰，得到其频率分布特征和时间变化幅值，其分解结果具有真实物理意义。可以保留并利用测量得到的所有数据对地震进行研究，使用与地震活动更为相关的分量有效的进行地震异常检测，弥补了上述技术去除太阳以及地磁活动较强时的数据，减少数据使用量，并使用内剩余数据直接进行异常检测的不足。

附图说明

图1为基于非负矩阵分解的卫星磁场数据地震异常检测方法流程图；

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

图3为卫星磁场Y分量数据原始曲线；

图4为卫星磁场Y分量原始数据与CHAOS-6模型数据曲线，残差曲线以及一阶差分曲线；图4(1)中实线为原始数据曲线，点划线为模型数据曲线，由于原始数据幅值较大且模型数据与其相差很小，图4(2)将图4(1)其进行部分放大，图4(3)为原始数据减模型数据得到的残差数据，图4(4)为残差数据的一阶差分数据；

图5为卫星Y分量磁场差分数据短时傅里叶变换时频幅值矩阵结果图；

图6为非负矩阵分解频率分布特征矩阵W结果图；

图7为非负矩阵分解权重矩阵H结果图；

图8为非负矩阵分解结果时域重构结果图，图8(1)为一阶差分数据，在研究区域内外均存在明显奇异值；图8(2)为W1和H1对应的时域分量，只在研究区域内存在明显奇异值；图8(3)为W2和H2对应的时域分量，不存在明显奇异值，在整条轨道内均匀分布；图8(4)为W3和H3对应的时域分量；

图9为异常轨道判断结果图，图9(1)为4月9日轨道3，没有点超出阈值，故不被判定为异常轨道；图9(2)为4月10日轨道1，存在超出阈值的点，且仅在研究区域内出现，同时通过检查其标志位没有显示已知电离层活动，可被判定为异常轨道；

图10为每日异常轨道累计与背景拟合直线结果图；

图11为每日异常轨道累计曲线偏离背景拟合直线程度指数曲线。

具体实施方式

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

参见图1所示，一种基于非负矩阵分解的卫星磁场数据地震异常检测方法，，该方法包括：

a、读取卫星的磁场数据，并根据标志位去除无效数据；

b、原始数据是去除无效数据之后的数据减去CHAOS-6磁场模型数据，并对结果求一阶差分；

c、对差分数据进行短时傅里叶变换构造磁场数据对应的非负时频幅值矩阵；

d、利用非负矩阵分解方法对磁场数据的时频幅值矩阵是的是非负时频矩阵进行分解，将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离；

e、根据能量比选出因地震产生的局部影响分量，并通过超限率的方法对该分量进行异常轨道判断；

f、对每日异常轨道个数进行累计，并通过其偏离背景拟合直线的程度检测地震异常；

g、输出偏离程度指数曲线结果图。

步骤a，读取卫星的磁场数据包括北向X分量B_x，东向Y分量B_y和垂直Z 分量B_z三个分量的数据。

步骤b，原始数据减去CHAOS-6磁场模型数据，并对结果求一阶差分，是用原始Y分量磁场数据减去CHAOS-6磁场模型计算的背景场，首先利用原始数据减去CHAOS-6磁场模型计算得到的背景场，去除地磁场幅值中较大的背景，得到残差数据，再对残差数据取一阶差分去除剩余较小背景值得到磁场的数据的变化量，其结果可用于短时傅里叶变换等后续操作。

i＝1，2，3，...L-1

其中B_y为原始Y分量磁场数据，B_m为CHAOS-6磁场模型计算的背景场， t_i为第i个数据对应的时刻，L为每条轨道数据的长度，dB为减模型求一阶差分后的预处理结果。

步骤c，对差分数据进行短时傅里叶变换构造磁场数据对应的非负时频幅值矩阵包括：

通过短时傅里叶变换对预处理后的磁场数据进行时频变换得到其时频矩阵；

取其幅值构成的时频幅值矩阵作为后续非负矩阵分解的输入矩阵。非负矩阵分解要求输入矩阵中的内每一个元素均为非负值，由于Y分量磁场数据为单通道数据无法构成矩阵，且预处理结果中存在负值，故通过短时傅里叶变换对数据进行时频变换，构造非负时频幅值矩阵，作为非负矩阵分解的输入矩阵。且短时傅里叶变换具有相关的物理意义，其结果反映的是数据的时频特性，即数据频率的幅值和相位分布随时间的变化，此处利用数据频率幅值分布随时间的变化，构成非负矩阵进行分解，后续可根据数据的时频幅值特性进行分解。

设预处理后每条轨道数据的时间序列为：

X＝[x₁，x₂，x₃，x₄，...，x_L-1] (2)

其中L为每条轨道数据的长度。

短时傅里叶变换的离散形式为：

n＝0，1，2，3，...N-1

k＝0，1，2，...L-1

其中X为数据的时间序列，g[l]为窗函数，*表示共轭，窗函数的长度为M，步长为S，N为窗的总个数，L为离散频点总个数，最后得到M×N的复数矩阵为其时频矩阵：

矩阵中的每一个元素为一复数z_mn＝a_mn+jb_mn，包含幅值信息和相位信息，其幅值为

相位为

取其幅值信息得到M×N的非负时频幅值矩阵为：

步骤d，利用非负矩阵分解方法对磁场数据的时频幅值矩阵进行分解，为利用非负矩阵分解方法根据数据的时频特性对数据的非负时频幅值矩阵进行自适应分解得到对应的频率分布特征矩阵和频率分布特征矩阵(W矩阵的列向量对应H矩阵行向量)对应的权重矩阵，将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量分离。

非负矩阵分解方法可以将任意一高维非负矩阵分解为两个低维的非负矩阵 W和H，使得W与H的乘积近似的逼近原高维矩阵V，其中W称为基矩阵，H称为权重矩阵。基矩阵W可以看成原矩阵V的一组基底，其列向量表征了数据的特征；而权重矩阵H为矩阵V在基向量矩阵空间的投影，权重矩阵的行向量表示基矩阵W对应列向量随时间变化的幅值。其中原始矩阵V中的元素均为非负数，其分解矩阵W和H中的元素也均为非负数，这使得分解结果更具有物理意义。基于该特性可将磁场数据对应的非负时频矩阵作为原始矩阵进行分解，得到数据的几个时频分布特征及其在时域上的幅值。将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离，不需要去除太阳和地磁活动较强的天的数据，同时可以得到与地震活动更为相关的分量。

非负矩阵分解的数学模型为：

V≈WH (6)

其中V∈R^≥0为一M×N的非负矩阵，W∈R^≥0为一M×R的非负矩阵， H∈R^≥0为一R×N的非负矩阵，分解特征个数R的选择通常应该满足 R＜MN/(M+N)，可以根据数据的已知信息或者经验信息进行设置。

本发明基于K-L散度制定非负矩阵分解的目标函数，同时对权重矩阵H添加行列式约束，使分解结果相对唯一。其中基于广义K-L散度的目标函数为：

i＝1，2，3...M

j＝1，2，3，...N

其中V_ij表示矩阵V的第i行、j列元素，[WH]_ij表示W矩阵和H矩阵相乘得到矩阵的第i行、j列元素。

由于V可以分解出无数个W和H，为使分解结果相对唯一，利用最小行列式准则对权重矩阵H进行约束，通过权重矩阵H的向量张成空间的体积最小进行约束。权重矩阵H的向量张成空间的体积的求解公式为：

则添加最小行列式约束的K-L散度非负矩阵分解目标函数为：

i＝1，2，3，...M

j＝1，2，3，...N

其中det表示求矩阵的行列式，α为平衡参数。

对该目标函数采用乘性迭代方法进行优化求解，通过交替更新变量求得最优解。利用梯度下降法计算更新公式：

其中令μ_w和μ_h分别为：

最终更新迭代公式为：

其中

表示矩阵元素的点乘。

标准化公式为：

由于添加了最小行列式约束不能保证W和H的非负性，故对W和H矩阵中的负值取0，保证其非负性。

上述迭代公式满足迭代次数越多目标函数越收敛。设置迭代停止条件为目标函数小于设置的阈值Th，同时设置迭代次数的上限为N_iter，即若迭代N_iter次后目标函数的值依旧不小于所设置的阈值则强行停止迭代。由于不同轨道得到的时频幅值矩阵的数值各不相同，故应根据不同轨道的数值设置其阈值Th：

其中ε为误差因子。

每条轨道分离得到的W矩阵为：

H矩阵为：

步骤e，根据能量比选出因地震产生的局部影响分量，通过超限率的方法对该分量进行异常轨道判断，是通过计算非负矩阵分解得到的H矩阵的各个行向量中研究区域内能量与整条轨道能量之比，选择能量比最大的分量HH作为因地震产生的局部影响分量，再选用可以反应数据能量平均水平的均方根作为阈值参数，根据HH的值是否大于整条轨道均方根的k倍通过经验选取的进行异常点判断，其中k通过经验选取，认为异常点只出现在研究区域内，且对应标志位显示异常点不是由于已知电离层活动引起的轨道为异常轨道。由于地震事件对数据的影响属于局部影响，故只对轨道中的地震研究区域内的数据产生一定的影响，此时数据能量主要集中在地震影响区域；而太阳活动和地磁活动属于全球影响事件，故会对整条轨道的数据产生影响，此时数据应较为均匀的分布在整条轨道上。所以选取研究区域内能量与整条轨道能量之比最大的分量，该分量数据能量主要集中在地震研究区域内，最有可能反映了地震对局部区域的影响。

其中能量比的计算公式为：

i＝1，2，3，...R

j＝1，2，3，...N

其中P_s表示数据在研究区域内的起点位置，P_e表示数据在研究区域内的终点位置。计算可得到R个能量比结果[η₁，η₂，η₃，...η_R]，选出最大能量比对应的 H行向量HH。

通过超限率法对选择分量HH进行异常轨道判断，选用可以反应数据能量平均水平的均方根作为阈值参数，异常数据应该远远超过该值，且持续一段时间。由于原始数据进行短时傅里叶变换时已经实现了对数据的平滑处理，故可以直接使用HH的值是否大于整条轨道的均方根的k倍进行异常点判断。计算长度为N的时间序列X_i的均方根公式为：

则阈值为kX_rms，大于阈值的点为异常点，同时认为异常点只出现在研究区域内，且对应标志位显示该异常不是由于已知电离层活动引起的轨道为异常轨道。

步骤f，对每日异常轨道个数进行累计，并通过其偏离背景拟合直线的程度检测地震异常，是首先累计每日异常轨道的个数得到累计结果，再使用每日所有轨道个数的最小二乘直线拟合作为背景，通过每日异常轨道累计结果偏离背景拟合直线的程度来显示地震异常的变化过程。因为对所有符合研究条件的轨道进行异常轨道判断之后，得到的结果无法直观的显示其与地震之间的联系，故对每日异常轨道个数进行累计，并使用每日所有轨道个数的最小二乘直线拟合作为背景，利用每日异常轨道个数累计结果偏离背景直线的程度反映地震对数据的影响，偏离程度越大表示其影响越大。

对每日异常轨道个数进行累计，结果为：

N_ano＝[N_ano(1)，N_ano(2)，N_ano(3)，...N_ano(d)] (23)

其中d为所研究的天数。

再对每日所有轨道个数进行累计，结果为：

N_nom＝[N_nom(1)，N_nom(2)，N_nom(3)，...N_nom(d)] (24)

由于每日所有轨道个数与每日异常轨道个数相差较多，通过离差标准化对两者进行归一化处理，其中离差标准化公式为：

计算得到每日异常轨道个数进行累计归一化结果为：

MN_ano＝[MN_ano(1)，MN_ano(2)，MN_ano(3)，...MN_ano(d)] (26)

每日所有轨道个数进行累计归一化结果为：

MN_nom＝[MN_nom(1)，MN_nom(2)，MN_nom(3)，...MN_nom(d)] (27)

对MN_nom进行最小二乘线性拟合，设拟合直线为：

y＝a₀+a₁x (28)

其中参数a₀和a₁估计求解公式为：

其中x为[1，2，3...d]，d为天数，x_i为x中的第i个元素；y为MN_nom，y_i为 MN_nom中的第i个元素；x_iy_i为x中的第i个元素与MN_nom中的第i个元素相乘；经计算得到直线方程参数

和

通过每日异常轨道个数累计曲线偏离背景拟合直线的程度来显示地震变化过程，偏离线性拟合程度计算公式为：

i＝1，2，3，...d

其中Y_i为MN_ano中的第i个元素；

ρ是表示每日异常轨道累计曲线偏离背景拟合直线程度的指数，越大表示数据因受地震影响产生的异常越大。

步骤g，输出偏离程度指数曲线结果图，是将表示每日异常轨道累计结果偏离背景拟合直线程度的指数随时间变化的曲线作为最终地震异常检测结果进行输出。由于ρ是表示每日异常轨道累计曲线偏离背景拟合直线程度的指数，越大表示数据因受地震影响产生的异常越大，故可以表示地震活动的异常变化情况，作为地震异常检测结果进行输出。

针对2016年4月16日发生的7.8级厄瓜多尔地震，以Swarm卫星群Alpha星的Y 分量磁场数据(1Hz)为例。根据Dobrovilsky′s公式R＝10^0.43M(M为地震的震级)选取厄瓜多尔地震的研究区域，以震中为中心±20°经度和纬度的矩形区域为研究区域，研究时间选取震前60天至震后30天，即2016年2月16日至5月16 日。厄瓜多尔地震震中、地震影响区域以及研究区域如图2所示，其中五角星为震中，圆形区域为地震影响区域，正方形区域为研究区域。

a、读取并录入Swarm卫星群Alpha星2016年1月17日至5月16日Y分量磁场数据B_y(1Hz)，并根据数据集中的标志位去除由仪器本身测量有误导致的无效数据。根据卫星数据测量特殊性，即测量值随时间和空间共同变化，以及后续要将因地震产生的局部影响和太阳以及地磁活动导致的全局分量进行分离，故使用地磁纬度在+50°～-50°范围内的且经过地震研究区域的轨道作为研究对象进行数据处理。

以2016年4月8日轨道2为例，Y分量磁场数据的原始曲线如图3所示，可看出磁场原始数据幅值非常大，无法观察到其变化情况，很难直接对其进行分解或异常提取等处理。

b、利用原始Y分量磁场数据减去CHAOS-6磁场模型计算得到的背景场，去除地磁场中幅值较大的背景，并对其结果取一阶差分，去除剩余较小背景值得到磁场数据的变化量。

原始数据减去模型磁场并对其进行差分的公式为：

i＝1，2，3，...L-1

其中B_y为原始Y分量磁场数据，B_m为CHAOS-6磁场模型计算的背景场， t_i为第i个数据对应的时刻，L为每条轨道数据的长度，dB为减模型求一阶差分后的预处理结果。

以2016年4月8日轨道2为例，其原始数据与模型数据曲线，原始数据与模型数据局部放大曲线，原始数据减模型数据的残差数据曲线，一阶差分数据曲线分别如图4中的(1)(2)(3)(4)所示。图4(1)中实线为原始数据曲线，点划线为模型数据曲线，由于原始数据幅值较大且模型数据与其相差很小，故图4(2)将图4(1)其进行部分放大，图4(3)为原始数据减模型数据得到的残差数据，图4(4)为残差数据的一阶差分数据。由结果可看出原始数据幅值非常大，通过减去磁场模型可以去除其幅值较大的背景场，但其残差数据依旧存在一定的背景幅值，需要进一步去除。故对残差数据求一阶差分，去除剩余背景幅值并得到磁场数据的变化量，得到的最终结果背景幅值非常小，且奇异值较为突出。

c、对预处理后的磁场数据进行短时傅里叶变换，得到其时频矩阵，并取其幅值构成非负时频幅值矩阵，作为后续非负矩阵分解的输入矩阵，利用数据的时频幅值特性进行分解。

设预处理后每条轨道数据的时间序列为：

X＝[x₁，x₂，x₃，x₄，...，x_L-1] (2)

其中L为每条轨道数据的长度，根据每条轨道长度不同取值不同。

短时傅里叶变换的离散形式为：

n＝0，1，2，3，...N-1

k＝0，1，2，...L-1

其中x[l]表示X为数据的时间序列，g[l]为窗函数，为减小频谱泄露选用汉宁窗，*表示共轭，窗函数的长度为M，取M＝64，步长为S，取S＝16，N为窗的总个数，L为离散频点总个数。最后得到M×N的复数矩阵为其时频矩阵：

矩阵中的每一个元素为一复数z_mn＝a_mn+jb_mn，包含幅值信息和相位信息，其幅值为

相位为

取其幅值信息得到M×N的非负时频幅值矩阵为：

以2016年4月8日轨道2为例，其短时傅里叶变换的时频幅值矩阵结果图如图5所示，其横坐标为地理纬度，纵坐标为频率，灰度表示其幅值大小。由此可见，该步骤可以将单通道磁场数据转变成M×N的矩阵，同时矩阵均为非负数，可作为非负矩阵分解的输入矩阵，并且具有其物理意义反映了信号的时频幅值特性，后续处理可根据该特性进行分解。

d、利用非负矩阵分解方法对磁场数据的时频幅值矩阵进行分解，是利用非负矩阵分解方法根据数据的时频特性对数据的非负时频幅值矩阵进行自适应分解得到对应的频率分布特征矩阵和其对应的权重矩阵，将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量分离。

非负矩阵分解方法可以将任意一高维非负矩阵分解为两个低维的非负矩阵 W和H，使得W与H的乘积近似的逼近原高维矩阵V，其中W称为基矩阵，H称为权重矩阵。基矩阵W可以看成原矩阵V的一组基底，其列向量表征了数据的特征；而权重矩阵H为矩阵V在基向量矩阵空间的投影，其行向量表示基矩阵W对应列向量随时间变化的幅值。其中原始矩阵V中的元素均为非负数，其分解矩阵W和H中的元素也均为非负数，这使得分解结果更具有物理意义。基于该特性可将磁场数据对应的非负时频矩阵作为原始矩阵进行分解，得到数据的几个时频分布特征及其在时域上的幅值。将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离，不需要去除太阳和地磁活动较强的天的数据，同时可以得到与地震活动更为相关的分量。

非负矩阵分解的数学模型为：

V≈WH (6)

其中V∈R^≥0为一M×N的非负矩阵，W∈R^≥0为一M×R的非负矩阵， H∈R^≥0为一R×N的非负矩阵，分解特征个数R的选择通常应该满足R＜MN/(M+N)，根据经验和已知信息，取R＝3。

基于K-L散度制定非负矩阵分解的目标函数，同时对权重矩阵H添加行列式约束，使分离结果相对唯一。其中基于广义K-L散度的目标函数为：

i＝1，2，3...M

j＝1，2，3，...N

其中V_ij表示矩阵V的第i行、j列元素，[WH]_ij表示W矩阵和H矩阵相乘得到矩阵的第i行、j列元素。

由于V可以分解出无数个W和H，为使分解结果相对唯一，利用最小行列式准则对权重矩阵H进行约束。通过矩阵H的向量张成空间的体积最小进行约束。矩阵H的向量张成空间的体积的求解公式为：

则添加最小行列式约束的K-L散度非负矩阵分解目标函数为：

i＝1，2，3...M

j＝1，2，3，...N

其中det表示求矩阵的行列式，α为平衡参数，取α＝0.005。

对该目标函数采用乘性迭代方法进行优化求解，通过交替更新变量求得最优解。利用梯度下降法计算更新公式：

其中令μ_w和μ_h分别为：

最终更新迭代公式为：

其中表示矩阵元素的点乘。

标准化公式为：

由于添加了最小行列式约束不能保证W和H的非负性，故对W和H矩阵中的负值取0，保证其非负性。

该迭代公式满足迭代次数越多目标函数越收敛。设置迭代停止条件为目标函数小于设置的阈值Th，同时设置迭代次数的上限为N_iter，取N_iter＝500，即若迭代N_iter次后目标函数的值依旧不小于所设置的阈值则强行停止迭代。由于不同轨道得到的时频幅值矩阵的数值各不相同，故应根据不同轨道的数值设置其阈值Th：

其中ε为误差因子，取ε＝0.05。

每条轨道分离得到的W矩阵为：

H矩阵为：

以2016年4月8日轨道2为例，其分解结果W矩阵和H矩阵分别如图6和图7所示，其中W1、W2、W3显示的是该轨道的三个频率分布特征，H1、H2、H3为其对应的幅值，H结果图中虚线之间为研究区域内数据。其中H1的较大幅值只出现在研究区域内，能量也主要分布在研究区域内，该分量应为因地震产生影响的局部分量；H2的幅值波动不大，能量均匀分布在整条轨道中，有可能是背景分量；H3的最大幅值出现在研究区域外，且研究区域内也存在幅值较大的值，该分量可能为太阳或地磁活动产生影响的全局分量。为从时域显示其分离效果，将分解结果重构到时域结果如图8所示。其中图8(1)为一阶差分数据，在研究区域内外均存在明显奇异值；图8(2)为W1和H1对应的时域分量，只在研究区域内存在明显奇异值；图8(3)为W2和H2对应的时域分量，不存在明显奇异值，在整条轨道内均匀分布；图8(4)为W3和H3对应的时域分量，最大奇异值出现在研究区域外，研究区域内也存在幅值较小的奇异值。从分解结果可知利用非负矩阵分解方法对数据的非负时频幅值矩阵进行分解，可以得到其频率分布特征和时间变化幅值，将因地震产生的局部影响分量和由于太阳活动以及地磁活动产生的全局影响分量进行分离，并利用与地震活动更为相关的分量进行地震异常检测。

e、计算H矩阵行向量地震研究区域与整条轨道的能量比，选取能量比最大的分量作为因地震产生的局部影响分量，并利用超限率的方法对该分量进行异常轨道判断。

能量比的计算公式为：

i＝1，2，3

j＝1，2，3，...N

其中P_s表示数据在研究区域内的起点位置，P_e表示数据在研究区域内的终点位置。经计算得到3个能量比结果[η₁，η₂，η₃]，选出最大能量比对应的H行向量HH。

计算行向量HH的均方根值：

取k＝5，设其阈值为kHH_rms，HH中超过阈值的点可认作异常点，得到异常点之后判断其是否只出现在研究区域内的，且对应标志位显示该异常不是由于已知电离层活动导致的异常，若仅存在满足上述条件的异常点，则该轨道为异常轨道。

以2016年4月8日轨道2为例，其分解结果H分别如图7所示，经计算其中H1的能量比为0.924；H2的能量比为0.624；H3的能量比为0.4095。根据研究区域数据占总数据的长度之比计算，数据均匀分布时其能量比应该为0.4，故 H1为因地震影响产生的局部分量，H3应为太阳或者地磁活动影响产生的全局分量，H2为背景分量或其他区分量。经统计所有轨道的最大能量比均值为0.57，较平均值增长45％，且所有异常轨道最大能量比均值为0.80；第二大能量比均值为0.39，与平均值0.4基本相同，最小能量比均值为0.23，较平均值下降42％。由此可见非负矩阵分解分解结果H的分量基本可以反应地震影响产生的局部分量，太阳和地磁活动影响产生的全局分量，在一定程度上将地震局部分量与太阳和地磁活动全局分量分离。

以2016年4月9日轨道3和4月10日轨道1为例，异常轨道判断示意图如图9所示，图中实线为HH分量数据，点划线为其阈值，虚线中间为研究区域内数据。图9(1)为4月9日轨道3，没有点超出阈值，故不被判定为异常轨道；图9(2)为4月10日轨道1，存在超出阈值的点，且仅在研究区域内出现，同时通过检查其标志位没有显示已知电离层活动，可被判定为异常轨道。

f、对每日异常轨道个数进行累计，结果为：

N_ano＝[N_ano(1)，N_ano(2)，N_ano(3)，...N_ano(d)] (23)

其中d为所研究的天数，d＝121。

再对每日所有轨道个数进行累计，结果为：

N_nom＝[N_nom(1)，N_nom(2)，N_nom(3)，...N_nom(d)] (24)

对两条曲线进行离差标准化实现归一化处理，利用离差化公式：

计算得到每日异常轨道个数进行累计归一化结果为：

MN_ano＝[MN_ano(1)，MN_ano(2)，MN_ano(3)，...MN_ano(d)] (26)

每日所有轨道个数进行累计归一化结果为：

MN_nom＝[MN_nom(1)，MN_nom(2)，MN_nom(3)，...MN_nom(d)] (27)

对MN_nom进行最小二乘线性拟合，设拟合直线为：

y＝a₀+a₁x (28)

其中参数a₀和a₁估计求解公式为：

其中x为[1，2，3...121]，y为MN_nom，经计算得到直线方程参数

和

其中x为[1，2，3...d]，d为天数，x_i为x中的第i个元素；y为MN_nom，y_i为 MN_nom中的第i个元素；x_iy_i为x中的第i个元素与MN_nom中的第i个元素相乘；经计算得到直线方程参数

和

通过每日异常轨道个数累计曲线偏离背景拟合直线的程度来显示地震变化过程，偏离线性拟合程度计算公式为：

i＝1，2，3，...121

其中Y_i为MN_ano中的第i个元素；

ρ是表示异常轨道累计曲线偏离背景拟合直线程度的指数，越大表示数据因受地震影响产生的异常越大。

每日异常轨道累计曲线与背景拟合直线归一化结果如图10所示，其中竖直虚线为4月16日7.8级厄瓜多尔地震。由结果可知从震前一段时间开始每日异常轨道个数累计曲线开始不正常的偏离背景拟合直线，并在地震发生附近达到最大偏离，震后几天开始恢复，最后回归平静。每日异常轨道累计曲线偏离背景拟合直线程度指数曲线如图11所示，其中黑色矩形内为地震前检测到的异常，由图1结果可以更直观的显示数据异常情况与地震的关系，每日异常轨道累计结果在震前20天左右开始出现不正常现象，震前10天开始出现大幅增长，并于震后两天达到最大，之后逐渐减小直至震后20天左右恢复正常水平。该结果显示该方法有效的检测到了目标地震产生的异常及异常随地震活动的变化，包括震前异常的开始，地震两天后异常的顶峰，以及之后的异常恢复。

步骤g，将表示每日异常轨道累计结果偏离背景拟合直线程度的指数随时间变化的曲线作为最终地震异常检测结果进行输出。

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