基于值域的多约束GNSS单频单历元定姿新算法

吴洪涛，赵修斌，庞春雷，张良，冯波

空军工程大学 信息与导航学院，西安 710077

单频单历元姿态测量由于所具有实时性高、低成本、不受整周跳变影响的优势,已成为当前航空工程领域的研究热点,但单频单历元解算模型强度太弱，整周模糊度的难度较大。利用多天线可以显著增强单频单历元解算模型的强度[1]，实现载体航向角、俯仰角和横滚角的全姿态测量[2]。Giorgi和Teunissen利用坐标旋转矩阵建立多天线定姿模型，将LAMBDA算法扩展为多基线多约束LAMBDA (MC-LAMBDA)算法[3-5]。龚昂将利用解析方法求解最优旋转矩阵，进一步简化了MC-LAMBDA算法，提高了算法效率[6]。李青松将坐标旋转矩阵的单位正交约束条件融入了模糊度浮点解和协方差矩阵的求解过程，提高了单历元固定的成功率和正确率[7]。这些方法都是基于LAMBDA算法从模糊度域确定搜索空间，解算效率随模糊度维数增加快速下降，计算量随天线数量增加显著增长。

除了上述算法外，模糊度解算算法还有部分模糊度法[8-11]，长短基线法[12-13]、模糊度函数法[14-15]等。部分模糊度法基于分组的思想，将模糊度分为基本组和剩余组，直接对基本组模糊度进行求解，但分步求解降低了理论成功率，且一般很难解析给出模糊度分组算法。长短基线法利用无模糊度短基线辅助长基线解算模糊度，计算效率高，但成倍地增加了硬件成本。模糊度函数法通过将多维模糊度的搜索转化为三维坐标域的搜索，计算效率与模糊度维数无关，但在多天线条件下建立坐标域搜索涉及到多基线之间的相对位置转换，实现较为复杂。

文献[16]提出的基于值域的模糊度搜索方法，在基线长度已知的条件下，将对模糊度域的搜索转化为对姿态域的搜索，方法简单，易于实现。首先，对多天线姿态测量算法进行了分析；然后，将姿态约束融入解算模型中，建立了约束条件下值域搜索模型，提高了模糊度单历元搜索效率，推导了基于最优姿态解近似估计的模糊度固定算法。该算法较标准迭代算法计算效率大大增加，且具有与标准算法相当的性能，能够实现模糊度的单频单历元可靠固定。

1 多天线定姿算法及分析

1.1 多天线定姿模型

基线矢量在载体坐标系和东北天(ENU)坐标系之间的转换关系为

a′=Ra

(1)

式中：a为载体坐标系下的基线矢量；a′为ENU坐标系下的基线矢量；R为由载体坐标系到ENU坐标系的坐标旋转矩阵。R可由姿态角表示为

R=

(2)

式中：Cξ=cosξ；Sξ=sinξ；θ、β和γ分别为载体的航向角、俯仰角和横滚角；θ∈[0,2π]、β∈[-π,π]、γ∈[-π,π]。设

(3)

R中的元素满足单位正交条件,即

(4)

理论上只要得到R前两行元素,即可解得全部姿态角。

(5)

β=-arcsinr31

(6)

(7)

当存在多条基线时，式(1)可写为矩阵形式，n+1个天线对应n条相互独立的基线。设天线平面与载体坐标系重合，则

(8)

A′=RA

(9)

式中:A为载体坐标系下的基线矢量矩阵；A′为ENU坐标系下的基线矢量矩阵。多天线双差载波相位观测方程可写为[17]

[Φ1Φ2…Φn]=

(10)

式中：G为双差视线向量矩阵；Φ1、Φ2、…、Φn分别为基线1、基线2、…、基线n的双差载波相位观测向量；N1、N2、…、Nn分别为基线对应的双差整周模糊度向量；λ为载波波长。

对式(10)进行向量化运算可得

(11)

式中：⊗为Kronecker积；vec(·)为向量化运算。同理可得到多天线双差伪距观测方程为

(12)

式中：P1、P2、…、Pn分别为基线1、基线2、…、基线n的双差伪距观测向量。

令

则式(11)和式(12)有[17]

(13)

P=BM

(14)

当得到正确的模糊度空间之后，利用最小二乘原理固定模糊度，代价函数为

(15)

式中：QΦ和QP分别为双差载波相位观测量和双差伪距观测量的协方差矩阵。

