基于機器學習的空間翻滾目標實時運動預測

余敏，羅建軍，王明明

1. 西北工業大學 深圳研究院，深圳 518057 2. 西北工業大學 航天動力學國家重點實驗室，西安 710072

空間機器人被認為是執行空間在軌服務任務(例如維修失效衛星)最有效的手段之一[1]，為保證在軌服務任務的安全性和可靠性，對空間目標進行快速、準確的運動狀態預測成為亟需，而空間目標大都呈現翻滾運動模態[2]。

現有關于空間目標的運動預測主要是基于物理模型的方法，需要先基于視覺信息對目標進行重構與建模。借助視覺敏感器和激光雷達傳感信息，通過多臺同步相機成像[3]，結合圖像處理和點云數據處理技術對空間目標進行三維重構[4-6]，基于目標歷史觀測運動狀態的分幀序列進行特征提取，獲得目標的構型和運動狀態測量信息，進而利用基于物理模型(例如非線性濾波)對目標的未來運動狀態進行預測。由于空間自由漂浮的翻滾目標幾乎不受外力/力矩作用[7]，給定目標的動力學參數和初始運動狀態，通過動力學積分可以得到目標的未來運動狀態。基于物理模型的預測方法的成功依賴于準確估計目標動力學參數和理想的建模假設。Hirzinger等[8]利用擴展卡爾曼濾波方法預測一個慣性參數已知的目標的線性運動，但只能預測未來較短時間內的運動狀態。文獻[9] 利用無損卡爾曼濾波預測無人機的運動狀態。Greenspan等[10]采用霍夫變換預測目標位姿，但其更注重計算速度，沒有考慮預測效果。Hillenbrand和Lampariello[3]利用線性最小二乘方法對自由漂浮目標進行辨識和運動估計，提出一種有效的長周期目標運動預測方法，但由于旋轉運動的非線性特性，該方法對目標旋轉運動的預測效果不佳。借助視覺反饋系統，Aghili[11]利用卡爾曼濾波估計自由漂浮目標的動力學參數和狀態，考慮了目標的不確定特性和數據噪聲。上述方法雖然能夠對空間目標進行預測，但是預測條件和前提假設較為嚴苛，很少考慮運動預測的計算效率。

近年來，機器學習的興起為目標運動預測問題提供了一種新的解決思路。類似于人的感知行為，從過去的經驗中預測未來的事件，機器學習方法通過對大量數據的學習，尋找過去事件與未來事件的隱含關系。Amalia等[12]利用部分可觀測馬爾科夫過程預測動態環境的狀態遷移，在擁擠的環境中為機器人進行自主導航。Sung等[13]從目標的觀測數據中通過聚類分析提取其運動模式，利用隱馬爾科夫過程預測其未來運動狀態，但是隱馬爾科夫過程預測的準確率相對較低。Peng和BAI[14]利用支持向量機、神經網絡和高斯過程回歸[15](Gaussian Process Regression, GPR)3種機器學習方法進行了軌道預測的研究，并對比分析了3種方法的優劣，發現支持向量機的預測效果不如神經網絡和GPR的預測效果。事實上，當神經網絡存在無限個隱層節點時，該神經網絡等價于高斯過程回歸[16]。GPR是一種基于貝葉斯理論和統計學習理論發展起來的監督式機器學習方法。Heravi和Khanmoh ammadi[17]利用GPR對運動目標進行長周期運動預測。Kim等[18]提出一種魯棒自回歸GPR，預測擁擠環境中行人的運動。由于GPR是一種靈活的非參數推斷方法，非參數特性直接導致其計算量較大的缺陷。為提高計算效率，相關研究者們提出了多種近似GPR方法，其中Snelson[19]提出的稀疏偽輸入高斯過程(Sarse Peudo-input Gaussian Pess, SPGP)回歸方法，被認為是最有效的近似GPR方法之一。

本文的創新之處在于提出一種兼顧計算效率和預測精度的啟發式機器學習方法，在較少數據驅動的情形下，仍然能夠精確實時地預測未來有限時域內空間翻滾目標的運動狀態。基于SPGP回歸，用較少的偽數據集代替目標的真實觀測數據，結合先驗知識和數據知識，不斷更新學習結果，進而預測目標未來有限時域內的運動狀態(包括位置和姿態)。為保證良好的目標預測效果，利用馬爾科夫鏈蒙特卡洛法代替原SPGP中的共軛梯度優化方法，優化偽數據等信息，避免由于學習數據較少，造成優化過程陷入局部極小值，在提高計算效率的同時，保證空間翻滾目標長期運動預測的精確性。

