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

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

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

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

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

1 矩陣求逆

1.1 廣義逆

廣義逆是線性代數(shù)中針對矩陣的一種運算[15]。一個矩陣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，需要一個合適的m×n矩陣G，即

A·G·A=A

(3)

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

G·A·G=G

(4)

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

廣義逆矩陣理論可用于推導(dǎo)線性方程組的解。對于上述任何一種廣義逆陣都可以用來判斷線性方程組是否有解，若有解時列出其所有的解。若以下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)解存在的條件是當(dāng)且僅當(dāng)Ag·b為其中一個解，也就是當(dāng)且僅當(dāng)A·Ag·b=b。

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

1.2 SVD分解求逆

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

M·x=λ·x

(7)

式中：λ是矩陣M的一個特征值；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的奇異值分解。

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

M=V·Σ-1·U*

(9)

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

1.3 Cholesky分解求逆

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

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

A=L·L*

(10)

式中：L是一個下三角矩陣且所有對角線元素均為正實數(shù)；L*是矩陣L的共軛轉(zhuǎn)置矩陣。每一對稱正定矩陣都有一個唯一的Cholesky分解過程。Cholesky分解的計算時間復(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 精密單點定位

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

2.1 觀測方程

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

(15)

(16)

(17)

其中，

(18)

(19)

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

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

(20)

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

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

2.2 參數(shù)估計

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

其中，

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

3 實驗分析

國際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)品綜合與服務(wù)中心、1個運行管理控制中心和通信網(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ù)的采集，并進行數(shù)據(jù)合理性檢驗和數(shù)據(jù)預(yù)處理，并將數(shù)據(jù)發(fā)送至數(shù)據(jù)中心備份。中國電子科技集團公司第五十四研究所長期穩(wěn)定運行一個監(jiān)測評估中心，研發(fā)多頻多系統(tǒng)GNSS監(jiān)測接收機，并裝配16個iGMAS跟蹤站，iGMAS跟蹤站名稱及全球分布情況如表1所示。

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

為了分析不同矩陣求逆方法在GNSS精密單點定位參數(shù)估計中的應(yīng)用效果，采用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基礎(chǔ)上開發(fā)的GNSS Data Processor軟件。在動態(tài)PPP數(shù)據(jù)處理實驗中，將接收機坐標參數(shù)作為白噪聲進行估計。

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

接下來分析SVD分解求逆與Choleskey分解求逆的計算速度。分別采用以上這兩種方案處理BJF1站同一天的觀測數(shù)據(jù)。在方案1中，采用SVD分解法處理PPP卡爾曼濾波中法方程求逆問題。在方案2中，則使用Choleskey分解法，求解相應(yīng)的逆矩陣。統(tǒng)計兩種方案中每個歷元所需要的時間，包括讀取文件、數(shù)據(jù)預(yù)處理和矩陣運算，其中矩陣求逆是主要項。兩種方案的計算時間分別如圖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個歷元所需時間較長，通過處理日志判斷該歷元觀測值粗差較多，需要多次矩陣求逆運算，導(dǎo)致該歷元計算所需時間較長。

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

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

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

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

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

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

4 結(jié)語

矩陣分解求逆在工程應(yīng)用中具有重要價值。論述了矩陣求逆的基本理論和方法，介紹矩陣廣義逆的性質(zhì)及其在線性方程解中的作用。重點論述了兩種矩陣分解求逆方法，包括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分解法應(yīng)用于矩陣求逆的時效性更好，動態(tài)定位精度統(tǒng)計結(jié)果也間接證明了Cholesky分解法矩陣求逆精度更高。這是因為PPP卡爾曼濾波法方程為對稱正定矩陣，在矩陣求逆時Cholesky分解法更為適用。