何瑞珠，劉 成，黃 康

(1. 中國科學(xué)院國家天文臺，北京 100012；2. 中國科學(xué)院大學(xué)，北京 100049)

衛(wèi)星定位中的一種分步加權(quán)解算方法*

何瑞珠1,2，劉 成1,2，黃 康1

(1. 中國科學(xué)院國家天文臺，北京 100012；2. 中國科學(xué)院大學(xué)，北京 100049)

衛(wèi)星定位中，當(dāng)可視衛(wèi)星數(shù)目多于4顆時常采用加權(quán)最小二乘(Web Login Server, WLS)算法，對各衛(wèi)星解算權(quán)重進行重新評估而獲得最優(yōu)解。然而由于受到多重因素的影響，權(quán)矩陣W的構(gòu)造與確定一直是各類加權(quán)算法中的重點和難點。從線性測量方程組出發(fā)，通過研究迭代解算過程中用戶等效偽距測量誤差對坐標位置誤差的傳遞與放大規(guī)律，提出了一種新的加權(quán)最小二乘解算方法，及其權(quán)矩陣的具體構(gòu)造與實現(xiàn)方法，從而對各未知數(shù)進行分步加權(quán)與分離解算。通過全球定位系統(tǒng)(Global Positioning System, GPS)實測實驗，對方法的可行性和精度水平進行了分析與驗證。結(jié)果表明，使用該分步加權(quán)方法進行定位，解算結(jié)果準確度更高、穩(wěn)定性更好。

衛(wèi)星定位；線性方程組；權(quán)矩陣；加權(quán)最小二乘算法

衛(wèi)星定位系統(tǒng)的精度主要取決于兩方面：一是觀測衛(wèi)星在空間的分布情況，通常稱為衛(wèi)星分布幾何圖形；二是各觀測量的精度。為評價定位結(jié)果，常采用精度因子(Dilution of Precision, DOP)的概念，也稱為衛(wèi)星圖形強度因子或精度衰減因子。由此，最終定位解算結(jié)果可表示為

(1)

GDOP及σUERE的改變均能對定位誤差造成影響。當(dāng)可視衛(wèi)星數(shù)目增加時GDOP值隨之單調(diào)遞減，從而能夠獲得更優(yōu)的定位精度。當(dāng)衛(wèi)星數(shù)目多于4顆時，常采用最小二乘算法(Least Square method, LS)進行解算。若各測量值之間相互獨立且符合高斯正態(tài)分布，則最小二乘算法能夠使得各觀測偽距具有最小殘差平方和，從而獲得最優(yōu)估計解。然而，在顧及電離層延遲、對流層延遲及多路徑誤差等情況下，各實際測量值之間往往不是相互獨立，也不是同分布的，甚至并不符合正態(tài)分布。此時，若仍沿用最小二乘算法則難以保證其最優(yōu)性。這種情況下，可以引入某個加權(quán)矩陣W，改變不同衛(wèi)星之間的權(quán)重，更合理地估計最優(yōu)位置解。這種方法稱為加權(quán)最小二乘估計算法(Weighted Least-Square estimate, WLS)[1-4]。

加權(quán)最小二乘估計算法更符合實際定位情況，能夠?qū)Ω饔^測衛(wèi)星進行更合理的衡量，獲得更高的解算精度。權(quán)矩陣W最常用的方法是作為一個對角矩陣W=diag(ω1,ω2, …,ωn)，對角元素ωi(i=1, 2, …,n)為定位解算中第i顆衛(wèi)星的權(quán)重系數(shù)，通過將偽距誤差方差對仰角作近似而獲得，隨著仰角的減小而單調(diào)遞增[1]。加權(quán)最小二乘估計算法求解中采用這樣一個協(xié)方差使得低仰角衛(wèi)星的權(quán)重降低，因為那些衛(wèi)星由于典型的多徑特征和殘留的電離層、對流層誤差而噪聲較大[1,5]。目前有一些關(guān)于加權(quán)新方法的有益探討和研究，文[3]利用模糊數(shù)學(xué)方法確定加權(quán)矩陣，文[6]提出將灰色關(guān)聯(lián)分析法與廣義延拓數(shù)學(xué)方法相結(jié)合構(gòu)造權(quán)矩陣。然而，由于σUERE取決于偽距測量精度、衛(wèi)星仰角及衛(wèi)星星歷數(shù)據(jù)質(zhì)量等多種因素，其間還存在錯綜復(fù)雜的相互作用，因此權(quán)矩陣的構(gòu)造與確定一直是一個難題[2-3,5]。本文提出了一種新的加權(quán)最小二乘算法，通過分步構(gòu)造加權(quán)矩陣的方法對各未知數(shù)進行分離解算。并且，通過全球定位系統(tǒng)實測實驗對方法的可行性和精度水平進行了分析與驗證。

