本發(fā)明涉及航天器相對運(yùn)動狀態(tài)偏差傳播分析技術(shù)，具體涉及一種解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法。

背景技術(shù)：

航天器軌道偏差傳播分析在空間態(tài)勢感知各功能任務(wù)(如目標(biāo)追蹤與數(shù)據(jù)關(guān)聯(lián)、碰撞預(yù)警、傳感器導(dǎo)引、機(jī)動檢測等)中具有重要應(yīng)用。隨著空間目標(biāo)不斷增加，而觀測設(shè)備有限，大多數(shù)空間目標(biāo)不能連續(xù)跟蹤，因此需要對其軌道及偏差進(jìn)行長時間預(yù)報，以引導(dǎo)下一次跟蹤或?qū)赡艿呐鲎策M(jìn)行預(yù)警。

近年來，衛(wèi)星編隊飛行在理論研究和工程實(shí)踐中都受到廣泛關(guān)注，因為通過多個航天器編隊飛行，可以低成本、高效率的完成很多航天任務(wù)，如構(gòu)成合成孔徑雷達(dá)、太空高精度干涉測量、對地觀測等。編隊衛(wèi)星的構(gòu)型設(shè)計及保持需要知道衛(wèi)星的絕對運(yùn)動狀態(tài)及衛(wèi)星間的相對運(yùn)動狀態(tài)，由于導(dǎo)航設(shè)備誤差、星間通信鏈路時延、及觀測盲區(qū)等因素影響，使獲得的編隊衛(wèi)星的相對運(yùn)動狀態(tài)存在不確定性。由于編隊衛(wèi)星間距離相對較近，如果不對存在于相對運(yùn)動狀態(tài)中的偏差進(jìn)行精確監(jiān)控，極有可能引起編隊衛(wèi)星相互碰撞。因此，偏差傳播分析結(jié)果作為碰撞概率計算的重要依據(jù)，對編隊、集群飛行衛(wèi)星的碰撞預(yù)警及規(guī)避任務(wù)尤為重要。

軌道動力學(xué)中經(jīng)典的偏差演化分析方法有線性協(xié)方差分析方法和montecarlo仿真方法。線性協(xié)方差分析是通過將動力系統(tǒng)線性化的方式來進(jìn)行偏差傳播，該方法使用簡單、計算量小，但是對強(qiáng)非線性系統(tǒng)或長時間偏差傳播問題計算誤差較大。montecarlo仿真雖然能通過統(tǒng)計的方法得到軌道偏差的高精度分布屬性，但需要進(jìn)行大量抽樣仿真，計算量大。近年來，不少學(xué)者提出了一些高精度半解析或數(shù)值的偏差演化分析方法，有效兼顧了偏差演化分析中的計算精度與計算效率。然而，這些方法多是針對衛(wèi)星絕對狀態(tài)偏差的傳播，且需要進(jìn)行一定量的軌道積分、依然有較大的計算量，對衛(wèi)星編隊飛行任務(wù)不太適用。

特別地，對編隊飛行衛(wèi)星的偏差分析與碰撞預(yù)警多是在軌進(jìn)行的。由于星載計算機(jī)的計算能力相對較弱，因此需要偏差傳播分析方法能解析計算，以便星載計算機(jī)能快速、準(zhǔn)確地進(jìn)行碰撞預(yù)警。因此，如何基于解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析，以高精度地對衛(wèi)星相對運(yùn)動偏差進(jìn)行長時間預(yù)報，已經(jīng)成為一項亟待解決的關(guān)鍵技術(shù)問題。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明要解決的技術(shù)問題：針對現(xiàn)有技術(shù)的上述問題，提供一種基于考慮j2攝動的非線性相對運(yùn)動方程，能夠用于兩相距較遠(yuǎn)衛(wèi)星的相對運(yùn)動狀態(tài)偏差的長時間、高精度解析預(yù)報，獲得的偏差信息可用于編隊衛(wèi)星碰撞概率計算及碰撞預(yù)警，設(shè)計方法正確合理、對實(shí)際工程任務(wù)的適用性好的解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法。

為了解決上述技術(shù)問題，本發(fā)明采用的技術(shù)方案為：

一種解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法，實(shí)施步驟包括：

1)為編隊飛行的兩衛(wèi)星指定主衛(wèi)星與從衛(wèi)星，并輸入初始時刻的參考衛(wèi)星絕對軌道狀態(tài)、兩衛(wèi)星標(biāo)稱相對運(yùn)動狀態(tài)、輸入初始相對運(yùn)動狀態(tài)偏差的概率密度函數(shù)；

