基于SVD和Cholesky求逆方法的精密單點(diǎn)定位研究

王彬彬，易卿武,2，高 銘，盛傳貞，應(yīng)俊俊，楊建雷，趙精博

(1.中國電子科技集團(tuán)公司第五十四研究所 衛(wèi)星導(dǎo)航系統(tǒng)與裝備技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，河北 石家莊 050081;2.西安電子科技大學(xué) 電子工程學(xué)院，陜西 西安 710071;3.中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院 大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430071)

全球衛(wèi)星導(dǎo)航系統(tǒng)是一種利用衛(wèi)星測距信號為用戶接收機(jī)提供全天候和高精度的定位、導(dǎo)航、授時(shí)服務(wù)(Positioning,Navigation and Timing,PNT)的現(xiàn)代化系統(tǒng)，并且在重力反演、海面測高等現(xiàn)代大地測量中具有重要的研究價(jià)值[1-2]。精密單點(diǎn)定位(Precise Point Positioning,PPP)作為一種全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System,GNSS)高精度定位算法，其能夠精確地估計(jì)接收機(jī)的絕對坐標(biāo)[3]，已被廣泛應(yīng)用于地殼形變監(jiān)測、低軌衛(wèi)星定軌、空間大氣監(jiān)測和時(shí)間同步等領(lǐng)域[4-9]。PPP算法的難點(diǎn)在于數(shù)學(xué)模型中待估參數(shù)過多，達(dá)到5+N個(gè)，其中N為觀測衛(wèi)星數(shù)。因此，參數(shù)估計(jì)會(huì)對PPP的處理結(jié)果產(chǎn)生重要影響，尤其是矩陣求逆方法更是直接決定了PPP的定位精度和時(shí)效性。

矩陣求逆是矩陣論的重要內(nèi)容，矩陣論中的廣義逆能夠解決工程應(yīng)用中的大量離散數(shù)學(xué)問題，在數(shù)理統(tǒng)計(jì)、系統(tǒng)理論、優(yōu)化計(jì)算和控制論等多領(lǐng)域中有重要應(yīng)用[10]。矩陣分解是一種有效的廣義逆的求解方法，常用的矩陣分解求逆方法包括正交三角分解、Cholesky分解和奇異值分解(Singular Value Decomposition，SVD)等。相關(guān)研究者在大地測量及相關(guān)領(lǐng)域?qū)仃嚽竽娣椒ㄩ_展了較為深入的研究工作。在海洋數(shù)據(jù)同化中，SVD分解、LU分解、正交三角分解的計(jì)算效率不同[11]。基于即現(xiàn)場可編程門陣列(Field-Programmable Gate Array,FPGA)平臺(tái)的Cholesky分解求逆法具有較高的計(jì)算效率[12]。此外，矩陣分塊運(yùn)算也可采用矩陣分解求逆的方法[13]，其能夠提高平差解算的精度和可靠性[14]。

為了研究不同矩陣求逆方法在PPP算法中的具體計(jì)算效果，首先介紹矩陣求逆的基本概念，包括廣義逆、SVD分解求逆和Cholesky分解求逆。其次，針對高精度GNSS數(shù)據(jù)處理的復(fù)雜過程，論述PPP算法的實(shí)現(xiàn)步驟，包括建立觀測方程和參數(shù)估計(jì)方法。最后，結(jié)合iGMAS觀測站數(shù)據(jù)，分別統(tǒng)計(jì)采用SVD分解求逆法和Cholesky分解求逆法的動(dòng)態(tài)PPP定位精度和計(jì)算所需時(shí)間，用以驗(yàn)證兩種矩陣分解求逆方法在GNSS精密單點(diǎn)定位算法中的應(yīng)用效果。

1 矩陣求逆

1.1 廣義逆

廣義逆是線性代數(shù)中針對矩陣的一種運(yùn)算[15]。一個(gè)矩陣A的廣義逆叫作A的廣義逆矩陣，具有部分逆矩陣的特性。假設(shè)一矩陣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)

假設(shè)矩陣A不可逆，或者是n≠m，需要一個(gè)合適的m×n矩陣G，即

A·G·A=A

(3)

因此，G·y為線性系統(tǒng)A·x=y的解。同樣，對于m×n階的矩陣G

G·A·G=G

(4)

成立。那么，廣義逆陣可定義為假設(shè)一個(gè)n×m的矩陣A，m×n的矩陣G可以使A·G·A=A成立，矩陣G即為A的廣義逆矩陣。

