基于可观测性约束的行星着陆轨迹随机优化方法

摘要

本发明涉及一种基于可观测性约束的行星着陆轨迹随机优化方法，属于深空探测器导航与制导技术领域。本方法以基于单目视觉导航的深空着陆制导控制任务为背景，考虑有效控制与可靠估计间的对偶问题，通过扩展状态空间将系统不确定性作为部分代价引入二次型性能指标中，从而采用线性二次型控制技术给出随机优化反馈控制律，代入建立的探测器系统的扩展状态空间描述模型，实现对行星着陆轨迹的实时优化；避免了动态规划与基于搜索方法的运算负担，有效地克服了着陆过程中可观测性缺乏问题，提高系统的导航估计性能，使行星探测器导航制导控制整体性能得到保障，达到可靠着陆行星的最终目标。

说明书

技术领域

本发明涉及一种基于可观测性约束的行星着陆轨迹随机优化方法，属 于深空探测器导航与制导技术领域。

背景技术

行星着陆过程中，由于着陆器受到各种不确定性因素影响，其制导控 制律的设计必须建立在对系统状态最优估计的基础上，这使得精确的导航 定位显得尤为关键；加之着陆时间较短，深空环境中通信延迟时间较长， 基于地面站的导航制导控制模式不再适用，因此发展深空着陆自主实时精 确的导航制导控制策略成为目前国内外学者研究的热点。

在先技术[1]中（参见R.R.Sostaric,J.R.Rea.Powered descent guidance methods for the moon and mars.San Francisco,USA:American Institute of Aeronautics and Astronautics Inc,2005.），传统的制导 控制设计过程是基于确定性等价假设，其制导控制策略的选取与导航估计 不确定性互不影响，从而可将系统的导航制导控制设计分为两个独立的过 程，即利用导航滤波器从噪声污染测量中估计出系统状态；再假设这样的 估计状态即为系统的真实状态，利用传统的制导控制技术给出到达目标着 陆区域所需的控制量。由于这种方法简便易行且运算量低，因此成为以往 着陆段任务制导控制策略的选择。

然而，对于非线性的着陆器系统而言，其可观测性与状态轨迹间存在 非线性耦合，因此其可观测性会受到系统状态轨迹的影响。当观测能力受 到限制时，有限的观测能力会导致缺乏状态观测的估计性能不佳，进而弱 化制导控制性能，尤其对于具有不规则引力场特性的小行星着陆制导控制 而言，一个很小的控制偏差会导致最终着陆器与目标着陆点的大范围偏离， 达到可靠的着陆性能就变得极具挑战性。因此，如何利用尽可能少的导航 敏感器，通过优化观测轨迹，在满足着陆任务要求的同时，尽可能地减小 导航定位误差，从而提高导航制导控制系统的整体性能是行星着陆过程中 亟待解决的问题。

发明内容

本发明的目的在于针对行星着陆视觉导航过程中由于有限观测能力引 起的估计性能不佳问题，提出一种基于可观测性约束的行星着陆轨迹随机 优化方法。

本方法以基于单目视觉导航的深空着陆制导控制任务为背景，考虑有 效控制与可靠估计间的对偶问题，通过扩展状态空间将系统不确定性作为 部分代价引入二次型性能指标中，从而采用线性二次型控制技术给出随机 优化反馈控制律。具体包括如下步骤：

步骤1，建立探测器系统的扩展状态空间描述模型

$\stackrel{·}{\chi }=\left(\begin{array}{c}\stackrel{·}{x}\\ \stackrel{·}{s}\end{array}\right)=\Psi (\chi ,u)---\left(1\right)$

其中，探测器状态x=[r v]^T，r，v分别为探测器相对于预定着陆点的位 置矢量与速度矢量，u为探测器的控制输入；扩展状态变量$\chi ={\left(\begin{array}{cc}x& s\end{array}\right)}_{x=\hat{x}}^T,$其中变量s=[s₁₁ s₂₁ s₂₂…s₆₆]^T包含误差方差平方根阵S的所有非零元 素，ψ为探测器状态变量χ及控制输入u的函数向量。

步骤2，设计反馈控制输入u：