1 空間翻滾目標運動建模

如圖1所示，將空間翻滾目標簡化為一個長方剛體，為便于分析，取其上點A為抓捕點。目標的本體坐標系otxtytzt固連于目標，其原點位于目標質心。

圖1 空間翻滾目標及其上抓捕點AFig.1 Space tumbling target with grapple fixture A

定義目標的轉動慣量為Ι=diag(Ιx,Ιy,Ιz)，其翻滾角速度為ω=[ωx，ωy，ωz]T，空間翻滾目標的歐拉動力學方程在otxtytzt中的分量可表示為

(1)

(2)

式中：A(q)為本體坐標系到慣性坐標系的姿態變換矩陣；inv(A)表示對矩陣A求逆。

圖2 抓捕點A的運動軌跡Fig.2 Motion trajectory of grapple fixture A

由圖2可知，空間翻滾目標的運動狀態呈現極其復雜的非線性運動特性，已有的運動預測算法，例如最小二乘估計方法并不能有效預測目標的非線性運動[3]。本文所提基于機器學習的方法由于其本身非線性概率建模原則，能夠有效預測非線性運動，且不需要目標的動力學模型，僅需要目標運動狀態的觀測數據。

本節建立空間翻滾目標的真實動力學模型目的如下。

1) 闡述空間翻滾目標的運動形式，說明所提算法對此類非線性翻滾運動預測的適用性。

2) 利用建立的動力學模型生成目標的歷史觀測數據，作為學習算法中的訓練數據。

3) 利用真實模型的仿真數據對比預測算法的運動預測數據，分析預測誤差。

2 標準高斯過程回歸

高斯過程[15]是任意有限個具有聯合高斯分布的隨機變量的集合，由其均值函數m(x)和協方差函數k(x,x′)表示。一個高斯過程表述為

f(x)～GP(m(x),k(x,x′))

(3)

假設潛在函數f(X)具有零均值的高斯先驗分布，即

(4)

θ=argmax(lnp(y|X,θ))

(5)

y=f(X)+ε

(6)

式中：ε為零均值高斯白噪聲。式(6)可等效為如下似然函數:

(7)

給定新的測試輸入數據x*，那么對應輸出數據y*的預測分布為

(8)

式中：

(9)

其中：k*=k(x*,xn)為測試輸入數據x*和訓練數據xn的協方差向量。類似地，K**=K(x*,x*)。

3 啟發式稀疏偽輸入高斯過程回歸

為提高目標運動預測的計算效率，本文采用稀疏偽輸入高斯過程回歸方法，通過分析目標真實觀測數據的特點，利用啟發式優化方法優化真實數據，得到與真實數據具有相似分布特性的稀疏的偽訓練數據，實現對目標實時、準確的運動狀態預測。

3.1 稀疏偽輸入高斯過程回歸

(10)

式中：kx=k(xm,xn)

推導得到真實數據集的似然函數為

(11)

式中：

(12)

(13)

根據貝葉斯準則，得到偽輸出數據的后驗分布如下：

(14)

給定新的測試輸入數據x*，聯立式(10)和式(14)，計算預測為

(15)

式中：

在稀疏偽輸入高斯過程訓練時，一般采用共軛梯度法優化超參數以及偽輸入。然而，梯度優化法對隨機初始值敏感，優化過程容易陷入局部極小解，即機器學習中常見的過擬合現象。那么，較差的訓練結果加上訓練數據較少，所建立的概率模型無法充分學習目標的潛在運動規律，進而無法構建準確的預測分布信息，導致預測精度降低。對于所提的目標運動預測問題，為保證較高的計算效率，又不影響稀疏偽輸入高斯過程回歸的預測效果，本文采用啟發式優化算法聯合優化訓練過程中的超參數和偽輸入，克服由于對隨機初始值敏感的問題造成的預測效果不佳。

3.2 啟發式優化：馬爾科夫鏈蒙特卡洛法

馬爾科夫鏈蒙特卡羅(Markov Chain Monte Carlo, MCMC)法[20]是統計模擬領域中著名的啟發式優化方法，通過對狀態空間進行大量采樣和游走，能夠有效解決由隨機初始化造成的局部極小問題，本文利用MCMC優化訓練過程中的偽輸入以及超參數。

