李 燁，郁 豐，王群仰

（1. 航天恒星科技有限公司，北京 100086；2. 南京航空航天大學 航天學院，南京 210016）

立方星剩磁在軌辨識與主動補償技術

李 燁1，郁 豐2，王群仰1

（1. 航天恒星科技有限公司，北京 100086；2. 南京航空航天大學 航天學院，南京 210016）

立方星的姿態(tài)測量與控制系統(tǒng)常采用磁測磁控結合偏置動量輪的方案，整星剩磁干擾力矩是影響姿態(tài)控制精度的重要因素之一。提出了一種利用磁強計實現(xiàn)剩磁矩在軌辨識與利用磁力矩器實現(xiàn)剩磁矩主動補償?shù)男路桨福夯诖艔娪嬢敵龊托l(wèi)星姿態(tài)動力學建立了剩磁矩在軌辨識模型，并利用采樣濾波器（UKF）提高單磁強計條件下的辨識效果；把控制對象簡化成線性定常系統(tǒng)，分析了剩磁干擾力矩對姿態(tài)的影響數(shù)學模型，并針對磁力矩器和磁強計分時工作的特點，基于疊加性原理提出了基于角速度的剩磁矩主動補償算法。仿真研究表明，在1000 s內剩磁矩在軌辨識精度為0.001 A·m2量級，主動補償后，偏航角、滾動角與俯仰角控制誤差分別從4.3°、4.6°與2.1°均減少至0.4°以內。提出的方法為類似配置衛(wèi)星減少剩磁干擾力矩的影響提供了一種新思路。

剩磁；姿態(tài)測量與控制系統(tǒng)；在軌辨識；剩磁補償

立方星是一種質量與體積都極小的低成本衛(wèi)星，大量采用商用器件，具有很高的功能密度，在國際上廣泛應用于航天科學研究與大學教育[1-2]。隨著技術的進步，近幾年立方星已經呈現(xiàn)出實用化的趨勢，如行星實驗室公司研制的鴿群對地成像星座。立方星由于質量與功耗極其苛刻的限制，姿態(tài)測量與控制系統(tǒng)常選用能易于小型化的磁強計、偏置動量輪和磁力矩器作為星上的敏感器和執(zhí)行機構，甚至將磁強計復用成定軌敏感器，以提高衛(wèi)星的功能密度[3-4]。

但是衛(wèi)星在軌工作時，整星內部會形成一個基本穩(wěn)定的磁場，稱之為剩磁。剩磁不但對某些載荷有影響[5]，而且剩磁矩與地磁場相互作用形成剩磁干擾力矩，會嚴重影響姿態(tài)控制的精度。為了減少剩磁對星體擾動的影響，通常在地面開展整星磁補償實驗，即首先利用多臺標準磁強計測量并反演出衛(wèi)星的剩磁矩，然后用小型永磁體安裝在星體內部抵消剩磁矩，最后再次測量確認磁補償?shù)男ЧＮ墨I[6]利用該方法在轉子上安裝一個與定子永磁體方向相反、大小相同的永磁體，降低了飛輪的剩磁矩，但是該方法會增加衛(wèi)星研制過程中的各種成本，并且剩磁控制水平存在限度，不是非常適合立方星。

在軌磁補償是一種利用磁力矩器來抵消星體剩磁矩的技術方案，能夠避免地面復雜的標定和補償過程，在一些場合具有明顯優(yōu)勢。SNAP-1衛(wèi)星在入軌后才發(fā)現(xiàn)星體存在較大剩磁，嚴重影響了姿態(tài)確定與控制系統(tǒng)的正常工作。地面設計了一個二步迭代的方案用于估計星體剩磁矩[7]，它首先利用星上估計的角速度作為觀測值估計了剩磁矩，但是該角速度測量值也受剩磁矩的影響，所以利用估計出的剩磁矩再對角速度估計做修正，獲得更精確的角速度估計值作為下次計算的輸入。實際操作中共進行了3次迭代，使用了3天中的3組數(shù)據(jù)，估計精度大約在0.01 A·m2量級。利用磁力矩器補償剩磁矩時采用直接扣除方法，沒有考慮磁力矩器與磁強計分時工作帶來的問題，補償后的姿態(tài)控制精度提高至5°~7.5°。

