李阳平 黄 玲 王 珂 赵海峰*

①(中国科学院大学 北京 100086)

②(中国科学院空间应用工程与技术中心中国科学院太空应用重点实验室 北京 100094)

③(中国科学院空天信息创新研究院 北京 100094)

1 引言

月球是与地球最近的天体，但时至今日，人类对月球表面月壤结构的勘探能力和成果仍十分有限。在美国Apollo任务中[1]，宇航员仅在近地表提取了月壤与月岩样本，月球更深部的原状结构信息和深部样品无法获取。我国嫦娥探月工程项目在月面上采用了探地雷达系统进行探测[2–4]，科学研究结果表明探地雷达系统不仅能够对月球地下介质分布开展有效的探测，也能为月面深部钻采的场地优选和路径规划提供依据，增强月面钻探的能力，降低钻探取样的风险。由于探地雷达的高频电磁波在月壤及月岩砂砾介质中具有很强的散射特征和各向异性特征，探测回波复杂，并且数据有限，采用传统的数据处理与目标成像技术难以对月壤中的月岩体形态进行精确成像，无法高效地适用于未来月球探测任务中。本文针对月面浅层(5～10 m深)钻进采样的科学需求，在雷达探测数据处理引入机器学习算法，增强探月雷达探测数据实时地层成像的效率，最终为月球探测任务的月岩钻进避障和钻进策略优化提供一种可行的解决方案。

由于近年来计算机计算性能的提高，机器学习算法的爆发性发展，人工智能神经网络算法相对于传统的模型反演算法展现出了巨大的发展潜力和优势。研究成果表明人工智能神经网络在探地雷达信号处理中具有较高效准确的优势[5–7]。深度神经网络由输入层、隐藏层和输出层组成，输入数据的特征变量作为多变量线性组合被传递给隐藏层的节点，隐藏层的每个节点设置非线性激活函数，通过每层的线性组合传递至输出层，从而实现神经网络对复杂非线性关系的建模。

随着深度神经模型的不断优化和扩展，已有学者将深度学习应用于探地雷达回波信号分类。Gao等人[8]利用基于区域的快速卷积网络方法(faster R-ConvNet)进行路面缺陷检测，探测精度可达89.13%。Li等人[9]利用快速目标检测(You Only Look Once, YOLO)模型检测沥青路面的隐藏裂隙，并在小数据集上也取得良好效果。2020年，Khudoyaro等人[10]将深度3维卷积神经网络应用于3维探地雷达数据，识别地下物体，将城市道路分类为管道、孔洞、检修孔和土壤，其准确率达到了97%。然而，当前的研究方法需要花费大量人力进行特征标定，而对于探地雷达图像识别而言，目标的精确几何特征仍无法定量描述，因此只能将其输出特征分类为不同的地质结构，分类结果单一且粗糙。

为了获得更准确的掩埋物体特征，Liu等人[11]提出基于全波形反演的检测算法，能够准确恢复地下掩埋的圆柱体位置和直径，但是麦克斯韦方程组的计算需要消耗大量资源，无法实时获得结果。于是，Giannakis等人[12]将深度学习与探地雷达(Ground-Penetrating Radar, GPR)信号处理结合，实时预测掩埋物体的深度，但目前仅支持对规则圆形进行处理，无法针对复杂表面预测，地质环境也局限于混凝土与钢筋结构，难以处理月球环境。

因此，本文结合深度神经网络在非线性问题中的优势，搭建了一种能分析预测不规则介质上表面几何特征与位置的深度学习模型。首先利用边缘检测等数字图像处理方法获取月岩轮廓，构建月下地质结构环境，模拟仿真得到探地雷达回波信号数据集；其次，利用主成分分析法提取雷达回波数据的主要特征，提取后的特征能够加强信号的表现力，作为神经网络模型的输入数据；最后，根据反向传播算法，搭建了一种深度神经网络模型，在训练过程中自主调整神经元节点权重，实现对掩埋介质的上表面几何特征与位置的回归预测。

