杜建麗 徐 勁 楊 冬 曹志斌 馬劍波

（1 中國科學(xué)院紫金山天文臺(tái) 南京 210023）

（2 中國科學(xué)院空間目標(biāo)與碎片觀測(cè)重點(diǎn)實(shí)驗(yàn)室 南京 210023）

1 引言

空間碎片的軌道信息是空間態(tài)勢(shì)感知的重點(diǎn)關(guān)注要素, 由其形成的編目數(shù)據(jù)庫是太空資源開發(fā)利用、太空資產(chǎn)碰撞預(yù)警、太空垃圾清理等應(yīng)用的基礎(chǔ), 其重要性得到了廣泛的認(rèn)可. 近20 yr來, 全球空間碎片專用、兼用和可用監(jiān)測(cè)設(shè)備的數(shù)量大幅增長(zhǎng), 形成以美國空間監(jiān)視網(wǎng)為代表的多個(gè)監(jiān)測(cè)網(wǎng), 設(shè)備類型主要為無線電雷達(dá)、光學(xué)望遠(yuǎn)鏡和激光雷達(dá), 同一目標(biāo)被多個(gè)設(shè)備探測(cè)或跟蹤的概率大大增加. 利用多站點(diǎn)、多類型設(shè)備獲取的資料編目定軌是組網(wǎng)觀測(cè)要解決的問題之一, 由于多站點(diǎn)組網(wǎng)改善了測(cè)量幾何、增加了觀測(cè)頻率、可聯(lián)合高精度資料, 融合定軌的精度優(yōu)于單站定軌.

編目工作中的軌道計(jì)算有別于一般應(yīng)用中的軌道計(jì)算, 不僅要求軌道計(jì)算的精度足以實(shí)現(xiàn)對(duì)空間碎片運(yùn)動(dòng)態(tài)勢(shì)的精確掌控, 還要求有較高的軌道計(jì)算效率. 編目工作對(duì)軌道計(jì)算效率的最低要求是軌道計(jì)算的處理時(shí)間不能大于自動(dòng)運(yùn)行設(shè)備的數(shù)據(jù)采集時(shí)間, 以避免產(chǎn)生數(shù)據(jù)擠壓的瓶頸效應(yīng), 使編目過程難以持續(xù). 以編目目標(biāo)2萬、每天進(jìn)行一次軌道更新計(jì)算為例, 單個(gè)目標(biāo)的平均處理時(shí)間要求不大于4.32 s, 對(duì)編目定軌的計(jì)算效率要求很高. 當(dāng)前編目工作中的軌道計(jì)算廣泛采用分析方法, 例如美國采用SGP4 （Simplified General Perturbation）/SDP4 （Simplified Deep space Perturbation）分析模型, 國內(nèi)編目工作通常采用擬平均根數(shù)法等. 雖然分析方法較好地解決了編目工作的時(shí)效性問題, 但隨著觀測(cè)技術(shù)的發(fā)展, 它的模型精度并不能與當(dāng)前的觀測(cè)精度（雷達(dá)測(cè)距精度優(yōu)于20 m,光學(xué)測(cè)角精度優(yōu)于10′′）相匹配, 從而不能達(dá)到理想的定軌精度, 制約了相關(guān)應(yīng)用水平的提升, 因此數(shù)值方法在編目工作中具有潛在的應(yīng)用價(jià)值. 另一方面, 現(xiàn)代計(jì)算機(jī)運(yùn)算速度的不斷提升也使數(shù)值方法應(yīng)用于編目處理成為了可能. 數(shù)值法應(yīng)用于編目定軌, 應(yīng)轉(zhuǎn)變應(yīng)用于人造衛(wèi)星定軌時(shí)“一星一方案”的處理思路, 重點(diǎn)解決通用性、時(shí)效性差的問題, 具體的研究方向?yàn)閿U(kuò)展數(shù)據(jù)類型、提升計(jì)算速度、全軌道類型適用以及海量數(shù)據(jù)帶來的計(jì)算穩(wěn)健性要求等.

作者在以下方面做了大量有益探索, 包括調(diào)研國內(nèi)主流空間監(jiān)測(cè)設(shè)備情況（設(shè)備類型、安裝方式、資料類型等）,調(diào)研空間碎片數(shù)據(jù)處理現(xiàn)狀和改造要求（精度和速度等）, 針對(duì)數(shù)值法軌道改進(jìn)計(jì)算耗時(shí)的關(guān)鍵環(huán)節(jié)研發(fā)兼顧精度和效率的優(yōu)化算法,針對(duì)極性奇點(diǎn)和數(shù)學(xué)奇點(diǎn)進(jìn)行無奇點(diǎn)化處理, 針對(duì)多源數(shù)據(jù)先驗(yàn)精度失真進(jìn)行自適應(yīng)加權(quán)方法研究等, 以上調(diào)研結(jié)果和研究結(jié)果是形成SPODFMD（Space debris Precise Orbit Determination Fusing Multi-source Data）軟件架構(gòu)和定軌算法的基礎(chǔ).本文將首先介紹軟件的功能和特點(diǎn), 再介紹軟件集成的關(guān)鍵技術(shù), 最后利用多源數(shù)據(jù)進(jìn)行測(cè)試并根據(jù)測(cè)試結(jié)果總結(jié)軟件的性能.

