基于变分模态分解的卫星磁场数据时变背景场建立方法

摘要

本发明为一种基于变分模态分解的卫星磁场数据时变背景场建立方法。读取Swarm A星的磁场Y分量数据，减去CHAOS‑6磁场模型数据去除主磁场的影响，剔除无效轨道；利用变分模态分解将磁场Y分量分解成2个模态，将不同的特征分离到不同的模态中去；通过各个模态与相应地磁指数的对比，找出与地磁活动相关性高的模态并将其去除；选择研究时间范围、研究区域，以一个小时为时间间隔，将该区域内的轨道数据划分到24个当地时间范围；将该区域网格化，剔除受地震影响的轨道数据；计算磁场Y分量能量的均值作为背景值并拟合网格内24个当地时间的背景值。利用变分模态分解将信号的特征进行分离，去除地磁活动干扰，同时对信号滤除高斯噪声，建立精确的时变背景场。

说明书

技术领域

本发明属于电离层磁场数据时变背景场技术领域，具体地而言为一种基于变分模态分解的卫星磁场数据时变背景场建立方法。

背景技术

地震灾害是地球上经常发生的一种自然灾害，是不可避免的。近年来，越来越多的研究表明，电磁异常对短临地震十分敏感。相较地面观测，电磁卫星技术具有全天候、全球性等优势。1958年，前苏联发射了第一颗地磁测量卫星SPUTNIK－3，2004年，法国发射了DEMETER卫星，其目的旨在检测和研究与地震活动相关的电离层电磁扰动。国内外学者研究了DEMETER卫星观测到的各种数据，比如电场、磁场、离子密度、离子温度、电子密度、电子温度、高能粒子通量等，结果表明，在地震前后均有电磁异常现象。2013年，欧洲航天局发射了Swarm卫星，其目的旨在精确测量地球磁场。

目前对电磁背景场的构建和电磁异常提取展开了深入研究，但已有的研究仅考虑电磁背景场的年变化，昼夜变化，季节变化等，没有考虑电磁背景场的当地时间变化，而电磁背景场的当地时间变化也有较大的差异。因此，基于Swarm卫星磁场数据的时变背景场建立是可行的。

CN110794464A公开了一种基于非震动态背景场的地震TEC异常信息提取方法，考虑TEC自身的周期性变化、太阳活动等空间环境的干扰，利用支持向量机回归模型(SVR)构建太阳活动、地磁活动指数与TEC原始信号减去TEC自身周期成分的残差之间的回归模型，从而建立非震动态TEC背景场，以减少地震电离层TEC异常信息提取过程中其它非震因素带来的干扰和误差，提高地震TEC异常检测的准确率，降低虚警率。

CN109409230A公开了一种空地数据互增强的震前热异常自动分级方法，基于遥感数据获取兴趣区域待提取的亮温图像，通过前置滑窗选择与待提取亮温图像方差最小的待选背景场亮温图像作为背景场图像；对处理后的待提取亮温图像和背景场图像进行地表温度反演，基于权重配比得到均衡地表温度进而获取热红外异常提取图像。

袁媛等利用2013-2016年间历年同期5月和同季度(4-6月)的中分辨率成像光谱仪(MODIS)卫星热红外海表温度SST数据的均值构建背景场，从而得到研究区域在2013-2016年间历年同月和同季度的SST背景场分布状况，并基于建立的背景场提取中国台湾新北海域MS6.2地震的热红外异常。杨牧萍在剔除受地磁活动影响的数据后，基于网格化对DEMETER卫星的电场VLF功率谱数据求取网格中的中值作为背景值，并分析东北亚地区背景场的年变化、昼夜变化、季节变化。

DEMETER卫星成功发射以后，学者们做了较多的电离层背景场研究。但到目前为止，尚未有基于变分模态分解的卫星磁场数据时变背景场建立方法，即使是对卫星磁场数据建立背景场，也只对静磁时间段的数据进行分析，没有将地磁活动的影响去除后建立背景场。且已有的研究仅考虑电磁背景场的年变化，昼夜变化，季节变化等，没有考虑电磁背景场随当地时间的变化。

发明内容

