高度約束MHE算法在偽距單點定位中的應用

孫淑光，楊祥源，陳萬通,，張巨聯，劉慶，任詩雨

( 1. 中國民航大學 電子信息與自動化學院, 天津 300300；2. 中國民航大學 民航航班廣域監視與安全管控技術重點實驗室, 天津300300；3. 上海飛機設計研究院, 上海 201109；4. 上海航天電子技術研究所, 上海 201109 )

0 引 言

在進行偽距單點定位(SPP)時，數據往往采用后處理的方式，傳統的精密衛星軌道及鐘差產品存在時間延遲的問題，限制了SPP技術在實時場景中的應用[1]. 近年來，國際全球衛星導航系統(GNSS)服務(IGS)推出了一些具有高精度的GNSS精密定位產品，可以很大程度上消除衛星定位相關誤差的影響[2].因此，一些相關研究重新聚焦于SPP技術.

SPP是僅利用一臺接收機即可實現高魯棒性的一種定位技術，而定位精度低仍是限制其實際應用的首要因素. 為進一步提高SPP的定位性能，近年來，研究人員發現一些場景下存在的先驗信息可以用來增強單點定位性能.

文獻[3]在SPP觀測模型中加入了先驗高度，將其作為準觀測方程，用于加權最小二乘(WLS)計算，實驗結果表明水平定位精度提高了42.84%. 文獻[4]提出了先驗海平面的高度約束精密單點定位(PPP)方法，并將約束方程通過差分處理實現線性化，實驗結果表明附加高度約束的定位精度在垂直方向上提高了79%，三維方向上提高了67%. 文獻[5]采用道路幾何模型約束進行導航定位，將車輛的狀態投影在地圖段上，實驗結果表明約束后的定位精度提高了31.4%.

傳統的SPP解算方法通常采用最小二乘(LS)估計[6]. 而LS只考慮觀測模型，不涉及系統狀態模型，故所得定位結果通常顯得粗糙雜亂. 因此，在一些研究中，研究人員引入了卡爾曼濾波(KF)算法，以實現KF算法的各種修改形式[7]. 其中，平滑約束卡爾曼濾波方法(SCKF)和無跡卡爾曼濾波方法(UKF)的出發點是先將非線性約束進行線性化，然后將其實現為完美測量[8-9]. 結果表明：上述方案所獲得的估計值僅近似滿足于非線性約束，而且非線性比簡單的線性化過程會具有更好的匹配度[9]. 不同于KF算法，滾動時域估計(MHE)是把帶約束的估計問題轉化成優化問題，同時把施加在系統變量上的約束條件直接表述為優化問題中目標函數的約束，并通過在線滾動優化進行求解[10]，是在約束條件下處理狀態參數估計的最有效方法[11-12]. 同時MHE算法可能會用到KF所不包括的信息量，而且對非線性約束可以直接應用，無需線性化.

基于以上考慮，本文將非線性高度約束方程引入到約束MHE濾波算法中進行SPP的解算，從而實現定位精度的提高，獲得更為平滑的定位結果.

1 SPP數學模型

SPP是以偽距觀測為基礎的一種定位方法，其偽距觀測是指對所觀測衛星和用戶接收機之間距離的測量.

1.1 消電離層觀測模型

本文分析的是雙頻的消電離層(IF)的偽距觀測，其觀測方程為[13]

式中：Pifj為基于載波頻率fj(j=1,2) 的用戶接收機到衛星i之間的偽距觀測； ρi為接收機和衛星的天線相位中心之間的真實幾何范圍； dtr和 dti分別為接收機和衛星的時鐘偏移；c為真空中光速；Ti為對流層延遲；bfj和Bfj分別為接收機和衛星的頻率相關偽距硬件延遲； εi為觀測噪聲和未建模的殘余誤差.

在上述觀測模型中，電離層延遲誤差通過IF組合模型的方式來消除. 為提高SPP的定位精度， dti、dtr、bfj和Bfj項應完全消除，在此，利用IGS精密產品消除衛星位置和時鐘誤差，利用對流層模型修正對流層延遲誤差. 而對于相對論效應、地球自轉效應、地球固體潮和海洋潮汐引起的誤差，若采用同樣的方式對其修正，反而致使得到的定位結果精度不高，故在此SPP模型中可被忽略，不作考慮.

此外，衛星鐘差項中集成了與衛星相關的硬件延遲，此由精確的衛星時鐘誤差乘積進行修正，同時與接收機相關的硬件延遲被整合到了接收機時鐘誤差項中，表示為[14]

