一种基于非震动态背景场的地震TEC异常信息提取方法

摘要

本发明提出一种基于非震动态背景场的地震TEC异常信息提取方法，包括以下步骤：动态背景场建立、小波多尺度分析、支持向量机回归模型和异常信息提取方法；本发明通过非震动态背景场，在考虑了TEC自身的周期性变化，同时还考虑了太阳活动、地磁活动等空间环境的干扰，可以提取的地震TEC异常更为准确，同时可以减少地震电离层TEC异常信息提取过程中其它非震因素带来的干扰和误差，提高地震TEC异常检测的准确率，在地震期有较高的预警率，而在非震期虚警率较低，可以提高地震TEC异常信息提取效率。

说明书

技术领域

本发明涉及地震TEC信息提取领域，尤其涉及一种基于非震动态背景场的地震TEC异常信息提取方法。

背景技术

电离层是日地空间环境的一个重要组成部分，它的剧烈变动会对人类的生产生活产生巨大的影响。

现有的地震电离层TEC异常提取主要通过以下两类方法：一是滑动时窗法；二是滑动四分位法，但是以上两种方法存在较多的缺陷，比如：由于背景场构建过程中没有考虑TEC自身的周期性波动影响，导致提取的异常存在周期性信息；由于背景场中没有考虑太阳活动等非震因素引起的地球空间环境的扰动的影响，使得在没有地震发生的情况下也能大量提取出TEC异常，存在较高的虚警率；同时，在地震事件和其它非震扰动事件同时发生时，无法判别TEC异常来源。因此，为了解决太阳活动等引起地球空间环境扰动以及TEC自身周期性波动的非震动态TEC背景场，以减少地震电离层TEC异常信息提取过程中其它非震因素带来的干扰和误差，提高地震TEC异常检测的准确率，降低虚警率，本发明提出一种基于非震动态背景场的地震TEC异常信息提取方法，以解决现有技术中的不足之处。

发明内容

针对上述问题，本发明通过非震动态背景场，在考虑了TEC自身的周期性变化，同时还考虑了太阳活动、地磁活动等空间环境的干扰，可以提取的地震TEC异常更为准确，同时可以减少地震电离层TEC异常信息提取过程中其它非震因素带来的干扰和误差，提高地震TEC异常检测的准确率，在地震期有较高的预警率，而在非震期虚警率较低，可以提高地震TEC异常信息提取效率。

本发明提出一种基于非震动态背景场的地震TEC异常信息提取方法，包括以下步骤：

步骤一：动态背景场建立

在地震电离层TEC异常信息提取前，首先应找到地震造成的电离层TEC扰动与其他扰动因素，然后对归一化后的太阳在10.7cm辐射通量的值和TEC值进行傅里叶变换，分析两者频谱图，对太阳活动与TEC之间的关系进行分析，然后先剔除TEC数据的半年周期和年周期性成分，再利用SVR模型进行与非震因素的回归建模；

对于某一地点电子总量为一个整体，其受到其他扰动因素影响，在时间序列上表现不断波动，根据公式(1)进行计算可以得出主要的非震TEC背景场T_back；

T₀＝T₁+T₂+T₃+T₄+T_e(1)

公式(1)中，T₁为太阳周期性活动的影响，主要表现在长周期变化上；T₂为TEC自身稳定的年周期及半年周期成分波动；T₃为地磁和空间宇宙环境活动的影响；T₄为极端天气、海啸等其他非地震因素的TEC扰动；T_e为地震引起的TEC扰动；

其中，太阳周期性活动的影响T₁、TEC自身稳定的年周期及半年周期成分波动T₂和地磁和空间宇宙环境活动的影响T₃为TEC波动的主要来源，准确模拟出T₁、T₂和T₃的扰动情况，再根据计算公式(2)进行计算可以得出最后的非震TEC背景场T_back；

T_back＝T₁+T₂+T₃(2)

步骤二：小波多尺度分析

利用小波分析中的一维小波多尺度分析，分析公式如公式(3)所示：

公式(3)中，a，b均为常数，函数是由Ψ(t)函数伸缩和平移变换得到，a为伸缩因子，b为平移因子，因此a,大小组合决定了Ψ_a,b(t)的形式；

小波多尺度分解可以看做一组高、低通滤波器对信号进行层叠滤波的结果，其中高通滤波器输出各尺度下的高频细节成分，低通滤波器输出反映信号的低频近似成分，通过对TEC原始信号T₀多层分解得到不同周期的高频细节成分，TEC原始信号中7阶高频细节成分正好对应着半年周期成分，TEC原始信号中8阶高频细节成分正好对应着年周期成分，提取7阶高频细节成分与8阶高频细节成分，即为公式(1)中TEC自身周期成分T₂，而TEC原始信号T₀减去TEC自身周期成分T₂得到TEC残差T_rest；