与传统探地雷达数据处理方法相比，本文所提模型有以下优点：(1)本文针对月球环境，采用数字图像处理方法获取掩埋介质轮廓，并构建大量接近真实几何特征的月球地质模型；(2)本文采用的深度神经网络模型，其中隐藏层神经元节点的激活函数具有非线性特征，从而对非线性数据处理具有较好的预测结果；(3)本文对于雷达信号进行分类，利用训练完成后的数字化成像模型，使用时直接输入回波信号数据获得成像结果，计算消耗小，耗时低，能够满足未来航天任务设计的实时性要求。

2 基于数字图像的月球地质环境搭建方法

在探地雷达的人工智能神经网络训练中，训练数据通常来自两种情况：(1)相对有限的一组探地雷达真实数据；(2)基于简单规则形状模型而获取的仿真模拟数据。这些数据有如下缺点：(1)真实数据的数据量不足、环境或系统影响因素多；(2)简单模型获得的模拟仿真数据不能反映实际情况。在当前更快速的全波形正向求解器和改进的分析方法背景下，如何获得更逼真的月壤模型模拟数据是解决问题的关键。本文根据月面环境，基于电磁波在地下的传播特性，结合嫦娥3号和4号搭载的月球探测器参数，以及Apollo系列任务返回的月壤与月岩样本图像，构建了一种接近真实的月面地质环境的建模和回波数据分类方法。

2.1 探月雷达参数设计

探月雷达的工作原理是发射天线将超宽带的毫微秒电磁波脉冲向月面辐射，当遇到电性特征不同的物质(如熔岩层和溅射物)时，电磁波将会发生反射和散射现象。探月雷达的接收天线将接收到的信号进行放大、采样，随后对获得的探测信号进行处理、反演成像，最终可获得月球地下介质的结构分布。

为深入探究不同深度下复杂月球地质结构特征与探月雷达回波信号的关系，本文首先针对月球浅表层0.5 m深单个掩埋非规则形状月岩块体进行仿真模拟，取得预测结果后，再进行10 m深月球复杂地质环境模拟。其月面雷达探测系统的主要技术指标如表1所示。

表1 探月雷达的参数设计

针对月面下0.5 m的浅层月面雷达穿透深度，设计月面雷达的中心频率为1.2 GHz。在月球环境条件下，其探测距离向分辨率约为2 cm，该月面雷达系统可用于探测月面表层掩埋月岩的轮廓和分布位置。针对10 m深的探测目标，设计月面雷达中心频率为700 MHz。在月壤条件下，其探测距离向分辨率约为20 cm，该月面雷达系统可用于探测月壤厚度及月岩分层结构。

2.2 构建2维无规则介质正演模型

在正演仿真模型中的月壤中月岩介质物性参数和几何参数设计，则采用基于Apollo 14任务返回的岩石样本作为参考模型[13]，从而可使得仿真模型和月岩表面形态更接近真实的月面条件与环境，能更有效地检验方法的有效性。其岩石样本部分照片如图1所示。

图1 Apollo 14号返回样品图像

构建几何参数的过程中，首先通过提取真实月岩照片的月岩轮廓来构建2维gprMax仿真模型[14]，从而可精确地模拟月岩样品表面较为真实的粗糙形态特征。在2维月岩图像中，图像存储的大部分信息主要由边缘组成，表现为局部特征不连续，可将其理解为图像灰度变化剧烈的区域，因此利用梯度算子来进行边缘检测是效果较好的一种方式。利用Canny梯度算子[15]对图像进行边缘检测后，获得结果如图2(a)所示。

可以看出，图像边缘的噪声较多，中间存在许多的空穴，其与真实月岩表面形态仍存在较大的差别，因此采用膨胀腐蚀算法进行处理，然后进行孔洞填充，最后得到与真实月岩表面形态更为接近的月岩截面轮廓，如图2(b)所示。

图2 获取月岩轮廓

在获取月岩截面图像后，利用gprMax软件进行2维有限差分数值模拟。gprMax软件是一个模拟电磁波传播的开源软件[14]。它基于Yee算法，通过有限时域差分方法求解麦克斯韦方程组。在模型的构建中，还能导入便携式网络图形(Protabel Network Graphics, PNG)格式的图片，通过图片中不同颜色像素来定义其材料特征，由此可以构建出月岩和月壤二相介质的模型。