2 軟件的功能和特點(diǎn)

SPODFMD軟件能滿足絕大部分主流空間監(jiān)測(cè)設(shè)備獲取數(shù)據(jù)的融合需求, 并且具備多種特點(diǎn).

2.1 軟件功能（見附錄）

2.2 軟件特點(diǎn)

多資料:軟件支持14類設(shè)備的數(shù)據(jù)輸入和任意融合, 設(shè)備類型見附錄.

高效性:軟件的效率在定軌算法和編程技巧兩方面得到了優(yōu)化. 算法方面, 一方面采用了先積分再插值的策略替代逐點(diǎn)積分, 一個(gè)軌道周期積分的基點(diǎn)數(shù)隨偏心率增大由25遞增至最多80個(gè), 實(shí)際構(gòu)造的基點(diǎn)數(shù)只需要覆蓋星歷點(diǎn)即可, 由于插值效率很高,星歷的計(jì)算時(shí)間約等于插值基點(diǎn)的積分時(shí)間,相比逐點(diǎn)積分密集星歷節(jié)省了大量時(shí)間;另一方面,計(jì)算一次積分最為復(fù)雜的項(xiàng)是地球重力及其偏導(dǎo)數(shù)的計(jì)算, 其耗時(shí)程度隨階數(shù)的增加呈平方增長(zhǎng),M′etris系數(shù)變換法可以解決這一問題, 無論是計(jì)算重力勢(shì)還是重力勢(shì)的任意階偏導(dǎo)數(shù), 都可轉(zhuǎn)換為球諧級(jí)數(shù)的求和, 球諧項(xiàng)的系數(shù)可事先通過變換一次求出, 再次節(jié)省了計(jì)算時(shí)間. 此外, 在計(jì)算大氣密度相對(duì)于地固矩坐標(biāo)的偏導(dǎo)數(shù)時(shí)采用了解析表達(dá)式,相比有限差分求導(dǎo)需要多次計(jì)算大氣密度的效率要高. 編程技巧方面, 軟件采用了高效的排序方法、搜索算法、滑動(dòng)窗口法等進(jìn)一步提高計(jì)算效率.

通用性:軟件適用于低、中、高軌各種偏心率的全軌道空間碎片. 偏心率通用性主要通過非均勻插值基點(diǎn)構(gòu)造方法和多種積分器來實(shí)現(xiàn), 非均勻插值基點(diǎn)構(gòu)造方法是針對(duì)大偏心率軌道利用均勻基點(diǎn)插值時(shí)近地點(diǎn)和遠(yuǎn)地點(diǎn)星歷精度失衡這一狀況提出的, 實(shí)際上適用于任何偏心率軌道, 可與任一單步法積分器結(jié)合使用; 軟件提供了4種積分器以及3種組合積分器以適應(yīng)軌道偏心率的變化, 其中AC （Adams-Cowell）積分器適用于偏心率小于0.1的軌道, 而RKF （Runge-Kutta-Fehlberg）、RKNF（Runge-Kutta-Nystrom-Fehlberg）和Everhart適用于任何偏心率軌道. 軌道高度的通用性主要通過定軌流程實(shí)現(xiàn), 除6個(gè)軌道根數(shù)外, 高軌目標(biāo)同時(shí)改進(jìn)與光壓相關(guān)的物理量,低軌目標(biāo)同時(shí)改進(jìn)與大氣相關(guān)的物理量, 中軌目標(biāo)同時(shí)改進(jìn)與大氣和光壓相關(guān)的兩個(gè)物理量, 通過對(duì)空間碎片軌道高度的區(qū)分實(shí)現(xiàn)了軌道高度全覆蓋.

穩(wěn)健性:軟件的穩(wěn)健性主要體現(xiàn)在無奇點(diǎn)上.軟件集成了兩種地球重力勢(shì)及其偏導(dǎo)數(shù)算法, 分別為Balmino和M′etris系數(shù)變換法, 它們的共同特點(diǎn)是都進(jìn)行了無奇點(diǎn)變換, 計(jì)算兩極上空空間點(diǎn)的引力和偏導(dǎo)數(shù)均不會(huì)引起定義和計(jì)算的奇點(diǎn).DTM94 （Drag Temperature Model）大氣密度模型也具備無奇點(diǎn)變換的條件, 經(jīng)過與地球重力勢(shì)相同的處理后也達(dá)到了消除奇點(diǎn)的效果, 進(jìn)一步提升了軟件的穩(wěn)健性.

