基于格理論的模糊度快速解算方法

王守華 ，吳黎榮 *，紀(jì)元法 ，孫希延

（1. 廣西精密導(dǎo)航技術(shù)與應(yīng)用重點實驗室（桂林電子科技大學(xué)），廣西桂林541004；2. 衛(wèi)星導(dǎo)航定位與位置服務(wù)國家地方聯(lián)合工程研究中心（桂林電子科技大學(xué)），廣西桂林541004）

0 引言

隨著全球衛(wèi)星導(dǎo)航系統(tǒng)（Global Navigation Satellite System，GNSS）的融合和發(fā)展，各行各業(yè)對導(dǎo)航定位的精度、可靠性的需求越來越高，而高精度GNSS 定位通常是基于載波相位觀測值進(jìn)行的，這使得快速、準(zhǔn)確地解算模糊度整型值成為了該領(lǐng)域的核心研究問題［1］。由于整周模糊度參數(shù)的離散性，通常都是基于模糊度的實數(shù)解及協(xié)方差矩陣構(gòu)造出一定的搜索空間，再進(jìn)一步采用搜索方法來獲取模糊度的最優(yōu)解，而現(xiàn)有的技術(shù)方法大多在高維情況下對此解算緩慢且成功率低。國內(nèi)外學(xué)者均對此問題做了深入研究，先后提出了不同的模糊度解算方法，如模糊度函數(shù)法、最小二乘搜索法以及快速搜索法，這些方法都是直接采用搜索的技術(shù)獲得模糊度的整型值。文獻(xiàn)［2］中提出的最小二乘模糊度降相關(guān)（Least squares AMBiguity Decorrelation Adjustment，LAMBDA）算法的基本思想是通過整數(shù)變換降低模糊度的相關(guān)性，然后在整數(shù)變換的空間內(nèi)采用搜索技術(shù)獲得最優(yōu)模糊度整型值；文獻(xiàn)［3］中根據(jù)改進(jìn)的球形搜索算法進(jìn)行擴(kuò)展，提出了改進(jìn)的LAMBDA 算法（MLAMBDA），縮小了模糊度整型搜索空間，提高了模糊度的搜索效率；文獻(xiàn)［4］中將格理論中的格基規(guī)約算法應(yīng)用到模糊度實數(shù)解的降相關(guān)中，但由于當(dāng)時各方面條件的欠缺，并沒有通過實測數(shù)據(jù)來對算法理論進(jìn)行分析。

本文基于格理論對GNSS 模糊度解算算法進(jìn)行改進(jìn)，并提出了最近格點（Closest Lattice Point，CLP）搜索算法對模糊度整型值進(jìn)行搜索。該算法首先將模糊度搜索轉(zhuǎn)化為對格中已知格點的最近格點搜索問題，再根據(jù)格基規(guī)約改進(jìn)得出具有最小可能長度且相互正交的格基向量，最后采用CLP 搜索算法搜索出最優(yōu)的模糊度參數(shù)值。

1 GNSS相對定位數(shù)學(xué)模型

對于GNSS 高精度定位雙差模式下的觀測方程，可以用下面通用的GNSS系統(tǒng)模型［5］來概括：

式中：L表示GNSS 雙差觀測值；b為觀測模型中的實數(shù)未知參數(shù)，包括位置、對流層和電離層延遲參數(shù)；B 為其相應(yīng)的系數(shù)矩陣；a 為模糊度未知參數(shù)；A 為其相應(yīng)的系數(shù)矩陣；e 為觀測噪聲。

用最小二乘法對式（1）進(jìn)行模糊度參數(shù)最優(yōu)值估計，得：

式中：P 是觀測值的權(quán)重矩陣，通過忽略模糊度參數(shù)a 的整數(shù)特性直接求出方程中參數(shù)浮點解a^ 和b^，及其對應(yīng)的協(xié)方差矩陣Qa^和Qb^。其中式（2）中的最小化問題可轉(zhuǎn)化成如下問題：

