面向月球科研站资源支持的永久阴影区水冰红外遥感探测方法研究

摘要

水冰是支撑国际月球科研站建设和运行所需的重要资源。红外遥感是探测月球水冰的重要手段，但却难以直接用于无阳光直射的永久阴影区。提出了一种基于光照散射强度模型的阴影区光谱信息校正方法，并构建了基于光谱吸收特征的水冰探测方法，实现了对于永久阴影区水冰信号的探测。应用该方法对月球南极沙克尔顿（Shackleton）和诺比利（Nobili）撞击坑进行分析，发现它们内部的永久阴影区存在水冰信号。该研究有望为国际月球研究站的选址和建设提供借鉴。

关键词

月球水冰; 红外遥感; 永久阴影区; 国际月球科研站

我国于2016年最早提出建立国际月球科研站（International Lunar Research Station）的概念。国际月球科研站将配备能源、通信导航、远程运输、天地往返、地面支持和生命保障等设施，支持探测利用、月球观测、基础科学实验和技术验证等多目标科研活动^［1］。国际月球科研站的建立不仅有助于月球资源开发，还能为人类在外太空扩展提供实践和经验，为在月球、火星甚至更远的行星上建立永久居住点奠定基础。我国对于国际月球科研站的构建规划分为三个阶段^［2］。第一阶段为“探测”阶段，主要包含嫦娥六号和七号任务。嫦娥六号已于2024年发射，并首次从月球背面采集样品。嫦娥七号将于2026年前后发射，开展月球南极的环境与资源详查^［3］。第二阶段为“建造”阶段，主要包括嫦娥八号及后续的五次探测任务^［2］。其中嫦娥八号将于2028年前后发射，开展月球资源利用试验验证，构建月球科研站的基本型。第三阶段为“使用”阶段，开展日地月空间环境探测及科学试验，建成鹊桥通导遥综合星座，形成完善型月球科研站，并逐步升级到实用型、多功能的月球基地^［2］。

月球南极是我国未来几年与月球科研站相关探测任务的目标区域。由于月球的自转轴倾角相对于轨道面的倾斜角度相对较小，月球南极存在无阳光直射的区域，形成永久阴影区（Permanently Shadowed Region）^［4］。永久阴影区内部温度极低，可能会积累大量水冰和其他挥发性物^［5］。月球南极的水资源潜力是其成为我国国际月球科研站重点关注区域的重要原因^［6］。长期以来，遥感一直是探测月球水冰的主要手段，包括中子、雷达和红外等数据类型^［7］。通过热中子计数可较为准确地探测月球南极的氢含量及其分布，但探测分辨率一般较低^［8］。雷达回波极化信号中包含水冰信息，但却会受到月表粗糙度的显著影响，探测结果仍存在较大争议^［9-10］。水冰对红外光谱在多个频段内都具有明确的吸收作用，并形成特有的光谱曲线，这使得红外遥感数据在水冰探测方面具有其独特优势，成为月球水冰遥感探测的重要手段^［7， 11- 12］。然而，永久阴影区内部无阳光直射，传感器接收到的红外光谱信息来自于周边地物反射，难以直接用于水冰探测分析。Li等提出通过除以月球极区有阳光直射区域反射率均值的方法对永久阴影区的光谱信息进行校正^［12］。这种全局平均校正方法虽然简单易行，但未考虑局部地形、光照和月壤物理性质对于散射光的影响^［12］。此外，永久阴影区光谱数据的信噪比较低，噪声影响导致光谱数据在光谱域上剧烈波动^［13］。这些波动会掩盖实际光谱特征，不同区域和时间获取的光谱数据的波动程度也存在差异，增加了水冰信息提取的不确定性。

