一种基于蒙特卡洛采样的GNSS单频单历元姿态确定方法

摘要

一种基于蒙特卡洛采样的GNSS单频单历元姿态确定方法，其步骤如下：一：准备工作：主要是给出线性化的站际‑星际双差模型；二：建立多变量的GNSS姿态模型；三：用蒙特卡洛取样方法构造模糊度的概率分布函数；四：用LAMBDA算法搜索出模糊度的候选值；五：计算模糊度最优整数解并确定姿态；通过以上步骤，用蒙特卡洛方法求得了模糊度的期望和协方差，然后将其用于求解模糊度的LAMBDA方法中进行模糊度的解算，减少了计算模糊度的计算量，摆脱了对伪距测量量的依赖，从而可以简便快捷的确定姿态，减弱了对环境条件的要求。本发明在确定姿态的过程中，可以快速固定模糊度，从而可以快速确定姿态。

说明书

技术领域

本发明提供一种基于蒙特卡洛采样的GNSS的单频单历元姿态确定方法，它涉及一种基于蒙特卡洛采样的GNSS单频单历元模糊度固定和姿态估计方法，即一种对静止或者运动载体利用GNSS进行姿态确定的单频单历元载波模糊度固定和姿态估计方法，属于导航技术领域。

背景技术

全球导航卫星系统(GNSS,Global Navigation Satellite System)目前主要包括美国的全球定位系统(GPS)、俄罗斯的格洛纳斯系统(GLONASS)、中国的北斗导航系统(BeiDou)以及欧洲的伽利略系统(Galileo)。GNSS系统具有实时、高精度和全球连续覆盖等优势，可以为地面车辆、舰船、航空飞行器以及低轨卫星提供精密授时、定位和定姿等导航服务。近年来，基于GNSS的姿态确定技术在导航领域得到了广泛的关注。通过在载体平台上安装至少三个GNSS接收机，构建单差或双差观测方程，可以获得载体的完整姿态信息。GNSS姿态确定具有零漂移、低成本、功耗低等优点，对现有的定姿系统可实现很好的替代、完善或补充。

GNSS系统的基本观测量包含伪距和载波两种，二者的观测噪声分别在分米和毫米水平。高精度GNSS姿态确定主要依赖载波信号，尤其是对于中小尺寸(米级或分米级)的运动载体，比如中小型车辆、无人机、探空火箭、卫星等，必须借助载波信号。但是，载波信号存在整周模糊度问题。因此，如何快速准确的固定载波相位整数模糊度，是实现高精度GNSS姿态确定的核心和关键。根据模糊度解算过程中是否利用载体和GNSS卫星之间的相对运动信息，可以将姿态确定模糊度解算方法分为多历元法和单历元法两种。多历元法又称运动法，其基本原理是利用观测信息的时间累积获取模糊度浮点解和协方差矩阵，从而得到固定整数解。多历元法容易受到载波信号失锁和模糊度周跳的影响，不利于在复杂的导航环境中(信号遮蔽、高动态飞行等)使用。相比之下，单历元法可以不受信号失锁和周跳的影响。而且，单历元法本身的固定速度要快于多历元法。实现单历元模糊度固定，可以采取若干种途径，比如遍历搜索、多频信号组合和多接收机长-短基线配置等。对于低成本单频GNSS定姿，搜索法是最佳选择。其中，LAMBDA方法(最小二乘模糊度去相关调整，Least squares AMBiguity Decorrelation Adjustment)是目前效率最高、应用最广的一种方法。LAMBDA方法在本质上是一种整数最小二乘搜索法，通过Z变换使不同通道的模糊度的相关程度降到最低，同时仍保留其整数特性，因此具有更高的搜索效率和准确度。然而，应用LAMBDA方法固定模糊度需要同时利用伪距和载波测量。伪距的主要作用是缩小模糊度的搜索范围，而不是直接提供姿态解算。当伪距噪声较小时(5～30cm)，LAMBDA方法可以实现100％的模糊度固定成功率；当伪距噪声过大(>50cm)或者多径误差存在时，性能会随之下降。

本发明在LAMBDA方法的基础上提出一种基于蒙特卡洛采样的GNSS单频单历元模糊度固定和姿态确定的方法。该方法从模糊度概率分布的角度出发，利用事先已知的多天线几何分布信息对模糊度的搜索空间进行限定，从而摆脱对伪距信号的依赖，因此可以用于高伪距噪声(低端接收机、弱信号环境等)和高多径环境(城市街道、载体形状不规则等)的姿态确定。该方法适合于利用短基线和超短基线定姿的中小型(静止或运动)载体。

