導航星座整網(wǎng)自主定軌的動力學短弧段偽逆平差法

李灝霖，王 玲，黃文德，周一帆

(1. 湖南大學電氣與信息工程學院，長沙 410012；2. 國防科技大學機電工程與自動化學院，長沙 410073； 3. 中國科學院上海天文臺，上海 200030)

0 引 言

隨著我國“北斗三號”全球衛(wèi)星導航系統(tǒng)建設(shè)的推進，導航星座借助星間鏈路實現(xiàn)自主運行的研究也成為了當前導航領(lǐng)域的研究熱點[1-3]。

由于星間鏈路是存在于導航衛(wèi)星之間的相對測量，在沒有外部基準的情況下，僅通過星間觀測進行星座整網(wǎng)位置平差，會出現(xiàn)秩虧問題，無法實現(xiàn)星座的絕對定向和定位。為了解決秩虧問題，蔡志武等[4]提出了解決星座自主定軌基準秩虧問題的參數(shù)約束方法及星間測角資料輔助法。文獻[5]采用借助中繼星進行聯(lián)合定軌的策略，消除秩虧及軌道整體漂移問題。文獻[6-7]利用上注星歷作為升交點先驗信息，提供位置基準約束，進行自主定軌。文獻[8]中采用增加地面錨固站的方法提供基準。上述方法主要通過外部基準或者精確的地面先驗信息解決秩虧問題，當外部基準不可用或者先驗信息精度不滿足要求時存在一定局限性。另一方面，僅采用單歷元距離觀測進行星座整網(wǎng)平差時，無法對速度分量進行有效估計。文獻[9-10]對該現(xiàn)象進行了簡單闡述，并指出通過增加速度觀測或動力學處理方法解決，但沒有進行展開分析。

本文提出一種自由網(wǎng)偽逆平差理論和動力學短弧定軌技術(shù)相結(jié)合的自主定軌算法。采用大地測量平差中附加最小范數(shù)條件的最小二乘方法，即偽逆平差法解決整網(wǎng)位置解算起算數(shù)據(jù)不足帶來的秩虧問題；采用動力學短弧定軌方法，通過狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣建立觀測值與待估參數(shù)之間的關(guān)系解決速度估計問題。以典型Walker導航星座為例，仿真了全星座星間鏈路測量數(shù)據(jù)，對算法性能進行了校驗。

1 問題描述

1.1 導航衛(wèi)星整網(wǎng)定軌模型

Li,j=h(Xi,Xj)=

(1)

式中:ε為星間測量誤差。

(2)

(3)

(4)

式中:Vi,j為測量殘差，A為系數(shù)矩陣。僅采用單歷元觀測數(shù)據(jù)進行定軌時，A可展開為

(5)

其中，

(6)

通過以上步驟的線性化，已經(jīng)將觀測數(shù)據(jù)和待估量的形式統(tǒng)一到定軌的殘差模型中

V=AδX-l

(7)

殘差最小二乘原理通過求取最小化誤差的平方和尋找超定方程的解，即

(8)

其中，P為觀測權(quán)矩陣。根據(jù)式(8)可以得到法方程

NδX-U=0

(9)

結(jié)果為

δX=N-1U

(9)

式中:N=ATPA，U=ATPl。

1.2 整網(wǎng)位置解算的秩虧問題

在有地面站或錨固站支持的情況下，衛(wèi)星整網(wǎng)定軌以地面站作為定軌網(wǎng)絡(luò)中的基準，式(9)能夠得到一組唯一的最優(yōu)解。但是，在自主運行模式下，衛(wèi)星導航系統(tǒng)脫離地面運行。由式(5)的系數(shù)矩陣可知，當兩顆衛(wèi)星的測量為相對測量時，矩陣中對應行的這兩顆衛(wèi)星的線性化偏導數(shù)互為相反數(shù)。因而，整網(wǎng)觀測中，所有衛(wèi)星之間都存在著相關(guān)性。當沒有外部基準且星座的所有節(jié)點的位置均未知時，會因為缺少必要的起算數(shù)據(jù)使得系數(shù)矩陣A列秩虧，進而導致法方程系數(shù)矩陣N秩虧，無法通過常規(guī)計算對N求逆。因此，由于星間觀測對衛(wèi)星網(wǎng)整體運動的不可觀測性，其衛(wèi)星網(wǎng)只能構(gòu)成空間的相對位置約束，整個衛(wèi)星網(wǎng)缺乏一個絕對的參考基準，法方程系數(shù)矩陣秩虧，只通過經(jīng)典自由網(wǎng)平差無法得到最小二乘準則下的唯一解。