由于參數(shù)b 為實數(shù)，所以式（3）可以等價于右邊兩項最優(yōu)估計，因此有

求解載波相位的模糊度最優(yōu)值-a，最為經(jīng)典的是基于最小二乘搜索的LAMBDA算法和在此基礎(chǔ)上改進(jìn)的MLAMBDA算法。

2 模糊度搜索與最近格點搜索的等價關(guān)系

由格的定義可知，設(shè)向量I 為一組格基向量且線性無關(guān)，則格的空間形式［7-9］可以表示如下：

式中，Λ(·)為格的數(shù)學(xué)表達(dá)式，如果格的所有格點α 相同，則格中的基本單元的容積大小不會因格基I 的變換而變化。根據(jù)這一思想可以通過格基規(guī)約對GNSS 模糊度的搜索區(qū)域進(jìn)行搜索，但前提是構(gòu)建與模糊度參數(shù)相關(guān)的格。

由格的一個著名問題知，對于一個已知的格中的點y，該點也許不是格點，那它肯定包含在由格基構(gòu)成的實線性空間中，并在格中存在一個格點x滿足：

其中格點x 被稱為y 在格中的最近格點。對于求解這個最近格點的問題被稱為最近格點搜索問題，格可看作是呈網(wǎng)格狀以一個點一個點的形式存在于n 維空間，在整數(shù)單位二維情況下示意圖如圖1所示。

圖1 給定點y的最近格點x搜索示意圖Fig.1 Search diagram of the closest lattice point x for the given point y

由式（3）、（4）知，估計的模糊度參數(shù)最優(yōu)值的協(xié)方差矩陣Qa^是正定的，因此權(quán)重矩陣P的Cholesky分解可以為：

式中R為下三角矩陣，同時式（4）的右式可以轉(zhuǎn)化為：

由式（9）、（10）可知式（4）最終可變換成：

同時由式（7）的最小格點搜索問題知，對于給定的格中的點 y ∈ Rn，搜索出最接近 y 的格點 x，存在 α ∈ Zk，使得對于x ∈ Λ，x =|Iα|，有如下等式：

由式（11）和（12）可以得出，由最小二乘法搜索模糊度參數(shù)最優(yōu)值等價于格上的最近格點搜索問題，因此可以將模糊度搜索轉(zhuǎn)化為對格中已知格點的最近格點搜索問題。

3 改進(jìn)的模糊度降相關(guān)處理

基于最小二乘法對模糊度參數(shù)解算時，為了消除參數(shù)的相關(guān)性會進(jìn)行降相關(guān)處理。而對格中已知點的最近格點搜索時，往往采用格基規(guī)約，通過將給定的格基變換為另一個所需格基，且保持格的大小不變，搜索出具有最小可能長度且相互正交的格基向量。因此模糊度解算中的降相關(guān)處理等價于格基規(guī)約，現(xiàn)對格基規(guī)約的三個基本步驟：正交化、尺寸減小和格基變換進(jìn)行改進(jìn)。

由式（11）中矩陣P 的Cholesky 分解得出的下三角矩陣R進(jìn)行Gram-Schmidt正交化得到：

則由式（8）和（13）可得：

為了獲得更好的尺寸減小，作出以下改進(jìn)：格基中ii的順序通過與對應(yīng)的Gram-Schmidt 正交基元素i*i 變換來改變，為了得到對應(yīng)的正交基i*i 更短的值，可以通過對ii和ii+1進(jìn)行變換，連續(xù)向量ii和ii+1的格基變化條件如下：

通過格基規(guī)約，對矩陣R進(jìn)行規(guī)約，獲得接近正交和最小長度的格基，在GNSS 定位解算中，基于矩陣R 獲得的最小長度且正交的格基即為對應(yīng)的模糊度參數(shù)值，由這一步提高了模糊度搜索的效率和成功率。

