,,,

(1.航天工程大學(xué) 研究生院，北京 101416; 2.航天工程大學(xué) 光電裝備系，北京 101416)

0 引言

可觀測性這一概念是由Kalman為了解決線性系統(tǒng)的相關(guān)問題而提出的，如果系統(tǒng)的狀態(tài)能被過去的觀測唯一確定，則該系統(tǒng)為可觀測的[1]。但航天工程領(lǐng)域，各種系統(tǒng)均為非線性的，例如自主導(dǎo)航及濾波定軌系統(tǒng)。對(duì)于非線性時(shí)變系統(tǒng)的可觀測性，還沒有一個(gè)統(tǒng)一而嚴(yán)密的定義。

Lee和Dunn提出的李函數(shù)準(zhǔn)則是目前非線性系統(tǒng)可觀測性判定的常用準(zhǔn)則[2]。僅判斷系統(tǒng)的可觀測性在實(shí)際應(yīng)用中作用并不大，所以引入可觀測度的定義，能將系統(tǒng)的可觀測性大小定量的體現(xiàn)出來。針對(duì)這一問題，國內(nèi)外先后提出多種非線性系統(tǒng)的可觀測度的判斷方法，其中可觀測矩陣SVD分解[3-4]和基于條件數(shù)的可觀測度計(jì)算方法[5-8]最為常見，然而這兩種方案均未考慮量測噪聲對(duì)系統(tǒng)可觀測性造成的影響，李恒年教授及孫仲康教授在文獻(xiàn) [9-10]中均提出了將可觀測矩陣和量測噪聲協(xié)方差矩陣結(jié)合的可觀測度判定方法，所求得的可觀測度受噪聲影響較大。

在多站單測量體制的實(shí)際應(yīng)用條件下，對(duì)濾波定軌系統(tǒng)進(jìn)行可觀測性分析。針對(duì)該種弱觀測條件下可觀測矩陣條件數(shù)過大時(shí)可觀測度判斷不準(zhǔn)確的這一情況，分析了誤差形成的原因并提出了相關(guān)解決方案。建立衛(wèi)星動(dòng)力學(xué)狀態(tài)方程與多站測速的量測方程，并求出其觀測矩陣，采用容積卡爾曼濾波(CKF)算法作為定軌算法，利用Matlab軟件對(duì)多測站濾波定軌系統(tǒng)的可觀測度進(jìn)行仿真。通過比較不同測站組合情況下定軌精度與可觀測度的關(guān)系。仿真結(jié)果表明，改進(jìn)后可觀測度判斷方法相較于傳統(tǒng)的基于條件數(shù)判斷方法，能更準(zhǔn)確判斷系統(tǒng)的定軌精度與可觀測度的對(duì)應(yīng)關(guān)系。

1 系統(tǒng)描述

常用非線性濾波算法如擴(kuò)展卡爾曼濾波(EKF)、無跡卡爾曼濾波(UKF)以及容積卡爾曼濾波(CKF)，適用于非線性系統(tǒng)的最優(yōu)狀態(tài)估計(jì)。在衛(wèi)星軌道估計(jì)中常采用非線性濾波算法進(jìn)行實(shí)時(shí)軌道估計(jì)，是一種在時(shí)域上對(duì)軌道狀態(tài)矢量進(jìn)行最優(yōu)估計(jì)的方法，需要對(duì)系統(tǒng)建立狀態(tài)方程與量測方程，利用狀態(tài)方程進(jìn)行時(shí)間更新同時(shí)對(duì)利用量測方程進(jìn)行量測更新，最終給出最優(yōu)估計(jì)。對(duì)系統(tǒng)進(jìn)一步進(jìn)行可觀測性分析時(shí)，根據(jù)控制理論的要求，需要求出非線性系統(tǒng)的狀態(tài)轉(zhuǎn)移方程與觀測方程。

1.1 狀態(tài)方程

