CN104964931B - 一种火星表面含水矿物定量反演方法 - Google Patents

Abstract

本发明公开一种火星表面含水矿物定量反演方法，该方法包括：S0、获取火星高光谱图像；S1：计算所述火星高光谱图像各波段的信噪比，并选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演；S2：根据预设的光谱特征参数，建立含水矿物识别模型，得到含水矿物识别区；S3：根据S1中所述信噪比大于预设阈值的波段以及预设矿物光谱库中的反射率光谱，构建稀疏解混光谱库；S4：根据所述稀疏解混光谱库，对所述含水矿物识别区进行稀疏解混，得到火星表面含水矿物的丰度。本发明的方法可以克服由于火星表面含水矿物含量低、分布零散、背景矿物不确定的特点导致的传统混合像元分解算法很难对含水矿物进行有效的丰度反演的问题。

Description

一种火星表面含水矿物定量反演方法

技术领域

本发明涉及深空探测矿物定量化反演技术领域，具体涉及一种火星表面含水矿物定量反演方法。

背景技术

火星表面含水矿物的定量化可以对火星表面的矿物成因进行限定，从而对火星的地质演化分析提供支持，并有助于进一步揭示火星早期水环境的细节特征。含水矿物是指含水分子或H⁺、OH^-、H₃O⁺等离子的矿物，火星含水矿物的主要存在形式是层状硅酸盐(主要为Fe/Mg层状硅酸盐和Al层状硅酸盐)和含水硫酸盐等。火星表面含水矿物具有丰度含量低、分布地域零散、背景矿物(像元内除含水矿物以外的其他矿物)不确定或未知等特点。

为了探测火星地质岩矿分布信息，多个火星探测卫星搭载有高光谱探测器，主要包括：1996年MGS搭载的热辐射光谱仪(Thermal Emission Spectrometer，TES)、2003年MEX搭载的可见光及红外矿物制图光谱仪(Visible and Infrared Mineralogical MappingSpectrometer/Observatoire pour la Minéralogie,l'Eau,les Glaces et l'Activité，OMEGA)和2005年MRO搭载的紧凑型侦查成像光谱仪(Compact Reconnaissance ImagingSpectrometer for Mars，CRISM)。

2012年8月，“好奇号”火星车在Gale撞击坑中央峰Sharp山的西北侧，对Gale地区的水环境、地质、气候等进行了实地调查，卫星高光谱遥感图像的含水矿物反演和分析结果可以为“好奇号”行进路线规划和样品采集分析提供支持。

矿物定量反演方法主要有基于诊断吸收谱带的丰度反演、光谱解混和数理统计方法。谱带深度易受环境因素和光谱重建精度等方面的影响，因而具有一定的不确定性；数理统计方法最常用方法是回归分析，需要测量和分析大量的样品，这在火星表面矿物丰度反演中不易实施；混合像元光谱解混算法是目前矿物丰度反演最常用的方法。

传统的光谱解混算法一般分为端元提取和丰度反演两个步骤，但是由于火星表面含水矿物具有含量低、分布零散、背景矿物不确定或未知的特点，使得从火星高光谱图像上提取端元变得相对困难，从而导致传统的混合像元分解算法很难对含水矿物进行有效的丰度反演。

发明内容

本发明所要解决的技术问题是由于火星表面含水矿物具有含量低、分布零散、背景矿物不确定或未知的特点，使得从火星高光谱图像上提取端元变得相对困难，从而导致传统的混合像元分解算法很难对含水矿物进行有效的丰度反演的问题。

为此，本发明提出一种火星表面含水矿物定量反演方法，提高定量反演的精度和效率，该方法包括：

S0、获取火星高光谱图像；

S1：计算所述火星高光谱图像各波段的信噪比，并选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演；

S2：根据预设的光谱特征参数，建立含水矿物识别模型，得到含水矿物识别区；

S3：根据S1中所述信噪比大于预设阈值的波段以及预设矿物光谱库中的反射率光谱，构建稀疏解混光谱库；

S4：根据所述稀疏解混光谱库，对所述含水矿物识别区进行稀疏解混，得到火星表面含水矿物的丰度。

可选的，所述步骤S1，包括：

S11：将所述火星高光谱图像分割为多个预设大小的子区域，并计算各子区域中各波段的噪声强度；

