林佳偉，王 平

（1.北京控制工程研究所，北京100190；2.空間智能控制技術(shù)國家級重點實驗室，北京100190）

基于GTLS的零動量衛(wèi)星慣量矩陣在軌辨識

林佳偉1，2，王 平1

（1.北京控制工程研究所，北京100190；2.空間智能控制技術(shù)國家級重點實驗室，北京100190）

使用廣義總體最小二乘（GTLS，generalized total least squares）方法對零動量衛(wèi)星進行慣量矩陣在軌辨識.提出了GTLS算法的先驗最小距離解的定義：當測量信息不足以確定唯一解時，解空間中最接近先驗估計的解.給出了先驗最小距離解的算法，并應用于慣量矩陣在軌辨識.仿真結(jié)果表明了該辨識方法的有效性及先驗最小距離解相對于最小范數(shù)解的優(yōu)越性.

零動量衛(wèi)星；慣量矩陣；在軌辨識；廣義總體最小二乘；先驗估計

慣量矩陣參數(shù)是衛(wèi)星力學模型的重要參數(shù)，當前的工程實踐中衛(wèi)星的慣量矩陣參數(shù)是由衛(wèi)星結(jié)構(gòu)分析得到的.通過對衛(wèi)星結(jié)構(gòu)進行建模，預先估計帆板展開、推進劑消耗等引起的衛(wèi)星的結(jié)構(gòu)變化，從而獲得慣量矩陣參數(shù)的先驗估計.但是在衛(wèi)星的實際運行中預先估計可能發(fā)生偏差，例如：推進劑余量難以準確估計，液體在失重狀態(tài)下漂浮不定，太陽翼和天線等大型附件的展開不能為結(jié)構(gòu)模型所完全描述，從而使先驗估計偏離真值.以卡西尼衛(wèi)星為例，“干”星的慣量矩陣的不確定性為±10%［1-2］.通過在軌辨識，可以獲得比先驗估計更準確的慣量矩陣信息，建立準確的力學模型，并為控制器提供準確的參數(shù)，從而提高控制系統(tǒng)的性能.

當前國際上對慣量矩陣的辨識方法存在如下缺點：沒有采用完整的動力學方程；未考慮測量噪聲；未能保證估計的一致性［3］.文獻［1］提出了基于總體最小二乘（TLS，total least squares）的慣量矩陣在軌辨識方法，可以克服這些缺點.本文使用GTLS方法對該方法進行擴展，適用于更一般的模型，并具有更好的數(shù)值穩(wěn)定性.考慮到工程實際中對機動次數(shù)的限制可能導致測量信息不足，因此本文對GTLS算法進行了改進，以便于利用先驗信息.

1 估計方程

假設整個星體相對于慣性空間的總角動量（包括衛(wèi)星主體角動量和存儲在星上動量裝置的角動量）在初始時刻為0.沒有外力矩的作用時，根據(jù)動量守恒原理，在本體坐標系中有

式中，h為存儲在星上動量裝置的角動量，J為衛(wèi)星慣量矩陣，ω為衛(wèi)星主體角速度.

使用星上動量輪驅(qū)動機構(gòu)將h控制到某一個值，然后用動量輪轉(zhuǎn)速測量機構(gòu)測量h，同時用陀螺測量ω.由于測量噪聲的存在，方程組等號不成立，即h＋Jω≈0.為了抑制噪聲的影響，得到更好的估計，進行多次機動以獲得多組方程（一般而言為超定的非一致方程組）.

將方程轉(zhuǎn)化為ωTJ＝－h(huán)T.以慣量矩陣J作為待估計參數(shù)X，進行m次機動，相應的方程組為

從方程可以看出，要確定唯一的解，要求衛(wèi)星3軸都安裝動量裝置.另外還應預先估計整星相對于慣性空間的總角動量，推力器噴氣將整星角動量調(diào)為0.