根據(jù)衛(wèi)星在軌道運(yùn)行的動(dòng)力學(xué)規(guī)律，建立衛(wèi)星的軌道動(dòng)力學(xué)模型。地球是一個(gè)形狀、質(zhì)量分布皆非均勻的扁球體，對(duì)航天器的引力需考慮地球非球形引力，大氣阻力和太陽光壓、潮汐等其他攝動(dòng)力。非球形引力項(xiàng)中最大的是J2項(xiàng)，數(shù)量級(jí)遠(yuǎn)遠(yuǎn)大于其攝動(dòng)力，在計(jì)算中僅考慮該項(xiàng)。在地心慣性坐標(biāo)系中，單顆衛(wèi)星的運(yùn)動(dòng)方程可以表示為引力、地心距離與衛(wèi)星運(yùn)動(dòng)矢量的微分方程。在J2000坐標(biāo)系下建立狀態(tài)方程，如下所示：

(1)

在計(jì)算機(jī)仿真分析中，連續(xù)系統(tǒng)無法進(jìn)行時(shí)域上的狀態(tài)估計(jì)，于是將式(1)的連續(xù)方程通過數(shù)值計(jì)算方法離散化，常用方法有歐拉法、龍格庫塔離散方法等，再對(duì)系統(tǒng)進(jìn)行泰勒一階展開，變?yōu)?2)所示的離散非線性時(shí)變系統(tǒng)：

Xk+1=Φk+1,kXk

(2)

式中,

(3)

Xk+1為k時(shí)刻的航天器的狀態(tài)量，Xk為k時(shí)刻的狀態(tài)量，Φk+1,k為k時(shí)刻到k+1時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣。

1.2 量測方程

量測過程是通過地面站對(duì)航天器進(jìn)行測量，常用的測量元有測距、測速及測角。在實(shí)際應(yīng)用中，相較于包含完整測元的體制，僅采用單測角或單測速的系統(tǒng)軌道估計(jì)精度較低，更具有研究意義，因此本文采用單種測速值的測量系統(tǒng)，對(duì)其進(jìn)行可觀測性研究。

每個(gè)測站給出的徑向速度量測值通過航天器的狀態(tài)矢量表示出，由此構(gòu)成量測方程組。多個(gè)測站對(duì)同一航天器進(jìn)行觀測獲取多個(gè)測速值，在單站的情況下進(jìn)行擴(kuò)維運(yùn)算，構(gòu)成多站僅測速系統(tǒng)的量測方程：

(4)

利用泰勒一階展開對(duì)量測方程進(jìn)行處理，對(duì)系統(tǒng)狀態(tài)求偏導(dǎo)，獲得系統(tǒng)觀測矩陣[11]：

(5)

量測方程建立在站心坐標(biāo)系中，狀態(tài)方程建立在J2000坐標(biāo)系中，進(jìn)行濾波運(yùn)算時(shí)需要提前將數(shù)據(jù)進(jìn)行坐標(biāo)系轉(zhuǎn)換，在同一坐標(biāo)系下進(jìn)行軌道估計(jì)。

2 可觀測性分析

可觀性及可控性是線性系統(tǒng)的一種定性分析方法，指系統(tǒng)內(nèi)的狀態(tài)是否可以由輸出和輸入進(jìn)行分析，對(duì)于系統(tǒng)的運(yùn)行有重要研究意義。在控制理論研究中，線性系統(tǒng)可觀測性有明確的判斷準(zhǔn)則，如格拉姆矩陣判據(jù)。任何線性系統(tǒng)的定理推廣到非線性系統(tǒng)時(shí)，復(fù)雜程度大大增加。因?yàn)榉蔷€性系統(tǒng)中某一時(shí)刻的狀態(tài)無法由初始狀態(tài)推導(dǎo)，非線性系統(tǒng)可觀測性分析沒有統(tǒng)一的判斷準(zhǔn)則。

2.1 非線性系統(tǒng)的可觀測性分析

