本發(fā)明涉及航天器相對運動軌跡預(yù)報技術(shù)，具體涉及一種基于最小二乘的非線性相對軌跡預(yù)報方法。

背景技術(shù)：

航天器相對運動軌跡預(yù)報在航天器編隊、集群等飛行控制任務(wù)中具有重要應(yīng)用。當(dāng)兩編隊航天器飛離測控區(qū)域時，需要對其相對運動軌跡進(jìn)行預(yù)報，以判斷是否需要進(jìn)行構(gòu)型保持、碰撞規(guī)避等操作。

用于描述航天器相對運動的方程主要有基于近圓軌道的C-W方程和基于橢圓軌道的T-H方程。但這兩個線性方程僅適用于二體假設(shè)下近距離的相對軌跡預(yù)報，作為改進(jìn)，不少學(xué)者相繼研究了考慮二階非線性項、J2攝動項等的非線性相對運動方程，可以用于更長時間、更高精度的相對軌跡預(yù)報。但這些解析的相對運動方程仍然只是對實際飛行力學(xué)環(huán)境的一定近似，即存在模型誤差。

再者，已有研究表明，相對軌跡預(yù)報精度不僅與采用的方程有關(guān)，還受初始條件(如參考航天器位置、初始相對狀態(tài)等)影響。實際任務(wù)中，一般會有很多歷史數(shù)據(jù)(即預(yù)報的初始條件)可以利用。現(xiàn)有相對運動方程都是基于單個數(shù)據(jù)點進(jìn)行相對軌跡預(yù)報，且沒有哪一個方程能對不同的初始點都具有優(yōu)于其他方程的預(yù)報精度。

因此，若直接采用現(xiàn)有相對運動方程進(jìn)行相對軌跡預(yù)報，一方面歷史數(shù)據(jù)沒有得到充分利用；另一方面也不知道應(yīng)該選擇哪個初始點和那個相對運動方程，才能獲得最高的預(yù)報精度。本發(fā)明基于最小二乘思想，利用歷史數(shù)據(jù)估計出新的預(yù)報初始條件，可以有效解決以上兩個問題，最小化相對運動方程中的模型誤差與初始條件中的測量誤差傳播造成的預(yù)報誤差，以獲得最高的相對軌跡預(yù)報精度。

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

本發(fā)明要解決的技術(shù)問題：針對現(xiàn)有技術(shù)的上述問題，提供一種考慮了J2攝動項和二階非線性項，可用于相距較遠(yuǎn)航天器的長時間、高精度相對軌跡預(yù)報，具有設(shè)計方法正確合理、對實際工程任務(wù)的適用性好的基于最小二乘的非線性相對軌跡預(yù)報方法。

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

一種基于最小二乘的非線性相對軌跡預(yù)報方法，步驟包括：

1)將已有的n個歷史歷元時刻對應(yīng)的、在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示的主航天器和從航天器的相對運動狀態(tài)歷史軌道數(shù)據(jù)寫入觀察向量；

2)利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量；

3)采用最小二乘方法迭代求解使殘差向量的殘差和最小的最優(yōu)初始相對運動狀態(tài)

4)將所述最優(yōu)初始相對運動狀態(tài)帶入所述相對運動方程，向前預(yù)報到指定的歷元時刻t_f，獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果。

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

1.1)分別指定編隊飛行的兩航天器為主航天器與從航天器，進(jìn)行相對軌跡預(yù)報的初始時刻為t₀，主航天器初始軌道要素為E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，其中，E為主航天器初始軌道，t₀為相對軌跡預(yù)報的初始時刻，a為半長軸，e為偏心率，i為軌道傾角，Ω為升交點赤經(jīng)，ω為近拱點角距，f為真近點角；令從航天器相對主航天器的初始相對運動狀態(tài)為x(t₀)，其中x(t₀)在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示，該坐標(biāo)系原點為主航天器質(zhì)心，x軸沿主航天器地心矢徑方向，z軸沿軌道面法向，y軸與x、z軸構(gòu)成右手坐標(biāo)系；

1.2)獲取n個歷史歷元時刻從航天器相對主航天器的相對運動狀態(tài)，記為：[t_-n,x(t_-n)],[t_-(n-1),x(t_-(n-1))],…,[t_-1,x(t_-1)]，其中x(t_-n)表示時刻t_-n從航天器相對主航天器的初始相對運動狀態(tài)，x(t_-(n-1))表示時刻t_-(n-1)從航天器相對主航天器的初始相對運動狀態(tài)，x(t_-1)表示時刻t_-1從航天器相對主航天器的初始相對運動狀態(tài)；

