测绘学报 | 赵东升, 李旺, 李宸栋, 等：1 Hz GNSS电离层相位闪烁因子提取及在北极区域的验证

本文内容来源于《测绘学报》2021年第3期（审图号GS（2021）1180号）

1 Hz GNSS电离层相位闪烁因子提取及在北极区域的验证

赵东升¹,李旺¹,李宸栋²,唐旭³,张克非¹

1.中国矿业大学环境与测绘学院, 江苏 徐州 221116;
2.宁波诺丁汉大学理工学院, 浙江 宁波 315100;
3.南京信息工程大学遥感与测绘学院, 江苏 南京 210044

基金项目：国家自然科学基金(41730109；41874040；42074226)；江苏省自然科学基金(BK20200664；BK20200646)；中央高校基本科研业务费(2020QN30；2020QN31)；徐州市重点项目(KC19111)；江苏省优势学科建设工程项目(PAPD)；江苏省双创博士项目

摘要：电离层闪烁作为北极区域频发的一种天文灾害，会影响全球导航卫星系统(GNSS)时空服务，需对其进行有效的监测。监测电离层闪烁通常需要高采样频率(50 Hz)的电离层闪烁接收机，但其分布有限，难以提供较大区域(如北极区域)的全面监测。为此本文以1 Hz GNSS观测数据为基础，详细研究了利用大地测量趋势分离、精密单点定位和小波变换等技术提取载波相位观测值中电离层闪烁信号的相关经验参数和方法，构建了可用于监测电离层闪烁的相位闪烁因子，并利用加拿大高纬度北极电离层闪烁监测网络(CHAIN)提供的11个测站188 d的观测数据对其有效性和准确性进行了验证。结果表明，与常用的ROTI因子相比，本文所提出的闪烁因子与ROTI均能有效地探测到电离层闪烁的发生，但是本文提出的闪烁因子与电离层闪烁接收机直接给出的相位闪烁因子具有更好的相关性，尤其是在有强电离层闪烁发生的时段，这说明该闪烁因子具有更为优良的电离层闪烁监测能力。

关键词：全球导航卫星系统    电离层闪烁    相位闪烁因子    1 Hz采样频率数据    监测

引文格式：赵东升, 李旺, 李宸栋, 等. 1 Hz GNSS电离层相位闪烁因子提取及在北极区域的验证[J]. 测绘学报，2021，50(3)：368-383. DOI: 10.11947/j.AGCS.2021.20200454.