针对上述永久阴影区水冰红外遥感探测的难点，提出了一种基于光照散射强度模型的阴影区光谱信息校正方法，以实现对于光谱形态的还原；在此基础上，构建了一种基于光谱曲线吸收特征的水冰探测方法，以实现对于永久阴影区内部水冰信号的探测。利用所提出的方法对沙克尔顿（Shackleton）和诺比利（Nobili）撞击坑内部永久阴影区的水冰进行探测，相关结果有望为国际月球科研站的选址和建设提供借鉴。

1 研究区域与数据

1.1　研究区域

现有研究对月球南极永久阴影区的水冰进行了探测与分析，部分研究结果支持沙克尔顿撞击坑内部含有水冰，但在分布和含量等方面仍存在争议^［12， 14-17］。鉴于良好的资源潜力和科学价值，沙克尔顿撞击坑已受到未来月球水冰探测任务的关注。除了沙克尔顿撞击坑之外，诺比利撞击坑也是月球水冰探测的重点区域。美国极地漫游车（VIPER）预计将着陆于诺比利撞击坑附近，并开展水资源探测。诺比利撞击坑分布有大面积的永久阴影区，现有遥感探测结果表明其内部可能含有水冰^［12］。因此，本文选择将沙克尔顿和诺比利撞击坑作为研究对象。

图1  沙克尔顿和诺比利撞击坑周边地形和光照情况

Fig.1  Topography and illumination in the vicinity of Shackleton and Nobili craters

1.2　研究数据

月船一号搭载的月球矿物绘图仪（Moon Mineralogy Mapper， M3）提供了0.43μm~3μm范围内85个波段的光谱辐照度。本文使用OP2c时期覆盖月球南极80º以上区域的所有影像。

月球南极表面数字高程模型（Digital Elevation Model）数据来自于月球轨道激光高度计（Lunar Orbiter Laser Altimeter）生成的栅格计（Gridded Data Records）产品。基于该数据，利用本课题组自主研发的高精度太阳光照模型^［18］，可获得月球南极表面与M3影像拍摄时刻相对应的太阳光照数据。

永久阴影区是指在月球进动周期内无直射太阳光照的区域，即太阳光照率为零的区域。通过计算月球进动周期内太阳直射时间所占的比例可以获得月球南极的太阳光照率。月表最高温度数据来自于由月球勘测轨道飞行器（Lunar Reconnaissance Orbiter）上面搭载的月球红外辐射计（Diviner）所反演的月表温度^［19］。此外，本文还纳入了来自于美国地质调查局（United States Geological Survey）光谱库的实验室纯水冰光谱数据。

2 研究方法

2.1　光谱信息校正

永久阴影区内部无阳光直射，传感器接收到的红外光谱会受到周边反射源地物信息的影响，难以直接用于水冰探测。因此，本文提出了一种基于光照散射强度模型的阴影区光谱信息校正方法。该方法的基本原理如图2所示。对于光照区像元（如示意图中A点），其光谱反射率可表示为卫星载荷记录的光谱辐照度与月表太阳光谱辐照度的比值。对于永久阴影区内部的任一像元（如示意图中P点），其入射光主要来自于周边有直射光照区域（如示意图中的S₁， S₂， …， S_n）对于太阳光的散射。假设卫星传感器记录的光谱辐照度为$I\left(\lambda \right)$，周边散射进入永久阴影区的光谱辐照度之和为${\displaystyle \sum _{i=\mathrm{1}}^n}{I}_i\left(\lambda \right)$，那么校正后的光谱反射率可表示为$R\left(\lambda \right)=I\left(\lambda \right)/{\displaystyle \sum _{i=\mathrm{1}}^n}{I}_i\left(\lambda \right)$。

图2  永久阴影区光谱信息校正方法原理示意图

Fig.2  Schematic diagram of the principle of the spectral correction method for permanently shadowed regions

图3  研究方法流程图

Fig.3  Methodology flowchart

如图3所示，本文所提出的光谱信息校正方法包括周边散射源定位、散射辐照度求解、以及光谱信息校正三个核心步骤。