S12：根据所述各子区域中各波段的噪声强度，计算所述火星高光谱图像各波段噪声强度；

S13：计算所述火星高光谱图像各波段的信噪比；

S14：选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演。

可选的，所述步骤S11，包括：

将所述火星高光谱图像分割为多个w×h大小的子区域，计算各子区域中各波段每个像元的残差，公式如下：

其中，r_i,j,k为子区域中第i行j列第k波段的像元残差，x_i,j,k为子区域中第i行j列第k波段的像元值，其中，1≤k≤N，N为所述火星高光谱图像包含的波段个数，当1<k<N时，由x_i,j,k-1和x_i,j,k+1线性拟合得到；当k＝1时，由x_i,j,2和x_i,j,3线性拟合得到；当k＝N时，由x_i,j,N-1和x_i,j,N-2线性拟合得到；其中，w和h为预设正整数；

计算所述各子区域中各波段所有像元的残差的标准差，得到各子区域中各波段的噪声强度。

可选的，所述步骤S12，包括：

根据所述各子区域中各波段的噪声强度，确定各波段所有子区域的噪声强度中的最大值和最小值；

在各波段的噪声强度中的最大值和最小值之间，划分各波段的噪声强度为m个噪声区间，m为预设值；

对于每个波段，将包含子区域数量最多的噪声区间中噪声强度的均值作为该波段的噪声强度。

可选的，在所述步骤S2中，所述预设的光谱特征参数包括：在波长1900nm处水的吸收特征BD₁₉₀₀、在波长2100nm处水的吸收特征BD₂₁₀₀、在2300nm波长处反射率的下降特征D₂₃₀₀以及在波长2290nm处的波峰特征SINDEX；

所述BD₁₉₀₀＝1－((R₁₉₃₀)/(a×R₁₈₅₀+b×R₂₀₄₆))，其中，a＝1－b，b＝(1930－1850)/(2046－1850)；

所述BD₂₁₀₀＝1－((R₂₁₃₂)/(a×R₁₉₃₀+b×R₂₂₅₀))，其中，a＝1－b，b＝(2132－1930)/(2250－1930)；

所述D₂₃₀₀＝1－(R₂₂₉₀+R₂₃₃₀+R₂₃₃₀)/(R₂₁₄₀+R₂₁₇₀+R₂₂₁₀)；

所述SINDEX＝1－((a×R₂₁₂₀+b×R₂₄₀₀)/(R₂₂₉₀))，其中，a＝1－b，b＝(2290-2120)/(2400-2120)；

式中下标表示波长，R表示反射率。

可选的，在所述步骤S2中，所述含水矿物识别模型为：

(BD₁₉₀₀＞t1)or(BD₂₁₀₀＞t2)or(D₂₃₀₀＞t3)or(SINDEX＞t4)

其中，t1、t2、t3和t4为预设常数。

可选的，所述步骤S3，包括：

S31：根据S1中所述信噪比大于预设阈值的波段的光谱范围和波长位置，对预设矿物光谱库中的反射率光谱进行重采样，以使所述反射率光谱与所述信噪比大于预设阈值的波段的光谱范围和波长位置一致；

S32：去除经过S31处理的所述预设矿物光谱库中任意两条光谱角小于预设最小光谱角的光谱的其中一条，得到稀疏解混光谱库。

可选的，在所述步骤S31中，所述重采样的公式如下：

R′_c＝a×R_l+b×R_r

其中，R′_c为重采样后得到的波长c处的反射率，波长c为S1中所述信噪比大于预设阈值的波段的波长位置，R_l和R_r分别为所述预设矿物光谱库中在波长c左右两侧距离波长c最近的波长，即l和r处的反射率值，a＝(r-c)/(r-l)，b＝1-a。

可选的，在所述步骤S4中，所述进行稀疏解混的公式如下：

其中，y为所述火星高光谱图像的反射率，A为所述稀疏解混光谱库，λ表示预设的正则化参数，m为光谱库中的光谱个数，s_j为光谱库中与第j个光谱对应的端元丰度，采用SUnSAL算法进行矿物丰度的求解。