2)基于初始時刻的參考衛(wèi)星絕對軌道狀態(tài)、兩衛(wèi)星標(biāo)稱相對運(yùn)動狀態(tài)，根據(jù)考慮j2攝動的非線性相對運(yùn)動方程計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的一階及二階狀態(tài)轉(zhuǎn)移張量；

3)基于衛(wèi)星相對運(yùn)動狀態(tài)的一階及二階狀態(tài)轉(zhuǎn)移張量，計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階及二階狀態(tài)轉(zhuǎn)移張量；

4)判斷初始及衛(wèi)星偏差是否均為高斯分布，如果均為高斯分布則跳轉(zhuǎn)執(zhí)行步驟5)，否則，跳轉(zhuǎn)執(zhí)行步驟6)；

5)采用協(xié)方差分析方法，獲得衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值及協(xié)方差矩陣的解析非線性傳播方程，跳轉(zhuǎn)執(zhí)行步驟7)；

6)采用高斯和模型，獲得衛(wèi)星非線性相對運(yùn)動狀態(tài)偏差的概率密度函數(shù)的解析傳播方程，跳轉(zhuǎn)執(zhí)行步驟7)；

7)輸出衛(wèi)星時刻兩衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值、協(xié)方差矩陣及概率密度函數(shù)。

優(yōu)選地，步驟1)中輸入初始時刻的參考衛(wèi)星絕對軌道狀態(tài)具體是指輸入主衛(wèi)星初始時刻t0的絕對軌道狀態(tài)e0＝[a0,e0,i0,ω0,ω0,f0]，其中，e0為主衛(wèi)星初始時刻t0的初始軌道狀態(tài)，a0為主衛(wèi)星初始時刻t0的半長軸，e0為主衛(wèi)星初始時刻t0的偏心率，i0為主衛(wèi)星初始時刻t0的軌道傾角，ω0為主衛(wèi)星初始時刻t0的升交點(diǎn)赤經(jīng)，ω0為主衛(wèi)星初始時刻t0的近拱點(diǎn)角距，f0為主衛(wèi)星初始時刻t0的真近點(diǎn)角；且輸入初始時刻的兩衛(wèi)星標(biāo)稱相對運(yùn)動狀態(tài)為δx(t0)，δx(t0)在主衛(wèi)星當(dāng)?shù)剀壍雷鴺?biāo)系中表示，該坐標(biāo)系原點(diǎn)為主衛(wèi)星質(zhì)心，x軸沿主衛(wèi)星地心矢徑方向，z軸沿軌道面法向，y軸與x、z軸構(gòu)成右手坐標(biāo)系；輸入初始相對運(yùn)動狀態(tài)偏差為dx(t0)及其概率密度函數(shù)為p(t0,dx0)，若初始偏差為高斯分布，則其概率密度函數(shù)為pg(dx0；m0,c0)，m0為初始相對狀態(tài)偏差的均值矩陣，c0為初始相對狀態(tài)偏差的協(xié)方差矩陣；

優(yōu)選地，步驟2)中考慮j2攝動的非線性相對運(yùn)動方程的函數(shù)表達(dá)式如式(1)所示；

式(1)中，δx(tf)表示時刻tf從衛(wèi)星相對主衛(wèi)星的預(yù)報相對運(yùn)動狀態(tài)，φ(tf,t0)是從t0時刻到tf時刻的狀態(tài)轉(zhuǎn)移矩陣，δx(t0)表示時刻t0從衛(wèi)星相對主衛(wèi)星的預(yù)報相對運(yùn)動狀態(tài)，ψ(tf,t0)為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量，t(tf)為以緯度幅角θf為自變量的無量綱化坐標(biāo)到以時間tf為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，為從θ0緯度幅角到θf緯度幅角的狀態(tài)轉(zhuǎn)移矩陣，t(t0)為以緯度幅角θ0為自變量的無量綱化坐標(biāo)到以時間t0為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，為從θ0緯度幅角到θf緯度幅角的二階狀態(tài)轉(zhuǎn)移張量，t^-1(t0)為轉(zhuǎn)換矩陣t(t0)的逆矩陣，為kronecker張量積；

優(yōu)選地，步驟2)中計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的一階狀態(tài)轉(zhuǎn)移張量如式(2)所示，計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的二階狀態(tài)轉(zhuǎn)移張量如式(3)所示；

式(2)和式(3)中，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的從θ0緯度幅角到θf緯度幅角的一階狀態(tài)轉(zhuǎn)移張量，表示θf＝θ0時的中間變量表示的逆矩陣，表示θf＝θ0時的中間變量為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的從θ0緯度幅角到θf緯度幅角的二階狀態(tài)轉(zhuǎn)移張量，為kronecker張量積，和均為公式推導(dǎo)使用的中間變量，且中間變量和的計算函數(shù)式如式(4)所示；

