利用shearlet变换对探月雷达数据进行月壤结构的重建

张领, 刘财, 胡斌* , 张泽强, 徐晨阳,曾昭发, 王德利, 徐懿

1 吉林大学 地球探测科学与技术学院， 长春 130026 2 澳门科技大学 月球与行星科学国家重点实验室， 澳门 999078 3 中水东北勘测设计研究院有限责任公司， 长春 130061

0 引言

作为中国月球和深空探测的重要组成部分，嫦娥三号于2013年12月着陆月球雨海地区，对着陆区域进行了深入的探测(Xiao et al., 2015).嫦娥三号所携带的探月雷达通过发射高频电磁波脉冲对地下结构进行探测.为了实现探测目标，探月雷达设计了两个通道的天线：第一个通道(CH-1)用于探测浅月壳结构(中心频率为60 MHz，带宽为40～80 MHz)，第二个通道(CH-2)用于探测月壤的厚度与结构(中心频率是500 MHz，带宽为250～750 MHz)(Fang et al., 2014).

国家天文台最初发表了有关探月雷达的原始与处理结果(Su et al., 2014).Xiao等(2015)根据探月雷达的CH-1和CH-2数据对嫦娥三号着陆区域进行了地下地质结构解释，获得了深度360 m以内的9层地下结构.Zhang等(2015)基于探月雷达数据对玄武岩浆的几次覆盖进行了解释.Yuan等(2017)利用雷达数据对月球的平静期进行了研究.值得注意的是，雷达结果解释的前提是数据的正确性，然而，Li等(2018)仔细分析了CH-1数据并与在地面上采集的实验数据进行对比，认为CH-1数据中深部的反射同相轴可能并非是来自地下的信息而是由仪器本身所引起的.

根据CH-2的雷达数据，嫦娥三号着陆区基岩之上的月壤层主要分为三层：最上层是厚度小于1 m的风化层，中间是厚度小于5 m的溅射物层，最下面是厚度约为5 m的古月壤层(Fa et al., 2015).Zhang等(2019)和Lai等(2016)推论出着陆点上包括月壤与古月壤的多层结构及其接受风化与撞击坑溅射覆盖的演化过程.Zhang等(2018a)和Dong等(2017)根据CH-2数据计算了月壤的电性参数与铁钛含量.Hu等(2019)基于CH-2数据利用相似性约束来定位了月壤中的岩石块.Ding等(2020)利用CH-2雷达数据与地形约束了着陆区的介电特性，表明紫薇撞击坑溅射物形成的小山丘是由爱拉托逊纪月海玄武岩构成.

探月雷达第二通道拥有两个接收器天线，因此，可以得到两组CH-2数据(CH-2A和CH-2B).因为CH-2A数据被一种横向噪声所干扰，所以前人的文章都只利用了CH-2B数据.CH-2A数据中这种横向噪声表现为中深部的大量近水平多次波.这种噪声可能是由于接收天线本身以及月球上复杂的采集环境所引起的.

实际上，不仅CH-2A数据被这种水平横向噪声所干扰，CH-2B数据中也存在少量这类噪声.这类水平横向噪声严重影响了CH-2数据的深度利用与解释.为了解决水平噪声问题，我们引入了多尺度几何分析方法中的shearlet变换对噪声进行分析和抑制.与常规傅里叶分析和小波分析相比，多尺度几何分析(例如contourlet变换(Do and Vetterli, 2003)，curvelet变换(Candès et al., 2006)，shearlet变换(Easley et al., 2008)等)可以更好地处理高维信号.它从局部性、方向性和多尺度分析信号，以找到具有表面奇异性或线奇异性的高维信号的最佳或最稀疏表示.近年来，许多多尺度几何分析方法已应用于噪声抑制中(Zhang et al., 2010).这些多尺度和多方向分析方法利用噪声和信号在方向和尺度上的分布差异来抑制噪声(Liu et al., 2014).在这些方法中，shearlet变换具有更好的方向灵敏度和稀疏表示(Häuser and Steidl, 2014).

本文将shearlet变换应用于嫦娥三号探月雷达CH-2A和CH-2B数据的噪声分析和抑制中，并融合两者的结果对嫦娥三号着陆区的月壤结构进行重建.首先，设计了一套针对第二通道雷达数据的预处理流程.随后，利用月壤综合模型验证shearlet变换抑制水平横向噪声的效果.其次，利用shearlet变换对CH-2A和CH-2B数据进行处理，并将两套处理后的图像进行融合.最后，根据探月雷达处理结果对嫦娥三号着陆点的月壤结构进行重建，重建结构中的溅射物厚度与撞击坑溅射物经验公式结果得到了相互印证.