當前工程中的動量輪轉(zhuǎn)速測量常使用測頻率的方法，即通過測量一定時間內(nèi)的轉(zhuǎn)角脈沖數(shù)來測量輪子的平均轉(zhuǎn)速，由于每一組機動結(jié)束后h的值不變，所以延長測量時間可以提高測量精度.同理也可以得到較高精度的角速度測量值.

2 求解算法

2.1 用GTLS算法進行慣量矩陣參數(shù)估計

考慮矩陣A∈Rm×n中有一些列是精確已知的情

形，即A＝[A1A2]，A1∈Rm×n1是不受噪聲污染的已知量，A2∈Rm×n2受到噪聲污染.未受擾方程為

而A20和B0的量測值為

式中：下標0表示未受噪聲污染的真值；待估計參數(shù)的真值X0∈Rn×d，X10∈Rn1×d，X20∈Rn2×d；量測值A2和B的真值分別為A20和B0，受到加性噪聲ΔA2和ΔB的污染；矩陣[ΔA2ΔB]均值為0，其各行向量之間獨立同分布，記其協(xié)方差陣C進行Cholesky變換得C＝GTG；假定A1和A20列滿秩.

可見此時方程是廣義誤差變量回歸（GEIV，generalized errors-in-variable）模型［4］.GTLS方法是針對GEIV模型的求解方法，要求滿足如下準則［4］：

準則1.尋找和使得range?range，并且最小.

其中，‖‖F(xiàn)指Frobenius范數(shù)（簡稱為F范數(shù)），range（·）指的是矩陣的值域，range?range等價于有解.

文獻［4］中給出了GTLS算法，可用來求解X.具體步驟如下.

1）若n1＝n，則退化為最小二乘方法；否則對［A1A2B］進行QR分解得

其中，Q為正交陣，R11∈Rn1×n1為上三角陣，R22∈R（m－n1）×（n2＋d）.若n1＝0，則把[A]B賦給R22.

2）若協(xié)方差陣C給定且與單位陣不成比例（即不等于單位陣乘上某個標量），則進行Cholesky分解，得C＝GTG；若C與單位陣成比例，則退化為TLS算法.

3）對矩陣對R22和G進行廣義奇異值分解，得

4）確定矩陣對R22和G的秩p.它可以由未受擾方程的性質(zhì)確定，也可以選取一個很小的數(shù)R0，若σ1≥…≥σp＞R0≥σp＋1≥…≥σn2＋d，則認為σp＋1≈…≈σn2＋d≈0.

5）將［zp＋1…zn2＋d］賦給Z2，其中zi為Z的列向量，下標表示列數(shù).再由R11Z1＝－R12Z2計算Z1.

6）若d＞1且p＜n2，通過QR分解，使交化，得＝QzRz，并將Qz賦給；否則進入下一步.

8）若Γ非奇異，則GTLS有解，對方程＾XΓ＝－Y進行求解可得＾X；否則應減小p的值再重新計算，將p－ρ賦給p，其中ρ是σp的重數(shù).

2.2 對算法的討論

GTLS方法在兩方面擴展了TLS方法：

1）經(jīng)典的TLS方法只能應用于行誤差向量的協(xié)方差陣為單位陣的情形，而GTLS方法擴展到了一般矩陣.

2）GTLS方法可以處理測量系數(shù)矩陣A中的一些列不含誤差的情形.

當n1＝0且G與單位陣成比例時，GTLS算法退化為經(jīng)典TLS算法，可見文獻［1］所提方法是本方法的特例.當n1＝0但G與單位陣不成比例時，GTLS算法理論上等價于對方程進行縮放變換后的TLS算法，此時估計具有一致性［1］.由于GTLS算法采用了廣義奇異值分解，避開了對G直接求逆，當協(xié)方差陣為病態(tài)矩陣時可以提高數(shù)值穩(wěn)定性.而當n1≠0，即A中有一些列精確已知時，無法使用TLS方法求解.

3 先驗最小距離解

有時候測量信息不足以給出唯一解，此時GTLS算法的步驟4）中接近0的奇異值的個數(shù)大于3.當機動次數(shù)不夠、機動設計不好（例如幾次機動完全相同）、動量裝置自由度小于3時，都可能出現(xiàn)這種情況.