1.3 整網(wǎng)平差速度改正的問題

由式(1)可知,

僅利用單歷元距離觀測數(shù)據(jù)進行整網(wǎng)平差時，沒有對速度的直接或間接觀測數(shù)據(jù),觀測值對速度的偏導數(shù)均為0，即

因此，在對衛(wèi)星狀態(tài)的修正值δX的求解中，無法得到速度修正值，導致速度概略狀態(tài)無法修正。文獻[9]提出兩種求取速度修正值的方法：1)星間觀測僅采用測距信息，對衛(wèi)星位置狀態(tài)修正值進行求取，利用動力學方程獲得速度修正值。2)星間觀測增加測速信息。

目前，衛(wèi)星導航系統(tǒng)多采用基于時間到達原理的星間距離測量作為星間鏈路的基本觀測數(shù)據(jù)。同時，由于同軌道面衛(wèi)星相對速度小，利用多普勒頻移進行速度的相對測量精度較低，并不能滿足速度狀態(tài)更新的要求。因此，本文采用短弧段定軌，利用動力學模型求解狀態(tài)轉(zhuǎn)移矩陣，通過位置狀態(tài)間接估計速度狀態(tài)從而實現(xiàn)高精度定軌。

2 導航星座動力學短弧段整網(wǎng)平差定軌

2.1 偽逆平差方法

在起算數(shù)據(jù)不足的情況下，通過經(jīng)典自由網(wǎng)平差無法得到最小二乘準則下的唯一解。下面給出利用偽逆平差得到秩虧方程的解的方法。

設(shè)滿足法方程(9)的一個最小二乘解為

則該解的范數(shù)為

(10)

(11)

設(shè)矩陣N有一廣義逆N+滿足如下四個條件：

NN+N=N

(N+N)T=N+N

N+NN+=Ν+

(NN+)T=NN+

則稱廣義逆N+為N的偽逆?？梢姡瑐文嬉驗闈M足最小范數(shù)逆的兩個條件，所以為最小范數(shù)逆的一種，同樣能得到的最小范數(shù)解

(12)

同時，偽逆N+唯一，可通過矩陣SVD分解求得?？梢宰C明，δXm的協(xié)因數(shù)矩陣為

QX=N+NN+=N+

2.2 動力學短弧段整網(wǎng)平差

衛(wèi)星的動力學模型可以表示為

(13)

式中:f為衛(wèi)星的受力模型，對于中高軌道衛(wèi)星一般考慮地球非球形攝動，日月引力攝動，太陽光壓攝動等攝動因素，W為動力學模型噪聲。將式(13)展開得

ΔW=F(t)·δX(t)+δW

(14)

式(14)的一般解為

δX(t)=Φ(t,t0)·δX(t0)

(15)

其中，Φ(t,t0)通過求解如下變分方程得到

且一般無法求得解析解，只能通過數(shù)值積分進行求解。

本文采用的短弧定軌的基本方法是，將某觀測弧段上得到的PNum個歷元的觀測值通過狀態(tài)轉(zhuǎn)移矩陣歸算到某一歷元(例如該弧段的初始歷元)。若觀測采樣間隔固定為TStep，則弧段長度TP滿足

TP=(PNum-1)TStep

(16)

