基于大气电离层参数的震前卫星轨道异常识别方法

摘要

本发明公开了一种基于大气电离层参数的地震前兆卫星轨道异常识别方法，包括：原始电离参数的野值剔除及标准化处理；定义卫星轨道上电离参数偏度和峰度；计算卫星各轨道上参数的偏度和峰度值；使用t location‑scale分布拟合得到分布参数及拟合误差；在不同地震震级下，使用不同阈值划分概率密度函数，建立震前卫星异常轨道识别模型；通过同历史地震目录信息对比，以证明本方法的有效和准确性；使用综合评价方法得到面向多参数的轨道识别模型。本发明的优点是：针对大气层电离参数数据量巨大、信息维度高等特点，给出了一种地震前兆卫星轨道异常的识别方法，可用于卫星轨道异常检测和地震前兆识别。

说明书

技术领域

本发明涉及基于大气电离层参数的地震震前卫星轨道异常识别方法，是强震监测大气电离层参数时卫星轨道异常的识别方法，属于航天遥感、地球物理学与计算机科学的交叉领域。

背景技术

我们生活居住的地球每日都有近万次地震发生。作为世界上最严重的自然灾害之一，地震直接导致及间接引发的次生灾害对全球带来了极大的生命和财产损失。我国位于环太平洋地震带和欧亚地震带交汇，是一个地震灾害严重的国家。如何分析、识别地震前兆，提高地震预测精度，对降低地震灾害、减少国家人民损失有着极为重要的意义。

据有关统计，全球每年约发生500多万次地震。近年来全球再次进入地壳板块运动活跃期，屡屡发生高等级强震。仅2015年全球就发生了19次7级以上强震，2016年初至今半年内共发生7级以上强震6次。我国位于亚欧板块、印度洋板块和太平洋板块交界处，欧亚地震带和环太平洋地震带交汇处，地震发生频率高、强度大、分布广、震源浅，是一个地震灾害严重的国家。因此，分析地震前兆中自然环境因素的变化，提高地震预测水平，对降低地震灾害带来的损失具有巨大现实意义。

由于地震突发性强、发生频度高等特点，如何提高地震预报精度和可靠度始终是一个世界性难题。传统方法的空间局限性，无法适用于大面积地域的地震预报。然而随着卫星遥感技术的发展，气体反演、红外遥感等技术被应用于地震的宏观预测，其精度高、重复观察周期短、不受地面条件限制等特点，成为了研究断裂活动性和地震前兆的重要观测手段，极大提高了地震监测能力。

观测表明，在地质学、地球物理学和地球化学等多种地球前兆表现中，电磁异常反应是最为敏感的观测手段之一。在地球100公里至400公里高空，空气极其稀薄，这个高度的空间称为“电离层”。电离层的气体元素都呈游离状态出现，当地磁发生不规则变化时，高层气体元素的核外电子出现丢失或整合现象，空间占用减少，出现了电离层的向下压缩现象。目前研究已证实强震前震区上空的电离层厚度会出现往下压缩的现象，这种压缩会使“太阳风”的高能带电粒子沿着被压缩的电离层向下沉降，从而引起大气电离层参数发生一定程度的突变。从上个世纪70年代至今，已有数百位地震电磁学领域的专家学者研究探索了电离层电磁异常扰动与地震事件间的存在的关联。大量实践研究基础表明部分电磁参数的异常扰动可以作为地震前兆。

对地观测卫星呈周期性围绕地球采集地球表面的参数信息特征，每一条卫星轨道包含了一段完整的信息采集过程。传统的震前异常识别方法大多基于卫星采样点，若采样点附近出现了电离层参数异常(如采样点位于震中附近)，则将该采样点作为异常信息进行处理。但由于地震的偶发性，异常点在全部采样点中所占的比例往往较小，后续震前异常模式识别工作较难进展。而将包含一个或多个异常信息采样点的卫星轨道作为研究单位，即将该轨道标记为震前轨道异常进行研究，在不丢失异常采样点信息的前提下，大大提高了异常数据所占比例，更容易展开震前识别和分析。