2.1節(jié)中的原始GTLS算法在解空間中取出最小范數(shù)解，其物理意義不清楚，可能與真值區(qū)別很大.為了提高估計精度，考慮利用先驗信息得到更接近真值的估計.在此定義先驗最小距離解：當GTLS方法存在多解時，在其解空間中尋找一個與先驗估計距離最小的解.在這里用兩個解之差的F范數(shù)描述它們的距離.與最小范數(shù)解相比，先驗最小距離解具有明確的物理意義，一般而言更接近真值.本節(jié)中先給出求解算法，然后證明算法求出的解與先驗估計的距離最小.

3.1 討論前提與符號表示

對［A B］G－1進行奇異值分解，再對矩陣對［A B］和G進行廣義奇異值分解，兩組奇異值相同［6］.

在本節(jié)中，僅考慮n1＝0的情形，即所有列都包含誤差.2.1節(jié)步驟3）中，當n1＝0時，R22＝［A B］，對矩陣對［A B］和G進行廣義奇異值分解的符號與對矩陣對R22和G進行廣義奇異值分解的符號相同.

對［A B］G－1進行奇異值分解，得［A B］G－1＝UΣVT.其中U和V為［A B］G－1的左右奇異矩陣，分塊成V＝［V1V2］＝，U＝［U1U2］.V1和U1的列數(shù)為p，V11∈Rn×p，V22∈Rd×（n－p＋d）.考慮具有多重廣義奇異值σp＋1＝…＝σn＋d的情形，并且p＜n.當m≥n＋d時，Σ＝其中Σ1＝diag{σ1，σ2，…，σp}∈Rp×p，Σ2＝diag{σp＋1，…，σn＋d}＝σp＋1In＋d－p，In＋d－p表示單位陣，下標為維數(shù).記與σp＋1，…，σn＋d相對應的右奇異向量為νp＋1，…，νn＋d.

3.2 先驗最小距離解的求解算法

記X的先驗估計為Xp，當矩陣對［A B］和G有多重廣義奇異值σp＋1＝…＝σn＋d時，對[zp＋1…zn＋d]進行列變換得

式中，M是列初等變換陣，Y∈Rn×（n－p），Id為d維單位陣.

選取N1∈R（n－p）×d，使得‖Xp－S－YN1‖F(xiàn)取最小，這是一個最小二乘問題.

本小節(jié)只列出算法，下面證明所得的S＋YN1是解空間中與先驗估計的距離最小的解.

3.3 先驗最小解的證明

定理1.滿足GTLS準則1的‖［AG－1‖F(xiàn)的最小值為σn＋d.

證明.根據(jù)矩陣低秩逼近的Eckart-Young-Mirsky定理［7］：其中σ為E的奇異值，其下標按大小排列.若方程有解，則rankG－1｝≤n.由［A B］－當σn＋1＝…＝σn＋d時，

定理2.存在一個滿秩矩陣K，使得G［zp＋1…zn＋d］K＝［νp＋1…νn＋d］.

證明.根據(jù)奇異值分解的性質(zhì)［6-7］，［νp＋1…νn＋d］和G［zp＋1…zn＋d］的列向量都是相對于［A B］G－1的多重特征值的特征向量.νp＋1，…，νn＋d互相正交，是特征向量空間的正交基，而G［zp＋1…zn＋d］列滿秩，也是特征向量空間的一組基.所以K是特征向量空間中兩組基的轉(zhuǎn)換陣，它是滿秩的.

定理3.｛S＋YN｜?N∈R（n－p）×d｝是GTLS問題的解集，對于某一個給定的N，相應的

其中d維方陣R由K－1M進行QR分解而得，即K－1M＝QR，Q的各列向量為單位正交向量組.

證明.首先證明G－1的F范數(shù)為最小.

由K－1M＝QR得G－1［νp＋1…νn＋d］Q，將其代入的表達式得＝σp＋1U2QQT［νp＋1…νn＋d］T

