一种可识别电离层行波扰动和赤道等离子泡的探测指数

1.本发明属于gnss数据处理与电离层扰动探测技术领域，具体涉及一种可识别电离层行波扰动和赤道等离子泡的探测指数。


背景技术：

2.电离层作为日地空间环境的重要组成部分，是一个典型的非平稳、多尺度、非线性复杂系统。由于受到太阳辐射、地磁活动以及其它因素的影响，电离层电子密度经常会发生急剧性的不规则变化，即电离层扰动。电离层扰动现象能显著的影响通过该电离层区域的无线电传播路径，严重时甚至能导致卫星信号失锁、航天器失常甚至受损等情况。为了防止电离层扰动对卫星导航定位、无线电通信的严重影响，监测电离层扰动现象，对提高导航定位精度、保障航天活动安全等具有非常重要的意义。
3.在中低纬度地区，电离层行波扰动(tids)和赤道等离子泡(epbs)是两种典型的电离层扰动现象。tids是由于太阳耀斑、地磁风暴、日食、地震、飓风、火山爆发等各种事件引起电离层电子密度发生的准周期性波动。epbs是夜间电离层等离子体密度的消耗，与由电离层f层底部扰动源触发的广义rayleigh-taylor不稳定性有关，产生的等离子泡可以到达低纬甚至中纬度地区，并可以从几米蔓延到数百公里。tids和epbs都会使得通过该区域的无线电信号的振幅和相位发生快速变化，从而对全球卫星导航系统(gnss)安全性和功能产生影响，在某些极端情况下可能会导致卫星通讯中断。并且这两类事件发生在中低纬度地区，可能出现同一时空产生、共存、耦合甚至互相转化等情况。因此，快速探测并准确识别tids和epbs至关重要。
4.目前探测tids和epbs等电离层扰动事件的手段包括气辉成像仪、测高仪、高频雷达等。然而这类设备分布稀疏，观测区域固定并且价格昂贵，难以满足大规模监测预警的需求。gnss由于其覆盖全球、全天候等特点为大范围和高精度监测电离层扰动提供了解决方案。但基于gnss数据的电离层扰动探测方法目前还存在着一些新的问题与挑战。目前最常用的基于gnss的电离层空间环境扰动探测方法主要分为二维去趋势的总电子含量图(dtec map)和单站的指数方法两类。前者可以较为准确地表征出tids特征，但在缺乏足够密集的参考站数据情况下，难以用于tids传播检测。而后者目前最常见的roti指数方法能探测出epbs，但其对tids却不敏感，无法有效识别tids。另一方面，对于电离层扰动类型(如tids、epbs)的区分通常是对图像进行人工目视判断，目前还难以通过量化指标自动区别并识别这两种扰动类型。


技术实现要素：