式(4)中，a、b為公式推導(dǎo)使用的中間變量，表示t0時刻平均軌道要素偏差到tf時刻平均軌道要素偏差的一階傳遞矩陣，表示t0時刻軌道要素偏差的無量綱化矩陣，表示tf時刻軌道要素偏差的無量綱化矩陣，表示t0時刻平均軌道要素偏差到tf時刻平均軌道要素偏差的二階傳遞張量，q(θf)表示tf時刻無量綱軌道要素偏差到相對狀態(tài)的二階轉(zhuǎn)換張量，t(θf)表示tf時刻以緯度幅角為自變量的無量綱化坐標(biāo)到以時間為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，t^-1(θf)表示t(θf)的逆矩陣，π表示相對狀態(tài)分量次序變換矩陣，σ(θf)表示tf時刻軌道要素偏差到相對狀態(tài)的轉(zhuǎn)換矩陣，表示平均軌道要素偏差到吻切軌道要素偏差的轉(zhuǎn)換矩陣，為kronecker張量積。

優(yōu)選地，步驟3)計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階及二階狀態(tài)轉(zhuǎn)移張量的函數(shù)表達(dá)式如式(5)所示；

式(5)中，dxf為衛(wèi)星時刻的相對運(yùn)動狀態(tài)偏差，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階狀態(tài)轉(zhuǎn)移張量，dx0為初始的相對運(yùn)動狀態(tài)偏差，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的二階狀態(tài)轉(zhuǎn)移張量，為的第i行第j列元素，為的第i行第j列元素，為的第i行k列j維元素，是從t0時刻到tf時刻的狀態(tài)轉(zhuǎn)移矩陣，φ(tf,t0)是從t0時刻到tf時刻的狀態(tài)轉(zhuǎn)移矩陣，ψ(tf,t0)為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量，為t0時刻的標(biāo)稱相對運(yùn)動狀態(tài)δx0的第k個分量，為kronecker張量積。

優(yōu)選地，步驟5)中獲得的衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值及協(xié)方差矩陣的解析非線性傳播方程的函數(shù)表達(dá)式如式(6)所示；

式(6)中，mⁱ(tf)為衛(wèi)星相對狀態(tài)偏差的均值矩陣m(tf)的第i行元素，m^j(tf)為衛(wèi)星相對狀態(tài)偏差的均值矩陣m(tf)的第j行元素，c^ij(tf)為衛(wèi)星相對狀態(tài)偏差的矩陣，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階狀態(tài)轉(zhuǎn)移張量的第i行第a列元素，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階狀態(tài)轉(zhuǎn)移張量的第j行第b列元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第i行a列b維元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第j行b列c維元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第i行b列c維元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第j行c列d維元素，分別為初始相對狀態(tài)偏差dx0的前四階矩。

優(yōu)選地，初始相對狀態(tài)偏差dx0的前四階矩的計算函數(shù)表達(dá)式如式(7)所示；

式(7)中，分別為初始相對狀態(tài)偏差dx0的前四階矩，分別為初始相對狀態(tài)偏差的均值矩陣m0的a,b,c,d分量元素，分別為初始相對狀態(tài)偏差的協(xié)方差矩陣c0的ab分量元素、bc分量元素、cd分量元素、ad分量元素。

優(yōu)選地，步驟6)的詳細(xì)步驟包括：

6.1)用n個高斯分布概率密度函數(shù)去逼近初始相對狀態(tài)偏差的概率密度函數(shù)p(t0,dx0)獲得每一個子高斯分布的權(quán)重、均值和協(xié)方差矩陣；

6.2)采用衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值及協(xié)方差矩陣的解析非線性傳播方程，將每個子高斯分布的均值mi(t0)及協(xié)方差矩陣ci(t0)預(yù)報到衛(wèi)星時刻tf，預(yù)報過程中子高斯分布的權(quán)重保持不變，獲得衛(wèi)星時刻每個子高斯的均值mi(tf)及協(xié)方差矩陣ci(tf)，其中i＝1,2,…,n，n為高斯分布概率密度函數(shù)的數(shù)量；

6.3)針對每一個第i個子高斯分布，其中i＝1,2,…,n，n為高斯分布概率密度函數(shù)的數(shù)量，根據(jù)其權(quán)重系數(shù)ωi、均值mi(tf)及協(xié)方差矩陣ci(tf)，根據(jù)式(8)和式(9)所示獲得衛(wèi)星非線性相對運(yùn)動狀態(tài)偏差的概率密度函數(shù)的解析傳播方程計算衛(wèi)星相對狀態(tài)偏差的均值及協(xié)方差矩陣以及衛(wèi)星偏差的概率密度函數(shù)p(tf,dxf)；