模塊化:軟件采用模塊化設(shè)計(jì), 將功能獨(dú)立且有重復(fù)使用需求的函數(shù)全部采用模塊管理, 由此帶來了程序設(shè)計(jì)的靈活性, 例如動(dòng)力學(xué)模型共有24個(gè)相關(guān)模塊, 分別用于定軌/預(yù)報(bào)、低/中/高軌、Balmino/M′etris、1/2階積分器及以上4類配置選項(xiàng)的任意組合, 例如“定軌+低軌+M′etris+2階積分器”對(duì)應(yīng)一個(gè)單獨(dú)的模塊, 程序的修改、調(diào)用非常便利, 也便于實(shí)現(xiàn)更為復(fù)雜的功能.

靈活性:軟件的靈活性是指多資料這一特點(diǎn)在概念上的延伸, 軟件支持的14類設(shè)備是邏輯設(shè)備的概念, 邏輯設(shè)備可與實(shí)際設(shè)備一一對(duì)應(yīng), 也可與實(shí)際設(shè)備完全不同, 應(yīng)根據(jù)數(shù)據(jù)特點(diǎn)和融合需求靈活設(shè)置設(shè)備的類型. 例如一臺(tái)相控陣?yán)走_(dá), 由于測(cè)距精度遠(yuǎn)高于測(cè)角精度, 若僅使用測(cè)距資料與其他高精度測(cè)角設(shè)備獲取的數(shù)據(jù)融合, 調(diào)用軟件時(shí)僅需要將相控陣?yán)走_(dá)的設(shè)備類型設(shè)置為激光雷達(dá). 類似的場(chǎng)景還有將一臺(tái)激光測(cè)距和測(cè)角設(shè)備作為單基雷達(dá)三元素設(shè)備進(jìn)行輸入等, 采用邏輯設(shè)備的理念使得數(shù)據(jù)融合場(chǎng)景更加豐富、人性化.

智能化:智能化體現(xiàn)在軟件對(duì)物理參數(shù)改進(jìn)方面的處理, 除6個(gè)軌道根數(shù)外, 定軌時(shí)還嘗試同時(shí)改進(jìn)與目標(biāo)相關(guān)的物理參數(shù), 如果迭代發(fā)散或改進(jìn)值不合理則自動(dòng)固定其物理參數(shù)重新改進(jìn), 整個(gè)過程無需人工干預(yù).

信息全反饋: 軟件處理的所有環(huán)節(jié)都會(huì)反饋必要的信息, 通過17類反饋信息的值可追溯不合理的輸入、積分異常、堆棧異常等問題, 符合工程應(yīng)用對(duì)軟件的要求.

3 關(guān)鍵技術(shù)

SPODFMD軟件集成的關(guān)鍵技術(shù)主要有地球重力勢(shì)及其偏導(dǎo)數(shù)無奇點(diǎn)快速計(jì)算方法、無奇點(diǎn)DTM94大氣密度模型及其偏導(dǎo)數(shù)的解析表達(dá)方法、大偏心率軌道密集星歷精密快速計(jì)算方法以及穩(wěn)健的自適應(yīng)加權(quán)方法等. 這些算法雖然是為解決某一問題而提出, 但在開發(fā)過程中都兼顧了定軌對(duì)效率、精度、通用和穩(wěn)健的要求, 它們的集成就是從各個(gè)角度同時(shí)對(duì)定軌過程進(jìn)行優(yōu)化, 最終決定了SPODFMD軟件的效果.

3.1 大偏心率軌道密集星歷精密快速計(jì)算方法

編目工作的軌道計(jì)算普遍涉及密集星歷的產(chǎn)生和計(jì)算問題, 星歷點(diǎn)密集時(shí)逐點(diǎn)積分步長(zhǎng)不能得到充分伸展, 積分右函數(shù)計(jì)算次數(shù)大大增加, 嚴(yán)重降低數(shù)值方法的計(jì)算效率, 避免這一問題的有效技術(shù)手段就是運(yùn)用內(nèi)插方法. 目前產(chǎn)生密集星歷的多種內(nèi)插方法均采用了定步長(zhǎng)劃分插值基點(diǎn)的方式,這種方式適用于運(yùn)動(dòng)特征變化非常平緩的近圓軌道, 對(duì)于編目工作中經(jīng)常涉及的運(yùn)動(dòng)特征變化較為劇烈的大偏心率軌道并不適用, 易造成近地點(diǎn)附近星歷因插值基點(diǎn)稀疏導(dǎo)致插值精度低, 而遠(yuǎn)地點(diǎn)附近星歷又因插值基點(diǎn)密集而精度過高的失衡狀況.可通過增加插值基點(diǎn)的數(shù)量來改善這一失衡狀況,但相應(yīng)需增加數(shù)值積分時(shí)間,這又易降低計(jì)算效率.為兼顧大偏心率軌道整體的星歷計(jì)算精度和計(jì)算效率, 提出一種非均勻插值基點(diǎn)構(gòu)造方法[1–3], 具體為從已知星歷的初始基點(diǎn)開始遞次構(gòu)造插值基點(diǎn),待構(gòu)造的基點(diǎn)由上一個(gè)基點(diǎn)的星歷計(jì)算得出, 基點(diǎn)的星歷由數(shù)值積分得出, 基點(diǎn)構(gòu)造與數(shù)值積分計(jì)算交替進(jìn)行,直到將涵蓋所有插值點(diǎn)的基點(diǎn)構(gòu)造完畢.非均勻插值基點(diǎn)的構(gòu)造方法見（1）式,基點(diǎn)間隔是空間目標(biāo)地心距的函數(shù), 該方法可以與一般的單步積分法在星歷計(jì)算過程中協(xié)調(diào)使用, 當(dāng)與Hermite插值多項(xiàng)式結(jié)合使用時(shí)總效率最高.