发明内容

本发明的目的是基于大气电离层参数，发明一个地震震前卫星轨道异常的识别方法。基于该方法，可有效且准确地识别强震震前记录异常大气电离层参数的卫星轨道，实现了在不丢失异常采样点信息的前提下，扩大异常电离层参数所占比例，提高了强震震前异常识别的能力。

本发明的具体技术方案包括以下几个步骤：

步骤一：数据预处理。针对卫星检测得到的大气电离层参数原始数据，剔除数据测量过程中由于仪器故障等因素产生的野值；卫星通常装载了多个电离层参数的测量仪器，不同参数的量纲往往有所不同，为消除不同量纲属性对数值分析过程造成的影响，利用线性函数对原始数据进行归一化；

步骤二：基于给出的卫星轨道上电离层参数偏度和峰度的新定义，针对预处理后的电离数据计算卫星各轨道上参数的偏度和峰度值；

步骤三：使用t location-scale分布拟合计算得到的各轨道电离层参数的偏度和峰度，得到拟合参数；使用残差平方和SSE(Sum of Squares for Error)作为拟合误差的衡量标准；

步骤四：根据步骤三得到的参数拟合分布，在不同地震震级下划分概率密度函数尾部并记录划分阈值；根据得到的划分阈值序列，建立震前卫星轨道异常识别模型；

步骤五：针对步骤四得到的各电离层参数轨道异常识别模型，结合历史地震目录信息，使用残差平方和SSE和判定系数R²，以证明模型的有效和准确；

步骤六：分别建立不同电离层参数的轨道异常识别模型，根据不同参数模型判定系数R²定义各自权重，使用综合评价方法结合各子模型得到一个综合考虑多参数的识别模型。

本发明的有益效果是：

本发明所提的方法从卫星轨道异常的角度出发，针对监测得到的电离层参数信息进行地震震前的卫星轨道异常变化分析，提出了一种基于大气电离层参数的地震震前卫星轨道异常 识别方法，有效地解决了大气电离异常信息数量少，异常程度低等问题。本发明经多次真实卫星数据测试验证以及相关专家评定，能有效识别强震震前的电离轨道异常。

附图说明

图1是本发明方法的总体流程图。

具体实施方式

下面结合附图和相关算法，对本发明做进一步的说明。

本发明的总体流程如图1所示。

本发明针对大气电离数据，剔除原始参数中野值并做标准化处理；定义了卫星轨道上电离层参数偏度和峰度；计算卫星各轨道上参数的偏度和峰度值；使用t location-scale分布拟合得到分布参数及拟合误差；在不同地震震级下，使用不同阈值划分概率密度函数，建立震前卫星轨道异常识别模型；通过同历史地震目录信息对比，以证明本方法的有效和准确性；使用综合评价方法得到面向多参数的轨道识别模型。具体实施步骤如下：

1.野值剔除及数据标准化

针对本发明中存在的小部分异常数据，使用切比雪夫滤波器进行野值剔除：

$>G_n\left(\omega \right)=\frac{1}{\sqrt{1+{ϵ}^2{T}_n^2\left(\frac{\omega }{{\omega }_0}\right)}}.$>

式中n表示I型切比雪夫滤波器的阶数，ω₀表示滤波器截止频率，ε表示滤波器参数且满足|ε|＜1，T_n(x)＝cos(n·arccos(x))表示n阶切比雪夫多项式。

同样，由于卫星装载的多个测量仪器对不同测量参数的量纲有所不同，为消除不同量纲属性对数值分析过程造成的影响，使用线性函数对原始数据进行归一化。

$>x^{\prime }=\frac{x-x_{\min }}{x_{max}-x_{min}}.$>

式中x为原始数据值，x_min和x_max分别表示属性的最大值和最小值，x′为归一化后的数值，其范围在0到1之间。

其过程可用算法1来描述：

