基于三测站数据特征点的电离层扰动传播测量方法及系统与流程

本发明涉及电离层波动传播分析技术领域，特别涉及一种基于三测站空间关系以及监测数据特征点的电离层扰动传播测量技术方案。

背景技术：

地球高层大气由于受到太严辐射作用，原子、分子发生电离形成大量自由电子，从而在60～2000km高度范围内形成电离层。电离层扰动是指电离层结构偏离其常规形态的急剧变化，又称电离层骚扰。电离源的突变、非平衡态动力学过程、不稳定的磁流动力过程和某些人为因素等，都可引起电离层扰动。电离层扰动会对卫星通讯、GNSS定位等产生严重影响。而寻找电离层扰动传播源，研究其传播方向、速度、范围等，对于探索磁层、电离层、对流层之间相互作用，空间天气变化规律具有重要意义。目前国际上采用的电离层扰动传播速度提取方法主要有距离-时间图法、三测站干涉测量法等。其中距离-时间图法主要是在图上按照斜率计算，准确性较差，计算结果偏差较大；三测站干涉测量法需要较长的稳定扰动波形，且计算过程十分复杂，易受外界因素影响，在扰动波形较短时计算结果不准确。因此，本领域亟待更实用的技术方案出现。

技术实现要素：

为了解决现有方法中存在的问题，本发明提供一种实现简单、结果精确的电离层扰动传播测量技术方案。

本发明技术方案提供一种基于三测站数据特征点的电离层扰动传播测量方法，包括以下步骤，

步骤1，建立局部坐标系，包括根据已知的三个测站穿刺点的坐标，建立以其中一点A为原点的局部坐标系，经线方向为y轴，纬线方向为x轴，计算其余两个测站穿刺点B、C在该局部坐标系中的坐标(x_B,y_B)、(x_C,y_C)；

步骤2，选取监测数据特征点，得到异常到达A、B、C三个点的时间分别为t_A、t_B、t_C，得到异常扰动从A到B和从A到C的传播时间，分别记为时间差dt_AB＝t_B-t_A，dt_AC＝t_C-t_A；

步骤3，求异常传播路径的斜率k和B、C在扰动传播方向投影点D、E坐标(x_D,y_D)、(x_E,y_E)如下，

步骤4：求异常传播速度和方位角，得到测量结果如下，

异常传播速度大小为

异常传播方位角为

其中，α_AD、α_AE为D、E坐标(x_D,y_D)、(x_E,y_E)相应的异常传播方位角。

而且，步骤2中，选取监测数据特征点的实现方式为，从三个监测站的二阶差分时间序列取极小值或极大值，若三个曲线在该区域前后波形非常相似，选取该时域内极小值或极大值作为异常到达三个测站的时间。

本发明相应提供一种基于三测站数据特征点的电离层扰动传播测量系统，包括以下模块，

初始化模块，用于建立局部坐标系，包括根据已知的三个测站穿刺点的坐标，建立以其中一点A为原点的局部坐标系，经线方向为y轴，纬线方向为x轴，计算其余两个测站穿刺点B、C在该局部坐标系中的坐标(x_B,y_B)、(x_C,y_C)；

时差提取模块，用于选取监测数据特征点，得到异常到达A、B、C三个点的时间分别为t_A、t_B、t_C，得到异常扰动从A到B和从A到C的传播时间，分别记为时间差dt_AB＝t_B-t_A，dt_AC＝t_C-t_A；

投影模块，用于求异常传播路径的斜率k和B、C在扰动传播方向投影点D、E坐标(x_D,y_D)、(x_E,y_E)如下，

测量输出模块，用于求异常传播速度和方位角，得到测量结果如下，

异常传播速度大小为

异常传播方位角为

其中，α_AD、α_AE为D、E坐标(x_D,y_D)、(x_E,y_E)相应的异常传播方位角。

而且，时差提取模块中，选取监测数据特征点的实现方式为，从三个监测站的二阶差分时间序列取极小值或极大值，若三个曲线在该区域前后波形非常相似，选取该时域内极小值或极大值作为异常到达三个测站的时间。

本发明利用三个监测站的扰动出现的时间差以及三个监测站的空间几何关系提取电离层扰动在空间中传播速度的大小和方向。该技术方案利用测站穿刺点精确的几何关系，在电离层扰动时间较短的情况下，仅用少量特征点即可准确解算出电离层波动传播的速度大小和方向，并且在实验中得到了验证。本发明对于研究日地空间环境、电离层、磁层耦合、电离层行扰传播形态等具有重要科学意义和应用价值，具有重要的市场前景。

附图说明

图1是本发明实施例的技术方案原理图；

图2是本发明实施例应用示例的三个监测站以及震中分布图；

图3是本发明实施例应用示例的三个GPS站在电离层扰动期间的监测数据二阶差分时间序列图；

图4是本发明实施例应用示例的按照三站特征点法计算得到的电离层扰动传播的速度和方向结果图。

具体实施方式

下面通过实施例，对本发明的技术方案作进一步具体的说明。

本发明首先根据已知的三个测站穿刺点的坐标建立以其中一点为原点的局部坐标系，然后根据三个测站监测数据特征点出现的时间计算电离层波动到达三个测站之间的时间差；最后根据三个测站穿刺点几何关系以及波动到达三个测站穿刺点的时间差在局部坐标系中解算出电离层波动传播的速度大小和方向。

参见图1，本发明实施例提供的技术方案包括如下步骤：

步骤1：建立局部坐标系