（1）周边散射源定位：对于永久阴影区内部任一目标像元，基于地形数据利用视域分析方法可以获得周边可视区域。与此同时，基于实时光照模型获取与永久阴影区像元拍摄时刻一致的阳光直射区域。周边可视区域与阳光直射区域的交集即为周边散射源像元的分布区域。

（2）散射辐照度求解：对于永久阴影区像元（假设为P像元）周边所寻找到的任一散射源像元（假设为S像元）。假设其某一波段（$\lambda$）来自太阳直射的光谱辐照度为$I_{\mathrm{s}\mathrm{u}\mathrm{n}}\left(\lambda \right)$，由S像元散射进入P像元的光谱辐照度为$I_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}}\left(\lambda \right)$，那么二者有以下关系：

$$I_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}}\left(\lambda \right)=\tau \times I_{\mathrm{s}\mathrm{u}\mathrm{n}}\left(\lambda \right)\times r(\lambda ,i,e,g)$$

（1）

式中：$\tau$为S像元处太阳圆盘可见率，能够由光照模型获得。$r(\lambda ,i,e,g)$为S像元至P像元的光谱反射率，可由Hapke模型求解：

$$r\left(\lambda ,i,e,g\right)=\frac{\omega (\lambda )}{\mathrm{4}\mathrm{\pi }}\frac{{\mu }_{\mathrm{0}}}{{\mu }_{\mathrm{0}}+\mu }\left\{\left[\mathrm{1}+B\left(g\right)\right]P\left(g\right)+H\left({\mu }_{\mathrm{0}},\omega (\lambda )\right)H\left(\mu ,\omega (\lambda )\right)-\mathrm{1}\right\}$$

（2）

式中：$i$为S像元处太阳入射角，$e$为出射角，$g$为相角；${\mu }_{\mathrm{0}}$为入射角$i$的余弦值，$\mu$为出射角$e$的余弦值；$B(g)$和$P(g)$分别为后向散射函数和相位函数，$H\left({\mu }_{\mathrm{0}},\omega \right)$和$H\left(\mu ,\omega \right)$为描述多重散射过程的函数。本文中$B(g)$，$P(g)$，$H\left({\mu }_{\mathrm{0}},\omega \right)$，$H\left(\mu ,\omega \right)$的求解方法和参数设置与以往研究保持一致^［20］。$\omega (\lambda )$为$\lambda$波段在S像元处的单次散射反照率，可利用该像元所在位置阳光直射状态下M3影像记录的光谱反射率，通过反向求解Hapke模型获得。

（3）光谱信息校正：通过上述方式，可以得到永久阴影区像元P周边所有散射源像元（S₁， S₂， …， S_n）散射的光谱辐照度（$I_{Scatter}^{\mathrm{1}}\left(\lambda \right)$， $I_{Scatter}^{\mathrm{2}}\left(\lambda \right)$，…， $I_{Scatter}^n\left(\lambda \right)$）。那么P像元所在位置校正后的反射率（$R\left(\lambda \right)$）可表示为

$$R\left(\lambda \right)=I(\lambda )/\left({\displaystyle \sum _{k=\mathrm{1}}^n}{I}_{Scatter}^k\left(\lambda \right)\times {\Omega }_k\times cos(i_k)\right)$$

（3）

式中：$I(\lambda )$为卫星传感器记录的P像元所在位置得光谱辐照度；$I_{\mathrm{S}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}}^k\left(\lambda \right)$为P像元周边第k个散射源像元散射进入P像元的光谱辐照度；$i_k$为散射光线在P像元处的入射角；${\Omega }_k$为投影立体角。

通过上述方式可实现对于永久阴影区内部每个像元光谱信息的校正，并获得校正后的光谱反射率。

2.2　水冰探测