5.针对现有技术不足，本发明提供一种可识别电离层行波扰动和赤道等离子泡的探测指数iroti，能有效解决基于gnss探测电离层扰动目前遇到的诸多困难，为本领域技术人员提供探测并区分tids和epbs的方法。
6.为了达到上述目的，本发明提供的技术方案是一种可识别电离层行波扰动和赤道
等离子泡的探测指数，包括以下步骤：
7.步骤1，获取观测区域的gnss原始观测文件和广播星历文件；
8.步骤2，对gnss原始观测数据进行周跳探测，获得具有连续观测弧段的历元编号与卫星prn号；
9.步骤3，基于gnss原始观测数据构建无几何相位观测方程，获得包含伪距和相位硬件延迟偏差的电离层总电子含量stec；
10.步骤4，利用步骤2得到的连续观测弧段的历元编号和卫星prn号对步骤3得到的stec进行筛选，得到处于同一连续弧段的stec序列；
11.步骤5，利用处于同一连续弧段的stec序列计算电离层扰动指数rot指数，进而获得roti指数的时间序列；
12.步骤6，采用滑动平均法对stec进行去趋势处理，得到去趋势的总电子量dtec的时间序列；
13.步骤7，联合步骤5获得的roti指数和步骤6获得的dtec，计算iroti指数。
14.而且，所述步骤2中利用turbedit方法对步骤1获得的原始观测数据进行周跳探测，得到具有连续观测弧段的历元编号与卫星prn号，并将最短连续弧段设置为μ1历元以确保后续滑动平均获取去趋势的总电子量dtec能够正常进行。
15.而且，所述步骤3中gnss信号由伪距和载波相位观测值两部分组成，可以表示为：
[0016][0017]
式中，下标r和上标s分别表示gnss接收机和卫星的序列号，下标i表示频率，和为卫星与接收机之间沿视线(los)方向的伪距和相位观测值，表示卫星与接收机之间的几何距离，c是真空中的光速，dtr和dts分别代表gnss接收机和卫星的钟差，为沿los方向的对流层延迟，为频率i处的电离层延迟，b
r,i
和分别表示gnss接收机和卫星在频率i处的伪距硬件延迟，b
r,i
和分别表示gnss接收机和卫星在频率i处的相位小数偏差，λi为波长，为相位整周模糊度，和是观测噪声和其他各种误差的组合项。
[0018]
无几何相位观测方程通过不同频率的载波相位观测值差值得到，即：
[0019][0020]
式中，为无几何相位组合观测值，和分别为f1、f2频率下沿los方向的相位观测值，为沿给定los方向的电离层延迟，γi＝(f0/fi)2是比例因子，其中f0＝10.23mhz，λ1和λ2分别为f1、f2频率的波长，和分别为f1、f2频率的整周模糊度，表示多径效应、噪声误差和其他类型的误差。
[0021]
电离层总电子含量(stec)定义为从卫星到接收机沿los方向的电子密度的积分，频率为f的stec与电离层延迟i在长度单位上的关系可以描述为：
[0022]
[0023]
式中，stec值以tecu为单位表示，1tecu＝10
16
el/m2；为无几何相位组合观测值；γi＝(f0/fi)2是比例因子，其中f0＝10.23mhz；表示多径效应、噪声误差和其他类型的误差。stec的截止高度角设置为μ2，以削弱多路径效应的影响。
[0024]
而且，所述步骤5中rot指数为两个连续时间期间stec值的时间导数，即：
[0025][0026]
式中，rot值作为测量相位波动活动的参数，以tecu/min为单位表示；分别为卫星k在i-1、i时刻的电离层总电子含量，t
i-1
、ti为历元时间，k为卫星序列号。
[0027]
roti为特定时期内rot的标准差，具体计算方式如下：
[0028][0029]
式中，尖括号表示特定时期内的平均值。
[0030]
而且，所述步骤6中滑动平均法计算去趋势的总电子量(dtec)的公式如下：
[0031][0032]
式中，stec(i)为第i时刻的stec值，n为滑动窗口，stec(j)为第j时刻的stec值，j的取值范围为[i-n/2,i+n/2]。
[0033]
而且，所述步骤7中通过计算dtec值在特定时期内的积分作为判断扰动发生的指数idtec：
[0034][0035]
式中，|dtec(i)|为i时刻dtec的绝对值，δt为离散dtec的时间间隔，n为特定时期内dtec个数，idtec值以tecu
·
min为单位表示。
[0036]
结合idtec和roti提出了iroti指数，该指数考虑振幅和相位的变化来探测电离层扰动进而区分tids和epbs，具体计算公式可以表示为：
[0037]
iroti＝idtec
·
roti(8)
[0038]
式中，iroti值以tecu2为单位表示。
[0039]
当iroti＜μ3tecu2时，可以判断电离层扰动类型为tids，当iroti＞μ4tecu2时，可以认为扰动类型为epbs。
[0040]
与现有技术相比，本发明具有如下优点：
[0041]
在传统电离层扰动指数roti的基础上，结合去趋势的总电子量dtec的优势提出iroti指数，可以有效解决当前传统电离层扰动探测指数tids和epbs难以同时探测的问题，并通过阈值将二者识别区分。
附图说明
[0042]
图1为本发明iroti指数计算总体流程框架。
[0043]
图2为本发明实施例人工模拟tids和epbs信号的dtec时间序列图，其中，图2(a)为
tids的模拟信号，图2(b)为epbs的模拟信号。
[0044]
图3为本发明实施例2015年6月22日hkws站计算的roti、iroti和dtec值，其中，图3(a)为hkws站计算的roti，图3(b)为hkws站计算的iroti，图3(c)为hkws站计算的dtec。
[0045]
图4为本发明实施例2012年7月27日metu站计算的roti、iroti和dtec值，其中，图4(a)为metu站计算的roti，图4(b)为tids时期计算的iroti，图4(c)为epbs时期计算的iroti，图4(d)为tids时期计算的dtec和图4(e)为epbs时期计算的dtec。
[0046]
图5为本发明实施例tids时期和电离层平静时期roti和iroti的统计结果，其中，图5(a)为tids时期和电离层平静时期roti的统计结果，图5(b)为tids时期和电离层平静时期iroti的统计结果。
[0047]
图6为本发明实施例tids和epbs事件roti和iroti的统计结果。
具体实施方式
[0048]
本发明提供一种可识别电离层行波扰动和赤道等离子泡的探测指数，该指数在传统电离层扰动指数roti的基础上，结合了去趋势的总电子量(dtec)的优势，可以有效地探测出tids和epbs，并通过阈值将二者识别区分。
[0049]
下面结合附图和实施例对本发明的技术方案作进一步说明。
[0050]
如图1所示，本发明提供一种可识别电离层行波扰动和赤道等离子泡的探测指数，包括以下步骤：
[0051]
步骤1，获取观测区域的gnss原始观测文件和广播星历文件。
[0052]
步骤2，对gnss原始观测数据进行周跳探测，获得具有连续观测弧段的历元编号与卫星prn号。
[0053]
利用turbedit方法对步骤1获得的原始观测数据进行周跳探测，得到具有连续观测弧段的历元编号与卫星prn号，并将最短连续弧段设置为300历元以确保后续滑动平均获取去趋势的总电子量dtec能够正常进行。
[0054]
步骤3，基于gnss原始观测数据构建无几何相位观测方程，获得包含伪距和相位硬件延迟偏差的电离层总电子含量stec。
[0055]
gnss信号由伪距和载波相位观测值两部分组成，可以表示为：
[0056][0057]
式中，下标r和上标s分别表示gnss接收机和卫星的序列号，下标i表示频率，和为卫星与接收机之间沿视线(los)方向的伪距和相位观测值，表示卫星与接收机之间的几何距离，c是真空中的光速，dtr和dts分别代表gnss接收机和卫星的钟差，为沿los方向的对流层延迟，为频率i处的电离层延迟，b
r,i
和分别表示gnss接收机和卫星在频率i处的伪距硬件延迟，b
r,i
和分别表示gnss接收机和卫星在频率i处的相位小数偏差，λi为波长，为相位整周模糊度，和是观测噪声和其他各种误差的组合项(如伪距和载波相位观测的多路径效应)。
[0058]
由于伪距观测值更容易受到多路径效应和观测噪声的影响，从而导致电离层tec
精度下降。因此，利用无几何相位组合方法可以提取沿los方向的电离层观测值，比常用的相位平滑伪距方法分析得到的电离层扰动精度更高。无几何相位观测方程通过不同频率的载波相位观测值差值得到，即：
[0059][0060]
式中，为无几何相位组合观测值，和分别为f1、f2频率下沿los方向的相位观测值，为沿给定los方向的电离层延迟，γi＝(f0/fi)2是比例因子，其中f0＝10.23mhz，λ1和λ2分别为f1、f2频率的波长，和分别为f1、f2频率的整周模糊度，表示多径效应、噪声误差和其他类型的误差。
[0061]
电离层总电子含量(stec)定义为从卫星到接收机沿los方向的电子密度的积分，频率为f的stec与电离层延迟i在长度单位上的关系可以描述为：
[0062][0063]
式中，stec值以tecu为单位表示，1tecu＝10
16
el/m2；为无几何相位组合观测值；γi＝(f0/fi)2是比例因子，其中f0＝10.23mhz；表示多径效应、噪声误差和其他类型的误差。stec的截止高度角设置为30
°
，以削弱多路径效应的影响。
[0064]
无几何相位组合不能分离电离层延迟和相位的整周模糊度，但相位观测值在连续观测弧段的整周模糊度是不变的，不同频率信号的整周模糊度的差值是一个常量，在后续的扰动探测指数计算中可以消除，这也正是步骤2计算周跳探测的原因。
[0065]
步骤4，利用步骤2得到的连续观测弧段的历元编号和卫星prn号对步骤3得到的stec进行筛选，得到处于同一连续弧段的stec序列。
[0066]
步骤5，利用处于同一连续弧段的stec序列计算电离层扰动指数rot指数，进而获得roti指数的时间序列。
[0067]
rot指数为两个连续时间期间stec值的时间导数，即：
[0068][0069]
式中，rot值作为测量相位波动活动的参数，以tecu/min为单位表示；分别为卫星k在i-1、i时刻的电离层总电子含量，t
i-1
、ti为历元时间，k为卫星序列号。
[0070]
roti为特定时期内rot的标准差，具体计算方式如下：
[0071][0072]
本实施例将特定时期设定为5分钟，尖括号表示5分钟时间间隔内的平均值。
[0073]
步骤6，采用滑动平均法对stec进行去趋势处理，得到去趋势的总电子量(dtec)的时间序列。
[0074]
滑动平均法计算dtec的公式如下：
[0075][0076]
式中，stec(i)为第i时刻的stec值；n为滑动窗口，本实施例设置为30分钟；stec(j)为第j时刻的stec值，j的取值范围为[i-n/2,i+n/2]。
[0077]
步骤7，联合步骤5获得的roti指数和步骤6获得的dtec，计算iroti指数。
[0078]
图2展示了一组使用人工模拟tids和epbs信号的dtec时间序列图。图2(a)为tids的模拟信号，主要由平静时期的波状特征信号(w1)和扰动期的tids(d1)信号组成。图2(b)为epbs的模拟信号，主要由平静期的波状特征信号(w2)和扰动时期的epbs信号(d2)组成。从图2(a)可知，在发生tids的5分钟时间窗口(d1)的大部分时间内，dtec的时间导数值都相对接近。这说明发生tids时期的rot基本相同，从而导致roti处于较低水平，这就是传统的roti方法无法有效识别tids的原因。由于不论是tids还是epbs的发生时期dtec值总是大于平静时期的dtec值，因此，本发明通过计算dtec值在特定时期内的积分作为判断扰动发生的指数idtec(图2(a)中的灰色区域)：
[0079][0080]
式中，dtec(i)为i时刻的dtec值，δt为离散dtec的时间间隔，n为5分钟内dtec个数，idtec值通常以tecu
·
min为单位表示。
[0081]
对比图2(a)和图2(b)，tids和epbs(图2中灰色区域)的idtec值都处于较高水平，因此无法区分这两种形式的扰动。但从图2(b)可以看出，epbs发生时5分钟时间窗口内dtec时间导数值差异较大，说明发生epbs时roti处于较高水平。因此在探测tids和epbs的基础上结合roti可以进一步区分两者。
[0082]
本发明结合idtec和roti提出了iroti指数，该指数考虑振幅和相位的变化来探测电离层扰动进而区分tids和epbs，具体计算公式可以表示为：
[0083]
iroti＝idtec
·
roti(8)
[0084]
式中，iroti值以tecu2为单位表示。
[0085]
当iroti＜4.91tecu2时，可以判断电离层扰动类型为tids，当iroti＞7.28tecu2时，可以认为扰动类型为epbs。
[0086]
图3展示了2015年6月22日hkws站计算的roti、iroti和dtec结果。这天在香港地区观测到了由强地磁暴引起的大尺度tids事件。从图3(c)中可以看出，在22:00ut到22:30ut之间，dtec值呈现了不同振幅的波状结构(tids)。通常当roti值小于0.25tecu/min时，电离层可视为处于平静状态。但从图3(a)可以看出，发生tids时期该站的roti值并没有明显的变化，且均小于0.2tecu/min，这说明roti无法用来对tids进行检测。从图3(b)可知，与roti相比，iroti值在tids发生期间显著增加，并且可以与电离层平静期明显区分，这说明iroti可以用于探测tids。
[0087]
图4展示了2012年7月27日metu站上空观测到tids和epbs事件的roti、iroti和dtec值。从图4(d)可得，g12卫星的dtec值在11:40左右呈现了持续时间近1小时的波状特征(tids)。对比图4(a)和图4(b)，在发生tids的时期，iroti值从0.06tecu2上升到0.65tecu2，，而roti的整体趋势均小于0.25tecu/min，这说明roti不能检测到tids。在图4中左边三个箭
头所指的时间点中，通过图4(d)可以判断此时并没有发生较为明显的扰动情况，但此时roti却急剧上升至0.36tecu/min，而iroti相对稳定，保持在0.2tecu以下的区间。此时的roti值突然升高可能是由于中午剧烈的太阳活动(lt＝12.4h)导致局部电离层不稳定造成的。这表明iroti在tids探测中具有较好的抗干扰能力，可以较好地避免中午电离层活跃的影响。从图4(e)可以看出，所有卫星的dtec值在19点25分左右出现了不同程度的剧烈波动，并且在相应的时间段内图4(a)中的roti值也达到了0.89tecu/min，根据roti的值将其判断为epbs。在epbs事件发生时，图4(c)中iroti的值也与roti一样在同一时间有不同程度的增大，这说明iroti也可以用于检测epbs。
[0088]
图5展示了发生tids事件时期和电离层平静时期roti和iroti的统计结果。其中实心点和空心点分别代表tids事件和平静期相关指数的最大值。从图5(a)可得，电离层平静时期roti的最大值在0.06～0.31tecu/min，扰动期roti最大值在0.08～0.63tecu/min，与“平静”期有大面积重叠(图5(a)中的灰色部分)。因为大部分电离层正常变化的roti也能到达该区域，所以仅利用roti指标无法成功探测重叠区域的tids事件。图5(a)中绝大多数的tids事件都处于重叠区域，说明使用roti来监测tids存在一定的缺陷。而如图5(b)所示，平静时段iroti最大值在0.06～0.39tecu2之间，监测到tids事件的iroti最大值在0.22～7.28tecu2之间，两者的重叠面积很小(图5(b)中的灰色部分)，说明本发明提出的iroti指数相比于roti指数对tids的识别效果具有较大的改善。
[0089]
图6给出了选定时段内tids和epbs事件roti和iroti的统计结果，其中浅色点和深色点分别代表与epbs和tids事件相关指数的最大值。从图6可以看出，epbs事件发生时iroti的最大值在4.91～72.78tecu2区间内，总体上高于tids事件。灰色区域是tids和epbs事件的重叠区域。处于重叠区域的tids和epbs事件难以用iroti指数进行分类，属于模糊区间，该区间的电离层扰动事件可能存在同时发生tids和epbs等情况。模糊区间的上下界分别为iroti＝4.91tecu2和iroti＝7.28tecu2，可以被认为是区分tids和epbs的阈值区间。当iroti＜4.91tecu2时，可以判断电离层扰动类型为tids，当iroti＞7.28tecu2时，可以认为扰动类型为epbs。
[0090]
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施案例，做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

