林佳偉，王 平

（1.北京控制工程研究所，北京100190；2.空間智能控制技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京100190）

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

林佳偉1，2，王 平1

（1.北京控制工程研究所，北京100190；2.空間智能控制技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京100190）

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

零動(dòng)量衛(wèi)星；慣量矩陣；在軌辨識(shí)；廣義總體最小二乘；先驗(yàn)估計(jì)

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

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

1 估計(jì)方程

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

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

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

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

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

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

2 求解算法

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

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

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

而A20和B0的量測(cè)值為

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

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

準(zhǔn)則1.尋找和使得range?range，并且最小.

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

文獻(xiàn)［4］中給出了GTLS算法，可用來(lái)求解X.具體步驟如下.

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

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

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

3）對(duì)矩陣對(duì)R22和G進(jìn)行廣義奇異值分解，得

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

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

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

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

2.2 對(duì)算法的討論

GTLS方法在兩方面擴(kuò)展了TLS方法：

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

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

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

3 先驗(yàn)最小距離解

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

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

3.1 討論前提與符號(hào)表示

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

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

對(duì)［A B］G－1進(jìn)行奇異值分解，得［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.當(dāng)m≥n＋d時(shí)，Σ＝其中Σ1＝diag{σ1，σ2，…，σp}∈Rp×p，Σ2＝diag{σp＋1，…，σn＋d}＝σp＋1In＋d－p，In＋d－p表示單位陣，下標(biāo)為維數(shù).記與σp＋1，…，σn＋d相對(duì)應(yīng)的右奇異向量為νp＋1，…，νn＋d.

3.2 先驗(yàn)最小距離解的求解算法

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

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

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

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

3.3 先驗(yàn)最小解的證明

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

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

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

證明.根據(jù)奇異值分解的性質(zhì)［6-7］，［νp＋1…νn＋d］和G［zp＋1…zn＋d］的列向量都是相對(duì)于［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問(wèn)題的解集，對(duì)于某一個(gè)給定的N，相應(yīng)的

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

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

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

第3個(gè)等號(hào)是由于V的各列向量為單位正交向量組而得的.

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

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

定理4.S＋YN1是GTLS解集中的先驗(yàn)最小距離解.

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

3.4 討 論

若要求解先驗(yàn)最小距離解，并不需要對(duì)Z進(jìn)行歸一化（2.1節(jié)的步驟3）.這是因?yàn)闅w一化相當(dāng)于對(duì)其右乘一個(gè)滿秩矩陣，而從證明過(guò)程中可以看出，這么做并不影響解集的表達(dá)式.

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

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

4 仿真算例

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

［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］ 林佳偉，王平.零動(dòng)量衛(wèi)星慣量矩陣的在軌辨識(shí)［C］.中國(guó)航天可持續(xù)發(fā)展高峰論壇暨中國(guó)宇航學(xué)會(huì)第三屆學(xué)術(shù)年會(huì)，北京，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］ 張賢達(dá).矩陣分析與應(yīng)用［M］.北京：清華大學(xué)出版社，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—），男，福建人，博士研究生，研究方向?yàn)楹教炱髟谲壉孀R(shí)（e-mail：linjw0207＠gmail.com）.

V448.22

A

1674-1579（2009）05-0026-05