從當前弧段起始歷元t0開始積分，得到弧段內(nèi)t1,t2,…,tk,…,t(PNum-1)歷元的概略狀態(tài)及相對于起始歷元的狀態(tài)轉(zhuǎn)移矩陣Φ(t1,t0),Φ(t2,t0),…,Φ(tk,t0),…,Φ(t(PNum-1),t0)。同時，通過t0,t1,…,tk,…,t(PNum-1)歷元的觀測數(shù)據(jù)，可以得到PNum組系數(shù)矩陣A(t0),A(t1),…,A(tk),…,A(t(PNum-1))。利用狀態(tài)轉(zhuǎn)移矩陣，將弧段內(nèi)所有歷元的系數(shù)矩陣統(tǒng)一歸算至弧段起始歷元t0，即

而通過式(15)中Φ(t,t0)和狀態(tài)估計值的關(guān)系可知

(17)

所以

A(tk,t0)=A(tk)·Φ(t,t0)

(18)

由此，可以得到式(7)殘差的變換形式

(19)

式中:δX(t0)為弧段起始歷元整網(wǎng)所有衛(wèi)星狀態(tài)修正值。求出δX(t0)后，通過式(15)對弧段內(nèi)其他歷元的狀態(tài)修正值δX(t)進行求解。

2.3 整網(wǎng)短弧段偽逆平差定軌算法流程

根據(jù)短弧段偽逆平差原理，首先確定弧段起始歷元，并確定起始點的概略狀態(tài)X-(t0)。然后進行整網(wǎng)短弧段偽逆平差定軌。整網(wǎng)短弧段偽逆平差定軌算法為本文研究的核心部分，主要包括短弧段數(shù)據(jù)累積、整網(wǎng)平差兩部分，定軌算法及執(zhí)行步驟如下所示。

步驟1：短弧段數(shù)據(jù)累積，包括以下步驟：

步驟1.1：獲取弧段內(nèi)起始歷元的觀測數(shù)據(jù)L(t0)，如式(1)，令k=1。

步驟1.2：如式(3)，求取起始歷元的概略觀測L-(t0)，同時，將起始歷元觀測數(shù)據(jù)L(t0)在X-(t0)處進行泰勒展開，得到起始歷元系數(shù)矩陣A(t0)。

步驟1.3：通過軌道積分獲得tk歷元的概略位置X-(tk)以及該歷元相對于弧段起始歷元的狀態(tài)轉(zhuǎn)移矩陣Φ(tk,t0)。

步驟1.4：獲得tk歷元觀測數(shù)據(jù)L(tk)。

步驟1.5：同步驟1.2，求取tk歷元概略觀測L-(tk)和系數(shù)矩陣A(tk)。

步驟1.6：通過式(18)，利用狀態(tài)轉(zhuǎn)移矩陣Φ(tk,t0)將歷元系數(shù)矩陣歸算至起始歷元，得A(tk,t0)。

步驟1.7：判斷tk是否為弧段終止歷元，若不是，則k=k+1且返回步驟1.3。

步驟1.8：整合弧段內(nèi)所有歷元觀測數(shù)據(jù)和概略觀測可得矩陣L=[L(t0),L(t1), …,L(tk), …,L(t(PNum-1))]和矩陣L-=[L-(t0),L-(t1), …,L-(tk), …,L-(t(PNum-1))],從而可得殘差矩陣l=L-L-，以及整弧段系數(shù)矩陣A=[A(t0,t0),A(t1,t0), …,A(tk,t0), …,A(t(PNum-1),t0)]。

步驟2：整網(wǎng)平差，包括以下步驟：

步驟2.1：構(gòu)造如式(19)的最小二乘整網(wǎng)定軌模型，并得到法方程，如式(9)。

步驟2.2：對法矩陣N進行SVD分解，求得法矩陣的偽逆N+，并通過式(12)求得弧段內(nèi)起始歷元的狀態(tài)修正值δX(t0)。

步驟2.3：通過式(15)，利用δX(t0)和狀態(tài)轉(zhuǎn)移矩陣Φ(t,t0)求取弧段內(nèi)除起始歷元外的其他歷元的狀態(tài)修正值，可得到δX，并加到概略位置上，得到全弧段軌道狀態(tài)估計值。

