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

何瑞珠，劉 成，黃 康

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

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

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

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

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

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

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

(1)

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

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

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

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

(2)

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

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

(3)

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

(4)

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

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

(5)

(6)

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

ΔX＾，并定義如下：

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

(7)

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

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

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

(8)

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

＾-2GTΔρ.

(9)

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

(10)

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

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

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

(11)

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

(12)

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

(13)

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

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

2.1 基本思想

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

(14)

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

(15)

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

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

(16)

式中，Δx0、Δy0、Δz0與Δt0為各坐標(biāo)方向及接收機(jī)鐘差誤差的真實(shí)值。

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

2.2 實(shí)現(xiàn)方法

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

H=(GTWG)-1GTW.

(17)

當(dāng)有i顆觀測衛(wèi)星參與解算時(shí)，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)造一個(gè)權(quán)矩陣W使得(21)式第1個(gè)方程的右端滿足：

).

(23)

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

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

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

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

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

(24)

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

(25)

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

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

(26)

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

(27)

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

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

(28)

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

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

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

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

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

3.1 仿真實(shí)驗(yàn)

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

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

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

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

圖2 定位解算點(diǎn)分布情況示意圖。(a) 最小二乘算法定位點(diǎn)分布;(b) 分布加權(quán)算法定位點(diǎ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)計(jì)特性可知，在同等條件下進(jìn)行定位解算，無論是在各坐標(biāo)方向還是三維位置方向上，分步加權(quán)算法都能夠獲得更準(zhǔn)確的定位結(jié)果，提高三維定位精度0.7 m(均方誤差)，減小相對誤差約10%～20%。并且，分步加權(quán)算法在各坐標(biāo)軸方向和三維位置上的均方根誤差都比最小二乘算法更小，說明定位解的波動(dòng)較小，方法具有更高的穩(wěn)定性。

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

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

4 總 結(jié)

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

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

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

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

[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é)報(bào), 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é)報(bào), 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] 王世儒, 王金金, 馮有前, 等. 計(jì)算方法[M]. 第二版. 西安: 西安電子科技大學(xué)出版社, 2004.

[9] 席少林, 趙鳳治. 最優(yōu)化計(jì)算方法[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é)院知識(shí)創(chuàng)新工程重要方向項(xià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