廣義逆矩陣?yán)碚摽捎糜谕茖?dǎo)線性方程組的解。對于上述任何一種廣義逆陣都可以用來判斷線性方程組是否有解，若有解時(shí)列出其所有的解。若以下n×m的線性系統(tǒng)

A·x=b

(5)

有解存在。其中：向量x為未知數(shù)；向量b為常數(shù)。那么，廣義逆矩陣?yán)碚摰玫降慕鉃?/p>

x=Ag·b+[I-Ag·A]·w

(6)

式中：I是單位矩陣；參數(shù)w為任意矩陣，而Ag為A的任何一個(gè)廣義逆矩陣。線性系統(tǒng)解存在的條件是當(dāng)且僅當(dāng)Ag·b為其中一個(gè)解，也就是當(dāng)且僅當(dāng)A·Ag·b=b。

當(dāng)矩陣為方陣時(shí)，矩陣分解是一種有效的廣義逆矩陣求解方法，通過將一個(gè)矩陣拆解為多個(gè)矩陣的乘積。常用實(shí)現(xiàn)一些矩陣的快速運(yùn)算，包括三角分解、滿秩分解、QR分解、Jordan分解、Jacobi分解、SVD分解和Cholesky分解等。下面主要分析SVD分解和Cholesky分解在矩陣分解求逆中的方法和效果。

1.2 SVD分解求逆

SVD分解是主成分分析、智能計(jì)算、機(jī)器學(xué)習(xí)和自然語言處理等領(lǐng)域的基礎(chǔ)[16-17]。而SVD分解求逆則是通過構(gòu)造特征值和特征向量，對矩陣進(jìn)行奇異值分解，從而便于矩陣求逆運(yùn)算。對于矩陣M，其與特征值和特征向量的關(guān)系為

M·x=λ·x

(7)

式中：λ是矩陣M的一個(gè)特征值；x是對應(yīng)的特征向量。對于n維矩陣M，如果求出所有的特征值λ1,λ2,…,λn，以及對應(yīng)的特征向量{x1,x2,…,xn}，那么矩陣M可分解為

M=U·Σ·V*

(8)

式中：U為特征向量{x1,x2,…,xn}構(gòu)成的n維矩陣，V*是矩陣U的逆矩陣。Σ則是特征值λ1,λ2,…,λn構(gòu)成的對角線矩陣，特征值也被稱為奇異值。因此，式(8)為矩陣M的奇異值分解。

奇異值分解可以被用來計(jì)算矩陣的逆矩陣。若矩陣M的奇異值分解為M=U·Σ·V*，那么矩陣M的逆矩陣為

M=V·Σ-1·U*

(9)

式中，Σ-1是Σ的逆矩陣，可通過矩陣Σ中主對角線上非零元素都求倒之后再轉(zhuǎn)置得到。

1.3 Cholesky分解求逆

Cholesky分解求逆是指將一個(gè)正定的埃爾米特矩陣，即正定對稱矩陣分解成一個(gè)下三角矩陣與其共軛轉(zhuǎn)置之乘積。通過計(jì)算下三角矩陣的逆矩陣，快速求解正定對稱陣的逆矩陣[18]。這種分解方式在提高代數(shù)運(yùn)算效率、蒙特卡羅方法等場合中十分有用。實(shí)際應(yīng)用中，Cholesky分解求逆在求解線性方程組中的效率約為三角分解求逆的兩倍。

對于n×n的對稱正定矩陣A，則存在一個(gè)主對角線元素嚴(yán)格正定的下三角矩陣L，滿足

A=L·L*

(10)

式中：L是一個(gè)下三角矩陣且所有對角線元素均為正實(shí)數(shù)；L*是矩陣L的共軛轉(zhuǎn)置矩陣。每一對稱正定矩陣都有一個(gè)唯一的Cholesky分解過程。Cholesky分解的計(jì)算時(shí)間復(fù)雜度為O(n3/3)，正定對稱陣A的逆矩陣可由下三角矩陣求逆確定，表示為

A-1=(L·L*)-1=(L-1)*·L-1

(11)

L是通過Cholesky分解得到的下三角陣，假設(shè)P是L的逆矩陣，則有

L·P=

(12)

由式(12)可知，P也是下三角矩陣，其中主對角線的值為

(13)

矩陣P是下三角矩陣L的逆矩陣，因此，矩陣A經(jīng)過Cholesky分解后的逆矩陣為

A-1=(L·L*)-1=P*·P

(14)

2 精密單點(diǎn)定位