3 仿真結(jié)果與分析

3.1 仿真算例

利用STK仿真衛(wèi)星星座，Walker星座參數(shù)為24/6/1，所有衛(wèi)星軌道統(tǒng)一高度為25678 km，軌道偏心率為0，軌道傾角55°，近地點角距為0°，其他初值不同的軌道6根數(shù)如表 1所示。動力學模型中，除地球中心引力外，還考慮10×10階地球非球形攝動、日月引力攝動以及太陽光壓攝動。其中，重力場模型采用JGM-3全球重力場模型；太陽光壓模型采用球模型，面質(zhì)比取0.002 m2/kg，反射系數(shù)取1.2，考慮錐形地影模型。本文僅進行衛(wèi)星定軌，不進行鐘差估計。

仿真起始歷元，各衛(wèi)星位置狀態(tài)加入均方根為1 m正態(tài)分布的隨機誤差，速度加入均方根為10-6m/s正態(tài)分布的隨機誤差。同時，仿真的觀測數(shù)據(jù)在理論距離的基礎(chǔ)上加入均方根0.2 m的隨機誤差。

表1 星座軌道根數(shù)初值Table 1 Part of initial value of constellation

注：M表示對應衛(wèi)星編號的平近點角； RAAN表示升交點經(jīng)度。

3.2 結(jié)果分析

在進行整網(wǎng)定軌中以用戶測距誤差值(User range error,URE)作為定軌精度判定指標[11]，其定義為：

ΔURE(i)=

其中，ΔURE(i)表示第i顆衛(wèi)星的URE值，R(i)表示徑向方向誤差，C(i)表示鐘差誤差，本文不對鐘差進行估計，因而鐘差為0，T(i)表示跡向方向誤差，N(i)表示與徑向、跡向平面垂直方向的誤差，ΔUREM表示星座中n顆衛(wèi)星的平均URE值。

本文以15 min為一個歷元，采用6歷元一個弧段進行自主定軌。依據(jù)式(16)，對應弧段長度為75 min。

若不進行自主定軌，僅通過動力學模型進行軌道預報90天，誤差如圖 4所示，24顆衛(wèi)星ΔUREM約600 m。

圖5為采用本文自主定軌算法運行90天，不同歷元24顆衛(wèi)星ΔUREM。從圖 5可以看出，衛(wèi)星整網(wǎng)星座自主運行90天，定軌結(jié)果未發(fā)現(xiàn)明顯的發(fā)散過程，定軌精度穩(wěn)定，ΔUREM曲線下包絡(luò)精度基本控制在0.5 m以下。經(jīng)統(tǒng)計，24顆衛(wèi)星ΔUREM，90天的均值為0.54 m，標準差為0.23 m，其中最大值為1.98 m，最小值僅0.054 m。

從圖6可以看出，RTN三個方向誤差都比較穩(wěn)定，與ΔUREM的變化趨勢基本一致，T方向誤差變化比較明顯，是導致ΔUREM曲線下包絡(luò)變化的主要原因。

4 結(jié)束語

本文針對傳統(tǒng)靜態(tài)自由網(wǎng)平差進行自主定軌存在的秩虧問題，提出一種自由網(wǎng)偽逆平差理論和動力學短弧段定軌技術(shù)相結(jié)合的自主定軌算法。該算法在傳統(tǒng)最小二乘解算的通解基礎(chǔ)上增加最小范數(shù)約束得到唯一位置解，解決整網(wǎng)位置解算起算數(shù)據(jù)不足帶來的秩虧問題；同時，利用狀態(tài)轉(zhuǎn)移矩陣對系數(shù)矩陣進行歸算，通過距離觀測間接修正速度，解決單歷元平差無法對速度分量進行有效估計的問題。本文方法可應用于我國“北斗三號”全球衛(wèi)星導航系統(tǒng)的自主運行。