1、引　言

了解星体表面矿物类型、丰度和粒径分布是行星探测和科学的基本目标之一，这有助于理解星体地质演化历史和空间环境(Clark 等，1990)。月球表面主要有橄榄石，单斜辉石，斜方辉石，斜长石和钛铁矿等，Lucey(2004) 利用Clementine光谱数据首次对全月表的矿物分布进行了识别，Shuai等(2013)利用嫦娥一号IIM(Interference Imaging Spectrometer)数据对月球表面的矿物进行了定量填图，但是轨道数据的空间分辨率受到限制，更精细的矿物组分、含量和颗粒大小等信息在研究月球的演化历史中有重要的意义。获取月表矿物的丰度能够揭示岩石类型和地球化学性质，分析物质起源和演化过程，估算月球的资源；了解粒径分布能够为模拟月壤提供参考，服务于月球着陆器和巡视器的设计和测试，还可以为研究月表的成熟度/风化程度提供证据，帮助理解月表的演化历史。阿波罗任务从月球采集到了月壤并返回地球，但是所有样品都是在地球上的实验室进行分析的(Pieters 等，2002；Taylor 等，2001a，2001b)，嫦娥三号任务中的玉兔号月球车携带的成像光谱仪对月表土壤进行了实地光谱测量，实地采集到的月壤光谱为理解月球表面的物质成分、含量及粒径分布等提供了一个新的视角。

矿物的丰度和有效颗粒大小可以利用高光谱数据通过求解辐射传输模型求得，主要包括Hapke模型(Hapke，1981)和Shkuratov模型(Shkuratov 等，1999)。例如，在卫星尺度上，Lucey(2004)基于Clementine UV/VIS 光谱数据利用Hapke 辐射传输模型构建查找表获得了斜长石、单斜辉石、斜方辉石、橄榄石在全月表的分布；Yan等(2010)利用Hapke模型将Clementine和端元数据转化成单次散射反照率，基于线性解混模型对斜长石、单斜辉石、钛铁矿、斜方辉石、橄榄石等5类矿物等进行了丰度填图；Shuai等(2013)基于IIM高光谱遥感数据在考虑太空风化的基础上，利用Hapke模型和多端元解混模型得到了全月表典型造岩矿物的丰度分布。Li和Li(2011)利用牛顿迭代法求解Hapke模型得到月壤样品中各端元的丰度和有效颗粒大小，但是这样只能得到一个有效的颗粒大小值，每个像元中各端元的粒径分布更能约束月球演化历史和环境。Zhang等(2015)根据Hapke模型利用在750—1550 nm波段范围去除连续统后的玉兔号测量光谱，在假定矿物颗粒大小已知的情况下(设定矿物粒径为17 μm)进行了矿物的丰度反演；Ling等(2015) 利用修正的高斯模型(MGM)估算了玉兔号月球车行进路线中的矿物丰度。然而，利用高光谱数据同时得到月壤每个端元的丰度和粒径分布的研究十分有限。

本文首先提出联合辐射传输模型和稀疏光谱分解算法的方法来获得每个可观测像元上矿物丰度和粒径分布曲线，然后利用实验室测量数据验证此方法的有效性，最后得到玉兔号月球车着陆区域的矿物类型、丰度及粒径分布曲线。

2、数据及预处理     (2.1) 玉兔号测量数据

嫦娥三号着陆点在雨海盆地(Imbrium)的北部虹湾地区，经纬度为340.49°E，44.12°N(Xiao 等，2015)。玉兔号携带的成像光谱仪在高度约为1 m处测量月壤获得4条光谱(图1)，分为可见光—近红外(450—950 nm)和短波红外(900—2400 nm)两个部分，可见光—近红外谱段的光谱分辨率为2—7 nm，短波红外为3—12 nm(He 等，2014)，本文所用测量光谱被校正到入射角30°出射角0°的反射率(Zhang 等，2015)，并利用小波变换进行了光谱平滑(图2)。为避免月表热效应对反射率的影响，本研究只对0.45—1.8 μm的数据进行分析(Combe 等，2011；Zhang等，2015)。

    (2.2) 实验室测量光谱数据

本文所用实验室测量光谱数据主要包括两个部分：