步骤三：支持向量机回归模型

利用支持向量机回归模型，将SVM推广到非线性回归估计中，得到支持向量回归算法，是通过核函数将线性不可分样本数据T＝{(x₁,y₁),……(x_l,y_l)}∈(x×y)^l映射到线性高维特征空间，然后通过解决最优化问题得到回归函数f(x)＝ω·φ(x)+b，φ(x)为非线性映射函数；

采用ε不敏感损失函数和核函数的支持向量回归数学形式如公式(4)所示：

公式(4)中，C为惩罚参数；ξ_i和为模型的松弛变量；

当实际值和预测值的误差未超过预设的ε时，则认为在该时间点的预测是不存在损失的，再利用SVR构建太阳活动、地磁活动指数与T_rest之间的回归模型；

步骤四：异常信息提取方法

T1：建立数据样本集

将2004年-2011年的全球电子总量地图数据按世界时和地理空间位置分成不同的数据集，然后分别对每个数据集进行处理；

T2：小波多尺度分解

将T1中的每个数据集按时间序列排成一维向量，然后分别对每个一维向量进行小波多尺度分解，其中小波基为db8，小波分解得到每个数据集序列的7阶高频和8阶高频，然后用原始的一维向量序列T₀减去小波分解的7阶高频与8阶高频结果，即得到去除半年周期与年周期的TEC残差T_rest；

T3：SVR建模过程

模型建模训练过程中，将表征太阳活动和地磁活动的9个变量输入模型内；

混频数据匹配：SVR建模前需要将混频数据进行匹配，统一所有数据为2小时分辨率；

非震样本筛选：为了保证模型训练输入输出均为未受地震干扰的数据，SVR建模前需要对每个数据集进行非震样本筛选，选择不受地震干扰的数据输入模型进行训练；

输入输出标准化处理：在SVR回归建模之前，需要对回归建模数据指标进行标准化处理，选择最大最小0-1归一化方法进行标准化处理；

SVR建模：对归一化之后的指标数据，随机选取输入数据中80％的数据作为SVR模型的训练集，剩下20％作为SVR模型的测试集，SVR模型核函数选择RBF核函数，模型输出为模型参数、模型平均误差σ以及模型精度，模型平均误差σ为TEC值T_r′与实际观测的TEC值T_rest之间的差值；

T4：非震动态背景场构建

对于每一个时刻、每一地理位置，将其对应的表征太阳活动和地磁活动的9个变量输入到已经训练好的SVR模型，输出该时间点、该位置空间环境TEC分量的预测值T_r′，然后T_r′与该时间点、该位置小波分解得到小波7阶高频值H7和8阶高频值H8进行相加，即得到该时刻该空间位置的TEC非震动态背景值T_back；

然后分别计算所有时刻、所有地理位置的非震动态背景值，即得到整个时空的TEC非震动态背景场；

T5：基于动态背景场的地震TEC异常提取

原始TEC观测值T₀减去TEC非震动态背景值T_back即得到TEC残差T_rest，然后取SVR模型1.5倍的平均误差σ为阈值进行异常判定；

若T_rest大于1.5倍平均误差σ即认为该TEC残差为地震引起的正异常扰动；

当T_rest小于-1.5倍平均误差σ即认为该TEC残差为地震引起的负异常扰动；

当T_rest小于1.5倍平均误差σ且大于-1.5倍平均误差σ，即认为该TEC残差为正常偏差，将正常偏差范围内的TEC残差值赋值为0，最后输出TEC残差，即为地震TEC异常提取结果。

进一步改进在于：所述步骤一中的其他扰动因素包括太阳活动、地磁活动、空间环境以及极端天气。

进一步改进在于：所述步骤一的公式(1)中先利用小波多尺度分解提取T₂部分信息，然后根据太阳活动与地磁活动和TEC信号相关性，利用支持向量机回归算法构建太阳、地磁活动指数与T₁、T₃两部分信息之间的回归预测模型，最后整合三部分信息，再根据计算公式(2)得到最后的非震TEC背景场T_back。

T_back＝T₁+T₂+T₃>

而剩余成分TEC残差T_rest中，地震TEC扰动信号得到放大，通过分析TEC残差，这样可以更为准确的提取地震电离层TEC扰动，剩余成分TEC残差T_rest计算公式如公式(5)所示。