发明内容

(一)发明目的：本发明以LAMBDA固定模糊度方法为基础，针对短基线和超短基线姿态确定问题，提出了一种基于蒙特卡洛采样的GNSS单频单历元模糊度固定和姿态估计方法。

(二)技术方案

本发明一种基于蒙特卡洛采样的GNSS单频单历元姿态确定方法，其步骤如下：

步骤一：准备工作

首先，给出m+1个接收机对n+1个GNSS卫星的载波观测量给出站际-星际双差观测线性化方程如下所示：

式(1)中，Φ为m条基线向量的载波双差观测量组成的n×m阶矩阵，每一列向量表示一条基线对应的n个双差观测量；B为m条未知基线向量在GNSS参考坐标系下的坐标值组成的3×m阶矩阵；G是从卫星到接收机方向的双差单位向量在GNSS参考坐标系下的表示组成的n×3阶矩阵(由于各个副接收机之间距离很小，而卫星距离各副接收机很远，所以这里认为同一个卫星到任何一个副接收机方向的单位向量相同)；λ是载波波长；Z是双差整周模糊度组成的n×m阶矩阵，每一列向量表示一条基线对应的n个双差整周模糊度；V观测噪声组成的n×m阶矩阵；符号vec(·)表示把矩阵的列向量按列序号从小到大的顺序依次从上到下重新排列成一列，组成新的列向量；Q是vec(V)的协方差矩阵；取Q为：

其中σ是载波噪声的标准差；符号表示克罗内克积；

蒙特卡洛取样方法：表示概率分布函数p_x(x)的随机测量，是取样点，是各点对应的权重，并且N_s是取样次数，那么概率分布函数p_x(x)可近似表示为

$p_x\left(x\right)\approx \underset{i=1}{\overset{N_s}{\Sigma }}{w}_i\delta (x-x_i)---\left(3\right)$

其中，δ(·)表示狄拉克δ函数；

步骤二：建立多变量的GNSS姿态模型

定义本地坐标系：原点在主接收机位置，x轴沿着第一条基线方向，y轴在第一条基线和第二条基线决定的平面内并且垂直于x轴，z轴方向由右手定则决定；m条基线向量在本地坐标系中的表示已知，则由m条基线向量组成的3×m坐标矩阵F如下所示：

$\begin{array}{c}m\ge 3,F=\left(\begin{array}{ccccc}{f}_{11}& f_{21}& f_{31}& \mathrm{...}& f_{m1}\\ 0& f_{22}& f_{32}& \mathrm{...}& f_{m2}\\ 0& 0& f_{33}& \mathrm{...}& f_{m3}\end{array}\right)\\ m=2,F=\left(\begin{array}{cc}{f}_{11}& f_{21}\\ 0& f_{22}\\ 0& 0\end{array}\right)\end{array}---\left(4\right)$

其中，F中每一列向量表示一条基线向量在本地坐标系下的表示，即：在本地坐标系下，副接收机1的坐标为(f₁₁,0,0)，副接收机2的坐标为(f₂₁,f₂₂,0)，副接收机3的坐标为(f₃₁,f₃₂,f₃₃)，副接收机m的坐标为(f_m1,f_m2,f_m3)，其中各个坐标值的含义见说明书附图中的图2；

令R表示从本地坐标系到GNSS参考坐标系的坐标转移矩阵，则

B＝RF(5)

将其代入式(1)即可得到多变量的GNSS姿态模型：

其中，3阶标准正交方阵R和双差模糊度矩阵为未知量；将R中的参数用四元数表示，则(6)式可写为：

其中，q＝(q₀,q₄)^T，q₀＝(q₁,q₂,q₃)^T，并且满足q^Tq＝1，

$R\left(q\right)=\left(\begin{array}{ccc}{q}_1^2-q_2^2-q_3^2+q_4^2& 2(q_1{q}_2+q_3{q}_4)& 2(q_1{q}_3-q_2{q}_4)\\ 2(q_1{q}_2-q_3{q}_4)& -q_1^2+q_2^2-q_3^2+q_4^2& 2(q_2{q}_3+q_1{q}_4)\\ 2(q_1{q}_3+q_2{q}_4)& 2(q_2{q}_3-q_1{q}_4)& -q_1^2-q_2^2+q_3^2+q_4^2\end{array}\right)---\left(8\right)$

q为未知量；