PPP是一種使用單個(gè)接收機(jī)的單系統(tǒng)偽距和載波相位觀測數(shù)據(jù)，采用事后的精密衛(wèi)星軌道和鐘差產(chǎn)品，通過修正觀測方程中的系統(tǒng)誤差源獲取高精度接收機(jī)坐標(biāo)、接收機(jī)鐘差、大氣延遲和浮點(diǎn)模糊度的衛(wèi)星定位方法。PPP不需要建立同步參考站，能夠獨(dú)立完成高精度絕對定位。同時(shí)，可以計(jì)算接收機(jī)鐘差和大氣延遲量，在科研和工程領(lǐng)域都具有較高的應(yīng)用價(jià)值，也是目前GNSS數(shù)據(jù)處理領(lǐng)域的研究熱點(diǎn)之一。

2.1 觀測方程

采用無電離層組合觀測方程[3]，包括無電離層偽距觀測方程和無電離層載波相位觀測方程，其表達(dá)式分別為

(15)

(16)

(17)

其中，

(18)

(19)

無電離層相位組合中的模糊度參數(shù)是寬巷模糊度和窄巷模糊度的組合，并且包含硬件延遲。因此，其不具備整數(shù)特性，在參數(shù)估計(jì)時(shí)需當(dāng)作浮點(diǎn)數(shù)進(jìn)行處理。在參數(shù)估計(jì)前需要對無電離層組合觀測量進(jìn)行系統(tǒng)誤差修正，包括對流層延遲干分量、接收機(jī)端和衛(wèi)星端的天線相位中心偏差(Phase Center Offset,PCO)和天線相位中心變化(Phase Center Variation,PCV)修正、相對論效應(yīng)修正、相位纏繞誤差修正、地球自轉(zhuǎn)修正、固體潮修正、海潮修正和極移潮修正等。

對無電離層組合進(jìn)行幾何距離線性展開、斜向?qū)α鲗友舆t投影到天頂方向和系統(tǒng)誤差修正處理后，將觀測方程按照矩陣形式表達(dá)。若某一時(shí)刻有m個(gè)觀測衛(wèi)星，未知參數(shù)X的表達(dá)形式為

(20)

函數(shù)模型和隨機(jī)模型統(tǒng)稱為PPP數(shù)學(xué)模型，表達(dá)式分別為

根據(jù)無電離層組合的表達(dá)式和誤差傳播定律，天頂方向相位無電離層組合先驗(yàn)精度約為1 cm，偽距無電離層組合先驗(yàn)精度約為1 m。斜向觀測值精度與高度角有關(guān)，高度角越小先驗(yàn)誤差越大，對應(yīng)的權(quán)重也就越小。

2.2 參數(shù)估計(jì)

卡爾曼濾波是一種最優(yōu)化的自回歸參數(shù)估計(jì)算法，其基于一組觀測序列以及動(dòng)力學(xué)模型信息，采用遞歸算法求解狀態(tài)向量的最優(yōu)估值[19]。狀態(tài)方程式和觀測方程式的數(shù)學(xué)模型分別為

其中，

卡爾曼濾波算法將時(shí)間更新和觀測更新依次迭代運(yùn)行，就能夠得到所有歷元當(dāng)前時(shí)刻的參數(shù)估值，是一種局部最優(yōu)解。有別于最小二乘算法，卡爾曼濾波算法不需要存儲(chǔ)所有觀測時(shí)刻的數(shù)據(jù)，僅在當(dāng)前歷元進(jìn)行狀態(tài)預(yù)測和參數(shù)估計(jì)，未知參數(shù)相對較少，很大程度上節(jié)省了計(jì)算機(jī)內(nèi)存并提高計(jì)算效率。另外，狀態(tài)方程能夠基于之前歷元的參數(shù)估值預(yù)測下個(gè)歷元的參數(shù)值，適用于未知參數(shù)隨時(shí)間變化的觀測方程。

3 實(shí)驗(yàn)分析