式(8)和(9)中，m(tf),c(tf)為衛(wèi)星相對狀態(tài)偏差dxf的均值及協(xié)方差矩陣，ωi表示第i個子高斯分布的權(quán)重系數(shù)，n為高斯分布概率密度函數(shù)的數(shù)量，ω為所有高斯分布的權(quán)重系數(shù)之和，mi(tf)為tf時刻第i個子高斯分布的均值，ci(tf)為tf時刻第i個子高斯分布的協(xié)方差矩陣，m^t(tf)為mi(tf)的轉(zhuǎn)置矩陣，p(tf,dxf)為衛(wèi)星偏差的概率密度函數(shù)，pg(dxf；mi(tf),ci(tf))表示tf時刻第i個子高斯分布的概率密度函數(shù)。

優(yōu)選地，步驟6.1)中逼近初始相對狀態(tài)偏差的概率密度函數(shù)p(t0,dx0)采用的逼近公式如式(10)所示；

式(10)中，p(t0,dx0)為初始相對狀態(tài)偏差的概率密度函數(shù)，dx0為初始相對狀態(tài)偏差，pg(dx0；mi,ci)表示第i個子高斯分布的概率密度函數(shù)，mi和ci分別為相應(yīng)的均值和協(xié)方差矩陣，ωi表示第i個子高斯分布的權(quán)重系數(shù)。

優(yōu)選地，第i個子高斯分布的權(quán)重系數(shù)ωi的計算函數(shù)表達(dá)式如式(11)所示；

式(11)中，ωi表示第i個子高斯分布的權(quán)重系數(shù)，m0、c0為初始相對狀態(tài)偏差的均值及協(xié)方差矩陣，及為將一維標(biāo)準(zhǔn)正態(tài)分布pg(x；0,1)等方差地分割為n個子正態(tài)分布的權(quán)重、均值及標(biāo)準(zhǔn)差，diag{a}表示以a向量為對角線元素的矩陣，v和λ分別為對c0進(jìn)行特征值分解所得的特征向量矩陣與特征值矩陣，λ1,λ2,...,λn為c0的特征值，λj與υj分別為沿選定的j方向的特征值與特征向量，即λ的第j行j列元素和v的第j列向量，j選擇為初始偏差協(xié)方差矩陣中徑向位置或橫向速度的對應(yīng)值。

本發(fā)明解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法具有下述優(yōu)點(diǎn)：1、本發(fā)明通過計算相對狀態(tài)偏差的狀態(tài)轉(zhuǎn)移矩陣及張量，計算相對運(yùn)動狀態(tài)偏差均值、協(xié)方差矩陣、及概率密度函數(shù)的非線性傳播結(jié)果，考慮了航天器實(shí)際飛行環(huán)境中的主要攝動因素j2項及二階非線性項，可用于兩相距較遠(yuǎn)航天器的長時間、高精度相對運(yùn)動狀態(tài)偏差傳播分析。2、本發(fā)明利用編隊衛(wèi)星的初始相對運(yùn)動狀態(tài)均值及協(xié)方差矩陣，計算相對狀態(tài)及相對狀態(tài)偏差的狀態(tài)轉(zhuǎn)移矩陣及張量，預(yù)報相對運(yùn)動狀態(tài)偏差的均值、協(xié)方差矩陣、及概率密度函數(shù)，方法解析，計算速度快，可用于衛(wèi)星在軌碰撞預(yù)警及碰撞規(guī)避。

附圖說明

圖1為本發(fā)明實(shí)施例一方法的基本流程示意圖。

圖2為本發(fā)明實(shí)施例一方法和和精確的蒙特卡洛仿真的終端誤差分布對比圖。

圖3為本發(fā)明實(shí)施例一方法用的高斯和模型對均值和協(xié)方差的預(yù)報與蒙特卡洛仿真的樣本點(diǎn)分布對比圖。

圖4為本發(fā)明實(shí)施例二中的終端誤差均值及協(xié)方差矩陣的預(yù)報結(jié)果對比圖。

具體實(shí)施方式

實(shí)施例一：

如圖1所示，本實(shí)施例解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法的實(shí)施步驟包括：

1)為編隊飛行的兩衛(wèi)星指定主衛(wèi)星與從衛(wèi)星，并輸入初始時刻的參考衛(wèi)星絕對軌道狀態(tài)、兩衛(wèi)星標(biāo)稱相對運(yùn)動狀態(tài)、輸入初始相對運(yùn)動狀態(tài)偏差的概率密度函數(shù)；

2)基于初始時刻的參考衛(wèi)星絕對軌道狀態(tài)、兩衛(wèi)星標(biāo)稱相對運(yùn)動狀態(tài)，根據(jù)考慮j2攝動的非線性相對運(yùn)動方程計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的一階及二階狀態(tài)轉(zhuǎn)移張量；

