CN111983654A - 一种基于gnss的北极区域电离层相位闪烁因子构建方法 - Google Patents

Abstract

一种基于GNSS的北极区域电离层相位闪烁因子构建方法，用于北极区域天文灾害监测使用。首先选取位于北极区域的原始GNSS相位观测值，通过模型的方法改正其中的卫星钟、接收机相位中心、相对论、相位缠绕、对流层延迟和固体潮误差，利用精密单点定位方法进一步消除接收机钟差和未能用模型改正的对流层湿延迟，利用离散小波变换技术降低信号的噪声，降低噪声对电离层闪烁信号的提取；利用连续小波变换技术对上述信号进行时频分析，提取经验信号频带范围，并提取出电离层闪烁信号；设定平移窗口，对窗口内的电离层闪烁信号取标准差，完成电离层相位闪烁因子构建。其步骤简单，检测成本低，具有广泛的实用性。

Description

一种基于GNSS的北极区域电离层相位闪烁因子构建方法

技术领域

本发明涉及一种电离层闪烁因子构建方法，尤其适用于一种基于GNSS的北极区域电离层相位闪烁因子构建方法，用于北极区域低电离层闪烁天文灾害的监测、建模及预报。

技术背景

电离层闪烁是北极区域频发的一种天文灾害，是近年来空间科学面临的一大难题，给GNSS系统位置、导航与时间(PNT)服务的稳定性带来巨大挑战，例如2003年的万圣节电离层闪烁事件造成GNSS导航卫星中断服务，严重干扰了电力、银行系统以及北极航线的航行安全，造成经济损失。由此可见，电离层闪烁灾害对人类社会可持续发展构成了重大威胁，科学认识电离层闪烁发生规律，监测及预报电离层闪烁天文灾害对于防灾减灾具有重大现实意义。

电离层闪烁会对GNSS信号带来振幅或者相位的剧烈波动。电离层闪烁因子是对该影响的定量表征，可以反映出电离层闪烁的强弱，是实现电离层闪烁的监测、建模及预报的基础数据，也是实时修正电离层闪烁造成的定位误差的重要参数。为了测量电离层闪烁对GNSS信号的影响，通常需要电离层闪烁监测接收机(Ionospheric ScintillationMonitoring Receivers(ISMR))，其可以直接给出两类闪烁因子：振幅闪烁因子和相位闪烁因子。ISMR通常运行在50Hz的采样频率，导致其需要较大的存储空间，且价格较贵，进而限制了ISMR的布站数量。当前布设ISMR的电离层闪烁监测网络主要有欧洲航天局的Monitor项目、巴西的CIGALA/CALIBRA、ICEA和LISN项目和加拿大的CHAIN项目，共105个测站，其中仅有27个测站位于北极区域。由此可见，ISMR网络尚未对北极区域形成全覆盖，无法利用ISMR所提供的闪烁因子进行北极地区电离层闪烁的建模。相比于ISMR，采样频率通常在1Hz及以下的测地型接收机分布十分广泛，在北极区域布设测地型接收机的观测站超过300个，基本可以实现对北极陆地及近海的全覆盖，且这些测站的设立时间更久，可以提供更长时间序列的观测数据，有助于所建立的电离层闪烁模型吸收电离层闪烁变化的长周期项的影响。因此，构建基于测地型接收机的电离层闪烁因子，有利于提高所构建北极区域电离层闪烁模型的精度。

目前，基于低采样频率的测地型接收机开发的电离层闪烁因子主要有四个，即总电子含量变化率指数(ROTI)、沿弧段总电子含量变化率指数(AATR)亦被称做TEC起伏、无电离层组合残差的标准差(σ_IF)和基于C/N₀的闪烁因子(S_4c)，其中S_4c被验证与电离层振幅闪烁因子具有较强时空相关性，而前三种闪烁因子主要与相位闪烁因子相关。由于电离层闪烁在低纬度区域和北极区域的产生机理不同，这导致发生在低纬度区域的电离层闪烁主要以振幅闪烁为主，而北极区域主要以相位闪烁为主，利用相位闪烁因子可以较好反映出北极区域电离层闪烁信息，所以这里将重点讨论前三种闪烁因子。测地型接收机与ISMR主要区别在于ISMR所采用的接收机钟更精确、更稳定。过大的钟差会淹没载波观测值中的电离层闪烁信号。为了消除测地型接收机中波动的钟差，ROTI和AATR均采用的历元间差分方法，该方法虽然可以消除一定的接收机钟差，但却降低了闪烁因子的时间分辨率，导致无法获得每个历元的闪烁因子；而σ_IF采用了无电离层组合估计接收机钟差的方法，克服了历元间做差的不足。为了检验测地型接收机给出的电离层闪烁因子的准确性，当前已有的研究验证了ROTI和AATR在电离层闪烁期间与ISMR给出的相位和振幅闪烁因子在低纬度和极区具有较强相关性，AATR已被欧洲地球静止导航重叠服务(EGNOS)选为表征电离层活动的参数之一。Juan等利用低纬度的数据无法验证σ_IF与ISMR提供的相位闪烁因子的相关性，因为σ_IF给出的是相位闪烁在两个频点共同作用的结果，而在相位闪烁对每个频点的影响是不同的。针对σ_IF无法提供每个频点的相位闪烁因子的缺点，Nguyen等引入0.1Hz截至频率的巴特沃斯滤波剔除载波信号中电离层折射作用影响，进而获得每个频点的闪烁因子，该闪烁因子采用的0.1Hz截至频率是依据低纬地区大量观测值而确定的菲涅耳频率(Fresnelfrequency)的经验值，北极地区的菲涅耳频率通常更高，而采用更高截至频率的巴特沃斯滤波会极大的削弱所估计闪烁因子的幅值，导致该滤波方法难以适用于北极。由此可见，虽然利用上述研究成果可以使得测地型接收机应用到电离层监测当中，但是目前针对构建北极地区的电离层闪烁因子方法尚缺乏深入研究。