非線性系統(tǒng)可通過泰勒展開等近似為線性化系統(tǒng)，再采用線性系統(tǒng)的判斷準(zhǔn)則來判斷，但該種方法對(duì)犧牲了對(duì)系統(tǒng)的準(zhǔn)確描述使計(jì)算變得簡明。而直接對(duì)非線性系統(tǒng)的可觀測性判定在工程運(yùn)用中，通常采用由Lee和Dunn提出的李函數(shù)定理[2]：

對(duì)非線性系統(tǒng):

(6)

y(t)=h(x(t),t)

(7)

如果在t∈[t0,t1]上，對(duì)凸集S∈Rn上的任意x0，x1∈S，h(x(t;x),t)=h(x(t;x1),t)， 都有x1=x0，則系統(tǒng)在凸集S上是完全可觀測的。線性系統(tǒng)在S集上的可觀測性可以通過對(duì)初始狀態(tài)的可觀性分析推導(dǎo)出來，與其狀態(tài)向量的變化無關(guān)。非線性系統(tǒng)的可觀性不能由初始狀態(tài)推導(dǎo)，因此考慮由狀態(tài)方程f(·)和量測方程h(·)進(jìn)行推導(dǎo)。

定義：

(8)

(9)

式中，H(t)、Φ(t,t0)為h(·)和f(·)的雅各比矩陣。若M(x0)是正定的，則系統(tǒng)在凸集S上是完全可觀測的。

將該定理應(yīng)用在非線性離散時(shí)變系統(tǒng)中，定義可觀測性矩陣Γ為：

(10)

式中，N表示觀測次數(shù)，n表示狀態(tài)維數(shù)。

如果rankΓ(k0,k0+N-1)=n，即在N次觀測中，可觀測矩陣Γ的秩等于n，則系統(tǒng)在S上是完全可觀測的。將式(3)、(5)帶入(10)可以得到可觀測矩陣，對(duì)其求秩可以對(duì)系統(tǒng)的可觀測性進(jìn)行判斷，通過求取特征值或進(jìn)行奇異值分解，都可以對(duì)系統(tǒng)的可觀測度進(jìn)行判斷。

系統(tǒng)中的狀態(tài)噪聲與量測噪聲對(duì)系統(tǒng)的可觀測性也存在影響，對(duì)系統(tǒng)本身結(jié)構(gòu)的可觀測性描述可以不考慮噪聲的影響。考慮噪聲時(shí)的系統(tǒng)可觀測性及可觀測度判斷有其他的判斷方法。

2.2 可觀測度分析方法

可觀測性定性的反應(yīng)出系統(tǒng)的實(shí)際運(yùn)行能力，為了更清楚反應(yīng)系統(tǒng)可觀測性的大小，引入可觀測度的概念，對(duì)系統(tǒng)的可觀測性給出一種定量的分析方法。

矩陣的條件數(shù)定義為最大奇異值和最小奇異值的比值，可觀測矩陣的條件數(shù)反映可觀測矩陣的病態(tài)程度，可以理解為可觀測矩陣各行列之間的相關(guān)性，相關(guān)性強(qiáng)表明量測值中獲取的相似信息量多，有效信息量少；相關(guān)性弱則表明量測值中獲取的有效信息量多。這樣能直觀的反映出算法的穩(wěn)定性、收斂性及收斂速度，將條件數(shù)的大小作為可觀測度的判斷依據(jù)。

對(duì)可觀測矩陣Γ進(jìn)行式(11)所示的奇異值分解：

Γ=U∑V

(11)

式中，U、V為奇異值分解后的酉矩陣，∑為對(duì)角矩陣，對(duì)角線上的元素為可觀測矩陣Γ的奇異值。條件數(shù)cond(Γ)為最大奇異值和最小奇異值的商。

(12)