國際GNSS監(jiān)測評估系統(tǒng)(International GNSS Monitoring and Assessment System,iGMAS)是對GNSS運(yùn)行狀況和主要指標(biāo)進(jìn)行監(jiān)測和評估，生成高精度星歷、鐘差和全球電離層延遲模型等產(chǎn)品的平臺(tái)。目前，由30個(gè)跟蹤站、3個(gè)數(shù)據(jù)中心、7個(gè)分析中心、2個(gè)監(jiān)測評估中心、1個(gè)產(chǎn)品綜合與服務(wù)中心、1個(gè)運(yùn)行管理控制中心和通信網(wǎng)絡(luò)組成[20-22]。iGMAS跟蹤站完成北斗衛(wèi)星導(dǎo)航系統(tǒng)(BeiDou Navigation Satellite System，BDS)、全球定位系統(tǒng)(Global Positioning System，GPS)、全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System，GLONASS)和Galileo等GNSS的信號接收和測量、原始觀測數(shù)據(jù)的采集，并進(jìn)行數(shù)據(jù)合理性檢驗(yàn)和數(shù)據(jù)預(yù)處理，并將數(shù)據(jù)發(fā)送至數(shù)據(jù)中心備份。中國電子科技集團(tuán)公司第五十四研究所長期穩(wěn)定運(yùn)行一個(gè)監(jiān)測評估中心，研發(fā)多頻多系統(tǒng)GNSS監(jiān)測接收機(jī)，并裝配16個(gè)iGMAS跟蹤站，iGMAS跟蹤站名稱及全球分布情況如表1所示。

表1 iGMAS跟蹤站名稱及全球分布

為了分析不同矩陣求逆方法在GNSS精密單點(diǎn)定位參數(shù)估計(jì)中的應(yīng)用效果，采用iGMAS跟蹤站中的BJF1站在2020年第96天的北斗三號和北斗二號觀測數(shù)據(jù)，進(jìn)行PPP數(shù)據(jù)處理驗(yàn)證和處理性能分析。BJF1站位于北京，經(jīng)緯度為39.6 °N和115.9 °E，接收機(jī)主機(jī)類型為CETC-54-GMR-4016，接收機(jī)天線類型為LEIAR25.R4。數(shù)據(jù)處理平臺(tái)為聯(lián)想E580 PC機(jī)，處理器為Intel Core i5-8250 U，主頻為1.6 GHz。數(shù)據(jù)處理軟件是在開源庫GPSTk基礎(chǔ)上開發(fā)的GNSS Data Processor軟件。在動(dòng)態(tài)PPP數(shù)據(jù)處理實(shí)驗(yàn)中，將接收機(jī)坐標(biāo)參數(shù)作為白噪聲進(jìn)行估計(jì)。

PPP數(shù)據(jù)處理流程主要包括數(shù)據(jù)準(zhǔn)備、數(shù)據(jù)預(yù)處理和參數(shù)估計(jì)。準(zhǔn)備的數(shù)據(jù)文件包括精密星歷文件和觀測文件，其中精密星歷文件可以通過FTP協(xié)議從IGS網(wǎng)站獲取。將星歷數(shù)據(jù)文件轉(zhuǎn)換為SP3標(biāo)準(zhǔn)格式，觀測數(shù)據(jù)文件轉(zhuǎn)換為RINEX標(biāo)準(zhǔn)格式，完成數(shù)據(jù)準(zhǔn)備工作。數(shù)據(jù)預(yù)處理主要包括周跳探測、粗差處理和誤差源模型化修正。采用TurboEdit方法對周跳進(jìn)行探測[23]，出現(xiàn)周跳的相位觀測值通過調(diào)整權(quán)重進(jìn)行處理，使用經(jīng)驗(yàn)?zāi)Ｐ头謩e修正，包括衛(wèi)星天線和接收機(jī)天線誤差、相對論效應(yīng)、相位纏繞誤差、地球自轉(zhuǎn)、固體潮、海潮和極潮等誤差。參數(shù)估計(jì)部分包括建立函數(shù)模型和隨機(jī)模型，并采用卡爾曼濾波算法依次對每個(gè)歷元的觀測值進(jìn)行參數(shù)估計(jì)。卡爾曼濾波與權(quán)函數(shù)調(diào)整迭代運(yùn)行，直到所有粗差被處理完成為止。

接下來分析SVD分解求逆與Choleskey分解求逆的計(jì)算速度。分別采用以上這兩種方案處理BJF1站同一天的觀測數(shù)據(jù)。在方案1中，采用SVD分解法處理PPP卡爾曼濾波中法方程求逆問題。在方案2中，則使用Choleskey分解法，求解相應(yīng)的逆矩陣。統(tǒng)計(jì)兩種方案中每個(gè)歷元所需要的時(shí)間，包括讀取文件、數(shù)據(jù)預(yù)處理和矩陣運(yùn)算，其中矩陣求逆是主要項(xiàng)。兩種方案的計(jì)算時(shí)間分別如圖1和圖2所示。

圖1 動(dòng)態(tài)PPP卡爾曼濾波法方程SVD分解求逆計(jì)算時(shí)間

圖2 動(dòng)態(tài)PPP卡爾曼濾波法方程Cholesky分解求逆計(jì)算時(shí)間