技术特征：
1.一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于，包括以下步骤：步骤1，获取观测区域的gnss原始观测文件和广播星历文件；步骤2，对gnss原始观测数据进行周跳探测，获得具有连续观测弧段的历元编号与卫星prn号；步骤3，基于gnss原始观测数据构建无几何相位观测方程，获得包含伪距和相位硬件延迟偏差的电离层总电子含量stec；步骤4，利用步骤2得到的连续观测弧段的历元编号和卫星prn号对步骤3得到的stec进行筛选，得到处于同一连续弧段的stec序列；步骤5，利用处于同一连续弧段的stec序列计算电离层扰动指数rot指数，进而获得roti指数的时间序列；步骤6，采用滑动平均法对stec进行去趋势处理，得到去趋势的总电子量(dtec)的时间序列；步骤7，联合步骤5获得的roti指数和步骤6获得的dtec，计算iroti指数。2.如权利要求1所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤2中利用turbedit方法对步骤1获得的原始观测数据进行周跳探测，得到具有连续观测弧段的历元编号与卫星prn号，并将最短连续弧段设置为μ1历元以确保后续滑动平均获取dtec能够正常进行。3.如权利要求1所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤3中gnss信号由伪距和载波相位观测值两部分组成，可以表示为：式中，下标r和上标s分别表示gnss接收机和卫星的序列号，下标i表示频率，和为卫星与接收机之间沿视线(los)方向的伪距和相位观测值，表示卫星与接收机之间的几何距离，c是真空中的光速，dt
r
和dt
s
分别代表gnss接收机和卫星的钟差，为沿los方向的对流层延迟，为频率i处的电离层延迟，b
r,i
和分别表示gnss接收机和卫星在频率i处的伪距硬件延迟，b
r,i
和分别表示gnss接收机和卫星在频率i处的相位小数偏差，λ
i
为波长，为相位整周模糊度，和是观测噪声和其他各种误差的组合项。4.如权利要求3所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤3中无几何相位观测方程通过不同频率的载波相位观测值差值得到，即：式中，为无几何相位组合观测值，和分别为f1、f2频率下沿给定los方向的相位观测值，为沿给定los方向的电离层延迟，γ
i
＝(f0/f
i
)2是比例因子，其中f0＝10.23mhz，λ1和λ2分别为f1、f2频率的波长，和分别为f1、f2频率的整周模糊度，表示多径效应、噪声误差和其他类型的误差。
5.如权利要求4所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤3中电离层总电子含量(stec)定义为从卫星到接收机沿los方向电子密度的积分，频率为f的stec与电离层延迟i在长度单位上的关系可以描述为：式中，stec值以tecu为单位表示，1tecu＝10
16
el/m2；为无几何相位组合观测值；γ
i
＝(f0/f
i
)2是比例因子，其中f0＝10.23mhz；表示多径效应、噪声误差和其他类型的误差；stec的截止高度角设置为μ2，以削弱多路径效应的影响。6.如权利要求1所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤5中rot指数为两个连续时间期间stec值的时间导数，即：式中，rot值作为测量相位波动活动的参数，以tecu/min为单位表示；分别为卫星k在i-1、i时刻的电离层总电子含量，t
i-1
、t
i
为历元时间，k为卫星序列号；roti为特定时期内rot的标准差，具体计算方式如下：式中，尖括号表示特定时期内的平均值。7.如权利要求1所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤6中滑动平均法计算去趋势的总电子量(dtec)的公式如下：式中，stec(i)为第i时刻的stec值，n为滑动窗口，stec(j)为第j时刻的stec值，j的取值范围为[i-n2,i+n2]。8.如权利要求1所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤7中通过计算dtec值在特定时期内的积分作为判断扰动发生的指数idtec：式中，dtec(i)为i时刻dtec的绝对值，δt为离散dtec的时间间隔，n为特定时期内dtec个数，idtec值以tecu
·
min为单位表示。9.如权利要求7所述的一种可识别电离层行波扰动和赤道等离子泡的探测指数，其特征在于：步骤7中结合idtec和roti提出了iroti指数，该指数考虑振幅和相位的变化来探测电离层扰动进而区分tids和epbs，具体计算公式可以表示为：iroti＝idtec
·
roti(8)式中，iroti值以tecu2为单位表示；当iroti＜μ3tecu2时，可以判断电离层扰动类型为tids，当iroti＞μ4tecu2时，可以认为扰动类型为epbs。

技术总结
本发明公开了一种可识别电离层行波扰动和赤道等离子泡的探测指数。该指数在传统电离层扰动指数ROTI的基础上，结合了去趋势的总电子量(DTEC)的优势，可以有效地探测出TIDs和EPBs，并通过阈值将二者识别区分。并通过阈值将二者识别区分。并通过阈值将二者识别区分。


技术研发人员：任晓东 勒烜 张小红
受保护的技术使用者：武汉大学
技术研发日：2023.03.06
技术公布日：2023/9/14