1.3)將n個歷史歷元時刻從航天器相對主航天器的初始相對運動狀態(tài)歷史軌道數(shù)據(jù)寫入觀察向量Z＝[x(t_-n),x(t_-(n-1)),...,x(t_-1)]，其中，x(t_-n)表示時刻t_-n從航天器相對主航天器的相對運動狀態(tài)，x(t_-(n-1))表示時刻t_-(n-1)從航天器相對主航天器的相對運動狀態(tài)，x(t_-1)表示時刻t_-1從航天器相對主航天器的相對運動狀態(tài)。

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

2.1)利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，得到各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量其中，表示時刻t_-n從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-(n-1)從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-1從航天器相對主航天器的預(yù)報相對運動狀態(tài)；

2.2)根據(jù)式(2.2-1)計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量；

式(2.2-1)中，e表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量，Z表示觀察向量，Y表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，x(t_-n)表示時刻t_-n從航天器相對主航天器的相對運動狀態(tài)，x(t_-(n-1))表示時刻t_-(n-1)從航天器相對主航天器的相對運動狀態(tài)，x(t_-1)表示時刻t_-1從航天器相對主航天器的相對運動狀態(tài)；表示時刻t_-n從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-(n-1)從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-1從航天器相對主航天器的預(yù)報相對運動狀態(tài)。

優(yōu)選地，所述步驟2.1)中利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量時，如果僅已知主航天器的軌道半長軸a(t₀)或平均軌道角速度其中μ為地心引力常數(shù)，a為半長軸，則利用式(2.1-1)所示C-W方程計算各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量；如果已知主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，則利用式(2.1-2)所示考慮J2攝動的非線性相對運動方程計算各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量；

式(2.1-1)中，表示時刻t從航天器相對主航天器的預(yù)報相對運動狀態(tài)，Φ(t,t₀)是從t₀時刻到t時刻的狀態(tài)轉(zhuǎn)移矩陣，x₀表示t₀時刻兩航天器的初始相對運動狀態(tài)，τ＝n(t-t₀)，s＝sinτ，c＝cosτ，為主航天器的平均軌道角速度，t₀為相對軌跡預(yù)報的初始時刻，a為半長軸，μ為地心引力常數(shù)；

式(2.1-2)中，表示時刻t從航天器相對主航天器的預(yù)報相對運動狀態(tài)，Φ(t,t₀)是從t₀時刻到t時刻的狀態(tài)轉(zhuǎn)移矩陣，為從θ₀緯度幅角到θ緯度幅角的狀態(tài)轉(zhuǎn)移矩陣，x₀表示t₀時刻兩航天器的初始相對運動狀態(tài)，Ψ(t,t₀)為從t₀時刻到t時刻的二階狀態(tài)轉(zhuǎn)移張量，為從θ₀緯度幅角到θ緯度幅角的二階狀態(tài)轉(zhuǎn)移張量，T(t)為以緯度幅角θ為自變量的無量綱化坐標(biāo)到以時間t為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，T(t₀)為以緯度幅角θ為自變量的無量綱化坐標(biāo)到以時間t₀為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，T^-1(t₀)為轉(zhuǎn)換矩陣T(t₀)的逆，為Kronecker張量積。

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

3.1)定義最小二乘方法的目標(biāo)函數(shù)為式(3.1-1)所示的殘差和，使得式(3.1-1)所示殘差和最小的條件為式(3.1-2)；

式(3.1-1)中，J表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量的殘差和，Y表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，Z表示觀察向量；

式(3.1-2)中，J表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量的殘差和，x₀表示t₀時刻兩航天器的初始相對運動狀態(tài)，Y表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，e表示各歷史歷元時刻的預(yù)報結(jié)果和歷史軌道數(shù)據(jù)中歷史數(shù)據(jù)的殘差向量；

3.2)將相對運動方程、各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量帶入式(3.1-2)，取一階近似可得初始相對運動狀態(tài)的最優(yōu)估計值滿足式(3.2-1)所示的迭代條件進(jìn)行迭代求解；

式(3.2-1)中，k為迭代次數(shù)，表示第k+1次迭代的結(jié)果，表示第k次迭代的結(jié)果，的迭代初值x(t₀)為時刻t₀從航天器相對主航天器的初始相對運動狀態(tài)，F(xiàn)^k向量的表達(dá)式為F^k＝[(Φ^k)^TΦ^k]^-1(Φ^k)^T，Φ^k表示第k次迭代中的一階狀態(tài)轉(zhuǎn)移矩陣，由式(2.1-1)及(2.1-2)可知，Φ^k在每次迭代中為常量；