算法1：数据野值剔除及其标准化

输入：原始电离层参数D

输出：预处理后的电离层参数

01： D’＝[]；//数据初始化

02： for i＝l：k //共有k个不同的电离层参数

03： X＝D(：，i)； //取出数据D中的第i个属性

04： X’＝Chev(X)；//经过车比雪夫滤波器，实现野值剔除

05： X”＝(X’-min(X’))/(max(X’)-min(X’))；//数据标准化

06： D’＝[D’；X”]； //将预处理后的第i个电离层参数存入D’中

07： endfor

2.卫星轨道上电离层参数的偏度和峰度

本发明完成各电离层参数的预处理后，进一步定义并计算各卫星轨道上参数的偏度和峰度。针对各卫星轨道上电离检测点，每一个检测点对应的地磁纬度为lat，对应电离层参数为f(lat)，同一轨道上所有检测点地磁纬度的集合为LAT_set，其均值记为下面给出卫星轨道上电离层参数偏度的计算方法：

$>{\alpha }_3=\frac{E{(X-\bar{X})}^3}{{\left(E{\left(X-\bar{X}\right)}^2\right)}^{3/2}}=\frac{\Sigma f\left(lat\right){(lat-\overline{{\mathrm{LAT}}_{set}})}^3}{{(\Sigma f(lat\left){\left(lat-\overline{{\mathrm{LAT}}_{set}}\right)}^2\right)}^{3/2}}.$>

同样给出卫星轨道上电离层参数峰度的计算方法：

$>{\alpha }_4=\frac{E{(X-\bar{X})}^4}{{\left(E{\left(X-\bar{X}\right)}^2\right)}^2}=\frac{\Sigma f\left(lat\right){(lat-\overline{{\mathrm{LAT}}_{set}})}^4}{{(\Sigma f(lat\left){\left(lat-\overline{{\mathrm{LAT}}_{set}}\right)}^2\right)}^2}.$>

其过程可用算法2来描述：

算法2：计算各卫星轨道上电离层参数的偏度和峰度

输入：原始电离层参数D

输出：计算得到的偏度A3和峰度A4

01： A3＝[]；A4＝[]； //数据初始化

02： for i＝l：k//共有k个不同的电离层参数

03： START＝min(D(：，i).orbitnum)；//标记第i个电离层参数的初始轨道序号

04： STOP＝max(D(：，j).orbitnum)；//标记第i个电离层参数的终止轨道序号

05： X＝[]；//初始化

05： for j＝START：STOP //遍历所有轨道序号

06： for w＝1：length(D(：，i))//遍历所有轨道

07： if D(w，i).orbitnum＝＝j//属于轨道j

08： X＝[X；D(w，j)]； //存入同一类

09： endif

10： endfor

11： alpha3＝sum(X.para.*(X lat)^3)/(sum(X.para.*X.lat^2))^1.5；//计算偏度值

12： alpha4＝sum(X.para.*(X.lat)^4)/(sum(X.para.*X.lat^2))^2；//计算峰度值

13： A3＝[A3；alpha3]；A4＝[A4；alpha4]；//记录偏度、峰度

14： endfor

15：endfor

3.t location-scale分布拟合

针对计算得到的卫星轨道上电离层参数的偏度和峰度，本发明使用t location-scale分布对其进行拟合。下面给出t location-scale分布概率密度函数的计算表达公式：

$>f\left(x\right|\mu ,\sigma ,v)=\frac{\Gamma \left(\frac{v+1}{2}\right)}{\sigma \sqrt{\pi v}\Gamma \left(\frac{v}{2}\right)}{\left[\frac{v+{\left(\frac{x-\mu }{\sigma }\right)}^2}{v}\right]}^{-\frac{v+1}{2}}.$>

式中的μ为位置参数，σ为尺度参数，v为形状参数，Γ(□)为伽马函数。

针对给出的偏度和峰度数据，计算不同数值区间下出现的频数和频率，根据以上给出的t location-scale分布概率密度函数，使用非线性最小二乘法进行参数拟合：