在使用gprMax进行正演模拟中，首先需要确定仿真模型区域，并划分网格大小。本文所使用的环境设定为：CPU Intel Xeon Gold 5118, 256 GB内存，Tesla V100 16 GB显卡的服务器，由于数据量和计算时间的限制，本文将浅层2维月球地质模型区域设计在0.7 m×0.5 m的空间中，网格大小为0.005 m，然后将岩石位置随机分布在空间中的任意坐标。其中一个代表性模型如图3所示。

图3 模拟二相不规则月岩环境

在此模型中(图3)，月岩体底部距月壤表面的埋深为0.25 m。波形的中心频率设置为1.2 GHz。而时间窗口设置为7.5 ns。雷达的发射天线与接收天线的间距为2 cm，测线上的测点间距为1 cm，即雷达自模型最左侧开始行进，每隔1 cm进行电磁波的发射与接收，从而获得雷达探测的B-Scan数据。

为进一步研究含有多月岩块体地层的回波信号特征与数据处理方法，本文将2维月球地质模型区域设置为10 m×20 m。在模型中随机掩埋不规则形态和大小的月岩介质，生成不同岩石含量的地质环境，其中两个代表性模型如图4所示。

图4 多目标非均匀月岩地层

在此模拟环境中，模拟区域设置为20 m×10 m，网格大小为0.05 m，岩石大小分为3种：10 cm左右、30 cm左右和50 cm左右。控制大中小石头比例，均匀地将它们散布在月壤地层中。而当月岩个数逐渐增多时，需要控制岩石面积百分含量的范围。本文按照黑色像素点/白色像素点的比值，来计算岩石含量。

2.3 月壤与月岩的电性特征参数设置

构建月球地质环境的过程中，需要考虑的月岩电性特征主要包括介电常数、电导率、磁导率、磁损耗。Apollo系列任务返回的月壤与月岩样本均进行了介电常数、密度的测量。采用的月壤参数如表2，月壤参数基于450MHz雷达下测量[16]。

表2 月壤电性参数表

不规则形态月岩的电性参数设置如表3所示。

表3 月岩电性参数表

研究表明，介电常数的实数部分大约与频率无关，在1 MHz下的探测结果与450 MHz的结果非常接近。而损耗正切角通常是频率的复杂函数，不同频率下的结果各不相同[16]。因此，本实验中，只模拟了月岩与月壤的介电常数，而并未模拟损耗正切值。

3 基于深度神经网络的探月雷达反演算法

本文基于模拟月岩的几何与形态参数进行gprMax仿真模拟。在获得大量雷达信号数据集的前提下，提取数据中的主要特征，作为神经网络的输入数据。最后采用反向传播算法(Back Propagation, BP)，构建深度神经网络，对地质结构特征进行预测，绘制表面轮廓图像。

3.1 输入数据特征提取

神经网络模型的输入数据为探地雷达回波数据，与其对应的输出数据为此收发射天线垂直界面岩石的上表面各像素位置。然而，0.5 m深雷达单道数据的回波量为339维，如果全部作为深度学习的输入数据，则数据量过大，导致模型泛化性能差。因此需要提取雷达回波数据的主要特征。主成分分析(Principal Component Analysis, PCA)是机器学习研究中的一种常用的降维方法[17]。主成分分析法提取特征时构建一个n×m的变换矩阵，该矩阵能够将样本x(即339维的单道回波数据)映射到一个新的特征子空间中，并使得变换后的特征子空间维数m少于原始特征空间维数n，其主要思路为

在本文的研究中，解释的特征值比率累计和在95%以内，即尽可能降低信息损失。通过分析对比多参数、多模型条件下的训练结果，选择了最合适的特征维数，将0.5 m地质模型特征的339维回波数据降为39维主特征。

3.2 神经网络模型结构

在完成雷达观测数据的降维压缩之后，本文利用反向传播算法的思想，搭建神经网络模型。模型如图5所示。

图5 神经网络模型示意图