其中, ?Tl為第l個(gè)插值基點(diǎn)與待構(gòu)造的第l+1個(gè)基點(diǎn)的間隔;a為軌道半長(zhǎng)徑; γ為時(shí)間變換參數(shù); β為軌道周期變換因子;Rl為利用第l個(gè)基點(diǎn)星歷計(jì)算得出的地心距;Pt為以時(shí)間變量t計(jì)算的軌道運(yùn)動(dòng)周期;N為一個(gè)軌道周期內(nèi)插值基點(diǎn)的數(shù)量. 此外,γ和N為插值基點(diǎn)構(gòu)造之前預(yù)先確定的常量, γ由數(shù)值實(shí)驗(yàn)結(jié)果確定的最佳普適數(shù)值為0.3, β由軌道偏心率e和γ的取值決定, 在基點(diǎn)構(gòu)造過程中也是一個(gè)常量.

3.2 地球重力勢(shì)及其偏導(dǎo)數(shù)無奇點(diǎn)快速計(jì)算方法

地球重力勢(shì)通常表達(dá)為球坐標(biāo)的函數(shù), 這種表達(dá)方式引入了計(jì)算奇點(diǎn), 表現(xiàn)為兩點(diǎn): （1）位于兩極上空點(diǎn)位的經(jīng)度定義不明確; （2）求解重力勢(shì)相對(duì)于地固矩坐標(biāo)偏導(dǎo)數(shù)時(shí)分母出現(xiàn)緯度余弦,在兩極上空時(shí)引起分母為零. 從物理角度考慮,地球外部重力勢(shì)連續(xù)分布不存在局部奇點(diǎn), 上述奇點(diǎn)非本性奇點(diǎn)而是數(shù)學(xué)奇點(diǎn), 通過引入適當(dāng)?shù)臄?shù)學(xué)方法可以消除. Pines[4]首先引入了地球重力勢(shì)的一種無奇點(diǎn)表達(dá)形式, 將由勒讓德多項(xiàng)式表達(dá)的地球重力勢(shì)球諧函數(shù)轉(zhuǎn)換為由亥姆霍茲多項(xiàng)式表達(dá), 并引入地心距和方向余弦取代球坐標(biāo);Balmino等[5]對(duì)Pines[4]的結(jié)果作了進(jìn)一步的整理和完善, 他們的基本思想都是利用復(fù)合求導(dǎo)法則將重力勢(shì)相對(duì)于地固坐標(biāo)的一、二階偏導(dǎo)數(shù)轉(zhuǎn)化為相對(duì)于地心距和方向余弦的偏導(dǎo)數(shù), 經(jīng)轉(zhuǎn)化后的偏導(dǎo)數(shù)不會(huì)在分母上產(chǎn)生緯度余弦的因子, 因此是無奇點(diǎn)的. Hotine[6]發(fā)現(xiàn)任何球諧級(jí)數(shù)相對(duì)于矩坐標(biāo)的偏導(dǎo)數(shù)還可以表示為一個(gè)球諧級(jí)數(shù), 于是將地球重力勢(shì)相對(duì)于地固矩坐標(biāo)的偏導(dǎo)數(shù)計(jì)算轉(zhuǎn)變?yōu)橐环N系數(shù)變換, Bettadpur[7]做了進(jìn)一步的研究和改進(jìn). M′etris等[8]吸收了Bettadpur[7]和Pines[4]方法的結(jié)論, 對(duì)重力勢(shì)展開式中所有單個(gè)立體球諧項(xiàng)的Cunningham[9]結(jié)果進(jìn)行求和, 并對(duì)和式中的立體諧函數(shù)按度和階進(jìn)行歸并和整理, 得到了重力勢(shì)相對(duì)于地固矩坐標(biāo)偏導(dǎo)數(shù)的無奇點(diǎn)緊致公式, 偏導(dǎo)公式在形式上完全類似于原重力勢(shì)表達(dá)式, 其諧系數(shù)和部分變量與原系數(shù)和原變量之間存在變換關(guān)系, 需要變換的量為:μ →μ′,N1→N′1,N2→N′2,以及L →L′的變換,上述變量的具體定義及變換方法見下文（2）–（3）式.

