基于自適應算法的月球引力常數和J2系數確定

王 祎，杏建軍，鄭黎明，于 洋

(中南大學航空航天學院，長沙410083)

0 引言

月球引力場的高精度測量對研究月球演化和開展月球探測具有重要意義，因此成為國際深空探測領域研究的熱點。然而由于月球自轉周期和公轉周期同步，導致月球始終只有一個面朝向地球，因此月球背面引力場數據很難用傳統方法測量，這也是高精度月球引力場模型建立所面臨的最大技術難題。Konopliv 等［1］、Matsumoto 等［2］和 Goossens 等［3］國內外學者也對高精度月球引力場模型的建立做了眾多的研究:Konopliv采用了美國的GRAIL［4］的月球引力場測量數據恢復月球引力場，月球引力常數精度達到4.90280031 ×1012±4.4 ×105m3/s2，其中J2項攝動系數精度為 9.0880835 ×10－5±1.4 ×10－9;Matsumoto和 Goossens則是分析了日本的SELENE［5］計劃的月球引力場的測量數據后得到精度為4.90280080×1012±9×105m3/s2的月球引力常數;利用環月飛行的月球探測器之間衛衛跟蹤測量輔助地面深空探測網進行月球引力場確定，是現今直接測量月球背面引力場的唯一方法，日本SELENE計劃和美國GRAIL任務分別應用這種方法進行了月球背面引力場的測量，有效地提高了月球引力場測量的精度。

Psiaki［6］提出了一種基于星間相對位置測量信息自主確定月球引力場的方法。該方法利用兩顆環月編隊飛行的衛星構成一個獨立自主的測量系統，環月編隊飛行的衛星在不依賴地面深空探測網支持的條件下，利用編隊衛星星間相對位置的測量數據，采用擴展卡爾曼濾波器(Extend Kalman Filter，EKF)濾波算法，自主確定月球衛星軌道，同時修正月球引力場參數，得到更高精度月球引力場模型。

考慮到實時導航要求，Psiaki［6］采用實時性較強的EKF算法，且EKF算法在處理白噪聲激勵下的平穩或非平穩向量的狀態估計問題時能夠得到最優的結果［7］。然而在實際應用中，EKF的動力學模型和觀測模型不可能準確無誤，系統噪聲和測量噪聲特性不能提前精確確定，導致狀態參數估計精度的降低，甚至濾波器的發散［8］。

針對EKF算法中的系統噪聲或者測量噪聲難以精確確定導致濾波器精度降低的問題，本文提出一種自適應擴展卡爾曼濾波器(Adaptive Extend Kalman Filter，AEKF)。該算法將改進的Sage-Husa次優無偏極大估值器［9］與擴展卡爾曼濾波器相結合，在計算過程中對未知噪聲的統計特性進行實時估計和修正，并利用測量值來判斷濾波是否發散，當濾波發散時通過引入自適應衰減因子來修正預測狀態協方差矩陣，進而調整測量方程協方差矩陣和濾波增益矩陣，抑制由于干擾所引起的濾波發散，提高濾波器的精度［10］。

1 自適應濾波算法

1. 1 動力學模型與測量模型

雙星環月編隊飛行動力學方程和測量方程可表示為:

式中:Xk表示目標在k時刻系統的狀態向量，包括參考衛星在慣性空間的位置速度矢量、伴隨衛星在慣性空間的位置速度矢量、月球引力場常數和月球J2項攝動系數14個量;Zk為k時刻編隊衛星和參考衛星相對位置矢量在慣性空間的測量值;fk－1(Xk－1)為狀態轉移函數;Hk為測量轉移矩陣;假設動力學方程的噪聲wk為期望值為零，方差為Q的高斯白噪聲序列，測量噪聲vk是期望值為零，在區間(－1，1)均勻分布的噪聲序列，測量噪聲特性不滿足標準EKF的條件。

1. 2 EKF 濾波算法

由式(1)可得預測狀態值和預測測量值分別為:

計算狀態轉移矩陣Φk/k－1為:

預測狀態的協方差矩陣和濾波增益矩陣為［11－12］:

估計狀態值和協方差矩陣更新為:

1. 3 AEKF 濾波算法