式中： ζ 為衛星鐘差 dti與衛星硬件延遲Bpj的整合誤差； μ 為接收機鐘差 dtr與接收機的硬件延遲bpj的整合誤差. 當上述相關誤差和噪聲被忽略或校正后，式(1)可以寫為

1.2 SPP參數估計

與傳統的迭代算法不同，偽距觀測方程的閉式解方法不需要迭代運算. 通過計算參考衛星和其他衛星偽距方程之間的線性無關差，則可將線性化應用于式(3).

假設有n顆(n≥5 )按仰角降序排列得可見衛星，則可以獲得關于矩陣向量形式的未知參數(x,y,z,μ)的線性觀測公式，記為

式中：r=[x,y,z]T；em=[1,1,···,1]T，m=n-1 ； ν 為 觀測噪聲；l=[ε1/δ1,ε2/δ2,···,εm/δm]T

，且

式中， χi=[αi/δi,βi/δi,γi/δi] ，i=1,2,···,n-1 ，有[15]為后續應用滾動時域估計算法，采用正交變換的方法來消除式(4)中的接收機時鐘誤差項μ.假定P∈Rm×m是正交變換矩陣，則有0,···,0)T. 根據householder 變換矩陣的定義，并通過相應代數計算，矩陣P可以具體表示為

為消除接收機時鐘偏移參數 μ ，將式中的轉換矩陣P代入式(4)，可以獲得初始正交變換為

特別注意，由于只有第一個等式包含時鐘偏移項μ，故余下的部分可寫為

由此，模型的觀測方程可表述為

式中：y為觀測向量；H為向量r對應設計矩陣； ε 為誤差修正項.

因此，可以使用上述線性方程直接求解接收機位置參數解，該方程不包括接收機時鐘偏移項 μ .

為求解式(10)中的未知參數r，可以采用一般LS. 同時基于LS估計是以下最小化問題的解決方案，可表示為[16]

■■■■■■■

2 高度約束型MHE算法

2.1 高度約束

在假設先驗高度信息沒有誤差的情況下，具有非線性等式的高度約束方程可表示為

式中：D=diag[1/(a+h),1/(a+h),1/(b+h)] ，a、b分別為WGS-84參考橢球的長半軸和短半軸；h為相對于WGS-84參考橢球的高度[17].

2.2 滾動時域估計

MHE算法按所用信息量的不同可分為全信息MHE和近似MHE. 全信息MHE方法利用了當前和過去時刻的所有測量數據，測量結果更為準確，但另一方面會造成“滾雪球”效應，導致問題不可解；近似MHE方法則是為了避免上述問題的出現，通過限制問題的維數，確保了求解的可行性.

2.3 附加高度約束的滾動時域估計

在實際應用中，有很多的狀態方程都可以用來模擬接收機的狀態模型. 本文使用如下狀態方程來模擬接收機的狀態模型：

式中：T為連續兩個歷元的時間間隔；系統噪聲Q為：

式中：qx、qy、qz分別為不同分量上的過程噪聲的功率譜密度；為使觀測方程式(10)中的狀態向量rk與過程方程式(14)中的狀態向量Xk相同，將矩陣修改為包含有3個零列，如下所示：

因此，所得SPP觀測方程表示為

式中，R為觀測噪聲. 同理，將高度約束方程的相關矩陣D補充相應零列修改為

相應地，高度約束方程則可表示為

因此，高度約束MHE狀態估計問題的目標函數構造如下：

滿足約束

為限制優化問題的維數，引入了長度為J的固定時域(J=N-M+1 )，將全信息估計問題的計算劃分為兩部分：{1 ≤k≤M-1 } 和{M≤k≤N}，采用近似MHE算法將狀態估計問題轉為如下約束優化問題的目標函數

滿足約束

式中： {Xk} 為 {XM,...,XN} ；N-M+1 為 滾 動 時域窗口長度；其維度為為后驗估計協方差矩陣.

3 實驗結果與分析

為評估MHE算法在高度約束SPP的性能，下面分別進行了靜態實驗和動態實驗.

3.1 靜態實驗

本靜態實驗是基于AREG觀測站(站點位置：阿雷基帕，經度為-71.492 9°，緯度為-16.465 4°，海拔高度為2 489.336 m)根據IGS提供的24 h GNSS觀測數據，并以靜態接收機的真實高度作為先驗高度值，在MATLAB平臺上進行了SPP定位實驗. 實驗方案如下：

1)無約束LS的SPP實驗；

2)視窗長度為2 h，高度約束MHE的SPP實驗(MHE-2)；