当成功固定整周模糊度N之后，可根据式(13)解出M,由式(5)～式(7)可得到待求解的姿态角。

1.2 算法分析

多天线增强了单频单历元解算模型的强度，但模糊度维数随天线数量成倍增加，从模糊度域确定搜索空间的计算量将显著增长。基于值域的模糊度搜索算法，将对模糊度域的搜索转化为对姿态域的搜索，搜索过程与模糊度维数无关，适用于多天线姿态测量。但该算法没有将实际应用中的姿态约束信息融入解算模型，单历元模糊度搜索效率不高。一般通过对低成本惯性测量单元(Inertial Measurement Unit，IMU)输出数据进行解算，或水平仪、微机电系统传感器(Micro Electro Mechanical Systems，MEMS)等其他低成本传感器可以得到俯仰角、横滚角的强约束[18-19]。

当得到正确的模糊度空间之后，利用总体最小二乘原理固定模糊度。考虑俯仰角、横滚角的约束信息，代价函数可写为[20]

(16)

利用纯载波相位方程求解N条件下M的最小二乘解，往往不能满足单位正交约束条件。理论上需要在单位正交约束条件下，迭代求解N条件下坐标旋转矩阵M的最优解，但迭代求解会大大增加计算量。为了减小计算量，可以通过求解M最优解的近似解[7]，固定模糊度。近似解一方面需要具有与迭代解相当的性能，另一方面需要尽量减小计算量。

2 本文所提算法

2.1 多约束多天线模糊度搜索模型的建立

由式(2)M可以由θ、β、γ表示，即

M=M(θ,β,γ)

代入式(13)整理得

(17)

式中：下标i为第i条基线;j为第i条基线对应的第j个双差观测方程。Φij的误差通常在1/4周以内，若已知θ、β、γ,通过直接四舍五入取整即可确定整周模糊度:

(18)

式中：round(·)为四舍五入取整。整周模糊度可以看做姿态角的函数。基于此，通过直接对姿态角进行遍历搜索，解算出相应的模糊度浮点解，取整即可得到模糊度的搜索空间。

为了提高搜索效率，合理选择搜索步长成为关键。为了确保整周模糊度不被漏搜，并考虑到由于取整时四舍五入导致的模糊度变化，要求搜索步长使模糊度的每次变化不能超过0.5周。以航向角为例，设Δθ为航向角的搜索步长，即

(19)

当Δθ小于10°有sin(Δθ)≈Δθ,cos(Δθ)≈1代入式(19)整理得

i=1,2,…,n;j=1,2,…,t

(20)

式中：n为相互独立的基线数；t为双差观测方程个数。

文献[16]通过式(20)确定搜索步长时不考虑俯仰角、横滚角约束，而将约束信息融入式(20)中能进一步增大搜索步长，提高搜索效率。

2.2 姿态角约束条件下搜索步长的确定

设

(21)

(22)

(23)

由范数理论

(24)

设

则航向角搜索步长Δθ为

(25)

同理可得俯仰角、横滚角的搜索步长Δβ、Δγ为

(26)

(27)

其中类似式(24)推导过程

(28)

(29)

(30)

(31)

(32)

(33)

2.3 姿态解近似估计的模糊度固定

利用欧拉角对式(16)重新参数化后,代价函数为

(34)

式中：

L=F(N,η)

(35)

(36)

利用最小二乘原理有

Δη≈(KT(η0)QK(η0))-1KT(η0)QΔL

(37)

式中：ΔL=L-F(N,η0)。

(38)

所提算法流程图如图1所示。

图1 所提算法流程图Fig.1 Flow chart of proposed algorithm

3 实验结果及分析

为验证算法效果，分别进行了静态验证实验和动态实验。

3.1 静态实验条件及结果分析

分别进行了两组静态实验，利用GPS-703-GGG 天线和NovAtel OEM628板卡采集GPS系统L1频点数据，采样频率为1 Hz，实验地点为空军工程大学信息与导航学院科研楼顶，数据采用MATLAB2014A进行事后处理，电脑硬件配置为Intel(R) Core(TM) i3-2350M双核CPU，主频2.3 GHz，内存4 G。

第1组静态实验天线配置如图2所示，共采集63 488个历元的数据。基线在载体系下的基线矢量矩阵为

(39)

