衛星定位和精度因子的改進方法

陳燦輝 張曉林

(北京航空航天大學 電子信息工程學院,北京 100191)

衛星定位和精度因子的改進方法

陳燦輝 張曉林

(北京航空航天大學 電子信息工程學院,北京 100191)

在衛星導航定位系統中,在精度因子計算和采用最小二乘法進行定位求解時,傳統上采用測量矩陣直接求逆方法來進行.為了克服矩陣求逆帶來的計算量大和數值穩定性差的不足,利用測量矩陣的對稱正定性,提出了一種基于矩陣 UTDU分解的定位解算和精度因子計算方法.改進方法具有嚴格的數學理論基礎,保證了方法的正確性和有效性.數值分析結果表明,相對直接求逆的傳統方法而言,在定位解算時,該方法能降低約 60%的運算量,而在精度因子計算中,約能降低 36%的運算量.且改進方法能大大降低求解矩陣的條件數,提高了求解的數值穩定性.

衛星導航;最小二乘;解算;精度因子;矩陣分解

全球導航衛星系統(GNSS,Global Navigation Satellite System)是一種以空間衛星為基礎的無線電導航與定位系統,該系統能為全世界任何地方的用戶全天候、全時間、連續和實時地提供三維位置、速度和時間 (PVT,Position,Velocity and Time)信息.由衛星導航定位系統確定的位置和時間的精度取決于各種因素錯綜復雜的相互作用.粗略來講,基于偽距的定位精度可以表示為精度因子(DOP,Dilution of Precision)和偽距誤差的乘積[1-3].為了提高定位精度,必須選擇精度因子小的衛星星座進行定位,計算精度因子是定位解算中必不可少的過程.另外,在采用 GNSS進行導航定位時,由于可見衛星數很多,在接收機容量和運算速度等因素的限制下,一般情形是不可能采用所有可見衛星來進行定位的,這時,就要進行選星,以選擇星座可用性滿足設計要求或是星座幾何精度因子(GDOP,Geometric Dilution of Precision)最優的少數可見星來進行定位,這也就需要進行精度因子求解.傳統上,計算精度因子采用矩陣求逆方法,計算量大,特別是在多星座組合導航中動態條件下需要頻繁進行選星操作時,其計算量更是一個極為突出的問題.此外,在各種偽距定位算法中,最小二乘法是一種比較簡單、基本而又有著廣泛應用的重要方法[3],該方法依據線性化定位模型進行迭代求解.由于最小二乘法每步迭代都涉及矩陣的乘法和求逆運算,因而整個迭代過程要多次重復這一過程,這無疑會增大導航解算的計算量與存儲量,從而影響實時定位的精度.因而,為了保證實時處理的要求,對接收機處理器速度的要求就大大提高了,這也就大大加重了用戶接收機的負擔,使其成本上升.

為了解決上述問題,探索新的、適合于衛星導航定位系統的快速獲取定位信息和精度因子的方法就顯得迫切而重要,并具有極為重要的現實意義和應用價值.

為了實現導航方程的快速求解,人們提出了一系列方法[4-6].文獻[4]提出了一種定位求解的遞推方法,但它是針對單星座系統來進行算法設計的,對多星座組合系統而言,隨著可見衛星數的大幅增加其計算量也會增加;文獻[5-6]提出了導航方程的非迭代解法,但它們均是針對單衛星導航系統來進行算法設計的,在多星座組合系統中難以直接采用.在 GNSS中,多星座組合導航勢成必然,針對 GNSS中快速導航定位的需要,并考慮到在導航定位系統中最小二乘法應用的廣泛性,這里從傳統的最小二乘法出發,針對其不足,提出了新的基于實對稱矩陣 UTDU分解的改進算法.該算法一方面避免了由于矩陣求逆帶來的計算量大而使得定位解算實時性降低的問題,可實現 GDOP值和定位求解的快速計算,降低計算復雜度;另一方面,分解方法能有效降低求解矩陣的條件數,改善計算過程的數值穩定性.

1 定位求解與精度因子計算

在衛星導航定位系統中,基于偽距定位的測量方程如下[1-3]:

其中,y表示偽距預測值與測量值之差,y∈Rn,n表示可見衛星數;x表示狀態變量的增量,x∈Rm,m表示狀態量維數,有 m=3+nsys,nsys表示衛星系統個數,狀態變量的前三個元素是接收機三維位置坐標,其它元素表示接收機鐘差;ε是測量誤差矢量;H是 x和 y之間的線性關聯矩陣,也稱為方向余弦矩陣,H∈Rn×m.例如,對于雙星座組合導航定位系統,H具有如下形式:

