卫星定位中的一种分步加权解算方法*

何瑞珠，刘 成，黄 康

(1. 中国科学院国家天文台，北京 100012；2. 中国科学院大学，北京 100049)

卫星定位中的一种分步加权解算方法*

何瑞珠1,2，刘 成1,2，黄 康1

(1. 中国科学院国家天文台，北京 100012；2. 中国科学院大学，北京 100049)

卫星定位中，当可视卫星数目多于4颗时常采用加权最小二乘(Web Login Server, WLS)算法，对各卫星解算权重进行重新评估而获得最优解。然而由于受到多重因素的影响，权矩阵W的构造与确定一直是各类加权算法中的重点和难点。从线性测量方程组出发，通过研究迭代解算过程中用户等效伪距测量误差对坐标位置误差的传递与放大规律，提出了一种新的加权最小二乘解算方法，及其权矩阵的具体构造与实现方法，从而对各未知数进行分步加权与分离解算。通过全球定位系统(Global Positioning System, GPS)实测实验，对方法的可行性和精度水平进行了分析与验证。结果表明，使用该分步加权方法进行定位，解算结果准确度更高、稳定性更好。

卫星定位；线性方程组；权矩阵；加权最小二乘算法

卫星定位系统的精度主要取决于两方面：一是观测卫星在空间的分布情况，通常称为卫星分布几何图形；二是各观测量的精度。为评价定位结果，常采用精度因子(Dilution of Precision, DOP)的概念，也称为卫星图形强度因子或精度衰减因子。由此，最终定位解算结果可表示为

(1)

GDOP及σUERE的改变均能对定位误差造成影响。当可视卫星数目增加时GDOP值随之单调递减，从而能够获得更优的定位精度。当卫星数目多于4颗时，常采用最小二乘算法(Least Square method, LS)进行解算。若各测量值之间相互独立且符合高斯正态分布，则最小二乘算法能够使得各观测伪距具有最小残差平方和，从而获得最优估计解。然而，在顾及电离层延迟、对流层延迟及多路径误差等情况下，各实际测量值之间往往不是相互独立，也不是同分布的，甚至并不符合正态分布。此时，若仍沿用最小二乘算法则难以保证其最优性。这种情况下，可以引入某个加权矩阵W，改变不同卫星之间的权重，更合理地估计最优位置解。这种方法称为加权最小二乘估计算法(Weighted Least-Square estimate, WLS)[1-4]。

加权最小二乘估计算法更符合实际定位情况，能够对各观测卫星进行更合理的衡量，获得更高的解算精度。权矩阵W最常用的方法是作为一个对角矩阵W=diag(ω1,ω2, …,ωn)，对角元素ωi(i=1, 2, …,n)为定位解算中第i颗卫星的权重系数，通过将伪距误差方差对仰角作近似而获得，随着仰角的减小而单调递增[1]。加权最小二乘估计算法求解中采用这样一个协方差使得低仰角卫星的权重降低，因为那些卫星由于典型的多径特征和残留的电离层、对流层误差而噪声较大[1,5]。目前有一些关于加权新方法的有益探讨和研究，文[3]利用模糊数学方法确定加权矩阵，文[6]提出将灰色关联分析法与广义延拓数学方法相结合构造权矩阵。然而，由于σUERE取决于伪距测量精度、卫星仰角及卫星星历数据质量等多种因素，其间还存在错综复杂的相互作用，因此权矩阵的构造与确定一直是一个难题[2-3,5]。本文提出了一种新的加权最小二乘算法，通过分步构造加权矩阵的方法对各未知数进行分离解算。并且，通过全球定位系统实测实验对方法的可行性和精度水平进行了分析与验证。

1 卫星定位中的加权最小二乘算法

卫星导航定位中，测量得到用户至各观测卫星的伪距值后，即可构造得到地心地固坐标系(ECEF)下的观测方程：

(2)

式中，xi、yi和zi为第i颗卫星在地心地固坐标系下的坐标；ρi为第i颗卫星至用户的观测伪距值；x、y、z和Δt是待求量，分别为用户接收机位置坐标及接收机时钟与系统时钟之间的偏差。

方程组(2)经线性化处理后可得GΔX=Δρ,

(3)

式中，G为方向余弦矩阵：

(4)

其中，ei1、ei2、ei3各项表示由近似用户位置指向第i颗卫星的单位矢量方向余弦。

ΔX为近似解与真解的偏差向量：

(5)

(6)

可以看出，真正求解的是用户位置与时钟误差初值估计的校正量ΔX。记ΔX的最大似然估计为

ΔX＾，并定义如下：

ΔX＾=arg maxp(Δρ/ΔX),