3.3)判斷是否滿足指定的迭代條件，如果不滿足指定迭代條件，則跳轉(zhuǎn)執(zhí)行步驟3.2)；否則當(dāng)滿足指定迭代條件時，將最后一次迭代得到的第k+1次迭代的結(jié)果作為最終使得主航天器和從航天器之間相對軌跡預(yù)報殘差的殘差和最小的最優(yōu)初始相對運動狀態(tài)

優(yōu)選地，所述步驟3.3)中指定的迭代條件具體是指第k+1次迭代的結(jié)果與第k次迭代的結(jié)果之間的變化量小于或等于給定迭代誤差、或者迭代次數(shù)k達(dá)到給定的最大值K。

優(yōu)選地，所述步驟4)獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果如式(4-1)或式(4-2)所示；

式(4-1)中，x(t_f)表示指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果，Φ(t_f,t₀)是從t₀時刻到t_f時刻的狀態(tài)轉(zhuǎn)移矩陣，表示步驟3.3)中獲得的最優(yōu)初始相對運動狀態(tài)；

式(4-2)中，x(t_f)表示指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果，Φ(t_f,t₀)表示是從t₀時刻到t_f時刻的狀態(tài)轉(zhuǎn)移矩陣，Ψ(t_f,t₀)表示從t₀時刻到t_f時刻的二階狀態(tài)轉(zhuǎn)移張量，表示步驟3.3)中獲得的最優(yōu)初始相對運動狀態(tài)。

本發(fā)明基于最小二乘的非線性相對軌跡預(yù)報方法的優(yōu)點如下：

1、本發(fā)明基于最小二乘的非線性相對軌跡預(yù)報方法利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量，考慮了航天器實際飛行環(huán)境中的主要攝動因素J2項及二階非線性項，可用于兩相距較遠(yuǎn)航天器的長時間、高精度相對軌跡預(yù)報。

2、本發(fā)明利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量，采用最小二乘方法迭代求解使殘差向量的殘差和最小的最優(yōu)初始相對運動狀態(tài)將最優(yōu)初始相對運動狀態(tài)帶入相對運動方程，向前預(yù)報到指定的歷元時刻t_f，獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果通過解析的相對運動方程進(jìn)行相對軌跡預(yù)報，計算速度快。

附圖說明

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

圖2為本發(fā)明實施例方法的相對軌跡預(yù)報位置誤差對比示意圖。

圖3為本發(fā)明實施例方法的相對軌跡預(yù)報速度誤差對比示意圖。

具體實施方式

實施例一：

如圖1所示，本實施例基于最小二乘的非線性相對軌跡預(yù)報方法的步驟包括：

1)準(zhǔn)備歷史軌道數(shù)據(jù)：

將已有的n個歷史歷元時刻對應(yīng)的、在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示的主航天器和從航天器的相對運動狀態(tài)歷史軌道數(shù)據(jù)寫入觀察向量；

2)計算各歷史歷元時刻預(yù)報結(jié)果和歷史軌道數(shù)據(jù)之間的殘差：

利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量；

3)采用最小二乘方法迭代求解最優(yōu)初始相對運動狀態(tài)：

采用最小二乘方法迭代求解使殘差向量的殘差和最小的最優(yōu)初始相對運動狀態(tài)

4)帶入最優(yōu)初始相對運動狀態(tài)進(jìn)行相對運動軌跡預(yù)報：

將最優(yōu)初始相對運動狀態(tài)帶入相對運動方程，向前預(yù)報到指定的歷元時刻t_f，獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果。

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

1.1)分別指定編隊飛行的兩航天器為主航天器與從航天器，進(jìn)行相對軌跡預(yù)報的初始時刻為t₀，主航天器初始軌道要素為E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，其中，E為主航天器初始軌道，t₀為相對軌跡預(yù)報的初始時刻，a為半長軸，e為偏心率，i為軌道傾角，Ω為升交點赤經(jīng)，ω為近拱點角距，f為真近點角；令從航天器相對主航天器的初始相對運動狀態(tài)為x(t₀)，其中x(t₀)在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示，該坐標(biāo)系原點為主航天器質(zhì)心，x軸沿主航天器地心矢徑方向，z軸沿軌道面法向，y軸與x、z軸構(gòu)成右手坐標(biāo)系；