1 衛(wèi)星定位中的加權(quán)最小二乘算法

衛(wèi)星導(dǎo)航定位中，測量得到用戶至各觀測衛(wèi)星的偽距值后，即可構(gòu)造得到地心地固坐標系(ECEF)下的觀測方程：

(2)

式中，xi、yi和zi為第i顆衛(wèi)星在地心地固坐標系下的坐標；ρi為第i顆衛(wèi)星至用戶的觀測偽距值；x、y、z和Δt是待求量，分別為用戶接收機位置坐標及接收機時鐘與系統(tǒng)時鐘之間的偏差。

方程組(2)經(jīng)線性化處理后可得GΔX=Δρ,

(3)

式中，G為方向余弦矩陣：

(4)

其中，ei1、ei2、ei3各項表示由近似用戶位置指向第i顆衛(wèi)星的單位矢量方向余弦。

ΔX為近似解與真解的偏差向量：

(5)

(6)

可以看出，真正求解的是用戶位置與時鐘誤差初值估計的校正量ΔX。記ΔX的最大似然估計為

ΔX＾，并定義如下：

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

(7)

式中，p(Δρ/ΔX)是對于一個固定的ΔX而言其測量值Δρ的概率密度函數(shù)。

若測量誤差之間相互獨立，且均符合均值為0、方差為σ2的高斯分布，則(7)式變?yōu)?/p>

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

(8)

將‖Δρ-GΔX＾‖2對ΔX＾求微分，則可將(8)式變換為

＾-2GTΔρ.

(9)

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

(10)

使得測量值矢量Δρ和GΔX之間的殘差平方和最小。其中，GΔX是根據(jù)ΔX的估計值計算得到的測量值矢量。

由于在實際定位中，各測量偽距不僅不是同分布、等精度的，而且還是相互關(guān)聯(lián)的，因此需要考

慮更一般化的情況。此時，最大似然估計ΔX＾可以表示為

(11)

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

(12)

記R-1=W，則(12)式可表達為ΔX＾=(GTWG)-1GTWΔρ,

(13)

(13)式所作的估計即稱為加權(quán)最小估計(WLS)，W即為所使用的加權(quán)矩陣。

2 分步加權(quán)解算方法

2.1 基本思想

由上述討論可知，最小二乘算法的目標函數(shù)可表示為

(14)

式中，Δρi為各衛(wèi)星偽距測量誤差。加權(quán)最小二乘算法的目標函數(shù)則為

(15)

它通過權(quán)矩陣W改變了解算估計過程中各衛(wèi)星起的作用和比重，使其更符合實際情況。

然而，最小二乘與加權(quán)最小二乘算法的數(shù)學(xué)估計目標均是使得全局殘差和最小，即使得X、Y、Z

(16)

式中，Δx0、Δy0、Δz0與Δt0為各坐標方向及接收機鐘差誤差的真實值。

這樣，通過對各未知數(shù)依次分步進行相應(yīng)的加權(quán)計算，改變原有解算模式，達到對各未知數(shù)進行分離解算的目的。此時不再要求權(quán)矩陣W使得所有未知數(shù)的全局殘差最小，而是分別使得各個未知數(shù)各自的單一殘差最小化，并認為此時為該未知數(shù)的最優(yōu)估計。最后，將各坐標方向上依次所得的最優(yōu)估計值組合在一起，構(gòu)成一組新的矛盾解，作為用戶位置坐標的最優(yōu)估計。

2.2 實現(xiàn)方法

將(13)式中的(GTWG)-1GTW項作為一個整體，并記為

H=(GTWG)-1GTW.

(17)

當(dāng)有i顆觀測衛(wèi)星參與解算時，H為

(18)

則(13)式可轉(zhuǎn)化為

(19)

式中，Δρ1，Δρ2，…，Δρi為i顆衛(wèi)星的偽距測量誤差。將(19)式寫成線性化形式可得

(20)

設(shè)各衛(wèi)星的等效偽距測量中誤差分別為σρ1,σρ2,σρ3,…,σρi，于是由誤差協(xié)方差傳播率可得[7]

(21)

(22)

因此，可以構(gòu)造一個權(quán)矩陣W使得(21)式第1個方程的右端滿足：

).

(23)