首先給出馬爾科夫鏈及其平穩分布特性，作為MCMC優化算法的基礎。

3.2.1 馬爾科夫鏈及其平穩分布

具有馬爾科夫性質的隨機狀態變量b1,b2,…,bn可以構建一條馬爾科夫鏈，如下：

P(bt+1=b′|bt,bt-1,…)=P(bt+1=b′|bt)

(16)

式中：P為轉移概率矩陣。

給定初始狀態的概率分布π0和固定的轉移概率矩陣P，那么未來狀態的概率分布為

(17)

無論選取什么初始概率分布π0，經過足夠次數的狀態轉移，馬爾科夫鏈最終都會收斂于同一個穩定的狀態概率分布，此收斂現象是大多數馬爾科夫鏈(非周期馬爾科夫鏈)的共同行為。此收斂行為由轉移概率矩陣P決定，最終平穩的概率分布計算為

p(b)=Pn,n→∞

(18)

3.1.2 馬爾科夫鏈蒙特卡洛優化

由于馬爾科夫鏈能收斂到平穩分布，很自然的想法就如何構造一個轉移概率矩陣為P的馬爾科夫鏈，使得該馬爾科夫鏈的平穩分布恰好是p(b)。那么，無論從何種初始狀態出發，沿著該馬氏鏈轉移，最終都會收斂到期望的概率分布，同時得到期望概率分布的若干樣本。

馬爾科夫鏈的收斂性質主要由轉移概率矩陣P決定，因此以馬爾科夫鏈方式做大量采樣的關鍵在于如何構造P，以得到期望分布p(b)。為此，引入以下定理。

定理1(細致平穩條件)如果非周期馬爾科夫鏈的轉移概率矩陣P和概率分布π(b)滿足:

π(bi)pij=π(bj)pji，?i,j

(19)

則π(b)為該馬爾科夫鏈的平穩分布。

假設一個轉移概率矩陣為Q的馬爾科夫鏈，qij表示從狀態bi轉移到狀態bj的概率。通常情況下:

p(bi)qij≠p(bj)qji

(20)

即細致平穩條件不成立，此時狀態分布p(b)并非該馬爾科夫鏈的平穩分布。引入接受率α，并使得

p(bi)qijαij=p(bj)qjiαji

(21)

為使得式(21)成立，基于對稱性原則，取

αij=p(bj)qji,αji=p(bi)qij

此時，式(21)顯然成立。令

q′ij=qijαij

q′ji=qjiαji

式(21)可改寫為

p(bi)q′ij=p(bj)q′ji

(22)

原轉移概率矩陣為Q的馬爾科夫鏈就轉化成一個轉移概率矩陣為Q′的新馬爾科夫鏈，而新鏈恰好滿足細致平穩條件，其平穩分布為p(b)。

接受率αij的物理含義是在原馬氏鏈上，以狀態bi從qij的概率轉移到狀態bj時，以αij的概率接受這個轉移。注意，如果αij取值過小，那么在采樣過程中，馬爾科夫鏈將拒絕大量的跳轉，馬爾科夫鏈遍歷所有的狀態空間將耗費大量時間，收斂速度變慢。為提高采樣過程中的接受率，取

(23)|

即同比例放大接受率αij和αji。

對于稀疏偽輸入最優化求解問題，期望的目標函數是偽輸入關于真實數據的極大似然函數(式(13))，即MCMC算法中期望的平穩分布。不同于梯度優化算法在每次迭代優化中追求更優的解，MCMC算法在每次迭代過程中以一定的概率α接受次優解，從而具備從局部極小值中跳出的能力。因此，無論給定何種初始猜想值，經過足夠次數的隨機采樣后，MCMC算法將克服由于隨機初始猜想造成的敏感性問題最終找到全局最優解，保證稀疏偽數據的質量，進而構造更準確的預測分布信息。

3.3 空間翻滾目標運動預測

表1 MCMC-SPGP算法Table 1 MCMC-SPGP algorithm

采用多步直接預測策略，所提稀疏偽輸入高斯過程回歸的多步預測表達式為

MCMC-SPGP(yt,yt-1,…)

(24)

