Maximum likelihood estimation of irregularity parameters with spectrum fitting
-
摘要:
获得闪烁相关电离层不均匀体参数是建立闪烁预测模型的基础. 针对长期历史数据为单站观测的情况，提出利用闪烁谱拟合方法，通过最大似然估计（maximum likelihood estimation, MLE）获得多个闪烁相关参数的最优估计，并利用Monte Carlo方法对谱参数MLE结果进行分析. 结果表明：利用单指数谱模型可以获得谱参数的无偏最优估计，估计参数的90%置信度范围可以覆盖真实值；利用双指数谱模型可以获得谱参数的近似无偏最优估计，绝大部分估计参数的90%置信度范围可以覆盖真实值；两种模型估计的统计结果均符合预定分布. 利用实测数据对本文方法进行进一步验证，估计参数的统计分析结果符合预期分布，验证了该方法的有效性；仿真数据和实测数据分析结果同时还表明，该方法可适用于中等及强闪烁情况. 本文方法为利用长期历史数据建立闪烁预测模型提供了重要实现途径.
-
关键词:
- 电离层闪烁 /
- 功率谱 /
- 最大似然估计(MLE) /
- 全球导航卫星系统
Abstract:Parameters on ionospheric irregularity are essential for scintillation modeling. As scintillation is generally observed on a single site basis, it had been a hard work to retrieve various irregularity parameters. However, development of scintillation theory model provides a possible solution to this difficulty. Here we present a new method to estimate irregularity parameters with maximum likelihood estimation (MLE) based on scintillation spectrum fitting. Quality of the method is verified with Monte Carlo simulation. It shows that an unbiased optimal estimation could be achieved for single spectral index model with all results located within the 90% confidence interval, while a near unbiased optimization could be reached for double spectral index model with most of results in the 90% confidence interval. Statistic distribution of estimation results follows the theoretical one well, showing the viability of estimation. Test is also conducted with real measurements and the result is promising. Two kind of data tests show the method could work efficiently under modest and strong scintillation conditions. This would pave the way for scintillation modeling with commonly observation on single site.
-
0 引 言
电离层闪烁是影响GNSS及其应用的重要空间环境因素，对其进行监测预测是应对闪烁影响的有效措施. 闪烁影响不仅与造成闪烁的小尺度电离层电子密度不均匀体特征有关（包括扰动强度、谱指数、内外尺度等），还与用户的观测几何有关[1-2]. 直接利用用户测量结果如闪烁指数S4、电子总含量（total electron content, TEC）变化率指数（rate of TEC index, ROTI）等进行闪烁预测建模，会使模型变得复杂且预测性能有限[3-4]. 因此，已有的闪烁模型均采用对不均匀体特征进行建模，并通过信号传播模型获得接收机处信号变化特征. 如WBMOD模型中对不均匀体强度、谱指数、漂移速度等参数进行建模；全球电离层闪烁模型（Global Ionospheric Scintillation Model，GISM）通过对电子密度起伏变化建模[5-6]，并通过相位屏传播模型模拟产生闪烁信号[7].
利用闪烁观测数据对扰动强度、谱指数、漂移速度等不均匀体参数进行有效估计是建立闪烁预测模型的基础. 通常利用卫星信号对闪烁进行观测，随着GNSS应用的发展，利用GNSS信号进行电离层闪烁监测已成为一种重要手段. 对卫星信号进行高频采样观测并做谱分析，可以获得闪烁信号频谱. 理论分析表明，造成闪烁的电离层不均匀体空间谱特征，与接收的卫星闪烁信号的频谱特征具有相关性[8-9]. 因此可以通过对卫星闪烁信号频谱特征的深入分析，获得电离层不均匀体参数的估计值.
利用谱拟合方法对空间环境参数进行估计是行星际闪烁、X射线等天文学研究中的一种常用方法[10-12]. 但GNSS信号电离层闪烁测量须要建立合适的卫星闪烁信号功率谱模型，才能利用卫星信号谱分析结果进行闪烁参数估计.
闪烁理论表明，弱散射情况下卫星信号的相位谱和功率谱在高频端表现出幂率谱下降特征. C.Rino基于相位屏理论，建立了弱散射情况下的信号谱模型，给出卫星信号谱的表达式，并指出在高频端存在幂率谱近似Tscin/Tscinfpfp，其中Tscin为1 Hz处的相位谱强度，p为谱指数[13]. 基于此模型，一些学者通过对GPS卫星闪烁信号功率谱高频端一定频率范围内的谱系数进行拟合，实现对谱强度和谱指数的估计[9,14]. 但该方法在实际应用中受到很多限制. 为保证信号功率谱在高频端呈现出幂率谱变化特征，须选择一定强度的闪烁信号，并须要合理选择拟合的频率范围. 实际观测表明，卫星信号闪烁功率谱在一定强度闪烁下会出现扩展或增强[15]，在信号高频端很难看到明显的幂率谱变化特征. 同时，不均匀体的漂移速度变化、卫星自身运动，以及接收机观测噪声影响，难以确定一个固定的频率范围进行拟合[9,14]. 为此，C.S.Carrano等人[15-16]进一步利用C.Rino建立的强扰动下闪烁信号功率谱模型[17]进行直接拟合，对谱强度、谱指数等参数进行最优估计. 该方法的优势是选择最小截止频率和接收机高频观测噪声之间一段相对确定的频谱范围进行拟合，提高了谱拟合进行参数估计的可用性，但该方法缺乏对估计结果的合理评价. 后面分析可以看出，准确的谱模型能够保证参数估计结果符合特定分布，从而实现对结果置信度的估计. 近年来C.Rino等人针对GNSS应用提出了双指数闪烁谱模型，并通过数值进行模拟发现，双指数谱模型可以表现出不同闪烁情况下的信号功率谱变化特征[15,18]；进一步又利用双指数谱模型对天基测量数据进行了分析，验证了模型的有效性[19].
已有的长期闪烁观测数据主要为单站测量的卫星信标数据，如GPS信号、通信卫星信号等，因此本文利用闪烁功率谱拟合方法对闪烁相关参数进行最优估计，并利用仿真数据和实测数据对该方法的性能进行分析评估，统计分析结果表明了利用谱拟合方法估计闪烁参数的有效性.
1 电离层闪烁谱模型
C. Rino等人提出的双指数电离层闪烁谱模型已广泛应用到GNSS闪烁影响分析中[15,20-21]，双指数谱模型中定义了如下不均匀体相位谱密度函数模型[18]：
Φδϕ(q)=C′p{|q|−p1,q⩽ (1) 式中：C_p'为相位谱强度；q为空间波数；{q_0}为中断波数；{p_1}和{p_2}分别为低频部分和高频部分的相位谱指数. 双指数谱模型通过引入中断波数和两个谱指数，提高了闪烁谱分析的灵活性. 当两个谱指数相等时，双指数谱模型褪变为已有的单指数谱形式，此时中断波数可取任意值而不再具有实际意义.
利用菲涅尔半径\;{\rho _{\text{F}}}对空间波数进行归一化处理，可以获得更为统一的相位谱密度形式[18]
{{P}}\left( \mu \right) = \left\{ {\begin{array}{*{20}{c}} {{U_1}{\mu ^{ - {p_1}}},}&{\mu \leqslant {\mu _0}} \\ {{U_2}{\mu ^{ - {p_2}}},}&{\mu > {\mu _0}} \end{array}} \right. (2) 式中：谱强度参数为
{U_1} = C_p'\rho _{\text{F}}^{{p_1} - 1}\text{，}{U_2} = C_p'q_0^{{p_2} - {p_1}}\rho _{\text{F}}^{{p_2} - 1} (3) 两个谱强度参数满足{U_2} = {U_1}\mu _0^{{p_2} - {p_1}}；\;\mu \;{\text{ = }}\;q{\rho _{\text{F}}}为归一化波数；\;{\mu _0}\;{\text{ = }}\;{q_0}{\rho _{\text{F}}}为归一化中断波数；菲涅尔半径 \;{\rho _{\text{F}}} = \sqrt {{{\textit{z}} \mathord{\left/ {\vphantom {z k}} \right. } k}} ，{\textit{z}}为信号从不均匀体相位屏传播至接收机的距离，k为自由空间波数.
利用归一化的相位谱密度函数，接收信号强度的功率谱可以表示为[15,17-18]
I\left( \mu \right) = 2\int_0^\infty {\exp \left( { - \gamma \left( {\eta ,\mu } \right)} \right)\cos \left( {\eta \mu } \right){\text{d}}\eta } (4) 式中：\eta \;{\text{ = }}\;{r \mathord{\left/ {\vphantom {r {{\rho _{\text{F}}}}}} \right. } {{\rho _{\text{F}}}}}为利用菲涅尔半径\;{\rho _{\text{F}}}对空间间隔r进行归一化的结果；\gamma 为结构交互函数，定义为[18]
\begin{split} \gamma \left( {\eta ,\mu } \right) =& 16\int_0^{{\mu _0}} {{U_1}{\chi ^{ - {p_1}}}{{\sin }^2}\left( {\chi \eta /2} \right){{\sin }^2}\left( {\chi \mu /2} \right)\frac{{{\text{d}}\chi }}{{2{\text{π }}}}} {\text{ + }} \\ &{\text{ }}16\int_{{\mu _0}}^\infty {{U_2}{\chi ^{ - {p_2}}}{{\sin }^2}\left( {\chi \eta /2} \right){{\sin }^2}\left( {\chi \mu /2} \right)\frac{{{\text{d}}\chi }}{{2{\text{π }}}}} , \end{split} \chi 为积分变量.
实际观测中获得的是闪烁影响下信号起伏变化的时间序列，而时间序列变化由卫星信号穿越电离层不均匀体的相对运动引起. 因此，利用有效扫描速度{v_{{\text{eff}}}}可以进行空域谱和时域谱参数之间的转换，即\;\mu \;{\text{ = }}{{2{\text{π}}f} \mathord{\left/ {\vphantom {{2{\text{π }}f} {{f_{\text{F}}}}}} \right. } {{f_{\text{F}}}}}，其中菲涅尔频率{f_{\text{F}}} = {{{v_{{\text{eff}}}}} \mathord{\left/ {\vphantom {{{v_{{\text{eff}}}}} {{\rho _{\text{F}}}}}} \right. } {{\rho _{\text{F}}}}}. 进而获得闪烁影响下信号强度变化的时域谱形式为[21-22]
\begin{split} & I\left( {f;U,{p_1},{p_2},{\mu _0},{f_{\text{F}}}} \right) = \\ & 2\int_0^\infty {\exp \left( { - \gamma \left( {\eta ,\frac{{2{\text{π}}f}}{{{f_{\text{F}}}}}} \right)} \right)\cos \left( {\frac{{2{\text{π}}f\eta }}{{{f_{\text{F}}}}}} \right)} {\text{d}}\eta \end{split} (5) 由闪烁理论进一步可知，在已知信号强度功率谱情况下，存在下述关系式[15,18]：
S _4^2 = \frac{1}{{2{\text{π }}}}\int_{ - \infty }^\infty {I\left( \mu \right)} {\text{d}}\mu - 1 (6) 2 闪烁谱参数的最大似然估计(MLE)
利用闪烁影响下的信号强度时域谱模型，采用最小二乘估计和最大似然估计(maximum likelihood estimation, MLE)对实测信号强度谱进行拟合[16,22]，实现对信号强度谱模型中不均匀体参数的估计. 在已知信号概率分布情况下，利用MLE可以获得对信号参数的无偏估计.
定义\theta \;{\text{ = }}\left( {U,{p_1},{p_2},{\mu _0},{f_{\text{F}}}} \right)为信号强度谱中的参数. 对于时域频谱中的第i个频率分量{f_i}，假设实测数据的谱分析结果为I_i^m = {I^m}\left( {{f_i}} \right)，信号强度谱模型的预测结果为{I_i} = I\left( {{f_i};\theta } \right)，则二者的比值 {R_i} = {I^m}\left( {{f_i}} \right)/I\left( {{f_i};\theta } \right) 为一个随机变量. 从谱分析理论可知，傅里叶变换后某一特定频率处的频谱分量估计符合正态分布[23]. 若信号强度谱模型是准确的，则 {R_i} 的分布为{R_i}\sim {{\chi _d^2} \mathord{\left/ {\vphantom {{\chi _d^2} d}} \right. } d}，其中参数d\;{\text{ = }}\;2M，M为{I^m}\left( {{f_i}} \right)估计中的数据分段数. 对上式进行变换， d{R_i} 符合自由度为d的卡方分布，即d{R_i}\sim \chi _d^2[12,23-24].
对于功率谱的第i个频率分量{f_i}，给定功率谱模型估计结果{I_i}情况下，实际信号功率谱测量分量的概率密度函数为
p\left( {I_i^m|{I_i}} \right)\sim \frac{d}{{{I_i}}}\chi _d^2\left( {d\frac{{I_i^m}}{{{I_i}}}} \right) (7) 假定功率谱各频率分量上的测量是统计独立的，则在给定功率谱模型预测结果时，实测信号功率谱的联合密度函数为
p\left( {{I^m}|I} \right)\sim \prod\limits_{i = 1}^N {\frac{d}{{{I_i}}}\chi _d^2\left( {d\frac{{I_i^m}}{{{I_i}}}} \right)} (8) 因此，信号强度谱参数的似然函数为L\left( {\theta |{I^m}} \right) = p\left( {{I^m}|I\left( \theta \right)} \right). 式(8)可进一步表示为对数似然函数形式：
{\text{ln}}\left( L \right)\sim \sum\limits_{i = 1}^N {\frac{d}{{{I_i}}}\chi _d^2\left( {d\frac{{I_i^m}}{{{I_i}}}} \right)} (9) 通过使上述对数似然函数最大化获得对谱参数的最优估计值\hat \theta ，即为谱参数MLE.
统计理论进一步表明，在已知MLE结果\hat \theta 情况下，固定部分参数仅对剩余部分参数进行估计，则对数似然估计函数的差值\Delta S = {\text{ln}}\left( {L\left( \theta \right)} \right) - {\text{ln}}\left( {L\left( {\hat \theta } \right)} \right)符合自由度为\nu 的卡方分布，参数\nu 为待估部分参数个数[24-25]. 利用该性质可对\theta 中某个参数估计结果的置信度进行分析.
3 仿真分析与试验验证
Monte Carlo方法是一种统计试验方法，通过构造和系统性能近似的概率模型，对系统的随机性能进行模拟和统计分析[26]. 由于谱强度、谱指数等参数难以通过实际观测获得准确值，无法对谱拟合最优估计结果进行评价，因此利用Monte Carlo方法对谱参数MLE方法的性能进行分析. 设定合适的参数（包括谱强度、谱指数、漂移速度和中断频率），采用第2节中的双参数谱模型，利用相位屏方法可以对闪烁影响下的信号变化进行仿真[20-21,27]；利用仿真获得的闪烁信号进行谱分析，得到闪烁谱；对闪烁谱进行拟合，并利用MLE获得对闪烁谱参数的最优估计；对上述过程进行多次实现（本文为100次）；对多次MLE结果与闪烁谱参数预设值（作为准确值）进行比较，分析估计方法的性能.
理论上可对谱模型中的所有参数\theta \;{\text{ = }} \left( {U,{p_1},{p_2},{\mu _0},{f_{\text{F}}}} \right)进行估计. 由于漂移速度可以利用多站组网测量实现[3-4]，因此菲涅尔频率{f_{\text{F}}}可以很容易获得. 对漂移速度估计的有效性分析将会单独著文论述，这里仅对谱强度U和谱指数{p_1},{p_2}估计结果的有效性进行分析，分单指数谱和双指数谱两种情况进行. 其中，单指数谱可通过设置{p_1}={p_2}来实现，此时\;{\mu _0}可取任意值.
3.1 单指数谱估计结果分析
单指数谱分析中，设置谱强度U = 0.6，谱指数{p_1} = {p_2} = 3，菲涅尔频率{f_{\text{F}}}{\text{ = }}2\;{\text{Hz}}，利用相位屏方法对闪烁信号进行仿真，并对仿真的闪烁信号进行谱参数MLE.
图1给出了一次利用相位屏方法仿真获得的GNSS闪烁影响下的信号时间变化序列，长度为5 min. 由该信号变化序列计算获得的闪烁指数为{S _4} = 0.59.
![]()
图 1 仿真的单指数谱GNSS闪烁信号时间序列Fig. 1 Simulated GNSS signal time series under scintillation with one-component power law spectrum图2给出了对图1中闪烁信号时间序列的谱分析和谱拟合结果. 其中，蓝色曲线为仿真数据的谱分析结果. 谱分析采用Welch方法实现，将仿真数据分成5个数据段，每段长度1 min，对每段数据分别进行谱分析并进行平均获得最终的谱拟合结果.
![]()
图 2 仿真的GNSS闪烁信号单指数功率谱及MLE理论功率谱曲线Fig. 2 One-component power spectrum of simulated GNSS scintillation signals and that of MLE theory在谱拟合中进一步考虑噪声模型可以提高MLE的稳定性[10-11,22]. 因此，在理论闪烁谱基础上增加了简单的白噪声模型来进行参数估计. 图2中的粉色曲线为增加噪声模型后获得的功率谱最优估计结果，红色曲线为不考虑噪声时的理论谱模型结果. 可以看出，通过MLE获得的理论谱模型可较好地对实测功率谱进行拟合.
利用MLE方法对该仿真数据的功率谱进行拟合获得的最优参数估计结果为U = 0.55，{p_2} = 3.01，与理论值相近. 由式(6)可知，对信号功率谱进行积分可以得到闪烁指数. 利用最优谱拟合曲线计算闪烁指数，得到{S _4} = 0.55，与实际仿真数据计算结果接近. 上述结果说明了MLE方法在单指数谱形式下估计的有效性.
实测功率谱与理论谱模型的比值 d{R_i} 符合卡方分布[12,22-25]. 图3给出了利用图1数据获得的功率谱和估计理论谱比值的统计分布. 可以看出， d{R_i} 统计分布与卡方分布变化一致. 统计学中利用Kolmogorov-Smirnov检验（K-S检验）确定观测数据的潜在概率分布是否与假设分布存在显著差异，K-S检验门限值为0.05，当观测数据的概率分布符合已知概率分布时，K-S检验统计量小于预定门限值[28]. 对 d{R_i} 分布进一步进行Kolmogorov-Smirnov测试，得到K-S统计检验量为0.025，说明 d{R_i} 分布分析结果的正确性.
![]()
图 3 单指数谱仿真数据功率谱与MLE理论谱比值的统计分布Fig. 3 Statistic distribution of simulated data for one-component power spectrum and that of MLE theory图4给出了每次仿真数据的U和{p_2}最优估计结果及其置信度范围. 可以看出，对于每组数据，最优估计结果在参数真实值附近变化，估计结果的90%置信度范围均能涵盖参数真实值，说明了最优估计结果的有效性.
图5进一步给出了对谱强度U和谱指数p2所有仿真数据最优参数估计的统计分析结果. 可以看出，绝大部分估计参数落入椭圆内，椭圆中心（即最优参数估计的平均值）非常接近参数理论值，表明利用MLE方法对闪烁单指数功率谱进行最优拟合，可以实现对谱参数的无偏估计.
![]()
图 4 U和p2的MLE结果及其置信度范围Fig. 4 Estimated MLE results of U and p2 and its confidence intervals3.2 双指数谱估计结果分析
双指数谱分析中，设置谱强度U = 0.6，谱指数{p_1} = 2.5，{p_2} = 3.5，菲涅尔频率{f_{\text{F}}}{\text{ = }}2\;{\text{Hz}}，归一化中断波数\;{\mu _0}{\text{=}}5. 利用相位屏方法对闪烁信号进行仿真，并对仿真的闪烁信号进行谱参数MLE，获得的闪烁信号时间变化序列如图6所示，长度为5 min，计算相应的闪烁指数{S _4} = 0.55.
![]()
图 6 仿真的双指数谱GNSS闪烁信号时间序列Fig. 6 Simulated GNSS signal time series under scintillation with two-component power spectrum最优估计过程同单指数谱仿真，图7给出了图6中闪烁信号时间序列的谱分析和谱拟合结果，蓝色曲线为仿真数据的谱分析结果，粉色曲线为增加噪声模型后的功率谱最优估计结果，红色曲线为不考虑噪声时的理论谱模型结果. 可以看出，通过MLE获得的理论谱模型可较好地对实测功率谱进行拟合.
![]()
图 7 仿真的GNSS闪烁信号双指数功率谱及MLE理论功率谱曲线Fig. 7 Two-component spectrum of simulated GNSS scintillation signals, and that of MLE theory利用MLE方法，对该仿真数据的功率谱进行拟合获得的最优参数估计结果为U = 0.61，{p_1} = 2.497，{p_2} = 3.480，与理论值相近. 利用最优谱拟合曲线计算闪烁指数，得到{S _4} = 0.56，与实际仿真数据计算结果接近. 上述结果表明MLE方法在双指数谱形式下估计的有效性.
图8给出了每次仿真数据的U, {p_1}, {p_2}, \;{\mu _0}参数MLE估计结果及其置信度范围. 可以看出，最优估计结果仍在参数真实值附近变化，但变化范围要大于单参数谱情况，且绝大部分估计结果的90%置信度范围能涵盖参数真实值，表明在双指数谱情况下仍可以获得有效的最优估计结果.
![]()
图 8 U, p1, p2 , µ 0的MLE结果及其置信度范围Fig. 8 MLE results and its confidence intervals for U, p1, p2 and µ 0图9给出了对所有仿真数据3个最优参数估计的统计结果. 可以看出，绝大部分参数估计值落入椭圆内，椭圆中心与各参数理论值略有偏差，表明双参数谱情况下利用MLE方法仍可以获得闪烁功率谱参数的近似最优无偏估计.
![]()
图 9 U和p1, p2, µ 0的MLE结果及其协方差椭圆Fig. 9 MLE of p1, p2 and µ 0 and its covariance ellipses with respect to optimized U estimation图10给出了图7数据所得功率谱和最优估计理论谱比值的统计分布. 可以看出， d{R_i} 统计分布与卡方分布变化一致. 对 d{R_i} 分布进一步进行Kolmogorov-Smirnov测试，得到K-S统计检验量为0.033，验证了 d{R_i} 分析结果的正确性.
![]()
图 10 双指数谱仿真数据功率谱与MLE理论谱比值的统计分布Fig. 10 Statistic distribution of simulated data for two-component power law spectrum and that of MLE theory对每个参数估计时，会存在个别参数估计结果的置信度范围不能涵盖参数理论值的现象，如图8中红线所示，说明增加参数估计个数会在一定程度上降低参数估计的可信度. 但对绝大部分数据参数估计结果是有效的. 说明对实际观测数据进行分析时，可通过选择不同初始值及增加估计次数，提高对谱参数估计的有效性.
3.3 实测数据结果分析
图11为三亚2011-10-24闪烁影响下的GPS L1信号时间变化序列，原始数据采样频率50 Hz[4]. 图中截取了从13:21UT开始长度为5 min的闪烁数据，计算获得的闪烁指数{S _4} = 0.67，表明该事件为一次强闪烁事件[29].
对实测数据利用Welch方法进行谱分析，数据分成5个数据段，每段长度1 min. 对实测功率谱利用单参数功率模型进行MLE，得到U = 1.18, {p_2} = 3.59, {f_{\text{F}}}{\text{ = }}0.55，各参数90%置信度区间分别为[0.98, 1.46]、[3.38, 3.81]和[0.49, 0.63]. 利用最优谱拟合曲线计算得到{S _4} = 0.64，与实际计算值相近. 图12给出利用最优估计参数获得的信号强度考虑和没有考虑噪声的理论功率谱与实测功率谱，可以看出理论谱可以很好地拟合实际谱.
![]()
图 12 实测GPS信号的闪烁功率谱及其最优拟合结果Fig. 12 Spectrum of real GPS scintillation data, and that of theoretical models from optimal estimation进一步对实测功率谱与最优估计获得的理论谱进行分析，图13给出了实测功率谱和最优估计理论谱比值 d{R_i} 的统计分布. 可以看出， d{R_i} 统计分布符合自由度为10的卡方分布. 对 d{R_i} 分布进行Kolmogorov-Smirnov测试，得到K-S统计检验量为0.042，表明对 d{R_i} 分布分析的正确性.
![]()
图 13 实测功率谱与MLE理论谱比值的统计分布Fig. 13 Statistic distribution of real data and that of MLE spectrum4 结 论
不均匀体参数是闪烁建模预测的重要输入参量. 在双指数电离层闪烁谱模型基础上，本文提出针对GNSS单站常规闪烁观测，利用对闪烁信号功率谱进行最优拟合，实现对多个谱参数进行MLE的方法.
利用Monte Carlo方法分别对单指数谱和双指数谱情况进行MLE分析. 结果表明：利用MLE方法，在单指数谱模型情况下可以实现对谱强度和谱指数的无偏估计；在双指数谱模型情况下，仍可以实现对谱强度和谱指数的近似无偏估计. 将最优估计结果带入闪烁谱模型并计算闪烁指数S4，得到的理论结果与利用实际闪烁信号计算的结果非常接近；利用K-S检验对实际数据功率谱与理论谱比值进行分析， d{R_i} 统计分布符合预期的卡方分布. 单指数谱模型每次估计结果的90%置信度范围均能涵盖参数真实值，双指数谱模型绝大部分估计结果的90%置信度范围能涵盖参数真实值，个别参数估计结果的置信度范围不能涵盖参数理论值的现象， 说明随着估计参数的增多，会出现个别情况下估计结果偏差增大的现象. 针对这种情况，实际分析中可以通过对数据进行多次估计并进行平均加以消除. 利用实测数据对谱参数最优估计方法进一步分析，强闪烁情况下可以获得闪烁谱参数的有效估计. 同时，仿真数据和实测数据的分析结果表明，通过谱拟合进行闪烁谱参数最优估计的方法适用于中等强度闪烁和强闪烁情况. 上述结果表明，通过谱拟合对电离层不均匀体参数进行MLE，是一种有效的闪烁谱参数估计方法.
现有单站电离层闪烁观测中仅利用高频采样测量获得闪烁指数，本文实现的方法为进一步利用单站观测数据进行闪烁参数（包括谱强度、谱指数，甚至漂移速度）提取和建模提供了可能.
文中分析的强闪烁事件中没有发生GNSS信号失锁情况，因此可以获得较好的谱分析结果. 若存在信号失锁或干扰等情况，信号闪烁谱可能会发生畸变，进而影响谱拟合过程及估计的结果. 同样，对于弱闪烁事件也需进一步明确其频谱是否存在不利于谱拟合最优估计的情况. 对基于谱拟合进行闪烁参数最优估计的适用范围进行分析，将是后继工作的一个重点.
-
![]()
图 1 仿真的单指数谱GNSS闪烁信号时间序列
Fig. 1 Simulated GNSS signal time series under scintillation with one-component power law spectrum
![]()
图 2 仿真的GNSS闪烁信号单指数功率谱及MLE理论功率谱曲线
Fig. 2 One-component power spectrum of simulated GNSS scintillation signals and that of MLE theory
![]()
图 3 单指数谱仿真数据功率谱与MLE理论谱比值的统计分布
Fig. 3 Statistic distribution of simulated data for one-component power spectrum and that of MLE theory
![]()
图 4 U和p2的MLE结果及其置信度范围
Fig. 4 Estimated MLE results of U and p2 and its confidence intervals
![]()
图 6 仿真的双指数谱GNSS闪烁信号时间序列
Fig. 6 Simulated GNSS signal time series under scintillation with two-component power spectrum
![]()
图 7 仿真的GNSS闪烁信号双指数功率谱及MLE理论功率谱曲线
Fig. 7 Two-component spectrum of simulated GNSS scintillation signals, and that of MLE theory
![]()
图 8 U, p1, p2 , µ 0的MLE结果及其置信度范围
Fig. 8 MLE results and its confidence intervals for U, p1, p2 and µ 0
![]()
图 9 U和p1, p2, µ 0的MLE结果及其协方差椭圆
Fig. 9 MLE of p1, p2 and µ 0 and its covariance ellipses with respect to optimized U estimation
![]()
图 10 双指数谱仿真数据功率谱与MLE理论谱比值的统计分布
Fig. 10 Statistic distribution of simulated data for two-component power law spectrum and that of MLE theory
![]()
图 12 实测GPS信号的闪烁功率谱及其最优拟合结果
Fig. 12 Spectrum of real GPS scintillation data, and that of theoretical models from optimal estimation
-
[1] KINTNER P M, LEDVINA B M, DE PAULA E R, et al. Size, shape, orientation, speed, and duration of GPS equatorial anomaly scintillations: GPS equatorial anomaly scintillations[J]. Radio science, 2004, 39(2). DOI: 10.1029/2003RS002878
[2] LEDVINA B M, KINTNER P M, DE PAULA E R. Understanding spaced-receiver zonal velocity estimation[J]. Journal of geophysics research,2004,109(A10):A10306. doi: 10.1029/2004JA010489
[3] CARRANO C S, GROVES K M, RINO C L. On the relationship between the rate of change of total electron content index (ROTI), irregularity strength (CkL), and the scintillation index (S4)[J]. Journal of geophysics research, space physics,2019,124:2099-2112. doi: 10.1029/2018JA026353
[4] LIU D, YU X, WANG Y, et al. Validation of relationship between ROTI and S4 with array of GNSS measurements[C]// Proceedings of China Satellite Navigation Conference. Singapore: Springer, 2022.
[5] SECAN J A, BUSSEY R M, FREMOW E J, et al. An improved model of equatorial scintillation[J]. Radio science,1995,30:607-617. doi: 10.1029/94RS03172
[6] CERVERA M A, THOMAS R M, K. GROVES M, et al. Validation of WBMOD in the Southeast Asian region[J]. Radio science,2001,36:1559-1571. doi: 10.1029/2000RS002520
[7] BÉNIGUEL Y. Global ionospheric propagation model (GIM): a propagation model for scintillations of transmitted signals[J]. Radio science, 37(3): 1032.
[8] CARRANO C S. The roles of diffractive and refractive scattering in the generation of ionospheric scintillation[C]// ISR Seminar, 2014.
[9] BEACH T L. Global Positioning System studies of equatorial scintillations[D]. Cornell University, 1998.
[10] RUDDICK B, ANIS A, THOMPSON K. Maximum likelihood spectral fitting: the Batchelor spectrum[J]. Journal of atmospheric and oceanic technology,2000,17(11):1541-1555. doi: 10.1175/1520-0426(2000)017<1541:MLSFTB>2.0.CO;2
[11] BARRET D, VAUGHAN S. Maximum likelihood fitting of X-ray power density spectra: application to high-frequency quasi-periodic oscillations from the neutron star X-binary 4U1608-522[J]. The astrophysical journal, 2012, 746(2). DOI: 10.1088/0004-637X/746/2/131
[12] VAUGHAN S. A simple test for periodic signals in red noise[J]. Astronomy & astrophysics manuscript, 431, 391-403,DOI: 10.1051/0004-6361:20041453
[13] RINO C. A power law phase screen model for ionospheric scintillation: 1 weak scatter[J]. Radio science,1979,14:1135-1145. doi: 10.1029/RS014i006p01135
[14] GHAFOORI F. Modeling the impact of equatorial ionospheric irregularities on GPS receiver performance[D]. University of Calgary, 2012.
[15] CARRANO C S, RINO C. A theory of scintillation for two-component power law irregularity spectra: overview and numerical results[J]. Radio science,2015,51:789-813. DOI: 10.1002/2015RS005903
[16] CARRANO C S, VALLADARES C E, GROVES K M. Latitudinal and local time variation of ionospheric turbulence parameters during the conjugate point equatorial experiment in Brazil[J]. International journal of geophysics, 2012: 103963. DOI: 10.1155 /2012/103963
[17] RINO C. A power law phase screen model for ionospheric scintillation: 2 strong scatter[J]. Radio science,1979,14:1147-1155. doi: 10.1029/RS014i006p01147
[18] RINO C, BREITSCH B, MORTON Y, et al. A compact multi-frequency GNSS scintillation model[J]. Navigation, 2018: 1-7.
[19] RINO C, CARRANO C. On the characterization of intermediate-scale ionospheric structure[J]. Radio science,2018,53:1316-1327. doi: 10.1029/2018RS006709
[20] JIAO Y, RINO C L, MORTON Y T. Ionospheric scintillation simulation on equatorial GPS signals for dynamic platforms[J]. Navigation,2018,65(2):263-274. DOI: 10.1002/navi.231
[21] XU D, MORTON Y T, YANG R. A two-parameter multi-frequency GPS signal simulator for strong equatorial ionospheric scintillation: modeling and parameter characterization[J]. Navigation,2020,67(1):181-195. DOI: 10.1002/navi.350
[22] CARRANO C S, RINO C L, GROVES K M. Maximum likelihood estimation of phase screen parameters from ionospheric scintillation spectra[C]//The 15th International Ionospheric Effects Symposium, Alexandria, VA, May 9-11, 2017.
[23] PRIESTLEY M B. Spectral analysis and time series[M]. London: Academic Press, 1981.
[24] 陈希孺. 高等数理统计学[M]. 合肥: 中国科学技术大学出版社, 2020. [25] RICE J A. 数理统计与数据分析 [M]. 北京: 机械工业出版社, 2019. [26] 杜比 A. 蒙特卡洛方法在系统工程中的应用[M]. 西安: 西安交通大学出版社, 2007. [27] 刘钝, 於晓, 冯健, 等. 电离层闪烁对GNSS信号影响的相位屏仿真方法[J]. 电波科学学报,2016,31(4):632-638. doi: 10.13443/j.cjors.2015083102 LIU D, YU X, FENG J, et al. Simulating the impacts of ionospheric scintillation on GNSS signals with phase screen method[J]. Chinese journal of radio science,2016,31(4):632-638. (in Chinese) doi: 10.13443/j.cjors.2015083102
[28] 张德丰. MATLAB概率与数理统计分析[M]. 北京: 机械工业出版社, 2009. [29] ITU-R Recommendation, Ionospheric propagation data and prediction methods required for the design of satellite services and systems[S]. Geneva, 2012.
-
期刊类型引用(0)
其他类型引用(1)
下载:
计量
- 文章访问数: 123
- HTML全文浏览量: 14
- PDF下载量: 47
- 被引次数: 1
