多星混合資料的定軌?

王 歆

(1中國科學(xué)院紫金山天文臺南京210008)(2中國科學(xué)院空間目標(biāo)與碎片觀測重點(diǎn)實(shí)驗(yàn)室南京210008)

多星混合資料的定軌?

王 歆1,2?

(1中國科學(xué)院紫金山天文臺南京210008)(2中國科學(xué)院空間目標(biāo)與碎片觀測重點(diǎn)實(shí)驗(yàn)室南京210008)

在空間目標(biāo)光學(xué)觀測資料定軌中常常會遇到多個目標(biāo)的觀測被標(biāo)記為同一個目標(biāo)的情況,由于包含了多個目標(biāo)的數(shù)據(jù),定軌過程無法收斂或者完全錯誤．從極大似然估計(jì)角度,采用EM(Expectation Maximum)方法提出一種將軌道改進(jìn)和識別過程相互融合的處理方法,并在具體實(shí)現(xiàn)過程中給出一種穩(wěn)健估計(jì)方法．?dāng)?shù)值模擬表明方法簡便、有效、可行．

航天器,天體力學(xué),方法:統(tǒng)計(jì)

1 引言

隨著航天活動的深入開展與航天技術(shù)的發(fā)展,衛(wèi)星編隊(duì)越來越多地被采用,編隊(duì)飛行目標(biāo)在光學(xué)觀測中會出現(xiàn)在臨近天區(qū)．由于空間目標(biāo)預(yù)報誤差,在觀測中往往會采用提前擺位等方式來進(jìn)行空間目標(biāo)的捕獲,當(dāng)目標(biāo)先后經(jīng)過視場時,由于視運(yùn)動相近,以及光學(xué)觀測干擾因素較多,實(shí)踐中很難正確判斷待捕獲目標(biāo)．對于大視場空間目標(biāo)觀測望遠(yuǎn)鏡,編隊(duì)飛行的目標(biāo)甚至?xí)瑫r出現(xiàn)在視場中,觀測中也難以正確識別．

這種情況使得采集的觀測資料中,雖然都按照預(yù)報進(jìn)行了標(biāo)記,但不同站圈實(shí)際上跟蹤了不同的目標(biāo),形成了本文提出的多星混合資料,這里混合有“混淆”的含義．按照標(biāo)記目標(biāo)對這些資料進(jìn)行軌道改進(jìn)時,往往無法收斂.過去對于這類資料多采取人工處理,依賴經(jīng)驗(yàn)并沒有成熟的方法.現(xiàn)在觀測能力與數(shù)量大幅度提高后,此類情況也大幅度增加,依賴人工已無法滿足需求.即便采用軌道識別方法,對于編隊(duì)目標(biāo)由于目標(biāo)過于接近而無法區(qū)分[1].如何自動處理這類資料是一個棘手的問題.

這個問題表面上可歸為穩(wěn)健估計(jì)的范疇,也就是將標(biāo)記錯誤的資料作野值處理.然而該問題并不能簡單處理,一方面從資料屬性上只是屬于兩個目標(biāo),而并不是“壞”的觀測;另一方面也無法確保正確標(biāo)記的觀測量占到多數(shù),因此單純從穩(wěn)健估計(jì)角度是不能解決問題的.

本文從另一個角度對問題進(jìn)行了細(xì)致的分析,雖然資料中可能包含了多個目標(biāo),但根據(jù)標(biāo)記的目標(biāo)代號,通?？梢灾蕾Y料屬于哪一組編隊(duì)目標(biāo),而編隊(duì)的目標(biāo)集通常是已知的.在目標(biāo)集確定后,標(biāo)記錯誤等價于標(biāo)記缺失.本文從極大似然估計(jì)角度,將標(biāo)記作為待估量和軌道量一同考慮,采用EM方法給出了解決方案,實(shí)現(xiàn)了觀測資料的區(qū)分與多個目標(biāo)軌道同時改進(jìn),提高了資料的利用率.

2 問題的分析

2.1 軌道確定問題

簡要回顧一下軌道確定(軌道改進(jìn))問題.令→Y為觀測量,→X為軌道狀態(tài)量,有測量方程

以及狀態(tài)轉(zhuǎn)移關(guān)系:

其中0為t0時刻的狀態(tài)量.F、G為非線性函數(shù).若t0時刻參考軌道為,上述方程可在處線性化為:

其中表示,即t0時刻真值和參考軌道的差.由上述兩式可得到軌道改進(jìn)的方程:

(5)式就是軌道改進(jìn)的線性化方程,左邊通常稱為O?C.后續(xù)所有討論都針對(5)式,在敘述中我們不再區(qū)分原始觀測量和O?C量(統(tǒng)一稱作觀測量),以及狀態(tài)量和狀態(tài)改變量(統(tǒng)一稱作狀態(tài)量).