(1)用于方法验证的橄榄石和古铜辉石混合光谱以及蒙脱石光谱(图3)。这些光谱数据从RELAB(Reflectance Experiment Laboratory)光谱数据库中获得(http://www.planetary.brown.edu/relab/，2017)。由于缺乏已知矿物颗粒分布的造岩矿物数据，因此本研究利用可获得的已知尺寸分布的蒙脱石光谱进行粒径分布反演的验证；橄榄石和古铜辉石的混合物光谱(XT-CMP-31和XT-CMP-32)是在较好的控制条件下测量的，并且经过应用分析表明此样本具有较高的可靠性(Mustard和Pieters，1989)。

(2)用于玉兔号测量数据分析的端元光谱数据。月表的矿物主要有橄榄石、辉石、斜长石和钛铁矿等(Papike 等，1991)，Zhang 等(2015)对玉兔号月球车着陆区域的矿物分析表明辉石、橄榄石和斜长石的存在，Ling等(2015)对玉兔号测量数据分析结果表明还有钛铁矿的存在，Wu和Hapke(2018)认为玉兔号测量光谱特征较弱可能是太空风化影响和熔融玻璃的存在导致的。基于此，本研究利用从RELAB(Reflectance Experiment Laboratory)收集的橄榄石，单斜辉石，斜方辉石，斜长石，钛铁矿和熔融玻璃的光谱作为初始端元(Li和Li，2011；Shuai 等，2013)，这些矿物光谱均为月壤样品的实验室测量光谱(图4)。表1为所采用矿物的编号、粒径和折射率实部(Bell 等，1976；Egan和Hilgeman，1975；Lucey，1998)。

3、方法

本文提出联合辐射传输模型和稀疏光谱分解模型利用玉兔号测量光谱数据反演矿物丰度和粒径分布的方法(图5)。首先利用Hapke辐射传输模型推导出端元的光学常数，然后计算出不同粒径的矿物单次散射反照率，构建解混光谱库，最后利用稀疏分解算法反演矿物丰度及其粒径分布。

矿物在反射率空间的混合具有很强的非线性特征，而各组分的单次散射反照率ω_i可线性相加获得混合物的单次散射反照率ω_mix(Hapke，1981)

${\omega _{\rm{mix}}} = \sum\limits_{i = 1}^K {{F_i}{\omega _i}} $

(1)

    (3.1) 混合物单次散射反照率计算

Hapke等(1981)推导出了反射率与单次散射反照率之间的经验关系式

$r = \frac{\omega }{{4({\mu _0} + \mu)}}\left( {(1 + B(g))P(g) + H({\mu _0})H(\mu){\rm{ - }}1} \right)$

(2)

式中，r是矿物反射率，μ₀和μ分别是入射角和出射角的余弦值，B(g)表示在相位角g处的后向散射，P(g)为在相位角g处的相函数，H(μ₀)和H(μ)为多重散射函数(Hapke，2012)

$H(x) = {\left( {1 -(1 - \gamma)x\left({r_0} +\left(1 - \frac{{{r_0}}}{2} - {r_0}x\right)\frac{{\ln(1 + x)}}{x}\right)} \right)^{ - 1}}$

(3)

式中， $\gamma = \sqrt {1 - \omega } $ ， ${r_0} =(1 - \gamma)/(1 + \gamma)$ 后向散射B(g)可表示为

$B(g) = 1/(1 +(1/h)\tan(g/2)),\;\;h = - \frac{3}{8}\ln(1 - \varphi)$

(4)

式中， $\varphi $ 为填充因子，对于月表风化层，一般设置为0.41(Mustard和Pieters，1989)。

相函数P(g)采用勒让德多项式

$P(g) = 1 + b{P_1}(g) + c{P_2}(g)$

(5)

${P_1}(g) = \cos(g),{P_2}(g) = \frac{3}{2}{\cos ^2}(g) - \frac{1}{2}$

(6)

式中，b和c分别设置为−0.4和0.25(Lucey，1998；Mustard和Pieters，1989)。

    (3.2) 端元光谱库构建

端元矿物的单次散射反照率取决于颗粒大小和光学常数(折射率实部n和虚部k)。折射率实部在可见光—短波红外的变化不大(Roush，2003)，因此本文将n在所有波段视为一个常数(表1)，折射率虚部k可由Hapke辐射传输模型求出(图4)，在k的求解过程中，各矿物多采用的粒径为表1中各端元粒径范围的均值(Li和Milliken，2015)。

因此，给定粒径D，该粒径的矿物单次散射反照率可根据Hapke模型求得

$\omega = {S_e} +(1 - {S_e})\frac{{(1 - {S_i})}}{{1 - {S_i}\varTheta }}\varTheta $

(7)

式中，S_e为表面反射

${S_e} \approx \frac{{{{(n - 1)}^2} + {k^2}}}{{{{(n + 1)}^2} + {k^2}}} + 0.05$

(8)

S_i为内部反射

${S_i} \approx 1 - \frac{4}{{n{{(n + 1)}^2}}}$

(9)

Θ为内部传输系数

$\varTheta = \exp(- \alpha \left\langle D \right\rangle)$

(10)

式中， $\alpha = \dfrac{{4{\text{π}} k}}{{\textit{λ}}}$ ， $\left\langle D \right\rangle = \dfrac{2}{3}\left( {{n^2} - \dfrac{1}{n}{{({n^2} - 1)}^{3/2}}} \right)D$

月表太空风化的影响会增加SMFe(Submicroscopic metallic Fe)和熔融玻璃的含量，降低物质的粒径(Li和Li，2011)。可通过改变加入SMFe的宿主矿物的吸收系数考虑太空风化的影响(Hapke，2001)

$\alpha = \frac{{4{\text{π}} nk}}{{\textit{λ}}} + \frac{{36{\text{π}}{\textit{z}}{M_{{\rm Fe}}}\rho }}{{{\textit{λ}} {\rho _{{\rm Fe}}}}}$

(11)

${\textit{z}} = \frac{{{n^3}{n_{{\rm Fe}}}{k_{{\rm Fe}}}}}{{{{({n_{{\rm Fe}}}^2 - {k_{{\rm Fe}}}^2 + 2{n^2})}^2} + 4{n_{{\rm Fe}}}^2{k_{{\rm Fe}}}^2}}$

(12)

式中，n_Fe、k_Fe、ρ_Fe分别为SMFe的折射率实部、折射率虚部和密度，M_Fe为SMFe的质量分数，ρ为矿物密度。

每个粒径的矿物都被视为是一个端元，由于月球表面的矿物粒径一般小于200 μm(Taylor 等，2001a)，因此本文设定5—200 μm间隔为5 μm的粒径计算每个矿物的单次散射反照率，构成解混端元光谱库。

    (3.3) 稀疏光谱分解模型

稀疏光谱分解模型旨在从一个大的端元光谱库中找出最优的组合来最好的模拟测量光谱(Iordache 等，2011)。由于混合物中的端元个数相对于端元光谱来说是稀疏的，因此在经典线性解混模型上可加入稀疏约束，表达式为

$\mathop {\min }\limits_x {\left\| x \right\|_0}$

(13)

约束条件： ${\left\| {y - Ax} \right\|_2} \leqslant \delta $ ， $x \geqslant 0$ ， ${1^{\rm{T}}}x = 1$

可以将式(13)转化成最优化的形式

$\min \frac{1}{2}\left\| {Mx{\rm{ - }}y} \right\|_2^2 + {\textit{λ}} {\left\| x \right\|_1} + {l_{R_ + ^n}}(x)$

(14)

式中，M为端元光谱库，x是丰度，y为测量光谱，λ为正则化参数, ${l_{R_ + ^n}}(x)$ 为丰度非负项。本文利用变量分裂和增广拉格朗日的稀疏分解算法(SUnSAL)进行求解(Bioucas-Dias和Figueiredo，2010)，可以同时得到端元光谱及其丰度，省去了确定端元光谱的过程。SUnSAL算法(http://www.lx.it.pt/~bioucas/publications.html [2017-12-20])利用ADMM(Alternating Direction Method of Multipliers)方法(Gabay和Mercier，1976)进行式(14)的求解，ADMM可以作为一个变量分割的过程来推导，然后采用增广拉格朗日方法解决由此产生的约束问题。

4、结果分析     (4.1) 实验室测量数据验证

本文首先利用实验室测量光谱数据来验证方法的有效性，包括两个方面：(1)利用已知粒径分布的纯净蒙脱石光谱验证尺度分布反演结果(Robertson 等，2016)；(2)利用已知比例混合的古铜辉石和橄榄石来验证丰度反演结果(Mustard和Pieters，1989)。

(1)粒径分布反演验证。3种不同粒径(25—32 μm，38—45 μm和63—75 μm)的蒙脱石样品用来验证本文方法，图7为反演的粒径分布曲线以及真实的颗粒大小统计。

反演的粒径分布曲线与真实的颗粒大小统计一致。对于25—32 μm蒙脱石样品，本文模拟的平均粒径为31.20 μm，实际统计为31 μm；对于38—45 μm蒙脱石样品，本文模拟的平均粒径为43.84 μm，实际统计为38 μm；对于63—75 μm蒙脱石样品，本文模拟的平均粒径为60.40 μm，实际统计为62 μm。

(2)丰度反演验证。已知丰度比例的颗粒大小在45—75 μm古铜辉石和橄榄石混合物(RELAB ID：XT-CMP-31和XT-CMP-32)用来验证本文丰度反演方法的准确性。图8为模拟的丰度结果及粒径分布曲线。实验室矿物丰度反演精度显示(表2)，各矿物的平均丰度反演误差为5.16%，表明本文方法具有较高的反演精度。由于没有真实的粒径分布数据，无法直接验证粒径分布曲线反演结果，但本文反演出的粒径大部分落在45—75 μm(图8)。

    (4.2) 玉兔号测量光谱数据分析

从玉兔号实地测量光谱数据的反演结果可以看出(图9)，测量光谱和模拟的光谱能够很好的拟合，表3为玉兔号月球车4个观测点的矿物丰度。辉石、橄榄石、斜长石、熔融玻璃和钛铁矿的丰度在4个观测点的平均丰度分别为28.1%、4.5%、39%、28%、0.4%。从各矿物的粒径分布看，在玉兔号着陆区域辉石的平均粒径为166.02 μm，其中橄榄石为8.34 μm，斜长石为196.31 μm，熔融玻璃为44.21 μm，一定程度上表明它们对太空风化的响应存在差异(图10)。

从本研究与其他相关研究(Zhang 等，2015；Ling 等，2015)的对比看(表4)，在丰度值上有较大的差别，主要原因为：(1)本研究将熔融玻璃作为了一个端元，因为在月球表面熔融玻璃的含量相对较高，LSCC(Lunar Soil Characterization Consortium)月海样品的成分分析结果中，熔融玻璃的含量为52.84±7.79 %，最高可达70.27%(Pieters 等，2000)。(2) 本研究用的是反射率光谱，Zhang 等(2015) 在反演矿物含量时采用的是去完连续统后的光谱，因此在一定程度上弱化了反照率差异造成的影响，其只分析了橄榄石、单斜辉石、斜方辉石及斜长石的相对含量；Ling等(2015)对特征位置利用MGM方法进行了半定量的分析。(3)本研究将矿物颗粒大小也作为了自由变量，同时考虑了太空风化的影响。

5、结　论

嫦娥三号任务中玉兔号携带的成像光谱仪获取了非常有价值的月壤实地测量光谱数据，从数据中得到矿物的类型、含量以及粒径分布对于进一步的理解月球演化历史有极其重要的意义。本文提出联合辐射传输模型和稀疏光谱分解算法进行矿物丰度及粒径分布反演的方法，并成功用于了玉兔号测量光谱数据的信息提取。首先利用辐射传输模型计算出玉兔号测量光谱和不同颗粒大小矿物端元的单次散射反照率，然后基于稀疏光谱分解算法反演出各矿物端元的丰度及其粒径分布曲线。

通过对实验室测量光谱的分析表明，本文方法可以很好地反演出混合物中的矿物丰度以及各端元的粒径分布。着陆点是年龄为3.5—2.0 Ga的月海玄武岩单元(Zhang 等，2015)，Apollo采集的月海样本分析(Heiken 等，1991)显示具有大量的玻璃，16.3 % 的斜长石和18.6%的辉石及少量的橄榄石(~3.4 %)。本文玉兔号测量数据结果表明，着陆点附近主要的物质为斜长石、辉石、熔融玻璃和橄榄石，它们在4个观测点的平均丰度分别为39 %、28.1 %、28 %、4.5 %，与Apollo样品分析结果的相对含量比较一致。辉石平均粒径为166.02 μm，橄榄石为8.34 μm，斜长石为196.31 μm，熔融玻璃为44.21 μm，一定程度上表明了各矿物所受太空风化的程度不同(Pieter 等，1993)。另外，本文方法在利用月球轨道高光谱遥感数据进行矿物识别和定量反演研究中具有很大的应用价值。