永久阴影区内部光谱数据的信噪比很低，表现为在光谱域上的随机抖动。光谱抖动会掩盖实际的水冰吸收特征，现有研究采用曲线拟合的方式从抖动的光谱数据中提取水冰信息^［12］。本文在以往研究的基础上，构建了基于光谱吸收特征的水冰探测方法。如图3所示，该方法主要包括光谱曲线拟合和水冰探测分析两个部分。

（1）光谱曲线拟合：对于永久阴影区内部任一目标像元，首先对其校正后的光谱反射率进行归一化处理，以突出光谱吸收特征。之后，采用csaps三次样条曲线插值方法对归一化后的光谱反射率进行拟合。该方法在曲线拟合时兼顾了光谱曲线整体分布形态和局部变化特征，其基本思路如下：为了实现对于光谱曲线整体形态的最佳拟合，构造与原始光谱值（$\mathit{y}=\left\{y_{\mathrm{1}},y_{\mathrm{2}}...y_n\right\}$）的差异最小的目标样条函数（$f$），即（$\mathrm{m}\mathrm{i}\mathrm{n}({\displaystyle \sum _{j=\mathrm{1}}^n}{\left|y_j-f(x_j)\right|}^{\mathrm{2}})$）。与此同时，利用目标拟合函数的m阶导数平方积分引入粗糙度惩罚因子（$\int {\left[D^{m}f(s)\right]}^{\mathrm{2}}ds$），以实现对于光谱曲线局部变化特征的捕捉。二者之间的权重可根据光谱曲线的抖动程度自适应选取。该方法能够根据光谱数据抖动情况自适应地进行最优拟合，有助于从光谱噪声信号中凸显水冰吸收特征。

（2）水冰探测分析：基于拟合后的光谱曲线，通过光谱连续统分析寻找其发生光谱吸收谱段的中心位置。对于永久阴影区内部任一目标像元，将其所寻找到的光谱吸收中心与已知的水冰吸收谱段进行对比。若在已知的水冰吸收谱段范围内均发现了光谱吸收中心，且吸收特征明显（吸收强度 > 10%），那么认为该永久阴影区像元为潜在的水冰像元。参考以往研究^［12， 21］，本文纳入了三个水冰吸收谱段，它们的中心位置约在1.3 μm， 1.5 μm和2.0 μm附近，具体如表1和图4所示。

表1  本文所使用的三个典型水冰吸收谱段波长范围

Tab. 1  Three typical wavelength ranges of the water ice absorption spectrum used in this study

水冰吸收谱段中心位置/μm	水冰吸收谱段中心波长范围/μm	水冰吸收谱段两侧波长范围/μm
1.3	1.242~1.323	1.130~1.350
1.5	1.503~1.659	1.420~1.740
2.0	1.945~2.056	1.820~2.200

图4  实验室纯水冰光谱反射率曲线和水冰吸收谱段范围

Fig.4  Spectral curves of laboratory pure water ice and the spectral range of water ice absorption

对于所探测的潜在水冰像元，假设其光谱反射率所形成的光谱向量为$\mathit{M}=\left\{m_{\mathrm{1}},m_{\mathrm{2}},...,m_n\right\}$，实验室纯水冰光谱数据所形成的光谱向量为$\mathit{W}=\left\{w_{\mathrm{1}},w_{\mathrm{2}},...,w_n\right\}$，二者的光谱角（spectral angle， S_A）如下：

$$S\underset{}{{}_{\mathrm{A}}}=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left(\frac{\mathit{M}W}{\left|\mathit{M}\right|\left|\mathit{M}\right|}\right)$$

（4）

$S_{\mathrm{A}}$衡量了光谱向量之间的相似程度。$S_{\mathrm{A}}$等于0º表示两个光谱向量完全相似，$S_{\mathrm{A}}$等于90º表示两个光谱向量完全不相似。参考以往探究^［12］，潜在水冰像元的$S_{\mathrm{A}}$应小于30º。

为了进一步增加探测结果的可靠性，本文纳入了光谱信息散度（spectral information divergence， S_ID），其表达式如下：