相比于现有技术，本发明的火星表面含水矿物定量反演方法针对火星含水矿物具有丰度含量低、分布地域零散、背景矿物不确定或未知等特点，采用目标识别与稀疏解混相结合的方法进行含水矿物的定量反演，目标识别有助于缩小含水矿物的潜在分布范围，提高丰度反演的运算效率，增加含水矿物丰度反演的准确性；稀疏解混算法不需要对火星矿物类别和端元数量进行专家预判断，而是基于完备光谱库逐像元自动进行矿物类别识别和丰度反演。能够进行有效的火星表面含水矿物定量丰度反演。该方法可推广至深空探测各星体矿物定量反演。

附图说明

图1为本发明实施例提供的一种火星表面含水矿物定量反演方法流程图；

图2为本发明实施例提供的CRISM高光谱图像的信噪比示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本实施例公开一种火星表面含水矿物定量反演方法，该方法可包括以下步骤S0至S4：

S0、获取火星高光谱图像；

S1：计算所述火星高光谱图像各波段的信噪比，并选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演；

S2：根据预设的光谱特征参数，建立含水矿物识别模型，得到含水矿物识别区；

S3：根据S1中所述信噪比大于预设阈值的波段以及预设矿物光谱库中的反射率光谱，构建稀疏解混光谱库；

S4：根据所述稀疏解混光谱库，对所述含水矿物识别区进行稀疏解混，得到火星表面含水矿物的丰度。

在一个具体的例子中，所述步骤S1，包括图1中未示出的步骤S11至S14：

S11：将所述火星高光谱图像分割为多个预设大小的子区域，并计算各子区域中各波段的噪声强度。

S12：根据所述各子区域中各波段的噪声强度，计算所述火星高光谱图像各波段噪声强度。

S13：计算所述火星高光谱图像各波段的信噪比。

在本实施中，每个子区域包含的所有像元的平均值作为该子区域的信号强度；所有子区域信号强度的均值作为火星高光谱图像的信号强度；根据火星高光谱图像每个波段的信号强度和噪声强度，计算火星高光谱图像各波段的信噪比。

S14：选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演。

在一个具体的例子中，所述步骤S11，包括：

将所述火星高光谱图像分割为多个4×4大小的子区域，计算各子区域中各波段每个像元的残差，公式如下：

其中，r_i,j,k为子区域中第i行j列第k波段的像元残差，x_i,j,k为子区域中第i行j列第k波段的像元值，其中，1≤k≤N，N为所述火星高光谱图像包含的波段个数。当1<k<N时，由x_i,j,k-1和x_i,j,k+1线性拟合得到；当k＝1时，由x_i,j,2和x_i,j,3线性拟合得到；当k＝N时，由x_i,j,N-1和x_i,j,N-2线性拟合得到；

计算所述各子区域中各波段所有像元的残差的标准差，得到各子区域中各波段的噪声强度。

在一个具体的例子中，所述步骤S12，包括：

根据所述各子区域中各波段的噪声强度，确定各波段所有子区域的噪声强度中的最大值和最小值；

在各波段的噪声强度中的最大值和最小值之间，划分各波段的噪声强度为m个噪声区间，m为预设值；本实施例中m＝60；每个波段的噪声强度均划分为60个区间。

对于每个波段，将包含子区域数量最多的噪声区间中噪声强度的均值作为该波段的噪声强度。

在一个具体的例子中，图2示出了本实施例中S13得到的CRISM高光谱图像的信噪比，在S14中，选择信噪比大于100的的波段进行含水矿物的识别及定量反演。

在一个具体的例子中，在所述步骤S2中，所述预设的光谱特征参数包括：在波长1900nm处水的吸收特征BD₁₉₀₀、在波长2100nm处水的吸收特征BD₂₁₀₀、在2300nm波长处反射率的下降特征D₂₃₀₀以及在波长2290nm处的波峰特征SINDEX；

所述BD₁₉₀₀＝1－((R₁₉₃₀)/(a×R₁₈₅₀+b×R₂₀₄₆))，其中，a＝1－b，b＝(1930－1850)/(2046－1850)；

所述BD₂₁₀₀＝1－((R₂₁₃₂)/(a×R₁₉₃₀+b×R₂₂₅₀))，其中，a＝1－b，b＝(2132－1930)/(2250－1930)；

所述D₂₃₀₀＝1－(R₂₂₉₀+R₂₃₃₀+R₂₃₃₀)/(R₂₁₄₀+R₂₁₇₀+R₂₂₁₀)；