4 基于格理論的模糊度搜索算法

由第3章中已知點的最近格點搜索問題，令r0為以y為中心的n維粒子的平方半徑，得：

同時由于R 為下三角形矩陣，不等式有下面一組條件關(guān)系：

式中：Rj，k為矩陣R的第j行第k列元素；ak為對應(yīng)的第k個需要固定的整數(shù)值，由于模糊度參數(shù)a值的整數(shù)特性，如果對于i +1 ≤j ≤n，aj表示已經(jīng)固定的參數(shù)值，則變量ai（i = n - 1，n -2，…，1）可以在整數(shù)范圍[Li，Ui]取整數(shù)值，區(qū)間范圍如下：

式中：round(*)表示取最接近的整數(shù)值；yi為對應(yīng)的第i個觀測中心點；ak為第k個固定的整數(shù)值；Ri，j為矩陣R 中的第i行第j列元素。

而對于已經(jīng)固定的aj(j = i + 1，…，n)的點，如果它們不屬于以y為中心的半徑為 r0的球體中，那這些點與y之間的歐幾里得距離平方可表示為：

通過以上條件，本文提出了CLP搜索算法，并對整數(shù)范圍[Li，Ui]搜索順序分級別處理，其中對每個級別在區(qū)間的中點開始，以Z 字形順序搜索，同時在搜索的過程中根據(jù)式（21）的限制條件估計確認(rèn)是否繼續(xù)搜索。先把r0設(shè)置為無窮大，設(shè)區(qū)間中點為Mi可得如下等式：

算法搜索的具體步驟可總結(jié)如下：

2）搜索過程從搜索第一個ai（i = n）開始，設(shè)此級別為第1級，生成的第一個格點設(shè)為aB，得出aB后，將r0設(shè)置為距離d2(y，Ra)。

3）再按照步驟1）開始搜索第2 級，對ai（i = n - 1）按照Z字形要求搜索得到整數(shù)點對應(yīng)的最小搜索半徑ri；若ri＜r0，則更新搜索半徑，使r0= ri，再進(jìn)行下一步。

4）依次對第3，4，…，n 級進(jìn)行搜索，按Z 字形順序搜索并得出整數(shù)點對應(yīng)的最小半徑，再與上一級得到的搜索半徑進(jìn)行比較，若小于，則對其進(jìn)行更新，否則同理進(jìn)行下一級搜索比較，直至未能在第n 級找到新的整數(shù)值時，搜索過程停止，則最小搜索半徑的整數(shù)點對應(yīng)的為最優(yōu)整數(shù)解。

5 算例及仿真分析

5.1 數(shù)值模擬分析

為了驗證本文提出的CLP搜索算法對模糊度搜索的有效性，基于Matlab 平臺在不同案例情況對模糊度參數(shù)最優(yōu)值搜索時間進(jìn)行比較驗證，針對不同案例，真實的模糊度參數(shù)浮點解a^ 被構(gòu)造為：

其中randn(n，1)是Matlab內(nèi)置函數(shù)，用于生成具有標(biāo)準(zhǔn)正態(tài)分布的n個目標(biāo)向量。這些向量都根據(jù)真實模糊度參數(shù)特性構(gòu)建而成，并構(gòu)成如下7 種不同模擬案例，前四種基于模糊度參數(shù)協(xié)方差矩陣Qa^=LTDL，其中L是單位下三角矩陣，每個lij（對于i＞j）都是由randn生成的隨機(jī)數(shù)。案例特點具體如下：

案例1：D= diag(di)，di= rand，其中rand是Matlab 內(nèi)置函數(shù)，用于在（0，1）中生成均勻分布的隨機(jī)數(shù)。

案例2：D= diag(n-1，(n- 1)-1，…，1-1)。

案例3：D= diag(1-1，2-1，…，n-1)。