$u=-K(\chi -{\chi }_{\mathrm{ss}}){|}_{x=\hat{x}}=-K_x(\hat{x}-x_d)-K_s(s-s_{\mathrm{ss}})---\left(2\right)$

$K=W_u^{-1}{B}_e^T\Gamma =\left(\begin{array}{cc}{K}_x& K_s\end{array}\right)---\left(3\right)$

使其满足

$u^{*}=\arg \underset{u}{\min }{\int }_0^{t_f}\{{\chi }^T{W}_e\chi +u^T{W}_{u}u\}\mathrm{dt}---\left(4\right)$

其中，u^*为采用无限时长稳态增益状态反馈控制方法，在时间区间 [0,t_f)上寻找到的最优控制输入；χ_ss=[x_d s_ss]^T为扩展系统稳态解，x_d为期 望的探测器状态，s_ss为系统不确定状态元素的稳态解，K_x为探测器状态偏 差控制系数，侧重于稳定控制目标，K_s提供附加反馈控制，以削弱系统不 确定性，提高控制性能。

矩阵Γ通过求解代数Riccati方程得到：

$0=A_e^T\Gamma +\Gamma A_e+W_e-\Gamma B_e{W}_u^{-1}{B}_e^T\Gamma ---\left(5\right)$

其中，为状态估计；t_f为探测器着陆时 间；W_u为控制燃耗权重矩阵，$W_e=\left(\begin{array}{cc}{W}_x& 0\\ 0& W_s\end{array}\right),$W_x为状态偏差权重矩阵， 为扩展状态偏差权重矩阵。

步骤3，将步骤2得到的反馈控制输入u代入步骤1建立的探测器系统 的扩展状态空间描述模型，从而实现对行星着陆轨迹的实时优化。

有益效果

本发明方法作为一种行星探测器随机优化制导控制方法，将导航估计 不确定性引入并扩展系统状态，给出包含与估计方差相关的二次型性能指 标评价方法，避免了动态规划与基于搜索方法的运算负担；同时，利用着 陆状态可观测性与着陆轨迹之间的非线性耦合特性，合理优化观测轨迹， 有效地克服了着陆过程中可观测性缺乏问题；通过合理地优化着陆器下降 轨迹，大幅度提高系统的导航估计性能，使行星探测器导航制导控制整体 性能得到保障，达到可靠着陆行星的最终目标。

附图说明

图1为本发明方法的流程图；

图2为具体实施方式中的着陆点固连坐标系示意图。

具体实施方式

为了进一步说明本发明的目的和优点，下面结合附图和实施例对本发明 内容作进一步说明。

①在着陆点固连坐标系下建立探测器动力学方程

$\stackrel{·}{r}=\stackrel{·}{v}$

$\stackrel{·}{v}=u+g-a_e-a_k+a_{\Delta }$

其中，r，v分别为探测器相对于预定着陆点的位置矢量与速度矢量，u为探 测器的控制输入，g为行星引力加速度，a_e,a_k分别为由行星自旋所引起的 离心惯性加速度与科氏加速度，a_△为未建模加速度。

所述着陆点固连坐标系定义为∑^l:o_l-x_ly_lz_l，该坐标系原点o_l位于预定的 着陆点，o_lz_l轴与行星质心指向着陆点矢量o_io_l方向一致。o_lx_l沿经线的切线 方向指向南极方向，o_ly_l与o_lx_l、o_lz_l之间满足右手法则，如图2所示。其中， 坐标系∑ⁱ:o_i-x_iy_iz_i为行星中心惯性坐标系，其原点在行星的质量中心，z_i轴 沿行星最大惯量轴方向，x_i轴沿历元时刻行星最小惯量轴所指方向，y_i轴 与x_i轴、z_i轴之间满足右手法则。

②将式(6)中的加速度项-a_e-a_k+a_△作为系统过程噪声，得到系统动力 学模型：

$\stackrel{·}{x}=\left(\begin{array}{cc}0& I\\ 0& 0\end{array}\right)x+\left(\begin{array}{c}0\\ I\end{array}\right)\{u+g\left(E_{r}x\right)+\omega \}$

$=\mathrm{Ax}+B\{u+g\left(E_{r}x\right)+\omega \}$

