基于机器学习的空间翻滚目标实时运动预测

2021-03-26余敏罗建军王明明

航空学报 2021年2期

余敏，罗建军，王明明

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]，但这不是本文的重点研究内容。