发明内容

针对上述技术问题，提供了一种基于GNSS的北极区域电离层相位闪烁因子构建方法，可以给出GNSS每一频率上闪烁因子，利用电离层相位闪烁因子将分布较为广泛但采样频率较低的测地型接收机引入到北极区域电离层闪烁监测领域，能够解决因北极区域电离层闪烁接收机分布不足而导致的无法完备监测北极区域电离层闪烁天文灾害问。

为实现上述技术目的，本发明的基于GNSS的北极区域电离层相位闪烁因子构建方法，首先选取位于北极区域的低采样频率观测信号的原始GNSS相位观测值，结合国际GNSS服务(IGS)提供的精密星历和测站坐标进行站星几何距离改正，通过模型的方法改正其中的卫星钟、接收机相位中心、相对论、相位缠绕、对流层延迟和固体潮误差，利用精密单点定位方法进一步消除接收机钟差和未能用模型改正的对流层湿延迟，利用离散小波变换技术降低信号的噪声，降低噪声对电离层闪烁信号的提取；利用连续小波变换技术，对上述信号进行时频分析，提取出北极区域电离层闪烁对应的经验信号频带范围，并对频带范围内的信号进行连续小波的逆变换，提取出电离层闪烁信号；设定一定长度的平移窗口，并对窗口内的电离层闪烁信号取标准差，完成基于GNSS低采样频率观测信号的电离层相位闪烁因子构建。

具体步骤为：

步骤1：准确确定测站坐标和每个观测历元的卫星坐标，之后改正改正每颗卫星至测站之间的几何距离：首先可通过公开的位于北极区域的GNSS数据观测网络获得GNSS原始观测数据；利用CPRS在线解算软件通过精密单点定位静态解算方法获得测站坐标，通过IGS提供的精密星历文件中记载的卫星坐标获取每个观测历元的卫星坐标，时间分辨率为15分钟，采用二阶拉格朗日算法内插出每个观测历元时刻的卫星坐标，根据测站坐标和每个观测历元的截止高度角在15°以上的卫星坐标，通过欧几里得度量获得卫星至测站几何距离改正参数，根据GNSS观测方程改正北极区域低采样频率GNSS原始载波相位观测值中卫星至测站几何距离，获得修正几何距离之后的观测值残差；

步骤2：由于经过站星几何距离改正后的GNSS原始载波相位观测值依旧存在多种误差，因此采用模型改正的方式进行修正，具体包括固体潮、天线相位中心、卫星钟差、相对论、相位缠绕、对流层干延迟和部分湿延迟误差的修正：

步骤3：对流层湿延迟精细改正和测站接收机钟差改正：利用精密单点定位(PPP)技术，将未通过建模改正部分的对流层湿延迟和测站接收机钟差作为精密单点定位中的待求参数，对二者进一步改正：

利用模型改正后GNSS原始观测信号组成双频无电离层组合，将该组合作为精密单点定位的基础观测值，组成如下观测方程：

V＝Ax-(L-D)

x＝[δX，δY，δZ，cδt_R，δρ_z，wc，B¹…Bⁱ…B^N]^T

其中，A为系数矩阵，L为无电离层组合观测值，D为无电离层组合模型改正值，x为待解算参数包括测站三维坐标(δX，δY，δZ)、δρ_z,wc为天顶对流层湿延迟、cδt_R为接收机钟差，B为模糊度，该观测方程采用卡尔曼滤波静态解算方式进行求解，滤波中的状态转移矩阵Φ为：

Φ＝diag(1，1，1，0，1，1，…，1)

系统噪声向量的协方差矩阵Q为：

Q＝diag(0，0，0，9×10¹⁰，10^-4，0，…，0)

采用正向和反向运算相结合的方式提高天顶对流层湿延迟和测站接收机钟差的估计精度，并将反向运算的结果作为未能利用模型改正的天顶对流层湿延迟δρ_z,wc和测站接收机钟差cδt_R的估计结果；通过Niell投影函数并结合GNSS观测方程将获得的天顶对流层延迟δρ_z,wc改正到每一卫星的观测值中；

步骤4：通过历元间做差的方法准确确定接收机钟差出现钟跳的发生历元，当钟跳存在时，历元间做差后的载波相位观测值会存在毫秒级的阶跃，利用精密单点定位方法无法准确估计出接收机钟差变化的细节部分，因而对经过站星几何距离改正、模型改正、对流层湿延迟改正和接收机钟差初步改正的GNSS载波相位观测值进行钟差的进一步精细改正；