案例4：D= diag(2 00，200，200，0.1，…，0.1 )。

案例5：協(xié)方差矩陣Qa^=UDUT，U是通過模擬的模糊度參數(shù)隨機(jī)矩陣QR分解獲得的隨機(jī)正交矩陣，D= diag(di)，di=rand。

案例 6：Qa^=UDUT，U由與案例 5 相同的方式生成，其中U中的其他對角元素隨機(jī)分布在di和dn之間，n是Qa^的維數(shù)。因此Qa^的條件數(shù)是

案例7：Qa^=ATA，A =randn(n，n)。

其中，每個案例都是有模擬理論的，如案例4，因為GNSS中的協(xié)方差矩陣Qa^通常在第三條件標(biāo)準(zhǔn)偏差和第四條件標(biāo)準(zhǔn)差之間存在較大差距，所以模擬下把其元素前后設(shè)置差距較大。在所有的模擬案例實驗中，其他條件設(shè)置一致，且對參數(shù)的搜索成功率被默認(rèn)為是100%，所有呈現(xiàn)的結(jié)果都在50次獨立運(yùn)行中執(zhí)行。

表1 基于統(tǒng)計學(xué)原理統(tǒng)計出維度為40 時，CLP、LAMBDA和MLAMDA 搜索算法在模擬案例1～7 下50 次獨立運(yùn)行的平均搜索時間，而且每次運(yùn)行中都統(tǒng)計到估計出參數(shù)最優(yōu)值的搜索耗時。

表1 三種算法在模擬案例1～7上50次獨立運(yùn)行的平均搜索時間Tab.1 Average search time of three algorithms in 50 independent runs on simulation cases 1～7

由表1可知，模糊度維數(shù)為40時，CLP 搜索算法在模擬案例1～7 下50 次獨立運(yùn)行的平均搜索時間約為LAMBDA 的1/97、1/374、1/349、1/179、1/64、1/205、1/187，約為 MLAMBDA的1/6、1/5、1/6、1/6、1/7、1/7；且MLABDA 的平均搜索時間約為LAMBDA 的 1/16、1/78、1/56、1/23、1/11、1/27、1/26 倍。根據(jù)數(shù)值模擬結(jié)果可以得出對模糊度浮點解參數(shù)a^ 的最優(yōu)值搜索，CLP搜索算法比LAMBDA和MLAMBDA的耗時更短。

圖2 和圖3 分別展示出CLP 搜索算法在案例1～7 進(jìn)行50次獨立運(yùn)行的平均搜索時間和候選者數(shù)量，可見CLP 搜索算法搜索的平均時間和搜索的候選者平均數(shù)量是一一對應(yīng)的關(guān)系。其中當(dāng)維數(shù)為10～20 時，各種案例的搜索時間基本無差別，但隨著維數(shù)的增加，其他案例的搜索時間都呈指數(shù)級增長，只有案例3 的搜索時間在相對穩(wěn)定的范圍內(nèi)波動，且案例3 正好與CLP 搜索算法解算參數(shù)設(shè)置相一致，因此可以得出CLP 搜索算法即使在高維情況下，對于模糊度的搜索依舊高效穩(wěn)定。

5.2 實測數(shù)據(jù)實驗仿真