在傳統的EKF算法中，系統噪聲和測量噪聲為白噪聲，并且方差Q和R都是已知的，但是在實際情況中噪聲不一定為白噪聲，且統計特性Q和R未知，不準確的噪聲統計特性會影響狀態估計的精度和穩定性，甚至引起濾波發散。針對這一問題，Sage與Husa(1969)提出了對Q或R進行修正估計的方案，同時張常云等［13］在對 Sage－Husa的 AEKF的后續研究中證明了該算法不能同時修正Q與R，并分析了該法導致分散的原因。

本文提出了一種改進的AEKF算法，該方法在濾波過程中，利用改進的Sage-Husa估計器估計未知測量噪聲的統計特性，通過已知的Q可以準確地估計出R，并對濾波發散情況進行判斷和抑制，有效提高了濾波數值的抗干擾性，減小了狀態估計誤差。

由(1)式與(3)式可得:

根據文獻［14］，對(9)式兩邊同時取協方差可以得到預測誤差向量的協方差矩陣:

Ck的計算通常采用移動開窗估計法，則Ck的估計值可以為:

式中:N為窗口函數大小，從式(12)中可以看出，若N值過大，Ck不能反映觀測噪聲的瞬態特性，若N值選取過小，Ck易受當前預測殘差向量的影響，導致估計值的平穩性較差。

1. 4 濾波發散的抑制

由于傳統的EKF常常會出現發散的情況，本文提出一種基于協方差匹配判據的方法對濾波發散趨勢進行判斷，由

式中:tr(^Rk)表示矩陣 ^Rk求跡

若該式不成立，修正Pk(－)可得:

當滿足式(13)的條件時，就采用前面所述Sage-Husa算法;而當不滿足式(13)的條件時，就采用衰減因子計算公式［15－16］，得到一自適應加權系數 Ωk，再用它對Pk(－)進行修正，其中，Ωk由下式確定

式中:ρ＞0為衰減系數，能進一步提高濾波器快速跟蹤能力，通常取值在11左右，其值越小，則k時刻以前的信息所占的比例就越小，當前殘差向量的影響就越突出，該方法有很強的關于突變狀態的跟蹤能力，并且在濾波達到穩定時，仍保持對于緩變以及突變狀態的跟蹤能力。

當采用自適應濾波算法時，由于算法實時對測量數據的方差R進行估計，并且將方差R應用到EKF增益式(6)的計算中。這樣當測量數據誤差較大時，EKF增益的值就會降低，測量數據在狀態估計中的權重降低，導致狀態估計精度降低速度減慢，提高算法的收斂性。

1. 5 AEKF濾波算法實現

針對式(1)描述的自主月球引力場確定問題，AEKF具體步驟如下(如圖1):

1)狀態初始條件為

2)預測更新，對于給定的 Qk－1，Pk－1，根據式(2－4)，求狀態預測，測量預測和預測協方差。

3)由式(11－12)估計測量噪聲統計特性。

4)判斷是否符合式(13)的條件，如不符合按照式(14－15)修正Pk(－)，符合就不需要修正Pk(－)直接進入下一步。

5)測量更新，根據式(6－8)得到量測更新和狀態估計值。

2 仿真校驗

為了校驗本文AEKF的效果，在測量噪聲統計特性未知的情況下，分別采用標準的EKF算法和AEKF算法進行仿真，并且與蒙特卡洛仿真(Monte Carlo Algorithm，MCL)的結果進行比較。

2. 1 坐標系的選取

根據文獻［17］，建立本文所采用的坐標系:以OE－XEYEZE表示J2000地球慣性坐標系，原點為地心，基本面為地球平赤道面，XE軸指向歷元J2000春分點，ZE軸指向地球北極，根據右手螺旋定律由XE軸和ZE軸確定YE軸方向，然后以月心為原點，將J2000地球慣性坐標系平移至月心，得到月心J2000慣性坐標系OL－XYZ;以OS-X'Y'Z'表示軌道坐標系，以參考衛星的質心為坐標原點，X'軸由月心指向參考衛星的質心，Z'軸指向參考衛星軌道面法向方向，Y'軸由右手法則，垂直于X'軸和Z'軸所組成的平面。本文用 XR，YR，ZR和 XF，YF，ZF分別表示參考衛星和伴隨衛星在OL－XYZ坐標系中的三個方向矢量。

2. 2 仿真場景