步骤5：对经过上述步骤改正后的GNSS载波相位观测值利用离散小波变换的方法剔除观测噪声：采用现有的小波去噪经验信息：小波基选为多贝西极限相位小波，消失矩为2，小波分解层数为5层，降噪阈值确定方法为考虑信噪比的Stein无偏似然估计和固定阈值估计折中的原则，即信噪比很小时，按照Stein无偏似然估计处理时，信号噪声很大，在这种情况下，采用固定阈值处理，计算方法为

其中n为数据的长度；

步骤6：利用连续小波变换提取闪烁对应频带区间及构建闪烁因子：对GNSS载波相位观测值经过上述步骤改正后，其残差为电离层缓慢变化的部分和电离层闪烁部分，利用连续小波变换的方法将信号从时间域变换到时间-频率域时，通过时频分析，定位出突变信号即电离层闪烁信号的发生时刻及所在的特征频率区间，最终确定出Morse小波的对称参数和时间带宽积；利用以上参数对残差进行连续小波变换，并对位于特征频率区间的小波参数进行连续小波逆变换，获得电离层闪烁信号；以60s的时间间隔为滑动窗口，对电离层闪烁信号求标准差，获得基于GNSS低采样频率观测数据的电离层相位闪烁因子。

采用模型改正的方式对固体潮、天线相位中心、卫星钟差、相对论、相位缠绕、对流层干延迟和部分湿延迟误差的修正具体包括：

1)采用二阶简化潮汐模型针对GNSS载波相位观测值中的固体潮误差进行改正；

2)利用IGS提供的包含接收机和卫星绝对相位中心改正数据(即ANTEX文件)改正站的接收机天线相位中心(包括接收机天线参考点、接收机天线平均相位中心和接收机天线瞬时相位中心)和卫星天线相位中心(包括卫星天线平均相位中心和卫星天线瞬时相位中心)；

3)利用IGS提供的精密钟差文件改正截止高度角15°以上的每一个卫星的钟差；

4)针对相对论效应影响的改正包括两部分：

一部分是使得卫星钟频率呈现周期性变化，采用如下公式改正：

其中，r^s和v^s分别为惯性系下卫星的位置矢量和速度矢量，c为光速；

另一部分为Shapiro延迟，即GNSS信号穿过地球引力场时出现的时空弯曲，采用如下公式改正：

其中，r^s、r^r分别为卫星和测站的地心距离，

为测站至卫星的距离，μ为地球引力常数；

5)采用如下公式对相位缠绕误差(ω)进行改正：

ω＝δφ+2Nπ

其中，N为整数部分，δφ为小数部分，采用如下公式计算：

其中，

ζ＝ρ′₀·(d×d)

其中，d和d′分别为测站和卫星的有效偶极矢量；ρ′₀为卫星指向测站的单位矢量；

6)对流层延迟(δρ_trop)分为干分量(δρ_z,d)和湿分量(δρ_z,w)两部分，采用如下公式估计：

δρ_trop(E,H)＝＝δρ_z,d·M_d(E,H)+δρ_z，w·M_w(E)

其中，E为卫星高度角；H为海拔高，本文用大地高近似替代，单位为米；M_d和M_w为对流层干分量采用如下方式估计：

δρ_z，d＝a·e^-βH

其中，a＝2.3m，β＝0.116×10^-3。对流层湿分量采用如下方式估计：

δρ_z，w＝δρ_z0,w+Δρ_z,w

由于对流层延迟湿分量难以利用模型精确估计，因此仅近似将δρ_z0,w取0.1m，Δρ_z,w作为随机游走参数进行估计，对于未估计的对流层延迟湿分量作为精密单点定位中的参数进行解算。

对经过站星几何距离改正、模型改正、对流层湿延迟改正和接收机钟差初步改正的GNSS载波相位观测值进行钟差的进一步精细改正具体步骤为：

将GNSS原始载波相位信号经过修正后的残差组成无电离层组合观测值，由于该组合可以消除绝大部分电离层的影响，故而经过修正后的残差仅剩余未能通过PPP技术估算的钟差的细节部分和整周模糊度，由于在同一观测弧段内，所有历元的整周模糊度具有同一性，但周跳会改变这种同一性，进而影响钟差的精细估计，为此首先需要对经过步骤1-3修正后的残差进行周跳的探测与修复，尤其是在电离层较为活跃的条件下，难以做到周跳的准确修复，因此采用双频宽项和无几何关系组合的办法进行初次周跳探测，这可以修复多数的周跳；未修复或未正确修复的周跳通过对经过步骤1-3修正后的残差组成的无电离层组合观测值在时间域做差的办法进行探测，由于该组合观测值可以极大的降低电离层的影响，可以提供精度较高的检测值，可以较为准确的检测出小周跳的存在，但只采用这一组合，无法确定周跳的大小，因此采用将探测发生周跳的历元作为新弧段的开始，这样可以保证在每一观测弧段内整周模糊度是不变的，对于每颗卫星每观测弧段内的无电离层组合观测值进行历元间做差，以消除整周模糊度的影响，做差之后的残差为钟差在时间上的变化梯度，对所有高度角大于30°卫星的钟差变化梯度在时间上进行加权求积分可以获得接收机钟差，其中定权方法采用卫星高度角定权的方法，即P＝sin²α，P为权重，α为卫星高度角。

通过对通常采用的0.1Hz至0.4Hz的特征频率区间的小波系数进行小波逆变换从而获得电离层闪烁信号

进而利用以下公式获得基于低采样频率数据的电离层相位闪烁因子