(7)

式中，p(Δρ/ΔX)是对于一个固定的ΔX而言其测量值Δρ的概率密度函数。

若测量误差之间相互独立，且均符合均值为0、方差为σ2的高斯分布，则(7)式变为

=arg min‖Δρ-GΔX＾‖2,

(8)

将‖Δρ-GΔX＾‖2对ΔX＾求微分，则可将(8)式变换为

＾-2GTΔρ.

(9)

令(9)式为0并进行求解，可得 ΔX＾=(GTG)-1GTΔρ.

(10)

使得测量值矢量Δρ和GΔX之间的残差平方和最小。其中，GΔX是根据ΔX的估计值计算得到的测量值矢量。

由于在实际定位中，各测量伪距不仅不是同分布、等精度的，而且还是相互关联的，因此需要考

虑更一般化的情况。此时，最大似然估计ΔX＾可以表示为

(11)

ΔX＾=(GTR-1G)-1GTR-1Δρ.

(12)

记R-1=W，则(12)式可表达为ΔX＾=(GTWG)-1GTWΔρ,

(13)

(13)式所作的估计即称为加权最小估计(WLS)，W即为所使用的加权矩阵。

2 分步加权解算方法

2.1 基本思想

由上述讨论可知，最小二乘算法的目标函数可表示为

(14)

式中，Δρi为各卫星伪距测量误差。加权最小二乘算法的目标函数则为

(15)

它通过权矩阵W改变了解算估计过程中各卫星起的作用和比重，使其更符合实际情况。

然而，最小二乘与加权最小二乘算法的数学估计目标均是使得全局残差和最小，即使得X、Y、Z

(16)

式中，Δx0、Δy0、Δz0与Δt0为各坐标方向及接收机钟差误差的真实值。

这样，通过对各未知数依次分步进行相应的加权计算，改变原有解算模式，达到对各未知数进行分离解算的目的。此时不再要求权矩阵W使得所有未知数的全局残差最小，而是分别使得各个未知数各自的单一残差最小化，并认为此时为该未知数的最优估计。最后，将各坐标方向上依次所得的最优估计值组合在一起，构成一组新的矛盾解，作为用户位置坐标的最优估计。

2.2 实现方法

将(13)式中的(GTWG)-1GTW项作为一个整体，并记为

H=(GTWG)-1GTW.

(17)

当有i颗观测卫星参与解算时，H为

(18)

则(13)式可转化为

(19)

式中，Δρ1，Δρ2，…，Δρi为i颗卫星的伪距测量误差。将(19)式写成线性化形式可得

(20)

设各卫星的等效伪距测量中误差分别为σρ1,σρ2,σρ3,…,σρi，于是由误差协方差传播率可得[7]

(21)

(22)

因此，可以构造一个权矩阵W使得(21)式第1个方程的右端满足：

).

(23)

也即使得X坐标方向上的定位解算误差σx最小。此时，所构造的权矩阵W仅使得X坐标方向上达到局部最优，而其他未知数的估计值却可能并非最优；然而，在这一次的分步解算中，并不关心其他未知数的误差估计情况。

对未知数x作出最优估计后，同理再对其他各未知数依次进行解算即可。这种思想，就是先构造一个权矩阵W，改变解算过程中误差的分布和影响，使得X方向上的误差估计最小；再构造另一个权矩阵W，使得Y方向上的误差估计最小，依次类推。最后，将各未知数的最优估计值组合在一起，构成一组新的矛盾解，作为用户位置坐标的最优解。为获得更高精度的解算结果，该方法的解算过程并非一步完成，而是分为若干个步骤，故称为“分步加权解算方法”。

2.3 加权矩阵W的构造

2.3.1 多元函数极值求解构造方法

以5颗卫星时的解算情况为例，当i=5时，方向余弦系数矩阵G为

(24)

加权对角矩阵W为

(25)

在实际解算过程中，G为一个常数矩阵。将(24)、(25)式代入(17)式，可得到一个由ωi(i=1～5)构成的矩阵。并且，对应于H矩阵中每一个元素hji(i=1～5,j=1～4)，均可用一个由ωi(i=1～5)构成的代数式表达：

hji=fji(ω1,ω2,ω3,ω4,ω5) (i=1～5,j=1～4).

(26)

将(26)式得到的各项hji表达式及各卫星测量等效伪距误差值σρi(i=1, 2, …, 5)一起代入(21)式，可构造分别对应于未知数x、y、z和Δt的目标函数：

(27)