由圖1和圖2可知，在一天所處理的2 500個(gè)觀測歷元中，SVD分解法所需計(jì)算時(shí)間變化較大，均值為291 ms。而Cholesky分解法所需時(shí)間比較穩(wěn)定，均值為189 ms。兩種方案均在第188個(gè)歷元所需時(shí)間較長，通過處理日志判斷該歷元觀測值粗差較多，需要多次矩陣求逆運(yùn)算，導(dǎo)致該歷元計(jì)算所需時(shí)間較長。

不同的矩陣分解方法求逆精度不一致，這會(huì)通過卡爾曼濾波參數(shù)估計(jì)從而影響PPP定位結(jié)果。分別分析SVD分解求逆法和Cholesky分解求逆法對BJF1站PPP定位結(jié)果的影響，對比分析方案1和方案2中BJF1站定位精度。以采用事后精密星歷，靜態(tài)PPP處理模式得到的BJF1坐標(biāo)為參考值，將方案1和方案2中的定位結(jié)果分別與參考值進(jìn)行對比，并將位置定位偏差轉(zhuǎn)換到北方向(North,N)、東方向(East,E)、垂向(Up,U)方向，結(jié)果分別如圖3和圖4所示。

圖3 動(dòng)態(tài)PPP卡爾曼濾波SVD分解求逆定位結(jié)果

圖4 動(dòng)態(tài)PPP卡爾曼濾波Cholesky分解求逆定位結(jié)果

在采用SVD分解求逆法的方案1所中，BJF1站動(dòng)態(tài)PPP定位偏差在45 min后首次收斂到0.2 m以內(nèi)，并在之后的時(shí)間段內(nèi)變化較大，尤其是U方向，抖動(dòng)劇烈。而在使用Choleskey分解法求逆的方案2中，BJF1站動(dòng)態(tài)PPP定位偏差在12 min就收斂到 0.2 m以內(nèi)，整個(gè)時(shí)間段內(nèi)的定位偏差也表現(xiàn)穩(wěn)定。經(jīng)統(tǒng)計(jì)，采用SVD分解法的動(dòng)態(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。

由此可知，無論是在計(jì)算所需時(shí)間方面，還是在因矩陣求逆誤差而導(dǎo)致的定位偏差方面，Cholesky分解法均優(yōu)于SVD分解法，具體的SVD分解和Cholesky分解效果對比情況如表2所示。

表2 SVD分解和Cholesky分解效果對比

4 結(jié)語

矩陣分解求逆在工程應(yīng)用中具有重要價(jià)值。論述了矩陣求逆的基本理論和方法，介紹矩陣廣義逆的性質(zhì)及其在線性方程解中的作用。重點(diǎn)論述了兩種矩陣分解求逆方法，包括SVD分解求逆法和Cholesky分解求逆法。其中，SVD分解法可用于求解廣義逆和方陣求逆，而Cholesky分解法適用于正定對稱矩陣求逆。精密單點(diǎn)定位是一種GNSS高精度定位算法，參數(shù)估計(jì)是PPP數(shù)據(jù)處理的重要環(huán)節(jié)，可通過矩陣分解進(jìn)行快速求逆運(yùn)算。描述了PPP算法觀測方程，介紹觀測方程中各參數(shù)的幾何和物理含義，包括各種誤差項(xiàng)的具體處理策略，并論述一種卡爾曼濾波參數(shù)估計(jì)算法。

采用iGMAS中的BJF1站北斗二號和北斗三號觀測數(shù)據(jù)進(jìn)行PPP處理，驗(yàn)證SVD分解法和Cholesky分解法的矩陣求逆效果?；赟VD分解法的動(dòng)態(tài)PPP每個(gè)歷元數(shù)據(jù)處理所需時(shí)間均值為291 ms，N、E、U三維方向定位誤差分別為5.95 cm、8.04 cm、20.31 cm。使用Cholesky分解法的動(dòng)態(tài)PPP每個(gè)歷元數(shù)據(jù)處理所需時(shí)間為189 ms，三維方向定位誤差分別為2.65 cm、4.05 cm、8.04 cm。因此，相比于SVD分解法，Cholesky分解法應(yīng)用于矩陣求逆的時(shí)效性更好，動(dòng)態(tài)定位精度統(tǒng)計(jì)結(jié)果也間接證明了Cholesky分解法矩陣求逆精度更高。這是因?yàn)镻PP卡爾曼濾波法方程為對稱正定矩陣，在矩陣求逆時(shí)Cholesky分解法更為適用。