其中<·>表示以一定时间间隔内的期望，设60s为滑动窗口，因此构建的电离层相位闪烁因子的频率仍为1Hz。

选择Morse小波作为连续小波变换的小波基，Morse小波包括对称参数和时间带宽积两个参数，通过确定Morse小波的对称参数和时间带宽积参数从而准确提取出电离层闪烁信息：

设定Morse小波将位于3至6之间的对称参数以1为步长，位于20至60之间的时间带宽积以5为步长组成多种参数组合，利用不同参数组合的Morse小波对步骤1-5修正后的低采样频率的GNSS载波相位观测值进行连续小波变换，并对特征频率区间的小波系数进行连续小波逆变换，构建每种参数组合对应的相位闪烁因子；将由高频观测值得出的闪烁因子(σ_φ)与所构建的每种参数组合对应的闪烁因子进行比较，分析二者的相关性及残差，以能够提供高相关性且小残差值的参数组合更为适合对电离层闪烁信息的提取作为标准进行判断，最终确定出Morse小波的最佳对称参数为3，时间带宽积为45。

有益效果：

1.相比于已有适用于低采样频率数据的电离层闪烁因子，本发明所提出的方法可以给出GNSS每一频率上闪烁因子，便于综合分析电离层闪烁对GNSS不同频率造成的影响；

2.通过本发明构建出的电离层相位闪烁因子可以将分布较为广泛但采样频率较低的测地型接收机引入到北极区域电离层闪烁监测领域，可以助力解决因北极区域电离层闪烁接收机分布不足而导致的无法完备监测北极区域电离层闪烁天文灾害的技术难题。

相比于传统利用高采样频率(50Hz)的电离层闪烁监测接收机开展电离层闪烁监测研究，该方法可以将普通测地型接收机引入电离层闪烁监测领域，大幅度降低电离层闪烁监测研究的成本，结合测地型接收机分布广泛的特性，有利于对区域电离层闪烁的完备监测，具有广泛的实用性。

附图说明

图1为本发明实施的基于GNSS低采样频率观测信号的北极区域电离层相位闪烁因子构建方法的流程图。

图2为PRN25卫星的L1P观测值及该观测弧段的卫星高度角曲线图。

图3为利用PPP方法获得的钟差示意图。

图4(a)为对PRN25卫星的L1P信号改正后的残差示意图。

图4(b)为对PRN25卫星的L1P信号进一步精细修正接收机钟差后的残差示意图。

图5为利用离散小波变换方法剔除的噪声示意图。

图6为去噪后的残差的时频波普图。

图7(a)为使用本发明基于GNSS的北极区域电离层相位闪烁因子构建方法构建的相位闪烁因子示意图；

图7(b)为利用电离层闪烁接收机直接输出的相位闪烁因子示意图。

具体实施方式

下面结合附图对本发明的实施例做进一步说明：

如图1所示，本发明基于GNSS低采样频率观测信号的北极区域电离层相位闪烁因子构建方法，首先选取位于北极区域的低采样频率观测信号的原始GNSS相位观测值，结合国际GNSS服务(IGS)提供的精密星历和测站坐标进行站星几何距离改正，通过模型的方法改正其中的卫星钟、接收机相位中心、相对论、相位缠绕、对流层延迟和固体潮误差，利用精密单点定位方法进一步消除接收机钟差和未能用模型改正的对流层湿延迟，利用离散小波变换技术降低信号的噪声，降低噪声对电离层闪烁信号的提取；利用连续小波变换技术，对上述信号进行时频分析，提取出北极区域电离层闪烁对应的经验信号频带范围，并对频带范围内的信号进行连续小波的逆变换，提取出电离层闪烁信号；设定一定长度的平移窗口，并对窗口内的电离层闪烁信号取标准差，完成基于GNSS低采样频率观测信号的电离层相位闪烁因子构建。

具体步骤为：

步骤1：站星几何距离改正。为改正每颗卫星至测站之间的几何距离，需准确确定测站坐标和每个观测历元的卫星坐标。测站坐标利用CPRS在线解算软件通过精密单点定位静态解算方法获得。每个观测历元的卫星坐标通过IGS提供的精密星历文件中提供的卫星坐标获得，但其时间分辨率通常为15分钟，本发明采用二阶拉格朗日算法内插出每个观测历元时刻的卫星坐标。根据测站坐标和每个观测历元的卫星坐标，通过欧几里得度量获得站星几何距离改正参数，并在原始的载波相位观测值中改正站星几何距离。

步骤2：模型改正。GNSS原始载波相位观测值经过站星几何距离改正(步骤1)后，其中依旧存在多种误差，其中的固体潮、天线相位中心、卫星钟差、相对论、相位缠绕、对流层干延迟和部分湿延迟误差本发明采用模型改正的方式进行修正，具体如下。

1)固体潮误差采用二阶简化潮汐模型进行改正。

2)测站的接收机天线相位中心(包括接收机天线参考点、接收机天线平均相位中心和接收机天线瞬时相位中心)和卫星天线相位中心(包括卫星天线平均相位中心和卫星天线瞬时相位中心)利用IGS提供的包含接收机和卫星绝对相位中心改正数据(即ANTEX文件)进行改正。

3)卫星钟差利用IGS提供的精密钟差文件进行改正。

4)本发明对相对论效应的影响分成以下两部分，分别进行改正：

一部分是使得卫星钟频率呈现周期性变化，采用如下公式改正：