本文針對磁測磁控結合偏置動量的姿態(tài)測量與控制方案，提出了一種僅利用磁強計和磁力矩器實現(xiàn)整星剩磁矩在軌辨識與主動補償?shù)男路桨?。它基于磁強計和衛(wèi)星姿態(tài)動力學建立了剩磁矩在軌辨識模型，并利用 UKF算法提高單磁強計條件下的辨識精度和速度；針對磁強計和磁力矩器分時工作的特點，把控制對象簡化成線性定常系統(tǒng)，分析剩磁干擾力矩對姿態(tài)影響數(shù)學模型，基于疊加性原理提出了基于角速度的剩磁矩主動補償算法。研究表明，本文方法具有所需器件少、不依賴地面輔助等優(yōu)勢，并且剩磁矩辨識精度高、速度快，主動剩磁補償使得控制誤差能小于0.4°，綜合效果良好。

1 剩磁矩對控制系統(tǒng)的影響分析

某型立方衛(wèi)星的姿態(tài)確定與控制系統(tǒng)采用磁強計作為姿態(tài)敏感器，以偏置動量輪和三軸磁力矩器為控制執(zhí)行機構。其中，-Y軸安裝的偏置動量輪常速旋轉，實現(xiàn)星體的被動穩(wěn)定，X和Z軸磁力矩器控制衛(wèi)星俯仰通道的姿態(tài)，Y軸磁力矩器實現(xiàn)章進動控制。

值得指出的是，由于立方星體積小，磁強計與磁力矩器之間的安裝距離很近，磁力矩器工作時產生的磁矩會影響磁強計測量地磁場，所以星上的磁強計和磁力矩器是分時工作的，即在一個控制周期內的一小段時間中，磁力矩器不工作，此時磁強計完成對當?shù)氐卮艌龅臏y量。

式中，m是衛(wèi)星剩磁矩，B是地磁場在星體坐標系中的磁感應強度。

該剩磁干擾矩會對衛(wèi)星的姿態(tài)控制產生較為嚴重的擾動，在立方星上表現(xiàn)得更明顯，主要原因是立方星動量輪的質量和尺寸受到嚴格限制，所以動量輪的角動量普遍偏小，抗干擾能力更差，嚴重影響了立方星姿態(tài)控制精度的提高。

因為立方體衛(wèi)星部件不可避免地使用磁性元件，并且星上電路中的電流也會激發(fā)磁場，所以衛(wèi)星在軌運行時，立方體衛(wèi)星本體會存在一個大小、方向未知的剩余磁矩[8]。如果星上工作模式不變，該剩磁場一般較穩(wěn)定，所以微小衛(wèi)星除受自身的控制力矩作用外，本體的剩磁矩與地磁場相互作用的形成剩磁干擾力矩[9]，

2 剩磁矩在軌辨識技術

為有效減小星體剩磁矩對姿態(tài)控制的擾動影響，需要通過估計算法將星體剩磁矩在軌辨識出來。顯然，剩磁矩最終影響的是衛(wèi)星姿態(tài)，所以本文建立了一種與姿態(tài)耦合的剩磁矩在軌估計算法。

衛(wèi)星的姿態(tài)動力學方程為

磁強計經過地面標定與在軌標定后，可以認為得到較準確的測量值[10]，所以本文選用的磁強計測量模型為

式中，Bc是磁強計輸出，B是衛(wèi)星本體坐標系中地磁場強度的真實值，Bε是磁強計的測量噪聲。

將式(1)代入式(2)中，并將式(3)代入式中，則考慮了剩磁干擾的姿態(tài)動力學方程為