本发明所要解决的技术问题在于提供一种基于变分模态分解的卫星磁场数据时变背景场建立方法，利用变分模态分解去除地磁活动的影响，且变分模态分解的降噪功能可以对分解出的模态去除高斯噪声，同时弥补了以往只对静磁时间段的数据进行分析的不足，从而建立更加可靠的时变背景场。

本发明是这样实现的，一种基于变分模态分解的卫星磁场数据时变背景场建立方法，该方法包括：

a、读取SwarmA星的磁场Y分量数据，减去CHAOS-6磁场模型数据，剔除磁场Y分量有跳变的轨道并根据标志位剔除无效轨道，判断每条轨道相邻采样点对应的ap指数是否有跳变，当有跳变时，在跳变处将该轨道分割开；

b、利用VMD分解将经过步骤a处理的磁场Y分量分解成2个模态，此时不同的特征被分离到不同的模态中，同时分解过程中通过维纳滤波对每一个模态滤除高斯噪声；

c、计算步骤b分解的各个模态与地磁指数的相关系数，找出与地磁指数相关性高的模态并去除；

d、对剩余的模态选择研究时间范围以及研究区域，以一个小时为时间间隔，将研究区域内的轨道数据划分到24个当地时间范围；

e、将选定的区域网格化，剔除每个网格中受地震影响的卫星轨道数据；

f、计算每个网格中24个当地时间范围内磁场Y分量能量的背景值并将同一网格内24个当地时间的背景值进行拟合；

g、绘制网格化后区域内24个当地时间的背景图和单一网格内背景值拟合图。

进一步地，步骤a所述的读取Swarm A星磁场Y分量数据包括读取SwarmA星的矢量磁力仪测得的磁场的Y分量，即NEC中指向东向的分量磁场数据。

进一步地，步骤b所述的利用VMD分解将磁场Y分量分解成2个模态，将不同的特征分离到不同的模态中去，包括：对每条轨道减模型后的磁场残余值Y(t)利用VMD分解成2个固有模态函数IMF分量，表达式为：

其中：u

变分模态分解定义一种新的固有模态函数，表达式为：

其中，

对信号分解的固有模态函数u

u

估计解析信号的中心频率并将所有模态函数的频谱调制到相应的基频带；根据频域卷积定理将解析信号调制至基频带，其表达式如下：

对上式L