T_rest＝T₄+T_e>

进一步改进在于：所述步骤步骤四T1中使用的全球电子总量地图为欧洲定轨中心发布的全球电子总量地图，所述欧洲定轨中心发布的全球电子总量地图的时间分辨率为2小时，即每2个小时一张地图，每张地图经度覆盖范围从西经180°至东经180°，空间分辨率为5°；纬度覆盖范围从北纬87.5°至南纬87.5°，空间分辨率为2.5°，每张地图共包括5183个格网点。

进一步改进在于：所述步骤四T1中的数据集为同一位置、同一时刻的数据放到一个样本集形成的数据集，所述地理空间位置为格网点，所述地理空间位置数量与网格点数量一致，共5183个，所述世界时包括12个时刻。

进一步改进在于：所述步骤四T1中的数据集一共包括12×5183＝62196个数据集，每个数据集从2004年-2011年共有2920个数据。

进一步改进在于：所述步骤四T3中的混频数据匹配具体为：对于原时间分辨率高于2小时的输入指标数据，采用求和取平均处理，降低时间分辨率到2小时，如DST指数时间分辨率为1小时，则将其0时刻值与1时刻值求和得到第一个时刻的值，以此将DST的一天24个数据值降低为一天12个数据值；对于原时间分辨率低于2小时的输入指标数据，按照该数据实际代表的时间插值，如F10.7指数原时间分辨率为1天，则对于该数据集1天12个数据对应的F10.7指数均取为同一个数据。

进一步改进在于：所述步骤四T3中的非震样本筛选原则为：将时间上前后15天且空间上R＝10^0.43M范围内出现地震的数据剔除，其中R为孕震区范围，M为地震的震级。

进一步改进在于：所述步骤四T3中最大最小0-1归一化方法计算公式如公式(6)所示：

公式(6)中，X_min与X_max分别为该指标输入模型的最大值与最小值，X_i为归一化之前的值，为归一化之后的值。

进一步改进在于：所述步骤四T3具体过程为：建模过程中，先将训练集数据输入SVR回归模型，其中输入为表征太阳活动和地磁活动的9个变量，模型输出为T_rest，通过不断调整模型参数进行模型训练，找到能刻画输入、输出数据回归关系的最优模型，并保存最优SVR模型的参数，按照每个时刻、每个地理空间位置分成进行训练与建模，总模型数量为12×5183＝62196个，然后将测试集数据中表征太阳活动和地磁活动的9个变量输入最优SVR模型得预测的TEC值T_r′，T_r′即表示太阳活动等非震因素对TEC值影响的部分，然后计算出模型回归输出T_r′与测试集原始T_rest之间平均偏差和平均精度，得到模型的精度与模型平均误差σ。

本发明的有益效果为：通过本发明方法提出的非震动态背景场，在考虑了TEC自身的周期性变化，同时还考虑了太阳活动、地磁活动等空间环境的干扰，可以提取的地震TEC异常更为准确，同时可以减少地震电离层TEC异常信息提取过程中其它非震因素带来的干扰和误差，提高地震TEC异常检测的准确率，降低了虚警率，并且可消除太阳、地磁非震电离层干扰和TEC自身周期波动带来的系统误差，同时还能提高地震TEC异常探测能力，通过利用SVR模型建立9个能够表征太阳活动与地磁活动的空间环境指标与电离层TEC值之间的非线性映射关系，从而确立TEC非震期动态背景场，本发明方法在地震期有较高的预警率，而在非震期虚警率较低，可以提高地震TEC异常信息提取效率。

附图说明

图1为本发明方法中地震TEC异常信息提取方法实施流程示意图。

图2为本发明方法中TEC小波多尺度分解示意图。

图3为本发明方法中小波分解流程图。

图4为本发明方法中SVR建模流程示意图。

图5为本发明方法中TEC非震动态背景场构建流程示意图。

具体实施方式

为了使发明实现的技术手段、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

根据图1、2、3、4、5所示，本实施例提出了一种基于非震动态背景场的地震TEC异常信息提取方法，包括以下步骤：

步骤一：动态背景场建立

在地震电离层TEC异常信息提取前，首先应找到地震造成的电离层TEC扰动与其他扰动因素，其他扰动因素包括太阳活动、地磁活动、空间环境以及极端天气，然后对归一化后的太阳在10.7cm辐射通量的值和TEC值进行傅里叶变换，分析两者频谱图，对太阳活动与TEC之间的关系进行分析，然后先剔除TEC数据的半年周期和年周期性成分，再利用SVR模型进行与非震因素的回归建模；