式中，m×Bε表示是由磁強計測量噪聲引起的建模不確定度，Tε+m×Bε是系統(tǒng)噪聲。剩磁矩在短時間內變化很小，所以建模成

衛(wèi)星的姿態(tài)用四元數(shù)q表示，四元數(shù)微分方程為

另一方面，在衛(wèi)星本體坐標系中的磁感應強度B和地心慣性坐標系中的磁感應強度Bi之間的關系為

式中：Pk-1是k-1時的狀態(tài)協(xié)方差陣；R=BεB是量測噪聲陣；Q是系統(tǒng)噪聲陣，計算式為

由于式(10)中的m是待估量，在實際濾波中可以采用m當前的估計值來代替，并考慮到Tε和Bε是兩個獨立的隨機分布，則Q的具體計算式為

式中：n是狀態(tài)量的維數(shù)；λ=α2(n+k)-n，是采樣比例因子，其中的α來決定采樣點的分布，通常取為一個很小的正值，k≥0。

濾波器的時間更新為

濾波器的量測更新為

3 剩磁矩在軌主動補償技術

本節(jié)主要論述利用磁力矩器補償剩磁矩對星體擾動的方法。由于星上的磁力矩器和磁強計采用分時工作方案，所以在一個控制周期內，磁力矩器僅工作部分時間，但是衛(wèi)星的剩磁矩是和地磁場每時每刻相互作用的，所以磁力矩器不能簡單地施加一個與剩磁矩大小相等的反向補償磁矩。把控制周期的時長記為t，在一個控制周期內磁力矩器的有效時間為t?。

一般來說，衛(wèi)星在三軸對地穩(wěn)定工作模式時，有效載荷正常工作，所以該模式下對衛(wèi)星姿態(tài)控制精度的要求較高，所以剩磁矩在軌補償算法主要針對三軸對地穩(wěn)定模式來設計。此時偏置動量衛(wèi)星的姿態(tài)動力學解耦成俯仰通道和滾動偏航通道，分別為

式中，ωx、ωy和ωz是衛(wèi)星角速度的三軸分量，φ和Ψ是衛(wèi)星的滾動角和偏航角，ω0是軌道角速率，Ix、Iy和Iz是衛(wèi)星的三個主轉動慣量，Tx、Ty和Tz是控制力矩的三軸分量。

從上述兩個微分方程的形式上看，衛(wèi)星三軸對地穩(wěn)定時，俯仰通道和滾動偏航通道的姿態(tài)動力學方程是線性定常的，所以該系統(tǒng)滿足疊加性原理。

因此在分析一個控制周期內的剩磁矩補償方法時，可以不用考慮控制系統(tǒng)正常的控制過程，所以首先單獨考慮考慮剩磁干擾力矩對衛(wèi)星姿態(tài)與角速度的擾動情況。利用第2節(jié)方法估計出的剩余磁矩計算剩磁干擾力矩：

因此在一個控制周期末端，俯仰通道產生的擾動角速度為

俯仰通道產生的擾動角度為

將式(24)離散化，滾動偏航通道在一個控制周期內產生的擾動姿態(tài)及其變化率為

則俯仰通道的補償磁力矩為

在滾動偏航通道有

則滾動偏航通道的補償磁力矩為

計算出一個周期內的補償力矩后，利用測量得到的磁場強度 Bc，計算出對應補償磁矩。記補償磁力矩為，則補償磁矩為

4 仿真與分析

利用數(shù)學仿真來驗證上述方法的有效性。設定3U立方星運行在高度為 500 km的圓軌道上，軌道傾角42°；三個主軸慣量為[0.0375 0.0375 0.0075] kg·m2，整星質量 4.5 kg；衛(wèi)星在三軸對地穩(wěn)定時的角動量指向軌道法向，大小為0.004 Nms；磁力矩器單軸的最大控制磁矩為0.18 A·m2；假定磁強計經過均值處理后的測量噪聲為 10 nT(1σ)；假定磁強計經過標定后殘留的零偏為[100 -8060]nT ；大氣密度為 6 ×10-13kg·m3，剩余磁矩矢量設定為[0.01320.01320.0132] A·m2，即整星剩磁矩的模為0.023 A·m2。衛(wèi)星控制系統(tǒng)的控制周期是 2 s，一個控制周期內磁力矩器的有效作用時間是1.6 s，剩余的0.4 s是磁力矩器的休眠時間。