其中，r^s和v^s分别为惯性系下卫星的位置矢量和速度矢量，c为光速；

另一部分为Shapiro延迟，即GNSS信号穿过地球引力场时出现的时空弯曲，采用如下公式改正：

其中，r^s、r^r分别为卫星和测站的地心距离，

为测站至卫星的距离，μ为地球引力常数。

5)相位缠绕误差(ω)采用如下公式进行改正：

ω＝δφ+2Nπ

其中，N为整数部分，δφ为小数部分，采用如下公式计算：

其中，

ζ＝ρ′₀·(d′×d)

其中，d和d′分别为测站和卫星的有效偶极矢量；ρ′₀为卫星指向测站的单位矢量。

6)对流层延迟(δρ_trop)分为干分量(δρ_z，d)和湿分量(δρ_z，w)两部分，采用如下公式估计：

δρ_trop(E，H)＝δρ_z，d·M_d(E，H)+δρ_z，w·M_w(E)

其中，E为卫星高度角；H为海拔高，本文用大地高近似替代，单位为米；M_d和M_w为对流层干分量采用如下方式估计：

δρ_z，d＝a·e^-βH

其中，a＝2.3m，β＝0.116×10^-3。对流层湿分量采用如下方式估计：

δρ_z，w＝δρ_z0，w+Δρ_z，w

由于对流层延迟湿分量难以利用模型精确估计，本文仅近似将δρ_z0，w取0.1m，Δρ_z，w作为随机游走参数进行估计，对于未估计的对流层延迟湿分量作为精密单点定位中的参数进行解算。

步骤3：对流层湿延迟精细改正和测站接收机钟差改正。对流层湿延迟和测站接收机钟差无法通过模型精确改正(步骤2)，本发明利用精密单点定位(PPP)技术，将对流层湿延迟(未通过建模改正部分)和测站接收机钟差作为精密单点定位中的待求参数，对二者进行进一步的改正，具体如下。

对GNSS原始观测信号经过模型改正(步骤2)后，组成双频无电离层组合，将该组合作为精密单点定位的基础观测值，组成如下观测方程：

V＝Ax-(L-D)

x＝[δX，δY，δZ，cδt_R，δρ_z，wc，B¹…Bⁱ…B^N]^T

其中，A为系数矩阵，L为无电离层组合观测值，D为无电离层组合模型改正值，x为待解算参数包括测站坐标(δX，δY，δZ)、天顶对流层湿延迟(δρ_z,wc)、接收机钟差(cδt_R)和模糊度(B)。该观测方程可以采用卡尔曼滤波静态解算方式进行求解，滤波中的状态转移矩阵(Φ)为：

Φ＝diag(1，1，1，0，1，1，…，1)

系统噪声向量的协方差矩阵(Q)为：

Q＝diag(0，0，0，9×10¹⁰，10^-4，0，…，0)

为提高天顶对流层湿延迟和测站接收机钟差的估计精度，本发明采用正向和反向运算相结合的方式，并将反向运算的结果作为未能利用模型改正的天顶对流层湿延迟(δρ_z,wc)和测站接收机钟差(cδt_R)的估计结果。通过Niell投影函数，将获得的天顶对流层延迟(δρ_z,wc)改正到每一卫星的观测值中。

步骤4：测站接收机钟中存在钟跳时的钟差精细估计。当接收机钟差漂移到某一阈值时，接收机厂商会通过对接收机钟差中插入钟跳的方式，使得接收机内部时钟与GPS时同步精度控制在一定范围。钟跳的幅度通常为1ms，其对接收机观测值得影响幅度远大于其他观测误差，所以可以通过历元间做差的方法准确确定钟跳的发生历元。当钟跳存在时，利用精密单点定位方法无法准确估计出接收机钟差变化的细节部分，因而需要对经过站星几何距离改正(步骤1)、模型改正(步骤2)、对流层湿延迟改正和接收机钟差初步改正(步骤3)的GNSS载波相位观测值进行钟差的进一步精细改正，具体如下。

将GNSS原始载波相位信号经过步骤1-3修正后的残差组成无电离层组合观测值，由于该组合可以消除绝大部分电离层的影响，故而经过步骤1-3修正后的残差仅剩余未能通过PPP技术估算的钟差的细节部分和整周模糊度。在同一观测弧段内，所有历元的整周模糊度具有同一性，但周跳会改变这种同一性，进而影响钟差的精细估计，为此首先需要对经过步骤1-3修正后的残差进行周跳的探测与修复。通常情况下，尤其是在电离层较为活跃的条件下，难以做到周跳的准确修复，因此本发明采用双频宽项和无几何关系组合的办法进行初次周跳探测，这可以修复多数的周跳；未修复或未正确修复的周跳通过对经过步骤1-3修正后的残差组成的无电离层组合观测值在时间域做差的办法进行探测，由于该组合观测值可以极大的降低电离层的影响，可以提供精度较高的检测值，可以较为准确的检测出小周跳的存在，但只采用这一组合，无法确定周跳的大小，本发明采用将探测发生周跳的历元作为新弧段的开始，这样可以保证在每一观测弧段内整周模糊度是不变的。对于每颗卫星每观测弧段内的无电离层组合观测值进行历元间做差，以消除整周模糊度的影响，做差之后的残差为钟差在时间上的变化梯度，对所有高度角大于30°卫星的钟差变化梯度在时间上进行加权求积分可以获得接收机钟差，其中定权方法采用卫星高度角定权的方法，即P＝sin²α，P为权重，α为卫星高度角。