minΣ[y-f(x|μ，σ，v)]²。

式中y表示数据x出现的频率值，f(x|μ，σ，v)表示数据x出对应t location-scale分布的概率密度值。保证拟合值与真实值间的平方和最小，由此估计得到概率分布的参数μ，σ和v。

使用残差平方和SSE(Sum of Squares for Error)作为拟合误差的衡量标准，其计算公式如下：

$>SSE=\underset{x\in X}{\Sigma }{(x-x^{\prime })}^2.$>

式中x表示数据的真实值，x′表示对应的拟合值，X表示所有数据x构成的集合。误差平方和SSE越小，表明拟合参数更精确，拟合效果更好。

其具体过程如算法3所述：

算法3：t location-scale分布拟合

输入：计算得到的偏度A3和峰度A4

输出：分布拟合参数MU、SIGMA、NU，拟合残差平方和SSE

01：MU＝[]；SIGMA＝[]；NU＝[]；SSE＝[]；//数据初始化

02：for i＝3：4 //3表示偏度，4表示峰度

03：for j＝l：k //共有k个不同的电离层参数

04： mu＝0；sigma＝0；nu＝0；sse＝0； //数据初始化

05： [mu，sigma，nu，sse]＝dfittool(t location-scale，Ai(：，j))；//使用非线性最小二乘法分布拟合

06： MU＝[MU；mu]；SIGMA＝[SIGMA；sigma]；NU＝[NU；nu]；SSE＝[SSE；sse]；//数据存储

07：endfor

08：endfor

4.建立震前卫星轨道异常识别模型

针对t location-scale分布拟合结果，在不同地震震级下，划分概率密度函数尾部并记录划分阈值。对于一个给定的地震震级M₀，查询历史地震目录，以地震发生时刻t为横坐标，地震发生数量为纵坐标，绘制地震发生时间与数量的曲线图划分阈值从概率密度函数尾部两端开始逐步向分布中心靠拢，记录划分阈值为±x时，尾部所包含的所有轨道序号，以轨道时间t为横坐标，轨道总数量N_±x为纵坐标，绘制轨道异常时间与数量的曲线图N_±x＝g(t)。计算与N_±x间的判定系数R²

$>R^2=\frac{\Sigma {(N_{M_0}-\overline{N_{\pm x}})}^2}{\Sigma {(N_{\pm x}-\overline{N_{\pm x}})}^2}.$>

式中表示N_±x的均值。

记录给定震级M₀下最小R²对应的±x值作为此时的划分阈值。

根据得到的划分阈值，建立震前卫星轨道异常识别模型：

$>N_{M>M_0}={\mathrm{aM}}^{b}t+c.$>

式中M₀表示给定的地震震级，表示截止至时间t共发生震级大于M₀的地震数量，a，b和c为分别为模型的尺度参数、形状参数和位置参数，用于控制不同震级限制下的模型的识别精度。

其具体过程如算法4所述：

算法4：建立震前卫星轨道异常识别模型

输入：计算得到的分布拟合参数MU、SIGMA、NU，历史地震目录表EARTHQUAKE

输出：轨道识别模型参数A、B、C

01：A＝[]；B＝[]；C＝[]；THRESHOLD＝[]；//数据初始化

02：for i＝l：k //共有k个不同的电离层参数

03： for M＝5：0.1：8 //震级选择

04：DATA_M＝select(EARTHQUAKE，M)；//筛选震级大于M的历史震例

05： N_X＝[]；N_M＝[]； //数据初始化

06： for w＝5：0.1：M

07：N_X＝[N_X；count(DATA_M＞w)]； //计算并记录DATA_M表中震级大于w的震例数量

08：N_M＝[N_M；model(MU(i，w)，SIGMA(i，w)，NU(i，w))]；//计算并记录模型中震级大于w的震例数量

09： endfor