1.2)獲取n個歷史歷元時刻從航天器相對主航天器的相對運動狀態(tài)，記為：[t_-n,x(t_-n)],[t_-(n-1),x(t_-(n-1))],…,[t_-1,x(t_-1)]，其中x(t_-n)表示時刻t_-n從航天器相對主航天器的初始相對運動狀態(tài)，x(t_-(n-1))表示時刻t_-(n-1)從航天器相對主航天器的初始相對運動狀態(tài)，x(t_-1)表示時刻t_-1從航天器相對主航天器的初始相對運動狀態(tài)；

1.3)將n個歷史歷元時刻從航天器相對主航天器的初始相對運動狀態(tài)歷史軌道數(shù)據(jù)寫入觀察向量Z＝[x(t_-n),x(t_-(n-1)),...,x(t_-1)]，其中，x(t_-n)表示時刻t_-n從航天器相對主航天器的相對運動狀態(tài)，x(t_-(n-1))表示時刻t_-(n-1)從航天器相對主航天器的相對運動狀態(tài)，x(t_-1)表示時刻t_-1從航天器相對主航天器的相對運動狀態(tài)。本實施例中，觀察向量具體為6n×1的觀測列向量，通過前述步驟1.3)最終將已有的n個歷史歷元時刻對應(yīng)的、在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示的兩航天器相對運動狀態(tài)寫入一個6n×1的觀測列向量。

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

2.1)利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，得到各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量其中，表示時刻t_-n從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-(n-1)從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-1從航天器相對主航天器的預(yù)報相對運動狀態(tài)；

2.2)根據(jù)式(2.2-1)計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量；

式(2.2-1)中，e表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量，Z表示觀察向量，Y表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，x(t_-n)表示時刻t_-n從航天器相對主航天器的相對運動狀態(tài)，x(t_-(n-1))表示時刻t_-(n-1)從航天器相對主航天器的相對運動狀態(tài)，x(t_-1)表示時刻t_-1從航天器相對主航天器的相對運動狀態(tài)；表示時刻t_-n從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-(n-1)從航天器相對主航天器的預(yù)報相對運動狀態(tài)，表示時刻t_-1從航天器相對主航天器的預(yù)報相對運動狀態(tài)。

本實施例中，步驟2.1)中利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量時，如果僅已知主航天器的軌道半長軸a(t₀)或平均軌道角速度其中μ為地心引力常數(shù)，a為半長軸，則利用式(2.1-1)所示C-W方程計算各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量；如果已知主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，則利用式(2.1-2)所示考慮J2攝動的非線性相對運動方程計算各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量；

式(2.1-1)中，表示時刻t從航天器相對主航天器的預(yù)報相對運動狀態(tài)，Φ(t,t₀)是從t₀時刻到t時刻的狀態(tài)轉(zhuǎn)移矩陣，x₀表示t₀時刻兩航天器的初始相對運動狀態(tài)，τ＝n(t-t₀)，s＝sinτ，c＝cosτ，為主航天器的平均軌道角速度，t₀為相對軌跡預(yù)報的初始時刻，a為半長軸，μ為地心引力常數(shù)；

式(2.1-2)中，表示時刻t從航天器相對主航天器的預(yù)報相對運動狀態(tài)，Φ(t,t₀)是從t₀時刻到t時刻的狀態(tài)轉(zhuǎn)移矩陣，為從θ₀緯度幅角到θ緯度幅角的狀態(tài)轉(zhuǎn)移矩陣，x₀表示t₀時刻兩航天器的初始相對運動狀態(tài)，Ψ(t,t₀)為從t₀時刻到t時刻的二階狀態(tài)轉(zhuǎn)移張量，為從θ₀緯度幅角到θ緯度幅角的二階狀態(tài)轉(zhuǎn)移張量，T(t)為以緯度幅角θ為自變量的無量綱化坐標(biāo)到以時間t為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，T(t₀)為以緯度幅角θ₀為自變量的無量綱化坐標(biāo)到以時間t₀為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，反之量綱化坐標(biāo)到無量綱化坐標(biāo)的轉(zhuǎn)換矩陣T^-1(t₀)為轉(zhuǎn)換矩陣T(t₀)的逆，為Kronecker張量積。例如，表示用矩陣X的每一個元素與矩陣Y相乘，所得結(jié)果為mp×nq的矩陣Z。本實施例中利用相對運動方程將初始相對狀態(tài)x₀向后分別預(yù)報到每個歷史數(shù)據(jù)對應(yīng)的歷元時刻，即[t_-n,t_-(n-1),...,t_-1]，獲得預(yù)報的相對運動狀態(tài)向量：若僅知道主航天器的軌道半長軸a(t₀)(或平均軌道角速度利用C-W方程(2.1-1)預(yù)報；若主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]均已知，則利用考慮J2攝動的非線性相對運動方程(2.1-2)預(yù)報。