本文基于Matlab7.0軟件，在配置為因特爾奔騰的處理器(雙核2.8 GHz)與兩根海力士DDR3內存(2 GB)的計算機上對近月軌道飛行衛星編隊進行仿真，考慮到月球 J2項攝動的影響，以月球LP165P模型為參考對象，通過AEKF算法實現伴隨衛星和參考衛星在OL－XYZ坐標系中精確定軌，同時估計月球引力場常數和J2項攝動常數。本文仿真采用的伴隨衛星和參考衛星的具體初始軌道根數如表1，2所示:

表1 參考衛星初始軌道根數Table 1 Initial orbit elements of reference satellite

表2 伴隨衛星初始軌道根數Table 2 Initial orbit elements of follow satellite

根據文獻［18］，假設系統狀態方程噪聲協方差矩陣 Qn=diag(1 × 10－6，1 × 10－6，1 × 10－6，1 ×10－10，1 × 10－10，1 × 10－10，1 × 10－6，1 × 10－6，1 ×10－6，1 × 10－10，1 × 10－10，1 × 10－10，1 × 104，1 ×10－36)。將表1中的初始軌道根數轉化成OL－XYZ坐標系中的位置速度和LP165P中的月球引力常數與J2攝動常數，分別由式(1)得到參考衛星和伴隨衛星在OL－XYZ坐標系中的絕對位置和速度，然后根據式(1)，得到編隊衛星相對位置的測量值。

2. 3 仿真結果

分別對AEKF(衰減系數ρ=11窗口函數N=14)、先驗噪聲為 Rn=diag(1/3，1/3，1/3)的不同初始狀態值的50組傳統EKF和MCL仿真結果進行比較。

AEKF采用的初始狀態誤差為 ΔX=［9.76×102，7.56 ×102，9.03 × 102，9.82 × 10－1，1.02，8.95×10－1，9.48 × 102，1.05 × 103，1.06 × 103，1.06，8.51 ×10－1，8.92 ×10－1，5.09 ×108，8.71 ×10－7］，初始協方差矩陣都為 P0=diag(9.53×105，5.71×105，8.15 ×105，9.64 ×10－1，1.04，8.01 ×10－1，8.98×105，1.11 × 106，1.12 × 106，7.24 × 10－1，7.95 ×10－1，9.93 ×10－1，2.59 ×1017，7.59 ×10－13)。50 組傳統EKF參考衛星三個方向的初始位置估計誤差分布符合期望值是0，方差分別是9.53×105，5.71×105，8.15×105的高斯分布，三個方向的初始速度估計誤差分布符合期望值是0，方差分別是9.64×10－1，1.04，8.01 ×10－1的高斯分布;伴隨衛星三個方向的初始位置估計誤差分布符合期望值是0，方差分別是 8.98 ×105，1.11 ×106，1.12 ×106的高斯分布，三個方向的初始速度估計誤差分布符合期望值是 0，方差分別是 7.24 × 10－1，7.95 × 10－1，9.93×10－1的高斯分布;月球引力常數初始估計誤差分布符合期望值是0，方差是2.59×1017的高斯分布，J2項攝動力常數初始估計誤差分布符合期望值是0，方差是 7.59×10－13的高斯分布;50 組 Monte Carlo仿真的初始狀態誤差估計分布與傳統EKF初始狀態估計誤差分布一樣。

通過對環月編隊衛星飛行7430 s的相對測量數據進行處理(測量更新間隔為5 s)，表3～5給出了三種方法估計最終時刻編隊衛星位置速度，月球引力常數與J2項攝動力常數估計的一倍標準差。

表3 位置標準差比較Table 3 The standard deviation comparison of position in three-directions

表4 速度標準差比較Table 4 The standard deviation comparison of velocity in three-directions

表5 月球引力常數與J2項攝動力常數標準差比較Table 5 The standard deviation comparison of the lunar gravitational constant and the coefficient of the J2 perturbation