對于空間翻滾目標上任一抓捕點，其運動狀態可以表示為[ri,qj]T,i=x,y,z,j=1,2,3,4.其中，ri為位置信息，qj為姿態信息。根據式(24)，目標的長期預測表達式為

(25)

4 仿真試驗

4.1 利用Snelson數據驗證算法

為驗證所提MCMC-SPGP算法的正確性，首先利用Snelson驗證SPGP算法的原始數據測試所提算法，同時與標準高斯過程(Gaussian Process, GP)和SPGP的回歸效果進行對比，仿真結果如圖3所示。

圖3 3種算法的Snelson數據預測分布Fig.3 Predictive distributions for Snelson’s data of three algorithms

Snelson數據[21]包含200個訓練數據和301個測試數據。利用基于高斯過程的回歸模型，學習訓練數據與測試輸入數據之間的隱含規律，求解301個測試輸入數據對應的輸出數據的預測分布。GP算法利用完整的200個訓練數據去學習隱含規律，而SPGP和MCMC-SPGP算法中僅利用20個偽數據代替真實的200個訓練數據。

圖3中藍線是預測分布的均值，兩條紅線是標準偏差，品紅色星號為真實的測試輸出值。注意，在圖3的圖3(b)和圖3(c)中，上方的紅色加號表示隨機初始偽輸入x軸位置信息，其y軸信息無物理意義。由圖3可知，3種算法都能獲得良好的回歸效果，驗證了所提MCMC-SPGP算法的正確性。

為說明所提算法的優越性，從均方根誤差和計算效率2個方面，對圖3做進一步數據分析，分析結果見表2。

表2 3種算法的性能對比Table 2 Performance comparison of three algorithms

由表2可知，MCMC-SPGP算法的均方根誤差略高于GP和SPGP算法，但近似相等。在訓練階段，GP和SPGP算法均采用共軛梯度方法優化超參數等，GP算法的訓練時間是0.489 1 s，SPGP算法的訓練時間是0.879 7 s，因為SPGP算法中還需要優化偽輸入，因此其訓練時間高于GP算法。所提MCMC-SPGP算法的訓練時間為1.079 8 s，其優化超參數、偽輸入的時間與MCMC算法中設置的打靶次數有關，為保證良好的回歸效果，打靶次數設置為2 000次。在預測階段，GP算法的預測時間為0.008 4 s，SPGP和MCMC-SPGP算法的預測時間分別為0.000 8 s和0.000 7 s。稀疏偽輸入高斯過程的預測效率是標準高斯過程預測效率的10倍，進一步說明了所提MCMC-SPGP算法在保證良好的回歸效果的同時，也保證了較高的預測效率。

4.2 標準GP算法的目標預測結果

第1節給出了空間翻滾目標的運動模型，以此目標為運動預測的研究對象。目標的轉動慣量為I=[50.3,0，0;0，105.18，0;0，0，105.97] kg·m2，目標初始翻滾角速度為[-4，-2，-4](°)/s。所提基于機器學習的預測模型僅需要目標歷史運動的觀測數據，作為學習算法的輸入數據。給定目標的歷史運動狀態數據(包括位置和姿態，即[ri,qj]T,i=x,y,z,j=1,2,3,4)，利用標準GP算法進行目標的運動預測。

給定0～20 s內的目標運動狀態數據，采樣時間為0.1 s，則有201個訓練數據，利用標準GP算法預測后續20 s內目標的運動狀態，預測結果如圖4所示。

圖4的前2幅圖中，藍色曲線代表目標位置和姿態的真實信息，綠色曲線代表GP算法預測分布的均值。由圖4可以看出，在20～40 s的測試階段，預測均值與真實數據幾乎吻合，位置預測誤差在毫米級別，姿態預測誤差不超過0.01。此外，由于高斯過程回歸模型滿足一致性條件[12]，預測誤差隨預測時間逐漸增加，因此越遠離訓練數據，其預測誤差越大。上述仿真結果說明GP算法對空間翻滾目標的長期運動預測效果良好。

圖4 利用GP算法的翻滾目標運動預測結果Fig.4 Simulation results by GP algorithm for motion prediction of tumbling target

4.3 SPGP算法的目標預測結果

與4.2節中所有仿真條件相同，本小節利用SPGP算法預測目標的運動狀態。在SPGP算法中，給定目標歷史運動觀測的201個數據，基于極大化邊際似然函數原則，通過共軛梯度法優化得到20個偽數據集，作為訓練數據。預測結果如圖5所示。

