[發明專利]一種基于重力場梯度不變量的軌道要素估計方法

有效

• 申請號: 201710024671.5
• 申請日: 2017-01-13
• 公開（公告）號: CN107065025B
• 公開（公告）日: 2019-04-23
• 發明（設計）人: 陳培;孫秀聰;張鍵
• 申請（專利權）人: 北京航空航天大學
• 主分類號: G01V7/00
• 分類號: G01V7/00
• 代理公司: 暫無信息
• 代理人: 暫無信息

1.一種基于重力場梯度不變量的軌道要素估計方法，其特征在于：其步驟如下：

步驟一：準備工作

1)地球引力勢函數

在初步研究航天器運動規律時，把地球看作理想的球體，它對航天器產生的引力指向地心，其大小F與航天器的質量m成正比，而與地心至航天器的距離r的平方成反比：

上式中，F為航天器所受地球引力，G為地球引力常量，m為航天器質量，r為地心至航天器的幾何距離，F除以m就是引力加速度g，把兩個常數合并μ＝GM，則得到下式：

上式中，μ成為地球引力常數；

此情況下地球引力場成為中心引力場，把引力加速度表示成引力勢的梯度，則中心引力場的引力勢函數是：

實際的地球并不是球對稱的，具有扁平度、梨形和赤道變形，因而引力加速度是距離r、緯度和經度λ的函數，而r，λ是在地球固連坐標系S_e中位置的球坐標，為了嚴格地而且方便地描述地球引力加速度的分布，把它表示成地球引力勢函數U的梯度：

上式中grad(*)為求梯度算子；

以下給出，在只考慮地球非球形J₂項影響下的引力場表達式：

上式中Re＝6378.1363km，為地球赤道半徑；

2)開普勒軌道根數

在慣性空間觀察航天器軌道，需要定義軌道要素；軌道要素，又稱軌道根數，是決定開普勒軌道的運動特征的一組常數；在考慮軌道攝動后，軌道根數不再是常數，可以作為狀態變量；

軌道根數一共有6個，表示及含義如下：

a：橢圓軌道半場軸；

e：橢圓軌道偏心率；

Ω：升交點赤經；航天器軌道由南向北穿越赤道的點稱為升交點B；赤道平面內，由春分點向東轉到升交點B的角度稱為升交點赤經，有效范圍是0至360度；

i：軌道傾角；軌道平面與赤道平面的夾角，有效范圍是0至180度；

t_p：近地點時刻；航天器通過近地點的時刻；

ω：近地點幅角；軌道平面內，有升交點轉到近地點的角度，有效范圍是0至360度；

他們的作用分別是：軌道半長軸a和偏心率e確定軌道的大小和形狀；升交點赤經Ω和軌道傾角i確定軌道平面在慣性空間的取向；近地點幅角ω確定拱線在軌道平面內的位置；紀元時刻的平近點角M(t₀)確定時間的起點；

3)軌道攝動方程

在攝動力作用下航天器的軌道不是開普勒軌道，其軌道要素是隨時間變化的：a(t)，e(t)，Ω(t)，i(t)，t_p，ω(t)；

以下給出軌道攝動方程：

上式中，與前文已出現的符號具有相同含義；f_r表示攝動加速度在地心軌道坐標系的徑向分量；f_u表示攝動加速度在地心軌道坐標系的橫向分量；f_h表示攝動加速度在地心軌道坐標系的副法向分量；p為半通徑；

步驟二：理想重力梯度張量在東北天(East-North-Up,ENU)坐標系下的分解和特征值

理想球體重力場模型如公式(3)所示，式中r有如下表示：

上式中，x，y，z分別是衛星位置向量r在地球固連坐標系上的三個分量投影；

重力梯度張量可以看成是引力加速度的梯度，引力加速度又是引力場的梯度，因而直接在ENU坐標系下分解：

上式中，下標E，N，U表示在ENU坐標系下分解的不同元素；

通過式(8)可以方便的得到重力梯度矩陣的三個特征值：

從式(9)可以知，三個特征值兩兩線性相關，只有一個特征量；

步驟三：求解測量歷元J2模型重力場張量在ENU坐標系下的分解和特征值

J2模型重力場張量由式(5)給出，采用求梯度算子求解J2模型重力場張量在ENU坐標系下的分解：

其中，各個分量形式如下：

