发明内容
发明目的:针对北斗二代特点,克服传统单频整周模糊度(LAMBDA)法实时性不高的缺陷和周跳探测和修复灵敏度不高的缺陷,本发明提供一种短基线下惯导辅助的北斗单频整周模糊度求解方法。
技术方案:为实现上述目的,本发明采用的技术方案为:
一种短基线下惯导辅助的北斗单频整周模糊度求解方法,包括如下步骤:
步骤一、信息初始化,采集北斗原始观测量并计算北斗双差观测量,为有效消除北斗双差观测量的大气延时误差,将北斗***输出的北斗原始观测量按照北斗星座类型的轨道高度分组差分,从而得到北斗双差载波相位观测值以及北斗卫星到接收机的双差单位矢量,构造北斗双差载波观测方程;
步骤二、利用惯导输出姿态矩阵、平台误差角信息以及天线配置信息估算基线矢量地心地固坐标(ECEF坐标)以及相应的基线矢量方差阵矩阵;
步骤三、将步骤二得到的惯导估算基线矢量代入步骤一中的北斗双差载波观测方程,利用递推最小二乘法求解整周模糊度浮点解及其协方差矩阵;
步骤四、根据步骤三得到的整周模糊度浮点解及其协方差矩阵利用改进最小二乘模糊度去相关方法固定整周模糊度整数解;
步骤五、利用惯导估算基线矢量结合北斗三差载波观测值一起构造检验量判断周跳是否发生,若发生周跳则利用检验量取整的方式估算周跳值并对整周模糊度进行修复。
所述步骤一中,为有效消除北斗双差观测量的大气延时误差,将北斗***输出的北斗原始观测量按照北斗星座类型的高度进行分组差分:MEO卫星一组,卫星数量为n,取仰角最大的MEO卫星为基准;IGSO卫星和GEO卫星一组,卫星数量为m,取仰角最大的IGSO/GEO卫星作为基准;构造北斗双差载波观测方程如下:
其中,为第j颗MEO卫星到基准MEO卫星(在MEO卫星组中标记为第1颗MEO卫星)的双差载波观测相位,为第j颗IGSO卫星/GEO卫星到基准IGSO卫星/GEO卫星(在IGSO卫星和GEO卫星组中标记为第1颗IGSO卫星/GEO卫星)的双差载波观测相位,为MEO卫星到载体之间的单位矢量,为GEO卫星/IGSO卫星到载体之间的单位矢量,为第j颗MEO卫星到基准MEO卫星的双差整周模糊度,为第j颗IGSO卫星/GEO卫星到基准IGSO卫星/GEO卫星的双差整周模糊度,为第j颗MEO卫星到基准MEO卫星的双差观测噪声,为第j颗IGSO卫星/GEO卫星到基准IGSO卫星/GEO卫星的双差观测噪声,λ为北斗B1频率信号的载波波长,为基线矢量。
所述步骤二中,利用惯导输出姿态矩阵、平台误差角信息以及天线配置信息估算基线矢量地心地固坐标及相应的基线矢量方差矩阵为:
其中,XI为惯导估算基线矢量,XG为北斗天线构成的基线在机体坐标系内的向量;为理想姿态矩阵,为惯导实际输出的机体坐标系到地理坐标系的姿态矩阵,为地理坐标系到地心地固坐标系的转换矩阵,φ为平台误差角向量,cov[φφT]为平台误差角方差阵,I为单位矩阵,cov[XI]为惯导估算基线矢量的协方差矩阵。
所述步骤四中,根据步骤三得到的整周模糊度浮点解及其协方差矩阵利用改进最小二乘模糊度去相关法固定整周模糊度整数解,包括:
(41)将整周模糊度浮点解协方差矩阵中的正对角线元素先进行顺序排序,再对排序后的协方差矩阵进行LDLT分解,得到整周模糊度去相关变换矩阵;
(42)利用整周模糊度去相关变换矩阵使整周模糊度去相关,然后结合整周模糊度的搜索条件固定整周模糊度整数解。
所述步骤四中,整周模糊度去相关变换矩阵获取方法如下:
(411)确定排序矩阵P,首先将矩阵P设定为与矩阵维数相同的零矩阵,然后将整周模糊度协方差矩阵中的正对角线元素Qii按升序排列,若Qii在排序后位于矩阵对角线元素的第m位,则令矩阵P中的元素Pmi=1,否则为零,计算得到矩阵P中所有元素的值后,计算排序后的矩阵其中i=1、2、…,m=1、2、…;
(412)然后对排序后的矩阵QP进行LDLT整数分解,得到下三角矩阵L后代入到Q=[L]-1QP[LT]-1中进行计算;
(413)基于(411)和(412)两个步骤,不断地进行迭代,直到矩阵[L]-1为或近似为单位阵结束;最终可以得到整周模糊度去相关变换矩阵为:Z=[L1]-1P1[L2]-1P2…[Lk]-1Pk。
所述步骤五中,利用惯导估算基线矢量集合北斗三差载波观测量一起构造检验量判断周跳是否发生,若发生周跳则利用检验量取整的方式估算周跳值并对整周模糊度进行修复,其中构造的检验量为:
T=(B(tk)·△b(tk+1,tk)-△y(tk+1,tk))/λ
其中,△y(tk+1,tk)=y(tk+1)-y(tk),y(tk)为tk时刻的北斗三差载波观测量;△b(tk+1,tk)=b(tk+1)-b(tk),b(tk)和b(tk+1)用惯导估算基线矢量代替, B(tk)为tk时刻对应分组的设计矩阵,或
若检验量T的绝对值大于1/2,则发生周跳,用检验量T四舍五入取整来计算发生周跳的数目。
有益效果:本发明提供的短基线下惯导辅助的北斗单频整周模糊度求解方法,针对北斗导航***特点将北斗原始观测量按照北斗星座类型的轨道高度进行分组差分,距离地面高度相近的卫星星座作为一组,能够有效消除北斗双差观测量中的大气延时误差,并且加入惯导辅助的方式,利用惯导输出姿态阵估算基线矢量信息辅助求解整周模糊度浮点解,求解北斗双差载波观测方程的方式为递推最小二乘法,相比于传统的最小二乘法模糊度去相关方法,可以提高整周模糊度的解算速率和精度;此外,利用惯导估算基线矢量和北斗三差载波观测值一起构造检验量进行周跳的探测和修复,相比于传统方法简单且十分灵敏;本方法适用于单频北斗导航***对高动态载体姿态进行定位或定姿。
具体实施方式
下面结合附图对本发明作更进一步的说明。
如图1所示为一种单频北斗导航***固定整周模糊度的求解方法,方法整体的流程如图1所示。首先,将北斗接收机输出的北斗原始观测量按照北斗星座类型的轨道高度进行分组差分得到北斗双差观测量和北斗双差载波观测方程,同时存储惯导输出姿态矩阵。接着,进行整周模糊度整数解的固定:首先利用惯导输出姿态矩阵估算基线矢量的地心地固坐标以及相应的基线矢量方差矩阵,再将惯导估算基线矢量代入北斗双差载波观测方程并利用递推最小二乘法计算整周模糊度浮点解及其协方差矩阵,接着利用改进最小二乘模糊度去相关法固定整周模糊度整数解,并进行整周模糊度ratio值检验。最后,利用惯导估算基线矢量以及北斗三差载波观测量一起构造检验量进行周跳的探测和修复。
下面结合具备惯导辅助的B1频信号北斗***为例进行说明。
短基线下惯导辅助的北斗单频整周模糊度求解方法,包括如下步骤:
一、初始信息获取
(1)为有效消除北斗双差观测量的大气延时误差,将北斗***输出的北斗原始观测量按照北斗星座类型的轨道高度分组差分:MEO卫星一组,卫星数量为n,选其中仰角最大的一颗MEO卫星作为基准;IGSO卫星和GEO卫星一组,卫星数量为m,选其中仰角最大的一颗IGSO卫星/GEO卫星作为基准,构造北斗双差载波观测方程如下:
其中,为第j颗MEO卫星到基准MEO卫星(在MEO卫星组中标记为第1颗MEO卫星)的双差载波观测相位,为第j颗IGSO卫星/GEO卫星到基准IGSO卫星/GEO卫星(在IGSO卫星和GEO卫星组中标记为第1颗IGSO卫星/GEO卫星)的双差载波观测相位,为MEO卫星到载体之间的单位矢量,为GEO卫星/IGSO卫星到载体之间的单位矢量,为第j颗MEO卫星到基准MEO卫星的双差整周模糊度,为第j颗IGSO卫星/GEO卫星到基准IGSO卫星/GEO卫星的双差整周模糊度,为第j颗MEO卫星到基准MEO卫星的双差观测噪声,为第j颗IGSO卫星/GEO卫星到基准IGSO卫星/GEO卫星的双差观测噪声,λ为北斗B1频率信号的载波波长,为基线矢量。
(2)利用惯导输出姿态矩阵、平台误差角和天线配置信息估算基线矢量地心地固坐标及相应的基线矢量方差矩阵:
式中,XI为惯导估算基线矢量(ECEF坐标,地心地固坐标),XG为北斗天线构成的基线在机体坐标系内的向量,为理想姿态矩阵(机体系到地理系的转换矩阵),为惯导实际输出的姿态矩阵,为地理坐标系到地球坐标系的转换矩阵,φ为平台误差角向量,cov[φφT]是平台误差角方差阵(由组合滤波中的方差矩阵获得),I为单位矩阵,cov[XI]为惯导估算基线矢量的方差矩阵。
二、整周模糊度的确认
(3)将由步骤(2)得到的惯导估算基线矢量XI代入步骤(1)得到的北斗双差载波观测方程中,利用递推最小二乘法求整周模糊度浮点解及其协方差矩阵。
(31)将惯导估算基线矢量XI代入北斗双差载波观测方程中,变为:
将上述两式均简化为R=-λa+v;其中,或a为对应分组的整周模糊度矢量,v为对应分组的测量噪声。
(32)将用惯导估算基线矢量XI来代替,R=-λa+v左侧R的协方差矩阵为:
(33)对整周模糊度矢量a进行递推最小二乘估计,得到整周模糊度浮点解为:
整周模糊度浮点解的协方差矩阵为:
其中,为整周模糊度浮点解的递推解,Qk为整周模糊度浮点解的协方差矩阵的递推解,PL为北斗原始观测量的方差矩阵。
(4)利用排序白化算法对整周模糊度浮点解进行整数Z变换,使得整周模糊度浮点解去相关。
(41)确定排序矩阵P,首先将矩阵P设定为与矩阵维数相同的零矩阵,然后将整周模糊度协方差矩阵中的正对角线元素Qii按升序排列,若Qii在排序后位于矩阵对角线元素的第m位,则令矩阵P中的元素Pmi=1,否则为零,计算得到矩阵P中所有元素的值后,计算排序后的矩阵其中i=1、2、…,m=1、2、…;
举例来说,初始时 首先设定 然后将中的正对角线元素进行升序排序,得到顺序为Q22、Q33、Q11,那么:位于第m=1位的元素i取值为2,则P12=1;位于m=2位的元素i取值为3,则P23=1;位于第m=3位的元素i取值为1,则P31=1;最终得到 计算得到
(42)对排序后的矩阵QP进行LDLT整数分解,得到下三角矩阵L后代入到Q=[L]-1QP[LT]-1中进行计算。
(43)基于(41)和步骤(42)两个步骤,不断地进行迭代,直到矩阵[L]-1为或近似为单位矩阵结束;最终得到整周模糊度去相关变换矩阵为:Z=[L1]-1P1[L2]-1P2…[Lk]-1Pk,则Z变换后的整周模糊度浮点解为其中k为迭代次数,Li为第i次迭代得到的L,Pi为第i次迭代得到的P。
(5)利用Z变换使整周模糊度浮点解去相关后,利用序贯最小二乘的方法搜索整周模糊度整数解。
(51)确定整周模糊度整数解的搜索范围:对Z变换后的整周模糊度浮点解进行取整,得到最接近的整数矢量z0。将z0中的某一个元素保持不变,该元素以外的其他对应元素取次接近于的整数值,得到一个矢量c,矢量c中的元素包括该不变元素以及其他元素取次接近整数的值并保持位置对应;重复该步骤直至z0中所有元素均有且只有一次保持不变的计算机会,得到对应z0中元素个数的矢量c;将所有得到的矢量c分别带入到公式中得到计算值,取其中次小的计算值作为搜索边界χ2的值,
举例来说,Z变换后的若整周模糊度那么z0=[1,3,6],分别取z0中的一个元素保持不变,得到三个c值,分别为[1,4,5]、[2,3,5]和[2,4,6],将三个c值分别带入到中进行计算,以得到χ2。
(52)将整周模糊度整数解的搜索条件表示为其中z为待搜索的整周模糊度整数解。
将n和m统一使用s进行表示,对进行LDLT分解,可以变为:
式中,lji是LDLT分解中矩阵L的元素,di是LDLT分解中矩阵D中正对角线上的元素,而且表示的方差,表示已知的情况下对zi的估计,将整周模糊度整数解的搜索条件表示为
(53)对χ2进行序贯条件方差调整,可以得到每一个整周模糊度整数解zj对应的搜索边界,具体方法如下:
第s-1个整周模糊度整数解的搜索边界为:第s-2个整周模糊度整数解的搜索边界为: ……,第一个整周模糊度整数解的搜索边界为: 搜索顺序是zs-1,zs-2,…,z1。
(6)将搜索得到的整周模糊度整数解z作为整周模糊度整数解的最优解通过Z变换的逆变换得到整周模糊度的最优解为
(7)将获得的整周模糊度的最优解通过ratio法进行检验,判断是否成立,若成立则整周模糊度的最优解成功固定,其中Ω1为最小残差平方和;Ω2为次小残差平方和;ratio为阈值,由经验确定。
三、周跳的探测和修复
(8)利用惯导估算基线矢量XI结合北斗三差载波观测值构造检验量,根据检验量大小判断周跳是否发生,若发生则利用检验量取整的方式估算周跳值并进行修复。
(81)将步骤(1)中的北斗双差载波相位载波观测方程变换为:
y=Bb-λa+ε
其中,或B为对应分组的设计矩阵或b为基线矢量,λ为北斗B1频率信号的载波波长,a为对应分组的整周模糊度矢量,v为对应分组的测量噪声。
(82)将两个相邻时刻的北斗双差载波观测方程相减以得到北斗三差观测方程:
△y(tk+1,tk)=y(tk+1)-y(tk)=B(tk+1)b(tk+1)-λa(tk+1)+ε(tk+1)-[B(tk)b(tk)-λa(tk)+ε(tk)]
其中,y(tk)表示tk时刻的北斗三差载波观测量,B(tk)表示tk时刻的设计矩阵,b(tk)表示tk时刻的基线矢量,a(tk)表示tk时刻的整周模糊度,ε(tk)表示tk时刻的测量噪声。
将相邻时刻载体到卫星的单位矢量看作不变,B用B(tk)来代替,将北斗三差载波观测方程变换为:
△y(tk+1,tk)=B(tk)·△b(tk+1,tk)-λ△a(tk+1,tk)+△ε(tk+1,tk)
其中,△b(tk+1,tk)=b(tk+1)-b(tk)为两个时刻基线矢量之差,b(tk)和b(tk+1)可利用惯导估算基线矢量来替代,其中
(83)设计检验量为T=(B(tk)·△b(tk+1,tk)-△y(tk+1,tk))/λ,若检验量的值在0附近,说明未发生周跳;若检验量T的绝对值大于1/2,则发生周跳,用检验量T四舍五入取整来计算发生的周跳数目。
至此,一种短基线下基于惯导辅助的北斗卫星单频整周模糊度解算方法流程结束。
以上所述仅是本发明的优选实施方式,应当指出:对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。