2022, Vol. 54
2. 哈尔滨工程大学青岛船舶科技有限公司，山东 青岛 266000
2. Harbin Engineering University Qingdao Ship Science and Technology Co., Ltd., Qingdao 266000, Shandong, China
卫星在围绕地球运动过程中除受地球引力影响，还受其他摄动力影响，导致其偏离预定轨道。为使卫星在预定轨道上运行，通常采用机动的方式调整卫星轨道。而全球导航定位用户在进行精密定位时依赖高精度、高可靠的轨道产品，当卫星发生机动时，需要对机动卫星单独处理以生成高精度轨道信息，即使使用广播星历定位，也必须对卫星机动进行监测，否则存在定位风险。除此之外，当卫星轨道发生机动时，该卫星服务中断，降低了卫星的连续性和可用性，从而影响卫星导航系统的定位、导航和授时服务性能[1]。为尽可能减少卫星机动不可用时间以及实现准确切换定轨策略目的[2]，需要对卫星的机动周期进行实时监测。欧洲定轨中心(The center for orbit determination in Europe，CODE)发布GPS卫星机动信息，但至少延时1 d，无法满足实时应用需求，而美国海岸警卫导航中心(U.S.coast guard navigation center，NAVCEN)发布的GPS卫星机动信息，普通用户无获取权限。广播星历的健康标志位同样会给出卫星预警信息，但此预警信息并不仅仅针对卫星机动，且广播星历的健康标志会提前标志机动开始时间，滞后标志机动结束时间，极大减少了卫星可用时间。除此之外，广播星历还有可能出现错误或丢失情况[3]。
为能准确地给出卫星机动信息，许多学者[4-6]开展了相关研究。由于利用卫星星历信息和观测量信息进行卫星机动探测及周期确定不需要额外测量设备，故被广泛应用。部分学者[7]仅利用广播星历进行卫星机动探测，基于广播星历的轨道互差探测BDS轨道机动，但该方法主要应用于事后处理。基于广播星历的卫星机动探测方法受广播星历更新频率限制，导致卫星机动探测准确度不高。Huang等[8]提出了一种结合广播星历和伪距的卫星机动探测方法，该方法能显著提高探测实时性。为进一步提高卫星机动探测及周期确定的准确度，有学者提出利用更高测量精度的载波相位进行卫星机动探测。Qiao等[9]提出了利用载波相位三差残差探测BDS卫星机动周期，该方法能准确探测卫星机动周期，其通过线性回归确定卫星周期，但其机动探测结果可能受基线选取的影响[10]。除此之外，还有学者利用轨道监测数据探测卫星机动，杜兰等[11]利用中国区域定位系统(China area positioning system，CAPS)的轨道监测数据确定GEO卫星的轨道机动，但普通用户没有获取该数据的权限。
小波变换是数字信号处理领域中一种对信号进行时频分析和处理的重要方法，可有效在时域、频域对信号进行局部化分析。本文基于载波相位三差模型进行卫星机动探测及周期确定，考虑差分测量噪声的影响，提出利用离散小波变换进行卫星机动分析，基于多分辨分析对载波相位三差残差进行分解，通过分析不同频率特性的近似系数和细节系数，提取卫星机动特性，综合载波相位三差及小波变换的优势，实现卫星机动快速探测及周期确定。
1 离散小波变换的卫星机动探测 1.1 载波相位三差模型在t时刻载波相位观测量方程:
| $ \begin{aligned} \lambda \phi(t)= & r(t)+T(t)+c\left(d t_u(t)-d t^s(t)\right)- \\ & I(t)-\lambda N+\delta(t) \end{aligned} $ | (1) |
式中：r(t)为卫星与接收机间的距离，T(t)、I(t)分别为对流层和电离层延迟误差，dtu(t)、dts(t)分别为接收机钟差及卫星钟差，N为整周模糊度，λ为对应的波长，δ(t)为测量噪声，c为光速。
在周跳探测策略上，采用电离层残差法和M-W组合联合探测；在误差处理策略上，采用无电离层组合消除电离层延迟；Saastamoinen模型及Niell投影函数补偿对流层延迟；采用星间和站间差分分别消除与接收机及卫星有关的误差；为进一步消除整周模糊度，将观测方程进行历元间差分：
| $ \begin{aligned} & \lambda^{\prime}\left(\phi_{M N}^{j k}\left(t_k\right)-\phi_{M N}^{j k}\left(t_{k-1}\right)\right)= \\ & r_{M N}^{j k}\left(t_k\right)-r_{M N}^{j k}\left(t_{k-1}\right)+\delta\left(t_k\right) \end{aligned} $ | (2) |
式中：M、N为不同测站，j、k为不同卫星，λ′为经过无电离层组合后的波长，
基于载波相位三差观测方程，提取载波相位三差残差量：
| $ \begin{aligned} \delta\left(t_k\right)= & \lambda^{\prime}\left(\phi_{M N}^{j k}\left(t_k\right)-\phi_{M N}^{j k}\left(t_{k-1}\right)\right)- \\ & r_{M N}^{j k}\left(t_k\right)+r_{M N}^{j k}\left(t_{k-1}\right)+\eta\left(t_k\right) \end{aligned} $ | (3) |
式中η(tk)为与卫星机动相关的变量。
在卫星未发生机动的情况下，η(tk)为0，由式(3)得到的载波相位三差残差相对较小；卫星机动发生后，η(tk)不再为0，广播星历不能很好地预测卫星实际轨道，载波相位三差残差会逐渐增大，机动结束后，卫星轨道又趋于平稳。利用载波相位三差残差在卫星机动发生及结束时刻的变化特性，实现卫星机动探测并确定机动周期。需要注意的是，虽然式(2)中的载波相位三差消除了大部分观测误差，但无电离层组合以及三次差分的操作都会增加测量噪声，如果此时不考虑测量噪声的影响，将导致不能快速准确地探测卫星机动周期，故在利用载波相位三差残差模型时必须考虑测量噪声的影响。
1.2 基于离散小波变换的卫星机动信息确定传统基于载波相位三差模型的卫星机动探测方法采用多历元线性回归减弱测量噪声，但多历元累积的方法增加了机动探测反应时间，且该方法的机动探测结果可能受基线选取的影响。考虑到机动信号和载波相位三差残差噪声进行小波变换后在各尺度上所表现的特性不同，本文利用小波变换对载波相位三差残差去噪，消除噪声对机动信号的影响，以此探测其在机动开始时刻和结束时刻的局部变化特性，实现卫星机动快速探测及周期确定。
根据式(3)，载波相位三差残差包含与卫星机动相关的变量和噪声两部分，可将其简写为δ(t)=η(t)+ε(t)，其中η(t)为与卫星机动相关的变量，ε(t)为载波相位三差残差噪声，且近似服从零均值正态分布[12]。不妨设噪声ε(t)均值为0、方差为σ2，用ψ(t)、Wε(m, n)分别表示选取的小波函数以及ε(t)的小波变换，则载波相位三差残差噪声小波变换后的平均功率Wε(m, n)为
| $ \begin{aligned} & E\left\{\left|W_{\varepsilon}(m, n)\right|^2\right\}= \\ & \iint E\{\varepsilon(u) \varepsilon(v)\} \psi_m(n-u) \psi_m(n-v) \mathrm{d} u \mathrm{~d} v= \\ & \sigma^2\|\psi\|^2 / m \end{aligned} $ | (4) |
式中：‖ · ‖2为平方积分，m为尺度因子，n为平移因子，u、v为积分变量。从式(4)可得载波相位三差残差噪声ε(t)小波变换后的平均功率与尺度因子m成反比，这与机动变量η(t)小波变换后各尺度上表现的特性不同[13]，小波变换正是通过多分辨分析，最终实现去噪目的。
小波变换的多分辨分析建立在空间概念上，在不同的尺度下，构造尺度函数向量空间V和小波函数向量空间W，载波相位三差残差δ(t)在尺度空间V做卷积得到表示低频信息的近似系数，在小波空间W做卷积得到表示高频信息的细节系数。考虑到机动期间db1小波函数与载波相位三差残差的阶跃特性相似性较高[14]，选取db1小波函数作为小波基，对载波相位三差残差δ(t)进行分解:
| $ \begin{aligned} \delta(t)= & \sum\limits_n c_{j, n} \phi_{j, n}(t)+\sum\limits_{0<m \leqslant j} \sum\limits_n d_{m, n} \psi_{m, n}(t)= \\ & \sum\limits_n c_{j, n} \phi_{j, n}(t)+\sum\limits_n d_{1, n} \psi_{1, n}(t)+\cdots+ \\ & \sum\limits_n d_{j, n} \psi_{j, n}(t) \end{aligned} $ | (5) |
式中：j为分解层数，最大分解层数N满足2N≤Length(δ)，其中Length(δ)表示载波相位三差残差δ(t)长度，
理论上，分解层数越高去噪效果越好，但分解层数过高可能会丢失载波相位三差残差的机动信息[15]，考虑到第1层受噪声影响较为严重[16]，所以分解层数一般大于1。由近似系数分解成下一层分解层数的近似系数和细节系数[17]，如果近似系数和细节系数探测结果相同，表明此时由近似系数分解的细节系数已体现机动信息，继续分解有可能丢失机动信息，故本文以基于近似系数与细节系数的卫星机动探测结果为依据，机动开始时刻相同的最低分解层数作为卫星机动探测分解层数。小波变换将大部分能量集中在少数、幅值比较大的小波系数上，而载波相位三差残差的噪声分布在各个尺度上，其小波系数幅值较小，考虑到噪声近似服从零均值高斯分布[12]，利用Donoho[18]提出的阈值公式作为阈值, 即
| $ T=\lambda \sqrt{2 \ln (n)} $ | (6) |
式中：n为三差残差长度; λ为噪声标准方差，其估计公式[19]：λ=median|ri|/0.674 5, 其中γi为小波分解第1层分解所对应的小波系数，median为求中值。考虑细节系数随分解层数的增加而衰减，故对细节系数的噪声标准方差λ做进一步处理[16]：λ′=λ/In(J+1)，其中，J为当前分解层数。
离散小波变换能有效消除噪声对三差残差的影响，放大卫星机动起止时刻三差残差的变化幅值，根据卫星机动未发生时的三差残差数据确定的判断阈值，可有效实现卫星机动探测。基于离散小波变换进行卫星机动探测及周期确定时，综合利用载波相位三差和小波变换优势，实现快速准确的卫星机动信息探测。
2 结果分析为了验证本文所提方法性能，本文设计了两个实验，利用不同基线长度测站的观测数据，分别针对GPS和BDS卫星进行卫星机动探测及周期确定，并与载波相位三差的线性回归方法进行比较。对于GPS卫星，利用IGS的METG、REDU、TLSE、WTZ3、WTZZ及DYNG共6个测站，其中DYNG作为主测站，与其他测站构成5条基线。对于BDS卫星，利用IGS的GMSD、JFNG、MIZU、PTGG及ANMG共5个测站，其中ANMG作为主测站，与其他测站构成4条基线。实验所选择的观测数据见表 1，数据采样率为30 s。
|
表 1 实验所选取的卫星观测数据 Tab. 1 Satellite observation data selected in the experiment |
对于实验1，根据CODE所公布的卫星机动信息，选取真实发生的G26机动卫星，取其机动时观测数据，构建载波相位三差探测模型，利用离散小波变换进行卫星机动探测及周期确定。图 1、2分别给出了基于线性回归和利用本文所提方法的机动探测及周期确定结果。图 1中，图 1(a)~(d)分别给出了5条基线所对应的载波相位三差残差曲线、以及利用2~4个历元做线性回归的机动探测曲线。红色实线为机动探测周期结果，红色虚线为判断阈值[1, 7, 20]。图 2中，图 2(a)、(c)、(e)分别给出了分解层数为3层的近似系数d3、d2及d1曲线；图 2(b)、(d)、(f)分别给出了分解层数为3层的细节系数a3、a2及a1曲线。红色实线为探测机动周期结果，红色虚线为判断阈值，其判断阈值按式(6)进行确定。
|
图 1 基于线性回归方法探测GPS卫星机动及周期确定 Fig. 1 GPS satellite maneuver detection and period determination based on linear regression method |
|
图 2 基于离散小波变换探测GPS卫星机动及周期确定 Fig. 2 GPS satellite maneuver detection and period determination based on discrete wavelet transform |
表 2给出了CODE所公布的GPS卫星G26的卫星机动信息、传统基于线性回归方法的机动探测结果，以及利用本文所提方法的机动探测结果，其中，线性回归方法的探测结果选取的是反应最敏感的TLSE-DYNG基线结果。考虑到GPS卫星机动时间一般不超过2 min[21]，本文在利用线性回归方法时最多选取了4个历元。
|
|
表 2 GPS机动探测及周期确定结果 Tab. 2 GPS maneuver detection and period determination results |
以表 2中CODE公布的机动开始时间为真值，从结果上看，采用离散小波变换探测机动开始时间滞后真值1 s；采用2~4个历元做线性回归探测机动开始时间分别滞后真值1 s、提前真值29 s和提前真值59 s。虽然CODE没有公布卫星机动的结束时间，但由图 1(b)可知，利用2个历元做线性回归探测的机动周期为一个点，与机动周期为一段时间的常识有矛盾之处，故基于线性回归存在提前确定卫星机动结束时间的可能。对于线性回归方法而言，线性回归选取历元较少时，斜率易受噪声影响导致误判风险；线性回归选取历元较多时，能有效降低因噪声导致的误判，但所需探测反应时间也会相应增加，如图 1及表 2结果所示，故基于线性回归的方法需在可靠性和时效性之间做一定权衡。离散小波变换能有效消除测量噪声对载波相位三差残差影响，放大载波相位三差残差的局部特性，实现平滑载波相位三差残差目的，且其单历元处理方式不需要额外增加探测反应时间。综上所述可得，相比于线性回归方法，基于离散小波变换方法探测卫星机动更快速、准确。
对于实验2，综合精密星历缺失性、广播星历健康标志和轨控指令信息，选取真实发生的C01卫星，取其机动时观测数据，建立载波相位三差探测模型，利用离散小波变换进行机动探测及周期确定。图 3、4分别给出了基于线性回归和利用本文所提方法的机动探测及周期确定结果。图 3(a)~(d)分别给出了载波相位三差残差曲线、以及以利用10、20和30个历元做线性回归的机动探测曲线，红色实线为探测机动周期结果，红色虚线为判断阈值。图 4(a)、(c)、(e)、(g)、(i)分别给出了分解层数为5层的近似系数d5、d4、d3、d2、d1曲线；图 4(b)、(d)、(f)、(h)、(j)分别给出了分解层数为5层的细节系数a5、a4、a3、a2、a1曲线。红色实线为探测机动周期结果，红色虚线为判断阈值。
|
图 3 基于线性回归方法探测BDS卫星机动及周期确定 Fig. 3 BDS satellite maneuver detection and period determination based on linear regression method |
|
图 4 基于小波变换探测BDS卫星机动及周期确定 Fig. 4 BDS satellite maneuver detection and period determination based on discrete wavelet transform |
表 3给出了关于C01卫星的轨控指令、传统基于线性回归方法的探测结果，以及利用本文所提方法的探测结果。其中，基于线性回归方法探测机动及周期的结果选取的是反应最敏感的ANMG-GMSD基线结果。
|
|
表 3 BDS机动探测及周期确定结果 Tab. 3 BDS maneuver detection and period determination results |
由表 3结果可知，利用小波变换和线性回归方法得到的BDS卫星机动周期都在轨控指令的机动时间段内。由于未查到BDS卫星公布的机动信息，所以无法给出探测结果的准确度，但就结果而言，BDS探测结果与GPS探测结果保持一致。通过对图 3对比，利用线性回归探测卫星机动时其机动探测结果可能受基线选取的影响。图 4表明，利用小波变换探测卫星机动时其不同基线探测到的卫星机动发生时刻及周期结果一致。综上所述，相比于传统的线性回归方法，离散小波变换在进行卫星机动快速探测及周期确定时具有一定优势。
3 结论1) 相比于传统基于线性回归的方法，基于离散小波变换探测GPS卫星机动开始时刻偏离真值1 s，优于传统方法。
2) 不同基线探测到的卫星机动发生时刻及周期结果一致，单历元处理方式不需要额外增加探测反应时间，能更及时准确的探测卫星机动并确定其周期。
| [1] |
YE Fei, YUAN Yuanbin, ZHANG Baocheng, et al. A three-step method for determining unhealthy time period of GPS satellite orbit in broadcast ephemeris and its preliminary applications for precise orbit determination[J]. Remote Sensing, 2019, 11(9): 1098. DOI:10.3390/rs11091098 |
| [2] |
QIN Zhiwei, HUANG Guanwen, ZHANG Qin, et al. A method to determine BeiDou GEO/IGSO orbital maneuver time periods[J]. Sensors, 2019, 19(12): 2675. DOI:10.3390/s19122675 |
| [3] |
HENG Liang, GAO G X, WALTER T, et al. GPS signal-in-space integrity performance evolution in the last decade[J]. IEEE Transactions on Aerospace and Electronic Systems, 2012, 48(4): 2932. DOI:10.1109/TAES.2012.6324670 |
| [4] |
QIN Wu, HUANG Hengli, ZHANG Yiyi, et al. Precise orbit determination for BeiDou GEO/IGSO satellites during orbit maneuvering with pseudo-stochastic pulses[J]. Remote Sensing, 2019, 11(21): 2587. DOI:10.3390/rs11212587 |
| [5] |
ZHOU Peiyuan, DU Lan, LI Xiaojie, et al. Near real-time BDS GEO satellite orbit determination and maneuver analysis with reversed point positioning[J]. Advances in Space Research, 2019, 63(5): 1781. DOI:10.1016/j.asr.2018.10.046 |
| [6] |
JU Bing, GU Defeng, HERRING T A, et al. Precise orbit and baseline determination for maneuvering low earth orbiters[J]. GPS Solutions, 2017, 21(1): 53. DOI:10.1007/s10291-015-0505-x |
| [7] |
YE Fei, YUAN Yunbin, TAN Bingfeng, et al. A robust method to detect BeiDou navigation satellite system orbit maneuvering/anomalies and its applications to precise orbit determination[J]. Sensors, 2017, 17(5): 1129. DOI:10.3390/s17051129 |
| [8] |
HUANG Guanwen, QIN Zhiwen, ZHANG Qin, et al. A real-time robust method to detect BeiDou GEO/IGSO orbital maneuvers[J]. Sensors, 2017, 17(12): 2761. DOI:10.3390/s17122761 |
| [9] |
QIAO Jing, CHEN Wu. Beidou satellite maneuver thrust force estimation for precise orbit determination[J]. GPS Solutions, 2018, 22(2): 42. DOI:10.1007/s10291-018-0705-2 |
| [10] |
陈倩. 缩短北斗卫星机动不可用时间的定轨方法研究[D]. 北京: 中国科学院大学, 2019 CHEN Qian. BDS orbit determination method during the period of satellite maneuver[D]. Beijing: University of Chinese Academy of Sciences, 2019 |
| [11] |
杜兰, 张中凯, 李晓杰, 等. GEO卫星东西机动监视与移动窗口轨迹拟合[J]. 测绘学报, 2014, 43(3): 233. DU Lan, ZHANG Zhongkai, LI Xiaojie, et al. Station-keeping maneuver monitoring and moving-window ground track fitting of GEO satellites[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(3): 233. DOI:10.13485/j.cnki.11-2089.2014.0034 |
| [12] |
ZHU Yongxing, ZHENG Kang, CUI Xianqiang, et al. Preliminary analysis of the quality and positioning performance of BDS-3 global interoperable signal B1C&B2a[J]. Advances in Space Research, 2021, 67(8): 2483. DOI:10.1016/j.asr.2021.01.045 |
| [13] |
PUGLIANO G, ROBUSTELLI U, ROSSI F, et al. A new method for specular and diffuse pseudorange multipath error extraction using wavelet analysis[J]. GPS Solutions, 2016, 20(3): 499. DOI:10.1007/s10291-015-0458-0 |
| [14] |
CHUI C K, HAN Ningning. Wavelet thresholding for recovery of active sub-signals of a composite signal from its discrete samples[J]. Applied and Computational Harmonic Analysis, 2020, 52: 1. DOI:10.1016/j.acha.2020.11.003 |
| [15] |
ZHONG P, DING X L, ZHENG D W, et al. Adaptive wavelet transform based on cross-validation method and its application to GPS multipath mitigation[J]. GPS Solutions, 2008, 12(2): 109. DOI:10.1007/s10291-007-0071-y |
| [16] |
CHEN Yong, CHENG Yanan, LIU Huanlin. Application of improved wavelet adaptive threshold de-noising algorithm in FBG demodulation[J]. Optik, 2017, 132: 243. DOI:10.1016/j.ijleo.2016.12.052 |
| [17] |
BAILI J, LAHOUAR S, HERGLI M, et al. GPR signal de-noising by discrete wavelet transform[J]. NDT & E International, 2009, 42(8): 696. DOI:10.1016/j.ndteint.2009.06.003 |
| [18] |
DONOHO D L. De-noising by soft-thresholding[J]. IEEE Transactions on Information Theory, 1995, 41(3): 613. DOI:10.1109/18.382009 |
| [19] |
DONOHO D L, JOHNSTONE I M. Adapting to unknown smoothness via wavelet shrinkage[J]. Journal of the American Statistical Association, 1995, 90(432): 1200. DOI:10.1080/01621459.1995.10476626 |
| [20] |
YANG Y. Robust estimation of geodetic datum transformation[J]. Journal of Geodesy, 1999, 73(5): 268. DOI:10.1007/s001900050243 |
| [21] |
DAI Xiaolei, LOU Yidong, DAI Zhiqiang, et al. Precise orbit determination for GNSS maneuvering satellite with the constraint of a predicted clock[J]. Remote Sensing, 2019, 11(16): 1949. DOI:10.3390/rs11161949 |