讓衛(wèi)星原有的控制系統(tǒng)正常工作，并運行于三軸對地穩(wěn)定狀態(tài)，仿真時間為8000 s，得到衛(wèi)星的姿態(tài)、角速度等參數(shù)如圖 2~4所示。在仿真時間內，X、Z兩軸的角速率在 0附近，最大角速率誤差小于 0.05 (°)/s，Y軸的角速率是軌道角速率；在剩磁等環(huán)境干擾力矩作用下，角速度控制誤差存在一定變化，這在姿態(tài)角控制誤差曲線上更為明顯，在仿真時間內，偏航角最大誤差小于4.3°，滾動角誤差小于4.6°，俯仰角誤差小于2.1°。滾動偏航通道的控制誤差要大于俯仰通道，這是因為俯仰通道時刻有較好的控制能力，而滾動偏航通道的控制能力與在軌地磁場的分布有關，所以更容易受到環(huán)境力矩的干擾。從另一方面評價滾動偏航通道的控制效果，繪制角動量與軌道法向的夾角，如圖4所示，該夾角最大誤差約為5.3°。

圖2 角速度Fig.2 Angular velocity

圖3 姿態(tài)角Fig.3 Attitude angle

圖4 角動量與軌道法向的夾角Fig.4 Angle between angular momentum and the normal of orbit

然后，再次運行該控制軟件，同時執(zhí)行剩磁矩在軌辨識算法模塊，運行EKF和UKF兩個濾波器，系統(tǒng)噪聲陣為 diag(1× 10-10, 1×1 0-10, 1× 10-10)Nm2，測量噪聲陣為 diag(400,400,400)nT2，狀態(tài)誤差協(xié)方差陣的初值設為 diag[(5°)2,(5°)2,(5°)2,(0.1°/ s)2,(0.1°/s)2, (0.1°/s)2,(0.01A·m2)2,(0.01A·m2)2,(0.01A·m2)2]，剩余磁矩估計值初值設為 0；UKF設置的其余參數(shù)有α= 0.005， k=0， β=3。EKF和UKF算法對三個軸向的剩磁矩估計誤差如圖5和6所示。

圖5 EKF剩余磁矩估計值Fig.5 Remanence estimation using EKF

由上述兩個仿真圖對比得出，利用單磁強計完成三軸剩磁矩的估計，采用UKF的收斂速度和估計精度都比EKF的效果更佳，EKF算法經過3000 s才完成明顯的收斂過程，并且后面的估計曲線仍有較大幅度的波動；UKF在1000 s內就完成了收斂，估計精度也明顯提高。在1000s時，UKF對剩磁矩的估計值為[0.0122 0.0119 0.0133] A·m2，相應的估計誤差是[-0.0010 -0.0013 0.0010] A·m2；對1000~8000 s內的估計值統(tǒng)計分析，估計均值為[0.0122 0.0124 0.0129] A·m2，估計均方差為[0.00018 0.00016 0.00022] A·m2。

圖6 UKF剩余磁矩估計值Fig.6 Remanence estimation using UKF

最后，再一次運行該控制軟件，執(zhí)行UKF剩余磁矩在軌辨識算法模塊1000 s，利用1000 s時的辨識結果對剩磁矩分別進行角速度補償，將補償值疊加到原控制系統(tǒng)的磁力矩器輸出中。