2.2 極大似然估計(jì)

對于本文討論的多星定軌問題,令采集了屬于n個目標(biāo)的m圈觀測資料,觀測量為共m圈,表示第i圈資料;每圈又包含若干個資料其中pi是第i圈資料個數(shù).狀態(tài)量為其中表示第j個目標(biāo)的狀態(tài)量．定軌還需要知道資料和目標(biāo)的對應(yīng)關(guān)系,引入變量當(dāng)Zij=1時,表示對應(yīng)于目標(biāo).每圈必然并且僅屬于一個目標(biāo),有:

軌道確定問題本質(zhì)是最優(yōu)估計(jì)問題,考慮極大似然估計(jì),對數(shù)似然函數(shù)為：

其中p(·)表示概率密度函數(shù),軌道改進(jìn)即求:

假設(shè)

則

對于一般情形,和的對應(yīng)關(guān)系Zij是已知的,可按照j將Zij=1的圈分為1組,即按照目標(biāo)分為n組:

其中Gk={i,Zik=1}.由于各Lk之間無關(guān),最大化L等價于最大化每個Lk.假設(shè)測量誤差為正態(tài)分布:

其中是第i圈的測量誤差:

則有:

Lk最大時滿足:

等價于逐個目標(biāo)分別求解．(19)式的解也就是單目標(biāo)定軌的極大似然估計(jì)[2].

若令:

上式中···表示每個Zik值重復(fù)pi次.每個目標(biāo)的解可寫為統(tǒng)一形式:

形式上和單目標(biāo)軌道改進(jìn)一樣,只是增加了權(quán)重→Zk.

2.3 缺失變量模型與EM方法

對于本文討論的多星混合資料定軌問題,由于標(biāo)號是錯誤的,因此可看作是未知的,即缺失變量.假設(shè)的分布函數(shù)為q(),有:

多星混合資料定軌問題的似然函數(shù)為:

可分解為:

其中右邊第2項(xiàng)KL(q||p)為Kullback-Leibler散度(KL divergence),其性質(zhì)為KL(q||p)≥0,當(dāng)且僅當(dāng)時取等號.顯然L為似然函數(shù)L的下界[3]．

由于缺失變量的存在,直接求解(25)式的最大值是十分困難的,可通過EM方法求解[4].EM方法的求解思路是,由于L與無關(guān),因此L最大時,KL(q||p)=0,等價于最大化L的下界L.而

EM方法為一迭代過程,每次迭代分為兩個步驟:(1)期望步驟(E步驟,E Step):固定old,得到;(2)最大化步驟(M步驟,M Step),固定E步得到的q,在完整模型下求得new,此時由于發(fā)生了變化,因此KL(q||p)＞0,將new作為old重新回到E步驟.迭代收斂時就得到了最大化L的.

3 問題的解決

根據(jù)上述思路,本文提出的問題得以解決.利用、估計(jì),由此得到完整的觀測量,化為多個一般定軌問題.由(10)～(11)式,根據(jù)Bayes公式可得:

其中πj理解為第j個目標(biāo)的先驗(yàn)概率,即獲得的m圈資料中屬于j目標(biāo)的先驗(yàn)比例.由此可得到的數(shù)學(xué)期望.E步驟為:

將γij作為Zij代入(13)式:

可見M步驟最大化的是似然函數(shù)相對于的數(shù)學(xué)期望.因此EM方法是最大化似然函數(shù)的期望.L依然按照下標(biāo)j分組,最大化L依然等價于最大化每個Lk.將γij作為Zij代入(23)式可得到解.

多星混合資料定軌過程總結(jié)如下:(1)已知資料和精度矩陣以及觀測所屬目標(biāo)的軌道初值;(2)由(28)式根據(jù)和計(jì)算;(3)根據(jù),按照一般軌道改進(jìn)過程,以1作為權(quán)矩陣,按目標(biāo)逐一求得改進(jìn)量判斷改進(jìn)量是否收斂,不收斂返回第2步.

4 穩(wěn)健估計(jì)

上述雖然實(shí)現(xiàn)了多星混合資料定軌的完整過程,但在具體算法實(shí)現(xiàn)中還需要解決一個問題.上述過程中雖然Zij∈{0,1},但γij∈[0,1],由(23)式可見,對于某個目標(biāo)的定軌中將會使用所有的資料,分別以γij為權(quán)重,在定軌的前幾次迭代中,由于識別尚不完全準(zhǔn)確,會出現(xiàn)同一圈資料以不同權(quán)重屬于不同目標(biāo)的情況,即0＜γij＜1,意味著對于某目標(biāo)的軌道改進(jìn)中會遇到大偏差的觀測量.由于最小二乘估計(jì)是不穩(wěn)健的,這將導(dǎo)致求解結(jié)果出現(xiàn)大偏差,從而影響后續(xù)識別過程,最終導(dǎo)致求解失敗.因此在前幾次迭代中必須使用穩(wěn)健方法,來消除因識別尚不十分精確帶來的影響.