SPODFMD集成了Balmino和M′etris兩種計(jì)算方法. 軟件集成的Balmino方法對(duì)原公式中分開處理的地球中心引力項(xiàng)、帶諧項(xiàng)和田諧項(xiàng)作了統(tǒng)一處理, 表達(dá)簡(jiǎn)潔且便于計(jì)算機(jī)程式設(shè)計(jì), 其不足之處在于需要遞推計(jì)算亥姆霍茲多項(xiàng)式及其一、二階偏導(dǎo)數(shù), 影響計(jì)算效率. M′etris方法實(shí)際為Hotine算法的一種推廣, 是迄今為止能夠?qū)崿F(xiàn)地球重力勢(shì)二階以上高階偏導(dǎo)數(shù)的唯一方法, 無論是計(jì)算重力勢(shì)還是重力勢(shì)的任意階偏導(dǎo)數(shù), 其計(jì)算形式是統(tǒng)一的, 均為球諧級(jí)數(shù)的求和, 只是球諧項(xiàng)的系數(shù)不同. 鑒于這些系數(shù)可事先通過對(duì)應(yīng)的變換求得, 大大節(jié)省了計(jì)算時(shí)間, 當(dāng)不考慮地球重力勢(shì)諧系數(shù)的變化時(shí)（如忽略潮汐的影響）, 至少可節(jié)省30%的計(jì)算時(shí)間[10].

設(shè)U為重力勢(shì)或其對(duì)地固矩坐標(biāo)的任意階偏導(dǎo)數(shù), 且有形式:

原形式中,μ為地心引力常數(shù);n和m分別為球諧級(jí)數(shù)展開式所截取的階和度,N1和N2分別為n的下限和上限;L為z坐標(biāo)偏導(dǎo)的階次, 原形式中其值為零;Cnm和Snm為歸一化的球諧系數(shù),Hnm為歸一化的赫姆霍茲多項(xiàng)式, 且有γm= cosm ?cosmλ,σm= cosm ?sinmλ,λ和?分別為目標(biāo)點(diǎn)的經(jīng)度和緯度. 則U相對(duì)于矩坐標(biāo)（x,y,z）求任意階（i,j,k）的偏導(dǎo)數(shù)公式為:

3.3 自適應(yīng)加權(quán)方法

具有不同測(cè)量精度的多類型資料融合時(shí)需重點(diǎn)考慮權(quán)的分配問題, 常見的做法是根據(jù)資料的先驗(yàn)精度加權(quán), 先驗(yàn)精度失真易造成軌道改進(jìn)質(zhì)量的下降. 考慮到精度失真是一種常見現(xiàn)象, 提出了一種穩(wěn)健的自適應(yīng)加權(quán)方法, 自適應(yīng)權(quán)是在固定權(quán)軌道改進(jìn)收斂的基礎(chǔ)上利用資料殘差重新分配權(quán)重直至收斂的過程, 其特點(diǎn)在于: （1）自適應(yīng)加權(quán)定軌利用了固定權(quán)軌道改進(jìn)收斂的結(jié)果作為初根數(shù), 因此啟動(dòng)前會(huì)默認(rèn)先執(zhí)行固定權(quán)定軌, 這樣處理可避免直接使用自適應(yīng)權(quán)時(shí)（自適應(yīng)權(quán)由初軌預(yù)報(bào)值與資料值的殘差確定）過度依賴初根數(shù)的精度, 這種依賴性使得軌道改進(jìn)偏離真值的概率大大增加, 甚至出現(xiàn)發(fā)散現(xiàn)象;（2）自適應(yīng)權(quán)定軌過程中不再剔除觀測(cè)資料, 利用固定權(quán)軌道改進(jìn)最后一次迭代剔除野值后的全部資料參與自適應(yīng)定軌, 這樣處理既對(duì)精度失真資料的權(quán)進(jìn)行了一定程度的糾偏, 又保證了自適應(yīng)迭代的穩(wěn)健性, 通常自適應(yīng)權(quán)定軌僅需要少量的迭代即可收斂, 不會(huì)顯著增加計(jì)算時(shí)間.

3.4 無奇點(diǎn)DTM94大氣密度模型及其偏導(dǎo)數(shù)的解析表達(dá)方法