第3個等號是由于V的各列向量為單位正交向量組而得的.

［νp＋1…νn＋d］T［νp＋1…νn＋d］、QTQ、都為單位陣，得

所以，YN＋S是＾AX＝＾B的解，即?N∈R（n－p）×d，S＋YN都存在一個相對應的［Δ＾AΔ＾B］，滿足2.1節(jié)的準則1.

定理4.S＋YN1是GTLS解集中的先驗最小距離解.

證明.定理3中指出｛S＋YN｜?N∈R（n－p）×d｝是GTLS問題的解集，而在3.2節(jié)的算法中，N1使‖Xp－S－YN1‖F(xiàn)取最小，即S＋YN1是解空間中與先驗估計Xp距離最小的解.

3.4 討 論

若要求解先驗最小距離解，并不需要對Z進行歸一化（2.1節(jié)的步驟3）.這是因為歸一化相當于對其右乘一個滿秩矩陣，而從證明過程中可以看出，這么做并不影響解集的表達式.

一般認為奇異值的個數(shù)為min｛m，n＋d｝，當m＜n＋d時，對［A B］G－1的奇異值分解UΣVT進行擴維得［U 0］VT，其中Σ2＝diag｛σp＋1，…，σn＋d｝，新增加的σm＋1，…，σn＋d為0，它們對應的左奇異向量為0.由于本文中只討論σp＋1＝…＝σn＋d的情形，此時Σ2＝0，對應測量信息不足的情形.代入證明過程中，結(jié)論依然成立.

實際上，無論是TLS算法還是GTLS算法，都可以看成是結(jié)構(gòu)總體最小二乘（STLS，structured total least squares）算法的特殊情形［8-10］.但是STLS方法一般情況下只有數(shù)值近似解法，無法保證求出指標函數(shù)的全局最小點，而TLS和GTLS算法由于模型的特殊性，可以求出使指標函數(shù)全局最小的解，從而保證了估計的一致性.由結(jié)構(gòu)總體最小二乘算法的性質(zhì)可知［9］，在高斯噪聲的條件下，GTLS是GEIV模型的極大似然估計.極大似然估計在漸近意義上能接近方差下界，有時被認為是最優(yōu)的估計，這就是選擇GTLS方法的原因.

4 仿真算例

首先驗證基于GTLS的慣量矩陣參數(shù)估計方法的有效性.慣量矩陣參數(shù)的真值為kg·m2.仿真條件為：角速度測量誤差的均值為0，標準差為0.003（°）／s；衛(wèi)星本體坐標系3軸的3個動量輪中，有兩個可以獲得精確的測量值，第3個動量輪的轉(zhuǎn)速測量包含高斯噪聲，噪聲均值為0，標準差為0.3 rad／min.3軸上的動量輪裝置的容量的上限為10 kg·m2／s.仿真1×103次，根據(jù)仿真結(jié)果計算估計誤差陣的F范數(shù)的均值.仿真結(jié)果如表1所示.

表1 估計誤差陣的F范數(shù)的均值

從表1中可以看出，隨著測量次數(shù)不斷增大，估計誤差越來越小.

當機動次數(shù)m為64時，進行1×104次仿真，得到估計結(jié)果的均值，它與真值的差為可以看出估計的均值與真值十分接近.

可見，仿真結(jié)果表明了基于GTLS的慣量矩陣參數(shù)估計方法的有效性.

下面比較最小范數(shù)解與先驗最小距離解.將先驗估計看成是真值加上一個加性噪聲矩陣，所有矩陣元素都是互不相關(guān)的零均值高斯噪聲，均方誤差為200 kg·m2.假設只進行兩次機動，所有測量都包含誤差，噪聲條件不變.

仿真1×104次，先驗估計、最小范數(shù)解、先驗最小距離解的誤差陣的F范數(shù)的均值分別為584.6、2 952.3和341.2這3種估計的均值與真值的差分別為

從仿真所得的估計誤差陣的F范數(shù)和均值可以看出，先驗最小距離解要優(yōu)于先驗估計和最小范數(shù)解.