在實(shí)際應(yīng)用及仿真中發(fā)現(xiàn)，式(11)、(12)直接計(jì)算出的條件數(shù)并不適用于多站單測量體制的濾波定軌算法。傳統(tǒng)的條件數(shù)計(jì)算中將不同測站獲取的同種體制的量測值直接帶入可觀測矩陣Γ計(jì)算，所得到的可觀測度僅僅是量測信息本身之間數(shù)值大小上的相關(guān)性，當(dāng)兩個(gè)處于不同位置的測站所獲取的量測值在僅僅在數(shù)值上相似時(shí)，條件數(shù)會(huì)猛然增大，可知其并不能正確體現(xiàn)不同測站所包含有效信息的多少，定軌精度的好壞與所得到的條件數(shù)大小并不存在對(duì)應(yīng)關(guān)系。

(13)

(14)

3 仿真分析

為了研究可觀測度是否能作為濾波定軌算法精度的參考依據(jù)，采用基于容積卡爾曼濾波的濾波定軌算法，在衛(wèi)星工具包(Satellite Tool Kit，STK)中建立場景模型，衛(wèi)星模型采用太陽同步軌道衛(wèi)星(Sun-synchronous satellite)，軌道估計(jì)數(shù)據(jù)在軌道估計(jì)算法中采用高精度軌道預(yù)報(bào)(high precision orbit propagation, HPOP)，在條件數(shù)計(jì)算中采用僅考慮J2項(xiàng)的軌道數(shù)據(jù)，更利于狀態(tài)方程的雅各比矩陣的計(jì)算同時(shí)仿真出的曲線平滑，最終利用MATLAB進(jìn)行算法驗(yàn)證。

對(duì)模型中隨機(jī)放置的兩個(gè)測站和隨機(jī)放置的三個(gè)測站在不同幾何布設(shè)的情況下計(jì)算優(yōu)化方法所得的條件數(shù)，同時(shí)進(jìn)行定軌算法仿真以得出不同測站組合下的定軌誤差結(jié)果，濾波步長為1 s，條件數(shù)計(jì)算時(shí)長為100 s，定軌算法計(jì)算時(shí)長為170 s，在仿真模型中F1、F2、F3、F4分別代表四個(gè)隨機(jī)測站，四個(gè)測站的經(jīng)緯度如表1所示。

表1 四個(gè)測站的經(jīng)緯度

仿真初始軌道值在模型真實(shí)軌道值的基礎(chǔ)上位置量加入1 km的誤差，速度量加上100 m/s的誤差，參考真實(shí)初始軌道值為：

初始協(xié)方差矩陣為：

狀態(tài)噪聲協(xié)方差矩陣為：

選擇通過均方根誤差判斷定軌精度，根據(jù)數(shù)據(jù)的量級(jí)分為速度均方根誤差和位置均方根誤差。

Eposition=

(15)

Evelocity=

(16)

圖1 雙測站的條件數(shù)變化曲線

圖2 雙測站軌道估計(jì)位置均方根誤差

圖3 雙測站軌道估計(jì)速度均方根誤差

圖4 三測站的條件數(shù)變化曲線

圖5 三測站軌道估計(jì)位置均方根誤差

圖6 三測站軌道估計(jì)速度均方根誤差

由圖1和圖4可以看出，由于系統(tǒng)是時(shí)變的，所以系統(tǒng)的條件數(shù)也是時(shí)變的，結(jié)合系統(tǒng)的條件數(shù)與圖2、圖3以及圖4、圖5的定軌精度誤差結(jié)果可以看出，條件數(shù)越小，定軌精度越高，算法的收斂速度更快。為了更直觀的分析條件數(shù)與定軌精度之間的關(guān)系，如表2所示，將傳統(tǒng)方法所求的各系統(tǒng)的條件數(shù)均值及改進(jìn)方法所求的條件數(shù)均值與均方誤差值進(jìn)行對(duì)比，比較不同方法求得的條件數(shù)大小與系統(tǒng)定軌精度的關(guān)系。

表2 定軌精度與兩種可觀測度方法所求得條件數(shù)的比較

圖7 兩種方法所得條件數(shù)與估計(jì)精度的變化趨勢

