附加基線矢量約束的雙站協同精密定位方法

陳萬通, 尚正輝, 汪竹青

(中國民航大學 天津市智能信號與圖像處理重點實驗室,天津 300300)

0 引 言

近年來全球衛星導航系統(GNSS)精密單點定位(PPP)技術越來越受到GNSS領域的廣泛關注. PPP利用測碼偽距和載波相位觀測值以及國際GNSS服務(IGS)提供的高精度的衛星星歷和衛星鐘差產品,同時應用精化的誤差模型和誤差改正產品改正定位過程中的多種誤差,進而實現國際地球參考框架下的高精度單點定位[1]. PPP相比于傳統的單點定位和實時動態(RTK)定位方式,具有定位精度高,不受相對定位作用范圍限制,操作方便,成本低廉等優點. 由于上述優點,PPP具有廣泛的工程應用價值,現已全面應用到工程測繪、交通、大氣水汽反演、精密授時等多個領域.

PPP通常采用Kalman濾波進行位置及其他參數的解算. Kalman濾波進行參數解算時不僅要包含必要的觀測模型,同時引入了動態模型. 因此,Kalman濾波的精度不僅取決于觀測值的精度,還依賴于建立的動態模型和給定的初值是否準確. 然而,由于接收機運動狀態具有很多的隨機不確定性,很難通過精確建模來描述,初始信息也很難精確獲取,這給Kalman濾波遞推計算帶來了一定困難,不僅可能降低精度,濾波甚至可能會發散. 改善Kalman濾波的性能通常從建立更準確的動態模型和采用特殊的外部先驗信息入手. 許多學者提出了抗差濾波[2]、自適應選權濾波[3]等改進模型來控制觀測異常和動態模型誤差的影響,另外,國內外學者對于帶約束的濾波算法進行了一定的研究,文獻[4]分析了線性系統下三種約束濾波算法,并利用車輛跟蹤仿真實驗比較了三者的性能;文獻[5]將制約車輛運動的道路信息引入Kalman濾波模型,提高了濾波結果的精度. 但是,上述研究都主要集中于抗差自適應算法以及對殘差的研究,即使考慮外部約束信息也多是在目標跟蹤領域基于計算機進行二維仿真,將考慮附加約束信息的Kalman濾波與精密定位相結合并以實測衛星觀測數據進行驗證的研究較少. 鑒于此,本文討論了一種利用先驗基線約束信息對Kalman濾波進行改進的雙站協同PPP定位方法,該方法在濾波解算過程中充分考慮了狀態參數之間的函數關系,增強了模型的強度,不同程度上提高了定位解算的精度.

1 PPP數學模型

1.1 觀測模型

傳統PPP算法采用基于雙頻偽距和載波相位消電離層(IF)組合的觀測方程[6]:

(1)

(2)

對于上述模型,電離層延遲項已經采用觀測值IF組合的方式加以消除. 對流層延遲可分為天頂對流層干延遲和天頂對流層濕延遲,其中干延遲部分采用經驗模型來改正,濕延遲部分通常作為參數來估計. 通常,與衛星相關的硬件延遲項會合并到衛星鐘差項中,可通過IGS提供的精密鐘差進行改正;與接收機相關的硬件延遲項會合并到接收機鐘差中,作為一個整體進行參數估計. 相位觀測方程中的相位延遲項(包含硬件相位延遲項和初始相位偏差項)會被相應的模糊度參數吸收,使得模糊度參數失去整數特性[7]. 相位纏繞、固體潮、大洋潮汐、相對論效應、地球自轉等引起的誤差根據相應的模型進行修正. 最終式(1)和式(2)經過誤差改正后通常簡寫為:

(3)

(4)

(5)

(6)

1.2 Kalman濾波參數估計

PPP中通常采用Kalman濾波算法進行參數估計,若待估的狀態序列為{xk,k=0,1,…,n},觀測值為{zk,k=0,1,…,n},則Kalman濾波模型可由如下線性化后的動態模型和觀測模型表示[8]

xk=Fk-1xk-1+wk-1,wk-1～N(0,Qk-1),

(7)

zk=Hkxk+vk,vk～N(0,Rk-1).

(8)

式中:F為狀態轉移矩陣；H為設計矩陣；w和Q為系統噪聲及其協方差矩陣,v和R為觀測噪聲及其協方差矩陣. 在靜態定位中,狀態向量主要包含接收機的位置(X,Y,Z)、接收機鐘差dt、對流層延遲trop以及IF組合模糊度(Ni,i=1…，n)四類基本參數. 即:

(9)

狀態轉移矩陣F在靜態定位中可設為單位陣.

設計矩陣H可表示為

(10)

式中:α、β、γ為星站矢量的方向余弦;c為真空中的光速;Mwet為對流層延遲投影函數.其余為模糊度項系數.