在低軌目標(biāo)編目定軌中需要考慮大氣阻力攝動(dòng), 對(duì)于軌道特別低或者面質(zhì)比特別大的低軌目標(biāo), 大氣阻力攝動(dòng)已與地球主要帶諧項(xiàng)J2項(xiàng)攝動(dòng)相當(dāng), 因此為了保證定軌迭代過程的收斂, 還必須計(jì)算大氣阻力相對(duì)于目標(biāo)位置和速度矢量的偏導(dǎo)數(shù). SPODFMD采用了DTM94熱大氣模型[11], 其為三維模型具有精度高的優(yōu)點(diǎn), 但由于也采用球坐標(biāo)作為基本變量, 因此模型本身及其相對(duì)于地固矩坐標(biāo)的偏導(dǎo)數(shù)也存在極性奇點(diǎn); 另一方面, 目前除Harris-Priester大氣模型已有成熟的偏導(dǎo)數(shù)解析表達(dá)式外, 其他大氣模型均因過于復(fù)雜難以求導(dǎo),利用差分計(jì)算偏導(dǎo)數(shù)是常用的替代方法, 然而當(dāng)模型不連續(xù)或不可微時(shí), 差分計(jì)算結(jié)果失真.

鑒于DTM94原模型中存在勒讓德多項(xiàng)式, 同樣借鑒Pines[4]消除地球重力勢(shì)球諧表達(dá)式奇點(diǎn)所用的方法, 用目標(biāo)地心距以及方向余弦取代球坐標(biāo)作為DTM94大氣模型的基本變量, 變換后的大氣密度模型是無奇點(diǎn)的, 且使得尋求大氣密度相對(duì)于空間目標(biāo)地固矩坐標(biāo)的解析表達(dá)成為可能.借鑒Pines[4]提出的復(fù)合求導(dǎo)法則, 提出了大氣密度ρ相對(duì)于地固矩坐標(biāo)的表達(dá)公式[12]:

其中,j= 1,2,3,ρ=ρ（x,y,z）為x、y、z的連續(xù)可微函數(shù), 且有X1=x,X2=y,X3=z,r、θ1、θ2、θ3分別為目標(biāo)的地心距和地固坐標(biāo)系下X、Y、Z軸上的方向余弦.相比差分求取偏導(dǎo)數(shù)的方法,解析表達(dá)式具有準(zhǔn)確、穩(wěn)健和高效的優(yōu)勢(shì).

4 測(cè)試與分析

多源觀測(cè)數(shù)據(jù)有兩層含義, 一是指來源于多個(gè)測(cè)站, 二是指多種數(shù)據(jù)類型, 兩者滿足其一即可稱為融合定軌, 以下將通過算例來測(cè)試軟件的性能,主要關(guān)注指標(biāo)為定軌精度和定軌時(shí)間.

4.1 單站資料定軌

單站資料指一臺(tái)“邏輯設(shè)備”同時(shí)獲取的測(cè)軌資料,例如雷達(dá)可同時(shí)獲取測(cè)角和測(cè)距三元素資料,有的雷達(dá)還能獲取包括測(cè)速在內(nèi)的四元素資料. 同時(shí)獲取的三元素或者四元素都是單站資料, 大部分單站資料都是由單臺(tái)實(shí)際設(shè)備獲取的, 也有部分單站資料涉及到收發(fā)分置的兩臺(tái)實(shí)際設(shè)備, 例如雙基雷達(dá)的距離和資料以及無線電時(shí)差和頻差資料, 此時(shí)一臺(tái)“邏輯設(shè)備”由兩臺(tái)實(shí)際設(shè)備構(gòu)成. 包含同類型數(shù)據(jù)的單站資料, 例如激光測(cè)距、光電測(cè)角這類資料由于僅包含一種元素或者精度相同的兩種元素, 定軌時(shí)等權(quán)分配權(quán)重即可, 而三元素和四元素資料包含了不同精度的測(cè)角、測(cè)距和測(cè)速3類數(shù)據(jù),即便來源于同一設(shè)備, 也涉及不同數(shù)據(jù)類型的融合問題, 空間碎片多源觀測(cè)數(shù)據(jù)融合定軌軟件采用馬爾可夫估計(jì)理論, 利用數(shù)據(jù)的先驗(yàn)精度確定各元素的權(quán)重, 解算結(jié)果穩(wěn)定準(zhǔn)確. 以下通過例子來說明單站資料的定軌能力, 各算例的數(shù)據(jù)類型、時(shí)間跨度和數(shù)據(jù)量等關(guān)鍵信息列于表1, 利用這些資料的定軌性能主要指標(biāo)列于表2, 采用的主要?jiǎng)恿W(xué)模型和公共配置參數(shù)見表3.

表1 單設(shè)備（SE）獲取的數(shù)據(jù)資料Table 1 Data from Single Equipment （SE）

表2 SE定軌算例性能Table 2 Performance of SE orbit determination cases

表3 SE定軌的力學(xué)模型及公共參數(shù)Table 3 Force models and public parameters of SE orbit determination

