M-Cholesky分解法快速求解GNSS整周模糊度

李克昭 田晨冬

1 河南理工大學測繪與國土信息工程學院, 河南省焦作市世紀路2001號,454003 2 北斗導航應用技術協同創新中心,鄭州市科學大道62號,450001

目前應用最廣泛的模糊度固定方法是基于整數最小二乘原則在模糊度域中確定模糊度的最小二乘模糊度去相關算法LAMBDA[1]。LAMBDA算法主要包括3個步驟:模糊度降相關、模糊度搜索、模糊度檢驗。其中,模糊度降相關涉及大量矩陣運算,耗時較長,因此有必要對其降相關過程進行研究,提高計算效率。國內外學者對降相關過程進行了大量研究:Liu等[2]基于上下三角過程構造聯合去相關算法;Xu[3]使用Cholesky分解法構建整數高斯矩陣,在高維情況下分解得到條件數最小的浮點解協方差矩陣;陳樹新[4]提出一種以降低協方差矩陣條件數為準則的模糊度去相關算法;Liu等[5]克服了Z變換可能帶來的矩陣病態分解,減少了20%的解算時間;安向東等[6]根據頻間偏差率范圍,采用搜索方法和線性模型去除相位頻間偏差對寬窄巷模糊度的影響;邱蕾等[7]提出在位移值較大的情況下通過多種載波相位組合方法解算短基線GPS整周模糊度。自從Hassibi等[8]將格理論中的規約算法Lenstra-Lenstra-Lovász(LLL)應用于模糊度降相關以來,諸多學者從格基規約理論出發分析LLL的降相關性能:Xu[9]提出一種性能優于整數高斯去相關算法和LLL算法的逆整數Cholesky去相關方法;范龍[10]基于分塊Gram-Schmidt正交化算法改進LLL算法;謝愷等[11]基于系統旋轉的Houscholder正交變換對LLL算法進行改進;盧立果[12]采用分塊降維方法改進LLL算法。

上述方法中,Cholesky去相關法效果較好,但也存在一定的缺陷:1)L、D、LT子矩陣的形成過程過于復雜;2)在實際編程計算中,L、D、LT矩陣元素需要較多的存儲矩陣;3)各矩陣元素只能按順序求出;4)未能有效利用各子矩陣元素之間的關系。基于此,本文提出一種改進的M-Cholesky分解法,用合成矩陣取代L、D、LT矩陣,簡化了L、D、LT子矩陣的關系;之后,利用高斯算法中的四角規則計算各子矩陣元素,提升分解時的解算效率,可減少約40%的計算元素個數。

1 基于Cholesky分解的模糊度去相關模型

GNSS觀測方程的模型求解可以等價于求解下式的最小值問題:

(1)

由于模糊度分量之間通常具有很高的相關性,因此搜索效率較低。為了使求解過程更高效,可以通過整數高斯變換(Z變換)將其轉化為新的最小二乘問題:

(2)

(3)

式中,Z為幺模矩陣,即矩陣的行列式為1。矩陣Z的元素均為整數。

(4)

(5)

式中,I為n維單位矩陣,η為對L矩陣的第i行第j列的元素取整,ei和ej為單位坐標向量。

2 M-Cholesky分解算法

傳統Cholesky分解模型復雜,且需要4個數組分別存放Q、LT、D、L矩陣。本文通過研究矩陣間各元素的相應關系,引入合成矩陣來減少計算所需的數組,并結合高斯算法中的四角規則提出一種改進的M-Cholesky分解法。

以4階矩陣Q為例,給出其Cholesky分解:

(6)

將LT、D、L矩陣整合得到的合成矩陣與Q矩陣各元素進行比較,合成矩陣與Q矩陣的對應元素關系如式(7)中的合成矩陣所示:

(7)

可以看出: 1)利用式(7)進行Cholesky分解時,可在Q矩陣的基礎上直接計算得到合成矩陣,僅需要1個數組即可存放原本需要3個數組存儲的內存;2)合成矩陣中的元素lij與lji關于主對角線元素對稱相等,因此在進行Cholesky分解時,可以僅計算對角線元素以上的lij元素,lji=lij。

在上述分析的基礎上,利用四角規則給出M-Cholesky算法中合成矩陣元素的計算步驟:

1)基于dii對第i行元素進行規約化處理,得到lij元素;

3)將第i行的lij賦值給第i列的lji元素;

4)對(i-1),(i-2),…,1行依次循環,完成所有計算。

