王彬彬,易卿武,2,高 銘,盛傳貞,應俊俊,楊建雷,趙精博
(1.中國電子科技集團公司第五十四研究所 衛(wèi)星導航系統(tǒng)與裝備技術國家重點實驗室,河北 石家莊 050081;2.西安電子科技大學 電子工程學院,陜西 西安 710071;3.中國科學院精密測量科學與技術創(chuàng)新研究院 大地測量與地球動力學國家重點實驗室,湖北 武漢 430071)
全球衛(wèi)星導航系統(tǒng)是一種利用衛(wèi)星測距信號為用戶接收機提供全天候和高精度的定位、導航、授時服務(Positioning,Navigation and Timing,PNT)的現(xiàn)代化系統(tǒng),并且在重力反演、海面測高等現(xiàn)代大地測量中具有重要的研究價值[1-2]。精密單點定位(Precise Point Positioning,PPP)作為一種全球衛(wèi)星導航系統(tǒng)(Global Navigation Satellite System,GNSS)高精度定位算法,其能夠精確地估計接收機的絕對坐標[3],已被廣泛應用于地殼形變監(jiān)測、低軌衛(wèi)星定軌、空間大氣監(jiān)測和時間同步等領域[4-9]。PPP算法的難點在于數(shù)學模型中待估參數(shù)過多,達到5+N個,其中N為觀測衛(wèi)星數(shù)。因此,參數(shù)估計會對PPP的處理結果產(chǎn)生重要影響,尤其是矩陣求逆方法更是直接決定了PPP的定位精度和時效性。
矩陣求逆是矩陣論的重要內容,矩陣論中的廣義逆能夠解決工程應用中的大量離散數(shù)學問題,在數(shù)理統(tǒng)計、系統(tǒng)理論、優(yōu)化計算和控制論等多領域中有重要應用[10]。矩陣分解是一種有效的廣義逆的求解方法,常用的矩陣分解求逆方法包括正交三角分解、Cholesky分解和奇異值分解(Singular Value Decomposition,SVD)等。相關研究者在大地測量及相關領域對矩陣求逆方法開展了較為深入的研究工作。在海洋數(shù)據(jù)同化中,SVD分解、LU分解、正交三角分解的計算效率不同[11]?;诩船F(xiàn)場可編程門陣列(Field-Programmable Gate Array,FPGA)平臺的Cholesky分解求逆法具有較高的計算效率[12]。此外,矩陣分塊運算也可采用矩陣分解求逆的方法[13],其能夠提高平差解算的精度和可靠性[14]。
為了研究不同矩陣求逆方法在PPP算法中的具體計算效果,首先介紹矩陣求逆的基本概念,包括廣義逆、SVD分解求逆和Cholesky分解求逆。其次,針對高精度GNSS數(shù)據(jù)處理的復雜過程,論述PPP算法的實現(xiàn)步驟,包括建立觀測方程和參數(shù)估計方法。最后,結合iGMAS觀測站數(shù)據(jù),分別統(tǒng)計采用SVD分解求逆法和Cholesky分解求逆法的動態(tài)PPP定位精度和計算所需時間,用以驗證兩種矩陣分解求逆方法在GNSS精密單點定位算法中的應用效果。
廣義逆是線性代數(shù)中針對矩陣的一種運算[15]。一個矩陣A的廣義逆叫作A的廣義逆矩陣,具有部分逆矩陣的特性。假設一矩陣A∈m×n及另一矩陣Ag∈m×n,若Ag滿足A·Ag·A=A,則Ag即為A的廣義逆矩陣。廣義逆也稱為偽逆,是針對非可逆矩陣求解出與可逆矩陣類似特性的求逆方法。任意矩陣都存在廣義逆陣,若某一矩陣存在逆矩陣,逆矩陣即為其唯一的廣義逆陣[7]??紤]以下線性方程
A·x=y
(1)
式中:A為n×m的矩陣;y∈(A),是A的列空間。若矩陣A為可逆矩陣,則x=A-1·y是方程式的解。若矩陣A為可逆矩陣,則
A·A-1·A=A
(2)
假設矩陣A不可逆,或者是n≠m,需要一個合適的m×n矩陣G,即
A·G·A=A
(3)
因此,G·y為線性系統(tǒng)A·x=y的解。同樣,對于m×n階的矩陣G
G·A·G=G
(4)
成立。那么,廣義逆陣可定義為假設一個n×m的矩陣A,m×n的矩陣G可以使A·G·A=A成立,矩陣G即為A的廣義逆矩陣。
廣義逆矩陣理論可用于推導線性方程組的解。對于上述任何一種廣義逆陣都可以用來判斷線性方程組是否有解,若有解時列出其所有的解。若以下n×m的線性系統(tǒng)
A·x=b
(5)
有解存在。其中:向量x為未知數(shù);向量b為常數(shù)。那么,廣義逆矩陣理論得到的解為
x=Ag·b+[I-Ag·A]·w
(6)
式中:I是單位矩陣;參數(shù)w為任意矩陣,而Ag為A的任何一個廣義逆矩陣。線性系統(tǒng)解存在的條件是當且僅當Ag·b為其中一個解,也就是當且僅當A·Ag·b=b。
當矩陣為方陣時,矩陣分解是一種有效的廣義逆矩陣求解方法,通過將一個矩陣拆解為多個矩陣的乘積。常用實現(xiàn)一些矩陣的快速運算,包括三角分解、滿秩分解、QR分解、Jordan分解、Jacobi分解、SVD分解和Cholesky分解等。下面主要分析SVD分解和Cholesky分解在矩陣分解求逆中的方法和效果。
SVD分解是主成分分析、智能計算、機器學習和自然語言處理等領域的基礎[16-17]。而SVD分解求逆則是通過構造特征值和特征向量,對矩陣進行奇異值分解,從而便于矩陣求逆運算。對于矩陣M,其與特征值和特征向量的關系為
M·x=λ·x
(7)
式中:λ是矩陣M的一個特征值;x是對應的特征向量。對于n維矩陣M,如果求出所有的特征值λ1,λ2,…,λn,以及對應的特征向量{x1,x2,…,xn},那么矩陣M可分解為
M=U·Σ·V*
(8)
式中:U為特征向量{x1,x2,…,xn}構成的n維矩陣,V*是矩陣U的逆矩陣。Σ則是特征值λ1,λ2,…,λn構成的對角線矩陣,特征值也被稱為奇異值。因此,式(8)為矩陣M的奇異值分解。
奇異值分解可以被用來計算矩陣的逆矩陣。若矩陣M的奇異值分解為M=U·Σ·V*,那么矩陣M的逆矩陣為
M=V·Σ-1·U*
(9)
式中,Σ-1是Σ的逆矩陣,可通過矩陣Σ中主對角線上非零元素都求倒之后再轉置得到。
Cholesky分解求逆是指將一個正定的埃爾米特矩陣,即正定對稱矩陣分解成一個下三角矩陣與其共軛轉置之乘積。通過計算下三角矩陣的逆矩陣,快速求解正定對稱陣的逆矩陣[18]。這種分解方式在提高代數(shù)運算效率、蒙特卡羅方法等場合中十分有用。實際應用中,Cholesky分解求逆在求解線性方程組中的效率約為三角分解求逆的兩倍。
對于n×n的對稱正定矩陣A,則存在一個主對角線元素嚴格正定的下三角矩陣L,滿足
A=L·L*
(10)
式中:L是一個下三角矩陣且所有對角線元素均為正實數(shù);L*是矩陣L的共軛轉置矩陣。每一對稱正定矩陣都有一個唯一的Cholesky分解過程。Cholesky分解的計算時間復雜度為O(n3/3),正定對稱陣A的逆矩陣可由下三角矩陣求逆確定,表示為
A-1=(L·L*)-1=(L-1)*·L-1
(11)
L是通過Cholesky分解得到的下三角陣,假設P是L的逆矩陣,則有
L·P=
(12)
由式(12)可知,P也是下三角矩陣,其中主對角線的值為
(13)
矩陣P是下三角矩陣L的逆矩陣,因此,矩陣A經(jīng)過Cholesky分解后的逆矩陣為
A-1=(L·L*)-1=P*·P
(14)
PPP是一種使用單個接收機的單系統(tǒng)偽距和載波相位觀測數(shù)據(jù),采用事后的精密衛(wèi)星軌道和鐘差產(chǎn)品,通過修正觀測方程中的系統(tǒng)誤差源獲取高精度接收機坐標、接收機鐘差、大氣延遲和浮點模糊度的衛(wèi)星定位方法。PPP不需要建立同步參考站,能夠獨立完成高精度絕對定位。同時,可以計算接收機鐘差和大氣延遲量,在科研和工程領域都具有較高的應用價值,也是目前GNSS數(shù)據(jù)處理領域的研究熱點之一。
采用無電離層組合觀測方程[3],包括無電離層偽距觀測方程和無電離層載波相位觀測方程,其表達式分別為
(15)
(16)
(17)
其中,
(18)
(19)
無電離層相位組合中的模糊度參數(shù)是寬巷模糊度和窄巷模糊度的組合,并且包含硬件延遲。因此,其不具備整數(shù)特性,在參數(shù)估計時需當作浮點數(shù)進行處理。在參數(shù)估計前需要對無電離層組合觀測量進行系統(tǒng)誤差修正,包括對流層延遲干分量、接收機端和衛(wèi)星端的天線相位中心偏差(Phase Center Offset,PCO)和天線相位中心變化(Phase Center Variation,PCV)修正、相對論效應修正、相位纏繞誤差修正、地球自轉修正、固體潮修正、海潮修正和極移潮修正等。
對無電離層組合進行幾何距離線性展開、斜向對流層延遲投影到天頂方向和系統(tǒng)誤差修正處理后,將觀測方程按照矩陣形式表達。若某一時刻有m個觀測衛(wèi)星,未知參數(shù)X的表達形式為
(20)
函數(shù)模型和隨機模型統(tǒng)稱為PPP數(shù)學模型,表達式分別為
根據(jù)無電離層組合的表達式和誤差傳播定律,天頂方向相位無電離層組合先驗精度約為1 cm,偽距無電離層組合先驗精度約為1 m。斜向觀測值精度與高度角有關,高度角越小先驗誤差越大,對應的權重也就越小。
卡爾曼濾波是一種最優(yōu)化的自回歸參數(shù)估計算法,其基于一組觀測序列以及動力學模型信息,采用遞歸算法求解狀態(tài)向量的最優(yōu)估值[19]。狀態(tài)方程式和觀測方程式的數(shù)學模型分別為
其中,
卡爾曼濾波算法將時間更新和觀測更新依次迭代運行,就能夠得到所有歷元當前時刻的參數(shù)估值,是一種局部最優(yōu)解。有別于最小二乘算法,卡爾曼濾波算法不需要存儲所有觀測時刻的數(shù)據(jù),僅在當前歷元進行狀態(tài)預測和參數(shù)估計,未知參數(shù)相對較少,很大程度上節(jié)省了計算機內存并提高計算效率。另外,狀態(tài)方程能夠基于之前歷元的參數(shù)估值預測下個歷元的參數(shù)值,適用于未知參數(shù)隨時間變化的觀測方程。
國際GNSS監(jiān)測評估系統(tǒng)(International GNSS Monitoring and Assessment System,iGMAS)是對GNSS運行狀況和主要指標進行監(jiān)測和評估,生成高精度星歷、鐘差和全球電離層延遲模型等產(chǎn)品的平臺。目前,由30個跟蹤站、3個數(shù)據(jù)中心、7個分析中心、2個監(jiān)測評估中心、1個產(chǎn)品綜合與服務中心、1個運行管理控制中心和通信網(wǎng)絡組成[20-22]。iGMAS跟蹤站完成北斗衛(wèi)星導航系統(tǒng)(BeiDou Navigation Satellite System,BDS)、全球定位系統(tǒng)(Global Positioning System,GPS)、全球衛(wèi)星導航系統(tǒng)(Global Navigation Satellite System,GLONASS)和Galileo等GNSS的信號接收和測量、原始觀測數(shù)據(jù)的采集,并進行數(shù)據(jù)合理性檢驗和數(shù)據(jù)預處理,并將數(shù)據(jù)發(fā)送至數(shù)據(jù)中心備份。中國電子科技集團公司第五十四研究所長期穩(wěn)定運行一個監(jiān)測評估中心,研發(fā)多頻多系統(tǒng)GNSS監(jiān)測接收機,并裝配16個iGMAS跟蹤站,iGMAS跟蹤站名稱及全球分布情況如表1所示。
表1 iGMAS跟蹤站名稱及全球分布
為了分析不同矩陣求逆方法在GNSS精密單點定位參數(shù)估計中的應用效果,采用iGMAS跟蹤站中的BJF1站在2020年第96天的北斗三號和北斗二號觀測數(shù)據(jù),進行PPP數(shù)據(jù)處理驗證和處理性能分析。BJF1站位于北京,經(jīng)緯度為39.6 °N和115.9 °E,接收機主機類型為CETC-54-GMR-4016,接收機天線類型為LEIAR25.R4。數(shù)據(jù)處理平臺為聯(lián)想E580 PC機,處理器為Intel Core i5-8250 U,主頻為1.6 GHz。數(shù)據(jù)處理軟件是在開源庫GPSTk基礎上開發(fā)的GNSS Data Processor軟件。在動態(tài)PPP數(shù)據(jù)處理實驗中,將接收機坐標參數(shù)作為白噪聲進行估計。
PPP數(shù)據(jù)處理流程主要包括數(shù)據(jù)準備、數(shù)據(jù)預處理和參數(shù)估計。準備的數(shù)據(jù)文件包括精密星歷文件和觀測文件,其中精密星歷文件可以通過FTP協(xié)議從IGS網(wǎng)站獲取。將星歷數(shù)據(jù)文件轉換為SP3標準格式,觀測數(shù)據(jù)文件轉換為RINEX標準格式,完成數(shù)據(jù)準備工作。數(shù)據(jù)預處理主要包括周跳探測、粗差處理和誤差源模型化修正。采用TurboEdit方法對周跳進行探測[23],出現(xiàn)周跳的相位觀測值通過調整權重進行處理,使用經(jīng)驗模型分別修正,包括衛(wèi)星天線和接收機天線誤差、相對論效應、相位纏繞誤差、地球自轉、固體潮、海潮和極潮等誤差。參數(shù)估計部分包括建立函數(shù)模型和隨機模型,并采用卡爾曼濾波算法依次對每個歷元的觀測值進行參數(shù)估計。卡爾曼濾波與權函數(shù)調整迭代運行,直到所有粗差被處理完成為止。
接下來分析SVD分解求逆與Choleskey分解求逆的計算速度。分別采用以上這兩種方案處理BJF1站同一天的觀測數(shù)據(jù)。在方案1中,采用SVD分解法處理PPP卡爾曼濾波中法方程求逆問題。在方案2中,則使用Choleskey分解法,求解相應的逆矩陣。統(tǒng)計兩種方案中每個歷元所需要的時間,包括讀取文件、數(shù)據(jù)預處理和矩陣運算,其中矩陣求逆是主要項。兩種方案的計算時間分別如圖1和圖2所示。
圖1 動態(tài)PPP卡爾曼濾波法方程SVD分解求逆計算時間
圖2 動態(tài)PPP卡爾曼濾波法方程Cholesky分解求逆計算時間
由圖1和圖2可知,在一天所處理的2 500個觀測歷元中,SVD分解法所需計算時間變化較大,均值為291 ms。而Cholesky分解法所需時間比較穩(wěn)定,均值為189 ms。兩種方案均在第188個歷元所需時間較長,通過處理日志判斷該歷元觀測值粗差較多,需要多次矩陣求逆運算,導致該歷元計算所需時間較長。
不同的矩陣分解方法求逆精度不一致,這會通過卡爾曼濾波參數(shù)估計從而影響PPP定位結果。分別分析SVD分解求逆法和Cholesky分解求逆法對BJF1站PPP定位結果的影響,對比分析方案1和方案2中BJF1站定位精度。以采用事后精密星歷,靜態(tài)PPP處理模式得到的BJF1坐標為參考值,將方案1和方案2中的定位結果分別與參考值進行對比,并將位置定位偏差轉換到北方向(North,N)、東方向(East,E)、垂向(Up,U)方向,結果分別如圖3和圖4所示。
圖3 動態(tài)PPP卡爾曼濾波SVD分解求逆定位結果
圖4 動態(tài)PPP卡爾曼濾波Cholesky分解求逆定位結果
在采用SVD分解求逆法的方案1所中,BJF1站動態(tài)PPP定位偏差在45 min后首次收斂到0.2 m以內,并在之后的時間段內變化較大,尤其是U方向,抖動劇烈。而在使用Choleskey分解法求逆的方案2中,BJF1站動態(tài)PPP定位偏差在12 min就收斂到 0.2 m以內,整個時間段內的定位偏差也表現(xiàn)穩(wěn)定。經(jīng)統(tǒng)計,采用SVD分解法的動態(tài)PPP在N、E、U三維方向的定位精度分別為5.95 cm、8.04 cm、20.31 cm,而使用Cholesky分解法的定位精度則分別是2.65 cm、4.05 cm、8.04 cm。
由此可知,無論是在計算所需時間方面,還是在因矩陣求逆誤差而導致的定位偏差方面,Cholesky分解法均優(yōu)于SVD分解法,具體的SVD分解和Cholesky分解效果對比情況如表2所示。
表2 SVD分解和Cholesky分解效果對比
矩陣分解求逆在工程應用中具有重要價值。論述了矩陣求逆的基本理論和方法,介紹矩陣廣義逆的性質及其在線性方程解中的作用。重點論述了兩種矩陣分解求逆方法,包括SVD分解求逆法和Cholesky分解求逆法。其中,SVD分解法可用于求解廣義逆和方陣求逆,而Cholesky分解法適用于正定對稱矩陣求逆。精密單點定位是一種GNSS高精度定位算法,參數(shù)估計是PPP數(shù)據(jù)處理的重要環(huán)節(jié),可通過矩陣分解進行快速求逆運算。描述了PPP算法觀測方程,介紹觀測方程中各參數(shù)的幾何和物理含義,包括各種誤差項的具體處理策略,并論述一種卡爾曼濾波參數(shù)估計算法。
采用iGMAS中的BJF1站北斗二號和北斗三號觀測數(shù)據(jù)進行PPP處理,驗證SVD分解法和Cholesky分解法的矩陣求逆效果?;赟VD分解法的動態(tài)PPP每個歷元數(shù)據(jù)處理所需時間均值為291 ms,N、E、U三維方向定位誤差分別為5.95 cm、8.04 cm、20.31 cm。使用Cholesky分解法的動態(tài)PPP每個歷元數(shù)據(jù)處理所需時間為189 ms,三維方向定位誤差分別為2.65 cm、4.05 cm、8.04 cm。因此,相比于SVD分解法,Cholesky分解法應用于矩陣求逆的時效性更好,動態(tài)定位精度統(tǒng)計結果也間接證明了Cholesky分解法矩陣求逆精度更高。這是因為PPP卡爾曼濾波法方程為對稱正定矩陣,在矩陣求逆時Cholesky分解法更為適用。