输入层有N个训练样本，每个训练样本中包含回波数据经过主成分分析后的39维主特征，对于训练样本集X=[x1,x2,···,xk,···,xn]来说，预测输出为Y=[y1,y2,···,yk,···,yn]，目标输出则为T=[t1,t2,···,tk,···,tn]。初始权重ω的选取则一般采用统计学上随机抽取的思想，在算法中随机产生大量的初始值，逐次进行迭代优化。在BP算法中，由于函数大多输出值在[–1,1]，因而初始权重选择一般也在这个范围内。采用逐步搜索法来解决初始权值选择的随机性。先将初始权值区域N等分，在N个区域内再重复上述步骤，直到误差函数不再减小，认为找到了最优解。只要区域选择得足够小，这种方法能够有效避免局部极小问题。算法核心就是将复杂的雷达回波信号识别转化为回归问题，从而预测不规则介质的上表面几何特征与位置。

其模型的主要参数如表4所示。

表4 神经网络模型结构参数

数据输入输入层后，需要衡量整个模型是否找到了满足误差最小的最优解，总误差采用MSE函数，其表达式为

在本文中，经过模拟实验，隐藏层选择过少会导致模型拟合度降低，而隐藏层过多时模型泛化能力较差，因此隐藏层被确定为8层。其中，第1层与第2层的激活函数为ReLU[18]，神经元共64个。第3层和第4层激活函数为Sigmoid[19]，神经元共64个。接下来第5层～第8层的激活函数都为ReLU，而第5层与第6层的神经元个数为32个，第7与第8层的神经元个数为16个。其结构如表5所示。

表5 神经网络结构

3.3 实验结果

根据探地雷达电磁波的工作方式，在前期实验中，为了研究神经网络算法对月岩几何形状预测的效果，因此首先将不规则月岩形状掩埋物体的介电常数设置为与月壤差异较大的数值，以便逐步深入研究方法的可行性和效果。利用gprMax对图3所示的模型结构进行仿真，掩埋的月岩上表面深度为0.11，由此得到一道A-Scan，与多道A-Scan组成的截面图像B-Scan。其结果如图6所示。

图6 不规则月岩雷达回波信号

数据样本的丰富度影响模型是否能预测此条件下的任意月岩形态，本文首先随机生成了8000个模拟月岩环境，由于数据量不足，因此训练后的神经网络模型预测精度非常低。增加数据量至26000个模拟月岩地质环境。图7是不同数据量对训练预测结果的影响。

图7 不同数据量对训练结果影响

获取足够的数据量后，利用人工神经网络开始进行学习预测，模型的评价标准则由R-square确定系数[20]来确定。其表达式为

其中，y1, y2,···,yn表示测量值，f1, f2,···,fn表示预测值。R-square取值在0～1，越接近1时则表明回归模型表现效果越好。由此确定系数，得到较大差异介电常数不规则形状掩埋月岩体上表面的位置预测结果，R-square为0.93，如图8所示。

图8 较大差异介电常数不规则形状掩埋月岩体上表面轮廓的预测对比图

由此可见，此方法预测的上表面较为清晰，误差较小，误差范围一般在10 cm以内。最后，从预测结果中随机选择3个地质环境，对其上表面轮廓进行绘制，其结果如图9所示。

图9 较大差异介电常数不规则形状掩埋月岩体上表面轮廓预测结果的误差直方图

图10(a)的上表面形态预测准确率为93.76%，图10(b)的上表面预测准确率为96.04%，而图10(c)的上表面预测准确率是95.13%左右。

图10 较大差异介电常数不规则形状掩埋月岩体上表面轮廓预测示例

在完成较大差异介电常数不规则形状掩埋月岩体模型的仿真与成像研究，并取得良好效果之后，进一步将模型月岩介电常数介质材料变换为表3所示(Apollo号返回的月岩样本材料数据资料)。模拟仿真如图11所示。

由图11可知，月岩样本材料掩埋介质的上表面位置预测的R-square为0.78，与大介电常数掩埋介质的预测准确率有一定差别。

图11 月岩材料不规则月岩上表面轮廓预测

在获取单个岩石的上表面轮廓预测结果后，针对含有多月岩块体地层的模型(对应图4(a)与图4(b))，其回波信号B-Scan图像如图12所示。