其中，x=[r v]^T为探测器状态，E_r＝[I 0]为系数矩阵，A=[[0 0]^T[I 0]^T] 为系统矩阵，B=[0 I]^T为输入矩阵，ω=-a_e-a_k+a_△为系统过程噪声。

③采用基于单目视觉矢量测量的光学相对导航方案，选择着陆点固连 坐标系下探测器的位置矢量作为观测量，得到观测模型为

$z_k=h\left(x_c\left(t_k\right)\right)+{\upsilon }_k$

$={\left(\begin{array}{cc}f\frac{r_{\mathrm{cy}}}{r_{\mathrm{cx}}}& f\frac{r_{\mathrm{cz}}}{r_{\mathrm{cx}}}\end{array}\right)}^T+{\upsilon }_k$

其中，r_c＝C_clr＝[r_cx r_cy r_cz]^T为探测器相对于目标着陆点的位置矢量在相机 坐标系下的投影，C_cl为着陆点固连坐标系与相机坐标系之间的转换矩阵， υ_k为光学导航观测噪声，f为相机焦距。

所述相机坐标系为探测器自身加载的有效载荷。

④将探测器状态x描述为分布在状态空间中的粒子云x_n(n＝1,2,…,N， N为粒子数)。粒子云为围绕探测器飞行过程中的期望状态x_d，随机估计得 到；一个飞行期望状态对应一个粒子云。

期望状态x_d与其对应的粒子云中每个粒子x_n之间距离的平方和作为优 化的代价函数。此代价函数反映了粒子与期望状态间的距离以及其分散程 度。用二次型表达式描述探测器着陆过程中的某一期望状态x_d所对应时刻 的控制性能为：

$\mathrm{\delta J}=\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\{{(x_n-x_d)}^T{W}_x(x_n-x_d)+u^T{W}_{u}u\}---\left(9\right)$