由表3～5可以看出:AEKF算法是收斂的，并且能有效地提高編隊衛星軌道確定的精度。在7430 s的時間內，將初始X軸方向的位置誤差的標準差從976 m(1σ)提高到 51.1 m(1σ)，并且使球引力常數精度和月球J2項攝動系數精度分別達到了4.43 ×108m3/s2(1σ)和 9.12 ×10－8(1σ);MCL仿真與AEKF算法得到的結果基本一致，校驗了AEKF方法的正確性和有效性;傳統EKF是發散的，說明與傳統EKF假設不符合的測量噪聲使得EKF不能有效地收斂，而AEKF算法由于對測量噪聲的統計特性進行了實時修正，因此有效減少干擾影響，抑制發散從而提高了參數估計精度;EKF與AEKF方法運行一次的計算時間分別為 206.45 s和211.06 s，AEKF運行一次的時間較長，這是由于AEKF算法利用濾波器殘差對測量噪聲的統計特性進行了自動修正，增大了計算量從而增加了計算時間;該方法得到的月球引力常數精度為4.902801056 ×1012±4.43×108m3/s2(1σ)，歸一化的J2項攝動常數為 －9.08 ×10－5±9.12 ×10－8(1σ)，雖然與 Alex S.Konopliv 等［1－3］人的結果相比得到的引力場模型精度較低，但是由于不需要例如深空探測網(美國GRAIL計劃)，中繼衛星(日本SELENE計劃)等其他設施的輔助，能自主確定月球引力場，因此大大降低了任務實施的成本。

圖2 均勻分布噪聲下的EKF算法位置估計誤差Fig.2 The estimation position error using EKF under uniformly distributed random measurement noisy

圖3 高斯分布噪聲下的EKF算法位置估計誤差Fig.3 The estimation position error using EKF under gaussian distributed random measurement noisy

對傳統EKF進一步分析，圖2是在測量噪聲是均勻分布噪聲，同時改變編隊衛星的初始估計誤差(50組)情況下，采用EKF算法時參考衛星在OLXYZ坐標系下X軸方向的位置估計誤差的曲線圖;圖3是在測量噪聲是高斯分布噪聲，同時改變編隊衛星的初始估計誤差(50組)情況下，采用EKF算法時參考衛星在OL－XYZ坐標系下X軸方向的位置估計誤差的曲線圖。比較圖2與圖3，可以看出，與EKF假設相同的測量噪聲(高斯白噪聲)，并且噪聲統計特性準確知道的情況下，傳統EKF運行的很好，結果可以快速收斂，但當測量噪聲與傳統EKF假設不相同時(非高斯白噪聲)，傳統EKF濾波精度不高，甚至引起濾波發散。

圖4是在測量噪聲是均勻分布噪聲，同時改變編隊衛星的初始估計誤差(50組)情況下，采用AEKF算法時參考衛星在OL－XYZ坐標系下X軸方向的位置估計誤差的曲線圖。比較圖2與圖4可以看出:噪聲不是高斯白噪聲時，傳統EKF算法的狀態估計誤差比AEKF算法的狀態估計誤差大，甚至引起濾波發散。這是因為傳統EKF在濾波過程中一般假設測量噪聲是高斯白噪聲，但實際中測量噪聲不一定是高斯白噪聲;AEKF算法中測量噪聲的統計特性會隨著濾波器殘差的大小自動修正，因此估計的精度和抗干擾性比較好。

3 結論

針對系統噪聲或者測量噪聲未知情況下月球的引力場確定問題，提出了一種新的基于殘差信息調整方差的自適應EKF算法，該算法基于編隊衛星相對位置測量能自主確定月球的引力場，同時能對可能出現的濾波算法發散情況進行抑制，減少外界干擾對月球引力場測量的影響。文中給出的仿真實例校驗了該方法的有效性，仿真結果表明本文提出的AEKF算法濾波效果要明顯優于傳統EKF算法，精度高，抗干擾性強，月球引力常數和月球J2項攝動常數測量精度分別達到了4.43×108m3/s2(1σ)和9.12 ×10－8(1σ)。

［1］ Konopliv A S，Park R S，Yuan D N，et al.The JPL lunar gravity field to spherical harmonic degree 660 from the GRAIL primary mission［J］.Journal of Geophysical Research，2013，118:1415－1434.

［2］ Matsumoto K，Hanada H，Namiki N，et al.A simulation study for anticipated accuracy of lunar gravity field model by Selene tracking data［J］.Advances in Space Research，2008，42(8):331－336.

［3］ Goossens S，Matsumoto K，Liu Q，et al.Lunar gravity field determination using SELENE same-beam differential VLBI tracking data［J］.Journal of Geodesy，2011，85(4):205 －228.

