0 引 言
深空探测雷达是当前探测月球和火星等深空目标的重要地球物理方法［1-2］。其原理是利用电磁波在地下传播的过程中遇到介电性质存在差异的界面时会产生反射、散射、绕射等,通过滤波、时频分析、偏移和反演等多种数据处理与成像方法对地下结构和介质属性进行解释。深空探测雷达的工作频率范围在几兆赫兹到上千兆赫兹之间,对应的分辨率从数百米至分米级别。由于具有采集速度快、设备便携、可穿透性等优势,美国、日本、欧洲以及我国均已研制了多款搭载在轨道器或巡视器上的深空探测雷达系统并成功应用［3-5］。
2003年,欧洲航天局发射的火星快车上搭载了火星地下和电离层探测先进雷达系统（Mars Advanced Radar for Subsurface and Ionosphere Sounding,MARSIS）［5］。MARSIS的工作频率为1.3～5.5 MHz之间的4个带宽为1 MHz的窄带,能够探测数百米至几千米深度的地下结构,但其垂直分辨率较低,约为150 m。2005年,美国宇航局发射的火星勘测轨道飞行器（Mars Reconnaissance Orbiter,MRO）上搭 载了SHARAD（Shallow Radar）［6］。SHARAD的工作频率为15～25 MHz,其垂直分辨率和沿轨水平分辨率分别为约15 m和300 m。Mouginot等提取了MARSIS雷达数据中的表面回波功率,并进行了地形与电离层衰减校正,构建首个火星全球表面介电常数分布图［7］。但由于MARSIS的数据覆盖率较低,其结果缺乏细节信息。与MARSIS雷达数据相比,SHARAD雷达数据的分辨率、信噪比和覆盖率都更高,因此火星轨道雷达数据的分析和解释主要使用SHARAD数据。Castaldo等基于MOLA地形数据,将火星极地地区作为纯净的水冰对SHARAD雷达数据进行了标定和反演,构建了火星全球介电常数图［8］。Lauro等利用MOLA（Mars Orbiter Laser Altimeter）数据筛选平坦区域回波后,反演了火星北极地区Boreales Scopuli区域的基底单元的介电常数［9］。但是目前的火星轨道雷达数据反演仍存在两个问题：1）由于处理过程复杂,雷达数据中有效占比最大的表面反射信息没有得到充分利用；2）尚未有火星样本返回,轨道雷达数据的反演缺乏准确标定。
与轨道雷达探测器不同,安装在火星着陆巡视器上的探地雷达系统可以对巡视路径下方的次表层结构进行高分辨率探测。在火星探测中,美国宇航局的“毅力号”火星车于2021年2月成功着陆在Jezero撞击坑。“毅力号”上搭载的火星地下探测雷达成像仪（Radar Imager for Mars’Subsurface Exploration,RIMFAX）的科学目标是探测火星过去的沉积层、寻找地下水冰以及探测可能保存的有机物质［10］。Casademont等使用双曲线拟合方法反演了“毅力号”火星车获得的RIMFAX雷达数据,计算结果显示Jezero撞击坑浅表层（5 m内）的相对介电常数平均值为8.9［11］。2021年5月,我国“天问一号”火星探测任务的“祝融号”火星车成功着陆在乌托邦平原南部,其搭载的火星次表层探测雷达系统（Mars Rover Penetrating Radar,RoPeR）旨在寻找次表层水冰、远古大洋存在证据及可能保存的生命迹象［4,12］。Liu等针对RoPeR高频数据提出了一种改进的双曲线拟合方法,通过解决反演过程中非唯一解提高了反演精度［13］。然而,由于火星车行驶范围非常有限,轨道雷达数据仍是分析大范围地下结构和介质属性的主要依据。例如,祝融号火星车的巡视距离仅为约1.9 km,其有限的探测范围远不足以揭示乌托邦平原南部的区域地质特征。
为解决上述难题,本文首先对SHARAD雷达数据进行处理,剔除了受电离层、卫星姿态和表面粗糙度影响的数据。之后,处理并反演了RoPeR低频数据用于标定SHARAD雷达数据的表面反射功率。最后,利用处理后的SHARAD雷达数据估算了乌托邦平原南部（20°N～30°N,100°E～120°E）的表面介电常数分布图并进行了地质解释。
1 数据处理与反演方法
本文提出的介电常数反演流程如图1所示,主要包括RoPeR数据处理、双曲线拟合、SHARAD数据处理以及SHARAD数据反演等。
图1 介电常数反演流程图
1.1 RoPeR数据处理
RoPeR低频通道通过一对单极子天线发射和采集雷达数据,具体参数如表1所示。其工作频率为15～95 MHz。本文处理了RoPeR低频通道前1 200 ns的数据用于反演和解释,具体的数据处理步骤如图2所示。
表1 RoPeR低频通道参数［4］
天线类型天线频率天线高度探测深度时窗单极子天线15～95 MHz 903 mm 100 m（冰）15 000 ns
图2 RoPeR低频通道数据处理步骤
1）自检数据道剔除：自检数据道指专门用于检查雷达工作状态的数据,在相应的数据标签文件中被标记为“Self_test”,需要剃除。
2）道间距归一化：早期的RoPeR低频通道数据的道间距为50 cm,但2021年8月17日之后的数据道间距调整为25 cm,本文通过插值处理使道间距统一为25 cm。
3）时基抖动校准：时基抖动是指雷达系统在进行数据采集时受采样时间影响产生的微小随机偏差。本文采用互相关算法将所有数据道的直达波进行对齐,以确保雷达数据解释的准确性。
4）去直流偏移：在使用模数转换器对模拟信号进行数字化时,输出通常包含一个小的直流分量。为消除RoPeR雷达数据中的直流偏移,本文分别对每一道原始数据进行减去该道平均值处理：
式中,A(n)是一道原始雷达数据,A′(n)是去除直流后的雷达数据,n表示单道数据中的第n个样本数,N是该道数据的总采样点数。
5）零时校正：根据数据说明,低频通道原始数据的212.5 ns位置对应实际的数据零时,因此该时间之前数据被删除。
6）背景去除：RoPeR的背景噪声主要来源于天线系统的振铃效应、地表的多次反射以及火星车的反射。本文采用滑动窗口滤波方法消除背景噪声：
式中,A′是原始雷达数据,A″是去背景后的数据,Num为滑动窗口中的数据道数,本文中Num取值为30。
7）带通滤波：采用梯形带通滤波器来提高雷达数据信噪比,该滤波器使用汉明窗来确保通带与阻带之间的平滑过渡,具体参数为：低频边界为15 MHz和50 MHz,高频边界为70 MHz和95 MHz。
8）增益：为了补偿信号随深度衰减的影响,采用包络增益曲线对雷达信号进行处理。包络增益曲线的计算步骤为：首先对雷达数据进行希尔伯特变换计算包络,之后将包络进行叠加和归一化,最后取其倒数作为增益曲线。
1.2 RoPeR数据反演
双曲线拟合方法可以根据雷达图像中的双曲线反射形态反演介电常数,进而解释地下介质的物质成分和地质成因［13］。本文推导了包含空气层折射的双曲线拟合方法,用于RoPeR雷达数据反演（见图3）。根据菲涅尔定律,可以推导出以下方程：
图3 电磁波传播过程的示意图
式中,θ1和θ2分别是电磁波入射角和反射角,ε0和εr分别是真空和地下介质的相对介电常数,且ε0=1。当探地雷达沿着测线移动时,探地雷达天线与埋藏物体之间的单程距离R为
式中,H是天线的高度,h是埋藏物体的深度。天线与埋藏物体之间的水平距离x可以表示为
基于公式（3）、（4）和（5）,天线到地下目标的单程距离R与两者水平距离x之间的双曲线关系为
根据公式（6）不同双曲线的顶点位置和次表层介电常数可以生成一系列双曲线,再通过手动对RoPeR低频通道中的双曲线进行特征点提取来拟合生成的双曲线。两者误差最小时对应的介电常数即为反演的次表层介电常数。
1.3 SHARAD雷达数据处理
SHARAD雷达数据采集高度通常在250～316 km的高空,因此其表面反射功率受电离层效应、卫星姿态（滚转角）和地表粗糙度的影响［14］。为了确保反演结果的可靠性和准确性,本文对SHARAD雷达数据进行以下处理：
1）表面反射功率提取：SHARAD雷达图像中的第一个强反射对应地表反射,其位置可以通过以下公式确定［7］：
式中,S（i）表示在传播时间i下的SHARAD信号强度。在此公式中,当∣S（i）∣达到峰值时,C（i）达到最大值,且此时的时间索引i表示地表回波的位置。
图4（a）为一道SHARAD雷达数据,图4（b）为根据公式（7）计算得到的C值。对比两者发现C在18.9 μs处最大,该位置对应该道数据的表面反射功率为33.14 dB。之后,将该算法用于轨道编号1911202的雷达数据,如图4（c）所示,该轨数据共包含20 964道信号。计算得到的地表反射位置如图4（d）所示,结果表明该算法能够准确定位SHARAD雷达数据中的表面反射位置,进而自动提取对应的表面反射功率。
图4 SHARAD表面反射自动提取示意图
2）电离层影响剔除：电离层的高电子密度导致雷达信号显著衰减。为了避免电离层影响,本文仅使用夜间轨道雷达数据。
3）卫星姿态影响剔除：卫星姿态影响SHARAD发射电磁波在火星表面的入射角,为保证数据的一致性,本文使用翻滚角在90°±1°范围内的雷达数据。
4）地表粗糙度影响剔除：地表粗糙度是决定SHARAD表面反射功率的重要参数,然而当前的火星表面高程数据的分辨率不能定量表征地表粗糙度的影响。本文计算了SHARAD雷达数据的-3 dB脉冲宽度,剔除该宽度大于150 ns的数据［9］,如图5所示。
图5 SHARAD雷达数据的-3 dB宽度计算示意图
1.4 SHARAD雷达数据反演
根据Porcello等（1974年）提出的方法,地表回波功率Ps可以表示为［15］
式中,Pt是传输功率,G是天线增益,λ是雷达波长,H是飞行器高度,Rs是地表反射系数,Ls是地表粗糙度损失项。
根据RoPeR反演结果和SHARAD在该区域的表面反射功率,定义校准系数K：
式中,Ps0、Rs0、Ls0和Hs0分别是SHARAD在祝融号着陆区的表面回波功率、表面反射系数、地表粗糙度损失项和卫星轨道高度。经地表粗糙度影响剔除后,Ls0可以忽略,则校准系数K简化为
在得到校准常数K后,其他位置处的地表反射系数可根据SHARAD表面反射功率Ps和轨道高度H根据下列公式计算得到：
地表反射系数Rs与次表层相对介电常数εr存在关系：
2 结果和解释
2.1 RoPeR结果
图6（a）为RoPeR低频通道的原始数据,图6（b）为经过数据处理之后的数据,处理步骤包括自检数据道剔除、道间距归一化、时基抖动校准、去直流偏移、零时校正、背景去除、带通滤波以及增益。经过处理后的雷达数据信噪比明显提升,从图中可以看出0～100 ns深度范围内的雷达反射比较杂乱,100～500 ns深度范围内包含了多个较连续的倾斜反射。之后,本文采用了1.2节的双曲线拟合算法计算了“祝融号”火星车巡视路径下方地层的介电常数,共选取42道双曲线,前250 m所选取的双曲线位置及介电常数在图6（b）和（c）中列出。
图6 前250 m的RoPeR低频通道雷达图
对反演的介电常数数值进行插值,用相对介电常数4.9进行时程转换后得到如图7（a）所示的介电常数剖面图。从图中发现,“祝融号”火星车巡视路径下方地层的相对介电常数在3～7.8之间,且介电常数值随着深度的增加逐渐增大,对应地下介质的密度也逐渐增大。根据图7（b）,将着陆区地下结构划分为4个主要分层：表层（0～10 m）具有较低的相对介电常数,主要为火星表面风化层或松散的堆积物；第二层（10～30 m）具有适中的相对介电常数（4～5）,可能是古代大洋的沉积物；第三层（30～50 m）显示出较高的相对介电常数（5～6）,可能为破碎的基岩；第四层（50～85 m）具有最高的相对介电常数（6～7.8）,表明其为较密实的基岩。
图7 RoPeR低频通道结果
根据前人研究结果,对于火星北半球低地单元,SHARAD表面反射包含了约10 m深度范围内的次表层参数信息［16］。因此,本文使用图7（a）中0～10 m范围内的平均介电常数（4.32）来计算校准系数K。图8（a）为“祝融号”火星车的巡视路径,该区域内1836701轨道的SHARAD雷达数据距离最近。如图8（a）所示,经过数据处理后,该轨数据中有5个数据点满足计算要求,这些数据点的反射功率在26.83～26.99 dB之间,平均值为26.89 dB。根据公式（10）,校准系数K的值为1.67×1012。
图8 SHARAD数据位于CTX图像中的位置
本文利用RIMFAX数据对计算得到的校准系数K进行验证。轨道6598501的SHARAD数据经处理后有4个数据道距离“毅力号”的行驶区域最近,如图8（b）所示,这些数据道的表面反射功率在55.1～55.9 dB之间,平均反射功率为55.6 dB。根据公式（10）和计算的校准常数K,反演的相对介电常数为8.12～8.26,平均值为8.2。Casademont等反演得到的次表层平均相对介电常数约为8.9［11］。两者之间的绝对误差为-0.7,证明了本文所提出数据处理和反演方法的有效性。
2.2 乌托邦平原南部介电常数图
本文的研究区域为20°N～30°N,100°E～120°E,包含391轨SHARAD雷达数据。图9（a）为原始数据直接提取的表面反射功率,从图中可以看到存在部分轨道反射功率明显偏高或者偏低的现象。图9（b）～（d）分别为经过使用电离层、卫星姿态和地表粗糙度数据剔除处理后的表面反射功率,经过处理后的轨道反射功率变得更加均匀。经过处理后,数据道从2 233 849个减少至87 962个,如图9（d）所示。
图9 SHARAD雷达数据处理过程
使用处理后的地表回波功率（图9（d））反演研究区域的次表层介电常数,结果如图10（a）所示。反演的介电常数值在3.5～4.5范围内,平均值为4.0,与前人研究成果一致［8］。图10（b）为在图10（a）基础上进行插值和平滑获得,图中黑线代表等高线,蓝线代表海岸线［17］,表现出了更多的细节。具体地,研究区域南部的介电常数值较高,而北部区域的值较低。该差异与研究区域的火星地质单元比较符合,其南部为晚西方纪过度单元,而北部为晚西方纪低地单元（图10（c））［18］。晚西方纪过度单元由泥流、河流、湖泊的沉积物和火山岩组成,而晚西方纪低地单元由河流、湖泊、海洋沉积物以及崩积物组成。因此南部地层的组成介质相对于北部地层更粗糙且密度也更高,也就是介电常数的值要更高,与我们的反演结果一致。由于本文反演计算的次表层介电常数分布图仅能反映表面以下约10 m范围内的介质信息,其成因可能与乌托邦平原早期大洋的沉积物或后期的撞击、风化等表面重塑产物有关［19］。
图10 研究区域的介电常数图与地质图
3 结束语
本文提出了一种联合轨道雷达和巡视器雷达的火星表面介电常数反演方法,得到了乌托邦平原南部的表层介电常数分布图。结果表明研究区域的介电常数在3.5～4.5范围内,但是整体呈现南高北低的特征。本文结果中呈现的介电常数分布与乌托邦平原的演化历史相关,未来可用于分析火星古海洋的存在和演变过程。此外,本文提出的方法极大简化了火星轨道雷达表面反射信号的处理和反演过程并得到了验证,该方法可用于更大范围（如乌托邦平原其他区域或克律塞平原等）的表层介电常数反演,为寻找火星次表层水冰、研究火星过去地质历史等提供数据支持。
［1］丁泽刚,曾涛,张光伟,等.分布式地基雷达深空探测技术［J］.雷达科学与技术,2022,20（1）：28-33.
［2］龙腾,丁泽刚,曾涛,等.地基深空探测雷达研究进展与展望［J］.信号处理,2024,40（1）：56-72.
［3］LI Chao,ZHENG Yikang,WANG Xin,et al.Layered Subsurface in Utopia Basin of Mars Revealed by Zhurong Rover Radar［J］.Nature,2022,610:308-312.
［4］ZHOU Bin, SHEN Shaoxiang, LU Wei, et al. The Mars Rover Subsurface Penetrating Radar Onboard China’s Mars 2020 Mission［J］. Earth and Planetary Physics,2020,4（4）:345-354.
［5］JORDAN R, PICARDI G, PLAUT J, et al. The Mars Express MARSIS Sounder Instrument［J］. Planetary and Space Science,2009,57（14）:1975-1986.
［6］SEU R, PHILLIPS R J, BICCARI D, et al. SHARAD Sounding Radar on the Mars Reconnaissance Orbiter［J］.Journal of Geophysical Research:Planets,2007,112（E5）:E05S05.
［7］MOUGINOT J, POMMEROL A, KOFMAN W, et al. The 3-5MHz Global Reflectivity Map of Mars by MARSIS/Mars Express: Implications for the Current Inventory of Subsurface H2O［J］.Icarus,2010（2）:612-625.
［8］CASTALDO L, MEGE D, GURGUREWICZ J, et al. Global Permittivity Mapping of the Martian Surface from SHARAD［J］. Earth and Planetary Science Letters, 2017, 462:55-65.
［9］LAURO S E, MATTEI E, SOLDOVIERI F, et al. Dielectric Constant Estimation of the Uppermost Basal Unit Layer in the Martian Boreales Scopuli Region［J］. Icarus, 2012（1）:458-467.
［10］HAMRAN S-E, PAIGE D A, AMUNDSEN H E F, et al. Radar Imager for Mars’Subsurface Experiment-RIMFAX［J］.Space Science Review,2020（8）:128.
［11］CASADEMONT T M, EIDE S, SHOEMAKER E S, et al.Rimfax Ground Penetrating Radar Reveals Dielectric Permittivity and Rock Density of Shallow Martian Subsurface［J］. Journal of Geophysical Research: Planets,2023（5）:e2022JE007598.
［12］杨佳琦,邵芸,卞小林,等.基于祝融号火星车次表层探测雷达的乌托邦平原浅表层结构分析［J］.中国科学院大学学报（中英文）,2025,42（1）：116-125.
［13］LIU Renrui, XU Yi, CHEN Ruonan, et al. An Improved Hyperbolic Method and Its Application to Property Inversion in Martian Tianwen-1 GPR Data［J］.IEEE Trans on Geoscience and Remote Sensing,2023,61:5104014.
［14］CAMPBELL B A, MORGAN G A, BERNARDINI F, et al. Calibration of Mars Reconnaissance Orbiter Shallow Radar（SHARAD）Data for Subsurface Probing and Surface Reflectivity Studies［J］.Icarus,2021（5）:114358.
［15］PORCELLO L J, JORDAN R L, ZELENKA J S, et al.The Apollo Lunar Sounder Radar System［J］. Proceedings of the IEEE,1974,62（6）:769-783.
［16］CROCI R,SEU R,FLAMINI E,et al.The Shallow Radar（SHARAD）Onboard the NASA MRO Mission［J］. Proceedings of the IEEE,2011,99（5）:794-807.
［17］CLIFFORD S M, PARKER T J. The Evolution of the Martian Hydrosphere: Implications for the Fate of a Primordial Ocean and the Current State of the Northern Plains［J］.Icarus,2001（1）:40-79.
［18］TANAKA K L, SKINNER J A, DOHM J M, et al. Geologic Map of Mars［M］.［S.l.］：Astrogeology Research Program（USGS）,2014.
［19］LI Jianhui, LIU Hai, MENG Xu, et al. Ancient Ocean Coastal Deposits Imaged on Mars［J］. Proceedings of the National Academy of Sciences of the United States of America,2025,122（9）:e2422213122.