步骤5：利用离散小波变换，剔除观测噪声影响。载波相位观测值中的观测噪声会干扰对其中电离层闪烁信息的提取，影响闪烁估计的准确性，因而需要对经过步骤1-4改正后的GNSS载波相位观测值去噪。本发明采用离散小波变换去噪的方法剔除观测噪声，通过实验确定了以下小波去噪的经验信息：小波基选为多贝西极限相位小波，消失矩为2，小波分解层数为5层，降噪阈值确定方法为考虑信噪比的Stein无偏似然估计和固定阈值估计折中的原则，即信噪比很小时，按照Stein无偏似然估计处理时，信号噪声很大，在这种情况下，采用固定阈值处理，计算方法为

其中n为数据的长度。

步骤6：利用连续小波变换提取闪烁对应频带区间。对GNSS载波相位观测值经过步骤1-4改正后，其残差主要包括两部分，即电离层缓慢变化的部分和电离层闪烁部分。本发明对该残差进行连续小波逆变换，变换的参数如下：小波基选为morse小波，其对称参数(symmetry parameter)设为3，时间带宽积(time-bandwidth product)设为60。通过连续小波变换可以获得经过步骤1-4改正后残差的时频谱图，其中电离层缓慢变化的部分对应低频部分，电离层闪烁部分主要为中频部分。本发明确定北极区域电离层闪烁对应的频带区间的经验值为0.003Hz至0.06Hz。

步骤7：对位于电离层闪烁经验频带范围内的信号进行连续小波逆变换，获得电离层闪烁对载波相位观测值得原始影响信号。以60s的时间间隔为滑动窗口，对窗口内的原始影响值求标准差，获得测地型接收机的电离层相位闪烁因子

具体实施例一、

利用加拿大高纬北极电离层监测网络(Canadian High Arctic IonosphericNetwork,CHAIN)的ARC监测站数据，对本发明做进一步详细说明。

ARC监测站位于北极圈内，设有Septentrio公司的PolaRx5S接收机，搭配PolaNtGG天线。接收机被设定为可以接收L1和L2频点的信号，采样频率为50Hz，可以直接提供电离层闪烁相位(σ_Φ)和振幅(S₄)闪烁因子。通过CHAIN提供的下载链接可以公开获得闪烁数据和采样频率为1Hz的GPS观测数据。

选取2020年第80天数据作为示例说明本发明的使用方法。

步骤1：将第80天ARC站观测数据上传至CPRS在线解算软件，采用静态精密单点定位的方式获得该站的准确坐标为北纬73°0'14.77639"、西经85°1'33.86753"、大地高为23.066m。选取12时至15时的PRN25卫星观测该时段L1频点的观测数据对本发明的步骤做示例说明，该观测数据和卫星高度角如图2所示。结合CPRS给出的精确坐标和IGS给出的精密星历文件，可以精确计算卫地距离。

步骤2：利用模型估算卫星钟、接收机相位中心、相对论、相位缠绕、对流层延迟和固体潮误差。

步骤3：利用精密单点定位技术估算未建模的对流层湿延迟和接收机钟差。

步骤4：接收机钟差精细修正。如图3所示，由于存在接收机钟跳，利用PPP获得的钟差难以精准的估算钟差的细微变化。若不顾及此细微误差，对原始L1P信号进行以上改正(包括步骤1、2、3中所涉及的改正)后，其结果较为波动，难以提取出电离层闪烁误差，如图4(a)所示。对GNSS原始载波相位信号经过步骤1、2和3修正之后的残差组成无电离层组合观测值，进一步精细估计接收机钟差。经过该钟差修正后，其载波相位残差结果，如图4(b)所示，精度有一定的提高，便于提取其中的电离层闪烁误差。

步骤5：利用离散小波变换技术对修正步骤1、2、3和4之后的信号残差进行去噪。提取出的观测噪声结果如图6所示，可以看出该方法可以去除0.5cm以内的观测噪声。

步骤6：对步骤5给出的去噪后的残差进行连续小波变换，其时频波普图如图7所示，可以看出，电离层闪烁发生的时间主要在13点45分左右，其对应的频带为0.003Hz至0.016Hz。

步骤7：对经验频带内的信号的小波系数进行连续小波逆变换，获得受电离层闪烁影响信号，并以60s的时间间隔为滑动窗口，对窗口内的值求标准差，获得接收机在低采样频率下的电离层闪烁相位闪烁因子

如图7(a)所示。对比利用高采样频率直接获得的电离层闪烁相位闪烁因子，如图7(b)可知，该方法可以较为有效的探测到电离层闪烁的发生，同时可以有效避免由于较低卫星高度角而产生的多路径的影响。

Claims (6)

1.一种基于GNSS的北极区域电离层相位闪烁因子构建方法，其特征在于：