1 数据与方法原理

1.1 探月雷达数据预处理流程

玉兔号月球车携带着探月雷达对地下地质结构进行探测，总共在月表上行走了114.8 m的距离(如图1所示).图1中的底图由嫦娥三号着陆器搭载的降落相机所拍摄.在对探月雷达数据进行进一步的分析与解释之前，需要对雷达数据进行预处理.

图1 玉兔车行进路径图Fig.1 The Yutu′s path on the Moon

根据数据的质量、采集环境以及采集参数设置，本文设计了针对嫦娥三号探月雷达数据的预处理流程.经过预处理后，得到了两套高分辨率的第二通道数据(图2).

图2 探月雷达第二通道数据预处理流程及结果Fig.2 The pipeline chart of CH-2 LPR data processing. N103—N209 denote the positions

N103—N209表示航点对应位置.尽管针对CH-2数据进行了一系列的常规方法处理(例如，带通滤波与去背景等)，CH-2A与CH-2B仍然被一种横向噪声所干扰(图2中黄框部分).这种明显的近水平的横向多次波严重影响了中深部数据的解释.这种横向噪声很可能来源于接受天线仪器本身和周围采集环境的共同干扰.为了解决这些横向噪声的干扰，本文引入了多尺度几何分析方法中的shearlet变换来分析与抑制噪声.

1.2 Shearlet变换

Shearlet变换作为一种有效的多尺度与多方向的方法，可以基于多分辨率分析准确地生成shearlet系数(Yang et al., 2014).与使用仿射系统理论的方法相比，shearlet变换具有更好的稀疏性和更简单的数学结构(Guo et al., 2004, 2006; Heil, 2006).并且，即使在各向异性条件下，shearlet系统也能够区分奇异分量并表征多维信号(Häuser et al., 2014).shearlet变换的连续仿射系统可以描述如下：

{TtDMψ,M∈G,t∈Rn}，

(1)

(2)

膨胀矩阵Mas是两个操作的结合Mas=SsAa，其中S代表剪切操作，A表示尺度操作：

(3)

其中s和a分别代表剪切因子和尺度因子.

此外，对于ξ=(ξ1,ξ2)∈R2,ξ1≠0，ψ定义为

(4)

(5)

t∈R2}.

(6)

考虑到公式(4)满足连续shearlet变换的要求，定义f的连续shearlet变换如下：

SH(f)(a,s,m)∶=〈f·φa,s,m〉

(7)

图3表示出了频域中的shearlet基函数的划分.当在尺度a，角度s和位置t上基函数有所不同时，Aa将会把数据分解为具有不同频带的若干子集.大尺度shearlet结果描述数据细节信息，小尺度shearlet结果表示趋势.同理，Ss按方向划分数据，s用来控制倾角范围.

图3 (a) 不同参数的水平shearlet基函数分频示意图； (b) 不同参数的垂直shearlet基函数分频示意图Fig.3 The frequency division of the horizontal shearlets (a) and vertical shearlets (b) for different parameters

2 模型验证

为了验证shearlet变换去除横向噪声的有效性，通过时间域有限差分方法(FDTD)模拟探月雷达合成数据(Zhang et al., 2018b).月壤模型和雷达模拟结果如图4所示，雷达模拟的有关参数如表1所示.将背景噪声和低切滤波的高斯随机噪声添加到探月雷达数据中(图5)，形成含噪的合成数据.本文的目标是衰减背景噪声(低频，水平方向)和随机噪声(中高频)，因此设计了包括直流分量消除、去背景和带通滤波在内的常规处理方法组合，与shearlet变换去噪效果进行对比.

图4 (a)月壤模型； (b) 雷达模拟结果Fig.4 (a) The lunar regolith model； (b) The simulation result.

图5 (a) 噪声； (b) 加噪数据Fig.5 (a) Noise； (b) Noisy data

表1 正演模拟参数Table 1 Parameter of simulation

图6展示了常规处理手段方法和shearlet变换方法的对比结果图.整体上，两者对大部分噪声都有较好的抑制作用.但是在常规处理结果中仍存在一些残留噪声，这是因为有用信号和噪声信号之间存在频带的重叠.而shearlet变换能将数据映射到几个子域中，这些子域可以稀疏地表示数据，使得噪声和信号更好的分离，从而得到更好的处理效果.图7通过频谱可以看出，shearlet变换能够更加真实地还原原始信号.

图6 (a) 常规处理结果； (b) shearlet变换处理结果Fig.6 Denoised data, (a) the conventional method and (b) the shearlet method