算例1–5是較為典型的雷達(dá)獲取空間碎片的數(shù)據(jù)條件, 數(shù)據(jù)的時(shí)間跨度約2–5 d, 每天能探測(cè)1–2次, 整體來說符合空間碎片稀疏性的特點(diǎn). 稀疏數(shù)據(jù)條件下軌道改進(jìn)收斂通常需要較多的迭代次數(shù),由表2的統(tǒng)計(jì)結(jié)果可知以上算例在SPODFMD軟件中的迭代次數(shù)為3–8次, AC、Everhart和RKF積分器的平均處理時(shí)間分別為0.91 s、1.62 s和1.78 s,積分器的計(jì)算效率以AC 10最佳, Everhart 11和RKF 8 （7）遞次降低. AC積分器的高效眾所周知, 然而可變步長(zhǎng)的單步法積分器Everhart和RKF的計(jì)算效率接近AC積分器也是因?yàn)椴捎昧讼确e分再插值的策略, 而星歷精度的保障則主要?dú)w功于關(guān)鍵技術(shù)——非均勻插值基點(diǎn)構(gòu)造方法的應(yīng)用. 該方法可與任何一種單步法積分器結(jié)合使用, 達(dá)到兼顧各種偏心率軌道星歷精度和計(jì)算效率的目的, 這一效果不論是定步長(zhǎng)積分器AC還是變步長(zhǎng)積分器逐點(diǎn)積分均不能實(shí)現(xiàn).

數(shù)據(jù)條件是影響定軌精度的主要因素, 利用同一時(shí)刻獲取多維資料中的不同元素及其組合定軌, 單站算例1–3僅利用了雷達(dá)測(cè)距元素的軌道改進(jìn)收斂且定軌精度達(dá)到了百米甚至以下量級(jí), 僅利用了雷達(dá)測(cè)角元素的軌道改進(jìn)精度在公里級(jí)別, 同時(shí)利用測(cè)角和測(cè)距三元素的軌道改進(jìn)精度不低于僅利用高精度測(cè)距元素的定軌精度;單站算例4–5針對(duì)高軌目標(biāo)的雷達(dá)觀測(cè)增加了測(cè)速元素, 僅利用測(cè)角或測(cè)距元素的定軌誤差分別為25 km和18 km,增加高精度測(cè)速元素的融合定軌誤差分別降至7 km和5 km.以上算例體現(xiàn)了高精度觀測(cè)資料在單站軌道改進(jìn)中發(fā)揮了決定性作用, 在各元素先驗(yàn)精度較為準(zhǔn)確的情況下, 融合多維資料的軌道改進(jìn)效果一般不遜于單維資料.

4.2 多源數(shù)據(jù)定軌

多源數(shù)據(jù)融合要求同一時(shí)間段內(nèi)能獲取到來自于不同設(shè)備的觀測(cè)數(shù)據(jù), 當(dāng)組網(wǎng)觀測(cè)或者協(xié)同觀測(cè)時(shí)會(huì)產(chǎn)生多源數(shù)據(jù)融合定軌的需求. 以下用典型多站場(chǎng)景來測(cè)試軟件的性能和融合效果, 表4列出了7個(gè)融合場(chǎng)景的數(shù)據(jù)特點(diǎn), 表5為對(duì)應(yīng)場(chǎng)景的定軌性能關(guān)鍵指標(biāo), 表6為主要的動(dòng)力學(xué)模型及公共配置參數(shù).

表4 多設(shè)備（ME）獲取的數(shù)據(jù)資料Table 4 Data from Multi Equipment （ME）

表5 ME定軌算例性能Table 5 Performance of ME orbit determination cases

表6 ME定軌的力學(xué)模型及公共參數(shù)Table 6 Force models and public parameters of ME orbit determination

表7 SPODFMD軟件功能簡(jiǎn)述Table 7 Brief introduction of the software SPODFMD’s function

融合算例1–7展示了多源數(shù)據(jù)融合定軌的效果, 表5顯示高軌目標(biāo)37210利用1臺(tái)天基光學(xué)和3臺(tái)地基光學(xué)融合的精度達(dá)到了433 m, 相比利用單站雷達(dá)四元素資料7 km的定軌精度提高了一個(gè)量級(jí)以上; 高軌目標(biāo)38091利用1臺(tái)天基和1臺(tái)地基光學(xué)的融合精度優(yōu)于5 km; 低軌目標(biāo)16908、7646和46469、1328和7647的融合精度均優(yōu)于50 m.雖然依據(jù)以上少量算例不能直接得出多站融合精度高于單站的結(jié)論, 單從數(shù)值上分析, 融合的精度普遍高于單站, 這個(gè)提升主要由數(shù)據(jù)幾何條件改善引起, 數(shù)據(jù)質(zhì)量和數(shù)量是弱相關(guān)因素. 融合算例1–7利用AC積分器的平均計(jì)算時(shí)間為1.34 s, 雖然多站融合定軌處理的數(shù)據(jù)量相比單站有所增加, 但更加豐富的數(shù)據(jù)使得迭代次數(shù)減少, 總體而言計(jì)算時(shí)間僅有少許增加, 仍處于秒級(jí)水平.

4.3 動(dòng)力學(xué)模型測(cè)試