其中,axj,ayj,azj表示第 j顆衛星的方向余弦.

采用最小二乘法進行定位求解時,可得 x與y的關系:

稱 M=HTH為測量矩陣,M∈Rm×m.

幾何精度因子 GDOP的表達式為

其中,VGDOP為 GDOP之值;tr(·)為求矩陣之跡.

由式(3)、式(4)可見,不管是進行定位求解還是進行精度因子的計算,如果采用常規方法,則對測量矩陣 M進行求逆運算是不可避免的.由于進行定位求解時需要進行迭代運算,且在選星過程中,也要對多種衛星組合方案的 GDOP值進行求解運算,其計算量很大.特別是在多星座組合導航定位系統中,因為測量矩陣 M的階數會隨著星座個數的增加而增加,而矩陣求逆中乘法和加法的運算量約與其階數的三次方成正比,因此在采用常規方法進行計算時,其計算負荷更顯突出,也會占用較多的時間,這對動態用戶特別是高動態用戶接收機而言是一個非常嚴峻的挑戰.另外,從式(4)還可以看出,如果只是進行精度因子的計算,并不需要知道測量矩陣的逆矩陣中的全部元素,只需求出其主對角線元素即可,這就促使人們采取一些改進算法來達到該目的,以獲得運算量的降低.文獻[7]提出了一種計算 GDOP的閉合公式,但該方法只適合單星座 4星狀態下的求解計算.文獻[8]以矩陣特征值與特征多項式的關系為基礎提出了一種可有效降低 GDOP求解運算量的方法,但它是針對單星座系統進行的,并且需要計算矩陣 M的行列式.文獻[9]以方向余弦矩陣 H的 QR分解為依據提出了一種計算 GDOP的改進算法,可知,QR分解的計算量很大,使用Gram-Schmidt方法對一個 n階方陣進行 QR分解大概需要 n3次乘法(是 LU分解的 3倍),并有大約相同次數的加法[8].而對一個 n×m維矩陣進行 QR分解,在衛星系統數確定、即 m一定的情形下,其乘法和加法運算次數均與 n2成正比,從而,當可見衛星數較多,即 n較大時,該方法的運算量也比較大.另外,在衛星導航定位系統中采用最小二乘法進行用戶速度測量時,也采用與式(3)同樣的方式進行求解,即速度與位置求解的測量矩陣是相同的.對于有相同系數矩陣的問題,LU分解是一種有效處理方法[10].因此,綜合考慮以上因素,本文以 LU分解為基礎,提出了一種進行定位求解和精度因子計算的改進方法.同時,考慮到測量矩陣 M的對稱性,這里實際采用的是基于對稱矩陣的 LU分解,即 UTDU分解.

2 衛星定位與 DOP的改進方法

在衛星導航定位系統的實際使用中,一般來說,測量矩陣 M是對稱正定矩陣[3],由線性代數理論可知,它可進行 UTDU分解,且分解是唯一的[9],即有

其中,U為單位上三角矩陣;D為對角矩陣;它們的階數與矩陣 M相同,即 U,D∈ Rm×m.矩陣 U,D形式如下:

設矩陣 M中的元素為 mij(i,j=1,2,…,m),通過推導,可得其 UTDU分解算法:

從而,定位求解方程式(3)可變換為

其中,b=HTy.求出 b后,按回代求解方式[8]可獲得 x的值,這樣,即可避免迭代求解過程中的矩陣求逆運算.對于矩陣 UTDU分解的計算量,根據式(8),可推導出如下結論:

m階實對稱矩陣 UTDU分解算法中,乘法(含除法,下同)的運算次數為(m3/3+m2/2-5m/6),加法的運算次數為(m3-m)/6.

而 m階方陣求逆算法中,其乘法和加法的運算次數均約為 m3.例如,對雙星座組合導航系統而言,有 m=5,矩陣 M的 UTDU分解算法中乘法次數為 50,加法次數為 20,而求逆計算時,乘法和加法次數均約為 125次.顯然,采用 UTDU分解方式可節省約 60%的乘法和 80%的加法運算,其計算量的改善值是很大的.而采用式(9)的回代求解方法中(不考慮 HTy的計算量),其乘法運算次數為 m2,加法運算次數為(m2-m),與采用求逆后再進行直接求解的方法相比,并不會增加計算量.由此可見,采用 UTDU分解方法進行求解可有效降低運算量.

