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

(1. 湖南大學(xué)電氣與信息工程學(xué)院，長(zhǎng)沙 410012；2. 國(guó)防科技大學(xué)機(jī)電工程與自動(dòng)化學(xué)院，長(zhǎng)沙 410073； 3. 中國(guó)科學(xué)院上海天文臺(tái)，上海 200030)

0 引 言

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

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

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

1 問(wèn)題描述

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

Li,j=h(Xi,Xj)=

(1)

式中:ε為星間測(cè)量誤差。

(2)

(3)

(4)

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

(5)

其中，

(6)

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

V=AδX-l

(7)

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

(8)

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

NδX-U=0

(9)

結(jié)果為

δX=N-1U

(9)

式中:N=ATPA，U=ATPl。

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

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

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

由式(1)可知,

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

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

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

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

2.1 偽逆平差方法

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

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

則該解的范數(shù)為

(10)

(11)

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

NN+N=N

(N+N)T=N+N

N+NN+=Ν+

(NN+)T=NN+

則稱廣義逆N+為N的偽逆。可見(jiàn)，偽逆因?yàn)闈M足最小范數(shù)逆的兩個(gè)條件，所以為最小范數(shù)逆的一種，同樣能得到的最小范數(shù)解

(12)

同時(shí)，偽逆N+唯一，可通過(guò)矩陣SVD分解求得。可以證明，δXm的協(xié)因數(shù)矩陣為

QX=N+NN+=N+

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

衛(wèi)星的動(dòng)力學(xué)模型可以表示為

(13)

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

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

(14)

式(14)的一般解為

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

(15)

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

且一般無(wú)法求得解析解，只能通過(guò)數(shù)值積分進(jìn)行求解。

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

TP=(PNum-1)TStep

(16)

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

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

(17)

所以

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

(18)

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

(19)

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

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

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

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

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

步驟1.2：如式(3)，求取起始?xì)v元的概略觀測(cè)L-(t0)，同時(shí)，將起始?xì)v元觀測(cè)數(shù)據(jù)L(t0)在X-(t0)處進(jìn)行泰勒展開(kāi)，得到起始?xì)v元系數(shù)矩陣A(t0)。

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

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

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

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

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

步驟1.8：整合弧段內(nèi)所有歷元觀測(cè)數(shù)據(jù)和概略觀測(cè)可得矩陣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：對(duì)法矩陣N進(jìn)行SVD分解，求得法矩陣的偽逆N+，并通過(guò)式(12)求得弧段內(nèi)起始?xì)v元的狀態(tài)修正值δX(t0)。

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

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

3.1 仿真算例

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

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

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

注：M表示對(duì)應(yīng)衛(wèi)星編號(hào)的平近點(diǎn)角； RAAN表示升交點(diǎn)經(jīng)度。

3.2 結(jié)果分析

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

ΔURE(i)=

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

本文以15 min為一個(gè)歷元，采用6歷元一個(gè)弧段進(jìn)行自主定軌。依據(jù)式(16)，對(duì)應(yīng)弧段長(zhǎng)度為75 min。

若不進(jìn)行自主定軌，僅通過(guò)動(dòng)力學(xué)模型進(jìn)行軌道預(yù)報(bào)90天，誤差如圖 4所示，24顆衛(wèi)星ΔUREM約600 m。

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

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

4 結(jié)束語(yǔ)

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