由表2和圖7可以看出，傳統(tǒng)方法所求得的條件數(shù)的大小無法準(zhǔn)確反映系統(tǒng)的定軌精度，并不適用于條件數(shù)大，量測值之間存在不同數(shù)值意義的弱觀測系統(tǒng)。改進(jìn)方法將不同測站坐標(biāo)作為向量進(jìn)行投影處理，代入可觀測矩陣所求條件數(shù)，仿真結(jié)果數(shù)值上大于傳統(tǒng)方法，但是對(duì)系統(tǒng)可觀測性的判斷準(zhǔn)確度高。因此，實(shí)驗(yàn)結(jié)果表明改進(jìn)方法所求得的條件數(shù)大小能正確反映出定軌精度的變化，當(dāng)條件數(shù)較大時(shí)，系統(tǒng)的定軌精度較差；條件數(shù)較小時(shí)，系統(tǒng)的定軌精度較高。

4 結(jié)論

研究了基于多站測速的濾波定軌算法的可觀測性及可觀測度判斷方法。傳統(tǒng)可觀測度判定方法中條件數(shù)的計(jì)算方法其結(jié)論對(duì)于多站單測量值的弱觀測系統(tǒng)不適用，改進(jìn)方法將不同測站獲取的量測值看作矢量進(jìn)行坐標(biāo)系轉(zhuǎn)換，將測站坐標(biāo)值帶入計(jì)算，給出一種適應(yīng)多站組合測量弱觀測條件下的基于條件數(shù)的可觀測度判定方法。數(shù)值仿真結(jié)果驗(yàn)證了改進(jìn)方法相對(duì)于傳統(tǒng)方法，有更好的估計(jì)效果，更能準(zhǔn)確作為系統(tǒng)可觀測度判定的方法。

參考文獻(xiàn)：

[1] Kalman R E. A New Approach to Linear Filtering and Prediction Problems.[J]. J.basic Eng.trans.ASME, 1960; 82D(1):35-45.

[2] Lee T, Dunn K, Chang C. On observability and unbiased estimation of nonlinear systems[J]. System Modeling and Optimization, 1982: 258-266.

[3] 周衛(wèi)東, 蔡佳楠, 孫 龍,等. GPS/SINS超緊組合導(dǎo)航系統(tǒng)可觀測性分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2013, 39(9):1157-1162.

[4] 黃翔宇, 崔平遠(yuǎn), 崔祜濤. 深空自主導(dǎo)航系統(tǒng)的可觀性分析[J]. 宇航學(xué)報(bào), 2006, 27(3):332-337.

[5] 寧曉琳, 房建成. 航天器自主天文導(dǎo)航系統(tǒng)的可觀測性及可觀測度分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2005, 31(6):673-677.

[6] 劉 準(zhǔn), 陳 哲. 條件數(shù)在系統(tǒng)可觀測性分析中的應(yīng)用研究[J]. 系統(tǒng)仿真學(xué)報(bào), 2004, 16(7): 1552-1555.

[7] 曲 毅, 劉 忠, 薛 鋒. 基于TDOA的雙站自適應(yīng)濾波算法及可觀測性研究[J]. 系統(tǒng)仿真學(xué)報(bào), 2007, 19(22):5222-5225.

[8] 常曉華, 崔平遠(yuǎn), 崔祜濤. 一種深空自主導(dǎo)航系統(tǒng)可觀測性分析方法[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2010, 42(11):1681-1685.

[9] 黃 普, 錢 山, 李恒年. 空間目標(biāo)地基雷達(dá)跟蹤可觀測性分析[J]. 宇航動(dòng)力學(xué)學(xué)報(bào), 2015, 5(3).

[10] 孫仲康,周一宇,何黎星.單/多基地有源無源定位技術(shù)[M].北京：國防工業(yè)出版社,1996.

[11] 李恒年. 航天測控最優(yōu)估計(jì)方法[M]. 北京：國防工業(yè)出版社, 2015.