赵东升，李 旺，李宸栋，唐 旭，张克非

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

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

为了测量电离层闪烁对GNSS信号的影响，通常需要电离层闪烁监测接收机(ionospheric scintillation monitoring receivers，ISMR)，它可以直接给出两类闪烁因子：振幅闪烁因子(S4)和相位闪烁因子(σφ)。由于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，余下历元的湿分量作为随机游走参数进行估计，对于未改正的对流层延迟湿分量作为精密单点定位中的参数进行解算。

(1)

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

VIF=Ax-(LIF-DIF)

(2)

(3)

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

Φ=diag[111011…1]

(4)

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

Q=diag[0009×101010-40…0]

(5)

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

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

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

(6)

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

(7)

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

(8)

为避免误读、误会，释义，便成了首要功课。“三杂”之杂，为多种多样、丰富多彩意，决非杂乱之杂……大凡杂文家，一定会有相应的知识储备、文体实践和人生阅历——如此之“三”以“三杂”谓之，不算离奇、搞怪吧？

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

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

P=sin2E

(9)

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

(10)

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

(11)

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

1.4 离散小波变换去噪

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

(12)

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

level=max(5,[log2n])

(13)

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

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

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

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

(14)

(15)

2 数据简介

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

本文选择以上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具有较强的地磁风暴，该日的观测数据将是本文重点研究对象。

3 利用小波变换提取电离层闪烁的经验参数确定

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

图1 本文选用的11个测站的分布Fig.1 Distribution of the 11 stations selected in this paper

图2 2020年年积日80至112日的地磁活跃程度Fig.2 Degree of geomagnetic activity on DOY 80 to DOY 112 of 2020

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)可以看出，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 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.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。

在大体确定最佳参数组合潜在范围的基础上，本文进一步利用全部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和图7给出的较优参数组合范围中确定最佳参数组合，本文进一步分析了每个测站的最优和次优参数组合，并计算了最优和次优参数组合在线性拟合斜率、相关性和残差的均方根误差的指标差，计算方法如下

线性拟合斜率指标差=

(16)

相关性指标差=最优参数组合相关性-

次优参数组合相关性

(17)

图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

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的影响。值得一提的是，本文所提出的闪烁因子构建方法中由于引入了小波变换技术，可以较为容易地克服未正确探测及修复的周跳。

图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

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

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

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因子均能准确地探测出电离层闪烁的发生，可以用于电离层闪烁的定性监测研究。

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)的情况，需进一步统计分析。

图8 σφf,wavelet、ROTI及σφ对chuc测站电离层闪烁的探测结果对比Fig.8 Ionospheric scintillation detection results provided by σφf,wavelet、 ROTI and σφ utilizing the observations from chuc station

为进一步比较σφ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%的概率获得与σφ相关性更高的闪烁因子，即给出的闪烁幅值更为准确。

图9 有闪烁发生时段σφf,wavelet、ROTI与σφ(参考值)的相关性比较Fig.9 Comparison of the correlation provided by σφf,wavelet、 ROTI with σφ (the reference)

5 结 论

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

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

图10 无闪烁发生时段，σφf,wavelet、ROTI与σφ(参考值)的相关性比较Fig.10 Comparison of the correlation provided by σφf,wavelet、ROTI with σφ (the reference)

图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 σφ