0 引言
月球探测是人类进行太空探测的开端，太空探测可加深人类对地球、月球和太阳的认识，带动一系列基础科学的创新[1-2]。在月球探测技术中，探月雷达是其中一种重要方法，也是探地雷达的一种。由于月球表面没有液态水，介质的电阻率很高，电磁波穿透能力强，雷达可以探测到月面以下数百米深的次表层[3]。雷达不仅能够实现月表的巡视探测，而且能在地球基地和绕月卫星上对月球表面进行大范围遥感测绘[4]。因此，在月球、火星等其他外星球探测中，雷达已经成为一种用于浅表层和次表层地质结构探测的有效工具，在现阶段的月球探测中发挥着重要作用[1]。
我国发射的嫦娥五号探测器已于2020年12月31日于月球表面着陆，嫦娥五号探测器负责嫦娥三期工程“采样返回”任务，为我国实施载人登月做技术储备[5-6]。嫦娥五号月壤结构雷达探测仪采用非屏蔽的高频天线，其工作环境相比嫦娥三号更加复杂，天线耦合、着陆器地板和其他金属构件的强反射以及底板与地面之间的多次反射信号等都会对月面以下目标的反射信号形成强烈干扰，甚至淹没目标信号[7]。根据嫦娥五号的任务安排，在返回雷达数据后，需要快速地对复杂探月雷达信号进行数据处理并实现高分辨率成像，及时为月壤钻取任务探测钻头下方区域2 m深度范围内的月壤结构和月岩分布。因此，应进一步研究数据处理方法，抑制杂波和多次波的干扰，才能保证数据解译的准确性；并进一步研究适合强电磁干扰环境下的探月雷达快速高精度成像算法，对算法保证计算效率的同时，兼顾成像精度。
基于此，本文开展月壤分层结构的建模、电磁响应计算与分析以及嫦娥五号探月雷达数据的快速高精度成像研究，并通过地面试验进行验证，为月壤、月岩物理参数的准确分析提供基础，为嫦娥五号探测器的月壤取样提供科学依据，也为进一步的月球以及深空探测提供方法和技术。
1 嫦娥五号月壤结构探测仪
嫦娥五号着陆器上搭载一套由12个Vivaldi天线组成的多发多收探月雷达系统，如图1所示。该雷达系统由中国科学院空天信息研究院团队研制，其中心频率为2 GHz。嫦娥五号探月雷达将在静止的状态下对着陆点下2 m深度范围内的月壤结构进行精细探测，为后续钻取采样过程的顺利实施提供月壤分层结构及可能存在的月岩等信息，并通过分析现场探测数据与返回样品实验室分析数据之间的联系，为研究月壤形成的原因和演化机理提供数据支撑。
图1 嫦娥五号月壤结构探测仪立体结构示意图(蓝色区域表示雷达探测仪天线和月壤钻头的布局)
由于嫦娥五号探测器所搭载的载荷众多、系统复杂，内部电子器件和金属构件对探月雷达数据干扰严重。此外，由于着陆器整体优化的需要，配置给雷达系统的空间有限，因此天线设计采用尺寸较小的空气耦合Vivaldi天线，天线布局方式如图2所示。由于雷达天线距月球表面上方仅有90 cm(见图3)，月表以下的目标反射信号受到天线耦合信号、着陆器金属底板和其他构件反射信号，以及地面多次反射信号等的严重干扰，信噪比低。
(a) 立面图
(b) 平面图
图2 嫦娥五号月壤结构探测仪天线布置图(单位：mm)
图3 嫦娥五号雷达探测仪的复杂干扰环境示意图
考虑到嫦娥五号着陆器底板和其他金属构件等对月壤结构反射信号的严重干扰，有必要对复杂环境下探月雷达的数据处理方法进行优化，研究适合嫦娥五号月壤结构雷达探测仪的数据处理流程，在复杂干扰环境中提取有效的弱目标信号，为月表实测雷达数据的正确解译提供依据。此外，根据嫦娥五号的任务安排，从探月雷达数据传回地球到完成月壤采样只有短短3个小时，给探月雷达数据的快速处理与成像提出了很大挑战。因此，研究探月雷达数据处理和快速成像方法，是成功实现嫦娥五号月壤结构雷达探测仪科学探测任务的关键问题之一。
2 探月雷达数据处理方法
2.1 系统校准
由于嫦娥五号月壤结构探测仪的数据成像是通过发射天线辐射波场和接收天线的逆时波场的互相关进行成像。在实际仿真模拟过程中，如果将天线甚至整个着陆器模型纳入到仿真模型中，将大大增大计算量，计算时间无法承受。实际计算过程需要在仿真模型中将天线近似成点源，这需要对系统进行精确校准满足近似条件。这些系统校准包括子波测试、点扫描组件及电缆延时校正、天线相位中心校准和背景干扰去除。
子波是指由安装在着陆器上的雷达天线实际辐射到月壤介质中的电磁波波形，并用于模拟发射天线的辐射波场。本文采用如图4所示的实验方法采集各个发射天线的子波，将整个着陆器在实验室内吊高至离地面约3 m高的位置，在地表放置吸波材料，然后在着陆器钻头的正下方架设一个接收天线，采集不同发射天线向下辐射的电磁波。探月雷达系统的天线阵列单元采用的Vivaldi天线，Vivaldi是行波天线。传输的激励信号从馈点沿着天线元件传播时会被延迟。在仿真计算中，应将天线简化为理想源点，该源点可看作球面辐射场的起源，通常该点被称作天线视在相位中心。
图4 天线子波测试实验示意图
月壤结构探测仪的天线阵列搭载在着陆器的底板上，而且着陆器大部分组件都是由金属构成，着陆器底板和其他金属组件将对天线辐射的电磁波产生很强的干扰，若不对这些背景干扰进行处理，来自月壤内部的目标回波信号将被完全淹没在这些背景干扰中，很可能造成成像失败。因此必须预先采集这些背景干扰，在数据处理过程中进行去除。在实验室内采用的背景干扰去除实验方法是将着陆器放置在探测区域的上方，采集包含着陆器的背景干扰、环境干扰和目标回波信号的数据，然后在地表放置一层20 cm厚的吸波材料，再次采集背景和环境干扰信号。两次数据相减，能够抑制背景和环境干扰信号，从而增强目标回波信号，提高信号干扰比。
在实际月壤表面环境中，通过实验室方式在月表铺设吸波材料测量背景干扰的方法是不现实的，因此将搭载月壤结构探测仪的着陆器放置在航天五院大型微波暗室以及在绕月飞行阶段开机采集背景干扰数据，储存后进行实测数据的背景去除处理。
2.2 月壤分层结构速度分析
系统校准后，通过速度谱分析法对嫦娥五号探月雷达探测器所采集的多偏移距数据进行分析，获取月壤分层结构不同月壤层的介电常数和厚度[8-9]。速度谱分析方法的基本原理如图5所示。图5(a)为地下单个水平反射界面进行多偏移距测量得到探月雷达剖面示意图，通过给定的试速度生成一系列的双曲线，然后沿着这些双曲线叠加不同偏移距通道的信号能量，可得到试速度与零偏移距双程走时的速度谱如图5(b)所示。由图5(b)可知：如果选择的速度非常接近真实的速度，则叠加的能量将在速度谱上有一个最大的峰值，可以选择相应的速度和零偏移距双程走时，并认为其为地下介质的电磁波波速和电磁波到达水平反射界面的双程走时，从而计算出水平反射界面的深度。
(a) 多偏移距探月雷达剖面图
(b) 速度谱示意图
图5 探月雷达速度谱计算示意图
如果地下介质中含有多个水平反射界面，在共中心点剖面上会出现多个反射波，则共中心点速度谱上出现多个能量峰值，根据速度谱上能量峰值提取的速度是相应的叠加速度。当收发天线距离很小时，叠加速度近似等于均方根速度，为了计算地下各层介质的层速度，可应用Dix公式将均方根速度转化为层速度。Dix公式(1)中第i层层间速度vint,i的计算表达式为
(1)
式中，vst,i+1为第i+1层的均方根速度，ti+1为从地表到i+1层的双程走时，vst,i为第i层的均方根速度，ti为地表到i层的双程走时。层速度是与地层岩性密切相关的速度，在各种速度中占有十分重要的地位，计算了各层的层速度之后，结合各层的双程到达时间，则各层的厚度都可计算出来，第i层介质的厚度为
di=vint,i(ti+1-ti)
(2)
2.3 逆时偏移成像
由于嫦娥五号的月壤结构探测仪固定在着陆器上，采用静止探测的工作方式,仅由其所装配的12个Vivaldi天线采集132道雷达数据，地下目标回波信息量很少，采用常规的基于路径追踪原理的常规雷达成像算法很难满足高分辨率成像需求。逆时偏移技术是近十年发展起来的一种最新的基于波动方程的偏移成像技术，已经在地震勘探领域取得了巨大的成功[10-11]。将逆时偏移技术应用于雷达成像，是基于以下的电磁波波动方程：
(3)
式中，
为电场强度，μ为磁导率，ε为介电常数。
相对其他偏移方法，逆时偏移算法采用双程电磁波动方程对波场进行延拓避免了对波动方程的近似。它没有倾角限制，并可以利用转换波、棱镜波或多次反射波进行成像，获得更精确的振幅等动力学信息，实现保幅成像，还可以更好地对复杂速度场进行更细化更精确的估计。逆时偏移的主要流程包括激励源波场的正演模拟、接收波场的逆时外推和成像条件的应用，其基本流程如图6所示。
图6 逆时偏移算法的基本流程
接收波场的逆时外推是波的逆时传播过程，沿时间轴的负方向递推，每完成一个时间步长的外推，则根据一定的成像条件，将该时刻的波场与预先存储的正演模拟中对应时刻的波场进行一次成像，得到的成像值累加到成像剖面，直到波场外推过程完成就可以得到最终的偏移成像剖面。常用的成像条件为叠前逆时偏移互相关成像条件，其数学表达式为
(4)
式中，S为激励源波场，R为接收波场，t和Txi分别表示时间和炮数。这种成像条件实际是利用发射波场和接收波场的互相关。
2.4 频率域逆时偏移成像
由于嫦娥五号工作任务要求在3个小时内完成钻取采样任务，希望探月雷达数据处理与成像的时间控制在分钟量级。为提高计算效率，课题组提出了一种分层介质的格林函数的频域逆时偏移算法[12]，选取了零时滞互相关通常用作成像条件：
(5)
式中，
为成像域中网格的位置，(x,z)和(x,y,z)分别对应二维和三维情况，
为探月雷达的重建图像，
为源波场，
为接收器波场，它是通过外推时间反转由对应源的产生的接收记录波形，T为记录的探月雷达数据的时间窗口。探月雷达通常有多个源或源位置有多个以便逆时偏移处理，因此最终重建图像是通过叠加所有源的每个逆时偏移图像生成的。
通过计算探月雷达信号的因果关系，公式(5)可以按照文献[13]的卷积形式进行重写：
(6)
由于时域中的卷积与频域中频谱的相乘相对应，方程(6)可以写成：
(7)
式中,
和
分别为源波场和接收波场，和为源波场和接收波场对应的傅里叶变换形式。
在成像区域中的波场的频谱，即
和
是直接通过分层介质格林函数与源频谱以及相应的接收器记录的探月雷达数据的频谱乘积来计算，其表达式如下：
(8)
(9)
式中,
为分层介质格林函数，
为源在位置rs处激发的源波形的频谱，
为接收器在位置
记录的探月雷达数据的频谱。
分层介质格林函数是成像域中每个网格的每个频率的二阶张量。探月雷达通常在垂直宽边模式下工作，并使用线性极化天线进行发送和接收。这样，探月雷达仅记录电场的一个极化分量，则
(10)
算法流程如图7所示，由于上式只有两个参数，只需计算格林函数
的一个分量，如gyy，可以大幅度减少计算分层介质格林函数的成本。
图7 分层介质格林函数的频域逆时偏移成像算法流程
3 仿真试验
3.1 月壤分层结构
为验证速度谱分析法可使嫦娥五号能够准确获取月壤层状介质信息，通过正演GPRMax2.0正演模拟软件建立了如图8(a)所示的仿真模型，采用时域有限差分方法(FDTD)模拟电磁场的传播。该分层模型包含三层月壤介质，相对介电常数从上至下分别为2，2.5和3，电导率均为0.000 01 S/m。天线简化为点源进行模拟，共12个，其空间布局和高度都设置成与嫦娥五号探测器一致。3个方向的网格的尺寸均为7.5 mm，馈电脉冲采用中心频率为2 GHz的Ricker子波。
将正演模拟后得到的探月雷达数据仿真结果进行天线延时矫正、增益、道间均衡等预处理，预处理后得到的结果如图8(b)所示，图中展示了 1-11号天线作为发射天线时，所采集的110道信号，可清楚地分辨月壤表面、第一层月壤界面和第二层月壤界面的雷达反射信号，随着收发天线偏移距的增加，反射信号的双程走势沿双曲线增长。
由该正演数据计算得到的速度谱如图8(c)所示，可以提取月壤表面、第一层月壤界面和第二层月壤界面的雷达反射信号的叠加速度及双程零偏移距走时，根据公式(1)、(2)，可计算得到月表上方空气层、第一层月壤介质和第二层月壤介质的电磁波速度及厚度，如表1所示。与其真实值比较接近，最大误差仅为11.2%。
图8 月壤分层结构的正演模拟及速度反演
表1 月壤分层结构的介电常数及厚度估算结果
参数空气层月壤层1月壤层2实际速度/(m·ns-1)0.30.210.19计算速度/(m·ns-1)0.3080.2290.188误差/%2.79.01.1实际厚度/cm955050计算厚度/cm94.955.652.3误差/%0.111.24.6
3.2 月壤中孤立月岩目标成像
为验证笔者课题组所提出的频域逆时偏移算法准确性，本文使用多输出多输入探月雷达系统进行二维数值仿真实验进行分析，并将计算结果与基于时域有限差分的时域逆时偏移算法进行比较。建立月壤中孤立月岩模型，如图9(a)所示，地下层(风化石)和圆柱形岩石物体的相对介电常数分别设定为2.5和5，对应电导率分别设定为10-5 S/m和10-4 S/m。为与嫦娥五号月球探测雷达系统中使用的天线阵列相类似，本模型同样轮流使用12个天线作为发射天线，剩余的11个天线作为接收天线，从而可得到132道合成探月雷达数据。使用中心频率为2 GHz的雷克子波作为源激励，且将接收天线的接收信号中的直达波信号滤除。
二维成像区域的尺寸为1.6 m×1.8 m。仿真图像区域中的场点的空间采样间隔在x和z方向均为4 mm。对于频域逆时偏移成像，源频谱和接收频谱均以20 MHz的频率步进从20 MHz到4 GHz进行采样。
图9(b)、(c)分别为时域逆时偏移和频域逆时偏移使用半空间初始模型重建的二维图像结果。除源位置处的存在相位差之外，两个重建图像之间没有区别，导致源位置处的相位差的原因是时域有限差分方法不能精确地模拟非常接近源点的区域中的电场。被掩埋的岩石和地表都可以清晰地成像，岩石的顶部和底部界面可被准确识别，岩石的顶部位置大致在深0.3 m处与实际的位置相当吻合。相比时域逆时偏移算法，频域逆时偏移算法消耗的内存不到五十分之一，计算时间成本只有百分之一。由于可以提前计算分层介质格林函数，因此可以在快速地根据嫦娥五号采集的探月雷达数据生成二维地表图像，满足嫦娥五号的任务要求。
图9 时域逆时偏移和频域逆时偏移效果对比
4 地面验证试验
本文通过建立火山坑地面模型试验对提出的频域逆时偏移成像算法进行验证，如图10所示，实验室采用7 m×3 m×2.5 m的火山灰填充坑模拟月壤，火山灰的相对介电常数约为2.5，其介电性能与月壤比较接近，电导率可忽略不计。选取深1 m处的大理石板模拟月岩结构，放置的大理石板的尺寸为60 cm×60 cm×3 cm，相对介电常数大约为8，在灰坑的底部(2.5 m深处)，放置整块金属板。
(a) 嫦娥五号月壤结构探测着陆器的实际尺寸模型
(b) 模型试验结构示意图
图10 火山坑地面模型试验
装备实际尺寸模型的嫦娥五号着陆器放置在火山坑表面以检测地下目标，并对实测数据进行三维逆时偏移成像处理。首先根据估计的表面反射信号到达时间，从接收记录的132道探月雷达数据中去除直达波。频域逆时偏移和时域逆时偏移算法均采用半空间速度模型。由于本文只对比地下区域的成像结果，因此频域逆时偏移算法不考虑计算地面上的分层介质格林函数。相反，由于时域有限差分方法必须从源点外推波场，时域逆时偏移方法须考虑地上部分。地下区域的空间采样步长均为4 mm，共有440×220×340个网格。源频谱和接收频谱采用20 MHz的频率步进从20 MHz至4 GHz采样。
由两种逆时偏移算法获得的三维地下成像结果如图11所示，由于两种算法成像结果的数值仅代表对比度的不同，因此对地下区域的数值进行了归一化处理，通过不同角度剖面所得图像对比可知，两种算法重建的三维图像基本没有区别，且两种算法的成像结果的位置尺寸与实际情况相当吻合。
(a) 时域逆时偏移
(b) 频域逆时偏移
图11 探月雷达实测数据进行三维成像结果的比较
从图11可以看出，其中y方向的成像分辨率比x方向的成像分辨率稍差，这是由于在y方向上部署的只有两个Vivaldi天线所产生的孔径有限。表2分别给出了频域逆时偏移和时域逆时偏移的计算成本。三维时域有限差分模拟计算成本相对较高，其消耗了较大的存储空间和近三天的计算时间。由于分层介质格林函数可以预先计算并保存在磁盘中，因此实际频域逆时偏移算法的三维成像所需的计算时间约为1 h，可以满足嫦娥五号任务的要求。
表2 时域逆时偏移与频域逆时偏移计算成本比较
成像算法CPU核/线程内存/MB时间/s二维频域逆时偏移1/12025二维时域逆时偏移1/1140010800三维频域逆时偏移1/112003570三维时域逆时偏移1/15400227000
5 结束语
本文提出了一种适合嫦娥五号探月雷达的数据处理方法，并实现月壤结构的三维高精度反演成像。通过对嫦娥五号月壤结构探测仪系统进行子波测试、天线相位中心校准和背景干扰去除等校正，有效地去除了嫦娥五号月壤结构雷达探测仪在探测过程的干扰杂波；利用速度谱分析方法获取了月壤分层结构模型的介电常数和厚度，计算值与真实值最大计算误差仅为11.2%，为逆时偏移成像提供了初始模型。地面验证试验结果表明：同等精度的频域逆时偏移算法在对层状结构进行成像时，计算时间仅为时域逆时偏移的1.5%，占用内存仅为其22.2%，计算效率有显著的优越性。综上所述，本文所提的频域逆时偏移算法可满足嫦娥五号搭载的探月雷达系统在分钟量级时间内对着陆器下方的月球表面可能的岩石进行高精度图像的要求，为嫦娥五号钻取采样任务提供了重要的信息支持，并可为未来的月球探测任务提供技术方法支撑。
[1] 丁春雨，封剑青，郑磊，等. 雷达探测技术在探月中的应用[J]. 天文研究与技术，2015，12(2)：228-242.
DING Chunyu, FENG Jianqing, ZHENG Lei, et al. A Review of Applications of Radar-Detection Techniques in Lunar Exploration[J]. Astronomical Research and Technology, 2015, 12(2):228-242.(in Chinese)
[2] 欧阳自远. 我国月球探测的总体科学目标与发展战略[J]. 地球科学进展，2004(3)：351-358.
[3] 郑磊，苏彦，郑永春，等. 地基雷达技术及其在太阳系天体探测中的应用[J]. 天文学进展，2009，27(4)：373-382.
[4] 尚可，钱荣毅. 雷达探测技术在月球科学探测研究中的进展[J]. 地球科学前沿，2017，7(2)：158-166.
[5] XIAO Y，SU Y，DAI S，et al. Ground Experiments of Chang’e-5 Lunar Regolith Penetrating Radar[J]. Advances in Space Research，2019，63(10)：3404-3419.
[6] LU W，LI Y，JI Y，et al. Ultra-Wideband MIMO Array for Penetrating Lunar Regolith Structures on the Chang’e-5 Lander[J]. Electronics，2020，10(1)：1-7.
[7] LI Y，LU W，FANG G，et al. The Imaging Method and Verification Experiment of Chang’e-5 Lunar Regolith Penetrating Array Radar[J]. IEEE Geoscience and Remote Sensing Letters，2018，15(7)：1006-1010.
[8] LIU H，TAKAHASHI K，SATO M. Measurement of Dielectric Permittivity and Thickness of Snow and Ice on a Brackish Lagoon Using GPR[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2014，7(3)：820-827.
[9] LIU H，SATO M. In Situ Measurement of Pavement Thickness and Dielectric Permittivity by GPR Using an Antenna Array[J]. NDT & E International，2014，64：65-71.
[10] LIU H，LONG Z，TIAN B，et al. Two-Dimensional Reverse-Time Migration Applied to GPR with a 3-D-to-2-D Data Conversion[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2017，10(10)：4313-4320.
[11] 朱尉强，黄清华. 探地雷达衰减补偿逆时偏移成像方法[J]. 地球物理学报，2016，59(10)：3909-3916.
ZHU Weiqiang, HUANG Qinghua. Attenuation Compensated Reverse Time Migration Method of Ground Penetrating Radar Signals[J]. Chinese Journal of Geophysics, 2016,59(10):3909-3916.(in Chinese)
[12] LIU H，LONG Z，HAN F，et al. Frequency-Domain Reverse-Time Migration of Ground Penetrating Radar Based on Layered Medium Green’s Functions[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2018，11(8)：2957-2965.
[13] XU K，ZHOU B，MCMECHAN G A. Implementation of Prestack Reverse Time Migration Using Frequency-Domain Extrapolation[J]. Geophysics，2010，75(2)：61-72.