所述SINDEX＝1－((a×R₂₁₂₀+b×R₂₄₀₀)/(R₂₂₉₀))，其中，a＝1－b，b＝(2290-2120)/(2400-2120)；

式中下标表示波长，R表示反射率。如R₁₉₀₀表示1900nm波长处的反射率。

在一个具体的例子中，在所述步骤S2中，所述含水矿物识别模型为：

(BD₁₉₀₀＞t1)or(BD₂₁₀₀＞t2)or(D₂₃₀₀＞t3)or(SINDEX＞t4)

其中，t1、t2、t3和t4为预设常数。

本实施例中，设定t1＝t2＝t3＝0.02,t4＝0.015。

在一个具体的例子中，所述步骤S3，包括图1中未示出的步骤S31至S32：

S31：根据S1中所述信噪比大于预设阈值的波段的光谱范围和波长位置，对预设矿物光谱库中的反射率光谱进行重采样，以使所述反射率光谱与所述信噪比大于预设阈值的波段的光谱范围和波长位置一致；

S32：去除经过S31处理的所述预设矿物光谱库中任意两条光谱角小于预设最小光谱角的光谱的其中一条，得到稀疏解混光谱库。本实施例中最小光谱角设定为2.5度。

在一个具体的例子中，在所述步骤S31中，所述重采样的公式如下：

R′_c＝a×R_l+b×R_r

其中，R′_c为重采样后得到的波长c处的反射率，波长c为S1中所述信噪比大于预设阈值的波段图像的波长位置，R_l和R_r分别为所述预设矿物光谱库中在波长c左右两侧距离波长c最近的波长，即l和r处的反射率值，a＝(r-c)/(r-l)，b＝1-a。

在一个具体的例子中，在所述步骤S4中，所述进行稀疏解混的公式如下：

其中，y为所述火星高光谱图像的反射率，A为所述稀疏解混光谱库，λ表示预设的正则化参数(本实施例中设定为10^-5)，m为光谱库中的光谱个数，s_j为光谱库中与第j个光谱对应的端元丰度，采用SUnSAL算法进行矿物丰度的求解。

上述实施例中的火星表面含水矿物定量反演方法针对火星含水矿物具有丰度含量低、分布地域零散、背景矿物不确定或未知等特点，采用目标识别与稀疏解混相结合的方法进行含水矿物的定量反演，目标识别有助于缩小含水矿物的潜在分布范围，提高丰度反演的运算效率，增加含水矿物丰度反演的准确性；稀疏解混算法不需要对火星矿物类别和端元数量进行专家预判断，而是基于完备光谱库逐像元自动进行矿物类别识别和丰度反演。能够进行有效的火星表面含水矿物定量丰度反演。该方法可推广至深空探测各星体矿物定量反演。

本领域的技术人员能够理解，尽管在此所述的一些实施例包括其它实施例中所包括的某些特征而不是其它特征，但是不同实施例的特征的组合意味着处于本发明的范围之内并且形成不同的实施例。

虽然结合附图描述了本发明的实施方式，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下做出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。

Claims (8)

1.一种火星表面含水矿物定量反演方法，其特征在于，包括：

S0、获取火星高光谱图像；

S1：计算所述火星高光谱图像各波段的信噪比，并选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演；

S2：根据预设的光谱特征参数，建立含水矿物识别模型，得到含水矿物识别区；

S3：根据S1中所述信噪比大于预设阈值的波段以及预设矿物光谱库中的反射率光谱，构建稀疏解混光谱库；

S4：根据所述稀疏解混光谱库，对所述含水矿物识别区进行稀疏解混，得到火星表面含水矿物的丰度；

S2中，所述预设的光谱特征参数包括：在波长1900nm处水的吸收特征BD₁₉₀₀、在波长2100nm处水的吸收特征BD₂₁₀₀、在2300nm波长处反射率的下降特征D₂₃₀₀以及在波长2290nm处的波峰特征SINDEX；

所述BD₁₉₀₀＝1－((R₁₉₃₀)/(a×R₁₈₅₀+b×R₂₀₄₆))，其中，a＝1－b，b＝(1930－1850)/(2046－1850)；

所述BD₂₁₀₀＝1－((R₂₁₃₂)/(a×R₁₉₃₀+b×R₂₂₅₀))，其中，a＝1－b，b＝(2132－1930)/(2250－1930)；