3)基于衛(wèi)星相對運(yùn)動狀態(tài)的一階及二階狀態(tài)轉(zhuǎn)移張量，計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階及二階狀態(tài)轉(zhuǎn)移張量；

4)判斷初始及衛(wèi)星偏差是否均為高斯分布，如果均為高斯分布則跳轉(zhuǎn)執(zhí)行步驟5)，否則，跳轉(zhuǎn)執(zhí)行步驟6)；

5)采用協(xié)方差分析方法，獲得衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值及協(xié)方差矩陣的解析非線性傳播方程，跳轉(zhuǎn)執(zhí)行步驟7)；

6)采用高斯和模型，獲得衛(wèi)星非線性相對運(yùn)動狀態(tài)偏差的概率密度函數(shù)的解析傳播方程，跳轉(zhuǎn)執(zhí)行步驟7)；

7)輸出衛(wèi)星時刻兩衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值、協(xié)方差矩陣及概率密度函數(shù)。

本實(shí)施例解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法通過基于考慮j2攝動的非線性相對運(yùn)動方程，能夠用于兩相距較遠(yuǎn)衛(wèi)星的相對運(yùn)動狀態(tài)偏差的長時間、高精度解析預(yù)報，獲得的偏差信息可用于編隊衛(wèi)星碰撞概率計算及碰撞預(yù)警，具有設(shè)計方法正確合理、方法解析計算量小、對實(shí)際工程任務(wù)的適用性好的優(yōu)點(diǎn)。

本實(shí)施例中，步驟1)中輸入初始時刻的參考衛(wèi)星絕對軌道狀態(tài)具體是指輸入主衛(wèi)星初始時刻t0的絕對軌道狀態(tài)e0＝[a0,e0,i0,ω0,ω0,f0]，其中，e0為主衛(wèi)星初始時刻t0的初始軌道狀態(tài)，a0為主衛(wèi)星初始時刻t0的半長軸，e0為主衛(wèi)星初始時刻t0的偏心率，i0為主衛(wèi)星初始時刻t0的軌道傾角，ω0為主衛(wèi)星初始時刻t0的升交點(diǎn)赤經(jīng)，ω0為主衛(wèi)星初始時刻t0的近拱點(diǎn)角距，f0為主衛(wèi)星初始時刻t0的真近點(diǎn)角；且輸入初始時刻的兩衛(wèi)星標(biāo)稱相對運(yùn)動狀態(tài)為δx(t0)，δx(t0)在主衛(wèi)星當(dāng)?shù)剀壍雷鴺?biāo)系中表示，該坐標(biāo)系原點(diǎn)為主衛(wèi)星質(zhì)心，x軸沿主衛(wèi)星地心矢徑方向，z軸沿軌道面法向，y軸與x、z軸構(gòu)成右手坐標(biāo)系；輸入初始相對運(yùn)動狀態(tài)偏差為dx(t0)及其概率密度函數(shù)為p(t0,dx0)，若初始偏差為高斯分布，則其概率密度函數(shù)為pg(dx0；m0,c0)，m0為初始相對狀態(tài)偏差的均值矩陣，c0為初始相對狀態(tài)偏差的協(xié)方差矩陣；

本實(shí)施例中，步驟2)中考慮j2攝動的非線性相對運(yùn)動方程的函數(shù)表達(dá)式如式(1)所示；

式(1)中，δx(tf)表示時刻tf從衛(wèi)星相對主衛(wèi)星的預(yù)報相對運(yùn)動狀態(tài)，φ(tf,t0)是從t0時刻到tf時刻的狀態(tài)轉(zhuǎn)移矩陣，δx(t0)表示時刻t0從衛(wèi)星相對主衛(wèi)星的預(yù)報相對運(yùn)動狀態(tài)，ψ(tf,t0)為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量，t(tf)為以緯度幅角θf為自變量的無量綱化坐標(biāo)到以時間tf為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，為從θ0緯度幅角到θf緯度幅角的狀態(tài)轉(zhuǎn)移矩陣，t(t0)為以緯度幅角θ0為自變量的無量綱化坐標(biāo)到以時間t0為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，為從θ0緯度幅角到θf緯度幅角的二階狀態(tài)轉(zhuǎn)移張量，t^-1(t0)為轉(zhuǎn)換矩陣t(t0)的逆矩陣，為kronecker張量積；

本實(shí)施例中，步驟2)中計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的一階狀態(tài)轉(zhuǎn)移張量如式(2)所示，計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的二階狀態(tài)轉(zhuǎn)移張量如式(3)所示；