根據分解結果可得

可以證明,單位上三角矩陣的逆仍為單位上三角矩陣,故可令

經推導,可得到計算 U-1中元素的算法,有

在該算法中,乘法和加法的運算次數均為(m3/6-m2/2+m/3).

由式(4)可知,計算 GDOP時只需要求出矩陣 M-1的主對角線元素,根據式(10),即可推導出求解精度因子的算法:

設矩陣 M-1的主對角線元素為 δi(i=1,2,…,m),則有

該計算式中,乘法運算次數為(m2-m),加法運算量為(m2-m)/2.根據式(4)即可得到GDOP的計算式,有

對于其它類型的精度因子,由 δi的值即可根據各精度因子的定義[1]進行計算,如:

其中,VHDOP為水平精度因子之值;VVDOP為垂向精度因子之值;VPDOP為位置精度因子之值;VTDOP為時間精度因子之值.由此可知,在進行精度因子計算時,只有 GDOP是與測量矩陣 M的逆的跡相對應,而對其它精度因子,只有得到矩陣 M-1的各主對角線元素之后才能進行計算.因此,在衛星導航系統中,也就難以采用通用的求矩陣之逆的跡的方法來進行各精度因子計算.

由上面的分析可知,基于 UTDU分解,通過式(12)、式(13)得到 m階矩陣 M-1的主對角線元素的方法,其計算量為(m3+2m2-3m)/2次乘法,(m3-m)/3次加法.通過矩陣求逆方式獲取m階矩陣 M-1的主對角線元素約需要 m3的乘法和基本相同次數的加法運算,因此采用本文所述改進方法能有效降低精度因子的計算量,特別是在多星座組合導航定位系統中,矩陣 M的階數 m是隨著星座數的增加而增加的,從而其計算量的降低量將更為明顯.例如,在雙星座組合導航定位系統中,m=5,采用本文所述改進方法計算矩陣M-1的主對角線元素時,需要 80次乘法和 40次加法運算;而采用直接求逆方式時,則需要約 125次乘法和加法,以乘法來衡量,降低了約 36%的運算量.當然,采用直接求逆方式來進行精度因子求解,125次乘法運算量并不大,采用改進方法節省 45次乘法運算似乎意義也不大.的確,對于計算一次 DOP而言,其改進意義是不很明顯,但在諸如選星求解等過程中,需要在一次求解中進行成百上千次 DOP計算的情形下,其計算量的改善意義就可得到充分體現.例如,如果在某次選星過程中需要進行 100次 DOP計算(事實上,這是很普遍的,通常的選星方法其 DOP計算次數會比這大得多),則改進方法可節省約 4 500次乘法運算和約 8500次加法運算,其改善值很大.

由以上分析過程可知,在采用 UTDU分解方式進行 PVT求解和精度因子計算時,求解過程中需要形成 3個矩陣 U,U-1和 D.在實際操作過程中,為了降低存儲空間,事實上,并不需要定義 3個矩陣來存儲它們,而只需定義一個矩陣即可,例如,就采用矩陣 U來存儲這 3個矩陣,此時,矩陣U中的元素排列如下:

3 性能分析

基于 UTDU分解的方法具有嚴格的數學理論基礎,因此,本文所提改進方法的正確性和有效性是毋庸置疑的.下面對該方法的復雜性等性能進行進一步分析討論.

3.1 復雜性

由前述分析可知,對于基于 m階矩陣的UTDU分解方法,只需按式(8)計算出矩陣 D的主對角元素和矩陣 U的主對角線以上元素,其求解運算量是 O(m3/3)階的,而直接求逆方法的求解運算量是 O(m3)階的,顯然,改進方法求解更為簡潔.圖 1所示是在采用雙系統定位時,在不同可見衛星數和不同迭代次數情形下,采用本文所述 UTDU分解法相對采用式(3)所示傳統的直接求逆法相比其乘法運算量的改善百分比曲線.由圖 1可見,新方法能節省 60%以上的運算量.

圖 1 PVT求解中UTDU分解方法相對傳統直接求逆法運算量改善百分比曲線

根據前述分析可知,在 DOP求解中,在雙系統條件下,采用對測量矩陣 M進行直接求逆時約需 125次乘法運算,而采用 UTDU分解方法時只需約 80次乘法運算,可節省約 36%的運算量.文獻[9]提出采用矩陣的 QR分解來計算 DOP的改進算法,本文對基于 UTDU分解方法與基于 QR分解方法的運算量進行了比較,圖 2所示為UTDU分解法相對 QR分解法乘法運算量的改善百分比曲線.由圖 2可見,UTDU分解法比 QR分解法的計算量小得多,且隨著可見衛星數的增加,改善量會大幅增加,當可見衛星數在 13顆以上時,運算量的改善量超過 90%.