$E[{(x-x_d)}^T{W}_x(x-\left(x_d\right)+u^T{W}_{u}u)$

其中，W_x为状态偏差权重矩阵，W_u为控制燃耗权重矩阵。

定义为x的状态估计，估计误差且满足

$\left(\begin{array}{c}E\left[\tilde{x}\right]=0\\ E\left[{\tilde{x}}^T\tilde{x}\right]=P\end{array}\right)---\left(10\right)$

则式(9)整理为

$\mathrm{\delta J}=E[{(\hat{x}+\tilde{x}-x_d)}^T{W}_x(\hat{x}+\tilde{x}-x_d)+u^T{W}_{u}u]$

$={(\hat{x}-x_d)}^T{W}_x(\hat{x}-x_d)+u^T{W}_{u}u+E\left[{\tilde{x}}^T{W}_x\tilde{x}\right]---\left(11\right)$

$={(\hat{x}-x_d)}^T{W}_x(\hat{x}-x_d)+u^T{W}_{u}u+\mathrm{tr}\left[W_{x}P\right]$

式(11)中前两项为典型的二次型，如果能将第三项描述为二次型形式， 则可以利用现有的数学工具对问题进行求解。注意到误差方差阵P为非负 定矩阵，由矩阵理论可知，对称的非负定矩阵P∈R^6×6能分解成P＝SS^T，其 中S∈R^6×6为下三角阵，由P唯一确定，即矩阵P的下三角分解平方根阵为

设W_x为一对角阵，且W_x(i，i)=w_ii(i=1,2,…,6)，则式(11)中最后一项整 理为

$\mathrm{tr}\left[W_{x}P\right]=\mathrm{tr}\left[W_x{S}^{T}S\right]=\underset{i=1}{\overset{6}{\Sigma }}\left\{w_{\mathrm{ii}}\underset{j=1}{\overset{i}{\Sigma }}{s}_{\mathrm{ij}}^2\right\}---\left(13\right)$

式(11)改写为状态估计偏差控制输入u和误差方差平方根阵S非 零元素的二次型函数：

将式(14)对时间积分，得到随机优化的性能指标：

$J={\int }_0^{t_f}\{{(\hat{x}-x_d)}^T{W}_x(\hat{x}-x_d)+u^T{W}_{u}u+\underset{i=1}{\overset{6}{\Sigma }}\{w_{\mathrm{ii}}\underset{j=1}{\overset{i}{\Sigma }}{s}_{\mathrm{ij}}^2\left\}\right\}\mathrm{dt}---\left(15\right)$

式中前两项对应于确定性等价控制代价，第三项对应于系统不确定性代 价。

⑤定义扩展状态变量$\chi ={\left(\begin{array}{cc}x& s\end{array}\right)}_{x=\hat{x}}^T,$其中变量s=[s₁₁ s₂₁ s₂₂…s₆₆]^T包含 误差方差平方根阵S的所有非零元素，为便于描述，设期望状态x_d=0，则 式(15)所描述的二次型性能指标可简化表示为

$J={\int }_0^{t_f}\{{\chi }^T{W}_e\chi +u^T{W}_{u}u\}\mathrm{dt}---\left(16\right)$

其中，$W_e=\left(\begin{array}{cc}{W}_x& 0\\ 0& W_s\end{array}\right),$为扩展状态偏差权重矩 阵。

⑥深空探测器下降轨迹优化问题的扩展状态空间描述模型为：

$\stackrel{·}{\chi }=\left(\begin{array}{c}\stackrel{·}{x}\\ \stackrel{·}{s}\end{array}\right)=\Psi (\chi ,u)---\left(17\right)$

在时间区间[0,t_f)上寻找最优控制输入u^*，满足

$u^{*}=\arg \underset{u}{\min }{\int }_0^{t_f}\{{\chi }^T{W}_e\chi +u^T{W}_{u}u\}\mathrm{dt}---\left(18\right)$

由于涉及协方差传播和扩展空间描述的复杂性，寻求非线性控制系统的 最优解十分困难。这里采用无限时长稳态增益状态反馈控制技术获得优化 问题的次优解。

具体过程为：求解满足代数Riccati方程的矩阵Γ：

$0=A_e^T\Gamma +\Gamma A_e+W_e-\Gamma B_e{W}_u^{-1}{B}_e^T\Gamma ---\left(19\right)$

其中，$A_e=\frac{\partial \Psi (\chi ,u)}{{\partial }_{\chi }}{|}_{x=\hat{x}},B_e=\frac{\partial \Psi (\chi ,u)}{\partial u}{|}_{x=\hat{x}},$则反馈输入形式为：

$u=-K(\chi -{\chi }_{\mathrm{ss}}){|}_{x=\hat{x}}=-K_x(\hat{x}-x_d)-K_s(s-s_{\mathrm{ss}})---\left(20\right)$

$K=W_u^{-1}{B}_e^T\Gamma =\left(\begin{array}{cc}{K}_x& K_s\end{array}\right)---\left(21\right)$

其中，χ_ss=[x_d s_ss]^T为扩展系统稳态解，K_x主要侧重于稳定控制目标，K_s提 供附加反馈控制用于削弱系统的不确定性，从而提高控制性能。

以小行星Eros433作为目标天体，其自旋角速度为1639.4°/day，名义半 径为16km，引力常数为4.4621×10⁵m³/s²，在仿真中采用其四阶引力场模型。 在着陆点固连坐标系下，探测器初始位置为[500 20 -50]^Tm，相对小行星 表面速度为[1.5 20]^Tms，预期着陆时间300s。考虑着陆器初始状态估计 误差，其位置初始估计误差与速度估计误差分别服从方差为100m²和 0.1m²/s²的随机分布时，利用先技术[1]对着陆器进行前馈控制，在预期着 陆时间内着陆器最终着陆速度高达3ms，在期望着陆时间内存在撞入目标 天体的危险。这主要是由于探测器状态的可观测性受系统状态轨迹的影响， 当观测能力受到限制时，导航估计误差将持续增大，且其制导控制性能会 随着导航误差的增大被大幅度弱化。考虑导航估计误差的存在，在行星着 陆控制过程中加入本发明设计的反馈控制输入对着陆轨迹进行随机优化， 探测器的估计误差随着下降轨迹的变化在不断减小，末端着陆器真实状态 与前馈制导所产生的路径点基本相符，其最终着陆位置偏差在5m范围内， 速度偏差在0.05ms左右。通过反馈控制对下降轨迹的随机优化作用使导航 与制导控制性能得到了整体提高。