式(2.1-2)可表示為如式(2.1-3)所示分量形式；

式(2.1-3)中，表示t時刻預(yù)報相對運動狀態(tài)的第i個分量，Φ_i,a(t,t₀)表示t₀時刻到t時刻的一階狀態(tài)轉(zhuǎn)移矩陣的第i行、第a列元素，Ψ_i,ab(t,t₀)表示t₀時刻到t時刻的二階狀態(tài)轉(zhuǎn)移張量的第i模塊、第a行、第b列元素，x_a(t₀)表示初始t₀時刻相對運動狀態(tài)的第a個分量，x_b(t₀)表示初始t₀時刻相對運動狀態(tài)的第b個分量，Φ_ij(t,t₀)表示t₀時刻到t時刻的一階狀態(tài)轉(zhuǎn)移矩陣的第i行、第j列元素，T_il(t)表示以緯度幅角θ為自變量的無量綱化坐標(biāo)到以時間t為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣的第i行、第l列元素，表示從θ₀緯度幅角到θ緯度幅角的一階狀態(tài)轉(zhuǎn)移矩陣的第l行、第m列元素，T_mj(t₀)表示以緯度幅角為自變量的無量綱化坐標(biāo)到以時間為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣的第m行、第j列元素，Ψ_ijk(t,t₀)表示t₀時刻到t時刻的二階狀態(tài)轉(zhuǎn)移張量的第i模塊、第j行、第k列元素，表示從θ₀緯度幅角到θ緯度幅角的二階狀態(tài)轉(zhuǎn)移張量的第l模塊、第m行、第n列元素，T_nk(t₀)表示t₀時刻以緯度幅角為自變量的無量綱化坐標(biāo)到以時間為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣的第n行、第k列元素。式(2.1-3)采用愛因斯坦求和助記，省略對啞變量的求和符號，即下標(biāo)a表示狀態(tài)x(t₀)的第a個分量，N＝6為狀態(tài)變量的維數(shù)。啞變量a,b,i,j,k,m,n,l取值均為集合{1,2,…,N}。

狀態(tài)轉(zhuǎn)移矩陣與二階狀態(tài)轉(zhuǎn)移張量可用式(2.1-4)求解。

式(2.1-4)中，表示從θ₀緯度幅角到θ緯度幅角的一階狀態(tài)轉(zhuǎn)移矩陣的第i行、第j列元素，表示狀態(tài)轉(zhuǎn)移矩陣的第i行、第l列元素，表示狀態(tài)轉(zhuǎn)移矩陣的第l行、第j列元素，表示狀態(tài)轉(zhuǎn)移張量的第i模塊、第l行、第m列元素，表示狀態(tài)轉(zhuǎn)移矩陣的第m行、第k列元素，表示狀態(tài)轉(zhuǎn)移張量的第i模塊、第l行、第m列元素，表示狀態(tài)轉(zhuǎn)移矩陣的第m行、第j列元素，表示狀態(tài)轉(zhuǎn)移矩陣的第n行、第k列元素，啞變量a,b,i,j,k,m,n,l取值均為集合{1,2,…,N}。

狀態(tài)轉(zhuǎn)移矩陣與二階狀態(tài)轉(zhuǎn)移張量可用式(2.1-5)求解：