圖 2 DOP求解中UTDU分解方法相對QR分解法運算量改善百分比曲線

3.2 矩陣條件數

矩陣條件數常常用來衡量矩陣的病態性.考慮到矩陣范數的等價性,為簡單計,取矩陣范數‖·‖∞來進行條件數的計算,從而,對任意可逆矩陣 A,其條件數計算式可表示為

在衛星導航定位系統的實際使用中,有時會出現測量矩陣 M為病態的情形.采用 UTDU分解法進行求解時可避免對矩陣 M的直接運算.采用分解法按式(9)所示形式進行回代求解時,其回代過程分為 3步:

式中,bT,bx為中間變量;矩陣 D為對角陣.顯然,第 2個回代求解方程 Dbx=bT其實就是簡單的除法運算,從而,考慮分解矩陣的條件數時,只需考慮矩陣 U和 UT的條件數即可.顯然,矩陣 U和UT的病態性是相同的,故只需計算矩陣 U的條件數.

下面,以具體實例來分析測量矩陣 M及UTDU分解后矩陣的條件數.

例如,在某雙星座組合導航定位系統中,某時刻得到的測量矩陣為

經計算,可知其條件數為 fcond(M)=4.546 5×106.顯然,此時矩陣 M的病態性較嚴重.而在采用本文所述分解法后,經計算,矩陣 U的條件數為 fcond(U)=202.75,相對原矩陣條件數而言,降低了 4個量級,矩陣 U已經是良態矩陣了.這就說明,通過分解方法能獲得條件數小得多的求解矩陣,可有效改善求解矩陣的病態性,提高求解的數值穩定性.

3.3 仿真結果

為了進一步檢驗基于 UTDU分解的改進方法的有效性,本文以我國的北斗衛星導航系統(COMPASS,BeiDou/COMPASSNavigation Satellite System)和美國的全球定位系統(GPS,Global Positioning System)組成的雙星座組合導航系統為例進行了仿真分析.在實際使用中,用于表征偽距測量誤差的用戶等效距離誤差(UERE,User Equivalent Range Error)可以近似表示為零均值高斯隨機變量[1],統計表明,GPS單頻接收機的典型偽距測量誤差的標準差約為 6m[2],即用戶等效距離誤差的標準偏差 σUERE=6m.故在仿真中,以隨機方式在偽距上施加了服從高斯分布 N(0,62)的誤差.仿真結果表明,采用 UTDU分解方法和采用直接求逆方法進行 PVT求解和精度因子計算時,兩者的求解結果完全一致,但如前所述,分解方法的計算量獲得了較大改善.仿真結果統計表明,采用 UTDU分解方法進行 PVT解算時,位置偏差的均方差約為 6.11m,與輸入的偽距測量誤差基本相當,沒有放大.圖 3所示是仿真結果的位置誤差分布圖.

圖 3 位置偏差

由圖 3可見,在絕大部分情形下,位置偏差不超過 12m(即 2σUERE).統計表明,位置偏差不超過 2σUERE的時刻約為 96%,同時,解算的速度偏差也獲得了類似的仿真結果,這就進一步表明了本文所述方法的正確性和有效性.另外需要說明的是,在測量矩陣嚴重病態時,因為分解方法有效降低了求解矩陣的病態性,所以,在數據精度有限的情況下,它能獲得更為穩定的結果.當然,因為改進方法也只是最小二乘法的一種求解方法,最小二乘法的本質特征決定了其難以獲得比測量誤差更高的定位精度,要想獲得更精確的定位結果,就需要采用合適的濾波方法.

4 結束語

在衛星導航定位系統中,進行精度因子計算和采用最小二乘法進行 PVT求解時,傳統方法是基于測量矩陣的直接求逆來進行的.因為矩陣直接求逆的計算量較大,且在實際使用中,不可避免地會出現測量矩陣病態性較嚴重的情形,為了克服傳統方法的不足,利用測量矩陣的對稱正定性,以矩陣 UTDU分解為基礎,提出了一種衛星定位求解和精度因子計算的改進方法,該方法有效降低了求解的計算量,并能有效降低求解矩陣的條件數,有效避免了在定位求解中因矩陣直接求逆帶來的計算量大和數值穩定性差的問題.本文所提出的改進方法具有嚴格的數學理論基礎,保證了方法的正確性和有效性.仿真結果也進一步驗證了該方法的正確有效性.目前,該改進方法已經應用于所開發的 COMPASS和 GPS雙星座兼容接收機中,所得結果滿足了項目的要求.

