在卫星导航系统中, 多径误差与系统周围环境有密切关系, 所以不能通过差分技术将其消除, 它是制约卫星导航系统定位精度提升的因素之一。多径干扰不仅影响着单点定位的精度, 而且还限制了高精度单点定位的收敛速度。在北斗卫星导航系统中, 码多径干扰同样也制约着北斗用户获得更高精度的定位结果。此外, 北斗系统有着与其他卫星导航系统不同的星座结构——5颗地球同步轨道卫星(GEO)、5颗倾斜地球同步轨道卫星(IGSO)和5颗中轨道卫星(MEO)。因为GEO卫星和IGSO卫星首次应用在卫星导航系统中, 所以关于这两种卫星有关的多径效应的研究成果相对较少。目前, 国际上主要应用GEO卫星作为星基增强系统(satellite-based augmentation system, SBAS)的组成部分, 如美国的广域增强系统(WAAS)、俄罗斯的差分校正和检测系统(SDCM)、日本的多功能星基增强系统(MSAS)以及印度的GPS辅助静止轨道增强导航系统(GAGAN)。文献[1]研究了WAAS系统中GEO卫星多径衰落特性, 并指出GEO卫星的伪码观测量遭受着比GPS卫星伪码观测量更为严重的多径干扰。文献[2-3]的研究成果表明, 在北斗系统中, GEO卫星的多径衰落幅度较之与IGSO和MEO卫星衰落幅度更大。
近年来, 有关抑制GEO卫星多径效应的研究成果层出不穷, 比如运用小波变换技术将多径误差信号中的低频分量滤除, 然后将其应用在之后几天的多径误差修正中。同时还发现北斗GEO多径的周期大致是一个恒星日(23小时56分04秒), 因此恒星日滤波器也被作为抑制GEO多径的一种有效手段[4]。尽管如此, 针对北斗GEO卫星多径产生机理的研究相对较少, 也就是说北斗GEO卫星的多径信号是发生在接收机周围还是发生在卫星表面目前尚不明晰。文献[5]在研究了GPS系统中SVN49的L1频点伪距异常现象后, 发现此异常现象来自于卫星发射天线的耦合电路网络, 并通过数学模型描述了卫星发射信号。因为伪距异常现象与卫星高度角相关, 所以在此数学模型中, 高度角是其中一个参数, 在卫星运动到不同的位置时, 信号的幅度是不同的。但此模型不适合于描述北斗GEO卫星的伪距异常现象, 因为GEO卫星与地面接收机保持相对静止。所以本文提出了卫星表面的信号反射模型, 并运用此模型推导出卫星多径曲线。同时运用“Morlet”小波分别对真实多径信号和建模多径信号进行分析。此外文章中采用数值优化算法计算模型中的反射参数。最后通过分析KARR站和XMIS站接收到的GEO卫星数据表明, 此模型能够很好地对真实多径信号做出预测。
1 GEO卫星表面多径发射模型北斗GEO卫星的结构如图 1所示, 其有效载荷的尺寸是1.8 m×2.2 m×2.5 m[6]。截止到2016年3月30日, 所有在轨运行的GEO卫星工作在DFH-3A通信平台上, 并由长征3号丙发射升空。由图 1所示, GEO卫星的导航天线阵列是由19个天线单元组成。
正如前文所提到的, André Hauschild博士的研究成果表明[5], SVN49卫星的伪距异常来源于L1频点的导航信号在发射天线的耦合单元发生了反射, 使得用户接收到的信号中固有的存在一个大约30 ns延迟的多径信号。这种伪码异常现象也在Mauna Loa(Hawaii)和Weilheim(Germany)等地区被观测到。André Hauschild博士还通过数学模型解析式来描述此异常现象, 但是, 此模型对于北斗导航系统的GEO卫星码多径异常现象是不适用的。因此, 本文提出了基于卫星表面的多径反射模型(见图 2), 同时关于建模的合理性也在文中进行了阐述。
在图 2中, S表示GEO卫星导航天线阵列中的一个单元。假设射频信号从天线发射出来后, 在距离发射源S点h米处的反射面Rs上发生反射, 此反射面的法向量与水平线lhorizons的夹角是θ。此外, ld和lm分别代表直射信号路径和反射信号路径, α是卫星的高度角。GEO卫星距离地面大致是35 786 km, 而射频信号在卫星发射天线周围的反射距离只有10 m左右, 所以反射路径与直射路径的夹角γ非常小, 也就是说直射信号的方位角和高度角与反射信号的方位角和高度角几乎相等。所以在计算路径差lm-d时, 计算难度从三维空间几何降低到平面几何, 那么lm-d的表达式如下
(1) |
在多径环境中, 接收机很难将同时收到的直射信号与反射信号分离开。假设接收机的射频前端带宽无限宽, 且没有噪声干扰, 那么在码相位跟踪阶段, 本地的码与中频信号的相关峰不再成等腰三角形, 其数学表达式为
(2) |
式中, Rd(·)为直射信号伪随机噪声码的相关函数, αn、τn以及ϕn分别代表第n条多径信号相对于直射信号的振幅、时延以及载波相位。由于直射信号相关峰的对称性被多径信号破坏, 这就使得延迟锁定环(DLL)不能将本地信号与直射信号的码相位对齐, 这就是造成伪码误差的根本原因。
非相干超前减滞后功率鉴相器被广泛地应用于实际工程当中, 所以本文以此鉴相器为例推导出因多径干扰所产生的码跟踪误差, 数学表达式为
(3) |
式中, E(t)和L(t)分别表示超前与滞后支路的相关幅值。假设接收机实际收到的信号s(t)是直射波及其一路反射波地叠加, 则s(t)可表达为
(4) |
式中,p(·)代表为赋值是±1的伪随机序列码, A和ω0分别指代直射信号的幅度和角速度, 参数α、τ和ϕ的定义与公式(2) 相同。本地接收机中超前支路和滞后支路地码片序列可表示为
(5) |
式中,ϕc是多径信号产生的载波相位误差, δP是锁相环(PLL)的跟踪误差, δD是延迟锁定环(DLL)的跟踪误差, Δ是超前-滞后支路之间的码间距, ωl是本地载波NCO的角频率。为了剥离输入信号s(t)的载波分量, 本地载波sE(t)和sL(t)同时与s(t)相乘, 被滤除高频分量后的混频结果分别如下
(6) |
式中,ωe=ωl-ω0是输入信号与本地信号之间的载波频率差异, Re{·}代表对花括号中的内容取实部。假设PLL进入锁定状态, 并且载波相位误差ωe=0, 这样(6) 式可以被重写为
(7) |
将(7) 式代入到(3) 式中, 可得
(8) |
在公式(8) 中, 等号右边的前两项可以被看作是码多径误差的期望值, 在最后一项中, α和方括号中的多项式决定了码多径误差的幅值, 余弦项cos(ϕm)反应了伪码误差的变换趋势, 其中ϕm是多径信号的载波相对相位, 它与多径时延有关, 随时间而变化。所以多径信号的相对相位ϕm可以表示为
(9) |
式中,λ是信号载波的波长。
3 实测数据仿真本节将对2016年1月12日分别在KARR站(位于西澳大利亚)以及XMIS站(位于圣诞岛, 澳大利亚海外领土)接收到的真实卫星数据给予分析(数据由IGS网站提供), 这两个观测站周围环境开阔, 无高大建筑物。通过“码减载波”的方法提取北斗GEO卫星的多径误差信息, 并提出基于交叉验证的小波降噪算法分离多径误差的低频分量。最后由数值优化算法计算反射模型参数。
3.1 码多径误差北斗GEO码多径误差随时间变化, 并且以恒星日为周期呈现规律性变化, 其多径观测量的表达式如下
(10) |
式中
(11) |
式中, MP*是多径观测量, ρ*和ϕ*分别代表码观测量和载波观测量, 下标i和j表示不同频点, 线性因子α是不同载波频率的比值, λ*是相应频点的载波波长, N*是相应频点载波的周整模糊度, l*是相应频点载波的多径误差。此外, B*主要由周整模糊度线性组合后产生的误差, ε*包括载波多径线性组合误差以及其他更小的误差ξ, 例如硬件延迟以及码观测量和载波观测量中的噪声误差。理论上, 载波多径误差最大不会超过载波波长的1/4, 它比码多径误差低大约2个数量级, 故可以忽略不计。假设接收到的载波观测量没有发生周跳, 那么B*就保持一个常数, 它可以被认为是多径观测量MP*的期望。但在实际工程中, 载波周跳是很难避免的, 所以在分析MP*时, 必须先要对周跳进行检测和恢复。
3.2 基于交叉验证小波降噪算法交叉验证(cross-validation, CV)算法是统计学中一种将数据样本分类为小子集的实用方法, 其基本思想是将原始数据信息分为训练子集(training set)和验证子集(validation set), 用训练子集中对模型进行训练, 再通过验证子集检验模型的准确性[7]。CV算法常用于人工智能、机器学习、模式识别、数据挖掘等领域[8-9]。下面将以二次回归模型和三次样条插值模型为例, 简单地介绍交叉验证算法的基本原理, 以及引出文本提出的基于2-fold cross validation小波去噪算法的思想。
本文选择Symlet小波基函数对北斗多径信号进行分解, 算法步骤如下:
Step1 利用FFT对原始多径误差信号s(n)进行频谱分析;
Step2 对原始多径误差信号进行K层小波分解;
Step3 对k(k=1, 2, 3, …K)层到K层的细节系数以及第K层的近似系数, 进行重构信号s′k(n); 。
Step4 对重构信号进行FFT计算, 如果重构信号缺少原始信号的中的频谱分类, 则K=K+1, 并回到Step2;
Step5 对重构信号s′k(n)进行cross-validation验证, 选择奇数采样点s′k(nodd)作为训练集, 并利用三次样条插值函数对信号进行拟合, 令拟合信号为fk(n)。利用验证集s′k(neven)检验拟合信号fk(n), 并计算下式
(12) |
Step6 如果C(k)的值达到最小就可以停止计算, 并定义当前分解层数K为最佳分解层数, 且将k层及k层以上的细节系数与近似系数进行重构计算获取的s′k(n)信号输出, 得到降噪信号。
根据以上算法对KARR和XMIS站点的GEO卫星多径误差进行降噪处理, 处理结果如图 3所示。图中浅色曲线表示原始多径误差, 深色曲线表示降噪后的多径误差低频分量。
本文的研究对象是来自卫星星体上的多径信号, 那么为了尽可能地摆脱来自接收机天线周围的反射信号, 有意地选择了KARR和XMIS处于开阔环境中的观测站点。此外, 一个有利的条件就是两观测站都装配了Trimble NetR9接收机以及TRM 59 800.00 NONE天线, 这使得在分析结果更具有客观性和说服力。图 3中的曲线表明, 同一观测站不同GEO卫星的多径误差衰落特性是不同的。在图 3a)中C01、C04和C05卫星对应的多径误差衰落频率要高于其他2颗GEO卫星对应的多径曲线。类似的现象也出现在XMIS站(见图 3b))。特别的, 以KARR观测站的分析结果为例, C01、C02、C04以及C05卫星的多径误差波形分别在14 h、11 h、15 h和11 h处呈现出一定的对称性。同样的, 在XMIS站提取的多径误差信号中, 这种对称性也出现在C04以及C05号卫星中。此外, 观察两观测站对应卫星号的多径序列可以发现, 同一颗卫星在同一时间段内的码多径衰落特性有明显的区别。考虑到2个观测站所处的地形地貌较为相似, 如果导致北斗卫星GEO码多径异常现象的原因与GPS系统中SVN49卫星类似, 是因为射频信号在天线耦合网络电路中出生发射, 那么在不同观测站提取到的相同GEO卫星的多径误差应该具有类似的衰落特性。这也就是说, 本文提出的卫星表面反射模型是合理的, 导航信号沿着不同的反射路径传输到不同地区的观测站, 会导致同一颗卫星在不同地区呈现出的多径衰落的不一致现象。
3.3 计算模型反射参数为了计算反射模型中的反射参数, 本文提出了数值优化算法, 其优化目标使得实测多径数据的小波变换矩阵F与模型多径数据的小波变换矩阵G的相关系数最大。此外, 目标函数还应满足以下2个约束条件:反射角度θ不应超过180°, 反射距离h不应超过星体尺寸一半的宽度(包含太阳能电池帆板)。北斗GEO卫星最大宽度为17.7 m[6], 考虑到卫星天线的分布情况, 令反射距离h小于11 m。优化算法的数学表达式如下
(13) |
式中
(14) |
在公式(13) 中fmn和gmn分别代表矩阵F和G中的元素, F和G分别表示矩阵F和G的期望, tr(·)表示对矩阵求迹, 同时令
(15) |
式中,x(k)是离散序列, s和τ分别是小波基函数进行伸缩和平移变换因子。小波变换的过程可以理解为, 首先对小波基函数进行伸缩变换, 之后对离散序列与伸缩变换后的小波基函数进行卷积运算, 完成后再对小波基进行新伸缩变换, 之后再做卷积运算, 如此往复, 直到小波基函数完成所有的伸缩变换。
因为“Morlet”小波基函数具有非正交性, 而且是由Gaussian调制的指数复小波, 很利于对时间序列进行分析, 其数学表达式为
(16) |
式中,β决定了小波基函数的包络, 且当β增大时, 小波基函数的时间分辨率降低, 频率分辨率提高。如果β趋近于0, “Morlet”小波函数具有很高的频率解析力, 因为它的包络趋近于余弦函数, 反之当β趋近于无穷大时, “Morlet”小波函数近似于Dirac函数, 具有较好的时间分辨率。通常情况下, 取β=1。此外, ω0=6对小波分析结果的时间与频率的定位有一个较好的平衡性[10]。将(16) 式带入到(15) 式中, 可以得到
(17) |
那么矩阵F和G可通过(17) 式计算获得小波矩阵, 其同时包含着多径信号的频率成分, 以及频率出现的时间信息。那么根据以上数值优化算法, 可计算出反射模型中的参数h和θ。
以KARR站C05号卫星在2016年1月12日的数据为例, 进行多径误差反射参数的计算。图 4a)是通过“Morlet”小波基对实测多径信号低频分量的时频分析结果, 图 4b)是由优化算法计算出的仿真多径误差曲线的时频结果, 经计算得到的模型反射参数h和θ分别为9 m和121°。比较2幅图中对应的频率槽, 在相应的时间段内所包含的频率分量是接近的, 通过计算2幅图对应矩阵的相关系数, 其相似性高达0.954 3。类似地, 图 5显示的是实测多径误差和仿真多径误差的时域曲线。可以看出, 仿真多径误差曲线具有和实测多径误差曲线相似的波动特性, 在0 h到6 h和13.5 h到24 h这2个时间段内, 有较为明显的波动, 在6 h至13.5 h时段内, 波动较为平缓。
利用本文提出的卫星星体反射模型和优化算法对KARR站点的C04号卫星以及XMIS站点的C03号卫星的多径误差信号进行分析和反射参数计算, 结果如下:
图 6分别比较了KARR C04和XMIS C03 GEO卫星实测多径误差波形与仿真多径误差曲线。可以看出通过反射模型计算获得的多径误差曲线与原始多径误差曲线的波动特性基本一致。也就是说, 本文提出的卫星形体反射模型可以较好地对实测环境中的多径误差进行拟合。
4 结论北斗GEO卫星码多径观测量伴随着一个缓慢波动的误差, 与卫星高度角呈现出相关性, 并且以恒星日为周期不断地重复。文献[5]针对GPS系统的SVN49号卫星的码观测量异常现象的结论是, 射频信号在天线耦合网络中发生了发射, 形成了一个约30ns延迟的多径信号, 因此在发射源对直射信号形成了污染, 并且不同地区都观测到了类似的码多径异常现象。但是在对北斗GEO卫星数据进行分析后发现不同地区接收同一颗卫星的多径观测量呈现出不同的衰落特征, 所以说导致北斗GEO卫星码多径异常的原因与SVN 49是不同的。所以本文提出了基于实测数据的卫星多径反射模型, 并通过数值优化算法计算出了反射多径参数, 为了验证此模型的正确性, 文中通过小波变换算法对多径误差数据进行分析。通过时域数据对比和小波变化结果相关性分析发现, 本文提出的基于实测数据的卫星多径反射模型可以很好地对北斗GEO卫星实测多径误差信号进行预测, 这对抑制GEO卫星多径误差, 提高北斗导航系统的定位精度有着重要意义。
[1] | Wang G, de Jong K, Zhao Q, et al. Multipath Analysis of Code Measurements for Beidou Geostationary Satellites[J]. GPS Solutions, 2015, 19(1): 129-139. DOI:10.1007/s10291-014-0374-8 |
[2] | Wang M, Chai H, Liu J, et al. BDS Relative Static Positioning over Long Baseline Improved by GEO Multipath Mitigation[J]. Advances in Space Research, 2016, 57(3): 782-793. DOI:10.1016/j.asr.2015.11.032 |
[3] | Zhang F, He H, Tang B, et al. Analysis of Signal Characteristic and Positioning Performance Affected by Pseudorange Multipath for COMPASS[C]//China Satellite Navigation Conference(CSNC) 2013 Proceedings, 2013 |
[4] | Ma X, Shen Y. Multipath Error Analysis of COMPASS Triple Frequency Observations[J]. Positioning, 2015, 5(1): 12-21. |
[5] | Hauschild A, Montenbruck O, Thoelert S, et al. A Multi-Technique Approach for Characterizing the SVN49 Signal Anomaly, Part 1:Receiver Tracking and IQ Constellation[J]. GPS Solutions, 2012, 16(1): 19-28. DOI:10.1007/s10291-011-0203-2 |
[6] | Montenbruck O, Schmid R, Mercier F, et al. GNSS Satellite Geometry and Attitude Model[J]. Advances in Space Research, 2015, 56(6): 1015-1029. DOI:10.1016/j.asr.2015.06.019 |
[7] | Kabashima Y, Obuchi T, Uemura M. Approximate Cross-Validation Formula for Bayesian linear Regression[C]//201654th Annual Auerton Corterance on Communication, Control, and Computing, 2016, 27-30 |
[8] | Ma L, Li M, Gao Y, et al. A Novel Wrapper Approach for Feature Selection in Object-Based Image Classification Using Polygon-Based Cross-Validation[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(3): 409-413. DOI:10.1109/LGRS.2016.2645710 |
[9] | Helal M A, Haydar M S, Mostafa S. Algorithms Efficiency Measurement on Imbalanced Data using Geometric Mean and Cross Validation[C]//2016 International Workshop on Computtional Intelligence(IWCI), 2016:12-13 |
[10] | Grinsted A, Moore J C, Jevrejeva S. Application of the Cross Wavelet Transform and Wavelet Coherence to Geophysical Time Series[J]. Nonlin Processes Geophys, 2004, 11(5/6): 561-566. DOI:10.5194/npg-11-561-2004 |