10： R2＝sum((N_M-mean(N_X))^2)/sum((N_X-mean(N_X))^2)；//计算判定系数R2

11： if R2＞max_r2

12： max_r2＝R2；temp＝w；//记录当前划分阈值

13： endif

14：THRESHOLD＝[THRESHOLD；w]；//记录划分阈值

15： [a，b，c]＝cftool(THRESHOLD，exp，2)；//二阶指数拟合

16： endfor

17：A＝[A；a]；B＝[B；b]；C＝[C；c]；

18：endfor

5.检验模型的有效和准确性

针对得到的震前卫星轨道异常识别模型，查询历史地震目录信息，分别筛选统计截止至时间t所有发生震级大于M₀地震的数量。以地震发生时刻t为横坐标，地震发生数量N_y为纵坐标，绘制地震发生时间与数量的曲线图N_t＝f(t)。计算与震前卫星轨道异常识别模型间的残差平方和SSE和判定系数R²，其计算公式如下：

$>SSE=\Sigma {(N_{M>M_0}-N_t)}^2,$>

$>R^2=\frac{\Sigma {(N_{M>M_0}-\overline{N_t})}^2}{\Sigma {(N_t-\overline{N_t})}^2}.$>

式中表示N_t的均值。残差平方和SSE和判定系数R²越小，表示得到的模型精度越高，轨道异常识别效果越好，更能反映震前异常状况。

其具体过程如算法5所述：

算法5：建立震前卫星轨道异常识别模型

输入：轨道识别模型参数A、B、C，历史地震目录表EARTHQUAKE

输出：残差平方和SSE、判定系数R2

01：SSE＝[]；R2＝[]；//数据初始化

02：for i＝3：4//3表示偏度，4表示峰度

03：for j＝l：k//共有k个不同的电离层参数

04：N_X＝[]；N_M＝[]；//数据初始化

05：for M＝5：01：8//震级选择

06：N_X＝[N_X；count(EARTHQUAKE＞w)]；//计算并记录EARTHQUAKE表中震级大于w的震例数量

07：N_M＝[N_M；model(A(i，M)，B(i，M)，C(i，M))]；//计算并记录模型中震级大于M的震例数量

08：endfor

09：SSE＝[SSE；sum((N_X-N_M)^2)]；//计算并记录SSE

10：R2＝sum((N_M-mean(N_X))^2)/sum((N_X-mean(N_X))^2)；//计算并记录判定系数R2

11：endfor

12：endfor

6.使用综合评价方法得到面向多参数的轨道识别模型

针对不同的电离层参数，已经得到多个含有不同参数值得轨道识别模型。使用一种综合评价的方法结合以上多个子模型，以获得一个面向多参数的轨道识别模型。

现有k个电离层参数分别建立轨道异常识别模型各模型分别计算得到一个判定系数R_i²。使用综合评价方法设定各子模型的权重

$>w_i=\frac{R_i^2}{\underset{i=1}{\overset{k}{\Sigma }}{R}_i^2},$>

得到最后的综合模型为

$>N_{M>M_0}^C=\underset{i=1}{\overset{k}{\Sigma }}{w}_i{N}_{M>M_0}^i.$>

其具体过程如算法6所述：

算法6：建立面向多参数的轨道识别模型

输入：轨道识别模型参数A、B、C，判定系数R2

输出：最终模型参数PARA

01：PARA＝[]；W＝[]；N_C＝[]； //数据初始化

02：for i＝l：k //共有k个不同的电离层参数

03：W＝[W；R2(i)/sum(R2)]；//计算并记录子模型权重

04：endfor

05：for M＝5：0.1：8 //震级选择

06：temp＝0；//数据初始化

07：for i＝l；k //共有k个不同的电离层参数

08：temp＝temp+W(i)*model(A(i，M)，B(i，M)，C(i，M)，M)； //子模型加权求和

09：endfor

10：N_C＝[N_C；temp]；

11：endfor

12：PARA＝cftool(N_C，exp，2)；//二阶指数拟合。