為了進(jìn)一步驗證CLP 搜索算法在實際場景下的效果，采用了兩組不同時間段同一基線的GNSS 觀測數(shù)據(jù)，使用Bootstrapping成功率作為模糊度解算的可靠性指標(biāo)，按照上述相對定位數(shù)學(xué)模式對20 顆衛(wèi)星進(jìn)行單歷元測試。其中實測實驗條件為使用ublox M8T 型低成本接收機(jī)在某大學(xué)校園采集GNSS 觀測數(shù)據(jù)，數(shù)據(jù)包含單頻全球定位系統(tǒng)（Global Positioning System，GPS）和北斗導(dǎo)航系統(tǒng)（BeiDou navigation Satellite system，BDS）的載波相位、偽距觀測量，如圖4為一天中可觀測的衛(wèi)星數(shù)量（圖中最上方曲線為GPS+BDS 衛(wèi)星數(shù)）。其中，采集數(shù)據(jù)的兩個觀測站點固定不動，且觀測環(huán)境相對較好，兩站之間的距離約為100 m，位置坐標(biāo)均采用靜態(tài)后處理軟件進(jìn)行長時間解算得到。CLP 降相關(guān)處理時設(shè)置δ= 1，使得各算法基于Qa^=LTDL分解的降相關(guān)性能一致。數(shù)據(jù)采集時間段為 2019 年 3 月 18 日的 11 時 00 分到 13 時 00 分，共計7200歷元，采樣歷元間隔為1 s。

在實際長時間解算中，算法存在抗干擾性差和不穩(wěn)定性等問題。為了分析比較算法在對模糊度解算速度性能上的優(yōu)勢，選取其中兩組不同實測數(shù)據(jù)實驗得出三種算法的模糊度最優(yōu)值的搜索時間和固定時間，結(jié)果如圖5和圖6所示。從圖中可以看出，CLP 搜索算法的搜索時間穩(wěn)定在0.01 s，LAMBDA 和 MLAMBDA 則分別為 0.052 s 和 0.03 s，且不穩(wěn)定；當(dāng)對模糊度固定解算時，增加去相關(guān)處理，固定時間相比搜索時間略有增加，搜索次數(shù)隨著歷元數(shù)的增加而迅速減小而后平穩(wěn)，且都能明顯看出CLP 搜索算法比LAMBDA 和MLAMBDA算法解算得更快。

圖4 可觀測衛(wèi)星數(shù)Fig.4 Number of observable satellites

圖5 三種算法在兩組不同數(shù)據(jù)上的模糊度最優(yōu)值搜索時間Fig.5 Search time for the optimal ambiguity by three algorithms on two different sets of data

圖6 三種算法在兩組不同數(shù)據(jù)上模糊度最優(yōu)值固定時間Fig.6 Fixed time for the optimal ambiguity by three algorithms on two different sets of data

6 結(jié)語

本文針對未來各GNSS 系統(tǒng)兼容與互操作，多頻多模高維模糊度在常規(guī)方法下解算效率低的問題，基于格理論，對GNSS模糊度解算算法進(jìn)行改進(jìn)，并提出了CLP搜索算法對模糊度整型值進(jìn)行搜索。該方法首先將模糊度搜索轉(zhuǎn)化為對格中已知格點的最近格點搜索問題，再根據(jù)格基規(guī)約改進(jìn)得出具有最小可能長度且相互正交的格基向量，最后采用CLP 搜索算法搜索出最優(yōu)的模糊度參數(shù)值。通過模擬實驗和實測數(shù)據(jù)實驗驗證得出，所提的CLP 搜索算法理論上相較經(jīng)典的LAMBDA/MLAMBDA 算法對最優(yōu)值模糊度參數(shù)解算效率更高、更可靠，且每一個參數(shù)搜索時間穩(wěn)定在0.01 s，即使在高維情況下，CLP搜索算法的搜索依然穩(wěn)定可靠。

基于格理論的模糊度解算方法的提出打破了常規(guī)解算思想，給未來多頻多模系統(tǒng)的高精度解算提供了新的思路。本文主要針對GPS 和BDS 雙系統(tǒng)情況進(jìn)行解算，而對于多頻多模GNSS 系統(tǒng)將進(jìn)一步研究；同時在長基線情況下由于雙差模型難以消除存在強(qiáng)相關(guān)性的參數(shù)之間的影響，致使分離模糊度參數(shù)耗費(fèi)時間長且不可靠，解算結(jié)果是否依然有效可靠還需進(jìn)一步研究。