式(2.1-5)中，表示公式推導(dǎo)使用的中間變量，沒有物理意義，T(t)表示t時刻以緯度幅角為自變量的無量綱化坐標(biāo)到以時間為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，Π表示相對狀態(tài)分量次序變換矩陣，Σ(t)表示t時刻軌道要素偏差到相對狀態(tài)的轉(zhuǎn)換矩陣，表示平均軌道要素偏差到吻切軌道要素偏差的轉(zhuǎn)換矩陣，表示t₀時刻平均軌道要素偏差到t時刻平均軌道要素偏差的一階傳遞矩陣，表示t₀時刻軌道要素偏差的無量綱化矩陣，表示公式推導(dǎo)使用的中間變量的第i模塊、第j行、第k列元素，A表示公式推導(dǎo)使用的中間變量，沒有物理意義，B表示公式推導(dǎo)使用的中間變量，沒有物理意義，A_il表示中間變量A的第i行、第l列元素，B_lj表示中間變量B的第l行、第j列元素，B_mk表示中間變量B的第m行、第k列元素，表示t₀時刻平均軌道要素偏差到t時刻平均軌道要素偏差的二階傳遞張量的第l模塊、第j行、第k列元素，Q_ilm(t)表示t時刻無量綱軌道要素偏差到相對狀態(tài)的二階轉(zhuǎn)換張量的第i模塊、第l行、第m列元素，表示t時刻軌道要素偏差的無量綱化矩陣。等帶上標(biāo)“bar”的矩陣或張量表示與平均軌道根數(shù)對應(yīng)的值，即將主航天器對應(yīng)時刻的平均軌道根數(shù)帶入矩陣D，Γ，H表達(dá)式計算獲得的值，反之沒有上標(biāo)“bar”的矩陣或張量表示與吻切軌道根數(shù)對應(yīng)的值。以上表達(dá)式中有對應(yīng)初始時刻t₀及任意時刻t的量，分別帶入t₀時刻或t時刻下主航天器的軌道參數(shù)，即可計算對應(yīng)的矩陣及張量。T為以緯度幅角θ為自變量的無量綱化坐標(biāo)到以時間t為自變量的量綱化坐標(biāo)的轉(zhuǎn)換矩陣，反之量綱化坐標(biāo)到無量綱化坐標(biāo)的轉(zhuǎn)換矩陣T^-1為T的逆；Π為相對狀態(tài)分量次序變換矩陣；Σ為軌道要素偏差到相對狀態(tài)的轉(zhuǎn)換矩陣；為初始時刻平均軌道要素偏差到任意時刻平均軌道要素偏差的一階傳遞矩陣；D為平均軌道要素偏差到吻切軌道要素偏差的轉(zhuǎn)換矩陣；Γ為軌道要素偏差的無量綱化矩陣；P為無量綱軌道要素偏差到相對狀態(tài)的一階轉(zhuǎn)換矩陣；Q無量綱軌道要素偏差到相對狀態(tài)的二階轉(zhuǎn)換張量；H為初始時刻平均軌道要素偏差到任意時刻平均軌道要素偏差的二階傳遞張量。需要說明的是，前述矩陣T，T^-1，Π，Γ，P及張量Q、H均為已知矩陣，其詳細(xì)表達(dá)式參見文獻(xiàn)[1]:Sengupta P,Vadali S R,Alfriend K T.Second-order state transition for relative motion near perturbed,elliptic orbits[J].Celestial Mechanics and Dynamical Astronomy,2006,97(2)∶101-129。前述矩陣Σ，D同樣也均為已知矩陣，其詳細(xì)表達(dá)式參見文獻(xiàn)[2]:Gim D W,Alfriend K T.State Transition Matrix of Relative Motion for the Perturbed Noncircular Reference Orbit[J].Journal of Guidance,Control,and Dynamics,2003,26(6)∶956–971。

利用相對運動方程將初始相對狀態(tài)x₀向后分別預(yù)報到每個歷史數(shù)據(jù)對應(yīng)的歷元時刻，即[t_-n,t_-(n-1),...,t_-1]，獲得預(yù)報的相對運動狀態(tài)向量：若僅知道主航天器的軌道半長軸a(t₀)(或平均軌道角速度利用C-W方程(1a)預(yù)報；若主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]均已知，則利用考慮J2攝動的非線性相對運動方程(2.2-2)預(yù)報。

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

3.1)定義最小二乘方法的目標(biāo)函數(shù)為式(3.1-1)所示的殘差和，使得式(3.1-1)所示殘差和最小的條件為式(3.1-2)；

式(3.1-1)中，J表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量的殘差和，Y表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，Z表示觀察向量；

式(3.1-2)中，J表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量的殘差和，x₀表示t₀時刻兩航天器的初始相對運動狀態(tài)，Y表示各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量，e表示各歷史歷元時刻的預(yù)報結(jié)果和歷史軌道數(shù)據(jù)中歷史數(shù)據(jù)的殘差向量；

3.2)將相對運動方程、各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量和觀察向量之間的殘差向量帶入式(3.1-2)，取一階近似可得初始相對運動狀態(tài)的最優(yōu)估計值滿足式(3.2-1)所示的迭代條件進(jìn)行迭代求解；

式(3.2-1)中，k為迭代次數(shù)，表示第k+1次迭代的結(jié)果，表示第k次迭代的結(jié)果，的迭代初值x(t₀)為時刻t₀從航天器相對主航天器的初始相對運動狀態(tài)，F(xiàn)^k向量的表達(dá)式為F^k＝[(Φ^k)^TΦ^k]^-1(Φ^k)^T，Φ^k表示第k次迭代中的一階狀態(tài)轉(zhuǎn)移矩陣，由式(2.1-1)及(2.1-2)可知，Φ^k在每次迭代中為常量；