圖5 利用SPGP算法的翻滾目標運動預測結果Fig.5 Simulation results by SPGP algorithm for motion prediction of tumbling target

圖5中，紅色加號表示隨機初始偽輸入x軸位置信息，其y軸信息無物理意義。圖5的前七幅圖中，藍色曲線代表目標位置和姿態的真實信息，綠色曲線代表SPGP算法預測分布的均值。顯然，單純利用SPGP算法的預測效果較差，特別是對于目標的位置預測，位置預測的最大誤差為0.6 m，姿態預測誤差達到了0.1。這是由于在訓練過程中，利用共軛梯度法優化偽輸入和超參數時，對隨機初始猜想敏感，容易陷入局部極小解，從而導致SPGP算法在稀疏數據條件下，無法充分學習到訓練數據與測試數據之間的潛在聯系。

4.4 MCMC-SPGP算法的目標預測結果

與4.2節中所有仿真條件相同，本小節利用所提MCMC-SPGP算法預測目標的運動狀態。在MCMC-SPGP算法中，給定目標歷史運動觀測的201個數據，基于極大化邊際似然函數原則，通過MCMC算法優化得到20個偽數據集，作為訓練數據。在MCMC優化過程中，打靶次數設置為2 000次。預測結果如圖6所示。

圖6 利用MCMC-SPGP算法的翻滾目標運動預測結果Fig.6 Simulation results by MCMC-SPGP algorithm for motion prediction of tumbling target

由圖6可以看出，本文所提MCMC-SPGP算法對空間翻滾目標的長期運動預測效果良好。在20～40 s的測試階段，預測均值與真實數據基本吻合，位置預測誤差不超過0.02 m，姿態預測誤差不超過0.02。為說明所提算法的高效性，從訓練時間和預測時間2個方面，對4.2節、4.3節和4.4節仿真算例做進一步數據分析，分析結果見表3。

表3 3種算法的計算時間對比

由表3可知，在訓練階段，GP和SPGP算法均采用共軛梯度方法進行超參數等的優化，訓練時間分別為2.017 7 s、4.582 6 s。本文所提MCMC-SPGP算法的訓練時間為0.932 4 s，明顯優于GP和SPGP算法的訓練時間。在預測階段，GP算法的預測時間為0.044 0 s，SPGP和MCMC-SPGP算法的預測時間分布為0.004 2 s和0.003 6 s。雖然SPGP算法預測效率較高，但預測效果不好。而本文所提MCMC-SPGP算法不僅顯著提高了預測效率，同時保證了良好的預測效果。

圖7的前3幅圖中，品紅色實心點代表含噪聲的訓練數據。2條正紅色曲線代表標準偏差，灰色區域是置信度為95%的置信區間。藍色曲線代表目標位置的真實信息，綠色曲線代表MCMC-SPGP算法預測分布的均值。由圖7可以看出，當數據含有噪聲時，所提算法的預測效果依舊良好，位置預測誤差不超過0.04 m，姿態預測誤差不超過0.05，且真實數據都落在置信區間內，驗證了所提算法的魯棒性。顯然，隨著噪聲數量級的增大，預測效果會降低，很多文獻中也給出了更詳盡的應對噪聲的處理方法[22]，但這不是本文的重點研究內容。

圖7 考慮噪聲利用MCMC-SPGP算法的 翻滾目標運動預測結果Fig.7 Simulation results by MCMC-SPGP algorithm with noise for motion prediction of tumbling target

5 結 論

本文提出了一種純數據驅動方式的空間翻滾目標運動預測方法，總結如下：

1) 基于監督式機器學習領域中的稀疏偽輸入高斯過程回歸方法，給定目標運動的歷史觀測數據，實現了對空間翻滾目標運動狀態的實時長期預測，且對數據噪聲具有一定的魯棒性。

2) 利用啟發式馬爾科夫鏈蒙特卡洛優化方法，在訓練過程中聯合優化偽輸入和超參數，解決了隨機初始猜想造成的敏感性問題，避免了稀疏偽輸入高斯過程回歸方法對于目標長期運動預測問題的過擬合現象，提高了目標運動預測的預測精度。

[21] GAO Y S. Github[EB/OL].(2013-11-13)[2020-04-27].http:11gitub.com/tyanshuaicao/gp_cholqr.