式(2)和式(3)中，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的從θ0緯度幅角到θf緯度幅角的一階狀態(tài)轉(zhuǎn)移張量，表示θf＝θ0時的中間變量表示的逆矩陣，表示θf＝θ0時的中間變量為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)的從θ0緯度幅角到θf緯度幅角的二階狀態(tài)轉(zhuǎn)移張量，為kronecker張量積，和均為公式推導(dǎo)使用的中間變量，且中間變量和的計算函數(shù)式如式(4)所示；

式(4)中，a、b為公式推導(dǎo)使用的中間變量，表示t0時刻平均軌道要素偏差到tf時刻平均軌道要素偏差的一階傳遞矩陣，表示t0時刻軌道要素偏差的無量綱化矩陣，表示tf時刻軌道要素偏差的無量綱化矩陣，表示t0時刻平均軌道要素偏差到tf時刻平均軌道要素偏差的二階傳遞張量，q(θf)表示tf時刻無量綱軌道要素偏差到相對狀態(tài)的二階轉(zhuǎn)換張量，t(θf)表示tf時刻以緯度幅角為自變量的無量綱化坐標(biāo)到以時間為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，t^-1(θf)表示t(θf)的逆矩陣，π表示相對狀態(tài)分量次序變換矩陣，σ(θf)表示tf時刻軌道要素偏差到相對狀態(tài)的轉(zhuǎn)換矩陣，表示平均軌道要素偏差到吻切軌道要素偏差的轉(zhuǎn)換矩陣，為kronecker張量積。本實(shí)施例中，等帶上標(biāo)“bar”的矩陣或張量表示與平均軌道根數(shù)對應(yīng)的值，即將主航天器對應(yīng)時刻的平均軌道根數(shù)帶入矩陣d，γ，h表達(dá)式計算獲得的值，反之沒有上標(biāo)“bar”的矩陣或張量表示與吻切軌道根數(shù)對應(yīng)的值。式(4)中有對應(yīng)初始時刻t0及任意時刻tf的量，分別帶入t0時刻或tf時刻下主航天器的軌道參數(shù)，即可計算對應(yīng)的矩陣及張量。將θf＝θ0帶入式(4)，即可計算及其中in為n維單位矩陣。需要說明的是，前述矩陣t，t^-1，π，γ，p及張量q、h均為已知矩陣，其詳細(xì)表達(dá)式參見文獻(xiàn):senguptap,vadalisr,alfriendkt.second-orderstatetransitionforrelativemotionnearperturbed,ellipticorbits[j].celestialmechanicsanddynamicalastronomy,2006,97(2):101-129。前述矩陣σ，d同樣也均為已知矩陣，其詳細(xì)表達(dá)式參見文獻(xiàn):gimdw,alfriendkt.statetransitionmatrixofrelativemotionfortheperturbednoncircularreferenceorbit[j].journalofguidance,control,anddynamics,2003,26(6):956–971。

本實(shí)施例中，步驟3)計算用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階及二階狀態(tài)轉(zhuǎn)移張量的函數(shù)表達(dá)式如式(5)所示；

式(5)中，dxf為衛(wèi)星時刻的相對運(yùn)動狀態(tài)偏差，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階狀態(tài)轉(zhuǎn)移張量，dx0為初始的相對運(yùn)動狀態(tài)偏差，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的二階狀態(tài)轉(zhuǎn)移張量，為的第i行第j列元素，為的第i行第j列元素，為的第i行k列j維元素，是從t0時刻到tf時刻的狀態(tài)轉(zhuǎn)移矩陣，φ(tf,t0)是從t0時刻到tf時刻的狀態(tài)轉(zhuǎn)移矩陣，ψ(tf,t0)為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量，為t0時刻的標(biāo)稱相對運(yùn)動狀態(tài)δx0的第k個分量，為kronecker張量積。

本實(shí)施例中，步驟5)中獲得的衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值及協(xié)方差矩陣的解析非線性傳播方程的函數(shù)表達(dá)式如式(6)所示；

式(6)中，mⁱ(tf)為衛(wèi)星相對狀態(tài)偏差的均值矩陣m(tf)的第i行元素，m^j(tf)為衛(wèi)星相對狀態(tài)偏差的均值矩陣m(tf)的第j行元素，c^ij(tf)為衛(wèi)星相對狀態(tài)偏差的矩陣，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階狀態(tài)轉(zhuǎn)移張量的第i行第a列元素，為用于解析預(yù)報衛(wèi)星相對運(yùn)動狀態(tài)偏差的一階狀態(tài)轉(zhuǎn)移張量的第j行第b列元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第i行a列b維元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第j行b列c維元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第i行b列c維元素，為從t0時刻到tf時刻的二階狀態(tài)轉(zhuǎn)移張量的第j行c列d維元素，分別為初始相對狀態(tài)偏差dx0的前四階矩。