式中，{u

引入拉格朗日乘子λ和惩罚因子α求约束变分模型的最优解：

<λ(t)，f(t)-∑

利用交替方向乘子算法迭代搜索约束变分问题的鞍点,具体求解过程是交替更新

其中，

在使用交替方向乘子算法搜索最优解时，其停止条件为：

其中，ε是迭代收敛值。

进一步地，步骤c所述的计算各个模态与地磁指数的相关系数，找到与地磁活动指数相关性大的模态并去除包括：

步骤b变分模态分解后，计算每条轨道各个模态的平均能量作为该轨道每个模态的能量值；计算各个模态的能量值的公式为：

其中，i＝1,2代表第i个固有模态；j＝1,…N代表第j个数据点；N代表该轨道的总数据点数。

计算所有轨道分解出的2个模态的能量值与对应的地磁指数ap在时滞为0时的相关系数：

其中

与ap指数相关系数高的模态主要反应地磁活动情况，相关系数低的模态表明与地磁活动关系不大，通过去除与地磁指数相关性高的模态，以去除地磁活动的影响。

进一步地，步骤d所述的选择研究时间范围、研究区域，是选取研究区域为地理纬度20°×地理经度20°，研究时间为三年的SwarmA星的磁场Y分量，将该区域内的轨道采样点划分到24个当地时间范围，以一个小时为时间间隔，分别是LT0(00：00-00：59)，LT1(01：00-01：59),.....,LT22(22：00-22：59),LT23(23：00-23：59)，当地时间LT＝世界时间UTC+α/15，α是地理经度。

进一步地，步骤e区域网格化，剔除每个网格中受地震影响的卫星轨道包括：将区域划分成地理经度5°×地理纬度5°的网格，对每个网格找出研究时间范围内五级以上的地震，并去除震前7天到震后7天共15天的卫星轨道，以去除地震的影响。

进一步地，步骤f所述的计算每个网格中24个当地时间范围内磁场Y分量能量的背景值并将同一网格内24个当地时间LT的背景值进行拟合包括：

先对第二个固有模态的磁场Y分量平方求能量，公式为：

去除奇异值，再以能量的均值作为该网格在该当地时间LT的背景值，选择区域内的任一网格，将该网格内24个当地时间LT的背景值进行随当地时间变化的背景值拟合。

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

在本发明中，通过对SwarmA星磁场Y分量进行变分模态分解，能够将地磁活动与其他影响因素分离到不同的模态中，通过去除与地磁指数相关性高的模态以去除地磁活动干扰，避免后续建立的时变背景场受到地磁活动的影响。把数据根据当地时间划分到24个LT(当地时间)时间范围，利用均值能够有效计算24个LT(当地时间)时间范围内磁场的背景场，从而建立时变背景场。

本发明充分利用了Swarm A星静磁与非静磁时间段的磁场Y分量数据，将变分模态分解应用到磁场数据进行特征分离，去除地磁活动的干扰，弥补了以往背景场研究中只对静磁时间段的数据进行分析的不足，同时利用变分模态分解过程中的维纳滤波对每个模态去除高斯噪声实现降噪。且本发明通过背景场随当地时间的变化，建立了时变背景场，区别于以往只研究背景场的季节变化，年变化等。同时可以将该方法推广到全球时变背景场的建立。从而建立全球卫星磁场数据时变背景场，使后续地震前兆异常的提取不受背景场变化的影响，提取结果更精确。

附图说明

图1为基于变分模态分解的卫星磁场数据时变背景场建立方法的流程图；

图2(a)和2(b)分别为Swarm A星磁场Y分量原始曲线和减模型后Swarm A星磁场Y分量减模型后残余值曲线；

图3为Swarm A星磁场Y分量残余值变分模态分解结果曲线，图3(a)和3(b)分别为变分模态分解后的第一个固有模态和第二个固有模态；

图4为2015年-2017年间所有轨道分解后的两个固有模态的平均能量与对应的ap指数的曲线图；

图5为选择的研究区域以及将其网格化后示意图；

图6为研究区域内每个网格的24个LT(当地时间)的背景值示意图。

图7为研究区域内某一网格24个LT(当地时间)的背景值的高斯拟合曲线示意图。

具体实施方式

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

参见图1所示，基于变分模态分解的卫星磁场数据时变背景场建立方法，首先读取SwarmA星的磁场Y分量数据，减去CHAOS-6磁场模型以去除主磁场的影响，再剔除磁场Y分量有跳变的轨道并根据数据的标志位Flags_B剔除有误数据；判断每条轨道相邻时间点对应的ap指数是否有跳变，当有跳变时，在跳变处将该轨道分割开；利用VMD分解将磁场Y分量分解成2个模态，将不同的特征分离到不同的模态中去，同时分解过程中通过维纳滤波对每一个模态滤除高斯噪声；通过各个模态的平均能量与地磁指数的对比，找到与地磁活动相关性高的成分并将其去除，只对剩余的模态进行分析；选择研究时间范围、研究区域，利用当地时间(LT)将该区域内的轨道的每个采样点划分到24个当地时间(LT)范围；将该区域网格化，找出研究时间范围内每个网格中五级以上的地震，针对每个地震去除震前7天到震后7天共15天的卫星轨道以去除地震对背景场的影响；对网格中每个当地时间范围内的数据点的磁场Y分量平方求能量，将这些能量值排序，去除首尾各20％的数据点，对剩余的磁场Y分量能量值求均值作为该网格在该当地时间(LT)的背景值；将同一网格内24个当地时间(LT)的背景值进行高斯拟合。具体包括：

a、读取Swarm A星的磁场Y分量数据，减去CHAOS-6磁场模型数据，剔除磁场Y分量有跳变的轨道并根据标志位剔除无效轨道，判断每条轨道相邻采样点对应的ap指数是否有跳变，当有跳变时，在跳变处将该轨道分割开；

b、利用VMD分解将磁场Y分量分解成2个模态，将不同的特征分离到不同的模态中去；

c、计算各个模态与地磁指数的相关系数，找出与地磁活动指数相关性高的模态并去除；

d、选择研究时间范围、研究区域，以一个小时为时间间隔，将该区域内的轨道采样点划分到24个当地时间(LT)范围；

e、将该区域网格化，剔除每个网格中受地震影响的卫星轨道；

f、计算每个网格中24个当地时间(LT)范围内磁场Y分量能量的背景值并将同一网格内24个当地时间(LT)的背景值进行拟合；

g、绘制网格化后区域内24个当地时间(LT)的背景图和单一网格内背景值拟合图；

进一步地，步骤a所述的读取SwarmA星磁场Y分量数据，是指读取Swarm A星的矢量磁力仪测得的磁场的Y分量，即NEC(North-East-Center)中指向东向的分量磁场数据，已经有相关研究表述，Y分量磁场能够更清晰的记录到岩石圈磁场异常变化；减去CHAOS-6磁场模型数据是为了去除主磁场的影响，卫星测量的磁场为地磁场，主要包括主磁场、岩石圈磁场、变化磁场和感应磁场四部分。地磁场幅值为10

进一步地，步骤b所述的利用VMD分解将磁场Y分量分解成2个模态，将不同的特征分离到不同的模态中去，主要包括：

对每条轨道减模型后的磁场残余值Y(t)利用VMD分解成2个固有模态函数(IMF)分量，用公式表示为：

其中：u

变分模态分解是以经典维纳滤波、希尔伯特变换和混频的变分问题求解为基础，通过迭代搜寻变分模型最优解自适应的将信号分解成多个有限带宽的固有模态函数。维纳滤波是一种基于最小均方误差准则的最优估计器。变分模态分解的降噪功能是通过应用维纳滤波在分解的过程中对每一个模态进行去除高斯噪声来实现的。因此，在用变分模态分解对结果进行分解时，同时可以对分解出的模态进行高斯滤波，滤除高斯噪声。

变分模态分解定义一种新的固有模态函数，表达式为：

其中，

对信号分解的固有模态函数u

u

估计解析信号的中心频率并将所有模态函数的频谱调制到相应的基频带；根据频域卷积定理将解析信号调制至基频带，其表达式如下：

对上式L

式中，{u

引入拉格朗日乘子λ和惩罚因子α求约束变分模型的最优解：

<λ(t)，f(t)-∑

利用交替方向乘子算法迭代搜索约束变分问题的鞍点具体求解过程是交替更新

其中，

在使用交替方向乘子算法搜索最优解时，其停止条件为：

其中，ε是迭代收敛值。

进一步地，步骤c所述的计算各个模态与地磁指数的相关系数，找到地磁活动指数相关性大的模态并去除主要包括：

VMD分解后，计算每条轨道各个模态的平均能量作为该轨道每个模态的能量值。计算各个模态的能量的公式为：

其中，i＝1,2代表第i个固有模态；N代表该轨道的数据点数。

地磁活动是地磁场的扰动变化程度，主要通过地磁指数反映。地磁指数通过地磁台站测量，是用于描述地磁扰动强度的一种分级指标。国际上常用的地磁指数主要有以下几种：(1)Kp指数。国际上选定12个标准地磁台站，将这些台站描述3小时地磁扰动强度的K指数扩展成28级得到Ks指数，Ks指数各个时段的均值即为Kp指数。(2)ap指数。每一级Kp指数可以转化成一定幅度的ap指数，称为行星性等效三小时幅度。(3)Dst指数。计算地磁赤道附近5个经度间隔大致相同的地磁台站每小时水平强度变化的平均值，即可得到Dst指数。(4)AE指数。用于描述极区磁亚暴强度，取决于极地附近地磁台站每小时的水平强度变化。磁暴是一种使所有地磁要素发生剧烈变化的全球性地磁扰动现象。

计算所有轨道分解出的2个模态的能量值与对应的ap指数在时滞为0时的相关系数：

其中

与地磁指数ap相关系数高的模态主要反应地磁活动情况，相关系数低的模态表明与地磁活动相关性不大，通过去除与地磁指数相关性高的模态，以去除地磁活动的影响。

进一步地，步骤d所述的选择研究时间范围、研究区域，是选取研究区域为20°×20°(地理纬度×地理经度)，研究时间为三年的Swarm A星的磁场Y分量。将该区域内的轨道采样点划分到24个当地时间(LT)范围，以一个小时为时间间隔，分别是LT0(00：00-00：59),LT1(01：00-01：59),.....,LT22(22：00-22：59),LT23(23：00-23：59)。LT(当地时间)＝UTC(世界时间)+α/15(α是地理经度)。

进一步地，步骤e所述的将该区域网格化，剔除每个网格中受地震影响的卫星轨道主要包括，将该区域划分成5°×5°(地理经度×地理纬度)的网格，以确保每个网格中有足够多的数据。对每个网格找出研究时间范围内五级以上的地震，并去除震前7天到震后7天共15天的卫星轨道，以去除地震对背景场的影响。

进一步地，步骤f所述的计算每个网格中24个当地时间(LT)范围内磁场Y分量能量的背景值并将同一网格内24个当地时间(LT)的背景值进行拟合主要包括：

先对第二个固有模态的磁场Y分量平方求能量，公式为：

去除奇异值，再对剩余的磁场Y分量能量值求均值，以能量的均值作为该网格在该当地时间(LT)的背景值。选择区域内的任一网格，将该网格内24个当地时间(LT)的背景值进行随当地时间变化的背景值拟合。

以下结合实例对本发明进行进一步详细说明。

以2015年-2017年间Swarm卫星群中的A星的矢量磁场Y分量(1Hz)为例。

a、读取Swarm卫星A星磁场Y分量数据，是读取Swarm A星的矢量磁力仪测得的磁场的Y分量，即NEC(North-East-Center)中指向东向的分量磁场数据；减去CHAOS-6磁场模型数据是为了去除主磁场的影响，得到Swarm A星的Y分量磁场残余值；剔除磁场Y分量有跳变的轨道，是判断一条轨道中相邻时刻测得的磁场Y分量是否有跳变。根据标志位剔除无效数据选择有效数据，是指卫星矢量磁场数据的质量主要通过标志位Flags_B进行筛选，剔除无效数据后，仅对有效数据进行分析。以2015年1月2日轨道15为例，该轨道的磁场Y分量原始曲线如图2(a)所示，减模型后的磁场Y分量残余值如图2(b)所示；可以看出，减CHAOS-6模型后，在数据中显示出低频趋势与高频抖动，数据在低频趋势上叠加了高频微弱的抖动。

b、先判断每条轨道相邻采样点对应的ap指数是否有跳变，当有跳变时，在跳变处将该轨道分割开，这是为了让每条轨道所有时刻对应的ap指数是同一个值，不会出现跳变，从而保证每条轨道的平均能量与该轨道的地磁活动指数是一一对应的关系。利用VMD分解将磁场Y分量分解成2个模态，将不同的特征分离到不同的模态中去，是指对每条轨道减模型后的的磁场残余值Y(t)利用VMD分解成2个固有模态函数(IMF)分量，用公式表示为：

其中：u

以2015年1月2日轨道15为例，该轨道VMD分解后的结果如图3所示，其中图3(a)是VMD分解后的第一个固有模态IMF1，图3(b)是VMD分解后的第二个固有模态IMF2。从图中可以看到第一个固有模态的幅值大，频率低，与减模型后的残余值的低频变化趋势相似，而第二个固有模态的幅值小，频率较高，与叠加在原始数据低频趋势上的高频抖动相似；通过变分模态分解，将A星的磁场Y分量按照频率分离到两个固有模态中。

c、计算所有轨道分解出的两个模态与地磁指数的相关系数，找出与地磁活动指数相关性高的模态并去除。VMD分解后，计算每条轨道各个模态的平均能量作为该轨道各个模态的能量值。计算各个模态的平均能量的公式为：

其中，i＝1,2代表第i个固有模态；N代表该轨道的数据点数。对每条轨道分解出的2个固有模态都计算其平均能量，在步骤b中，已经保证了每条轨道所有数据的ap指数都相同。研究所有轨道的各个模态的平均能量与对应的地磁指数ap指数曲线如图4所示，选用窗长为401个数据点，步长为1的滑动窗对IMF1、IMF2的平均能量和ap指数进行平滑处理。图4(a)、(b)中的黑色的线分别代表每条轨道分解的IMF1、IMF2的平均能量，图4(c)中的黑色的线代表对应的ap指数；图4(a)、(b)中的白色的线分别代表每条轨道分解的IMF1、IMF2的平均能量的滑动平均值，图4(c)中的白色的线代表ap指数滑动平均值。从图可以看到，第一个模态的平均能量的滑动平均变化趋势与地磁指数相近，这意味着，当地磁指数高、地磁活跃时，相应的第一个模态的平均能量也大，猜测第一个模态可能与地磁活动相关；而第二个模态的平均能量的滑动平均变化趋势则与地磁活动的变化趋势无明显相似。为了定量分析固有模态与地磁活动的关系，进行以下步骤。

计算所有轨道分解出的2个模态的平均能量值与对应的地磁指数ap在时滞为0时的相关系数：

其中

采用窗长为81个数据点，步长为1个数据点的滑动窗，将窗沿着数据进行滑动，分别计算每个窗内2个固有模态与相应的地磁指数的相关系数。2015-2017年间，所有有效轨道的第一个固有模态与相应的地磁指数的相关系数的中值达到85.90％，相关程度非常高，说明第一个固有模态主要反应地磁活动情况。而第二个固有模态与地磁指数的相关系数的中值为39.87％，相关程度不太高，因此第二个固有模态受到地磁活动的影响不大。经过变分模态分解，将受地磁活动大的信息分离到第一个固有模态中，因此接下来只对第二个固有模态进行分析从而建立背景场。

e、选择研究时间范围、研究区域，是选取研究区域5°N-25°N(地理纬度),35°E-55°E(地理经度)内2015年到2017年Swarm A星的磁场Y分量。利用当地时间(LT)将该区域内的轨道的每个采样点划分到24个当地时间(LT)范围，以一个小时为时间间隔，分别是LT0(00：00-00：59),LT1(01：00-01：59),.....,LT22(22：00-22：59),LT23(23：00-23：59)。LT(当地时间)＝UTC(世界时间)+α/15(α是地理经度)

f、将该区域网格化，划分成5°×5°(地理经度×地理纬度)的网格，以确保每个网格中有足够多的数据。如图5所示，其中的矩形是选择的研究区域以及将其网格化后的效果图。对每个网格找出2015年-2017年间五级以上的地震，对每个地震去除震前7天到震后7天共15天的卫星轨道，以去除地震对背景场的影响。

g、计算每个网格中24个当地时间(LT)范围内磁场Y分量能量的背景值是指对网格中每个当地时间范围内的数据点的磁场Y分量平方求能量，对磁场Y分量平方求能量的公式为：

将能量值从大到小排序，去除前后各20％的点。再对剩余的磁场Y分量能量值求均值，以均值作为该网格在该当地时间(LT)的背景值。绘制网格化后区域内24个当地时间(LT)的背景图，如图6所示，由图可知，当地时间为白天(LT6-LT17)的背景值远大于当地时间为夜晚(LT0-LT5，LT18-LT23)，且在当地时间为LT11(11:00-12:00)时，背景值最大，大致符合随当地时间先增大后减小的规律。将同一网格内24个当地时间(LT)的背景值进行拟合，图7为该研究区域内某一网格的24个LT(当地时间)的背景值的高斯拟合曲线示意图，拟合的确定系数为0.9506，拟合标准差为0.1235。从图中可以看出该区域24个LT(当地时间)的背景值大致符合高斯分布。

本发明基于变分模态分解的卫星磁场数据时变背景场建立方法，利用CHAOS-6模型去除地磁主磁场的影响；对Swarm A星的磁场Y分量数据进行变分模态分解，通过各个模态与地磁指数的对比，找到与地磁活动相关性高的模态并将其去除，只对剩余的模态进行分析；选择研究区域，研究时间，划分区域内轨道数据到24个LT(当地时间)时间段并通过均值定义每个时间段内的背景值，建立时变背景场。该方法弥补了以往只对静磁时间段的数据进行分析，不能充分利用Swarm卫星的A星非静磁时间段数据的不足，全面建立时变背景场。

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