图7 频谱分析结果(a) 原始数据； (b) 含噪数据； (c) 常规处理方法； (d) shearlet变换.Fig.7 Frequency spectrum of (a) original data, (b) noise, (c) conventional denoised data and (d) shearlet denoised data

本文所提出的shearlet变换方法不仅是一种抑制噪声信号的方法，还是信号分析的有利工具.准确的数据分析，有针对性的策略和稀疏的表示方式，使得shearlet变换带来更理想的去噪结果.

3 探月雷达数据分析、处理及融合

3.1 利用shearlet变换进行多尺度几何分析

为更好地说明shearlet变换分析的过程，本节以CH-2A数据中N203到N204航点的数据为例(图10a).嫦娥三号探月雷达数据，尤其是CH-2A数据，存在明显的横向水平噪声(80～150 ns)，这严重影响深部信息的获取.为了更准确地研究月壤深部结构信息，使用shearlet变换(a=4,s=4)来分析探月雷达数据，并研究噪声和信号在频率和方向上的分布情况.Shearlet变换按频率和方向将信号映射到不同的子集中.将相同比例或角度的子集相加，可获得具有不同频带和不同方向带的子信号(见图8和图9).

图8 CH-2A中N203—N204段数据不同尺度下的shearlet分量Fig.8 The shearlet coefficients of the selected portion (N203 to N204 in CH-2A data) in different scales

图9 CH-2A中N203—N204段数据不同角度下的shearlet分量Fig.9 The shearlet coefficients of the selected portion (N203 to N204 in CH-2A data) in different angles

图8显示了雷达数据三个主要尺度(频率)分量.可以看到信号的频带低于噪声，信号集中在中低频段，而噪声集中在中高频段，因此可以在低尺度上观察到月壤的结构.图9显示了雷达数据三个主要角度分量.水平噪声主要集中在水平角度分量中.总之，对于探月雷达的第二通道数据，抑制噪声和提高降结构识别能力的关键是消除数据信号中的水平分量.

3.2 利用改进的水平滑动滤波在shearlet域中去除噪声

从图9b中可以看到，在水平的分量中仍然有许多有用的信号(0～60 ns).本文认为水平shearlet分量中的水平横向噪声是背景场.因此，提出了一种基于shearlet变换的水平噪声抑制方法.该方法实际上是在传统的背景消除方法，即水平滑动滤波器的基础上，代替直接去除水平方向上的小波系数，使用水平滑动滤波来尽可能多地保留有用信号(曾昭发等，2010).水平滑动滤波的目的是去除水平方向上的慢速振荡分量，以突出显示快速振荡的分量.该方法的原理是确定窗口的水平宽度，对窗口中的所有道集求平均，并从窗口中间减去平均值，如公式(8)所示：

(8)

其中，x′(m,n)是处理后的结果，xc是窗口中间的数据，x(m,n)是原始数据，m和n是坐标，N是道集数，用于控制处理数据的平滑度.

然而，时空域中的水平滑动滤波仍然会破坏浅部的有用信号．因在shearlet域中，子集的大小与原始数据的大小相同，为了抑制shearlet域中水平子集分量的背景，shearlet域中的水平滑动滤波可以表示为

(9)

其中，c′s(m,n)是处理的shearlet系数，csc是窗口中间的shearlet系数，cs(m,n)是原始的shearlet系数.

图10显示了处理后的结果.显然，大多数的水平噪声都得到了有效抑制，有效信号(90～130 ns)得以恢复.根据处理结果，可以勾画出一个新的分界面(图10c中箭头所示).所提出的方法会稍微破坏浅部的有用信号(0～60 ns).去噪后的数据为研究月壤的结构提供了良好的图像.

图10 N203—N204段数据的噪声抑制(a) 原始数据； (b) shearlet变换分离出的噪声； (c) 去噪结果.Fig.10 Demonstration of horizontal noise attenuation for the selected portion (N203 to N204 in CH-2A data)

3.3 CH-2A与CH-2B处理结果融合

将shearlet变换去噪方法应用至第二通道的全部数据，图11a和图11b分别显示的是CH-2A与CH-2B数据shearlet变换去噪后的结果.经过处理后，CH-2A大部分的横向近水平噪声都得到了有效抑制(红框).特别是在10～30 m区域，有用信号(受噪声干扰)在一定程度上得以恢复.虽然CH-2A数据的处理结果中存在一些噪声残留，这是因为为了保留更多有用信号，只删除了最水平的shearlet系数.与CH-2A数据不同，CH-2B数据中的水平噪声虽然幅度较小，但是它仍然会影响深部有用信号.CH-2B数据经过处理后，深部数据信号更清晰(红框)，这对于确定月壤的下边界具有一定意义.