由于掩埋月岩在深度方向上互相交叠，出现反射、散射、多次回波等现象，因此回波信号十分杂乱。且回波数据特征空间维度较高，利用PCA提取主要特征后，降维的输入数据为56维。岩石含量为38%的预测结果如图12(b)所示，深度学习模型的预测准确率R-square值仅为0.48。

图12 含有多月岩块体地层模型雷达回波信号的B-Scan多道雷达回波数据图像

3.4 实验结果分析

对月岩上表面轮廓进行模型和预测结果对比分析，可见0.5 m深的大介电常数模型可以看出大部分区域拟合得较好，误差大约在10 cm以内，但下表面轮廓由于预测精度较低，因此无法进行拟合。其原因是： (1)受上表面强反射影响，下表面的反射信号回波幅值弱；(2)部分下表面的回波信号与上表面的2次回波信号重叠，难以区分上下表面回波信号，这也是探地雷达数据处理的共性问题。采用本文所提方法，其优势在于，在将模型训练好之后，只需要约为0.0032 s的时间即可获得预测结果，极大提高了数据处理效率，达到现场处理与解释的能力。分析结果还表明当月岩的介电常数与月壤接近时，预测精度有所下降，但仍可预测其表面轮廓的形态趋势。

如图13所示，在多月岩块体地层模型的仿真与模型形态预测中，雷达仿真数据复杂程度增大，回波数据包含多目标散射、多次回波信息。月岩表面的模型形态与预测结果形态之间的匹配关系下降较为严重，但其仍显示出一定的线性相关性，说明预测方法仍能反映出月岩表面形态的主要特征。为进一步提高预测精度，需要在后续研究中首先从大量回波数据中提取反射信息，压制目标间的多次反射与回波，之后再进行主成分分析，提取信号主要特征。另外，有限训练数据集无法满足复杂模型的回归拟合优化需求，需要增大多月岩复杂模型的仿真训练数据多样性和数据量，以使得模型泛化能力强，增强预测的稳定性。

图13 含有多月岩块体地层模型上表面轮廓预测结果

4 结束语

本文基于深度学习方法的优势，针对月壤下月岩形态探测实时处理的需求，将探地雷达回波信号的解释和几何重建问题转化成深度学习中的回归问题，并通过仿真分析来验证该思路的可行性。

(1)本文通过提取Apollo 14号返回的月岩的2维轮廓，构建月岩形状的2维模型，并根据Apollo返回样品的电性参数来设置模型的电学特性，从而更真实地模拟了月岩结构的真实形态。

(2)利用gprMax软件构建探月雷达观测数据集，针对获得的高维回波数据，采用主成分分析方法来降低模型输入数据的复杂程度，防止过拟合现象。然后，构建神经网络模型对月岩几何形态进行直接预测。

(3)在基于2维模拟数据建立的神经网络模型预测中，当月壤中掩埋的月岩具有较大介电常数材料时，上表面的预测准确率达到0.93，能较好地重构出上表面轮廓，与Giannakis等人[12]的混凝土钢筋规则圆形介质预测实验相比，本文工作已能够对不规则介质的上表面进行有效的预测重构；利用已训练完成的模型进行数据处理，预测时间约为0.003 s，能够达到实时性预测的要求。当月岩样本介电常数材料接近月壤时，上表面的预测准确率为0.78，表明月壤和月岩介电常数对比度较小时，模型预测精度下降。

(4)在含复杂多月岩块体分布的地层模型上表面轮廓预测中，由于回波信号复杂，深度学习模型的预测准确率仅为0.48，其原因可能是复杂月岩分布导致雷达信号散射叠加，在后续工作中，可以研究不同掩埋介质深度、不同岩石含量下的模型拟合效果，并对比其他神经网络方法，选择更优化的数学模型，提高预测准确率。

在后续研究中，需要对更多地质介质的物性参数数据库进行统计和分类，建立完整的基于数据驱动的神经网络地层结构成像方法的数据集和理论体系。此方法对地球或者地外地质探测工作有着非常大的应用价值和科学意义。