［4］ Yan J G，Baur O，Fei L，et al.Long-wavelength lunar gravity field recovery from simulated orbit and inter-satellite tracking data［J］.Advances in Space Research，2013，52(11):1919 －1928.

［5］ Tanaka S，Shiraishi H，Kato M，et al.The science objectives of the SELENE-Ⅱ mission as the post SELENE mission［J］.Advances in Space Research，2008，42(2):394－401.

［6］ Mark L P.Absolute orbit and gravity determination using relative position measurements between two satellites［J］.Journal of Guidance，Control，and Dynamics，2011，34(5):1285 －1297.

［7］ 尚琳，劉國華，張銳，等.基于BP神經網絡的自主定軌自適應Kalman濾波算法［J］.宇航學報，2013，34(7):926－931.［Shang Lin，Liu Guo-hua，Zhang Rui，et al.An adaptive Kalman filtering algorithm for autonomous orbit determination based-on BP neural network［J］.Journal of Astronautics，2013，34(7):926 －931.］

［8］ 席燕輝，葉志成，彭輝.一種基于自適應粒子濾波的多層感知器學習算法［J］.中南大學學報(自然科學版)，2013，44(4):1397 －1402.［Xi Yan-hui，Ye Zhi-cheng，Peng Hui.An algorithm for MLPs training based on adaptive particle filter［J］.Journal of Central South University(Science and Technology)，2013，44(4):1397 －1402.］

［9］ Sage A P，Husa G W.Adaptive filtering with unknown prior statistics［C］.Joint Automatic Control Conference，Boulder Colorado，USA，August，1969.

［10］ 石勇，韓崇昭.自適應UKF算法在目標跟蹤中的應用［J］.自動化學報，2011，37(6):755－759.［Shi Yong，Han Chong-zhao.Adaptive UKF method with applications to target tracking［J］.Acta Automatica Sinica，2011，37(6):755 －759.］

［11］ 陳宇波，宋迎春.非高斯噪聲下的車載GPS信號定位算法［J］.中南大學學報(自然科學版)，2010，41(4):1462－1466.［Chen Yu-bo，Song Ying-chun.A Bayes filter algorithm with non-Gaussian noises based on location of vehicular GPS［J］.Journal of Central South University(Science and Technology)，2010，41(4):1462 －1466.］

［12］ 傅薈璇，王宇超，孫楓.融合Kalman濾波的自適應帶寬Mean Shift算法［J］.中南大學學報(自然科學版)，2011，42(1):784 － 788.［Fu Hui-xuan，Wang Yu-chao，Sun Feng.Tracking algorithm based on Kalman and mean shift with adaptive bandwidth［J］.Journal of Central South University(Science and Technology)，2011，42(1):784 －788.］

［13］ 張常云.自適應濾波方法研究［J］.航空學報，1998，19(7):96－99.［Zhang Chang-yun.Approach to adaptive filtering algorithm［J］.Acta Aeronautica et Astronautica Sinica，1998，19(7):96 －99.］

［14］ 王志忠，朱建軍.非線性模型中方差和協方差分量的估計［J］.測繪學報，2005，34(4):288－293.［Wang Zhi-zhong，Zhu Jian-jun.Estimation of variance and covariance components in the nonlinear model［J］.Acta Geodaetica et Cartographica Sinica，2005，34(4):288 －293.］

［15］ Yang Y X，Gao W G.An optimal adaptive Kalman filter［J］.Journal of Geodesy，2006，80(4):177－183.

［16］ Hide C，Moore T，Smith M.Adaptive Kalman filtering for lowcost INS/GPS［J］.The Journal of Navigation，2003，56(1):143－152.

［17］ 李立濤，楊滌，崔祜濤.奔月軌道的一種求解方法［J］.宇航學報，2003，24(2):150－155.［Li Li-tao，Yang Di，Cui Hutao.One method of calculating cislunar transfer trajecotry［J］.Journal of Astronautics，2003，24(2):150 －155.］

［18］ 張艷.基于星間觀測的星座自主導航方法研究［D］.長沙:國防科學技術大學，2005.［Zhang Yan.Study on autonomous navigation of constellation using inter-satellite measurement［D］.Changsha:National University of Defense Technology，2005.］