3.3)判斷是否滿足指定的迭代條件，如果不滿足指定迭代條件，則跳轉(zhuǎn)執(zhí)行步驟3.2)；否則當(dāng)滿足指定迭代條件時，將最后一次迭代得到的第k+1次迭代的結(jié)果作為最終使得主航天器和從航天器之間相對軌跡預(yù)報殘差的殘差和最小的最優(yōu)初始相對運動狀態(tài)

本實施例中，步驟3.3)中指定的迭代條件具體是指第k+1次迭代的結(jié)果與第k次迭代的結(jié)果之間的變化量小于或等于給定迭代誤差、或者迭代次數(shù)k達(dá)到給定的最大值K。例如或k≤20；則結(jié)束迭代，獲得最優(yōu)初始相對運動狀態(tài)的估計值：

將基于最小二乘方法計算得到的新初始條件帶入原相對運動方程，向前預(yù)報到需要的歷元時刻，獲得具有較高精度的相對運動軌跡預(yù)報結(jié)果。

本實施例中，步驟4)獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果如式(4-1)或式(4-2)所示；

式(4-1)中，x(t_f)表示指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果，Φ(t_f,t₀)是從t₀時刻到t_f時刻的狀態(tài)轉(zhuǎn)移矩陣，表示步驟3.3)中獲得的最優(yōu)初始相對運動狀態(tài)；

式(4-2)中，xⁱ(t_f)表示指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果，Φ(t_f,t₀)表示是從t₀時刻到t_f時刻的狀態(tài)轉(zhuǎn)移矩陣，Ψ(t_f,t₀)表示從t₀時刻到t_f時刻的二階狀態(tài)轉(zhuǎn)移張量，表示步驟3.3)中獲得的最優(yōu)初始相對運動狀態(tài)。其中，方程(4-1)為基于線性C-W方程的預(yù)報，方程(4-2)為考慮J2攝動的非線性相對運動方程預(yù)報結(jié)果。仿真結(jié)果表明，對同樣的相對運動方程，采用作為初始條件比直接采用x(t₀)作為初始條件的相對軌跡預(yù)報精度要高。而對同樣的歷史數(shù)據(jù)，基于方程(4-2)的預(yù)報結(jié)果要比基于方程(4-1)的預(yù)報結(jié)果精度高。

綜上所述，本實施例基于最小二乘的非線性相對軌跡預(yù)報方法為一種基于最小二乘的非線性相對軌跡預(yù)報精度改進(jìn)方法，利用歷史軌道數(shù)據(jù)，基于最小二乘原理最小化預(yù)報值與歷史值的殘差和，對相對軌跡預(yù)報的初始條件進(jìn)行估計改進(jìn)，以提高相對軌跡預(yù)報的精度。其中，所用的相對運動方程可根據(jù)具體情況選擇線性C-W方程或考慮J2攝動的非線性相對運動方程，考慮了J2攝動項和二階非線性項、可用于相距較遠(yuǎn)航天器的長時間、高精度相對軌跡預(yù)報，具有設(shè)計方法正確合理、對實際工程任務(wù)的適用性好等優(yōu)點。

實施例二：

本實施例與實施例一基本相同，其主要不同點為：本實施例中，步驟2.1)中利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量時，僅已知主航天器的軌道半長軸a(t₀)或平均軌道角速度其中μ為地心引力常數(shù)，a為半長軸，且利用式(2.1-1)所示C-W方程計算各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量；步驟4)獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果如式(4-1)所示。

本實施例僅根據(jù)歷史相對運動狀態(tài)數(shù)據(jù)、初始相對運動狀態(tài)及主航天器半長軸(或平均軌道角速度)就可以進(jìn)行相對軌跡預(yù)報。預(yù)報過程為解析的線性C-W方程，通過最小二乘估計新的初始條件能最小化預(yù)報誤差，預(yù)報精度得到改進(jìn)；且基于最小二乘的迭代過程均為解析計算，能夠保證迭代過程快速收斂，計算效率高。參見圖2所示相對軌跡預(yù)報位置誤差對比可知，和現(xiàn)有技術(shù)采用的CW方程直接預(yù)報相比，本實施例(基于CW方程的改進(jìn))的相對距離誤差x更?。粎⒁妶D3所示相對軌跡預(yù)報速度誤差對比可知，和現(xiàn)有技術(shù)采用的CW方程直接預(yù)報相比，本實施例(基于CW方程的改進(jìn))的相對速度誤差x更小。

實施例三：