从GPS卫星到测站的卫星信号与电离层的交点称为穿刺点，三个测站在同一时刻会形成三个穿刺点。为了简化后续解算过程，本发明选取三个穿刺点中的任意一点作为原点，经线方向为y轴，纬线方向为x轴，建立局部坐标系。同时解算出其他两个点在局部坐标系中的坐标。图1中A、B、C三点为监测点，以A点为原点建立局部坐标系，B、C坐标为(x_B,y_B)、(x_C,y_C)。

步骤2：计算扰动的传播时间

假设电离层异常传播方向如图1中所示，异常平行向前推进。从B点向异常传播方向做垂线，交点为D点；从C点向异常传播方向做垂线，交点为E点。则异常到达B点的时间和到达D点的时间相同；同样，异常达到C点的时间和到达E点的时间相同。

具体实施时，可选取监测数据特征点确定时间，优选的实现方式为，从三个监测站的二阶差分时间序列取极小值或极大值，若三个曲线在该区域前后波形非常相似，选取该时域内极小值或极大值作为异常到达三个测站的时间。

A、B、C三个测站异常特征点出现时间已经确定，设异常到达A、B、C三个点的时间分别为t_A、t_B、t_C，则异常从A到B和从A到C的传播时间分别为

步骤3：求异常传播路径的斜率k和除原点以外其余两点在扰动传播方向投影点坐标

设异常到达D、E三个点的时间分别为t_D、t_E，异常从A到D和从A到E的传播时间分别为dt_AD、dt_AE，传播距离为S_AD、S_AE。因为t_B＝t_D，t_C＝t_E，则dt_AD＝dt_AB，dt_AE＝dt_AC。设异常速度大小为V，方位角为α。则

由于时间差dt_AD和dt_AE都有方向，即正负，能够表示异常传播的方向。而A点为原点，D、E在X轴的坐标x_D、x_E和在Y轴的坐标y_D、y_E也有方向，可得到

上式中，dt_AD＝dt_AB，dt_AE＝dt_AC均为已知，设异常传播方向直线的斜率为k，则过C点和过B点的垂线斜率为-1/k。

其余两点B、C电离层异常传播方向投影点分别为D点和E点，D点为直线AD和直线BD的交点，由AD和AB方程

可得D点坐标为

同理可得E点坐标为

公式(3)、(5)、(6)联立

最终求得，以A点为原点求出异常传播的斜率为

并且

根据(9)求得异常传播斜率后，即可将k回代到式(5)、(6)得到D、E坐标(x_D,y_D)、(x_E,y_E)。

步骤4：求异常传播速度和方位角

由于求得k后，D点和E点坐标都可以根据式(5)、(6)求出来，则异常传播速度大小为

异常传播方位角(从Y轴正方向起，依顺时针方向到异常传播方向线之间的水平夹角)α_AD、α_AE可由D点坐标和E点坐标求出，最终以A点为原点的求出的异常传播方位角为

同理，可以将B、C设为原点，计算异常传播的速度和方位角，实现方式相同。

采用以上方案进行实验，采用中国大陆构造环境监测网络西藏地区三个GPS测站的数据。研究对象为2015年尼泊尔地震引起的同震电离层扰动。

步骤1：建立局部坐标系

选取图2中三个GPS测站xzar、xzzf、xzrk中任意一个测站的穿刺点为原点，记为点A，其他两个点记为A点和B点；经线方向为y轴，纬线方向为x轴，建立局部坐标系，同时计算其余两个GPS站的穿刺点在该局部坐标系中的坐标。

步骤2：选取监测数据特征点，计算扰动的传播时间。

从图3中三个监测站的二阶差分时间序列可以看出，三个箭头位置均出现了极小值，且三个曲线在该区域前后波形非常相似，因此选取该时域内最小值作为异常到达三个测站的时间，得到异常到达A、B、C三个点的时间分别为t_A、t_B、t_C，则计算时间差dt_AB＝t_B-t_A，dt_AC＝t_C-t_A。

步骤3：求异常传播路径的斜率k和B、C在扰动传播方向投影点D、E坐标

步骤4：求异常传播速度和方位角

求得k后，D点和E点坐标都在步骤3已求出来，则异常传播速度大小为

异常传播方位角可由AD或AE求出，最终

最终计算结果如图4所示，V～＝930m/s，α(Azi)～＝56°，图中曲线为该时间段内三个GPS测站穿刺点的轨迹，箭头代表电离层扰动的方向。

具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方式实现相应系统。

本发明实施例相应提供一种基于三测站数据特征点的电离层扰动传播测量系统，包括以下模块，

初始化模块，用于建立局部坐标系，包括根据已知的三个测站穿刺点的坐标，建立以其中一点A为原点的局部坐标系，经线方向为y轴，纬线方向为x轴，计算其余两个测站穿刺点B、C在该局部坐标系中的坐标(x_B,y_B)、(x_C,y_C)；

时差提取模块，用于选取监测数据特征点，得到异常到达A、B、C三个点的时间分别为t_A、t_B、t_C，得到异常扰动从A到B和从A到C的传播时间，分别记为时间差dt_AB＝t_B-t_A，dt_AC＝t_C-t_A；

投影模块，用于求异常传播路径的斜率k和B、C在扰动传播方向投影点D、E坐标(x_D,y_D)、(x_E,y_E)；

测量输出模块，用于求异常传播速度和方位角。

各模块具体实现可参见相应步骤，本发明不予赘述。

需要强调的是，本发明所述的实施例是说明性的，而不是限定性的。因此本发明包括并不限于具体实施方式中所述的实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。