其中,步驟2)中參與計算的元素為:

(8)

由式(8)可以看出,每次參與合成矩陣計算的4個元素均位于矩形的4個角上,因此稱之為四角規則。在M-Cholesky算法計算過程中,所有計算元素轉換為相應的對角元素、交叉元素和消元元素時,最多只用到3個元素,因此可以將該計算過程視為規約化后的消元計算。M-Cholesky算法與Cholesky算法在計算過程中所需的計算元素個數如表1所示。

表1 M-Cholesky算法與Cholesky算法所需計算元素個數對比Tab.1 Comparison of numbers of computing elements required by M-Cholesky algorithm and Cholesky algorithm

由表1可見,M-Cholesky算法可減少約40%的計算元素,提高了計算效率,且不影響計算結果的精度。

3 實驗驗證

解算精度與時間效率是判斷一個模糊度解算算法是否可行的關鍵,但由于本文所提出的M-Cholesky算法僅改變了計算元素,并不影響定位精度,因此僅采用仿真數據與實測數據對Cholesky算法與M-Cholesky算法的運行時間進行比較。

本文基于Microsoft Visual Studio 2019,利用C語言進行M-Cholesky算法的程序設計,并進行數值模擬測試。所有計算均在Windows 10操作系統的電腦平臺上進行,處理器為16 GB內存的Intel(R) Core(TM) i7-8750H CPU@ 2.20 Hz。

3.1 仿真實驗

參照文獻[13]的仿真參數進行仿真實驗,仿真維數為5～40。為避免偶然情況的出現,每次仿真實驗取100次均值。2種算法的平均運行時間如圖1所示,所有情況下40維的平均運行時間如表2所示。

圖1 平均運行時間Fig.1 Average running time

表2 40維解算時間對比Tab.2 Comparison of 40 dimensional solution time

由圖1可見,在觀測維數為5～40的情況下,2種方法的計算時間雖然偶爾有波動,但整體上都隨維數的增加而增加,其中M-Cholesky算法計算時間總體低于Cholesky算法,且計算效率也隨維數的增加而提升。由表2可見,M-Cholesky算法的計算時間更短,相比于Cholesky算法,情況1～7中,M-Cholesky算法計算時間分別提升10.6%、13.5%、12.5%、16%、15.7%、13.5%、20.3%。

3.2 實測實驗

實測實驗分為3組:1)短基線組采用和芯星通UB4B0-MINI板卡采集的靜態觀測數據,基線長89.83 m,采樣間隔1 s,截取3 000個歷元對解算時間和結果進行分析;2)中長基線組采用武漢大學IGS數據中心下載的香港地區HKSL和HKWS站的IGS數據,2站構成的基線長度為42.51 km,采樣間隔30 s;3)長基線組采用武漢大學IGS數據中心下載的香港地區HKSL和武漢地區JFNG站的IGS數據,2站構成的基線長度為903.26 km,采樣間隔30 s。所使用的數據均為標準RINEX格式數據,包括星歷數據和觀測數據。

采用GPS、BDS、Galileo、GLONASS多系統數據進行整周模糊度融合解算,3種基線解算時間對比情況如圖2～4所示。

圖2 短基線解算時間對比Fig.2 Comparison of short baseline solution time

由圖2可見,Cholesky算法的解算時間大致在(1.6～1.8)×104μs之間,而本文提出的M-Cholesky算法解算時間大致在(1.5～1.7)×104μs之間,所需解算時間較少,且M-Cholesky算法的整周模糊度解算結果相對穩定。由圖3、4也可以看出,M-Cholesky算法在解算時間和解算穩定性上表現更好。需要說明的是,短基線實驗結果是從觀測到的17 570個歷元中截取的3 000個歷元,采樣間隔為1 s,而中長、長基線的觀測時段長于短基線,采樣間隔均為30 s,因此圖2解算時間曲線較圖3、4更平緩。

圖3 中長基線解算時間對比Fig.3 Comparison of medium baseline solution time

圖4 長基線解算時間對比Fig.4 Comparison of long baseline solution time

4 結 語

本文分析影響Cholesky算法計算效率的因素,引入四角規則,提出一種改進的M-Cholesky算法。實驗結果表明,相對于Cholesky算法,M-Cholesky算法求解模糊度的計算時間減少約15%,提高了整周模糊度固定的效率,且解算精度與Cholesky算法一致。

如何處理四角規則中分母為0或接近0時的特殊情況,是完善M-Cholesky算法需要進一步考慮的問題。