本實施例與實施例一基本相同，其主要不同點為：本實施例中，步驟2.1)中利用相對運動方程計算各歷史歷元時刻的預(yù)報相對運動狀態(tài)向量時，已知主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，且利用式(2.1-2)所示考慮J2攝動的非線性相對運動方程計算各歷史歷元時刻[t_-n,t_-(n-1),...,t_-1]的預(yù)報相對運動狀態(tài)向量；步驟4)獲得指定的歷元時刻t_f的相對運動軌跡預(yù)報結(jié)果如式(4-2)所示。

本實施例用于初始條件中有歷史相對運動狀態(tài)數(shù)據(jù)及初始時刻相對運動狀態(tài)，且主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]均已知的情況。因此在步驟2)及步驟4)中，采用考慮J2攝動的非線性相對運動方程(2.1-2)及(4-2)進(jìn)行相對軌跡預(yù)報。

令編隊飛行的兩航天器分別為主航天器與從航天器，已知進(jìn)行相對軌跡預(yù)報的初始時刻為t₀,主航天器初始軌道要素為E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，其中μ為地心引力常數(shù)，a為半長軸，e為偏心率，i為軌道傾角，Ω為升交點赤經(jīng)，ω為近拱點角距，f為真近點角。從航天器相對主航天器的初始相對運動狀態(tài)為x(t₀)，其中x(t₀)在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示，該坐標(biāo)系原點為主航天器質(zhì)心，x軸沿主航天器地心矢徑方向，z軸沿軌道面法向，y軸與x、z軸構(gòu)成右手坐標(biāo)系。

本實施例根據(jù)歷史相對運動狀態(tài)數(shù)據(jù)、初始相對運動狀態(tài)及主航天器初始軌道根數(shù)就可以進(jìn)行相對軌跡預(yù)報。預(yù)報過程為解析的考慮J2攝動的非線性相對運動方程，通過最小二乘估計新的初始條件能最小化預(yù)報誤差，預(yù)報精度得到改進(jìn)；且基于最小二乘的迭代過程均為解析計算，能夠保證迭代過程快速收斂，計算效率高。對同樣的歷史數(shù)據(jù)及初始條件，本實施例所得精度比實施例二要高。

參見圖2所示相對軌跡預(yù)報位置誤差對比可知，和現(xiàn)有技術(shù)采用J2非線性方程直接預(yù)報相比，本實施例(采用J2非線性方程的改進(jìn))的相對距離誤差x更??；參見圖3所示相對軌跡預(yù)報速度誤差對比可知，和現(xiàn)有技術(shù)采用J2非線性方程直接預(yù)報相比，本實施例(采用J2非線性方程的改進(jìn))的相對速度誤差x更小。

實施例四：

本實施例用于初始條件中沒有歷史相對運動狀態(tài)數(shù)據(jù)，僅知道初始時刻相對運動狀態(tài)，及主航天器的6個軌道根數(shù)E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]的情況。該情況下，直接采用考慮J2攝動的非線性相對運動方程(2.1-2)進(jìn)行相對運動軌跡預(yù)報。

令編隊飛行的兩航天器分別為主航天器與從航天器，已知進(jìn)行相對軌跡預(yù)報的初始時刻為t₀,主航天器初始軌道要素為E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]，其中μ為地心引力常數(shù)，a為半長軸，e為偏心率，i為軌道傾角，Ω為升交點赤經(jīng)，ω為近拱點角距，f為真近點角。從航天器相對主航天器的初始相對運動狀態(tài)為x(t₀)，其中x(t₀)在主航天器當(dāng)?shù)剀壍雷鴺?biāo)系中表示，該坐標(biāo)系原點為主航天器質(zhì)心，x軸沿主航天器地心矢徑方向，z軸沿軌道面法向，y軸與x、z軸構(gòu)成右手坐標(biāo)系。

將初始軌道要素為E(t₀)＝[a(t₀),e(t₀),i(t₀),Ω(t₀),ω(t₀),f(t₀)]及從航天器相對主航天器的初始相對運動狀態(tài)為x(t₀)直接帶入考慮J2攝動的非線性相對運動方程(2.1-2)，即可獲得任意時刻的相對運動狀態(tài)預(yù)報結(jié)果。

本實施例在沒有歷史相對運動狀態(tài)數(shù)據(jù)的情況下就可以進(jìn)行相對軌跡預(yù)報。預(yù)報過程為解析的考慮J2攝動的非線性相對運動方程，計算效率高。對同樣的歷史數(shù)據(jù)及初始條件，本實施例所得精度比實施例三要低。

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