此时，确定权矩阵对角元素ωi(i=1～5)的过程实际上是一个多元函数的极值求解问题。首先，对(27)式的第1个目标函数fx求极小值min(fx)，得到与之对应的一组权值ωi(i=1～5)，并生成权矩阵W。将此权矩阵W代入(17)式并重新进行坐标位置迭代解算，即可得到关于未知数x的最优估计。

同理，可以依次得到其余各变量各自的最优估计值。此过程中，W为一正定对角矩阵，各对角元素满足下列归一化要求：

(28)

这样，既把权矩阵数据映射到0～1的范围内进行处理，也为多元函数的极值求解提供了区间范围约束条件。当卫星数目较多时，无论多元函数是否为凸函数特性，均能够在约束条件和区间内求解得到其局域极值点。具体的极值求解算法，可采用极速下降法、单纯形法等常用的非线性求解算法[8-10]。

2.3.2 可能性模板矩阵构造方法

多元函数极值求解的方法能够从数学原理上严格确定权矩阵元素的值，但在实际运用中，对一个多元函数求极值的计算不仅繁琐，还需花费大量的计算时间，限制了方法的实用性。为此，可在求解过程中结合实际情况引入枚举组合的方法，对权矩阵中各对角元素的可能取值进行列举与排列组合。在解算中，根据函数值极小化的判定条件进行计算和判断，筛选出最优的一组权值组合。

diag((4/40)1/2,(4/40)1/2,(4/40)1/2,(4/40)1/2,(24/40)1/2)

diag((9/40)1/2,(14/40)1/2,(4/40)1/2,(4/40)1/2,(9/40)1/2)

diag((14/40)1/2,(14/40)1/2,(4/40)1/2,(4/40)1/2,(4/40)1/2)

diag((19/40)1/2,(4/40)1/2,(9/40)1/2,(4/40)1/2,(4/40)1/2)

…

矩阵的实际个数由列举组合的具体细化程序确定。

在分步加权解算时代入这些模板矩阵进行计算，并根据判定条件选出最合适的矩阵模板作为最终权矩阵W即可。这样，虽不能严格计算获得使目标函数极小化的权矩阵元素值，但却能够避免权矩阵的反复构造，减少计算量，提高解算速度，在算法准确性与实用性上取得更好的平衡。由此，分布加权解算算法的使用步骤可表示为图1。

3 仿真结果与分析

3.1 仿真实验

为验证方法的可靠性，在北京地区某地对全球定位系统卫星进行了连续观测，并选出了5颗具有最佳精度因子值的可视卫星进行绝对单点定位。定位数据输出时间间隔为1 s。利用输出数据在Matlab仿真软件中分别采用最小二乘算法和分步加权解算方法进行坐标解算，两种方法的定位点分布情况如图2。

从图2可以看出，分布加权算法相对于最小二乘算法有着明显的改善，减小了误差分布的范围。统计其误差特性结果，如表1。

图1 分步加权解算方法流程示意图

Fig.1 A schematic diagram of the method for multiple-step Weighted Least-Square estimate

图2 定位解算点分布情况示意图。(a) 最小二乘算法定位点分布;(b) 分布加权算法定位点分布

Fig.2 Distributions of position solutions for satellite positioning. (a) A distribution from the LS method; (b) A distribution from the new method

3.2 有效性分析

由表1中的误差统计特性可知，在同等条件下进行定位解算，无论是在各坐标方向还是三维位置方向上，分步加权算法都能够获得更准确的定位结果，提高三维定位精度0.7 m(均方误差)，减小相对误差约10%～20%。并且，分步加权算法在各坐标轴方向和三维位置上的均方根误差都比最小二乘算法更小，说明定位解的波动较小，方法具有更高的稳定性。

表1 两种解算方法定位结果误差特性Table 1 Statistics of errors of position solutions yielded by the LS method and the new method

在解算时间上，最小二乘算法单点定位解算平均耗时约为4.92×10-4s；分步加权算法由于有着对可能性模板矩阵进行数值计算及再次求解位置坐标的过程，单点定位解算的平均耗时约为3.10×10-3s。

4 总 结

本文从卫星定位线性测量方程组出发，通过研究迭代解算过程中伪距测量误差对待求未知数误差的放大规律，提出了一种加权最小二乘解算方法，及其权矩阵的具体构造与实现方法。该方法对测量方程组中的各个未知数进行分步加权与单独解算，在每一次的权矩阵构造过程中，并不要求测量方程组的全局伪距残差和最小，而是依次构造使得各个未知数达到独立最优估计的权矩阵W，对未知数进行分离解算；最后，将分步解算得到的各未知数最优估计值作为新的用户坐标位置。通过仿真计算与分析，表明该分步加权解算方法能够在各坐标轴方向和三维坐标方向上提高定位解的精度。