基于角速度補償算法的結果如圖7~圖9所示。在仿真時間內，X、Z兩軸的角速率在0附近，最大角速率誤差仍小于0.05 (°)/s；Y軸的角速率在1000 s后出現(xiàn)了一個擾動，這是由于剩磁補償加入后產生的。對比圖2和圖7，在3000 s區(qū)域附近，采取主動剩磁補償措施后的Y軸角速率比補償前平穩(wěn)。圖8顯示，主動磁補償技術的效果在姿態(tài)角控制誤差曲線上更為明顯，主動磁補后滾動角、偏航角和俯仰角均逐漸減小，在4000 s后可以認為基本收斂穩(wěn)定，在仿真時間內控制誤差小于0.4°，對比圖3能夠說明基于角速度的剩磁補償方案的有效性。類似地，角動量與軌道法向的夾角如圖9所示，該圖描述了在剩磁補償后，衛(wèi)星在章進動控制作用下滾動偏航通道的收斂過程。

與文獻[5]比較，本文算法取得了更好的效果，其原因有多個方面：1）文獻[5]采用了角速度作為觀測值，而對比圖2和圖7，以及圖3和圖8，可以發(fā)現(xiàn)剩磁干擾力矩很小，對角速度擾動影響較小，但是經過長期作用，姿態(tài)參數(shù)會產生更明顯的偏差，所以本文算法通過磁強計輸出作為觀測值，相當于間接引入了姿態(tài)作為觀測值，更利于估計剩磁矩；2）UKF算法具有更好的估計性能；3）文獻[5]算法采用了二步迭代估計，在角速度估計中不考慮剩磁的影響，在剩磁估計中未考慮角速度誤差，如果迭代次數(shù)不夠多，還是會影響估計精度；4）本論文磁補償算法推導了較精確的數(shù)學模型。為說明本文主動磁補償算法的性能，采用同樣的剩磁矩辨識值，利用文獻[5]中的剩磁矩直接補償方法的控制結果如圖10~12所示。

圖7 角速度Fig.7 Angular velocity

圖8 姿態(tài)角Fig.8 Attitude angle

圖9 角動量與軌道法向的夾角Fig.9 Angle between angular momentum and the normal of orbit

圖10 角速度Fig.10 Angular velocity

圖11 姿態(tài)角Fig.11 Attitude angle

圖12 角動量與軌道法向的夾角Fig.12 Angle between angular momentum and the normal of orbit

比較圖8與圖11以及圖9與圖12的控制效果，兩者的區(qū)別僅僅是剩磁補償方法不一樣，雖然已有文獻的補償方法較原系統(tǒng)也有較明顯的效果，但是姿態(tài)精度比本文方法有明顯下降，補償后的滾動角和偏航角誤差已達到1°左右，角動量與軌道法向的夾角在仿真時間內達到了1.3°，這充分說明了本文提出的主動剩磁補償方案的有效性。

5 結 論

本文提出了一種利用單磁強計實現(xiàn)剩磁矩在軌辨識與利用磁力矩器實現(xiàn)剩磁矩主動補償?shù)男路椒?。本方法所需器件少，不依賴地面輔助，并且剩磁矩辨識精度在0.001 A·m2量級，收斂速度小于1000 s，主動剩磁補償使得控制誤差能小于0.4°，綜合效果良好。所提方法可以作為類似配置衛(wèi)星在軌消除剩磁干擾的一個借鑒和參考。

（References）：

[1] Boshuizen C R, Mason J, Klupar P, et al. Results from the planet labs flock constellation[C]//The 28th AIAA/USU Conference on Small Satellites, Utah, USA, 2014: 1-8.