本實(shí)施例中，初始相對狀態(tài)偏差dx0的前四階矩的計算函數(shù)表達(dá)式如式(7)所示；

式(7)中，分別為初始相對狀態(tài)偏差dx0的前四階矩，分別為初始相對狀態(tài)偏差的均值矩陣m0的a,b,c,d分量元素，分別為初始相對狀態(tài)偏差的協(xié)方差矩陣c0的ab分量元素、bc分量元素、cd分量元素、ad分量元素。本實(shí)施例中，上標(biāo)a,b,c,d為矩陣及張量對應(yīng)的分量元素。對高斯分布，概率密度函數(shù)僅有均值及協(xié)方差矩陣確定，因此終端相對狀態(tài)偏差的概率密度函數(shù)可表示為pg(dxf；mf,cf)。需要說明的是，即使初始相對狀態(tài)偏差為高斯分布，經(jīng)非線性傳播后，終端偏差一般也為非高斯的。因為僅用均值和協(xié)方差矩陣無法完整描述非高斯偏差的概率密度函數(shù)，所以本發(fā)明用高斯和模型來預(yù)報非高斯分布偏差的概率密度函數(shù)。

本實(shí)施例中，步驟6)的詳細(xì)步驟包括：

6.1)用n個高斯分布概率密度函數(shù)去逼近初始相對狀態(tài)偏差的概率密度函數(shù)p(t0,dx0)獲得每一個子高斯分布的權(quán)重、均值和協(xié)方差矩陣；

6.2)采用衛(wèi)星相對運(yùn)動狀態(tài)偏差的均值及協(xié)方差矩陣的解析非線性傳播方程，將每個子高斯分布的均值mi(t0)及協(xié)方差矩陣ci(t0)預(yù)報到衛(wèi)星時刻tf，預(yù)報過程中子高斯分布的權(quán)重保持不變，獲得衛(wèi)星時刻每個子高斯的均值mi(tf)及協(xié)方差矩陣ci(tf)，其中i＝1,2,…,n，n為高斯分布概率密度函數(shù)的數(shù)量；

6.3)針對每一個第i個子高斯分布，其中i＝1,2,…,n，n為高斯分布概率密度函數(shù)的數(shù)量，根據(jù)其權(quán)重系數(shù)ωi、均值mi(tf)及協(xié)方差矩陣ci(tf)，根據(jù)式(8)和式(9)所示獲得衛(wèi)星非線性相對運(yùn)動狀態(tài)偏差的概率密度函數(shù)的解析傳播方程計算衛(wèi)星相對狀態(tài)偏差的均值及協(xié)方差矩陣以及衛(wèi)星偏差的概率密度函數(shù)p(tf,dxf)；

式(8)和(9)中，m(tf),c(tf)為衛(wèi)星相對狀態(tài)偏差dxf的均值及協(xié)方差矩陣，ωi表示第i個子高斯分布的權(quán)重系數(shù)，n為高斯分布概率密度函數(shù)的數(shù)量，ω為所有高斯分布的權(quán)重系數(shù)之和，mi(tf)為tf時刻第i個子高斯分布的均值，ci(tf)為tf時刻第i個子高斯分布的協(xié)方差矩陣，m^t(tf)為mi(tf)的轉(zhuǎn)置矩陣，p(tf,dxf)為衛(wèi)星偏差的概率密度函數(shù)，pg(dxf；mi(tf),ci(tf))表示tf時刻第i個子高斯分布的概率密度函數(shù)。

本實(shí)施例中，步驟6.1)中逼近初始相對狀態(tài)偏差的概率密度函數(shù)p(t0,dx0)采用的逼近公式如式(10)所示；

式(10)中，p(t0,dx0)為初始相對狀態(tài)偏差的概率密度函數(shù)，dx0為初始相對狀態(tài)偏差，pg(dx0；mi,ci)表示第i個子高斯分布的概率密度函數(shù)，mi和ci分別為相應(yīng)的均值和協(xié)方差矩陣，ωi表示第i個子高斯分布的權(quán)重系數(shù)。

本實(shí)施例中，第i個子高斯分布的權(quán)重系數(shù)ωi的計算函數(shù)表達(dá)式如式(11)所示；