文獻(xiàn)[5-6]使用最小一乘估計(jì)用于初軌計(jì)算,指出最小一乘估計(jì)一方面具有很高的崩潰點(diǎn),另一方面無需任何參數(shù)的設(shè)置,適用性比較好.這兩個特點(diǎn)非常有利于上述問題的解決,而且識別精確后,對于本文問題可以轉(zhuǎn)為最小二乘估計(jì),從而彌補(bǔ)最小一乘估計(jì)在解唯一性以及漸近效率上的不足.文獻(xiàn)[5-6]將最小一乘估計(jì)問題轉(zhuǎn)換為線性規(guī)劃問題求解．本文討論的是軌道改進(jìn)問題,已具有初始參考解,可采用迭代重加權(quán)方法(IRLS,Iteratively Reweighted Least Square)求解[7].對于定軌問題的最小二乘估計(jì)是:

最小一乘估計(jì)是:

IRLS方法也是一個迭代過程:

即以參考軌道計(jì)算的|O?C|?1作為權(quán)重,進(jìn)行最小二乘估計(jì),迭代收斂時

由于軌道改進(jìn)本身就是迭代過程,因此無需每次迭代中都通過迭代嚴(yán)格求解最小一乘估計(jì),軌道改進(jìn)迭代過程中同時更新參考解和權(quán)重,因此應(yīng)用IRLS方法只需在軌道改進(jìn)過程中增加權(quán)重|O?C|?1即可.

5 數(shù)值模擬

采用數(shù)值模擬方法對上述方法進(jìn)行了驗(yàn)證,選用的兩個目標(biāo)在歷元時刻直角坐標(biāo)系下的位置、速度見表1.兩個目標(biāo)位于相同軌道,平近點(diǎn)角差0.6°.

表1 模擬參考軌道Table 1 The reference orbits of simulation

模擬了3天2站6圈的觀測數(shù)據(jù),不同目標(biāo)選用不同測站的觀測,模擬數(shù)據(jù)類型為赤經(jīng)和赤緯(α,δ).力模型包含90×90階地球引力場、日月第三體攝動,太陽光壓和大氣阻力攝動.得到的資料情況如表2.兩個目標(biāo)相對同一個測站觀測角距＜2°.

表2 模擬圈次情況Table 2 Description of simulated passes

為了驗(yàn)證程序的正確性,首先采用了和模擬完全一樣的力模型,觀測資料沒有加誤差.改進(jìn)參數(shù)為目標(biāo)的狀態(tài)量,求解中前3次采用最小一乘估計(jì).沒有其它先驗(yàn)信息,取π1=π2=0.5.定軌中采用的初值見表3.

表3 定軌初值Table 3 The initial values of orbit determination

經(jīng)過6次迭代后,兩個目標(biāo)的位置改進(jìn)量都＜1m,得到的對應(yīng)關(guān)系與模擬一致.最終狀態(tài)量見表4,與模擬采用的參考值一致,表明方法和程序的正確性.

表4 無誤差資料的軌道改進(jìn)結(jié)果Table 4 The results of orbit improvement with the observations without errors

為了仿真實(shí)際情況,定軌中將地球引力場在20×20階處截斷,大氣阻力中系數(shù)CD由2.2變?yōu)?.3,作為模型誤差.觀測資料中增加了5′′的隨機(jī)差,其它參數(shù)保持不變.定軌初值不變,對應(yīng)關(guān)系求解正確,最終結(jié)果見表5.

表5 有誤差資料的軌道改進(jìn)結(jié)果Table 5 The results of orbit improvement with the observations with errors

在實(shí)際工作中,當(dāng)根據(jù)標(biāo)記目標(biāo)處理數(shù)據(jù)失敗時,是無法預(yù)先得知其中是否包含了多星,可能是只包含了編隊(duì)中的另一個目標(biāo).因此模擬了一個目標(biāo)的情況,僅保留目標(biāo)2的3圈資料,采用同樣的過程處理,依然選取π1=π2=0.5,結(jié)果見表6.目標(biāo)2得到了正確的結(jié)果,最終識別結(jié)果為γ11=γ21=γ31=0,γ12=γ22=γ32=1,即所有3圈都屬于目標(biāo)2.目標(biāo)1結(jié)果并不正確,是因?yàn)樵诘?次迭代的E步驟中由于初值誤差,得到γ31=0.07,γ32=0.93.從而導(dǎo)致目標(biāo)1在M步驟中定軌錯誤,第2次迭代開始已沒有任何屬于目標(biāo)1的資料,軌道也不再修正.這個現(xiàn)象在實(shí)際計(jì)算中是容易判斷的,并不會構(gòu)成任何混淆.

