天问一号火星探测任务是我国首次火星探测任务，2020年7月，用于执行任务的天问一号探测器在海南文昌发射场发射，这次火星探测任务完成了环绕探测和巡视探测两项任务。2021年5月15日，着陆探测器成功软着陆于火星乌托邦平原南部，预示着探测任务工程目标即将圆满完成，并将继续开展科学探测工作。火星探测器由环绕器和着陆巡视器组成，环绕器将环绕火星进行探测，获取火星全面、综合的信息，而着陆巡视器则着陆在火星表面，分为火星车和着陆平台2个部分，并由火星车来完成火星表面的探测工作。

“祝融号”火星车及着陆平台的结构如图 1所示，火星车车体采用箱板式构型。车体顶板上方安装了太阳翼和桅杆，桅杆顶部安装有一对导航相机。着陆巡视器到达火星表面之后，火星车依次完成桅杆展开、太阳翼展开、定向天线展开、车体抬升、导航相机环拍等动作，为自身补充能量、获取自身健康信息的同时完成对着陆环境的了解，也为驶离动作提供数据支持。

火星车驶离着陆平台是探测器系统着陆后进行的第一个关键动作，是后续火星表面区域巡视探测工作的基础，因此保证火星车驶离动作的安全有着十分重要的工程意义。驶离工作的风险主要来源于着陆平台驶离轨道在发射、飞行、着陆阶段产生的位移和形变，以及着陆器着陆点地形平整度不明，可能会产生着陆平台驶离轨道坡度角过大、2驶离轨道间异面角过大等情况，从而导致火星车在驶离过程中脱轨或倾翻。因此，在着陆巡视器到达火星表面并展开动作后，立即对包括驶离轨道在内的着陆点周围环境进行拍摄，进而计算驶离轨道异面角、驶离轨道坡度角，以支持火星车驶离决策，保证火星车能成功驶离着陆平台、到达火星表面并开展后续巡视探测任务。

基于导航相机图像计算驶离轨道异面角及驶离轨道坡度角的工作包含相机标定、立体匹配、点云获取、平面拟合等多项内容。关于深空探测器导航相机标定问题，国内外学者提出了很多方法，Gennery^[1-2]先后提出CAHVOR和CAHVORE相机模型。Zhang^[3]提出一种通过立体像对对极几何的探测车自标定算法。在嫦娥3号、4号任务中，文献[4-6]高精度地完成了月球车立体视觉系统的几何参数标定。Ma等^[7]完成了“玉兔2号”月球车的立体视觉系统在轨标定。在点云及特征点坐标获取环节，采用前方交会方法，计算得到目标地形、特征点在导航相机坐标系下的坐标^[8]。在平面拟合环节，孔利等^[9]提出结合RANSAC方法与特征值法的平面拟合方法。本文使用总体最小二乘算法^[10]拟合出2驶离轨道平面，提出结合RANSAC算法与总体最小二乘算法的平面拟合方法拟合出火星表面平面，并进一步计算平面间夹角。在天问1号火星探测任务中，利用上述方法基于导航相机图像对着陆平台驶离轨道及火星表面进行了平面拟合与角度计算，为火星车驶离提供了决策支持。

1 火星车在轨驶离决策方法流程

天问1号着陆巡视器着陆后，为使祝融号火星车安全驶离着陆平台，需确定驶离轨道坡度角与驶离轨道异面角。先构建地面标定实验场，通过严密的控制场相机标定方法获得导航相机的内参(焦距、主点位置及镜头畸变)和外参(左右导航相机相对位姿参数、导航相机相对于火星车的位姿参数)。在火星车着陆火星表面之后，导航相机对火星表面及驶离轨道区域拍摄立体图像，然后对立体图像进行核线矫正、立体匹配，选取特征点并利用前方交会算法获左导航相机坐标系下的点云坐标及特征点坐标。之后结合相机外方位元素，将上述点云坐标及特征点坐标纳入火星车本体坐标系下。再利用上述数据分别拟合出火星表面和2驶离坡道平面，并计算驶离轨道异面角与驶离轨道坡度角。最后将计算结果与工程值比对用于支持驶离决策。具体流程如图 2所示。