也即使得X坐標方向上的定位解算誤差σx最小。此時，所構(gòu)造的權(quán)矩陣W僅使得X坐標方向上達到局部最優(yōu)，而其他未知數(shù)的估計值卻可能并非最優(yōu)；然而，在這一次的分步解算中，并不關(guān)心其他未知數(shù)的誤差估計情況。

對未知數(shù)x作出最優(yōu)估計后，同理再對其他各未知數(shù)依次進行解算即可。這種思想，就是先構(gòu)造一個權(quán)矩陣W，改變解算過程中誤差的分布和影響，使得X方向上的誤差估計最?。辉贅?gòu)造另一個權(quán)矩陣W，使得Y方向上的誤差估計最小，依次類推。最后，將各未知數(shù)的最優(yōu)估計值組合在一起，構(gòu)成一組新的矛盾解，作為用戶位置坐標的最優(yōu)解。為獲得更高精度的解算結(jié)果，該方法的解算過程并非一步完成，而是分為若干個步驟，故稱為“分步加權(quán)解算方法”。

2.3 加權(quán)矩陣W的構(gòu)造

2.3.1 多元函數(shù)極值求解構(gòu)造方法

以5顆衛(wèi)星時的解算情況為例，當(dāng)i=5時，方向余弦系數(shù)矩陣G為

(24)

加權(quán)對角矩陣W為

(25)

在實際解算過程中，G為一個常數(shù)矩陣。將(24)、(25)式代入(17)式，可得到一個由ωi(i=1～5)構(gòu)成的矩陣。并且，對應(yīng)于H矩陣中每一個元素hji(i=1～5,j=1～4)，均可用一個由ωi(i=1～5)構(gòu)成的代數(shù)式表達：

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

(26)

將(26)式得到的各項hji表達式及各衛(wèi)星測量等效偽距誤差值σρi(i=1, 2, …, 5)一起代入(21)式，可構(gòu)造分別對應(yīng)于未知數(shù)x、y、z和Δt的目標函數(shù)：

(27)

此時，確定權(quán)矩陣對角元素ωi(i=1～5)的過程實際上是一個多元函數(shù)的極值求解問題。首先，對(27)式的第1個目標函數(shù)fx求極小值min(fx)，得到與之對應(yīng)的一組權(quán)值ωi(i=1～5)，并生成權(quán)矩陣W。將此權(quán)矩陣W代入(17)式并重新進行坐標位置迭代解算，即可得到關(guān)于未知數(shù)x的最優(yōu)估計。

同理，可以依次得到其余各變量各自的最優(yōu)估計值。此過程中，W為一正定對角矩陣，各對角元素滿足下列歸一化要求：

(28)

這樣，既把權(quán)矩陣數(shù)據(jù)映射到0～1的范圍內(nèi)進行處理，也為多元函數(shù)的極值求解提供了區(qū)間范圍約束條件。當(dāng)衛(wèi)星數(shù)目較多時，無論多元函數(shù)是否為凸函數(shù)特性，均能夠在約束條件和區(qū)間內(nèi)求解得到其局域極值點。具體的極值求解算法，可采用極速下降法、單純形法等常用的非線性求解算法[8-10]。

2.3.2 可能性模板矩陣構(gòu)造方法

多元函數(shù)極值求解的方法能夠從數(shù)學(xué)原理上嚴格確定權(quán)矩陣元素的值，但在實際運用中，對一個多元函數(shù)求極值的計算不僅繁瑣，還需花費大量的計算時間，限制了方法的實用性。為此，可在求解過程中結(jié)合實際情況引入枚舉組合的方法，對權(quán)矩陣中各對角元素的可能取值進行列舉與排列組合。在解算中，根據(jù)函數(shù)值極小化的判定條件進行計算和判斷，篩選出最優(yōu)的一組權(quán)值組合。

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)

…

矩陣的實際個數(shù)由列舉組合的具體細化程序確定。

在分步加權(quán)解算時代入這些模板矩陣進行計算，并根據(jù)判定條件選出最合適的矩陣模板作為最終權(quán)矩陣W即可。這樣，雖不能嚴格計算獲得使目標函數(shù)極小化的權(quán)矩陣元素值，但卻能夠避免權(quán)矩陣的反復(fù)構(gòu)造，減少計算量，提高解算速度，在算法準確性與實用性上取得更好的平衡。由此，分布加權(quán)解算算法的使用步驟可表示為圖1。

3 仿真結(jié)果與分析

3.1 仿真實驗