本发明首先选取位于北极区域的低采样频率观测信号的原始GNSS相位观测值，结合国际GNSS服务(IGS)提供的精密星历和测站坐标进行站星几何距离改正，通过模型的方法改正其中的卫星钟、接收机相位中心、相对论、相位缠绕、对流层延迟和固体潮误差，利用精密单点定位方法进一步消除接收机钟差和未能用模型改正的对流层湿延迟，利用离散小波变换技术降低信号的噪声，降低噪声对电离层闪烁信号的提取；利用连续小波变换技术，对上述信号进行时频分析，提取出北极区域电离层闪烁对应的经验信号频带范围，并对频带范围内的信号进行连续小波的逆变换，提取出电离层闪烁信号；设定一定长度的平移窗口，并对窗口内的电离层闪烁信号取标准差，完成基于GNSS低采样频率观测信号的电离层相位闪烁因子构建。

2.根据权利要求1所述的基于GNSS的北极区域电离层相位闪烁因子构建方法，其特征在于具体步骤为：

步骤1：准确确定测站坐标和每个观测历元的卫星坐标，之后改正改正每颗卫星至测站之间的几何距离：首先可通过公开的位于北极区域的GNSS数据观测网络获得GNSS原始观测数据；利用CPRS在线解算软件通过精密单点定位静态解算方法获得测站坐标，通过IGS提供的精密星历文件中记载的卫星坐标获取每个观测历元的卫星坐标，时间分辨率为15分钟，采用二阶拉格朗日算法内插出每个观测历元时刻的卫星坐标，根据测站坐标和每个观测历元的截止高度角在15°以上的卫星坐标，通过欧几里得度量获得卫星至测站几何距离改正参数，根据GNSS观测方程改正北极区域低采样频率GNSS原始载波相位观测值中卫星至测站几何距离，获得修正几何距离之后的观测值残差；

步骤2：由于经过站星几何距离改正后的GNSS原始载波相位观测值依旧存在多种误差，因此采用模型改正的方式进行修正，具体包括固体潮、天线相位中心、卫星钟差、相对论、相位缠绕、对流层干延迟和部分湿延迟误差的修正：

步骤3：对流层湿延迟精细改正和测站接收机钟差改正：利用精密单点定位(PPP)技术，将未通过建模改正部分的对流层湿延迟和测站接收机钟差作为精密单点定位中的待求参数，对二者进一步改正：

利用模型改正后GNSS原始观测信号组成双频无电离层组合，将该组合作为精密单点定位的基础观测值，组成如下观测方程：

V＝Ax-(L-D)

x＝[δX,δY,δZ,cδt_R,δρ_z,wc,B¹…Bⁱ…B^N]^T

其中，A为系数矩阵，L为无电离层组合观测值，D为无电离层组合模型改正值，x为待解算参数包括测站三维坐标(δX，δY，δZ)、δρ_z,wc为天顶对流层湿延迟、cδt_R为接收机钟差，B为模糊度，该观测方程采用卡尔曼滤波静态解算方式进行求解，滤波中的状态转移矩阵Φ为：

Φ＝diag(1,1,1,0,1,1,…,1)

系统噪声向量的协方差矩阵Q为：

Q＝diag(0,0,0,9×10¹⁰,10^-4,0,…,0)

采用正向和反向运算相结合的方式提高天顶对流层湿延迟和测站接收机钟差的估计精度，并将反向运算的结果作为未能利用模型改正的天顶对流层湿延迟δρ_z,wc和测站接收机钟差cδt_R的估计结果；通过Niell投影函数并结合GNSS观测方程将获得的天顶对流层延迟δρ_z,wc改正到每一卫星的观测值中；

步骤4：通过历元间做差的方法准确确定接收机钟差出现钟跳的发生历元，当钟跳存在时，历元间做差后的载波相位观测值会存在毫秒级的阶跃，利用精密单点定位方法无法准确估计出接收机钟差变化的细节部分，因而对经过站星几何距离改正、模型改正、对流层湿延迟改正和接收机钟差初步改正的GNSS载波相位观测值进行钟差的进一步精细改正；

步骤5：对经过上述步骤改正后的GNSS载波相位观测值利用离散小波变换的方法剔除观测噪声：采用现有的小波去噪经验信息：小波基选为多贝西极限相位小波，消失矩为2，小波分解层数为5层，降噪阈值确定方法为考虑信噪比的Stein无偏似然估计和固定阈值估计折中的原则，即信噪比很小时，按照Stein无偏似然估计处理时，信号噪声很大，在这种情况下，采用固定阈值处理，计算方法为

其中n为数据的长度；

步骤6：利用连续小波变换提取闪烁对应频带区间及构建闪烁因子：对GNSS载波相位观测值经过上述步骤改正后，其残差为电离层缓慢变化的部分和电离层闪烁部分，利用连续小波变换的方法将信号从时间域变换到时间-频率域，通过时频分析，定位出突变信号即电离层闪烁信号的发生时刻及所在的特征频率区间，最终确定出Morse小波的对称参数和时间带宽积；利用以上参数对残差进行连续小波变换，并对位于特征频率区间的小波参数进行连续小波逆变换，获得电离层闪烁信号；以60s的时间间隔为滑动窗口，对电离层闪烁信号求标准差，获得基于GNSS低采样频率观测数据的电离层相位闪烁因子。

3.根据权利要求2所述的基于GNSS的北极区域电离层相位闪烁因子构建方法，其特征在于采用模型改正的方式对固体潮、天线相位中心、卫星钟差、相对论、相位缠绕、对流层干延迟和部分湿延迟误差的修正具体包括：

1)采用二阶简化潮汐模型针对GNSS载波相位观测值中的固体潮误差进行改正；