1.1 火星车导航相机几何标定

导航相机的高精度几何标定，提供相机内方位元素、外方位元素以及畸变系数。首先，构建地面标定试验场，通过2台经纬仪，将控制场各控制点坐标归纳在火星车本体系下；利用经纬仪进行精密的控制点坐标量测后，可获得导航相机系与试验场系的相对位姿关系，将此关系作为导航相机外参数的初始值；然后，控制火星车俯仰轴和偏航轴转动，导航相机在不同位姿对标定试验场拍摄，在经过立体匹配后量测出特征点的像平面坐标；最后，将控制点三维坐标、像平面坐标、相机畸变参数和内参数作观测量，将控制点坐标、相机外参数、相机内参数和像平面坐标误差作未知量建立共线方程式，从而求解出各未知量。依据左右像片外方位元素即可间接计算火星车导航相机的相对方位元素。

1.2 立体图像的密集匹配

首先利用相机标定所得的导航相机相对方位元素进行立体图像核线校正, 并使用Sobel滤波^[11]在校正后的立体图像对之间获取可靠的支撑点并获取支撑点的视差值，再利用支撑点构建德劳内三角网，然后为获取更加密集的匹配点视差图，利用各支撑点的视差值在三角网中进行插值。同时，也可通过人工检查的方法选择符合前方交会要求的少量同名像点。

1.3 三维点云及特征点坐标

利用获取到的密集视差图，通过前方交会获得在左导航相机坐标系下的三维点云及少量驶离轨道面上特征点坐标。然后通过左导航相机外方位元素，将点云坐标转换到火星车本体系下，完成点云生成及少量驶离轨道面上特征点的三维坐标获取。

1.4 平面拟合及角度计算

本次火星探测任务所选择的着陆点为火星乌托邦平原南部，乌托邦平原是被火山熔岩填平的低矮平原，地形平缓，陨石坑较少。所以在计算着陆平台驶离轨道坡度角时，可以将小范围火星表面近似视为平面。

驶离轨道为工业制品，其表面平整，选取特征点后采用基于总体最小二乘算法的平面拟合方法即可拟合出轨道平面；火星表面虽相对平整，但不规则分布着散落的岩石块，为避免岩石块对平面拟合结果的影响，采用结合RANSAC算法与总体最小二乘算法的平面拟合方法基于点云数据进行平面拟合。

之后计算各平面法向量间夹角即为本文所需的轨道异面角和驶离轨道坡度角。最后，将计算得到的角度值与工程值进行对比，若计算角度均小于工程值，则火星车可以执行驶离动作。

1.4.1 基于总体最小二乘算法的平面拟合

总体最小二乘法是一种较先进的最小二乘法^[12]，总体最小二乘法认为回归矩阵存在干扰，在计算最小二乘解时考虑了这个因素，而在一般的最小二乘法中没有考虑该因素的影响。总体最小二乘法应用广泛，以下是基于该方法的平面拟合基本原理。

设平面方程为

$z=a x+b y+c .$

(1)

对同平面上的n个数据点{(x_i, y_i, z_i)，i=1, 2，…，n}。考虑到x, y, z这3个方向上的误差值，将式(1)改写为

$z+v_z=a\left(x+v_x\right)+b\left(y+v_y\right)+c .$

(2)

在式(2)中，v_x, v_y, v_z分别为x, y, z这3个方向上的误差改正数。

整理式(2)，则有

$\left(\boldsymbol{A}+\boldsymbol{E}_A\right) \boldsymbol{X}=\boldsymbol{L}+\boldsymbol{E}_L .$

(3)

其中：