SPODFMD軟件支持70階以內(nèi)地球重力勢(shì)及其偏導(dǎo)數(shù)的計(jì)算, 并可擴(kuò)充附加潮汐、地球反照輻射壓、行星引力等攝動(dòng)力. 針對(duì)低軌目標(biāo)7647和高軌目標(biāo)37210, 圖1和圖2分別展示它們利用簡(jiǎn)約動(dòng)力學(xué)模型和精密動(dòng)力學(xué)模型進(jìn)行軌道確定并預(yù)報(bào)軌道的位置誤差（為表述方便, 由簡(jiǎn)約動(dòng)力學(xué)模型定軌及預(yù)報(bào)的軌道稱為F軌道, 由精密動(dòng)力學(xué)模型定軌及預(yù)報(bào)的軌道稱為P軌道）, 在定軌期內(nèi)7647目標(biāo)的F軌道與P軌道基本重合, 預(yù)報(bào)至2 d兩者的偏差約50 m, 預(yù)報(bào)至6 d兩者的偏差約230 m; 高軌目標(biāo)37210整個(gè)定軌及預(yù)報(bào)期內(nèi)F軌道與P軌道幾乎重合. 在效率方面, 兩個(gè)例子均顯示利用了精密動(dòng)力學(xué)模型的軌道改進(jìn)時(shí)間相比利用簡(jiǎn)約動(dòng)力學(xué)模型增加了約30%的計(jì)算時(shí)間. 已知編目目標(biāo)的軌道更新頻率約1–2 d 1次, 短期內(nèi)簡(jiǎn)約動(dòng)力學(xué)模型的定軌和預(yù)報(bào)精度與精密模型的差距不顯著, 因此不推薦犧牲計(jì)算效率的精密動(dòng)力學(xué)模型用于空間碎片編目定軌. 基于這一認(rèn)知, SPODFMD軟件在測(cè)試算例時(shí), 均采用了簡(jiǎn)約的動(dòng)力學(xué)模型, 配置方法如下: 低軌目標(biāo)取地球引力場(chǎng)20×20階, 中高軌目標(biāo)取8×8階, 地球重力勢(shì)偏導(dǎo)數(shù)計(jì)算取4階, 章動(dòng)模型取4項(xiàng), 考慮日月引力, 低軌目標(biāo)考慮大氣阻力, 高軌目標(biāo)考慮太陽光壓, 采用以上簡(jiǎn)約模型即可取得表2和表5中的定軌效果. 圖1–2用于測(cè)試的精密動(dòng)力學(xué)模型的配置如下: 低軌目標(biāo)取地球引力場(chǎng)20×20階, 中高軌目標(biāo)取8×8階, 地球重力勢(shì)偏導(dǎo)數(shù)計(jì)算取4階, 章動(dòng)模型取108項(xiàng), 考慮日月、行星引力和地球反照輻射壓, 低軌目標(biāo)考慮大氣阻力和潮汐（海潮、固體潮、大氣潮）, 高軌目標(biāo)考慮太陽光壓.

圖1 低軌目標(biāo)分別利用簡(jiǎn)約動(dòng)力學(xué)模型和精密動(dòng)力學(xué)模型定軌并預(yù)報(bào)軌道的位置誤差Fig.1 Orbital position errors of a LEO object determined and predicted by using a simplified dynamic model and a precise dynamic model

圖2 高軌目標(biāo)分別利用簡(jiǎn)約動(dòng)力學(xué)模型和精密動(dòng)力學(xué)模型定軌并預(yù)報(bào)軌道的位置誤差Fig.2 Orbital position errors of a GEO object determined and predicted by using a simplified dynamic model and a precise dynamic model

5 總結(jié)和展望

SPODFMD軟件處理上述空間碎片典型多源融合算例的運(yùn)算環(huán)境為:Windows 7; Inter Core i5-2400 CPU@3.10GHZ, 若采用主流i7及以上處理器,運(yùn)算速度還有很大提升空間, 軟件的定軌效率可以確切地說達(dá)到了秒級(jí)的水平, 在單臺(tái)計(jì)算機(jī)上運(yùn)行一天可執(zhí)行約8萬目標(biāo)的編目定軌, 并行計(jì)算預(yù)計(jì)可管理數(shù)十萬空間碎片, 在很長(zhǎng)時(shí)間內(nèi)軌道計(jì)算能力將不會(huì)成為空間監(jiān)測(cè)組網(wǎng)的瓶頸. SPODFMD軟件還是一款純數(shù)值法定軌軟件, 數(shù)值法的高精度特點(diǎn)已得到了國際國內(nèi)同行的認(rèn)可. 軟件目前可針對(duì)14種主流設(shè)備獲取的數(shù)據(jù)進(jìn)行融合, 但尚不能涵蓋所有的空間碎片觀測(cè)設(shè)備, 下一步擬增加無線電測(cè)角測(cè)速、無線電時(shí)差頻差設(shè)備類型以更好地適應(yīng)多源數(shù)據(jù)融合的需求.