為驗證方法的可靠性，在北京地區(qū)某地對全球定位系統(tǒng)衛(wèi)星進行了連續(xù)觀測，并選出了5顆具有最佳精度因子值的可視衛(wèi)星進行絕對單點定位。定位數(shù)據(jù)輸出時間間隔為1 s。利用輸出數(shù)據(jù)在Matlab仿真軟件中分別采用最小二乘算法和分步加權(quán)解算方法進行坐標解算，兩種方法的定位點分布情況如圖2。

從圖2可以看出，分布加權(quán)算法相對于最小二乘算法有著明顯的改善，減小了誤差分布的范圍。統(tǒng)計其誤差特性結(jié)果，如表1。

圖1 分步加權(quán)解算方法流程示意圖

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

圖2 定位解算點分布情況示意圖。(a) 最小二乘算法定位點分布;(b) 分布加權(quán)算法定位點分布

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中的誤差統(tǒng)計特性可知，在同等條件下進行定位解算，無論是在各坐標方向還是三維位置方向上，分步加權(quán)算法都能夠獲得更準確的定位結(jié)果，提高三維定位精度0.7 m(均方誤差)，減小相對誤差約10%～20%。并且，分步加權(quán)算法在各坐標軸方向和三維位置上的均方根誤差都比最小二乘算法更小，說明定位解的波動較小，方法具有更高的穩(wěn)定性。

表1 兩種解算方法定位結(jié)果誤差特性Table 1 Statistics of errors of position solutions yielded by the LS method and the new method

在解算時間上，最小二乘算法單點定位解算平均耗時約為4.92×10-4s；分步加權(quán)算法由于有著對可能性模板矩陣進行數(shù)值計算及再次求解位置坐標的過程，單點定位解算的平均耗時約為3.10×10-3s。

4 總 結(jié)

本文從衛(wèi)星定位線性測量方程組出發(fā)，通過研究迭代解算過程中偽距測量誤差對待求未知數(shù)誤差的放大規(guī)律，提出了一種加權(quán)最小二乘解算方法，及其權(quán)矩陣的具體構(gòu)造與實現(xiàn)方法。該方法對測量方程組中的各個未知數(shù)進行分步加權(quán)與單獨解算，在每一次的權(quán)矩陣構(gòu)造過程中，并不要求測量方程組的全局偽距殘差和最小，而是依次構(gòu)造使得各個未知數(shù)達到獨立最優(yōu)估計的權(quán)矩陣W，對未知數(shù)進行分離解算；最后，將分步解算得到的各未知數(shù)最優(yōu)估計值作為新的用戶坐標位置。通過仿真計算與分析，表明該分步加權(quán)解算方法能夠在各坐標軸方向和三維坐標方向上提高定位解的精度。

該方法也存在著一些需要繼續(xù)研究和改進的問題，主要體現(xiàn)在：

(1)方法需要利用各觀測衛(wèi)星等效偽距測量誤差的先驗統(tǒng)計信息。這些信息會在衛(wèi)星星歷中給出，但在實際定位時卻通常不會嚴格一致，因此會給權(quán)矩陣的計算帶來誤差，從而影響算法的精度。

(2)由于方法復(fù)雜度更高，因此增加了計算耗時。雖然仍能滿足一般性的、非高動態(tài)實時單點定位需求，但在觀測衛(wèi)星數(shù)目較多時，仍然需要考慮和研究更為簡單與快速的權(quán)矩陣確定方法，以減少計算量，提高解算速度，使算法具有更好的實用性和可靠性。

[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, 談?wù)怪? 衛(wèi)星導(dǎo)航幾何因子的分析和仿真[J].電子學(xué)報, 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] 寧春林, 施滸立, 李圣明, 等. 一種構(gòu)造WDOP中加權(quán)矩陣的新方法[J]. 宇航學(xué)報, 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] 武漢大學(xué)測繪學(xué)院測量平差學(xué)科組. 誤差理論與測量平差基礎(chǔ)[M]. 武漢: 武漢大學(xué)出版社, 2003.

[8] 王世儒, 王金金, 馮有前, 等. 計算方法[M]. 第二版. 西安: 西安電子科技大學(xué)出版社, 2004.

[9] 席少林, 趙鳳治. 最優(yōu)化計算方法[M]. 上海: 上?？茖W(xué)技術(shù)出版社, 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

國家自然科學(xué)基金 (61001109)；中國科學(xué)院知識創(chuàng)新工程重要方向項目 (KGCX2-EW-4071) 資助.

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

何瑞珠，女，碩士. 研究方向：衛(wèi)星導(dǎo)航與通信. Email: heruizhu11@mails.ucas.ac.cn

TN967.1

A

1672-7673(2015)01-0036-08