对于某一地点电子总量为一个整体，其受到其他扰动因素影响，在时间序列上表现不断波动，根据公式(1)进行计算可以得出主要的非震TEC背景场T_back；

T₀＝T₁+T₂+T₃+T₄+T_e>

公式(1)中，T₁为太阳周期性活动的影响，主要表现在长周期变化上；T₂为TEC自身稳定的年周期及半年周期成分波动；T₃为地磁和空间宇宙环境活动的影响；T₄为极端天气、海啸等其他非地震因素的TEC扰动；T_e为地震引起的TEC扰动；

其中，太阳周期性活动的影响T₁、TEC自身稳定的年周期及半年周期成分波动T₂和地磁和空间宇宙环境活动的影响T₃为TEC波动的主要来源，准确模拟出T₁、T₂和T₃的扰动情况，再根据计算公式(2)进行计算可以得出最后的非震TEC背景场T_back；

T_back＝T₁+T₂+T₃>

先利用小波多尺度分解提取T₂部分信息，然后根据太阳活动与地磁活动和TEC信号相关性，利用支持向量机回归算法构建太阳、地磁活动指数与T₁、T₃两部分信息之间的回归预测模型，最后整合三部分信息，再根据计算公式(2)得到最后的非震TEC背景场T_back；

T_back＝T₁+T₂+T₃>

而剩余成分TEC残差T_rest中，地震TEC扰动信号得到放大，通过分析TEC残差，这样可以更为准确的提取地震电离层TEC扰动，剩余成分TEC残差T_rest计算公式如公式(5)所示：

T_rest＝T₄+T_e>

步骤二：小波多尺度分析

利用小波分析中的一维小波多尺度分析，分析公式如公式(3)所示：

公式(3)中，a，b均为常数，函数是由Ψ(t)函数伸缩和平移变换得到，a为伸缩因子，b为平移因子，因此a,大小组合决定了Ψ_a,b(t)的形式；

小波多尺度分解可以看做一组高、低通滤波器对信号进行层叠滤波的结果，其中高通滤波器输出各尺度下的高频细节成分，低通滤波器输出反映信号的低频近似成分，通过对TEC原始信号T₀多层分解得到不同周期的高频细节成分，TEC原始信号中7阶高频细节成分正好对应着半年周期成分，TEC原始信号中8阶高频细节成分正好对应着年周期成分，提取7阶高频细节成分与8阶高频细节成分，即为公式(1)中TEC自身周期成分T₂，而TEC原始信号T₀减去TEC自身周期成分T₂得到TEC残差T_rest；

步骤三：支持向量机回归模型

利用支持向量机回归模型，将SVM推广到非线性回归估计中，得到支持向量回归算法，是通过核函数将线性不可分样本数据T＝{(x₁,y₁),……(x_l,y_l)}∈(x×y)^l映射到线性高维特征空间，然后通过解决最优化问题得到回归函数f(x)＝ω·φ(x)+b，φ(x)为非线性映射函数；

采用ε不敏感损失函数和核函数的支持向量回归数学形式如公式(4)所示：

公式(4)中，C为惩罚参数；ξ_i和为模型的松弛变量；

当实际值和预测值的误差未超过预设的ε时，则认为在该时间点的预测是不存在损失的，再利用SVR构建太阳活动、地磁活动指数与T_rest之间的回归模型；

步骤四：异常信息提取方法

T1：建立数据样本集

使用全球电子总量地图为欧洲定轨中心发布的全球电子总量地图作为数据，欧洲定轨中心发布的全球电子总量地图的时间分辨率为2小时，即每2个小时一张地图，每张地图经度覆盖范围从西经180°至东经180°，空间分辨率为5°；纬度覆盖范围从北纬87.5°至南纬87.5°，空间分辨率为2.5°，每张地图共包括5183个格网点，再将将2004年-2011年的欧洲定轨中心发布的全球电子总量地图数据按世界时和地理空间位置分成不同的数据集，所述世界时共12个时刻，所述地理空间位置为格网点，所述地理空间位置数量与网格点数量一致，共5183个，数据集一共包括12×5183＝62196个数据集，每个数据集从2004年-2011年共有2920个数据，然后分别对每个数据集进行处理；

T2：小波多尺度分解

将T1中的每个数据集按时间序列排成一维向量，然后分别对每个一维向量进行小波多尺度分解，其中小波基为db8，小波分解得到每个数据集序列的7阶高频和8阶高频，然后用原始的一维向量序列T₀减去小波分解的7阶高频与8阶高频结果，即得到去除半年周期与年周期的TEC残差T_rest；