觀測噪聲的協方差矩陣Rk可根據衛星高度角模型對偽距和相位觀測值分別賦值,即:

(11)

過程噪聲的協方差矩陣Qk為

(12)

式中,σ2為噪聲對應的方差,可根據經驗賦值.

對于靜態PPP,位置參數和模糊度在過程噪聲矩陣中對應的元素為0,接收機鐘差對應的元素取105,對流層濕延遲參數對應的元素取10-9[9].

2 附加基線約束的雙站協同PPP模型

2.1 附加基線約束信息的Kalman濾波

在實際的參數解算過程中,待估的狀態參數間會存在一些已知的函數關系,這些先驗的函數關系可以作為約束用來改善Kalman濾波的穩定性,從而增強PPP的定位性能[10].狀態約束信息通常可表示為一系列的等式或者不等式,對于這些約束信息通常有以下幾種處理方式:

1)模型簡化法,通過引入等式狀態約束來消除一些系統模型參數,從而簡化模型;

2)偽觀測法,將等式約束視為零觀測噪聲的高權重的完美觀測值,使得觀測方程得到增廣,加強觀測空間的結構[11],從而將約束估計問題轉化為非約束估計問題. 批處理方法和序貫處理方法是偽觀測法中常見的計算手段;

3)投影法,將未約束的Kalman濾波的狀態估計投影到約束空間[12]. 事實上,上述三種約束處理方法具有等價性,前兩種方法可視為投影法的特殊情況[13]. 為了便于理解狀態約束對狀態估計及其方差-協方差矩陣的影響,以下采用更一般的投影法的形式來推導附加約束的Kalman濾波遞推公式.

系統的狀態方程及觀測方程如第1.2節的定義. 另外,等式狀態約束條件通常表示為如下線性的約束方程:

Dk·xk=dk,

(13)

式中:Dk為系數矩陣;dk為約束向量.

利用投影法處理帶約束的參數估計問題可表示為如下數學形式[12]:

(14)

上述問題最優解為

xc=x-W-1DT(DW-1DT)-1(Dx-d).

(15)

當過程噪聲和觀測噪聲為高斯的,可以取W的值為無約束條件下估值的誤差協方差矩陣,即W=(Pk)-1,此時估計值具有最小的估計誤差協方差.

采用投影法的約束Kalman濾波算法有兩種形式,可分為開環形式和閉環形式. 具體區別在于,開環形式中,直接利用約束對最后的狀態估計值進行改正,而閉環形式中,需要將每個歷元修正后的估計值反饋給下一個歷元. 兩種形式的流程框圖如圖1所示. 由于閉環投影具有更好的約束性能[4],本文中采取閉環投影的形式.

(a)開環形式

(b)閉環形式圖1 投影法的不同形式

最終得到附加狀態約束的Kalman濾波遞推公式如下:

(16)

Pk,k-1=FPk-1FT+Qk-1,

(17)

(18)

(19)

(20)

(21)

由上面的遞推方程可以看出,將約束條件引入Kalman濾波方程,其思路是在原有濾波的基礎上,利用約束條件對濾波遞推方程中狀態估值的表達式(19)作了修正,以改善濾波結果[10].

進一步研究約束修正后的估值的方差-協方差矩陣. 由式(21)可以得到:

(22)

(23)

根據誤差協方差的定義則有:

=Pk-JPk-PkJT+JPkJT.

(24)

根據J的定義,有如下性質PkJT=JPkJT,且JPk為正定矩陣,由式(24)可以進一步得到:

(25)

(26)

(27)

可以發現,經約束信息修正后的狀態估值的誤差協方差陣小于無約束條件下狀態估值的誤差協方差陣[14],即狀態約束信息可以改善參數估值的誤差-協方差矩陣,從而參數估值的精度得到了提高.

2.2 附加基線約束信息的PPP

考慮兩個觀測站同時進行靜態PPP的定位場景,則這兩個觀測站之間會構成一個基線矢量,這個基線矢量不僅包含了方向信息也包含了長度信息. 因此,如果提前根據相對定位或者其他手段得到兩測站的基線(相對位置坐標),則用各個方向上的相對位置坐標差除以基線長度可以獲得這個基線矢量的方向余弦,即方向信息. 另外,根據基線坐標可以計算得到兩個測站之間的距離信息. 上述方向和距離約束信息本質上提供測站之間的幾何約束關系,將此幾何約束信息融入到Kalman濾波解算就有可能提高每個測站的PPP定位性能.

估值相對應的基線向量的基線長度為

(28)

首先,任何一個矢量的方向都可以由它的方向角來決定,而向量的三個方向余弦分別為這三個方向角的余弦值,因此可用矢量的方向余弦表征矢量的方向性. 設α,β,γ為基線矢量與地心地固坐標系三個坐標軸的方向角,cosα,cosβ,cosγ為對應的三個方向余弦,根據方向余弦的定義:

(29)

則任意k時刻,這兩個測站位置參數滿足如下方向約束:

(30)

(31)

(32)

以X方向的基線方向約束方程為例,可寫出矩陣方程Dk·xk=dk形式表示的基線方向約束:

Dk=

(33)

(34)

其他各方向上的基線方向約束方程構建方法同理.將得到的三個方向上的基線方向約束方程帶入2.1節的Kalman濾波遞推方程則可以得到方向約束后的兩測站坐標.

此外,之前計算得到的兩個測站之間的距離,也可以作為約束對兩天線的位置進行修正. 兩天線之間的先驗距離為S0,則任意k時刻,這兩個測站位置參數滿足如下距離約束:

(35)

上述距離約束方程同樣為非線性方程,需在預測值處進行線性化,只保留一階項后可將上述方程整理為矩陣方程Dk·xk=dk的形式.有:

(36)

(37)

將上述方程帶入2.1節的Kalman濾波遞推方程則可以得到距離約束后的兩測站坐標.

3 數據處理與分析

3.1 數據處理策略

本文以開源的PPP程序PPPH[15]為框架,根據附加約束的Kalman濾波算法原理對程序改進后進行PPP的解算. 數據預處理采用Turbo-Edit方法進行周跳的探測. 誤差改正部分采用了IGS的精密產品和精化的誤差改正模型來進行修正. PPP解算中未進行模糊度固定,參數估值為浮點解. 具體處理策略如表1所示.

表1 數據處理策略

實驗過程中采用15 min采樣率的精密星歷和30 s采樣率的精密衛星鐘差,并采用9階拉格朗日內插對衛星星歷和鐘差進行插值加密,以獲得每一觀測時刻的衛星星歷和鐘差.

3.2 實驗結果分析

為了更加全面地評估附加基線約束的PPP算法的定位性能,實驗利用采集的真實GPS觀測數據進行定位解算. 實驗數據采集于山西省太原市中北大學校園內,接收機采用華測CHCNAV的接收機. 在相距位置約2.5 m的兩個點分別架設兩個天線,在同一時間段內對兩測站進行同步觀測,接收機采樣率為1 s. 觀測時間段為2020年1月2日6:39:26-9:18:15. 截取了其中三個不同時長的數據進行解算,分別為8:00-9:00的1 h觀測數據;7:30-9:00的1.5 h觀測數據;7:00-9:00的2 h觀測數據.

為了確定方向約束和距離約束的約束方程的系數矩陣和約束向量,首先將觀測數據利用RTKLIB軟件static模式解得的基線作為兩天線之間的約束. 解得的基線坐標為(-2.4597,-0.8604,-0.0028),基線長度為2.605 m. 以2.2節所述構建約束方程,得到約束方程后,依次采用以下四種方案進行PPP解算:

1) 采用不附加約束的Kalman濾波;

2) 采用附加基線方向約束的Kalman濾波;

3) 采用附加基線長度約束的Kalman濾波;

4) 采用同時附加基線方向約束和基線長度約束的Kalman濾波.

以兩個天線在各個方向的定位偏差的RMS值為評價指標,評價兩個天線在各個方向的定位精度. 各方向的RMS值由收斂時刻(CT)到觀測截止時刻的坐標偏差計算得到,以連續60 s內三維坐標偏差的2范數不超過0.15 m作為判斷收斂的條件. 得到如表2所示的三個不同時段內兩天線X、Y、Z方向上的RMS統計.

表2 兩個天線定位偏差的RMS統計

分析表2可知,相比于方案1,參數解算過程中附加約束后(方案2～方案4),PPP定位性能得到了不同程度的改善,收斂時間也得到加快. 具體體現在:使用8:00-9:00時間段內數據進行解算,方案2下,1號天線Y坐標,2號天線Y、Z坐標估計精度得到提高,尤其是2號天線Y方向坐標RMS值由0.131 1變為0.029 4,減小了77%;方案3下,1號天線X坐標,2號天線X、Y、Z坐標估計精度得到提高;方案4下,1號天線Y、Z坐標,2號天線Y、Z坐標估計精度得到提高,同樣2號天線Y方向坐標RMS值變化明顯,由0.131 1變為0.033 8,減小了74%. 使用7:30-9:00時間段內的數據進行解算,三種附加約束信息的解算方案兩天線各個方向上定位精度均得到提高. 使用7:00-9:00時間段內的數據進行解算,方案2下,1號天線Y坐標,2號天線Y、Z坐標估計精度得到提高,特別是2號天線Y方向坐標RMS值由0.114 6變為0.037 7,減小了67%;方案3下,1號天線Y坐標,2號天線X、Y、Z坐標估計精度得到提高;方案4下,1號天線Y坐標,2號天線Y、Z坐標估計精度得到提高,同樣2號天線Y方向坐標RMS值變化明顯,由0.114 6變為0.029 6,減小了74%.

