路 阳,何 伟,王 乐
(中国电子科技集团公司第二十研究所,西安 710068)
GBAS 在飞机起降的应用,属于高端应用,能有效地提高飞机起降引导精度,为飞机航行保驾护航。但在实际应用中,由于新增卫星的伪距精度差,采用最小二乘解算,会引起定位解算结果的变化,从而造成高度跳变。在本文中,采用有效的抗差估计方法,能有效地提高伪距收敛精度,缩短收敛时间,提升GBAS 差分定位在飞机航行引导中的可靠性。
在GBAS 差分定位中,地面站伪距差分修正量的正确与否直接影响到机载端接收机的定位精度。实际应用中差分量生成流程如图1所示。
图1 差分量生成流程
生成步骤如下:
(1)原始残差生成
一颗卫星(编号为i)通过星历外推可以得到t时刻卫星位置为(xi,yi,zi),差分卫星地面站基准接收机天线的位置为(xr,yr,zr),从差分卫星地面站接收机天线r到卫星i的几何距离(真距)rr(i)为:
对4 个接收机同时收星,卫星(编号为i)在t时刻的伪距测量值为ρr(i)(m),则含有残差的伪距差分修正量为ρc(i)(m):
式中,i代表卫星号,i=1,2,3…32;m=1,2,3,4,代表4个接收机号。
(2)原始校正差分量的生成
由于式(2)中计算公式所得的差分量含有本地钟差,虽然对机载差分定位没有影响,但为了减少空地传输字节量,一般采用如下公式,估算本地钟差并消除:
(3)融合差分量的生成
对4 台校正本地钟差后的差分修正量进行融合处理,生成最终的差分修正量:
机载接收机(编号为u)对共视卫星i的伪距测量值为ρu(i),则修正后的机载关于这颗星的伪距测量值为:
机载和地面有N颗共视卫星,机载接收机修正后的伪距定位方程为:
通过解算式(7)可以计算得到机载接收机的差分定位结果(x,y,z)。再通过坐标转换得到大地坐标系(Lon′ ,Lat′ ,Height′)。
当i号低仰角卫星刚刚跟踪上,伪距的跟踪精度还没稳定时,伪距测量值为ρr(i)中包含的有误差分量 Δρr(i),Δρr(i)精度收敛图如图2所示。在接收机刚跟踪上卫星时,产生的伪距误差达到3 m 以上,随着时间的推移,伪距误差逐渐减小,由实践经验可知通常在40 s 以后,能够稳定在厘米级。图3所示为地面伪距误差引起对应的机载端高度变化图。
图2 伪距精度收敛图
图3 机载端高度变化图
i号卫星伪距未达稳定时,对4 个接收机会产生错误的差分修正量ρc′(i)为:
采用上述方法估算本地钟差并消除:
对4 台校正本地钟差后的差分修正量进行融合处理,生成最终的差分修正量为:
机载接收机收到错误的差分修正量后,修正共视卫星i的伪距测量值,则得到带有偏差的修正后的伪距为:
机载和地面有N颗共视卫星,其中包含一颗含有错误伪距修正值的机载伪距值,机载接收机的伪距定位方程:
式(13)中,由于第一个方程式中的错误伪距,导致在解算方程组时得到误差较大的机载差分定位(x′,y′,z′),坐标转换得到错误的大地坐标系(Lon′ ,Lat′ ,Height′)。因此,当一颗低仰角卫星被接收机接收到,伪距测量值的跟踪精度还没到达稳定时,会产生错误的机载修正伪距,导致式(13)解算出错,产生高度跳变。从图2~图3也可以看出,只有当地面接收机的伪距跟踪精度稳定后,解算上述方程组才会得到正确的定位结果,从而机载的差分定位结果恢复到正确的定位值[1]。
由此看出,只有伪距精度达到稳定状态后,才能生成正确的伪距差分修正量。传统方案根据实践经验,当新增卫星出现60 s 后,再让其参与差分运算,但此方法效率不高,较为简单原始,缺乏理论依据。而本文采用的自适应抗差估计算法能够对新增卫星伪距精度作出灵活自主判断,将其应用到新增卫星伪距测量值的预处理中,使问题得到更好的解决。
2.1.1 传统最小二乘估计
设观测向量 L= [L1L2…Ln]T,Li(i=1,2,…,n)独立,未知参数向 X= [X1X2…Xm]T,其误差向量Δ,线性化后的观测方程为:
误差方程为:
经典最小二乘(LS)估计要求:
由式(17)可得:
式中,∑X为X 的协方差矩阵;为方差因子估值。
我们以伪距双差作为对伪距精度的估计,并假设双差值各历元间相互独立,由此可得式(18)~式(20)中权阵P 为n阶单位阵,系数矩阵A 是各元素值均为1 的n×1 维列向量,故:
2.1.2 极大似然性估计(M 估计)
设观测样本L1,L2,…,Ln,测值Li分布密度为f(X-Li),极大似然估计要求参数估值满足:
用ρ(·)代替 -Inf(·),则极大似然原理要求:
即M 估计转化为求上述极小化问题的解,对式(22)求导后可得:
式中,ψ是ρ的导函数。
对于M 估计,将式(25)改写,令:
将式(26)代入式(25)可得:
式中,iW取Tukey 的ψ函数所相应的值,即:
式中,C=6.0,
2.2.1 相关最小二乘估计
设各观测值L1,L2, …,Ln均服从正态分布,故其联合概率密度函数为:
式中,L 为n维测值向量;∑L为L 的协方差矩阵;为测值协方差矩阵的行列式。式(30)可简写为:
实际求解时真实误差 Δ = L -E(L)无法直接求得,故只能求测值残差 V = L -,其中为L的估计值,因此式(32)可转换为:
将式(15)代入式(34)后求导可得:
表示成矩阵形式,可得估计值[2]为:
2.2.2 相关IGG3 方案
相关IGG3 方案基于相关最小二乘估计,且根据测值的不同情况分别作出处理,即当测值偏差较小时(测值主部)采用效率高的LS 估计;当测值超出一定范围时采用降权估计;当个别测值明显过大时,采用淘汰法(零权估计)。由此可得如下相关等价权[3]:
式中,k0可取值为1.0~1.5,k1可取值为2.5~3.0。
伪距精度采用伪距双差进行估计,伪距双差变量X的近似真值为0 m,对其进行了10 次观测得到:L=[0.099268,0.071994,0.027277,0.099270,0.108213,0.072887,0.077807,0.061709,0.094797,0.028172]。将L的首个历元测值L(1)加入大约0.3 m 的粗差后变为L′,然后对L′分别用M 估计、相关最小二乘估计、相关IGG3 方案这三种抗差估计算法进行处理,并使各算法最终收敛时,每种抗差方法的最终收敛结果、迭代次数及最终收敛精度如表1所示。
表1 三种抗差算法收敛结果对照表
各抗差算法收敛结果对比图如图4所示。
图4 三种抗差算法收敛效果对比图
由以上试验结果可得如下结论:从最终收敛结果方面,相关IGG3 方案最优,其次为M 估计,最后为相关最小二乘估计;从迭代次数方面,相关IGG3 方案和相关最小二乘估计均为2 次,算法效率最高,而M 估计迭代5 次,效率最低;从最终收敛精度方面,IGG3 方案最大,说明其收敛速度最快,而M 估计的最终收敛精度最小,说明其收敛速度最慢,相关最小二乘估计收敛速度介于两者之间。
表2为某地面接收机静态下的连续28 个观测历元的伪距双差观测值[4]。从数据变化上看,各历元观测值围绕伪距双差的近似真值0 上下小幅波动,我们依旧将观测值首个历元加入大约0.3 m 的粗差后,分别采用三种抗差算法进行估计,试验结果如图5所示。
表2 接收机连续28个观测历元的伪距双差观测值
由图5可以明显得出如下结论:三种抗差算法均能较好剔除较大的伪距粗差(首个历元跳点被“过滤”掉),从估计值准确度方面,起始阶段具有独立性的M 估计相比于其他两个具有相关性的估计算法效果较好,但随着历元数不断增多,后期(从14 历元起)两个相关估计算法占优,具有与真值更为接近的估计,相比较而言,IGG3 方案比相关LS 更能有效剔除观测值明显异常的点,并对部分测值较高点采用降权估计,因此采用IGG3 方案更为合适;从估计算法的连续性方面,具有独立性的M 估计连续性最好,而两个相关估计算法由于需要运用到权阵P,但在某些连续历元相关性较差的时刻权阵P 无法生成,因此会导致这些时刻无法生成估计值,故连续性较差。
图5 三种方案抗差结果图
通过上述分析可知,任何算法都不会十全十美,正所谓“鱼与熊掌不可兼得”,必须在算法估值准确度和连续性之间做出取舍,因此实际应用中,将M 估计和IGG3 方案综合使用效果更好,即在连续历元相关性较好的时刻,权阵P 存在,故此时采用相关IGG3 方案,而在连续历元相关性较差的时刻,权阵P 不存在,故此时采用独立M 估计,具体试验结果如图6所示。
图6 M 与IGG3 组合抗差结果图
由图6可知,原始伪距双差在34 s 以后精度达到0.1 m 以下,而采用M 与IGG3 组合抗差估计算法后,伪距双差在13 s 以后就能达标,由此可见,该算法适用于接收机初始跟踪锁定阶段的伪距抗差校正,能极大缩短接收机伪距的收敛时间,并能适当提高伪距精度[5]。
通过以上分析可知,采用M 与IGG3 组合抗差估计算法并将其应用到新增卫星伪距测量值的预处理中,可极大提高接收机的收敛速度和精度。