所述D₂₃₀₀＝1－(R₂₂₉₀+R₂₃₃₀+R₂₃₃₀)/(R₂₁₄₀+R₂₁₇₀+R₂₂₁₀)；

所述SINDEX＝1－((a×R₂₁₂₀+b×R₂₄₀₀)/(R₂₂₉₀))，其中，a＝1－b，b＝(2290-2120)/(2400-2120)；

式中下标表示波长，R表示反射率。

2.如权利要求1所述的方法，其特征在于，所述步骤S1，包括：

S11：将所述火星高光谱图像分割为多个预设大小的子区域，并计算各子区域中各波段的噪声强度；

S12：根据所述各子区域中各波段的噪声强度，计算所述火星高光谱图像各波段噪声强度；

S13：计算所述火星高光谱图像各波段的信噪比；

S14：选择信噪比大于预设阈值的波段进行含水矿物的识别及定量反演。

3.如权利要求2所述的方法，其特征在于，所述步骤S11，包括：

将所述火星高光谱图像分割为多个w×h大小的子区域，计算各子区域中各波段每个像元的残差，公式如下：

$r_{i,j,k}=x_{i,j,k}-{\hat{x}}_{i,j,k}$

其中，r_i,j,k为子区域中第i行j列第k波段的像元残差，x_i,j,k为子区域中第i行j列第k波段的像元值，其中，1≤k≤N，N为所述火星高光谱图像包含的波段个数，当1<k<N时，由x_i,j,k-1和x_i,j,k+1线性拟合得到；当k＝1时，由x_i,j,2和x_i,j,3线性拟合得到；当k＝N时，由x_i,j,N-1和x_i,j,N-2线性拟合得到；其中，w和h为预设正整数；

计算所述各子区域中各波段所有像元的残差的标准差，得到各子区域中各波段的噪声强度。

4.如权利要求2所述的方法，其特征在于，所述步骤S12，包括：

根据所述各子区域中各波段的噪声强度，确定各波段所有子区域的噪声强度中的最大值和最小值；

在各波段的噪声强度中的最大值和最小值之间，划分各波段的噪声强度为m个噪声区间，m为预设值；

对于每个波段，将包含子区域数量最多的噪声区间中噪声强度的均值作为该波段的噪声强度。

5.如权利要求1所述的方法，其特征在于，在所述步骤S2中，所述含水矿物识别模型为：

(BD₁₉₀₀＞t1)or(BD₂₁₀₀＞t2)or(D₂₃₀₀＞t3)or(SINDEX＞t4)其中，t1、t2、t3和t4为预设常数。

6.如权利要求1所述的方法，其特征在于，所述步骤S3，包括：

S31：根据S1中所述信噪比大于预设阈值的波段的光谱范围和波长位置，对预设矿物光谱库中的反射率光谱进行重采样，以使所述反射率光谱与所述信噪比大于预设阈值的波段的光谱范围和波长位置一致；

S32：去除经过S31处理的所述预设矿物光谱库中任意两条光谱角小于预设最小光谱角的光谱的其中一条，得到稀疏解混光谱库。

7.如权利要求6所述的方法，其特征在于，在所述步骤S31中，所述重采样的公式如下：

R′_c＝a×R_l+b×R_r

其中，R′_c为重采样后得到的波长c处的反射率，波长c为S1中所述信噪比大于预设阈值的波段的波长位置，R_l和R_r分别为所述预设矿物光谱库中在波长c左右两侧距离波长c最近的波长，即l和r处的反射率值，a＝(r-c)/(r-l)，b＝1-a。

8.如权利要求1所述的方法，其特征在于，在所述步骤S4中，所述进行稀疏解混的公式如下：

$min\left\{\frac{1}{2}\right||y-{\mathrm{As}}_j|{|}_2^2+\mathrm{\&lambda;}\left|\right|s_j|{|}_1,j=1,2,...,m\}$

$\begin{array}{cccc}s.t.& s_j\&GreaterEqual;0,& j=1,2,...,m,& \underset{j=1}{\overset{m}{\&Sigma;}}{s}_j=1\end{array}$

其中，y为所述火星高光谱图像的反射率，A为所述稀疏解混光谱库，λ表示预设的正则化参数，m为光谱库中的光谱个数，s_j为光谱库中与第j个光谱对应的端元丰度，采用SUnSAL算法进行矿物丰度的求解。