3)視窗長度為4 h，高度約束MHE的SPP實驗(MHE-4).

表1分別計算出了不同方案下SPP定位估計過程中的處理時長. 時間跨度是由MATLAB中的tic和toc函數計算，其中tic功能自動記錄啟動時間，toc計算程序運行時的時間差，可見相對LS算法，MHE算法的處理時間較長，且隨MHE視窗大小的增加，程序運行時間也會急劇增加.

表1 不同方案的處理時長 s

圖1為SPP在WGS-84直角坐標系下的定位估計誤差. 利用接收機參考坐標與估計坐標之間的差值來計算誤差值. 考慮到MHE計算量較大，所以圖1分別給出了WGS-84直角坐標系下LS方法和前300個歷元的MHE方法的SPP的估計誤差結果. 可以看出，LS定位誤差波動起伏較大，結果粗糙雜亂；而附加高度約束的MHE方案抗擾動性能更強，定位結果更為平滑，明顯優于傳統無約束LS方案的定位結果.

圖1 三種方案在X、Y、Z方向上的定位誤差比較

結合表2，與LS方案相比，MHE-2方案在X、Y、Z三個方向分量上的均方根(RMS)值分別從1.11 m、3.23 m、1.01 m降低到0.21 m、0.91 m、0.23 m，整體定位精度提高了73.03%；同樣，MHE-4方案的RMS值分別降低到0.20 m、0.87 m、0.26 m，整體定位精度提高了81.46%.

表2 SPP定位誤差RMS值

3.2 動態實驗

為了進一步評估附加高度約束MHE的SPP定位性能，通過實時采集的動態GPS衛星觀測數據進行了動態SPP定位實驗.

如圖2所示本實驗數據于2021年6月2日采集于武漢的中國科學院測量與地球物理研究所導航樓頂，從GPS時間(CPST)的10:37:00開始，觀測時長約5 min. 實驗過程中的動態接收機沿著自身圓形軌道勻速轉動. 此外，在距離該數據采集地點約4.5 m的地方設有一個基準站(精確坐標已知)同時進行觀測獲取衛星觀測數據，然后用RTKLIB軟件的載波相位差模式進行解算動態接收機的位置坐標，以獲得接收機的精確定位數據. 圖3為實驗中動態接收機的運動軌跡坐標.

圖2 動態數據收集站點

圖3 動態接收機運動軌跡

本次動態實驗中，先驗高度值被設置為所有歷元的接收機定位的平均高度. 為了與靜態實驗相一致，并且較小視窗尺寸的MHE方法的定位精度的提高對實際應用也更有意義，因此，實驗中MHE的視窗尺寸也設置為2和4.

圖4為動態接收機做勻速圓周運動時，LS方法和前300個歷元的MHE方法的SPP定位誤差結果. 結果表明：與無約束LS相比，MHE兩種方案的結果具有較小的波動，且MHE-4方案的結果優于MHE-2方案. 同時也不難發現，隨時間推移，動態接收站工作趨于穩定及MHE-4方案觀測數據信息的增加，其定位結果在250 s后相較與LS和MHE-2結果存在明顯偏離，定位結果更為精確.

圖4 三種方案在X、Y、Z方向上的定位誤差比較

表3給出了前300個歷元SPP定位的RMS誤差. 由表3可知，與LS方案相比，MHE-2方案在X、Y、Z三個方向分量上的RMS值分別從3.81 m、7.31 m、4.14 m降低到2.53 m、4.70 m、3.11 m，整體定位精度提高了33.0%； MHE-4方案的RMS值分別降低到1.95 m、3.23 m、2.26 m，整體定位精度提高了52.4%.與MHE-2方案相比，MHE-4方案的整體定位精度提高了19.4%.

表3 SPP定位誤差RMS值

4 結 論

SPP技術是一項重要的定位技術，可為一些特定的應用場景提供高魯棒性的定位. 本文從國內外研究的先驗信息可以用來增強單點定位性能的基礎上，提出了附加高度約束的MHE算法來進一步提高SPP定位精度，并進行了靜態和動態實驗分析. 結果表明：相對于傳統的LS方法，附加高度約束的MHE算法的單點定位結果更加準確和平滑，驗證了該方法的有效性、可行性，可適用于更為廣泛的應用場景.

但同時MHE算法的估計精度與運行效率互相矛盾，兩者不可兼得，其中運行效率主要與視窗的大小有關. 在實際應用過程中，需要根據應用場景的不同，在定位估計精度和計算工作量之間折衷考慮視窗參數大小.