References)

[1]Elliott D Kaplan,Christopher J Hegarty.GPS原理與應用[M].2版.寇艷紅,譯.北京:電子工業出版社,2007:36-472 Elliott D Kaplan,Christopher J Hegarty.Understanding GPS:principles and applications[M].2nd ed.Translated by Kou Yanhong.Beijing:Publishing House of Electronics Industry,2007:36-472(in Chinese)

[2]Pratap Misra,Per Enge.全球定位系統:信號、測量與性能[M].2版.羅鳴,等,譯.北京:電子工業出版社,2008:22-69 Pratap Misra,Per Enge.Global positioning system:signals,measurements,and performance[M].2nd ed.Translated by Luo Ming,et al.Beijing:Publishing House of Electronics Industry,2008:22-69(in Chinese)

[3]謝鋼.GPS原理與接收機設計[M].北京:電子工業出版社,2009:96-153 Xie Gang.Principles of GPS and receiver design[M].Beijing:Publishing House of Electronics Industry,2009:96-153(in Chinese)

[4]常青,柳重堪,張其善.GPS的幾何精度因子和定位解的遞推算法[J].通信學報,1998,19(12):83-88 Chang Qing,Liu Zhongkan,Zhang Qishan.The recurrence algorithm for GDOP and positioning solution in GPS[J].Journal of China Institute of Communications,1998,19(12):83-88(in Chinese)

[5]Bancroft S.An algebraic solution of the GPS equations[J].IEEE Transactions on Aerospace and Electronic Systems,1985,AES-21(1):56-59

[6]Krause L O.A direct solution to GPS-type navigation equations[J].IEEE Transactions on Aerospace and Electronic Systems,1987,AES-23(2):225-232

[7]Zhu Jijie.Calculation of geometric dilution of precision[J].IEEE Transactions on Aerospace and Electronic Systems,1992,28(3):893-895

[8]Shing H Doong.A closed-form formula for GPSGDOP computation[J].GPS Solutions,2009,13(3):183-190

[9]陳小平,滕云龍,康榮雷,等.幾何精度因子改進算法研究[J].電子科技大學學報,2008,37(增刊):27-30 Chen Xiaoping,Teng Yunlong,Kang Ronglei,et al.Study of improved arithmetric of GDOP[J].Journal of University of Electronic Science and Technology of China,2008,37(suppl):27-30(in Chinese)

[10]Timothy Sauer.數值分析[M].吳兆金,等,譯.北京:人民郵電出版社,2010:67-212 Timothy Sauer.Numerical analysis[M].Translated by Wu Zhaojin,et al.Beijing:Posts＆Telecom Press,2010:67-212(in Chinese)

[11]Steven J Leon.線性代數[M].7版.張文博,張麗靜,譯.北京:機械工業出版社,2007:286-319 Steven J Leon.Linear algebra with applications[M].7nd ed.Translated by Zhang Wenbo,Zhang Li jing.Beijing:China Machine Press,2007:286-319(in Chinese)

(編 輯 :婁 嘉)

Improved method of satellite positioning and dilution of precision

Chen Canhui Zhang Xiaolin

(School of Electronics and Information Engineering,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)

In satellite navigation system,the traditional algorithm of solving dilution of precision(DOP)and satellite positioning based on least square method is the direct matrix inverse(DMI)method.In order to overcome the disadvantages of high computational burden and poor numerical stability of traditional DMI method,an improved method of satellite positioning and DOP was presented based on the matrix UTDU decomposition,which made use of the symmetric and positive definite performance of the measurement matrix.The correctness and validity of the new method can be guaranteed by the strict mathematical theory.The numerical results show that,in comparison with the traditional DMI method,the reduction of operational volume of positioning is about 60%and that of solving DOP is about36%by the proposed method.At the same time,the condition number of the solving matrix of the improved method has reduced considerably after decomposition and the numerical stability is significantly improved.

satellite navigation;leastsquare;solutions;dilution of precision(DOP);matrix decomposition

TN 967.1

A

1001-5965(2011)04-0472-06

2010-06-08

國防科工局航天民用專項資助項目;北京市重點學科基金資助項目(XK 100070525)

陳燦輝(1973-),男,湖南汨羅人,博士生,canhuich@yahoo.com.cn.