根据事先测定，载体在ENU坐标系中的航向角为277.424°，俯仰角为2.459°，横滚角为1.814°。实验中，给俯仰角、横滚角加入均方差为1°的零均值高斯噪声，模拟由其他传感器得到的俯仰角、横滚角信息。俯仰角、横滚角的搜索范围为±5°，航向角的搜索范围为±180°。

图2 第1组静态实验天线配置Fig.2 Configuration of antennas of first static experiment

第2组静态实验天线配置如图3所示，共采集1 716个历元的数据，基线在载体系下的基线矢量矩阵为

图3 第2组静态实验天线配置Fig.3 Configuration of antennas of second static experiment

(40)

实验中，利用俯仰角、横角滚先验信息由NovAtel SPAN-CPT惯性导航设备粗对准获得(标准差设为0.6°)。俯仰角、横滚角的搜索范围为±3°，航向角的搜索范围仍为±180°

3.1.1 模糊度搜索算法比较

为了对本文模糊度搜索算法与原算法进行比较，采用以下两种算法对数据进行解算，单历元确定模糊度搜索空间。

算法1原算法，即不考虑姿态角约束推导搜索步长确定模糊度搜索空间。

算法2本文所提算法，即利用姿态角约束推导步长确定模糊度搜索空间。

利用两种算法对两组静态实验数据进行解算，逐历元搜索步长和模糊度空间大小如图4所示，图4(a)为第1组实验结果，图4(b)为第2组实验结果。虚线代表算法1，实线代表算法2，平均搜索步长和模糊度空间平均大小如表1所示。由于第2组实验基线长度大于第1组实验，因此第2组实验姿态角搜索步长明显减小，模糊度空间也显著增加。但两组实验结果中算法2均优于算法1。

第1组实验中算法2与算法1相比，航向角平均搜索步长增大了34.5%，俯仰角平均搜索步长增大了150.6%，横滚角平均搜索步长增大了126.9%，模糊度空间减小了65.8%。

第2组实验中算法2与算法1相比，航向角平均搜索步长增大了41.3%，俯仰角平均搜索步长增大了59.1%，横滚角平均搜索步长增大了293.4%，模糊度空间减小了67.7%。

实验结果表明，本文所提算法充分利用了姿态角约束条件，单历元模糊度搜索效率明显提高。

图4 两种模糊度搜索算法实验结果Fig.4 Experimental results of two ambiguity search algorithms

表1 两种模糊度搜索算法实验结果对比Table 1 Comparison of experimental results of two ambiguity search algorithms

3.1.2 模糊度固定算法比较

根据事先精确测量的姿态角，逐历元反解得到正确的模糊度。在利用约束条件下的值域搜索模型得到正确的模糊度空间，采用以下4种方案对两组静态实验数据进行实验。

方案3本文所提算法,即基于最优条件姿态解近似估计，计算代价函数式(17)。

方案4标准迭代算法,即利用牛顿迭代法求解最优姿态解，计算代价函数式(17)。

4种方案计算的正确模糊度对应的代价函数值如图5所示，图5(a)为第1组前500历元实验结果,图5(b)为第2组前500历元实验结果。对4种方案计算的正确模糊度的平均代价函数值和固定模糊度的平均耗时进行统计，结果如表2所示。

图5 4种方案计算的前500历元正确模糊度对应的代价函数值Fig.5 Cost function value calculated from correct ambiguity for the first 500 epoch data by four schemes

表24种方案计算的正确模糊度平均代价函数值和固定模糊度的平均耗时
Table2Meanofcostfunctionvaluecalculatedfromcorrectambiguityandcomputation-timetofixambiguityoffourschemes

实验结果方案1方案2方案3方案4第1组实验正确模糊度的代价函数值0．0532．0170．2090．196固定模糊度平均耗时/s0．0070．0060．0450．226第2组实验正确模糊度的代价函数值0．0914．9332．6532．078固定模糊度平均耗时/s0．0470．0640．3141．301

方案1不考虑角度约束信息(式(16)中不包含角度约束信息加权部分)，因此得到正确模糊度的代价函数值最小。方案2没有充分利用角度约束信息，计算得到正确模糊度的代价函数值远大于方案4。方案3充分考虑了角度约束信息，且利用最优条件姿态角近似估计简化了迭代过程。因此，方案3与方案4相比，正确模糊度的代价函数值几乎一致，且计算时间大大减少。第2组实验中，方案4的模糊度固定时间大于1 s，几乎不能实时解算，而方案3的计算时间远远小于1 s。