綜合來看,通過在PPP解算中附加約束條件,能有效縮短收斂的時間. 不同觀測時長的數據,利用附加約束的PPP解算方案進行解算,收斂時間都得到了縮短. 尤其是方向約束條件(方案2及方案4),對于收斂時間改善效果較為明顯. 三種附加約束的解算方案,均能較大程度提高2號天線Y方向的定位精度. 另外,采用1 h和2 h觀測時長的數據,不同解算方案下,三個方向上定位精度變化不同,但總體上三維方向坐標估計偏差呈減小的趨勢,這點從收斂的判定條件和收斂時間縮短中也可以看出來. 總的來看,三種附加約束的PPP算法定位性能均滿足精密定位的要求,相比無約束定位算法,在不同方向上有不同的提升效果.

為了觀察約束條件對于參數估計的影響過程,將7:30-9:00這一時間段四種方案的解算結果同RTKLIB的PPP模式解得的位置坐標(視為真值)進行比較,可以得到三個方向上4種解算方案的坐標估計誤差變化,如圖2所示.

(a) 天線1-X方向

(b) 天線1-Y方向

(c) 天線1-Z方向

(d) 天線2-X方向

(e) 天線2-Y方向

(f) 天線2-Z方向圖2 不同解算方案坐標估計誤差

另外,將解得的兩天線位置坐標轉換到大地坐標系下,可以得到上述4種解算方案下的兩天線估計位置分布,如圖3所示.

(a) 方案1

(b) 方案2

(c) 方案3

(d) 方案4圖3 大地坐標系下兩天線位置分布散點圖

對比圖2中不同方向上4種解算方案的位置估計誤差圖可以發現,附加約束的PPP解算結果(方案2、方案3、方案4)的精度都不同程度優于無約束的PPP解算結果. 可以發現,不同方向定位結果收斂穩定之前,附加約束的解算方案誤差曲線相比無約束解算方案更加平緩,其對應的定位誤差更小,而估計結果穩定后,4種方案的誤差曲線差距并不明顯. 說明約束條件主要影響前期解算精度不高的時段,通過附加約束能夠較快使坐標估計值符合到真值處,降低定位誤差,從而提高收斂速度. 待解算收斂后定位達到較高精度,約束條件對于估值的作用削弱.

對比圖3中的天線位置散點圖可以發現,無約束時估計得到的天線位置分布較為散亂,具有很大的隨機性. 附加約束后(方案2、方案3、方案4),天線位置分布均產生變化,尤其是方案2及方案4附加方向約束的情況下,得到的天線位置分布具有很強的規律性,估計得到的天線位置分布在兩天線所在直線上. 這也體現了約束條件對于接收機位置估值產生了影響.

另外,綜合對比表2、圖2～3可以發現,當僅考慮距離約束(方案3)的情況下,其定位精度提升并不明顯,估計得到的天線位置分布仍近似于無約束的參數估計結果,較為離散. 而在附加基線長度后,另外附加約束方向約束(方案4),對于解算結果影響較大,其最終的接收機定位分布接近于僅考慮基線方向約束的情況,并且兩天線Y方向定位精度提高不大的情況得到了有效改善. 究其原因可能是因為PPP作為一種精密定位手段,該技術本身的定位精度很高(本實驗中無約束下為厘米級),而RTKLIB相對定位解算出的基線,其精度也在厘米級,因此通過基線坐標求平方和得到的距離約束條件對于兩個天線的坐標改正量較小. 而基線方向約束條件隱含的幾何約束效果對于兩個天線的坐標改正量較大. 在構建約束條件時,特別要注意先驗信息要足夠精確,其精度最好高于無約束定位結果1到2個數量級,可使得約束效果更為明顯. 如可通過采集更長時間的觀測數據或者采用雙差定位的方式,求得毫米級的先驗點坐標信息,從而獲取更準確的基線信息.

4 結束語

針對PPP中采用Kalman濾波進行數據解算時,其性能受限于動態模型和濾波初值的準確程度的情況,本文提出將待估參數之間的函數關系作為約束應用到濾波解算過程中以改善定位性能. 通過實測的GPS觀測數據進行實驗發現,附加天線之間所成基線隱含的方向信息和距離信息,可以不同程度地提高傳統PPP的定位精度,縮短收斂時間. 但是同樣發現,由于約束條件本身精度的限制,限制了約束條件對于參數估值的修正作用. 因此,首先需要保證事先獲得的約束信息有較高精度,從而建立準確的約束條件;另外在選取約束條件時,應更仔細地分析系統本身的特點,找到約束效果更強參數間的函數關系,如直線軌跡約束強于距離約束和平面約束,這是保證約束定位算法解算性能的前提.