图11 去噪后的CH-2数据(a) CH-2A； (b) CH-2B.Fig.11 Horizontal noise attenuation for LPR CH-2 data

尽管CH-2A数据受到噪声困扰，但经过处理后，嘈杂的CH-2A数据信噪比有了显著的提高，其中浅部的信号(30～50 ns)比CH-2B数据清晰.这些浅部信号对于后面将要解释到的撞击坑溅射物的厚度至关重要，因此为了全面利用CH-2A和CH-2B两个数据集的整体优势，本文借鉴了图像融合的思想.具体步骤如下(图12)：

(1)使用shearlet变换将探月雷达CH-2A和CH-2B数据映射到不同的稀疏域.

(2)将稀疏域系数分为低频尺度和高频尺度.

(3)对于低频尺度的系数，选择两个数据的平均值作为融合系数.

(4)对于高频尺度的系数，选择两个数据中的最大值作为融合系数.

(5)将融合系数映射回时空域，获得最终融合结果.

图12 CH-2A和CH-2B数据图像融合流程Fig.12 Image of LPR CH-2A and B data fusion flow chart

经过图像融合以后，获得一个深部信息与浅部信息更加均衡的结果(图13).

图13 CH-2A和CH-2B融合结果，时深转换所用的相对介电常数为3Fig.13 The stacked result of LPR CH-2A and B data. The relative dielectric constant is set as 3

4 结果解释

根据月壤的形成机制、嫦娥三号着陆区域特点以及前人对探月雷达的数据分析，对融合后的CH-2数据进行分层解释(图14)：

(1)基岩.基底岩石是在雨海纪晚期或埃拉托逊纪早期所形成的玄武岩覆盖层(McKay et al., 1991).

(2)古月壤层.在埃拉托逊纪和哥白尼纪，各类陨石撞击月球表面的玄武岩盖层形成古月壤层(McKay et al., 1991；欧阳自远，2005).图中的蓝点代表未彻底风化的岩石碎块或陨石碎块(Hu et al., 2019；Fa, 2020).

(3)溅射层.在哥白尼纪的某个时间点，嫦娥三号着陆点附近的撞击坑形成(图14)，其所形成的溅射物覆盖了其周围的月壤层(Xiao et al., 2015).图中的黄点代表溅射出的岩石碎块(Hu et al., 2019；Fa, 2020).

图14 嫦娥三号着陆区域月壤层结构解释重建图Fig.14 The interpretation of regolith structure on the CE-3 landing site

另一方面，嫦娥三号降落在一个直径约为450 m的撞击坑的附近，距撞击坑的边缘约50 m处.通过众多撞击坑及其溅射物厚度的研究及公式(Mcgetchin et al., 1973)可以计算出着陆点的溅射层的厚度为4 m左右，与上述的分层结果(图14)一致，也侧面验证了融合数据的有效性与分层结构的可信度.

5 结论

本文利用shearlet变换对探月雷达数据进行了噪声分析与去除，并对嫦娥三号着陆区的月壤结构进行重建.首先，在shearlet域中，观察横向噪声的分布，并分离出以噪声为主的shearlet分量，得到重建后的雷达数据.随后，为充分利用嫦娥三号第二通道两套数据的整体优势，对两套经过shearlet变换去噪后的雷达数据进行融合，得到综合探月雷达图像，有效地去除了水平噪声并增强了来自浅月表的有效信号.最后，根据探月雷达处理结果对嫦娥三号着陆点的月壤结构进行重建，重建结构中的溅射物厚度与撞击坑溅射物经验公式结果得到了相互印证.

在分析了CH-2A和CH-2B的融合数据之后，获得了嫦娥三号着陆区域的月壤结果与演化模型.其月壤层可大致分为两层：覆盖在基质玄武岩之上的古月壤层和由附近撞击坑溅射物形成的溅射层.通过溅射物厚度经验公式，验证了融合数据的有效性与分层结构的可信度.

嫦娥四号探月雷达到目前为止也取得了阶段性成果，包括利用第一通道对深部结构进行解释(Zhang et al., 2020; Lai et al., 2020)，以及利用第二通道对嫦娥四号着陆区域溅射层与月壤层的解释(Li et al., 2020; Lai et al., 2019).其中嫦娥四号利用的第二通道也均是CH-2B数据.本文方法为嫦娥四号第二通道的两套数据(CH-2A和CH-2B)的综合利用提供了依据，同时为后续的嫦娥五号以及天问一号的雷达数据处理提供了处理方法支撑.