4种方案单历元固定模糊度正确率，如表3所示。方案1仅考虑的基线约束，没有考虑角度约束信息，固定模糊度正确率最低。方案2计算最优条件解时利用角度约束条件不充分，固定模糊度正确率低于方案3。方案3与方案4由于充分利用了约束条件，固定模糊度正确率最高。

利用第1组实验数据，对4种方案在不同可视卫星条件下单历元固定模糊度的正确率进行了统计，如表4所示。

方案1、方案2固定模糊度的正确率随着可视卫星数量的减少而降低，方案3由于充分利用了约束条件，固定模糊度的正确率不受影响，与方案4相同。

实验结果表明，利用俯仰角、横滚角约束条件能够进一步提高模糊度固定的成功率和可靠性。本文所提算法相比直接法牺牲了一定的效率，但却有与标准算法相当的良好性能，且计算效率远远高于标准算法。

表3 4种方案固定模糊度正确率Table 3 Fixing ambiguity success rate of four schemes

表4不同可视卫星数量下4种方案固定模糊度正确率

Table4Fixingambiguitysuccessrateoffourschemesonconditionofdifferentvisiblesatellites

3.2 动态实验条件及结果分析

动态实验时，将3个天线固定在实验车顶，采集GPS系统L1频点的数据，采样频率为1 Hz。将NovAtel SPAN-CPT惯性导航设备水平固定在实验车的中轴线上，以初始化后的组合输出作为姿态角参考值(测姿精度俯仰角为0.02°，横滚角为0.02°，航向角为0.06°)，为了与静态实验保持一致。实验地点为西安市东郊，实验场景如图6所示，共采集896个历元的数据，其中前347 s实验车处于静止状态，用于标定3个天线的姿态测量系统与组合导航系统的安装误差。基线在载体坐标系下的基线矢量矩阵为

(41)

利用方案3对数据进行解算，为了与静态实验一致，俯仰角、横滚角先验信息的标准差仍设为0.6°，搜索范围为±3°，航向角的搜索范围为±180°。

采用方案3和方案4对数据进行单历元解算，确定姿态信息。

方案4不考虑姿态角约束推导搜索步长确定模糊度搜索空间，利用牛顿迭代法计算代价函数式(17)，单历元固定模糊度。

方案3考虑姿态角约束推导搜索步长确定模糊度搜索空间，利用最优条件姿态解近似估计计算代价函数式(17)，单历元固定模糊度。

图6 动态实验天线配置Fig.6 Configuration of antennas of dynamic experiment

两种方案得到模糊度空间大小和平均耗时如表5所示。方案4几乎不能实时输出。与方案4相比，方案3具有更小的搜索模糊度空间和更高的计算效率，能够做到实时输出。

以姿态角参考值作为理论真值，将解算得到俯仰角、方位角与真值对比，得到误差曲线如图7所示，姿态角均方误差(RMSE)如表6所示。

航向角的RMSE为0.170 2°，俯仰角的RMSE为0.390 2°，横滚角的RMSE为0.368 9°。

实验结果表明，利用本文所提算法解算的姿态结果外符合精度良好，姿态信息测量准确，模糊度得到了正确的固定。

表5 两种方案计算效率对比Table 5 Efficiency of two schemes

图7 动态实验测量结果Fig.7 Attitude results of dynamic experiment

表6 航向角、俯仰角和横滚角的 RMSETable 6 RMSE of heading, pitch and roll angles

4 结 论

1) 融入姿态约束的值域搜索模型相对不考虑姿态角约束，能够进一步增大搜索步长，提高搜索效率，减小模糊度空间。

2) 利用最优姿态解近似估计固定模糊度，具有与标准迭代算法相当的正确率，且计算效率远远高于标准迭代算法。

虽然利用俯仰角、横滚角约束能够提高值域搜索算法的效率，但当基线较长时，利用本文算法确定的模糊度空间仍较大，计算效率不高。从静态实验结果来看，在基线长度达到3.678 m时(第2组静态实验)，本文算法的实时性明显下降。因此使用本文算法进行单频单历元姿态测量时，基线长度一般不应大于4 m。下一步重点考虑如何有效利用航向角约束和伪距观测信息，进一步提高值域类算法在较长基线条件下的适用性。