$ \begin{gathered} \boldsymbol{A}=\left[\begin{array}{ccc} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ \vdots & \vdots & \vdots \\ x_n & y_n & 1 \end{array}\right], \boldsymbol{E}_A=\left[\begin{array}{ccc} v_{x_1} & v_{y_1} & 1 \\ v_{x_2} & v_{y_2} & 1 \\ \vdots & \vdots & \vdots \\ v_{x_n} & v_{y_n} & 1 \end{array}\right], \\ \boldsymbol{X}=\left[\begin{array}{c} a \\ b \\ c \end{array}\right], \boldsymbol{L}=\left[\begin{array}{c} z_1 \\ z_2 \\ \vdots \\ z_n \end{array}\right], \boldsymbol{E}_L=\left[\begin{array}{c} v_{z_1} \\ v_{z_2} \\ \vdots \\ v_{z_n} \end{array}\right] . \end{gathered} $

之后用矩阵奇异值分法进行解算，即

$\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{A}}&\mathit{\boldsymbol{L}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathop {{\boldsymbol{U}_1}}\limits_{m + 1} }&{\mathop {{\boldsymbol{U}_2}}\limits_{n - m - 1} } \end{array}} \right]\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}\\ \mathit{\pmb{0}} \end{array}} \right]{\mathit{\boldsymbol{V}}^{\rm{T}}} = {\mathit{\boldsymbol{U}}_1}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{T}}}.$

(4)

其中

$ \begin{gathered} \boldsymbol{V}^{\mathrm{T}}=\left[\begin{array}{ll} \boldsymbol{V}_{11} & \boldsymbol{V}_{12} \\ {\mathop {{\boldsymbol{V}_{21}}}\limits_{m} }& {\mathop {{\boldsymbol{V}_{22}}}\limits_{1} } \end{array}\right]\begin{array}{*{20}{c}} m\\ 1 \end{array}, \boldsymbol{U}_1=\left[\begin{array}{cc} {\mathop {{\boldsymbol{U}_{11}}}\limits_{m} }& {\mathop {{\boldsymbol{U}_{12}}}\limits_{1} } \\ \end{array}\right], \\ {{\mathit{\boldsymbol{ \boldsymbol{\varSigma}}}}}=\left[\begin{array}{cc} \mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_1 & \mathit{\pmb{0}} \\ {\mathop {{\mathit{\pmb{0}}}}\limits_{m} }& {\mathop {{\boldsymbol{\varSigma}_{2}}}\limits_{1} } \end{array}\right]\begin{array}{*{20}{c}} m\\ 1 \end{array}. \end{gathered} $

由此得到模型参数估值为

$\boldsymbol{X}=-\boldsymbol{V}_{12} \boldsymbol{V}_{22}^{-1}.$

(5)

残差矩阵为

$\left[\begin{array}{ll} \boldsymbol{E}_A & \boldsymbol{E}_L \end{array}\right]=\boldsymbol{U}_{12} \mathit{\boldsymbol{ \boldsymbol{\varSigma}}}\left[\begin{array}{ll} \boldsymbol{V}_{12}^{\mathrm{T}} & \boldsymbol{V}_{22}^{\mathrm{T}} \end{array}\right].$

(6)

本文方法用如下公式直接计算出参数：

$\boldsymbol{X}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}-{\&}^2 \boldsymbol{I}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{L} .$

(7)

式中：& 为增广矩阵[A L]的最小特征值。

1.4.2 结合RANSAC方法与总体最小二乘平差法的平面拟合方法

RANSAC方法是由Fischler和Bolles^[13]提出，通过计算数据集中每一个点到随机平面的距离，划分出内点与外点，再通过不断选择新的随机平面来最大限度地保留内点，从而消除异常点的影响。其优点是对于包含较多异常值的数据集鲁棒性高，但没有考虑到测点误差。而总体最小二乘法却注意到了3个方向上可能存在的误差，所以相对于传统平面拟合方法，将随机采样一致性算法与总体最小二乘法相结合，将会使平面拟合结果拥有更高的鲁棒性及准确性，算法流程如下：