綜上所述，所述的“建立多变量的GNSS姿态模型”，其建立的过程是将(4)式和(5)式代入(1)式得到(6)式，再将(8)式代入(6)式得到(7)式，(7)式即为建立的多变量GNSS姿态模型；

步骤三：蒙特卡洛取样方法构造模糊度的概率分布函数

根据(7)式，双差模糊度可以表示为：

$Z=\frac{1}{\lambda }(\Phi -GR(q)F-V)---\left(9\right)$

观测噪声遵从均值为0协方差为Q的正态分布；vec(V)的概率分布函数为：

$p_{vec\left(V\right)}\left(v\right)=n/(\nu ;\theta ,Q)\Delta =\frac{1}{\sqrt{2\pi \det Q}}{e}^{-/{\nu }^T{Q}^{-1}\nu }\left(10\right)$

如果已知四元数q有上界q^u和下界q^l，根据先验姿态信息，q遵从均匀分布：

$p_q\left(q\right)=R(q;q^{lu},q)\Delta ={\{}^{cq\in [q^{lu},q]}0q\notin [q^{lu},q]\left(11\right)$

其中，

$\underset{q\in [q^l,q^u]}{\int }cdq=1---\left(12\right)$

根据q^Tq＝1，可得

$p_q\left(q\right)=\{\begin{array}{cc}{c}^{\prime }& q_0\in [q_0^l,q_0^u],q_4=\pm \sqrt{1-{||q_0||}^2}\\ 0& q_0\notin [q_0^l,q_0^u],q_4=\pm \sqrt{1-{||q_0||}^2}\end{array}---\left(13\right)$

并且，

$\underset{||q_0||\le 1}{\underset{q_0\in [q_0^l,q_0^u]}{\int }}{c}^{\prime }{\mathrm{dq}}_0=1---\left(14\right)$

其中，是q^l的前三个分量，是q^u的前三个分量。在缺少先验姿态信息的情况下，可取q^l＝[1,1,1]^T，q^u＝[1,1,1]^T；根据概率分布函数p_q(q)和p_vec(V)(v)对q和v进行N_s次取样

$\begin{array}{c}{q}_i~p_q\left(q\right);v_i~p_{vec\left(V\right)}\left(v\right)\\ i=1,2,\mathrm{...},N_s\end{array}---\left(15\right)$

每次取样点对应的权重为

$w_i=\frac{1}{N_s},i=1,2,\mathrm{...},N_s---\left(16\right)$

根据(9)式可以得出模糊度的样本点则

$p_{vec\left(Z\right)}\left(z\right)\approx \underset{i=1}{\overset{N_s}{\Sigma }}{w}_i\delta (z-z_i)---\left(17\right)$

模糊度的期望和协方差为

$\bar{z}=\underset{i=1}{\overset{N_s}{\Sigma }}{w}_i{z}_i---\left(18\right)$

$P=\underset{i=1}{\overset{N_s}{\Sigma }}{w}_i(z_i-\bar{z}){(z_i-\bar{z})}^T---\left(19\right)$

这里得出的期望和方差将用于步骤四的LAMBDA算法中；

綜上所述，所述的“构造模糊度的概率分布函数”，是指先根据姿态转换矩阵R(q)中的四元数元素q和噪声的概率分布函数对q和噪声分别进行取样，由取样结果得到模糊度的样本，这样可以根据蒙特卡洛方法得到模糊度的概率分布函数；

步骤四：LAMBDA算法搜索出模糊度的候选值

由于LAMBDA算法是已有的现成方法，这里不做详细说明；通过LAMBDA算法可以搜索出N_c个模糊度的候选值

步骤五：计算模糊度最优整数解并确定姿态

由模糊度的候选值根据(9)式可以算得各候选值对应的姿态矩阵

${\hat{R}}_k=G^{+}(\Phi -\lambda {\hat{Z}}_k)F^{+},k=1,2,...,N_c---\left(20\right)$

X⁺表示矩阵X的伪逆。其中四元数的候选值可由下式得到

$\left(\begin{array}{c}{q}_0=\frac{1}{4q_4}\left(\begin{array}{c}{R}_{23}-R_{32}\\ R_{31}-R_{13}\\ R_{12}-R_{21}\end{array}\right)\\ q_4=\frac{1}{2}{(\mathrm{trR}+1)}^{\frac{1}{2}}\end{array}\right)---\left(21\right)$

其中{R_ij,i,j＝1,2,3}是矩阵R的元素；则模糊度的最终整数解和姿态矩阵为

$[\hat{q},\hat{Z}]=\underset{k=1,...,N_c}{min}||vec(\Phi -GR({\hat{q}}_k)F-\lambda {\hat{Z}}_k)||---\left(22\right)$

其中，在步骤四中所述的“LAMBDA算法搜索出模糊度的候选值”，其搜索的过程如下：

先不考虑模糊度的整数限制，直接用最小二乘法解出模糊度的浮点解，把上文得出的协方差矩阵构造转换矩阵Z，采用Z整数转换矩阵将模糊度和协方差矩阵进行转换，进行整周模糊度搜索，求得模糊度的整数解及其协方差矩阵，再对求得的模糊度和协方差矩阵进行Z逆变换，得到模糊度的整数解及其协方差矩阵；

通过以上步骤，用蒙特卡洛方法求得了模糊度的期望和协方差，然后将其用于求解模糊度的LAMBDA方法中进行模糊度的解算，减少了计算模糊度的计算量，摆脱了对伪距测量量的依赖，从而可以简便快捷的确定姿态，减弱了对环境条件的要求。

(三)优点

本发明提供的一种基于蒙特卡洛采样的GNSS的单频单历元姿态确定方法的优点在于：

①本发明对于基线的姿态确定和模糊度的解算不需要伪距测量量的辅助，很好地适用于在高伪距测量噪声，多路径条件下。

②本发明在确定姿态的过程中，计算模糊度方差的计算量小，从而减小总体计算量。

③本发明在确定姿态的过程中，可以快速固定模糊度，从而可以快速确定姿态。

附图说明

图1是本发明所述方法流程图。

图2是本地坐标系下各基线向量的示意图。

具体实施方式

本发明一种基于蒙特卡洛采样的GNSS的单频单历元姿态确定方法，见图1所示，其具体实施步骤如下：

步骤一：给出线性化的星际-站际双差模型

m+1个接收机对n+1个GNSS卫星的载波观测量的站际-星际双差观测线性化方程如下所示：

式(23)中，Φ为m条基线向量的载波双差观测量组成的n×m阶矩阵，每一列向量表示一条基线对应的n个双差观测量；B为m条未知基线向量在GNSS参考坐标系下的坐标值组成的3×m阶矩阵；G是从卫星到接收机方向的双差单位向量在GNSS参考坐标系下的表示组成的n×3阶矩阵(由于各个副接收机之间距离很小，而卫星距离各副接收机很远，所以这里认为同一个卫星到任何一个副接收机方向的单位向量相同)；λ是载波波长；Z是双差整周模糊度组成的n×m阶矩阵，每一列向量表示一条基线对应的n个双差整周模糊度；V观测噪声组成的n×m阶矩阵；符号vec(·)表示把矩阵的列向量按列序号从小到大的顺序依次从上到下重新排列成一列，组成新的列向量；Q是vec(V)的协方差矩阵。取Q为：

其中σ是载波噪声的标准差；符号表示克罗内克积。

步骤二：建立多变量的GNSS姿态模型

定义本地坐标系：原点在主接收机位置，x轴沿着第一条基线方向，y轴在第一条基线和第二条基线决定的平面内并且垂直于x轴，z轴方向由右手定则决定。m条基线向量在本地坐标系中的表示已知，则由m条基线向量组成的3×m坐标矩阵F如下所示：

$\begin{array}{c}m\ge 3,F=\left(\begin{array}{ccccc}{f}_{11}& f_{21}& f_{31}& \mathrm{...}& f_{m1}\\ 0& f_{22}& f_{32}& \mathrm{...}& f_{m2}\\ 0& 0& f_{33}& \mathrm{...}& f_{m3}\end{array}\right)\\ m=2,F=\left(\begin{array}{cc}{f}_{11}& f_{21}\\ 0& f_{22}\\ 0& 0\end{array}\right)\end{array}---\left(25\right)$

其中，F中每一列向量表示一条基线向量在本地坐标系下的表示，其中各个坐标值的含义见附图2。

令R表示从本地坐标系到GNSS参考坐标系的坐标转移矩阵，则B＝RF。那么，多变量的GNSS姿态模型可由式(23)得到：

其中，3阶标准正交方阵R和双差模糊度矩阵Z为未知量。将R中的参数用四元数表示，则(26)式可写为：

其中，q＝(q₀,q₄)^T，q₀＝(q₁,q₂,q₃)^T，并且满足q^Tq＝1，

$R\left(q\right)=\left(\begin{array}{ccc}{q}_1^2-q_2^2-q_3^2+q_4^2& 2(q_1{q}_2+q_3{q}_4)& 2(q_1{q}_3-q_2{q}_4)\\ 2(q_1{q}_2-q_3{q}_4)& -q_1^2+q_2^2-q_3^2+q_4^2& 2(q_2{q}_3+q_1{q}_4)\\ 2(q_1{q}_3+q_2{q}_4)& 2(q_2{q}_3-q_1{q}_4)& -q_1^2-q_2^2+q_3^2+q_4^2\end{array}\right)---\left(28\right)$

q为未知量。

根据(27)式，双差模糊度可以表示为：

$Z=\frac{1}{\lambda }(\Phi -GR(q)F-V)---\left(29\right)$

步骤三：蒙特卡洛取样方法构造模糊度的概率分布函数

观测噪声遵从均值为0协方差为Q的正态分布，则vec(V)的概率分布函数为：

如果已知四元数q有上界q^u和下界q^l，根据先验姿态信息，q遵从均匀分布：

$p_q\left(q\right)=R(q;q^{lu},q)\Delta ={\{}^{cq\in [q^{lu},q]}0q\notin [q^{lu},q]\left(31\right)$

其中，

$\underset{q\in [q^l,q^u]}{\int }cdq=1---\left(32\right)$

根据q^Tq＝1，可得

$p_q\left(q\right)=\{\begin{array}{cc}{c}^{\prime }& q_0\in [q_0^l,q_0^u],q_4=\pm \sqrt{1-{||q_0||}^2}\\ 0& q_0\notin [q_0^l,q_0^u],q_4=\pm \sqrt{1-{||q_0||}^2}\end{array}---\left(33\right)$

并且，

$\underset{||q_0||\le 1}{\underset{q_0\in [q_0^l,q_0^u]}{\int }}{c}^{\prime }{\mathrm{dq}}_0=1---\left(34\right)$

其中，是q^l的前三个分量，是q^u的前三个分量。在缺少先验姿态信息的情况下，可取q^l＝[1,1,1]^T，q^u＝[1,1,1]^T。

根据概率分布函数p_q(q)和p_vec(V)(v)对q和v进行N_s次取样

$\begin{array}{c}{q}_i~p_q\left(q\right);v_i~p_{vec\left(V\right)}\left(v\right)\\ i=1,2,\mathrm{...},N_s\end{array}---\left(35\right)$

每次取样点对应的权重为

$w_i=\frac{1}{N_s},i=1,2,...,N_s---\left(36\right)$

根据(29)式可以得出模糊度的样本点则

$p_{vec\left(Z\right)}\left(z\right)\approx \underset{i=1}{\overset{N_s}{\Sigma }}{w}_i\delta (z-z_i)---\left(37\right)$

模糊度的期望和协方差为

$\bar{z}=\underset{i=1}{\overset{N_s}{\Sigma }}{w}_i{z}_i---\left(38\right)$

$P=\underset{i=1}{\overset{N_s}{\Sigma }}{w}_i(z_i-\bar{z}){(z_i-\bar{z})}^T---\left(39\right)$

这里得出的期望和方差将用于步骤四的LAMBDA算法中。

步骤四：LAMBDA算法搜索出模糊度的候选值

由于LAMBDA算法是已有的现成方法，这里不做详细说明。通过LAMBDA算法可以搜索出N_c个模糊度的候选值

步骤五：计算各模糊度候选值对应的姿态矩阵并筛选出最优解

由模糊度的候选值根据(27)式可以算得各候选值对应的姿态矩阵

${\hat{R}}_k=G^{+}(\Phi -\lambda {\hat{Z}}_k)F^{+},k=1,2,...,N_c---\left(40\right)$

符号X⁺表示矩阵X的伪逆。其四元数的候选值可由下式得到

$\left(\begin{array}{c}{q}_0=\frac{1}{4q_4}\left(\begin{array}{c}{R}_{23}-R_{32}\\ R_{31}-R_{13}\\ R_{12}-R_{21}\end{array}\right)\\ q_4=\frac{1}{2}{(\mathrm{trR}+1)}^{\frac{1}{2}}\end{array}\right)---\left(41\right)$

其中{R_ij,i,j＝1,2,3}是矩阵R的元素。

模糊度的最终整数解和姿态矩阵为

$[\hat{q},\hat{Z}]=\underset{k=1,...,N_c}{min}||vec(\Phi -GR({\hat{q}}_k)F-\lambda {\hat{Z}}_k)||---\left(42\right)$