GPS 定位算法研究3常 青 高 晖 柳重堪 张其善 中国工程物理研究院科学技术基金资助项目 收稿日期:1998203209 文 摘 给出卫星信号接收时刻GPS 时的计算方法,然后在Gau ss -N ew ton 法的基础上得到一个正常情况下的完整的GPS 迭代定位算法。
最后讨论只有三颗可见卫星时的处理方法。
主题词 全球定位系统 导航 矩阵 算法1 引 言GPS 系统(Global Po siti on ing System )是美国国防部研制的第二代卫星导航与定位系统。
它能为全球的用户提供全天候、连续、实时的高精度位置、速度和时间信息,因此在军事、民事方面都具有重要的意义和应用价值。
目前,GPS 已在飞机、水面船只和陆地车辆的导航定位系统中得到了广泛的应用。
就GPS 定位的观测量而言,有伪距和载波相位之分。
就GPS 定位的定位方式而言,有单点定位和差分定位之分。
无论采用何种观测量和何种定位方式,其定位模型在求解上都归结为非线性最小二乘问题,需要迭代求解。
本文首先计算了由伪距单点定位模型得到的非线性最小二乘问题的目标函数的偏导数和H essian 矩阵,发现不仅目标函数的偏导数的计算非常简单,而且目标函数的H essian 矩阵中的二阶导数部分很小,可以忽略。
这表明把Gau ss -N ew ton 方法用于GPS 定位模型是非常有效的[1]。
然后给出卫星信号接收时刻GPS 时的计算方法,并在此基础上得到了一个完整的GPS 迭代定位算法。
该算法在每一步迭代中通过对信号传播时间、信号发射时刻、卫星钟差及卫星坐标的更新,来完成导航解的更新。
最后讨论只有三颗可见卫星时的处理方法。
2 伪距单点定位的数学模型某一时刻,接收机对卫星j 的伪距观测方程为Θζj =Θj+c (d t -d t j )+d j i on +d j trop +v(1)其中Θζj 为接收机对卫星j 的观测伪距,Θj=(x j -x )2+(y j -y )2+(z j -z )2为接收机与卫星j 之间的几何距离,(x ,y ,z )和(x j ,y j ,z j)分别为接收机和卫星j 在ECEF 坐标系中的坐标,d t 和d t j 分别为接收机和卫星j 的钟偏差,d j i on 和d j trop 分别为电离层和对流层时延引起的距离偏差,c 为光速,v 为观测噪声,j =1,2…n ,而n 为同步观测的卫星数,n ≥4。
在式(1)中,d t j ,d j i on 和d jtrop 可采用通用公式计算。
将上述三项移至式(1)的左边,包含未知参数x ,y ,z ,d t 的项留在式(1)的右边,并将左边记为Θλj ,cd t 记为b ,则式(1)可改写为Θλj =Θj+b +v(2)式(2)便是伪距单点定位的数学模型。
3 定位算法记X =[x ,y ,z ,b ]Tf i (X )=Θj+b -Θλj f (X )=[f 1(X ),f 2(X ),…f n (X )]T于是求解模型(2)归结为求解如下的非线性最小二乘问题。
m in{[f (X )]T f (X )}(3)可采用Gau ss -N ew ton 法求解式(3),函数f j (X )的梯度为x f j (X )=[ f j x , f j y , f j z , f j b]T=[-x j-xΘj,-y j-yΘj,-z j-zΘj,1]T=[-l j ,-m j ,-n j ,1]T(4)其中,l j ,m j ,n j 为点(x ,y ,z )到点(x j ,y j ,z j )的方向余弦。
于是式(3)中目标函数的导数g (X )为g (X )=2A (X )f (X )(5)其中A (X )=[ x f 1(X ), x f 2(X ),…, x f n (X )]=-l 1-l 2…-l n-m 1-m 2…-m n-n 1-n 2…-n n11…1(6)由于A (X )仅由接收机到卫星的方向余弦和数1组成,因此,当把Gau ss -N ew ton 方法用于GPS 定位模型时,目标函数偏导数的计算是非常简单的。
目标函数的H essian 矩阵G (X )为G (X )=2A (X )[A (X )]T+22nj =1f j (X )2X Xf j (X )(7)矩阵2X Xf j (X )的每一个元素都含有1Θj 。
以位于第一行第一列位置的元素为例,2f jx2=(Θj )2-(x j -x )2(Θj )3=1Θj +(x j -x )2(Θj )21Θj (8)其余元素也都具有类似的形式。
由于Θj 很大,达到两万公里,故矩阵 2X X f j (X )的每一个元素都很小。
由式(7)和式(8),并考虑到f j (X )与Θj之比是非常小的,便可以看出式(7)右边第二项的每一个元数都很小,将其忽略掉,得到G (X )≈T (X )=2A (X )[A (X )]T(9)式(9)说明,把Gau ss-N ew ton方法应用于GPS定位模型效果是好的。
设X k为Gau ss-N ew2 ton法第k步迭代的结果,则第k+1步的迭代公式为X k+1=X k+s k(10)其中,s k为在点X k处的Gau ss-N ew ton方向,它满足T(X k)s=-g(X k)(11)假设A(X k)为行满秩矩阵(这在实际中通常是可以满足的),则T(X k)为正定对称矩阵,使用Cho lesky分解算法,便可获得s k。
式(10)、(11)便是求解式(3)的迭代公式。
它从一个给定的初始值X0开始,当满足‖X k+1-X k‖<∆(12)时结束。
其中‖・‖为欧几里德距离,∆为事先给定的解的精度。
当算法起动以后,迭代的初始值X0可选为上一时刻迭代的结果。
由于卫星钟差和卫星坐标的计算都要涉及到卫星信号接收时刻GPS时的计算,下面给出卫星信号接收时刻GPS时的计算方法。
以N表示由接收机产生的时间标准T s的计数,则信号接收时刻GPS时的计算可分为以下三个步骤:1)将软件开始运行后从接收机获得的Grego rian时间(即目前通用的阳历)换算为GPS 星期和秒,并利用此时的T s计数N对N=0时的GPS星期数和秒数进行估计。
以y,m,d,h h, m m,s s代表软件开始运行后从接收机获得的Grego rian时间的年,月,日,时,分,秒,则从1980年1月6日零时至该时刻的年数为(y-1980)。
若(y-1980)能整除4,且m≤2,则闰日数l d 为l d=(y-1980)4否则闰日数为l d=[(y-1980)4]+1其中[・]表示取整数部分。
于是自1980年1月6日零时以来的总天数D为D=(y-1980)×365+d(m-1)+d+l d-6其中d(m-1)表示忽略闰年后前(m-1)个月的天数。
由上式可以计算出GPS星期数G w和秒数G s:G w=[D7], G s=(D◊7)×86400+h h×3600+m m×60+s s其中D◊7表示D除7后的余数。
上述算法是忽略了U TC(世界协调时)跳秒后的一种简化算法。
以z w、z s表示N=0时的GPS星期数和秒数,于是z w、z s的估计算法为:z w=G w,z s=G s-3600×T z-N×T s其中T z表示时区。
它指明了当地时间比U TC时间提前的小时数。
若z s<0,令z w=z w-1, z s=z s+604800若z s≥604800,令z w=z w+1, z s=z s-604800z w、z s的计算在软件开始运行时进行,利用它们可以求得任一信号接收时刻的GPS星期数和秒数。
2)利用信号接收时刻的T s 计数N ′计算信号接收时刻的GPS 星期数g w 和秒数g s 。
算法为:令g w =z w , g s =z s +N ′×T s若g s <0,令g w =g w -1, g s =g s +604800若g s ≥604800,令g w =g w +1, g s =g s -604800 3)对上一步的结果进行修正。
2)中求得的g w 和g s 实际上是接收机时间,它与GPS 时的关系为接收机时间-接收机钟偏差=GPS 时设估计的当前信号接收时刻的接收机钟偏差的星期数和秒数为R w 和R s ,接收机钟差变率为R r ,上一信号接收时刻对应的T s 计数为N ″,于是g w 、g s 的修正算法为:令g w =g w -R w , g s =g s -R s -d T R r其中d T =(N ′-N ″)T s (1-R r )若g s <0,令g w =g w -1, g s =g s +604800若g s ≥604800,令g w =g w +1, g s =g s -604800综上所述,求GPS 定位解的迭代步骤为:1)获得全部可见卫星的星历信息及Θζj 、d j i on 、d j trop 、g js 。
2)令k =0,卫星信号传播时间∃t j 的初始值∃t j0=0。
选取初始值X 0,允许误差∆>0。
3)按式g jk =g js -∃t jk 计算卫星信号发射时刻的GPS 时的秒数g j k 。
计算时,若g j k <0,令g jk =g j k +604800然后以g j k 计算卫星信号发射时刻的卫星坐标q j k (卫星坐标q j k 的计算见文献[4]～[6])。
考虑到ECEF 坐标系的旋转,按式p jk =co s Ηjk-sin Ηjk 0sin Ηjkco s Ηj k0001q jk对q j k 进行旋转。
其中p j k 为旋转后的卫星坐标,Ηj k =-∃t j k 8αe ,8αe 为W GS -84地球自转角速度。
4)以g jk 代入式d t j =a j f o +a jf 1(t -t j oc )+a jf 2(t -t j oc )2+∆tj计算卫星钟差d t j k 。
其中t joc 为时钟数据参考时间,a j f o ,a j f 1,a j f 2是导航电文子帧1中给出的多项式系数,∆t j 为相对论校正项,可以根据导航电文第2、3子帧中提供的卫星轨道参数计算。
假设用户为单频L 1用户,则上式应进一步修正为d t j =a j f o +a j f 1(t -t j oc )+a j f 2(t -t j oc )2+∆t j -t j G D其中t j G D 为时延差改正,在导航电文子帧1中提供给用户。
计算时,若g j k -t joc >302400,令g j k -t j oc =(g j k -t joc )-604800若g j k -t joc <-302400,令g jk -t j oc =(g jk -t joc )+604800 5)以X k 、p j k 计算Θj k 和A (X k )(Θj k =‖[X k ]3-P jk ‖,A (X k )按式(6)计算。