5 結(jié) 論

本文使用GTLS方法進行零動量衛(wèi)星的慣量矩陣在軌辨識，提出了GTLS算法的先驗最小距離解的定義和算法，并將其應用在慣量矩陣辨識中.與文獻［1］提出的基于TLS的慣量矩陣在軌辨識方法相比，本文提出的方法適用于更一般的模型，具有更強的數(shù)值穩(wěn)定性.文中對GTLS算法進行了改進，當測量信息不足時在解空間中取出與先驗估計距離最小的解，相對于文獻［4］所提的最小范數(shù)解，先驗最小距離解充分利用了先驗估計信息，有利于提高估計精度.在工程實際中經(jīng)常使用動量輪對衛(wèi)星進行單次的姿態(tài)機動，其測量信息不足以確定唯一的慣量矩陣參數(shù)，此時可以應用先驗最小距離解.

本文的研究對象是零動量衛(wèi)星，采用了簡單的剛體模型，并且沒有考慮外界干擾.在后續(xù)的研究工作中應考慮放松這些限制，擴大本方法的應用范圍.

［1］ Lee A Y，Wertz JA.In-flight estimation of the Cassini spacecraft’s inertia tensor［J］.Journal of Spacecraft and Rockets，2002，39（1）：153-155

［2］ Wertz J，Lee A Y.In-flight estimation of the Cassini spacecraft’s inertia tensor［C］.The 11thAAS／AIAA Space Flight Mechanics Meeting，Santa Barbara，California，F(xiàn)eb 11-15，2001

［3］ 林佳偉，王平.零動量衛(wèi)星慣量矩陣的在軌辨識［C］.中國航天可持續(xù)發(fā)展高峰論壇暨中國宇航學會第三屆學術(shù)年會，北京，2008

［4］ Van Huffel S，Vandewalle J.Analysis and properties of the generalized total least squares problem AX≈B when some or all columns in A are subject to error［J］.SIAM J.Matrix Anal，1989，10（3）：294-315

［5］ Van Huffel S，Vandewalle J.The total least squares problem：computational aspects and analysis［M］.Philadelphia：SIAM，1991

［6］ Golub G H，Van Loan C F.Matrix computations［M］.3 rd ed.Baltimore，Maryland：The John Hopkins University Press，1996

［7］ 張賢達.矩陣分析與應用［M］.北京：清華大學出版社，2004

［8］ Lemmerling P.Structured total least squares：analysis，algorithms and applications［D］.Department of Electrical Engineering，Katholieke Universiteit，Leuven，Begium，1999

［9］ Abatzoglou T，Mendel J，Harada G.The constrained total least squares technique and its application to harmonic superresolution［J］.IEEE Trans.Signal Process，1991，39（5）：1070-1087

［10］ De Moor B.Total least squares for affinely structured matrices and the noisy realization problem［J］.IEEE Trans.Signal Process，1994，42（11）：3104-3113

In-Orbit Identification of the Inertial Matrix of Zero Momentum Satellite

LIN Jiawei1，2，WANG Ping1
（1.Beijing Institute of Control Engineering，Beijing 100190，China；2.National Laboratory of Space Intelligent Control，Beijing 100190，China）

A generalized total least squaresmethod is adopted to in-orbit identification of the inertialmatrix of a zero momentum satellite.A prior minimum distance solution is defined：this solution closest to the prior estimate in solution space as ifmeasurement information is not enough to determine a unique solution.An algorithm of the priorminimum distance solution is proposed，and applied to in-orbit identification of an inertialmatrix.Simulation results validate feasibility of the identification method and advantage of the priorminimum distance solution over the minimum norm solution.

zero momentum satellite；inertial matrix；in-orbit identification；generalized total least squares；prior estimation

2009-05-28

林佳偉（1982—），男，福建人，博士研究生，研究方向為航天器在軌辨識（e-mail：linjw0207＠gmail.com）.

V448.22

A

1674-1579（2009）05-0026-05