本发明属于GNSS卫星导航定位和测绘科学技术领域，具体是指一种利用GNSS测速原理的双天线定向测姿的方法。
背景技术：
在GNSS卫星导航定位和测绘科学技术及其应用领域中，如车载、机载、船载导航定位、地下工程贯通测量以及武器系统机动发射等，均需要确定载体方位或姿态。传统的方法有磁罗盘、寻北陀螺仪、惯性导航系统。其中，磁罗盘实际测姿精度易受周边电磁场环境影响，寻北陀螺仪测姿受到地域限且其造价成本较高，惯性导航系统测姿误差会随时间积累且其造价非常高。卫星导航系统具有全球、全天候、抗干扰性能好的特点，广泛应用于授时、定位和测速。同时高精度GNSS设备迅速发展，使得以GNSS构建载体测姿系统成为全新的定向测姿技术途径。
目前，现有GNSS定向测姿系统均是通过测量天线间基线向量实现载体定向或测姿(简称为基线向量测姿)。因此，若需确定载体方位，至少需要安装两个GNSS天线；若需确定载体三个姿态角，至少需要安装三个GNSS天线。如：专利CN1632465A采用两个天线实现了卫星定向；专利CN104597471A采用任意布置同一平面的三个GNSS天线，并通过载波单差观测方差估算载体姿态矩阵；专利CN102998690A将三个GNSS天线布置成直角三角形，利用载波双差方程估算测姿矩阵。需指出的是，现有定向测姿方法利用GNSS载波相对定位原理获得高精度基线向量，必须处理整周模糊度固定和周跳探测修复等棘手难题。而且，对于动态载体，模糊度固定和周跳探测修复的实时理论算法仍属于开放性研究难题。
技术实现要素：
为解决上述现有难题，本发明提出了利用GNSS测速原理的双天线定向测姿的方法，首先避免了模糊度固定和周跳探测修复难题，其次降低了三维测姿对GNSS天线数量的要求，其三同时适用于静态和动态载体的定向测姿需求。本发明同时适用于静态和动态载体的定向(二维测姿)和定姿(三维测姿)，且天线数量和装置成本均较低。
本发明采取的技术方案如下：步骤一、在载体上安装GNSS双天线定向测姿装置，所述双天线定向测姿装置包括GNSS双天线、接收机、电动导轨、滑块、步进电机、运动控制器，所述GNSS双天线包括天线a和天线b，天线a和电动导轨设于载体上，天线b通过滑块固设于电动导轨上，调节天线a和天线b高度使之与电动导轨平面平行，运动控制器与天线b连接，通过运动控制器使天线b在电动导轨上运动；
步骤二、通过通讯数据线连接双天线定向测姿装置和计算机，实时接收GNSS载波、伪距和多普勒数据；
步骤三、GNSS载体位置伪距测量；
记GNSS天线a、天线b和第j号卫星在t时刻的伪距观测值分别记为则固定于载体的天线a、b的伪距观测方程分别为：
式中，表示天线a到卫星j的几何距离；δta(t)为天线a钟差，δtj(t)为卫星j钟差；为相对论效应误差；为电离层和对流层误差，为多路径效应等其他误差；C表示光速；
忽略GNSS天线标记，将伪距方程中的几何距离在天线近似坐标(x0，ｙ0，z0)处进行泰勒级数展开，取一次项，可得：
式中，为t时刻天线近似位置至卫星j方向上的方向余弦，(x0，ｙ0，z0)表示天线近似三维坐标，(Xj，Yj，Zj)为卫星j坐标，表示天线到卫星j的近似距离，dx(t)＝(δx，δy，δz)T为天线坐标改正值；
由此，可得线性化伪距定位观测方程：
式中，
当同步观测n颗卫星，n大于四颗，将上式写成矩阵形式并由最小二乘估计可得GNSS天线三维空间直角坐标改正数：
式中，表示伪距观测值的权阵；
通过空间直角坐标与大地坐标转换公式，将天线空间直角坐标(x，y，z)转换成大地坐标(B，L，H)；由于测姿应用中两个天线距离一般很小，采用GNSS天线a或b得到的经纬度(Ba，La)或(Bb，Lb)均可用来表示载体位置；空间直角坐标与大地坐标迭代变换公式如下：
式中，N为该点的卯酉圈曲率半径；e分别为椭球第一偏心率；
步骤四、GNSS天线多普勒测速；
记天线a、天线b和第j号卫星t时刻的多普勒观测值分别为和设固定在载体上天线a的速度为安装在导轨上的天线b的速度为则GNSS天线测速方程如下：
式中，表示t时刻卫星j的多普勒频移观测值；为伪距变化率；为天线a钟差变化率；为卫星j的钟差变化率；为卫星j与天线a相对论效应变化率；C为光速；为其他观测值噪声变化率等；
忽略天线下标，伪距变化率进行一阶泰勒级数展开，可表示为：
式中，表示卫星j的速度；表示天线a或天线b的速度；
由此，可得GNSS天线多普勒测速线性化方程：
式中，
当同步观测n颗卫星(大于四颗)，写成矩阵形式并由最小二乘估计可得GNSS天线速度估计公式：
式中，表示多普勒观测值的权阵；因此，通过本步骤可以获得天线a和天线b的速度
步骤五、GNSS天线站心系速度换算；
由于天线b速度是载体速度和相对导轨速度的合成，天线b导轨速度记为则天线b导轨速度为将天线a速度和天线b轨速换算到站心系速度：
式中，Na(t)，Ea(t)，Ua(t)分别表示天线a在站心系下三个方向速度；Norb(t)，Eorb(t)，Uorb(t)分别表示天线b导轨在站心系下三个方向速度；
或者
步骤六、动态载体三维姿态实时估计；
根据天线a载体的站心速度和天线b导轨的站心系速度可构建载体姿态矩阵，为姿态矩阵第一列向量，为姿态矩阵第二列向量，姿态矩阵第二列向量与第一列向量叉乘作为姿态矩阵第三列向量，从而得到动态载体的三维姿态估计公式：
式中，×表示向量外积，右箭头上标表示列向量单位化，θ为横滚角，h为俯仰角，α为方位角。
采用上述方案本发明取得的有益效果如下：本方案利用GNSS测速原理的双天线定向测姿的方法，与现有GNSS定向测姿方法相比，首先，本发明只需要GNSS双天线即可实现动态三维姿态测量，减少了测姿系统对GNSS天线数的要求，因而降低了系统成本。其次，本发明除了适用于动态载体的三维测姿需求，还适用于静态情况下的俯仰角和方位角的静态定向(二维测姿)需求，而且此时仅仅需要一个GNSS天线。再次，该发明避免了载波相位模糊度固定和周跳探测修复难题，适用于高动态载体定向测姿需求。
附图说明
图1为本发明方法流程图；
图2为直线往返式GNSS双天线定向测姿装置图；
图3为圆弧旋转式GNSS双天线定向测姿装置图；
图4为本发明实施例计算得到的方位、俯仰和横滚角结果图。
其中，①GNSS天线a，②GNSS天线b，③步进电机，④导轨滑块，⑤固定基座，⑥控制器与数据存储处理终端，⑦、⑧、⑨为通讯数据线，⑩直线导轨，载体连接螺栓，圆弧导轨。
具体实施方式
下面结合具体实施对本专利的技术方案作进一步详细地说明，本发明所述的技术特征或连接关系没有进行详细描述的部分均为采用的现有技术。
以下结合附图，对本发明做进一步详细说明。
如图1所示，一种利用GNSS测速原理的双天线定向测姿装置及方法，包括以下步骤：
第一步、在载体上安装GNSS双天线定向测姿装置。该装置包括以下主要部件：GNSS双天线、接收机、电动导轨、滑块、步进电机、运动控制器等。天线a和电动导轨直接安装在载体上，天线b通过滑块固连在电动导轨上，调节天线a、b高度使之与电动导轨平面平行，通过运动控制器使天线b在电动导轨上运动。
第二步、通过通讯数据线连接双天线定向测姿装置和计算机，实时接收GNSS载波、伪距和多普勒数据。
第三步、GNSS载体位置伪距测量
记GNSS天线a、b和第j号卫星在t时刻的伪距观测值分别记为则固定于载体的天线a、b的伪距观测方程分别为：
式中，表示天线a到卫星j的几何距离；δta(t)为天线a钟差，δtj(t)为卫星j钟差；为相对论效应误差；为电离层和对流层误差，为多路径效应等其他误差；C表示光速。
忽略GNSS天线标记，将伪距方程中的几何距离在天线近似坐标(x0，ｙ0，z0)处进行泰勒级数展开，取一次项，可得：
式中，为t时刻天线近似位置至卫星j方向上的方向余弦，(x0，ｙ0，z0)表示天线近似三维坐标，(Xj，Yj，Zj)为卫星j坐标，表示天线到卫星j的近似距离，dx(t)＝(δx，δy，δz)T为天线坐标改正值。
由此，可得线性化伪距定位观测方程：
式中，
当同步观测n颗卫星(大于四颗)，将上式写成矩阵形式并由最小二乘估计可得GNSS天线三维空间直角坐标改正数：
式中，表示伪距观测值的权阵。
通过空间直角坐标与大地坐标转换公式，将天线空间直角坐标(x，y，z)转换成大地坐标(B，L，H)。由于测姿应用中两个天线距离一般很小，采用GNSS天线a或b得到的经纬度(Ba，La)或(Bb，Lb)均可用来表示载体位置。空间直角坐标与大地坐标迭代变换公式如下：
式中，N为该点的卯酉圈曲率半径；e分别为椭球第一偏心率。
第四步、GNSS天线多普勒测速
记天线a、b和第j号卫星t时刻的多普勒观测值分别为和设固定在载体上天线a的速度为安装在导轨上的天线b的速度为则GNSS天线测速方程如下：
式中，表示t时刻卫星j的多普勒频移观测值；为伪距变化率；为天线a钟差变化率；为卫星j的钟差变化率；为卫星j与天线a相对论效应变化率；C为光速；为其他观测值噪声变化率等。
忽略天线下标，伪距变化率进行一阶泰勒级数展开，可表示为：
式中，表示卫星j的速度；表示天线a或b的速度。
由此，可得GNSS天线多普勒测速线性化方程：
式中，
当同步观测n颗卫星(大于四颗)，写成矩阵形式并由最小二乘估计可得GNSS天线速度估计公式：
式中，表示多普勒观测值的权阵。
因此，通过本步骤可以获得天线a和天线b的速度
第五步、GNSS天线站心系速度换算
由于天线b速度是载体速度和相对导轨速度的合成，天线b导轨速度记为则天线b导轨速度为将天线a速度和天线b轨速换算到站心系速度：
式中，Na(t)，Ea(t)，Ua(t)分别表示天线a在站心系下三个方向速度；Norb(t)，Eorb(t)，Uorb(t)分别表示天线b导轨在站心系下三个方向速度。
或者
第六步、动态载体三维姿态实时估计
根据天线a载体的站心速度和天线b导轨的站心系速度可构建载体姿态矩阵，为姿态矩阵第一列向量，为姿态矩阵第二列向量，姿态矩阵第二列向量与第一列向量叉乘作为姿态矩阵第三列向量，从而得到动态载体的三维姿态估计公式：
式中，×表示向量外积，右箭头上标表示列向量单位化，θ为横滚角，h为俯仰角，α为方位角。
本发明所述的电动导轨可以采用多种形式，两种典型形式为直线往返式(如图2所示)和圆弧旋转式(如图3所示)电动导轨。其中，直线往返式要求控制器具有自动往返功能，圆弧旋转式则无此要求。其次，当滑块为精密数控滑块，精确获得GNSS天线b(精密数控滑块)相对导轨速率数据后，则利用附约束参数估计方法可进一步提高定向测姿精度。再次，所述的双天线可以采用独立时钟双天线或者共时钟双天线技术、所述的多普勒可以采用原始多普勒观测值或者导出多普勒观测值、所述的GNSS可采用GPS、GALILEO、GLONSS、BDS或者组合技术。最后，本发明适用于静态载体定向(二维测姿)，也适用于动态载体测姿(三维测姿)。
以直线往返式电动导轨、独立时钟GPS双天线多普勒观测值的动态车载三维测姿为实施例(华测X91，采样率5Hz，载体平均时速约50km/h)，图4为本发明实施例计算得到的方位、俯仰和横滚角结果图。参见图4可知，车辆在略微下坡道路上基本沿直线行驶，所得方位角约47.75°、俯仰和横滚角都很小，与动态载体行驶过程中的姿态情况相一致。
以上对本发明及其实施方式进行了描述，这种描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构并不局限于此。总而言之，如果本领域的技术人员受其启示，在不脱离本发明创造宗旨的情况下能够想到的变化和优点都被包含在本发明中，均应属于本发明的保护范围。