2)利用IGS提供的包含接收机和卫星绝对相位中心改正数据(即ANTEX文件)改正站的接收机天线相位中心(包括接收机天线参考点、接收机天线平均相位中心和接收机天线瞬时相位中心)和卫星天线相位中心(包括卫星天线平均相位中心和卫星天线瞬时相位中心)；

3)利用IGS提供的精密钟差文件改正截止高度角15°以上的每一个卫星的钟差；

4)针对相对论效应影响的改正包括两部分

一部分是使得卫星钟频率呈现周期性变化，采用如下公式改正：

其中，r^s和v^s分别为惯性系下卫星的位置矢量和速度矢量，c为光速；

另一部分为Shapiro延迟，即GNSS信号穿过地球引力场时出现的时空弯曲，采用如下公式改正：

其中，r^s、r^r分别为卫星和测站的地心距离，

为测站至卫星的距离，μ为地球引力常数；

5)采用如下公式对相位缠绕误差(ω)进行改正：

ω＝δφ+2Nπ

其中，N为整数部分，δφ为小数部分，采用如下公式计算：

其中，

ζ＝ρ′₀·(d′×d)

其中，d和d′分别为测站和卫星的有效偶极矢量；ρ′₀为卫星指向测站的单位矢量；

6)对流层延迟(δρ_trop)分为干分量(δρ_z,d)和湿分量(δρ_z,w)两部分，采用如下公式估计：

δρ_trop(E,H)＝δρ_z,d·M_d(E,H)+δρ_z,w·M_w(E)

其中，E为卫星高度角；H为海拔高，本文用大地高近似替代，单位为米；M_d和M_w为对流层干分量采用如下方式估计：

δρ_z,d＝a·e^-βH

其中，a＝2.3m，β＝0.116×10^-3。对流层湿分量采用如下方式估计：

δρ_z,w＝δρ_z0,w+Δρ_z,w

由于对流层延迟湿分量难以利用模型精确估计，因此仅近似将δρ_z0,w取0.1m，Δρ_z,w作为随机游走参数进行估计，对于未估计的对流层延迟湿分量作为精密单点定位中的参数进行解算。

4.根据权利要求2所述的基于GNSS的北极区域电离层相位闪烁因子构建方法，其特征在于：对经过站星几何距离改正、模型改正、对流层湿延迟改正和接收机钟差初步改正的GNSS载波相位观测值进行钟差的进一步精细改正具体步骤为：

将GNSS原始载波相位信号经过修正后的残差组成无电离层组合观测值，由于该组合可以消除绝大部分电离层的影响，故而经过修正后的残差仅剩余未能通过PPP技术估算的钟差的细节部分和整周模糊度，由于在同一观测弧段内，所有历元的整周模糊度具有同一性，但周跳会改变这种同一性，进而影响钟差的精细估计，为此首先需要对经过步骤1-3修正后的残差进行周跳的探测与修复，尤其是在电离层较为活跃的条件下，难以做到周跳的准确修复，因此采用双频宽项和无几何关系组合的办法进行初次周跳探测，这可以修复多数的周跳；未修复或未正确修复的周跳通过对经过步骤1-3修正后的残差组成的无电离层组合观测值在时间域做差的办法进行探测，由于该组合观测值可以极大的降低电离层的影响，可以提供精度较高的检测值，可以较为准确的检测出小周跳的存在，但只采用这一组合，无法确定周跳的大小，因此采用将探测发生周跳的历元作为新弧段的开始，这样可以保证在每一观测弧段内整周模糊度是不变的，对于每颗卫星每观测弧段内的无电离层组合观测值进行历元间做差，以消除整周模糊度的影响，做差之后的残差为钟差在时间上的变化梯度，对所有高度角大于30°卫星的钟差变化梯度在时间上进行加权求积分可以获得接收机钟差，其中定权方法采用卫星高度角定权的方法，即P＝sin²α，P为权重，α为卫星高度角。

5.根据权利要求2所述的基于GNSS的北极区域电离层相位闪烁因子构建方法，其特征在于：通过对通常采用的0.1Hz至0.4Hz的特征频率区间的小波系数进行小波逆变换从而获得电离层闪烁信号

进而利用以下公式获得基于低采样频率数据的电离层相位闪烁因子

其中<·>表示以一定时间间隔内的期望，设60s为滑动窗口，因此构建的电离层相位闪烁因子的频率仍为1Hz。

6.根据权利要求2所述的基于GNSS的北极区域电离层相位闪烁因子构建方法，其特征在于：选择Morse小波作为连续小波变换的小波基，Morse小波包括对称参数和时间带宽积两个参数，通过确定Morse小波的对称参数和时间带宽积参数从而准确提取出电离层闪烁信息：

设定Morse小波将位于3至6之间的对称参数以1为步长，位于20至60之间的时间带宽积以5为步长组成多种参数组合，利用不同参数组合的Morse小波对步骤1-5修正后的低采样频率的GNSS载波相位观测值进行连续小波变换，并对特征频率区间的小波系数进行连续小波逆变换，构建每种参数组合对应的相位闪烁因子；将由高频观测值得出的闪烁因子(σ_φ)与所构建的每种参数组合对应的闪烁因子进行比较，分析二者的相关性及残差，以能够提供高相关性且小残差值的参数组合更为适合对电离层闪烁信息的提取作为标准进行判断，最终确定出Morse小波的最佳对称参数为3，时间带宽积为45。

Images