该方法也存在着一些需要继续研究和改进的问题，主要体现在：

(1)方法需要利用各观测卫星等效伪距测量误差的先验统计信息。这些信息会在卫星星历中给出，但在实际定位时却通常不会严格一致，因此会给权矩阵的计算带来误差，从而影响算法的精度。

(2)由于方法复杂度更高，因此增加了计算耗时。虽然仍能满足一般性的、非高动态实时单点定位需求，但在观测卫星数目较多时，仍然需要考虑和研究更为简单与快速的权矩阵确定方法，以减少计算量，提高解算速度，使算法具有更好的实用性和可靠性。

[1] Kaplan E D, Hegarty C J. Understanding GPS: principles and applications[M]. 2nded. UK: Artech House London, 1996.

[2] Sairo H, Akopian D, Takala J. Weighted dilution of precision as quality measure in satellite positioning[J]. IEEE Radar, Sonar and Navigation, 2003, 150(6): 430-436.

[3] Yang Y, Miao L J. GDOP results in all-in-view positioning and in four optimum satellites positioning with GPS PRN codes ranging[C]// PLANS 2004 Position Location and Navigation Symposium. 2004: 723-727.

[4] 丛丽, Ahmed I Abidat, 谈展中. 卫星导航几何因子的分析和仿真[J].电子学报, 2006, 34(12): 2204-2208. Cong Li, Ahmed I Abidat, Tan Zhanzhong. Analysis and simulation of the GDOP of satellite navigaion[J]. Acta Electronica Sinica, 2006, 34(12): 2204-2208.

[5] Misra P, Enge P. Global Positioning System-Signal, Measurements, and Performance[M]. 2nded. Lincoln: Ganga-Jamuna Press, 2006.

[6] 宁春林, 施浒立, 李圣明, 等. 一种构造WDOP中加权矩阵的新方法[J]. 宇航学报, 2009, 30(2): 526-531. Ning Chunlin, Shi Huli, Li Shengming, et al. A new method of constructing weighting matrix of WDOP[J]. Journal of Astronautics, 2009, 30(2): 526-531.

[7] 武汉大学测绘学院测量平差学科组. 误差理论与测量平差基础[M]. 武汉: 武汉大学出版社, 2003.

[8] 王世儒, 王金金, 冯有前, 等. 计算方法[M]. 第二版. 西安: 西安电子科技大学出版社, 2004.

[9] 席少林, 赵凤治. 最优化计算方法[M]. 上海: 上海科学技术出版社, 1983.

[10]Zaguskin O O. Solution of Algebraic and Transcendental Equations[M]. New York: Pergamon Press, 1961.

CN 53-1189/P ISSN 1672-7673

A Method of Multiple-Step Weighted Least-Square Estimate for Satellite Positioning

He Ruizhu1,2, Liu Cheng1,2, Huang Kang1

(1. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China; 2. University of Chinese Academy of Sciences, Beijing 100049, China, Email: heruizhu11@mails.ucas.ac.cn)

The method of Weighted Least-Square estimate (WLS) is often used for satellite positioning to reevaluate the weights for various linkable satellites (of more than 4) in the fitting objective function for deriving optimal position solutions. However, because of a number of complex factors construction and determination of the important weight matrixWare difficult for all algorithms involving weighting for satellites. Starting from linear equations for optimal position solutions, we have studied the patterns of propagations of errors of equivalent pseudo ranges into coordinate errors of position solutions. The propagation patterns involve accumulation of errors in the iteration processes leading to solutions. We subsequently propose a new WLS method together with the approaches of construction and practical calculation of the weight matrix. What is mainly modified by the new method from the position-solving process of the original WLS method is as follows. In the new method arguments whose values are not known beforehand are solved in separate steps, with different steps having independent weight matrices. The new method effectively avoids the weight-matrix calculation that is used to minimize the global residual (residual-squared sum) of all the arguments. Instead, it minimizes the residual of each of the arguments separately. The best estimates of all coordinates constitute the optimal position solution for the location of a user. We have carried out a GPS measurement experiment to test feasibility and accurateness of the new method. The test results show that the new WLS method can appreciably improve the accuracy and stability of a position solution in satellite positioning.

Satellite positioning; Linear equations; Weight matrix; Weighted Least-Square estimate

国家自然科学基金 (61001109)；中国科学院知识创新工程重要方向项目 (KGCX2-EW-4071) 资助.

2014-03-13；修定日期：2014-04-14

何瑞珠，女，硕士. 研究方向：卫星导航与通信. Email: heruizhu11@mails.ucas.ac.cn

TN967.1

A

1672-7673(2015)01-0036-08