式(11)中，ωi表示第i個子高斯分布的權(quán)重系數(shù)，m0、c0為初始相對狀態(tài)偏差的均值及協(xié)方差矩陣，及為將一維標(biāo)準(zhǔn)正態(tài)分布pg(x；0,1)等方差地分割為n個子正態(tài)分布的權(quán)重、均值及標(biāo)準(zhǔn)差，diag{a}表示以a向量為對角線元素的矩陣，v和λ分別為對c0進(jìn)行特征值分解所得的特征向量矩陣與特征值矩陣，λ1,λ2,...,λn為c0的特征值，λj與υj分別為沿選定的j方向的特征值與特征向量，即λ的第j行j列元素和v的第j列向量，j選擇為初始偏差協(xié)方差矩陣中徑向位置或橫向速度的對應(yīng)值。需要說明的是，及可通過優(yōu)化方法求解，也可直接調(diào)用現(xiàn)有的數(shù)據(jù)庫獲得，見文獻(xiàn)vittaldev,v.,andrussell,r.p.,spaceobjectcollisionprobabilityusingmultidirectionalgaussianmixturemodels[j].journalofguidance,control,anddynamics,2016,39(9):2163-2169。

綜上所述，本實(shí)施例解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法為一種考慮j2攝動、解析的衛(wèi)星非線性相對運(yùn)動偏差傳播分析方法，利用編隊衛(wèi)星的初始相對運(yùn)動狀態(tài)均值及協(xié)方差矩陣，計算相對狀態(tài)及相對狀態(tài)偏差的狀態(tài)轉(zhuǎn)移矩陣及張量，預(yù)報相對運(yùn)動狀態(tài)偏差的均值、協(xié)方差矩陣、及概率密度函數(shù)。其中，可根據(jù)初始偏差及終端偏差是否為高斯分布，選擇是否采用高斯和模型。該方法考慮了j2攝動項和二階非線性項、可用于相距較遠(yuǎn)航天器的長時間、高精度相對運(yùn)動狀態(tài)偏差傳播，具有設(shè)計方法正確合理，方法解析、計算速度快，對實(shí)際工程任務(wù)的適用性好等優(yōu)點(diǎn)。

參見圖2所示終端誤差分布對比可知，和精確的蒙特卡洛仿真相比，本方法用的高斯和模型的所有子高斯分布很好的包圍了的蒙特卡洛仿真的樣本點(diǎn)，說明本方法精確可靠。進(jìn)一步由圖3可知，本方法用的高斯和模型對均值和協(xié)方差的預(yù)報與蒙特卡洛仿真一致，其概率密度等高線與蒙特卡洛仿真的樣本點(diǎn)分布一致，且很好的表征了終端偏差的非高斯分布特性，說明本方法在非線性、非高斯偏差均值、協(xié)方差矩陣及概率密度函數(shù)的傳播分析中具有較高的精度。

實(shí)施例二：

本實(shí)施例與實(shí)施例一基本相同，其主要不同點(diǎn)為：本實(shí)施例中，步驟3)中初始及終端偏差均為高斯分布，因此不再進(jìn)行步驟6)的計算，僅需要采用步驟5)對非線性相對狀態(tài)偏差的均值及協(xié)方差矩陣進(jìn)行解析預(yù)報，終端偏差的概率密度函數(shù)由其均值及協(xié)方差矩陣即可確定。本實(shí)施例僅需要預(yù)報相對狀態(tài)偏差的均值及協(xié)方差矩陣，不需要采用高斯和模型預(yù)報概率密度函數(shù)，方法解析、計算效率高，實(shí)施簡單；計算效率高。參見圖4所示終端誤差均值及協(xié)方差矩陣的預(yù)報結(jié)果對比可知，本本實(shí)施例方法用的高斯和模型對均值和協(xié)方差矩陣的預(yù)報與蒙特卡洛仿真一致，相比已有的線性方法，本方法精度顯著提高。

實(shí)施例三：

本實(shí)施例與實(shí)施例一基本相同，其主要不同點(diǎn)為：本實(shí)施例中，不再進(jìn)行步驟2)的計算，將步驟3)中(6)式的矩陣直接用步驟1)中(1)式的矩陣替換。因為衛(wèi)星絕對狀態(tài)偏差可作為一個相對狀態(tài)，采用(1)式進(jìn)行預(yù)報，所以本實(shí)施例可用于解析地預(yù)報衛(wèi)星絕對狀態(tài)偏差的均值、協(xié)方差矩陣及概率密度函數(shù)。本實(shí)施例考慮了j2攝動項和二階非線性項、可用于單個航天器長時間、高精度的絕對狀態(tài)偏差傳播，方法解析、計算效率高。

以上所述僅是本發(fā)明的優(yōu)選實(shí)施方式，本發(fā)明的保護(hù)范圍并不僅局限于上述實(shí)施例，凡屬于本發(fā)明思路下的技術(shù)方案均屬于本發(fā)明的保護(hù)范圍。應(yīng)當(dāng)指出，對于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明原理前提下的若干改進(jìn)和潤飾，這些改進(jìn)和潤飾也應(yīng)視為本發(fā)明的保護(hù)范圍。