$$S_{\mathrm{I}\mathrm{D}}={\sum }_{i=\mathrm{1}}^n\left(m_i\mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{m_i}{w_i}\right)+w_i\mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{w_i}{m_i}\right)\right)$$

（5）

$S_{\mathrm{I}\mathrm{D}}$衡量了光谱向量在它们概率分布上的相对熵差异，能够捕捉到光谱信息分布上的细微变化。$S_{\mathrm{I}\mathrm{D}}$ 的数值范围是 ［0，+∞），其中 0 表示两个光谱向量完全相同。本文将水冰像元的$S_{\mathrm{I}\mathrm{D}}$限制在相对较小的范围（$S_{\mathrm{I}\mathrm{D}}$ < 8）。

3 结果与讨论

3.1　散射辐照度

图5展示了月球南极80º以上区域单次散射反照率（single scattering albedo）的空间分布情况。可以看到，单次散射反照率的整体分布呈现出下侧低于上侧、极区低于非极区的趋势。撞击坑周边及其内部区域的单次散射反照率通常也相对较低。由于剔除了阴影区和质量较差的数据，部分区域的单次散射反照率出现了缺失。若某一永久阴影区像元周边所寻找到的散射像元的单次散射反照率出现了缺失，该像元各波段的单次散射反照率由其他像元对应波段的均值替代，以确保散射辐照度的顺利求解。

图5  月球南极单次散射反照率（1-2.5 μm波段均值）

Fig.5  The single scattering albedo of lunar south (average of 1-2.5 μm bands)

图6展示了沙克尔顿和诺比利撞击坑内部永久阴影区像元所对应的周边散射源像元数量。首先，周边散射源像元的数量在空间分布上存在明显差异。沙克尔顿和诺比利撞击坑中心区域的地势低于周边，高程差异可达数公里（图1），这导致位于中心区域的永久阴影区像元难以与周边光照区域形成视线通路。因此，沙克尔顿和诺比利撞击坑内部各永久阴影区像元所寻找到的周边散射源像元数量存在空间差异，总体表现为永久阴影区中心区域低于四周区域的空间分布特点。其次，永久阴影区像元周边散射源像元的数量会随时间发生变化，这与太阳光照角度和范围的变化有关。

图6  永久阴影区内部像元对应的周边散射源像元数量

Fig.6  Number of scattered source pixels corresponding to the permanently shadowed pixels

图7展示了沙克尔顿和诺比利撞击坑内部永久阴影区像元的散射辐照度。散射辐照度的空间分布与周边散射源像元数量存在相似性。周边散射源像元数量多的永久阴影区像元所在位置通常具有较高的散射辐照度。然而，这两者的空间分布并不完全相同。这是由于散射辐照度还受到太阳位置、局部地形和月壤物理性质等因素的影响。此外，通过散射辐照度可以观察到永久阴影区表面的纹理变化，这与以往研究结果具有一定的相符性^［22， 23］。

图7  永久阴影区内部像元的光谱散射辐照度（1-2.5 μm波段均值， mW m^-2 μm^-1）

Fig.7  Spectral scattering irradiance (average of 1-2.5 μm bands, mW m^-2 μm^-1) for permanently shadowed pixels

3.2　水冰分布

图8展示了永久阴影区典型像元在校正前后光谱形态的对比结果。可以看到这些像元的光谱形态在校正后出现了较为明显的改变，部分像元光谱的局部剧烈波动情况得到了改善。更为重要的是，这些像元校正前的光谱形态与实验室纯水冰光谱存在较大差异，但校正后的光谱形态与实验室纯水冰光谱的形态（图4）显示出了更高的相似度。

图8  校正前（左侧）和校正后（右侧）的光谱反射率

Fig.8  Spectral reflectance before (left) and after (right) correction