[2] 廖文和. 立方體衛(wèi)星技術發(fā)展及其應用[J]. 南京航空航天大學學報, 2015, 47(6): 792-797. Liao Wen-he. A survey of Cubesat technology development and applications[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2015, 47(6): 792-797.

[3] 郁豐, 華冰, 吳云華, 等. 基于自適應卡爾曼濾波的簡化地磁定軌[J]. 中國慣性技術學報, 2014, 22(4): 519-524. Yu Feng, Hua Bing, Wu Yun-hua, et al. Simplified geomagnetic orbit determination based on adaptive Kalman filter[J]. Journal of Chinese Inertial Technology, 2014, 22(4): 519-524.

[4] Wang Xin-long, Zhang Qing, Li Heng-nian. An autonomous navigation scheme based on starlight, geomagnetic and gyros with information fusion for small satellites[J]. Acta Astronautica, 2014, 94: 708-717.

[5] 楚中毅, 任善永. 帶柔性伸桿小衛(wèi)星振動控制的半物理仿真實驗[J]. 宇航學報, 2013, 34(6): 748-754. Chu Zhong-yi, Ren Shan-yong. Semi-physical simulation on vibration control of small satellite with flexible manipulators[J]. Journal of Astronautics, 2013, 24(6): 748-754.

[6] 樂韻, 房建成, 湯繼強, 等. 磁懸浮反作用飛輪剩磁矩分析與補償方法研究[J]. 航空學報, 2011, 32(5): 881-890. Le Yun, Fang Jian-cheng, Tang Ji-qiang, et al. Research on analysis and compensation method of remnant magnetic moment for magnetically suspended reaction flywheel[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(5): 881-890.

[7] Steyn W H, Hashida Y. In-orbit attitude performance of the 3-axis stabilized SNAP-1 nanosatellite[C]//The 15th AIAA/USU Conference on Small Satellites. Utah, USA, 2001: 1-10.

[8] Yan Liang, Liu Delong, Jiao Zongxia. Magnetic field modeling based on geometrical equivalence principle for spherical actuator with cylindrical shaped magnet poles [J]. Aerospace Science and Technology, 2016, 49: 17-25.

[9] Lee D, Springmann J C, Spangelo S C, et al. Satellite dynamics simulator development using lie group variational integrator[C]//AIAA Modeling and Simulation Technologies Conference. Portland, Oregon, 2011: 547-566.

[10] Zhang Zhen, Jin Jin, Xiong Jian-ping. Research on attitude determination of micro-satellite with impact of remanence [C]//2013 International Conference on Optical Instruments and Technology: Optical Sensors and Application. Beijing, China, 2013: 1-10.

On-orbit estimation and compensation for CubeSats remanence

LI Ye1, YU Feng2, WANG Qun-yang1
(1. Space Star Technology Co. Ltd., Beijing 100086, China; 2. Nanjing University of Aeronautics and Astronautics, College of Astronautics, Nanjing 210016, China)

A CubeSat usually utilizes magnetometer, magnetorquer and momentum wheel to determination and control its attitude, and the remanence disturbance torque is one of factors that influence the performance of the attitude determination and control system. A new method for realizing the remanence’s on-orbit estimation and compensation is proposed which only utilizes data from the magnetometer and the magnetorquer. First, a remanence on-orbit estimation model is derived by employing the magnetometer outputs and satellite attitude dynamics, and an unscented Kalman filter is utilized to improve the estimation with a single magnetometer. Second, after simplifying the system into a linear time-invariant system, the remanence disturbance which affects the satellite attitude is mathematically modeled according to the superposition principle, and then a remanence compensation method with the criterion of removing angular velocity disturbance is presented in consideration of time-sharing between magnetometer and magnetorquer. Finally, simulations are made which show that the remanence’s estimation accuracy is better than 0.001 A·m2in 1000 s, and the heading, roll and pitch angles are reduced from 4.3°, 4.6° and 2.1° respectively to within 0.4°. The proposed method is a new way to reduce the remanence disturbance for satellites.

remanence; attitude determination and control system; on-orbit estimation; remanence compensation

V448.22

：A

2016-03-23；

：2016-05-28

國家自然科學基金資助項目（61203197）

李燁（1973—），男，高級工程師，從事導航、通信技術研究。E-mail: liyeal@sina.com

聯(lián) 系 人：郁豐（1980—），男，副研究員，從事微小衛(wèi)星技術、組合導航技術研究。E-mail: yufeng@nuaa.edu.cn

1005-6734(2016)03-0342-07

10.13695/j.cnki.12-1222/o3.2016.03.012