由式(10)求解矩陣的特征值λ₁,λ₂,λ₃，特征值與r,的關系如下：

可以驗證λ₁，λ₂，λ₃有如下關系：

λ₁+λ₂+λ₃＝0 (13)

因此，三個特征值中只有兩個獨立量；

步驟四：利用J2重力梯度矩陣特征值求解各測量歷元r和

1)計算r初值

使用J2模型求解的重力梯度張量矩陣特征值(12)，和理想球體重力場模型特征值表達式(9)計算r初值；

式(12)中的特征值λ₁和λ₂任選其一皆可，帶入式(9)，反解r則有如下初值表達式：

或

2)計算初值

由式(12)做如下計算：

從式(16)可以得到的表達式：

從式(17)可以得到的另一表達式：

以上的兩種表達式，任選其一便可；

3)利用更新r

利用式(12)中的λ₂表達式更新r，將r₀和帶入反解r³項，表達式如下：

4)利用r₁更新

方法一：

用式(12)中的三個特征值做如下計算：

將r₁和帶入上式(23)，反解上式中的項，如下：

方法二：

使用式(21)和式(22)做如下計算：

將r₁和帶入上式(25)，反解上式中的項，如下：

步驟五：利用r和計算初始軌道要素

1)計算半長軸和偏心率

在所有計算歷元的結果r中搜索最大值r_max和最小值r_min，則有以下關系：

r_a＝r_max,r_p＝r_min (27)

上式中，下標a表示遠地點，下標p表示近地點；

則半長軸：

偏心率：

2)計算軌道傾角

在此假設，通過步驟四所求得的值的正負號已得知，在所有計算歷元的結果中搜索最大值，

軌道傾角：

3)計算近地點幅角和初始真近點角

在2)中搜索到的r_min其對應的緯度記為則近地點幅角ω由下式計算：

初始歷元真近點角θ₀，由下式計算：

上式中，為第一個歷元所計算出的緯度；

步驟六：軌道要素平滑

將初始歷元五個軌道要素作為待估參數，如下：

x＝[a₀ e_x0 e_y0 i₀ u₀] (33)

上式中，下標“0”表示初始歷元參數，e_x,e_y,u表達式如下：

e_x＝ecosω

e_y＝esinω (34)

u＝ω+θ

待估參數x的初值，由步驟五給出；

軌道動力學模型由式(6)描述，半通徑由于J2項引起的攝動力加速度在地心軌道坐標系下的分解如下：

結合式(6)與式(35)，以及由步驟五提供的初值，便構成了已知初值的微分方程組求解問題，則通過數值積分方法可求得任意k歷元時刻的五個軌道要素，x_k＝[a_k e_xk e_yk i_k u_k]，則k歷元時刻的觀測方程如下：

上式中，為步驟四中所求得的各歷元衛星到地心的幾何距離與緯度，Δr,表示計算誤差，下標“k”表示第k個歷元；

簡單地假設Δr,是一個高斯噪聲，則k歷元測量協方差矩陣如下：

上式中，σ_λ為重力梯度不變量誤差，其誤差大小要根據測量儀器的測量噪聲與J2模型誤差大小選取；

將各個歷元的測量值z放置在一起寫成如下形式：

上式中，N為測量歷元的個數；各個歷元測量協方差矩陣寫成以下形式：

采用準西格瑪最小二乘批處理濾波估計狀態參數，其估計式如下：

上式中，下標j表示迭代次數，H_j表示利用計算出的測量量z的計算值與敏度矩陣；

式(40)中，求解H陣難度較大，將Z寫成關于x＝[a₀ e_x0 e_y0 i₀ u₀]的隱函數表達，如下：

Z＝f(a₀ e_x0 e_y0 i₀ u₀) (41)

通過普通解析法求偏導計算敏度矩陣H,計算量非常巨大，因此可以通過以下方法求解；

在第j次迭代完畢，得到狀態參數估計值則進行如下采樣：

上式中，L為x的維度，P₀為先驗協方差矩陣，表示P₀的第i列平方根；先驗協方差矩陣由下式給出：

將式(42)通過式(6)積分，和式(36)傳遞，計算測量值如下：

χ_i,j→γ_i,j,i＝0,...,L (44)

則敏度矩陣H可以近似通過下式得出：

上式中，表示P₀第i個對角線元素的平方根,i＝1,2,3,4,5。