T3：SVR建模过程

模型建模训练过程中，将表征太阳活动和地磁活动的9个变量输入模型内；

混频数据匹配：对于原时间分辨率高于2小时的输入指标数据，采用求和取平均处理，降低时间分辨率到2小时，如DST指数时间分辨率为1小时，则将其0时刻值与1时刻值求和得到第一个时刻的值，以此将DST的一天24个数据值降低为一天12个数据值；对于原时间分辨率低于2小时的输入指标数据，按照该数据实际代表的时间插值，如F10.7指数原时间分辨率为1天，则对于该数据集1天12个数据对应的F10.7指数均取为同一个数据；

非震样本筛选：为了保证模型训练输入输出均为未受地震干扰的数据，SVR建模前需要对每个数据集进行非震样本筛选，非震样本筛选原则为：将时间上前后15天且空间上R＝10^0.43M范围内出现地震的数据剔除，其中R为孕震区范围，M为地震的震级，选择不受地震干扰的数据输入模型进行训练；

输入输出标准化处理：在SVR回归建模之前，需要对回归建模数据指标进行标准化处理，选择最大最小0-1归一化方法进行标准化处理，最大最小0-1归一化方法计算公式如公式(6)所示：

公式(6)中，X_min与X_max分别为该指标输入模型的最大值与最小值，X_i为归一化之前的值，为归一化之后的值；

SVR建模：对归一化之后的指标数据，随机选取输入数据中80％的数据作为SVR模型的训练集，剩下20％作为SVR模型的测试集，SVR模型核函数选择RBF核函数，模型输出为模型参数、模型平均误差σ以及模型精度，模型平均误差σ为TEC值T_r′与实际观测的TEC值T_rest之间的差值，建模过程中，先将训练集数据输入SVR回归模型，其中输入为表征太阳活动和地磁活动的9个变量，模型输出为T_rest，通过不断调整模型参数进行模型训练，找到能刻画输入、输出数据回归关系的最优模型，并保存最优SVR模型的参数，按照每个时刻、每个地理空间位置分成进行训练与建模，总模型数量为12×5183＝62196个，然后将测试集数据中表征太阳活动和地磁活动的9个变量输入最优SVR模型得预测的TEC值T_r′，T_r′即表示太阳活动等非震因素对TEC值影响的部分，然后计算出模型回归输出T_r′与测试集原始T_rest之间平均偏差和平均精度，得到模型的精度与模型平均误差σ；

T4：非震动态背景场构建

对于每一个时刻、每一地理位置，将其对应的表征太阳活动和地磁活动的9个变量输入到已经训练好的SVR模型，输出该时间点、该位置空间环境TEC分量的预测值T_r′，然后T_r′与该时间点、该位置小波分解得到小波7阶高频值H7和8阶高频值H8进行相加，即得到该时刻该空间位置的TEC非震动态背景值T_back；

然后分别计算所有时刻、所有地理位置的非震动态背景值，即得到整个时空的TEC非震动态背景场；

T5：基于动态背景场的地震TEC异常提取

原始TEC观测值T₀减去TEC非震动态背景值T_back即得到TEC残差T_rest，然后取SVR模型1.5倍的平均误差σ为阈值进行异常判定；

若T_rest大于1.5倍平均误差σ即认为该TEC残差为地震引起的正异常扰动；

当T_rest小于-1.5倍平均误差σ即认为该TEC残差为地震引起的负异常扰动；

当T_rest小于1.5倍平均误差σ且大于-1.5倍平均误差σ，即认为该TEC残差为正常偏差，将正常偏差范围内的TEC残差值赋值为0，最后输出TEC残差，即为地震TEC异常提取结果。

通过本发明方法提出的非震动态背景场，在考虑了TEC自身的周期性变化，同时还考虑了太阳活动、地磁活动等空间环境的干扰，可以提取的地震TEC异常更为准确，同时可以减少地震电离层TEC异常信息提取过程中其它非震因素带来的干扰和误差，提高地震TEC异常检测的准确率，降低了虚警率，并且可消除太阳、地磁非震电离层干扰和TEC自身周期波动带来的系统误差，同时还能提高地震TEC异常探测能力，通过利用SVR模型建立9个能够表征太阳活动与地磁活动的空间环境指标与电离层TEC值之间的非线性映射关系，从而确立TEC非震期动态背景场，本发明方法在地震期有较高的预警率，而在非震期虚警率较低，可以提高地震TEC异常信息提取效率。

以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。