1) 从点云数据中随机抽取3点，判断是否共线，若是则重新抽取，若不是则给出这3个点所在的平面方程。

2) 计算数据集中剩余各点到上述平面的距离及标准偏差。标准偏差是一种数据值偏离算术平均值程度的度量，利用点到模型平面距离的标准偏差作为判断局内外点的准则，可以有效区分出局内点与局外点。

3) 设置合适的阈值，当标准偏差大于阈值时，该点为局外点；当标准偏差小于阈值时，该点为局内点，并统计局内点数量。

4) 重复上述步骤，比较选取局内点个数最多的数据集。

5) 采用基于总体最小二乘法的平面拟合方法对以上步骤选取出的数据集中的内点进行平面拟合，得到平面方程参数值。

在以上过程中，随机抽取次数可以根据数据集的大小由魏英姿和刘晓莉^[14]提出的方法计算得到。

2 实验验证 2.1 点云拟合平面仿真实验

设平面方程为

$ z=2 x+y+0 . $

本文使用matlab模拟点云数据进行平面拟合，从设定平面中选取300个随机点。为验证本文方法的有效性及不同异常值数量情况下的鲁棒性，分别加入100、500、800和1 000个异常值形成4组模拟点云数据。采用本文方法，结合随机采样一致性算法与总体最小二乘算法的平面拟合方法对上述模拟点云数据进行平面拟合，计算出估算值与真实值间的标准误差。

当点集中分别存在不同数量的异常点时，仿真数据及分别使用2种方法的平面拟合结果如图 3所示，得到的平面参数及标准偏差如表 1所示。

从表 1可以看出，当模拟点云中没有异常值时，2种方法均可以准确地估计出平面参数，当模拟数据中加入异常值时，随着异常值数量的增加，整体最小二乘平差法更容易受到异常值的干扰从而不能很好地估计出平面参数，而本文提出的结合了RANSAC算法与总体最小二乘平差法的平面拟合方法估计出的平面参数依旧可以保持良好的准确性，标准偏差均小于0.04，由此可知本文方法在平面拟合中的有效性及鲁棒性均优。

2.2 火星车在轨驶离支持 2.2.1 导航相机几何标定

在地面实验中，布置高精度的标定实验控制场。本次地面实验的三维控制场环境如图 4所示，在距离导航相机1.3~5.0 m范围内共布设33根测杆、121个控制点。导航相机标定结果如表 2及表 3所示。其中，镜头畸变模型由径向畸变参数k₁、k₂和切向畸变参数p₁、p₂组成，旋转矩阵是描述导航A相机系到导航B相机系下的旋转，位移向量是导航A相机系原点在导航B相机系下的三维坐标。

2.2.2 点云及特征点坐标获取

火星车及着陆平台着陆火星表面后，桅杆展开，导航相机获取驶离坡道及火星表面图像如图 5所示。利用导航相机标定结果，对获取到的立体图像对进行密集匹配，通过前方交会获得离散点云及特征点标。

2.2.3 平面拟合及角度计算

利用本文提出结合随机采样一致性算法和总体最小二乘算法的平面拟合方法对点云数据进行平面拟合。通过驶离坡道特征点，利用总体最小二乘算法估算出2驶离坡道平面参数。最后，通过平面法向量计算各平面间夹角，结果如表 4所示。

如表 4所示，2驶离坡道异面角为0.203°，远小于工程值(10°)，其与火星表面夹角分别为-18.947°与-19.154°，均远小于工程值(34°)，故判定火星车驶离着陆平台动作安全性高，无脱轨、倾翻危险，可执行驶离指令。

3 结语

本文以我国天问1号火星探测任务中的火星车在轨驶离着陆平台决策问题作为研究对象，提出基于导航相机图像的平面间角度计算方法，创新性地采用RANSAC算法与总体最小二乘算法，并利用仿真数据进行算法有效性、鲁棒性的验证。最终将此方法运用在火星探测任务中并圆满完成火星车驶离决策支持，进一步验证了本文方法的有效性。