表6 單目標(biāo)軌道改進(jìn)結(jié)果Table 6 The results of orbit improvement with the observations of single object

本文采用的穩(wěn)健方法也可單獨(dú)用于常規(guī)的軌道改進(jìn)中,保留了3圈目標(biāo)2和1圈目標(biāo)1的觀測資料(去除第1和第5圈資料)共142個資料,僅采用最小二乘估計(jì),無法完成軌道改進(jìn).前3次按本文加權(quán)方法采用最小一乘估計(jì)后再使用最小二乘估計(jì),并根據(jù)3倍殘差進(jìn)行資料剔除.最終正確剔除了屬于目標(biāo)1的33個資料,得到了正確的改進(jìn)結(jié)果,見表7.

表7 穩(wěn)健軌道改進(jìn)結(jié)果Table 7 The results of robust orbit improvement

最后模擬了3星混合資料的情況,新增目標(biāo)3的參考軌道和定軌初值見表8,目標(biāo)3和目標(biāo)2的觀測角距最大約為2°.將目標(biāo)3的2圈模擬資料和原先的6圈資料共同組成8圈3星混合資料,軌道改進(jìn)中取π1=π2=π3=0.33,最終求解結(jié)果正確,見表9.

表8 目標(biāo)3的模擬參數(shù)Table 8 The simulation parameters of the object 3

表9 3星軌道改進(jìn)結(jié)果Table 9 The results of orbit improvement with the three objects

6 結(jié)論與討論

本文針對光學(xué)觀測日常定軌中的多星混合資料問題,給出了有效的處理方法.不同于現(xiàn)有方法,方法將資料標(biāo)記和軌道量一同作為待估量,并融合了識別和定軌過程,識別不僅依賴初值,在軌道改進(jìn)過程中隨著軌道修正不斷動態(tài)調(diào)整.在軌道改進(jìn)中,采用IRLS方法實(shí)現(xiàn)最小一乘穩(wěn)健估計(jì),不改變原有迭代計(jì)算過程,權(quán)重選取也無需調(diào)整參數(shù),應(yīng)用有效、簡便.

本文沒有采用逐點(diǎn)識別,是考慮到現(xiàn)有跟蹤系統(tǒng)都通過航跡關(guān)聯(lián)實(shí)現(xiàn)閉環(huán)跟蹤.同一圈資料一般屬于同一個目標(biāo),按圈設(shè)置Zij更為合理.對于單點(diǎn)采集資料,例如搜索或巡天觀測方式,本文方法同樣適用,只需將每圈觀測看作只有一個采樣即可.

[1]吳連大.人造衛(wèi)星和空間碎片的軌道和探測.北京:中國科學(xué)技術(shù)出版社,2011:250-253

[2]Tapley B D,Schuts B E,Born G R.Statistical Orbit Determination.Amsterdam:Elsevier Academic Press,2004:190-194

[3]Bishop C M.Pattern Recognition and Machine Learning.New York:Springer,2006:430-443

[4]Dempster A P,Laird N M,Rubin D B.Journal of the Royal Statistical Society.Series B(Methodological),1977,39:1

[5]王歆.天文學(xué)報,2013,54:274

[6]Wang X.ChA&A,2013,37:455

[7]Scales J A,Gersztenkorn A.InvPr,1988,4:1071

Orbit Determination with Mixture Observations of Multiple Objects

WANG Xin1,2
(1 Purple Mountain Observatory,Chinese Academy of Sciences,Nanjing 210008)(2 Key Laboratory of Space Object and Debris Observation,Purple Mountain Observatory,Chinese Academy of Sciences,Nanjing 210008)

In the operational orbit determination with optical measurements of space objects,some observations of different objects are tagged as the same object.For this kind of data,the orbit improvement according to the tag is failed because of the composition of multiple objects.A method is proposed from the view of maximum likelihood,and it combines the orbit improvement and identi fication by employing the EM(Expectation Maximum)method.In the implementation of this method,a robust estimation is also given.Corresponding numerical simulations show that the method is feasible,effective,and convenient.

space vehicles,celestial mechanics,methods:statistical

P135;

A

2014-05-04收到原稿,2014-06-10收到修改稿

?國家自然科學(xué)基金項(xiàng)目(11373072)資助

?wangxin@pmo.ac.cn