ZHAO Dongsheng, LI Wang, LI Chendong, et al.     Extracting an ionospheric phase scintillation index based on 1 Hz GNSS observations and its verification in the Arctic region[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(3): 368-383.   DOI: 10.11947/j.AGCS.2021.20200454.

阅读全文：http://xb.sinomaps.com/article/2021/1001-1595/2021-3-368.htm

引言

全球导航卫星系统(GNSS)已经成为监测电离层的重要手段^[1-4]。电离层闪烁是电离层中的不规则体，是近地磁赤道和地磁极区域频发的一种天文灾害^[5]，会对穿越其中的GNSS信号造成衍射和散射，进而使得GNSS信号的振幅或者相位剧烈波动，降低GNSS信号的信噪比甚至失锁，给GNSS系统位置、导航与时间(PNT)服务的稳定性带来巨大挑战^[6-8]。电离层闪烁因子是电离层闪烁对GNSS信号影响的定量表征，可以反映出电离层闪烁的强弱，是实现电离层闪烁的监测、建模、预报的基础数据^[9]，对于科学认识电离层闪烁发生规律、修正或减弱电离层闪烁对GNSS的不利影响进而实现高精度定位具有重要意义。

为了测量电离层闪烁对GNSS信号的影响，通常需要电离层闪烁监测接收机(ionospheric scintillation monitoring receivers，ISMR)，它可以直接给出两类闪烁因子：振幅闪烁因子(S₄)和相位闪烁因子(σ_ϕ)。由于ISMR运行在50 Hz的采样频率，使得ISMR可以直接跟踪到由于电离层闪烁而产生的信号振幅和相位的变换，但是这也导致其需要较大的存储空间，且价格较贵，限制了ISMR的布站数量。当前布设ISMR的电离层闪烁监测网络主要有欧洲航天局的Monitor项目、巴西的CIGALA/CALIBRA、ICEA和LISN项目和加拿大高纬度北极电离层闪烁监测网络(CHAIN)，共约100个测站。相比于ISMR，采样频率通常在1 Hz及以下的接收机分布十分广泛，基本可以实现对陆地及近海的全覆盖，且这些测站的设立时间更久，可以提供更长时间序列的观测数据，有助于研究电离层闪烁变化的长周期项，提高电离层闪烁监测及建模的精度^[10]，但在1 Hz及以下采样频率的观测数据中，由于观测数据时间分辨率的降低，电离层闪烁误差极容易被其他误差(尤其是接收机钟差)所淹没，为此需要专门的方法构建基于低采样频率(1 Hz及以下)GNSS观测数据的电离层闪烁因子。

地磁北极区域的电离层闪烁以相位闪烁为主^[11]，因而本文将以相位闪烁因子作为主要研究对象。为构建基于低采样频率观测数据的相位闪烁因子，文献[12]基于无几何组合观测值提出了总电子含量(TEC)变化率指数(ROTI)，可以较好减弱接收机钟差的影响，已经成为电离层闪烁监测领域应用最为广泛的标准闪烁因子之一^[12-14]，ROTI可以给出每个观测历元的闪烁因子。在此基础上，文献[15]通过对一段时间内所有可视卫星获得的TEC值以卫星高度角定权的方式进行加权平均从而提出了另一种闪烁因子——TEC起伏(AATR)，它可以用来预报发生闪烁的时间段和区域，已经被欧洲地球静止导航重叠服务(EGNOS)选为表征电离层活动的参数之一^[14]。文献[16]将无电离层组合残差的标准差作为监测电离层闪烁的因子并提出了对应的修正无电离层组合中的接收机钟差方法。文献[17]利用1 Hz的斜路径总电子含量和小波变换技术提出了一种类相位闪烁因子，该因子被认为可以有效地预警欧洲部分北极区域的电离层闪烁事件。上述研究所构建的闪烁因子均基于GNSS组合信号，其所获得的闪烁因子为电离层闪烁对组合信号作用的结果，无法用于研究电离层闪烁对每一频点信号的影响。

为构建基于1 Hz GNSS载波相位观测值的电离层相位闪烁因子，首先，本文将利用大地测量趋势分离(geodetic detrending)、精密单点定位、无电离层组合、周跳探测与修复等技术剔除载波相位观测值中的几何距离以及固体潮、海潮、大气潮、极潮、天线相位中心、相对论、相位缠绕、接收机钟差和对流层延迟误差，获得载波相位观测值的残差；然后利用离散小波变换技术降低噪声对电离层闪烁信号进行提取；最后利用连续小波变换技术对去噪后的信号进行时频分析，对经验频带内的小波系数进行逆变换，提取出电离层闪烁信号，并设定滑动窗口，对窗口内的电离层闪烁信号取标准差，从而完成每个频点上电离层相位闪烁因子的构建。本文将基于CHAIN提供的大量观测数据，通过对比分析提出的电离层相位闪烁因子、ROTI与ISMR提供的相位闪烁因子的相关性和准确性，验证其有效性。

1  闪烁因子的构建方法

本文基于1 Hz GNSS观测数据，致力于构建可以提供每个频点的电离层相位闪烁因子，详细过程如下。

1.1  大地测量趋势分离

大地测量趋势分离的概念为利用公式、模型及公开参数等手段剔除GNSS载波相位中较为容易估计的误差^[16]。为提取GNSS载波相位观测值中的电离层闪烁误差，本文采用了如下方法进行大地测量趋势分离。

(1) 站星几何距离改正。为改正每颗卫星至测站之间的几何距离，需准确确定测站坐标和每个观测历元的卫星坐标。测站坐标利用CSRS- PPP在线精密单点定位软件通过静态解算方法获得。卫星坐标通过国际GNSS服务组织(IGS)提供的精密星历文件并采用二阶拉格朗日内插算法获得。根据测站坐标和每个观测历元的卫星坐标，通过欧几里得度量获得站星几何距离改正参数。

(2) 固体潮误差采用二阶简化潮汐模型进行改正，海潮、大气潮和极潮误差均采用国际地球自转服务(IERS)推荐的惯用模型改正^[18]。

(3) 测站的接收机天线相位中心(包括接收机天线参考点、接收机天线平均相位中心和接收机天线瞬时相位中心)和卫星天线相位中心(包括卫星天线平均相位中心和卫星天线瞬时相位中心)利用IGS提供的天线相位中心模型(即igs14.atx)改正。

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

(5) 相对论效应和相位缠绕误差采用文献[19]所述相关方法予以改正。

(6) 对流层延迟分为干分量和湿分量两部分，其中干分量采用文献[19]所述相关方法予以改正。由于难以利用模型对湿分量进行精密估计，本文采用如下较为简单的方法估计：将湿分量的初值设为0.1 m，余下历元的湿分量作为随机游走参数进行估计，对于未改正的对流层延迟湿分量作为精密单点定位中的参数进行解算。

GNSS载波相位观测值经过以上大地测量趋势分离后，每个频点信号(f)的残差可以表示为

(1)

式中，c为光速；M_w为对流层湿延迟投影函数；E为卫星高度角；B_f为硬件延迟误差；λ_f为波长；N_f为整周模糊度。由式(1)可以看出，大地测量趋势分离可以消除GNSS载波相位观测值的大部分误差的影响，但为提取出电离层闪烁误差(I_f^d)，仍需要消除或削弱载波相位观测值中未能通过大地测量趋势分离消除的接收机钟差(δt_R)，未通过模型改正的对流层误差(δρ_{z, wc})、对流层折射作用产生的误差(I_f^r)、周跳和观测噪声(∈_f)。

1.2  对流层湿延迟未能通过模型改正部分和测站接收机钟差的改正

对流层湿延迟和测站接收机钟差无法通过大地测量趋势分离技术精确改正。本文利用静态精密单点定位(PPP)技术，将对流层湿延迟(未通过建模改正部分)和测站接收机钟差作为精密单点定位中的待求参数，对二者进行进一步的改正，具体如下。

对式(1)中的GNSS原始观测信号经过大地测量趋势分离后的残差组成双频无电离层组合L_IF，将该组合作为精密单点定位的基础观测值，组成如下观测方程

(2)

(3)

式中，V_IF为残差；A为系数矩阵；D_IF为无电离层组合模型改正值；x为待解算参数包括测站坐标(δX，δY，δZ)、天顶对流层湿延迟(δρ_{z, wc})、接收机钟差(cδt_R)和同观测历元的m颗可视卫星的模糊度(N)。该观测方程采用卡尔曼滤波静态解算方式进行求解，滤波中的状态转移矩阵(Φ)为

(4)

系统噪声向量的协因数矩阵(Q)为

(5)

为提高天顶对流层湿延迟和测站接收机钟差的估计精度，本文采用正向和反向运算相结合的方式，并将反向运算的结果作为未能利用模型改正的天顶对流层湿延迟(δρ_{z, wc})和测站接收机钟差(c·δt_R)的估计结果。值得指出的是在估计接收机钟差和对流层湿延迟时，也可采用非组合精密单点定位技术^[20]。利用Niell湿延迟投影函数，可将天顶对流层延迟(δρ_{z, wc})换算到每一卫星斜路径上，进而改正该卫星观测值中的对流层延迟误差。如果接收机钟比较稳定，不存在钟跳时，利用PPP方法获得的接收机钟差可直接用来改正观测值中的钟差；但对于大多数测地型接收机，接收机厂商会在钟差漂移到某一阈值时，通过对其中插入钟跳的方式，使得接收机内部时钟与卫星钟同步精度控制在一定范围^[21]。此时利用PPP方法无法准确估计出接收机钟差变化的细节部分，需要对钟差进行进一步精细估计。

1.3  测站接收机钟中存在钟跳时的钟差精细估计及周跳探测与修复

GNSS载波相位观测值经过上文所提出方法修正后，其残差为

(6)

式中，δt_{R, cj}表示未能经过PPP方法修正的钟差。为提取该残差中的钟差，本文利用两个频点的残差组成无电离层组合观测值，可表示为

(7)

该组合可以较为完好地消除电离层折射作用的影响，但无法消除电离层衍射与散射作用(即电离层闪烁)和模糊度的影响。由于模糊度在同一观测弧段内为同一值，因此可以通过历元间做差的方法消除。忽略观测噪声的影响，历元间做差之后的残差可表示如下

(8)

式中，Δ表示历元间做差运算符；ΔI_IF^d表示电离层闪烁随时间的变化率，其量值通常较小，本文将其忽略；ΔN_IF表示周跳。因此在未发生周跳且忽略电离层闪烁时间变化率的假设下，上式表示的历元间做差之后的残差即为接收机钟差在时间上的变化率。

由式(8)可以看出，周跳会影响钟差的估计精度，故而在估计钟差之前需要进行周跳探测与修复。本文首先利用式(6)给出的两个频点的残差组成双频宽项组合和无几何关系组合联合进行初次周跳探测，该方法可以修复大多数的周跳，但通常情况下，尤其是在电离层较为活跃的条件下，难以准确修复全部周跳。为此需要对未准确修复的周跳进行进一步探测和修复。本文将初步修正周跳后的每个频点的残差组成的无电离层组合观测值并进行历元间差分，由于该组合观测值可以极大地降低电离层的影响，从而提供精度较高的检测值，能够较为准确的检测出小周跳的存在。但只采用这一组合，无法确定周跳的大小，导致无法修复探测到的周跳。为此本文将探测发生周跳的历元作为新弧段的开始，从而在该观测弧段内可以为钟差的估计提供无周跳影响的历元间差分无电离层组合观测值。需要指出的是单独用无电离层组合观测值进行周跳探测会存在不敏感周跳组问题^[22]，仍然难以探测全部周跳，但因为本文采用了小波变换技术，少量周跳不会影响最终闪烁因子的估计。

为避免多路径效应对估计接收机钟差的影响，本文设置卫星的截止高度角为20°。为综合利用每个历元内截止高度角以上的所有卫星观测值，本文采用卫星高度角定权的方法对各卫星的历元间差分无电离层组合观测值进行加权平均，其权重(P)为

(9)

加权平均后的历元间差分后的钟差可表示为

(10)

式中，m表示历元k中卫星的数量；s表示卫星的序号。历元k时刻的钟差可以表示为从初始历元至历元k的历元间差分后钟差的数值积分，即

(11)

式中，n表示在初始历元到历元k的积分变量；dn为积分步长，通常等于观测的历元间隔。本文中初始的历元钟差设为0 s。

1.4  离散小波变换去噪

GNSS载波相位观测值经过上文修正后，其残差可表示为

(12)

式中，影响整周模糊度N_f在同一弧段同一性的周跳误差已经在1.3中修正过。本文采用离散小波变换方法削弱观测噪声(ε_f)对提取电离层闪烁信息的干扰。式(12)中I_f^r+I_f^d+B_f+λ_fN_f是目标信号，其在时间域具有连续性，在频率域内为低频信号，故而在小波域，该信号产生的小波系数的模值往往较大；ε_f通常为高斯白噪声，其在时间域上不具有连续性，在频率域内为高频信号，因而其在小波域的小波系数的模值通常较小，因此可以利用离散小波变换的方法在小波域内区分出目标信号和噪声。

利用离散小波变换去噪的过程包括将式(12)提供的残差数据转换到小波域，并对小波系数进行阈值处理，然后将小波系数进行逆变换获得目标信号。具体如下：

(1) 分解。选用合适的小波基对数据进行多级分解，获得小波系数。去噪的小波基需要正交或双正交，本文选用较为常用的多贝西极限正交相位小波，消失矩设为2。分解层数(level)的确定方法如下

(13)

式中，[·]表示向下取整运算；n为数据长度。

(2) 小波系数阈值处理。阈值分为软阈值和硬阈值。相比于软阈值，硬阈值可以保持闪烁的幅值，故而本文按照硬阈值的原则对每一层的小波系数进行处理，即所有小于阈值的小波系数均赋值0，大于阈值的小波系数保持不变。本文采用的降噪阈值确定方法为考虑信噪比的Stein无偏似然估计和固定阈值估计折中的原则，即信噪比很小时，按照Stein无偏似然估计处理，而信号噪声很大时，采用固定阈值处理，其计算方法为，其中n为数据的长度。

(3) 重构。利用阈值处理后的多级小波系数进行小波重构，获取目标信号。

1.5  电离层闪烁信号提取及闪烁因子构建

GNSS载波相位观测值经过1.4所提出方法进一步修正后，观测噪声已经被大大削弱了，其残差可以表示为

(14)

本小节致力于将电离层闪烁信号I_f^d从电离层折射信号I_f^r中分离出来，基本思路如下：因为电离层折射信号和闪烁信号的频率不同，二者混杂后组成的混合信号不是稳态信号，故而可以利用小波变换的方法将信号从时间域变换到时间-频率域，通过时频分析并对特征频率区间的小波系数进行小波逆变换从而获得电离层闪烁信号，进而利用以下公式获得基于1 Hz数据的电离层相位闪烁因子(σ_{ϕf, wavelet})

(15)

式中，＜·＞表示以一定时间间隔内的期望，本文选用60 s为滑动窗口。值得指出的是，虽然本文借用了文献[16]提出的大地测量趋势分离技术和钟差精细估计方法，但本文在观测值定权和电离层闪烁信息提取方法上具有显著的不同。

2  数据简介

为确定构建所提出闪烁因子的相关经验参数及验证其对电离层闪烁的探测效果，本文选取了加拿大高纬度北极电离层闪烁监测网络(CHAIN)中的11测站^[23]，它们的分布如图 1所示，均位于极隙区至极光区内，其中chuc测站具有并址磁测站，可以提供时间分辨率为1 min的地磁场强度信息。每个测站均配备Septentrio PolaRxS Pro电离层闪烁监测接收机和Septentrio的PolaN GG天线。CHAIN中的接收机被设定为接收双频GPS信号，且可以利用采样频率为50 Hz的观测数据直接产生电离层相位闪烁因子σ_ϕ，本文将其视为评价1 Hz数据闪烁因子的参考值。CHAIN同时可以提供1 Hz的GPS观测数据，本文利用其来研究所提出电离层相位闪烁因子的有效性。

图 1 本文选用的11个测站的分布

Fig. 1     Distribution of the 11 stations selected in this paper      

图选项

本文选择以上11个测站于2009年DOY 230至DOY 260、2013年DOY 120至DOY 150、2015年DOY 90至DOY 110、2017年DOY 100至DOY 130、2018年DOY 220至DOY 250和2020年DOY 80至DOY 112，共6组188 d的观测数据作为研究对象，每组数据中均包括强弱不同的地磁活动，利用每组数据获得的结论相似，为避免重复，这里重点介绍最后一组的地磁强弱情况。电离层闪烁与地磁活动相关^[24]，因此本文选择Kp、Dst和ASY 3个参数来反映的地磁活动情况(如图 2所示)。Kp指数表征全球地磁风暴的幅度，其中Kp≥5表示发生地磁风暴。Dst指数表示地球表面南北极与极赤道的轴向对称扰动磁场。Dst中主要扰动为负值，即地磁场的减少，主要产生于磁层中的电子环流；而Dst的正向扰动主要是太阳风作用在磁层造成的。ASY表征中纬度地区地磁扰动，其有两个分量，即水平分量(偶极子-极方向，H)和垂直分量(东-西方向，D)，其中ASY-H的变化趋势被证实与AE指数具有强相关性。AE指数通常被用来表征极光区电子活跃程度，但由于近期的AE指数无法获得，故本文用ASY指数近似替代AE指数来表征研究区域的电子活跃程度。由图 2可以看出，2020年DOY 111具有较强的地磁风暴，该日的观测数据将是本文重点研究对象。

图 2  2020年年积日80至112日的地磁活跃程度

Fig. 2  Degree of geomagnetic activity on DOY 80 to DOY 112 of 2020      

图选项

3  利用小波变换提取电离层闪烁的经验参数确定如上文所述，本文利用连续小波变换的方法将电离层闪烁信号与电离层折射信号相分离，而恰当的选取小波变换的小波基和小波变换相关参数会提高该分离效果，进而提高所构建相位闪烁因子的准确性。下面将对电离层闪烁时频分析的小波基及对称参数(symmetry parameter) 和时间带宽积(time-bandwidth product)的选取进行详细阐述。

3.1  小波基的选取

可进行时频分析的连续小波变换通常有3种小波基，即广义Morse小波(generalized morse wavelets)、解析Morlet小波(analytic morlet wavelet)和Bump小波。本文采用受电离层闪烁影响较为严重的观测弧段(arcc测站PRN07卫星2020年DOY 111 15:00至19:00的观测数据)，对上述3种小波的时频分析能力进行比较，方法为：首先以1 Hz为采样频率、以该弧段的观测历元数14 400为信号长度(其中Morse小波的对称参数和时间带宽分别设置为3和60)分别建立以上3种小波的小波滤波器库，并利用Matlab提供的频响特性函数(freqz)绘制幅频响应图(图 3)；然后利用建立的滤波器库对该弧段观测数据进行小波变换，绘制时频谱图(图 4(a)—图 4(c)，观察其对电离层闪烁信息的响应能力；接着以0.1 Hz至0.4 Hz为特征频率区间，将对应的小波系数进行逆变换，确定每种小波基对应的电离层相位闪烁因子(图 4(d)—图 4(f)，依据其与参考值(σ_ϕ)的契合程度，尤其是在发生电离层闪烁时段(σ_ϕ≥0.2 rad)的契合程度，选取恰当的小波基。

图 3   3种小波基的幅频响应图

Fig. 3    Amplitude-frequency response diagrams of the three wavelets      

图选项

图 4   arcc测站PRN07卫星2020年DOY 111 15时至19时观测弧段的时频谱图和相应闪烁因子

Fig. 4  The time-frequency spectrum and corresponding scintillation index based on the observation arc from 15h to 19h of PRN07 collected on DOY 111 of 2020 at arcc station      

图选项

从幅频响应图(图 3)可以看出，Bump小波的频率响应区间最小(如图 3(a)所示)，而电离层闪烁通常跨越的频率区间比较大，使得利用Bump提取的电离层闪烁信号强度会明显低于电离层闪烁接收机高频观测数据给出的相位闪烁因子(如图 4(a)所示)，导致最终确定的闪烁因子远低于参考值(如图 4(d)所示)。如图 3(b)和图 3(c)所示，Morlet小波和Morse小波在全频率域均有响应，但Morse小波在对高频信号响应更佳。Morlet小波和Morse小波均能较好地反映出电离层闪烁信号的时频信息(如图 4(b)和图 4(c)所示)，但利用Morlet小波变换后所得到的电离层相位闪烁因子的幅值明显大于参考值(如图 4(e)所示)，进而导致过多的误报，而利用Morse小波获得相位闪烁因子与参考值基本相当(如图 4(f)所示)，且Morse小波可以通过调整其对称参数和时间带宽积参数而进一步提高所构建闪烁因子的准确性，所以本文选择Morse小波作为连续小波变换的小波基。

3.2  对称参数和时间带宽积的选取

对称参数和时间带宽积是利用Morse小波基时频分析的重要参数。Morse小波的对称参数通过解调偏斜(demodulate skewness)控制着小波在时间上的对称性；时间带宽积与小波的持续时间成正比，决定时域小波的中心窗口在峰值频率处可以容纳振荡的个数。二者共同决定着小波的形状，进而影响小波变换的效果，因此恰当的选择小波参数对于电离层闪烁信息的准确提取至关重要。

为确定最优Morse小波参数，进而准确提取1 Hz GNSS观测值中的电离层闪烁信息，首先分别利用有电离层闪烁(arcc测站2020年DOY 111 PRN07卫星16:00至18:00)和无电离层闪烁(arcc测站2020年DOY 111 PRN02卫星00:00至03:30)弧段以及两组距离较远的参数组合作特例分析，以确定最优参数组合的所在区间。分析方法如下：利用对称参数3和6、时间带宽积60和20组成两个参数组合(3, 60)和(6, 20)，分别对以上两种弧段的观测数据进行连续小波变换，并对特征频率区间(0.1~0.4 Hz)的小波参数进行连续小波逆变换，构建相应的相位闪烁因子；以参考值为依据对不同参数组合构建的闪烁因子进行3方面评估：相关性、线性拟合的斜率和残差的均方根误差(RMS)。相关性越高，说明二者的相似程度越高；线性拟合的斜率越接近于1，残差的均方根误差越小说明构建的闪烁因子越接近参考值。

特例分析的结果如图 5所示。由图 5(a)—图 5(d)可以看出利用相同参数组合作用于闪烁强度不同的数据获得的闪烁指数与参考值的符合程度不同，利用不同参数作用于同一弧段数据所得闪烁指数也不相同。为评价这两种参数组合参数表现优劣，进一步对由小波变换构建的闪烁因子与参考值做线性拟合并进行相关性分析和残差的均方根误差分析，其结果如图 5(e)—图 5(l)所示。可以看出，在有电离层闪烁发生时段，本文所提出方法构建的闪烁因子与参考值的相关性更高，其拟合斜率更接近于1，但由于在强闪烁发生时，闪烁幅值较大，导致残差的均方根误差较大；而在闪烁强度小于0.05 rad时，由于受到未能剔除的观测噪声影响，估计的闪烁指数呈现出较强的随机性，进而降低了其与参考值的相关性。利用参数组合(3, 60)获得的闪烁指数的拟合斜率小于1而参数组合(6, 20)获得拟合斜率大于1，因此通过以上特例分析可以看出，小波参数的选择会影响所估计的闪烁指数，最佳的参数组合应该位于对称参数3~6、时间带宽积20~60。

图 5 利用arcc测站2020年DOY 111的数据对组合参数(3, 60)和(6, 20)特例分析的结果

Fig. 5     Special case analysis of two parameter groups (3, 60) and (6, 20) observation arcs collected on DOY 111 of 2020 at arcc station      

图选项

在大体确定最佳参数组合潜在范围的基础上，本文进一步利用全部11测站188 d的观测数据选取最佳的参数组合，方法为：分别以1和5为步长选取位于3至6区间和20至60区间的对称参数和时间带宽积，组成多种参数组合；基于以上数据，利用多种组合参数的Morse小波构建相应的相位闪烁因子，并比较其与参考值线性拟合的斜率、相关性和残差的均方根误差。图 6和图 7分别统计了由188 d的数据获得的上述3种指标在有闪烁发生时段和无闪烁发生时段的变化情况，并以误差线的形式展示每种指标在对应参数组合下波动的标准差。可以看出，以上分析可以确定获得较好结果的参数组合所在的大致范围。在有闪烁发生时段，对称系数为3时可以给出更为优良的结果，其与时间带宽积45、50、55和60组成的参数组合具有相似且较为优良的表现，均可以使线性拟合斜率趋近于1，相关性达到92%而残差的均方根误差降至0.03 rad；而在无闪烁发生时段，依旧是对称系数3表现更佳，其与时间带宽积30、35、40、45组成的参数组合可以获得较好的结果，但其相关性要略低于有闪烁发生时段，笔者认为这主要是由于观测噪声在无闪烁发生时段带来的影响导致估计的闪烁因子呈现一定的随机性导致。

图 6 有闪烁发生时段，所构建闪烁因子与参考值之间的统计分析

Fig. 6     Statistical analysis of the reference and the proposed scintillation index      

图选项

图 7 无闪烁发生时段，所构建闪烁因子与参考值之间的统计分析

Fig. 7     Statistical analysis of the reference and the proposed scintillation index      

图选项

为了从图 6和图 7给出的较优参数组合范围中确定最佳参数组合，本文进一步分析了每个测站的最优和次优参数组合，并计算了最优和次优参数组合在线性拟合斜率、相关性和残差的均方根误差的指标差，计算方法如下

(16)

(17)

式中，表示绝对值运算符；残差指标差的计算方法同相关性指标差。当线性拟合斜率指标差为负值、相关性指标为正值、残差指标差为负值时，说明最优参数组合的对应指标更优，反之说明次优参数组合的相对应指标更优。表 1和表 2给出了以上分析的结果。可以看出，参数组合的指标结果随测站的变化而不同，但从统计的角度，(3, 45)和(3, 35)分别在有闪烁和无闪烁发生时段被选为最优或次优参数组合的次数最多。较大的时间带宽积(≥45)更适用于有闪烁发生时段，而较小的时间带宽积(30至45)在无闪烁发生时段可以获得更好的指标。值得注意的是，尽管适用于不同测站的最优及次优组合不同，但各参数组合的指标差相差不大，均在标准差以内(见图 6和图 7)，说明表 1和表 2所统计的参数组合均可在一定条件下对电离层闪烁进行较为精确的估计。在实际应用中难以做到在不同的电离层闪烁条件下选择不同的参数组合进行电离层闪烁监测，为此本文选择(3, 45)作为最佳参数组合。需要注意的是，这里仅仅通过试算的方法选取出本文所考虑的几组参数组合中的较优参数组合，不是全部参数组合中的最佳选择，但试验发现当确定对称参数后，由相临近的时间带宽积给出的电离层闪烁因子幅值相差不大，因此本文不进行更小范围内的参数选取研究。

表 1 有闪烁发生时段，各测站最优与次优参数组合及其指标差

Tab. 1 The optimal and sub-optimal parameter groups and their differences on each station during periods with scintillation

表 2 无闪烁发生时段，各测站最优与次优参数组合及其指标差

Tab. 2 The optimal and sub-optimal parameter groups and their differences on each station during periods without scintillation

4  与ROTI指数的比较为验证所提出的闪烁因子σ_{ϕf, wavelet}的有效性，本文利用全部11测站6组共188 d的观测数据，以CHAIN网络高频采样数据输出的相位闪烁因子σ_ϕ为参考值，对所提出的闪烁因子和ROTI指数进行对比。由于每组数据获得的结果类似，为避免重复，这里重点介绍由第6组(即2020年DOY 80至DOY 112)观测数据获得的结果。ROTI是当前最为常用的基于低采样频率观测数据的闪烁因子^[25]，被诸多学者用来监测电离层闪烁^[26-32]。ROTI表示斜路径总电子含量(STEC)的变化率，其计算方法可参考文献[12]。ROTI中的STEC直接由双频载波无几何关系组合获得，因此ROTI容易受到周跳的影响。当前已有方法难以探测并修复全部周跳，尤其对于低卫星高度角或有闪烁发生的弧段，为此本文将截止高度角设为20°并采用文中所述方法，以最大限度减弱周跳对ROTI的影响。值得一提的是，本文所提出的闪烁因子构建方法中由于引入了小波变换技术，可以较为容易地克服未正确探测及修复的周跳。

4.1  电离层闪烁的监测效果

为验证所提出闪烁因子对电离层闪烁的监测效果，本文利用chuc测站2020年DOY 110和DOY 111的观测数据，对比分析了σ_{ϕf, wavelet}、σ_ϕ和ROTI在每小时时间内所探测出的电离层闪烁历元的数量。选用该测站的原因是因为chuc测站有并址磁测站，而电离层闪烁与本地的地磁场活动直接相关，因而可以利用磁测站提供的地磁场数据，验证所探测出受电离层闪烁影响历元数量的正确性。需要指出的是σ_{ϕf, wavelet}与ROTI的时间分辨率与GPS数据的采样频率相同(本文中为1 s)，而σ_ϕ的时间分辨率为60 s，为方便对比，本文将σ_{ϕf, wavelet}与ROTI按照σ_ϕ中对应的时刻进行了重采样。本文采用经验值的方式确定电离层闪烁的阈值，即σ_{ϕf, wavelet}与σ_ϕ为0.2 rad、ROTI为0.1 TECU/min^[33]。当闪烁因子的幅值大于阈值时即判定发生电离层闪烁。

图 8展示了chuc测站在第110和111天的地磁场变化情况及每小时受电离层闪烁影响的历元数量。由图 8(a)可以看出DOY 111的地磁活动要强于DOY 110。由图 8(b)和8(c)可以看出3种闪烁因子所探测出发生电离层闪烁的历元均分布在地磁场较为活跃的时段，这在一定程度上再次证明了电离层闪烁与本地的地磁场活动直接相关的结论以及3种闪烁因子对电离层闪烁监测的有效性；σ_{ϕf, wavelet}与ROTI所探测出的发生电离层闪烁的历元数量与参考值所给出的数量基本相同(尤其是L2载波)，说明σ_{ϕf, wavelet}与ROTI均可以较为准确地监测电离层闪烁；尽管ROTI无法给出每个频点上的闪烁数量，但是其探测结果与σ_{ϕf, wavelet}在L1、L2频点探测结果均相同，这主要是由于在高纬度地区，电离层闪烁对GPS的L1、L2频点的影响几乎相同^[34]。综上所述，本文所提出的闪烁因子σ_{ϕf, wavelet}和ROTI因子均能准确地探测出电离层闪烁的发生，可以用于电离层闪烁的定性监测研究。

图 8 σ_{ϕf, wavelet}、ROTI及σ_ϕ对chuc测站电离层闪烁的探测结果对比

Fig. 8     Ionospheric scintillation detection results provided by σ_{ϕf, wavelet}、ROTI and σ_ϕ utilizing the observations from chuc station      

图选项

4.2  电离层闪烁的幅值估计准确度

由4.1中的试验可以证明σ_{ϕf, wavelet}和ROTI均可以较为准确地探测出电离层闪烁发生，但无法证明二者估计的电离层闪烁幅值的准确度。本文将相关性作为为衡量电离层闪烁因子估计的幅值准确度的指标。在σ_{ϕf, wavelet}和ROTI均能准确探测出电离层闪烁发生的前提下，二者与σ_ϕ的相关性越高，说明其对电离层闪烁估计的幅值越接近σ_ϕ。本文将全部11个测站188 d的观测数据分成有闪烁发生(σ_ϕ≥0.2 rad)和无闪烁发生(σ_ϕ＜0.2 rad)两种类别，分别研究由每颗卫星获得的σ_{ϕf, wavelet}、ROTI与σ_ϕ的相关性。

图 9和10分别展示了利用2020年DOY 80至DOY 112数据获得的有闪烁和无闪烁发生时段的相关性结果，其中，空白区域表示该卫星不可用或无闪烁发生。因在相同纬度附近测站的结果类似和篇幅限制，本文仅展示了arcc、sarc、repc和chuc这4个测站的结果。由图 9可以看出，σ_{ϕf, wavelet}与σ_ϕ的相关性在全部有闪烁发生的弧段均在0.8以上，且高于ROTI与σ_ϕ的相关性，说明在有闪烁发生时，σ_{ϕf, wavelet}的幅值更加接近σ_ϕ。由图 10可以看出，在无闪烁发生的时段，纬度较高的测站(arcc和sarc)所给出的σ_{ϕf, wavelet}、ROTI与σ_ϕ的相关性均要高于纬度较低的测站，这主要是因为arcc和sarc测站位于极隙区，而其余两个测站分布在极光区，极隙区的磁场线呈现开放状态，导致该区域的电离层更容易受到太阳风的影响而变得活跃，进而导致电离层闪烁的发生频率更高。尽管此时极隙区和极光区闪烁均没有超过阈值，但是极隙区的闪烁幅度要高于极光区的闪烁，使得由1 Hz采样频率观测数据估算出来的闪烁因子受观测噪声的影响更小，进而提高了σ_{ϕf, wavelet}、ROTI与σ_ϕ的相关性。图 10也可以看出，对于极隙区测站(arcc和sarc)，σ_{ϕf, wavelet}的相关性明显优于ROTI，但该图难以反映二者的相关性在其他区域测站(如repc、chuc)的情况，需进一步统计分析。

图 9 有闪烁发生时段σ_{ϕf, wavelet}、ROTI与σ_ϕ(参考值)的相关性比较

Fig. 9     Comparison of the correlation provided by σ_{ϕf, wavelet}、ROTI with σ_ϕ (the reference)      

图选项

图 10 无闪烁发生时段，σ_{ϕf, wavelet}、ROTI与σ_ϕ(参考值)的相关性比较

Fig. 10     Comparison of the correlation provided by σ_{ϕf, wavelet}、ROTI with σ_ϕ (the reference)      

图选项

为进一步比较σ_{ϕf, wavelet}与ROTI在无闪烁发生时段的电离层闪烁的幅值估计准确度，本文对全部的σ_{ϕf, wavelet}、ROTI与σ_ϕ的相关性按照下式进行如下统计分析

(18)

式中，corr表示相关性；s表示测站卫星；NS表示测站在一天内全部可视卫星数量；DOY表示年积日；P表示σ_{ϕf, wavelet}与σ_ϕ的相关性大于ROTI与σ_ϕ的相关性的概率。式(18)表示的分析方法为：将每个测站每天数据获得的σ_{ϕf, wavelet}与σ_ϕ的相关性减去ROTI与σ_ϕ的相关性，找出差值大于0的卫星的数量，并除以全部的卫星数量，若最终的结果大于50%，则说明在该测站该天，σ_{ϕf, wavelet}比ROTI有更大的可能性获得相关性更高的闪烁因子。全部11个测站于2020年DOY 80至DOY 112观测数据的统计结果如图 11所示，相比于ROTI，σ_{ϕf, wavelet}有72.9%的概率获得与σ_ϕ相关性更高的闪烁因子，即给出的闪烁幅值更为准确。

图 11 σ_{ϕf, wavelet}与σ_ϕ的相关性大于ROTI与σ_ϕ的相关性概率统计

Fig. 11     Statistics of the probability that the correlation provided by σ_{ϕf, wavelet} and σ_ϕ is greater than that provided by ROTI and σ_ϕ

图选项

5  结论

本文以1 Hz采样频率GNSS载波相位观测值为基础，通过大地测量趋势分离、精密单点定位、无电离层组合估计接收机钟差、周跳探测与修复、离散小波去噪等技术，剔除了载波相位观测值中的卫地几何距离以及固体潮、海潮、大气潮、极潮、天线相位中心、相对论、相位缠绕、对流层延迟、接收机钟差、周跳和观测噪声误差，并详细研究了利用连续小波变换技术提取电离层闪烁信号中的最优小波基及相关参数的经验值，进而构建了基于1 Hz观测数据的电离层相位闪烁因子。为了验证所提出闪烁因子的有效性和准确性，本文以高采样频率数据的相位闪烁因子为参考值，采用加拿大高纬度北极电离层闪烁监测网络(CHAIN)中11个测站188 d的1 Hz观测值对所提出闪烁因子的有效性进行验证，并与已经广泛使用的ROTI因子进行对比，结果表明在电离层闪烁发生时段，本文所提出的闪烁因子与ROTI均能有效地探测到电离层闪烁的发生，但本文提出的闪烁因子所给出的电离层闪烁的幅值更接近参考值，准确度更高，尤其在有强电离层闪烁发生的时段。

文献[35]指出接收机相位偏差存在受温度影响的短期变化，而本文为了简便，认为相位偏差在同一观测弧段不变化，未来的研究若能考虑到短期的温度变化，或可在一定程度上提高本文所构建相位闪烁因子的对电离层闪烁的监测精度。需要指出的是，构建该闪烁因子的原理也可用于测地型接收机的观测值，由于篇幅限制，笔者将在未来通过试验验证其可靠性，以拓展该闪烁因子的应用场景。

作者简介

第一作者简介：赵东升(1992—), 男, 博士, 讲师, 研究方向为GNSS定位、GNSS遥感、电离层闪烁监测与预报。E-mail: [e-mail]

通信作者：张克非. E-mail: [e-mail]

初审：张艳玲复审：宋启凡
终审：金   君