图9展示了利用本文方法探测到的水冰像元的空间分布。可以看到，沙克尔顿撞击坑内部的水冰像元主要分布在坑壁四周，并且集中分布于撞击坑的左右两侧。诺比利撞击坑内部的水冰像元数量少于沙克尔顿撞击坑，并且倾向于分布于撞击坑的左侧。结合温度数据可以发现绝大多数水冰像元分布于月表最高温度低于110 K的区域。这种低温环境有助于月表水冰的长期保存。此外，水冰探测结果会受到参数设置的影响。如图10所示，两个撞击坑内所探测到的水冰像元数量均会随着吸收强度阈值的增大而快速减少，但都会随着光谱信息散度阈值的增大而快速增加。

图9  永久阴影区内部所探测到的水冰像元（白点）的空间分布

Fig.9  Spatial distribution of water ice pixels (white points) detected in permanently shadowed regions

3.3　不确定性分析

本文所得到的水冰像元的空间分布与以往基于M3影像的探测结果具有一定的相似性（倾向于分布于坑壁四周）^［12］，但也存在差异。本文探测到的水冰像元数量更多且分布范围更广。造成这种差异的可能原因如下：

（1）水冰探测结果会受到探测方法的影响。以往研究虽然对阴影区光谱进行了校正，但方法较为简单（阴影区光谱反射率除以非阴影区光谱发射率的均值）^［12］。本文则为每一个永久阴影区像元单独构建校正模型，考虑了光谱波段、地形、几何形态和物理性质等因素，有助于提升水冰探测结果的可靠性。

（2）水冰探测结果受到数据处理方式的影响。以往研究对M3数据进行了重采样，并将所有影像拼接后统一处理^［12］。然而，考虑到光照情况会随着时间发生变化，本文对每张影像的每个像元逐一进行处理。由于永久阴影区内部M3数据的信噪比较低，光谱数据表现出剧烈波动（图8），这放大了数据处理方式对于光谱形态的影响。

（3）水冰探测结果会受到参数设定的影响。采用更为严格的参数阈值能够增加水冰探测结果的可信度，但会减少探测到的水冰像元数量（图10）。这种参数设定带来的不确定性会导致水冰像元空间分布的变化。

总之，水冰探测结果受到众多因素的影响，导致不同研究结果存在不一致性。考虑到探测结果的不确定性，应更关注水冰的整体分布趋势，而非局部水冰像元的确切位置。以沙克尔顿撞击坑为例，虽然基于不同遥感数据和方法的水冰分布结果存在差异，但多数研究表明其内部存在水冰信号^［12， 14， 16-17］。

图10  参数设置对于永久阴影区水冰探测结果的影响

Fig.10  Influence of parameter settings on the water ice detections in permanently shadowed regions

4 结论

月球水冰有望为国际月球科研站提供水资源，是月球探测任务的重要目标之一。红外遥感是探测月球水冰的重要手段，但却难以直接用于永久阴影区内部。针对该问题，本文提出了一种基于光照散射强度模型的光谱信息校正方法，能够还原永久阴影区内部的光谱特征。在此基础上，构建了一种基于光谱吸收特征的水冰探测方法，实现对于永久阴影区内部水冰信号的探测。应用该方法对月球南极沙克尔顿和诺比利撞击坑进行分析，发现它们内部的永久阴影区存在水冰信号。更为重要的是，所探测到的水冰信号主要分布于坑壁附近，与以往研究结果具有一致性。本文提出的研究方法和分析结果有望为国际月球研究站的选址和建设提供借鉴。

作者贡献声明

杨其全：模型构建，数据处理，论文撰写；

谢欢：学术指导，论文修改，项目资助；

冯永玖：学术指导，结果讨论，论文修改；

柳思聪：结果讨论，论文修改；

王超：结果讨论，论文修改；

刘成保：学术指导，论文修改；

黄潜：数据处理，结果讨论；

黄荣：结果讨论，论文修改；

盛丽艳：结果讨论，论文修改；

彭松：结果讨论，论文修改；

童小华：学术指导，项目构思，论文修改。