CN102721650B - 基于特征指标的矿物成分遥感信息提取方法及装置 - Google Patents

Abstract

本发明提供了一种基于特征指标的矿物成分遥感信息提取方法和装置，其中所述方法包括如下步骤：A、对遥感图像进行预处理；B、对经过预处理的遥感图像的像元的特征指标进行计算；C、建立归方程；D、根据实际样品反推出所述回归方程中a、b、c、d和e的值；E、根据步骤C得到的回归方程计算像元处的矿物成分含量。本发明可以应用于矿产勘查，对矿产勘查工作起到节省时间、节约人力物力的作用。

Description

基于特征指标的矿物成分遥感信息提取方法及装置

本发明受国家自然科学基金项目（编号：41102205）、中央级公益性科研院所基本业务费项目（编号：K1207）、地质大调查项目（编号：1212011120994）、国家重大基础研究计划项目（编号：2007CB411406-06）、国家科技支撑专题（编号：2006BAB07B00-7-2）联合资助。

技术领域

本发明涉及一种地质遥感信息分析方法和装置，特别是在遥感信息中分析矿物成分含量的方法及装置。

背景技术

矿产资源是影响国家安全和国家战略布局重要的战略资源。随着开采时间的延续，找矿难度不断增大，目前急需找矿勘查技术的创新。地球化学已经成为了一个系统科学，而随着技术进步和其他学科的交叉联系，地球化学的新的分支领域不断的产生，地球化学已经有了多个外延。遥感技术作为一门新兴的高新技术，与地球化学的结合形成了遥感地球化学这样一个分支。遥感地球化学的产生既丰富了遥感应用研究领域，又拓宽了地球化学研究领域。

遥感地球化学是以物质电磁波理论为基础，借助遥感技术获取数据，研究化学元素在地表或其它行星表面的分布、含量及迁移的科学，它的特点是快速、大范围获取数据。除能获取地球静态参数外，遥感地球化学还可以测量地球动态参数。

目前利用遥感进行找矿已经取得了大量的成果。这些主要是利用地球化学晕形成指示矿床或矿化带的蚀变粘土矿物、碳酸盐、氢氧化物、硫酸盐等矿物，可以绘制矿物成分分布图及丰度图，而地下油气藏在地表形成的地球化学场异常特征的波谱相应也能圈出勘探靶区。王晋年等利用航空红外细分光谱仪（FIMS）进行了遥感地球化学找矿研究。根据从FIMS图像上提取的反射率曲线，在2.175μm处有绢云母A1-OH吸收峰;在2.29μm处有绿泥石Mg-OH吸收峰。然后利用这2个波段的光谱吸收指数(SAI)进行矿物吸收鉴别分类，识别出蚀变玄武岩、蚀变凝灰岩、凝灰质砂岩、玄武岩、戈壁等5类岩石。将该技术推广到未知区，在新疆博孜阿特与博格特区找到2条稳定的金矿化蚀变带（王晋年，郑兰芬，童庆禧.成像光谱图像光谱吸收鉴别模型与矿物填图研究[J].环境遥感，1996，11(1):20-32.）。

遥感地球化学在从遥感信息中分析矿物成分含量时，主要利用波段反射率、DN值来反演矿物成分含量。因为矿物的吸收特征主要是分子键振动形成的吸收特征或者是电子跃迁形成的吸收特征。主要的吸收特征（主要指吸收谷的波长）与分子振动和电子跃迁有关，这与分子晶体结构有关，如羟基、碳酸盐基团的吸收特征，和铁离子、锰离子的吸收特征，而与矿物成分含量关系不大。但是矿物成分含量对分子吸收特征是有影响的，如高铝白云母和低铝白云母，其吸收特征有一定的差异性，铝含量的增加吸收谷向左偏移。在月球上的铁、钛元素含量与月壤光谱特征有关系，Lucey（Lucey Paul G.，Blewett David T.，and Hawke B.Ray.Mapping the FeO and TiO2 content of the lunar surface with multispectral imagery[J].Journal of Geophysical Research，1998，103(2):3679～3699）给出了经验公式。这就说明矿物成分含量对光谱有一定的影响。当矿物成分含量具有一定差异时，相应地在波段反射率或DN值上也会有差异体现出来，但还未发现矿物成分含量与波段反射率或DN值之间的直接的函数关系。因此，当遥感信息中波段反射率或DN值的波动较小时，直接利用波段反射率、DN值来反演矿物成分含量得到的结果往往是不准确的，有时甚至会产生较大的错误。

发明内容

为了提高从遥感信息中分析矿物成分含量的准确性，本发明提供了一种基于特征指标的矿物成分遥感信息提取方法和装置。

本发明的技术方案如下：

基于特征指标的矿物成分遥感信息提取方法，包括如下步骤：

A、对遥感图像进行预处理；

B、对经过步骤A处理的遥感图像的像元的特征指标进行计算，所述特征指标包括：

波段相关斜率SL：SL=（B_i-B_j）/（B_j-B_m）；

波段斜率指数SI：SI=e^((Bi-Bm)/Bi)；

波段相关吸收深度指数SD：SD=e^{((Bi-Bk)/(Bi-Bm))}；

波段相关多项式指数SP：SP=e^{((Bi-Bm)/Bi+(Bj-Bk)/Bm)}；

其中，B_i为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最大值对应的波段；B_j为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第二大值对应的波段；B_k为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第三大值对应的波段；B_m为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最小值对应的波段；

C、建立回归方程：

Con_ij=a×SL_ij+b×SI_ij+c×SD_ij+d×SP_ij+e

其中，Con_ij为像元处的矿物成分含量，SL_ij为所述像元的波段相关斜率，SI_ij为所述像元的波段斜率指数，SD_ij为所述像元的波段相关吸收深度指数，SP_ij为所述像元的波段相关多项式指数，a、b、c、d为系数，e为回归分析截距；

D、根据实际样品反推出所述回归方程中a、b、c、d和e的值；

E、根据步骤C得到的回归方程计算像元处的矿物成分含量。

优选的，所述基于特征指标的矿物相关遥感信息提取方法，还包括步骤F：基于步骤E得到的所述像元处的矿物成分含量形成彩色图像。

优选的，在步骤C中还包括对所述回归方程进行回归检验的步骤。

优选的，步骤A所述预处理包括如下步骤：

选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合形成第一优化数据，或对属于ETM数据的所述原始数据选择7、4和3波段组合形成第一优化数据。

优选的，步骤A所述预处理还包括如下步骤：在所述第一优化数据中选取研究区域：分割第一优化数据，判断分割区内每个波段的直方图是否为正态分布，如果是正态分布，则将该分割区作为研究区域；如果不是正态分布，则将所述分割区划分为多个次级区域，使得所述次级区域内的每个光谱波段的直方图为正态分布，并将所述次级区域作为所述研究区域。

优选的，步骤A所述预处理还包括如下步骤：将所述原始数据中不包含全部波段数据的边界区域的去除步骤。

优选的，在步骤A所述预处理中，对所述第一优化数据还按如下公式进行处理：

公式为：

$R=\frac{\mathrm{pi}\×L\×d^2}{\mathrm{Esun}\×\cos \left(A\right)}$

其中：R为行星反射率；pi=3.14;d表示所述原始数据成像当天的日地距离；A是太阳高度角；Esun值为大气层外相应波长的太阳光谱辐照度；L是辐亮度，可通过下式求出：

L=gain×DN+bias

其中：gain为增益，bias为偏移；

对于热红外数据需要温度和比发射率的分离，公式如下：

$T=\frac{c_2}{\mathrm{\λ}\ln (1+\frac{\mathrm{\ϵ}{c}_1}{{\mathrm{\λ}}^{5}R})}$

$\mathrm{\ϵ}=\frac{R{\mathrm{\λ}}^5(e^{\frac{c_2}{\mathrm{\λT}}}-1)}{c_1}$

其中：T为温度，λ为波长，ε为比发射率，c₁，c₂为常数，c₁=3.74818×10^-4Wμm²，c₂=1.43878×10⁴Kμm；R为光谱辐射亮度，可以通过下面公式计算出：

$R={\mathrm{LMIN}}_{\mathrm{\λ}}+\left(\frac{{\mathrm{LMAX}}_{\mathrm{\λ}}-{\mathrm{LMIN}}_{\mathrm{\λ}}}{\mathrm{QCALMAX}}\right)\mathrm{QCAL}$

QCAL是数据的实际辐射，LMIN_λ是QCAL=0时的光谱辐射值，LMAX_λ是在QCAL=QCALMAX的光谱辐射值，QCALMAX是数据的图像辐射值。

优选的，步骤A所述预处理还包括如下步骤：即对所述波段直方图上存在的偏离正态分布的值或者在所述波段直方图上有突出的值进行去除。

优选的，步骤A所述预处理还包括掩膜步骤：即进行数值设置，使得所述第一优化数据中不需要进行研究的数据被去除。

一种基于特征指标的矿物相关遥感信息提取装置，包括如下部件：

预处理部件：用于对遥感图像进行预处理；

特征指标计算部件：用于对经过所述预处理部件处理的遥感图像的像元的特征指标进行计算，所述特征指标包括：

波段相关斜率SL：SL=（B_i-B_j）/（B_j-B_m）；

波段斜率指数SI：SI=e^((Bi-Bm)/Bi)；

波段相关吸收深度指数SD：SD=e^{((Bi-Bk)/(Bi-Bm))}；

波段相关多项式指数SP：SP=e^{((Bi-Bm)/Bi+(Bj-Bk)/Bm)}；

其中，B_i为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最大值对应的波段；B_j为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第二大值对应的波段；B_k为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第三大值对应的波段；B_m为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最小值对应的波段；

回归方程部件：用于建立如下回归方程，

Con_ij=a×SL_ij+b×SI_ij+c×SD_ij+d×SP_ij+e

其中，Con_ij为像元处的矿物成分含量，SL_ij为所述像元的波段相关斜率，SI_ij为所述像元的波段斜率指数，SD_ij为所述像元的波段相关吸收深度指数，SP_ij为所述像元的波段相关多项式指数，a、b、c、d为系数，e为回归分析截距；

系数推演部件：用于根据实际样品反推出所述回归方程中a、b、c、d和e的值；

结果获取部件：用于根据步骤C得到的回归方程计算像元处的矿物成分含量。

优选的，所述基于特征指标的矿物相关遥感信息提取装置还包括图像形成部件：基于结果获取部件得到的所述像元处的矿物成分含量形成彩色图像。

优选的，所述回归方程部件还执行如下步骤：对所述回归方程进行回归检验。

优选的，所述预处理部件用于执行如下步骤：

选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合形成第一优化数据，或对属于ETM数据的所述原始数据选择7、4和3波段组合形成第一优化数据。

优选的，所述预处理部件还执行如下步骤：在所述第一优化数据中选取研究区域：分割第一优化数据，判断分割区内每个波段的直方图是否为正态分布，如果是正态分布，则将该分割区作为研究区域；如果不是正态分布，则将所述分割区划分为多个次级区域，使得所述次级区域内的每个光谱波段的直方图为正态分布，并将所述次级区域作为所述研究区域。

优选的，所述预处理部件还执行如下步骤：将所述原始数据中不包含全部波段数据的边界区域的去除。

优选的，所述预处理部件还执行如下步骤：对所述第一优化数据按如下公式进行处理：

公式为：

$R=\frac{\mathrm{pi}\×L\×d^2}{\mathrm{Esun}\×\cos \left(A\right)}$

其中：R为行星反射率；pi=3.14;d表示所述原始数据成像当天的日地距离；A是太阳高度角；Esun值为大气层外相应波长的太阳光谱辐照度；L是辐亮度，可通过下式求出：

L=gain×DN+bias

其中：gain为增益，bias为偏移；

对于热红外数据需要温度和比发射率的分离，公式如下：

$T=\frac{c_2}{\mathrm{\λ}\ln (1+\frac{\mathrm{\ϵ}{c}_1}{{\mathrm{\λ}}^{5}R})}$

$\mathrm{\ϵ}=\frac{R{\mathrm{\λ}}^5(e^{\frac{c_2}{\mathrm{\λT}}}-1)}{c_1}$

其中：T为温度，λ为波长，ε为比发射率，c₁，c₂为常数，c₁=3.74818×10^-4Wμm²，c₂=1.43878×10⁴Kμm；R为光谱辐射亮度，可以通过下面公式计算出：

$R={\mathrm{LMIN}}_{\mathrm{\λ}}+\left(\frac{{\mathrm{LMAX}}_{\mathrm{\λ}}-{\mathrm{LMIN}}_{\mathrm{\λ}}}{\mathrm{QCALMAX}}\right)\mathrm{QCAL}$

QCAL是数据的实际辐射，LMIN_λ是QCAL=0时的光谱辐射值，LMAX_λ是在QCAL=QCALMAX的光谱辐射值，QCALMAX是数据的图像辐射值。

优选的，所述预处理部件还执行如下步骤：对所述波段直方图上存在的偏离正态分布的值或者在所述波段直方图上有突出的值进行去除。

本发明的技术效果：

直接利用波段反射率、DN值来反演矿物成分含量没能把波段反射率或DN值的区别进行放大和细节区分，这导致反演结果不准确。如果在反演矿物成分含量时能够把波段反射率或DN值的区别进行放大，则对于矿物成分含量的信息提取就能够更准确。尽管波段反射率或DN值与矿物成分含量之间没有直接的函数关系，仍然可以利用统计学的方法来进行信息的提取。

本发明的矿物成分遥感信息提取方法或装置中，利用波段相关斜率、波段斜率指数、波段相关吸收深度指数和波段相关多项式指数这四个特征指标，以指数的方式将波段反射率或DN值的区别进行放大，从而得到了准确性更高的矿物成分含量反演结果。

附图说明

图1为本发明矿物成分遥感信息提取方法的一个实例的流程图。

图2为本发明矿物成分遥感信息提取方法的另一个实例的流程图。

图3为铜元素含量与ASTER对应光谱之间的关系。

图4为一个地区铜元素分布的地球化学图。

图5为现有技术反演的图3所示地区铜元素分布的地球化学图。

图6为本发明方法反演的图3所示地区铜元素分布的地球化学图。

图7为本发明矿物成分遥感信息提取装置的一个实例的原理图。

图8为本发明矿物成分遥感信息提取装置的另一个实例的原理图。

具体实施方式

在对本发明的技术方案进行详细描述之前，对使用的名词检查进行说明：

DN值（Digital Number）是遥感影像像元亮度值，记录的地物的灰度值。无单位，是一个整数值，值大小与传感器的辐射分辨率、地物发射率、大气透过率和散射率等有关。

以下结合附图对本发明的技术方案进行详细说明。

图1显示了本发明基于特征指标的矿物成分遥感信息提取方法的一个实例。

首先，对遥感图像进行预处理步骤是对遥感图像这一原始数据进行必要的处理，以简化后续数据处理步骤。例如，去除遥感图像中明显的有云遮挡部分、水体部分等不能获得有用结果的部分。

其次，进行特征指标计算，本发明名称中的特征指标就是指本步骤中计算的特征指标。前一步骤得到的预处理结果在本步骤得到进一步处理，即对遥感图像中像元的如下特征指标进行计算：

波段相关斜率SL：SL=（B_i-B_j）/（B_j-B_m）；

波段斜率指数SI：SI=e^((Bi-Bm)/Bi)；

波段相关吸收深度指数SD：SD=e^{((Bi-Bk)/(Bi-Bm))}；

波段相关多项式指数SP：SP=e^{((Bi-Bm)/Bi+(Bj-Bk)/Bm)}；

其中，B_i为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最大值对应的波段；B_j为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第二大值对应的波段；B_k为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第三大值对应的波段；B_m为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最小值对应的波段。

图3是上述样本波段光谱图的一个例子，是通过如下方法获得：在遥感图像所示区域内进行实地采样，选取铜元素含量不同的样本，并按采集样本的坐标做了ASTER图像对应位置的采样光谱，基本上是按1ppm铜含量递进，分别为9.8ppm、10.8ppm、11.9ppm、13.1ppm、14.1ppm、15.1ppm、16.1ppm和17.1ppm铜元素含量。可以看出，随着铜元素含量的增多，第7和第8波段与第4波段之间的差值在增大，也就是光谱曲线的斜率在逐步增大。但这种对应关系不是特别明显，因此采用上述特征指标，就可以把铜元素含量对光谱线斜率的影响进行放大。

第三，建立回归方程：

Con_ij=a×SL_ij+b×SI_ij+c×SD_ij+d×SP_ij+e

其中，Con_ij为像元处的矿物成分含量，SL_ij为所述像元的波段相关斜率，SI_ij为所述像元的波段斜率指数，SD_ij为所述像元的波段相关吸收深度指数，SP_ij为所述像元的波段相关多项式指数，a、b、c、d为系数，e为回归分析截距。

回归分析是处理变量的相关关系的一种数理统计分析。主要目的是在自变量和因变量之间没有严格、确定的函数关系中找出最能代表其关系的数学表达形式。在本步骤建立的回归方程就是利用这一原理在上述四个特征指标与矿物成分含量之间建立数学表达形式。

第四，反推回归方程系数。在这一步骤中，利用实地采集的样本的数据（所述样本已与遥感图像中的位置坐标对应）反推第三步骤中回归方程的系数，即a、b、c、d和e的值。理论上通过5个样本的数据就可以得到上述系数，但更多样本数据能够获得更准确的系数。

第五，得到矿物成分含量。在第四步骤中，已经得到了回归方程，在本步骤中，利用该回归方程对遥感图像中像元进行计算，得到该像元处矿物成分含量（例如铜的含量）。

图2显示了本发明矿物成分遥感信息提取方法的另一个实例。

第一步骤，遥感图像选择。

为了提高后续数据处理的准确性，必须对获得的遥感图像（原始数据）进行严格的筛选。要求遥感图像的时相尽可能地选择在植被发育弱、冰雪覆盖较少的季节，同时要求数据获取期间，该区域的天空云量较少。选择适合的遥感图像后，进一步地，要进行影像合成，一般对于ASTER数据一般选择6、3、1波段组合合成影像数据（即第一优化数据），ETM选择7、4、3波段组合合成影像数据（即第一优化数据）。这样可以明确看出植被、冰雪、云和水等地物。植被一般呈现绿色，冰雪是亮白色，或亮浅蓝色，云为亮白色，水为暗蓝色或暗黑色等。

第二步骤，研究区域选择。

这一步骤在所述第一优化数据中选取研究区域，研究区域作为后续数据处理的对象。研究区域选择的方法如下：分割第一优化数据，判断分割区内每个波段的直方图是否为正态分布，如果是正态分布，则将该分割区作为研究区域；如果不是正态分布，则将所述分割区划分为多个次级区域，使得所述次级区域内的每个光谱波段的直方图为正态分布，并将所述次级区域作为所述研究区域。

对这一步骤的方法进行进一步说明：采用滑动窗口直方图统计法，设计一个M×N的窗口，判断这个窗口内的各波段直方图是否近似为正态分布，也可以采用目估法。

图示目估法：设定一个M×N的窗口，这个窗口内遥感图像某一个波段像元x_j，k（j=1，m;k=1，n）区间为[x₀，x_n]，统计这个窗口内的直方图，公式如下：

其中i∈[0，n]，x_j，k=x_i为逻辑操作，目估直方图是否符合正态分布。

计算判断法：设定一个M×N的窗口，这个窗口内遥感图像某一个波段像元值为x_j，k（j=1，m;k=1，n），像元均值为

标准差为σ，利用偏度系数和峰度系数进行判断。

偏度系数满足公式：

$g=\frac{m\×n}{(m\×n-1)(m\×n-2)}\underset{j=1,k=1}{\overset{m,n}{\mathrm{\Σ}}}{\left(\frac{x_{j,k}-\stackrel{\‾}{x}}{\mathrm{\σ}}\right)}^3\∈[-{\mathrm{\ϵ}}_1,{\mathrm{\ϵ}}_1]$

其中ε₁为给定一个很小的正数。

峰度系数满足公式：

$f=\frac{m\×n(m\×n+1)}{(m\×n-1)(m\×n-2)(m\×n-3)}\underset{j=1,k=1}{\overset{m,n}{\mathrm{\Σ}}}{\left(\frac{x_{j,k}-\stackrel{\‾}{x}}{\mathrm{\σ}}\right)}^4-\frac{3{(m\×n-1)}^2}{(m\×n-2)(m\×n-3)}\∈[-{\mathrm{\ϵ}}_2,{\mathrm{\ϵ}}_2]$ 其中ε₂为给定一个很小的正数。

第三步骤，是否有边界区域。

在这一步骤中要判断第一优化数据中是否有边界区域。边界区域是指在缺失某个波段的信息的区域，或者说各波段信息不重合区域，这一区域往往处于遥感图像的边界。

第四步骤，去除边界区域。

如果在第三步骤中判定存在边界区域，在本步骤中需要去除边界区域。具体的方法就是对一个区域是否含有每一个波段信息进行判断，如果是，则附值为1；如果为否，则附值为0，生成一个二值图像，然后把每一个波段二值图像相乘形成一个新的二值图像，最后把每个波段与二值图像相乘，这样就去除了边界信息。具体公式如下：

$y_i=\underset{i=0}{\overset{n}{\∩}}(x_i>0)\·x_i$

其中n指所使用的遥感图像波段总数，i=1，…，n，x_i和y_i指分别指i波段去除波段前后的值。

第五步骤，进一步处理。

到第四步骤为止，获得的图像仍存在几何变形、遥感器增益和偏移等，在经本步骤处理后得到带有坐标信息的行星反射率图像，具体公式为：

$R=\frac{\mathrm{pi}\×L\×d^2}{\mathrm{Esun}\×\cos \left(A\right)}$

其中：R为行星反射率；pi=3.14；d表示所述原始数据成像当天的日地距离；A是太阳高度角；Esun值为大气层外相应波长的太阳光谱辐照度；L是辐亮度，可通过下式求出：

L=gain×DN+bias

其中：gain为增益，bias为偏移；

对于热红外数据需要温度和比发射率的分离，公式如下：

$T=\frac{c_2}{\mathrm{\λ}\ln (1+\frac{\mathrm{\ϵ}{c}_1}{{\mathrm{\λ}}^{5}R})}$

$\mathrm{\ϵ}=\frac{R{\mathrm{\λ}}^5(e^{\frac{c_2}{\mathrm{\λT}}}-1)}{c_1}$

其中：T为温度，λ为波长，ε为比发射率，c₁，c₂为常数，c₁=3.74818×10^-4Wμm²，c₂=1.43878×10⁴Kμm；R为光谱辐射亮度，可以通过下面公式计算出：

$R={\mathrm{LMIN}}_{\mathrm{\λ}}+\left(\frac{{\mathrm{LMAX}}_{\mathrm{\λ}}-{\mathrm{LMIN}}_{\mathrm{\λ}}}{\mathrm{QCALMAX}}\right)\mathrm{QCAL}$

QCAL是数据的实际辐射，LMIN_λ是QCAL=0时的光谱辐射值，LMAX_λ是在QCAL=QCALMAX的光谱辐射值，QCALMAX是数据的图像辐射值。

对于ASTER数据来说，由于存在三个不同传感器获取的数据，分别率不一样，首先需要通过插值合并的归一化处理，并生成包含所使用波段的单一文件，插值分为几个像元合并成一个像元，采用像元和均值法和一个像元分成直接分成几个像元。然后再对像元按第五步骤中的公式进行处理。

第六步骤，是否有“尖锐”信息。

在本步骤，对于前序步骤处理得到的研究区域的直方图进行判断，判断是否存在“尖锐”信息。所谓“尖锐”信息是指图像上有一些数量很少而又集中产生某种特定像元的信息，包括某一波段集中产生亮像元或暗像元的信息和大量集中的透射率小于5％的背景噪声。判断依据主要是每一个波段的近似正态分布的直方图上两端是否存在着少量偏离正态分布的值或者直方图上有突出的值。

第七步骤，去除“尖锐”信息。

在判断出“尖锐”信息后，采用如下方法进行去除：

对某一波段集中产生亮像元或暗像元的“尖锐”信息利用高端切割或低端切割来去除，对集中产生透射率小于5％的背景噪声采用波段比值法来去除，公式如下：

或

前一个是高端切割和低端切割公式，后一个是波段比值切割公式，其中，i=0，...，n，n指所使用的遥感图像波段总数，x_i和y_i指分别指i波段去除“尖锐”信息前和后的的波段值，a，b∈[1，…，n]，C_a，C_b是常数，x_a，x_b是原始a、b波段对应的值。这个公式的目的就是：给定一个约束的条件，使得这个条件大于或小于某一个数值的图像保留下来，其他的全部被赋值为零。

第八步骤，生成掩膜。

掩膜的作用是确定哪些数据需要参加计算，哪些不需要。本步骤生成一个掩膜，其实就是生成一个二值图像，0值表示不需要参与计算的数据，1值表示需要参与计算的数据。生成掩膜一般采用波段的逻辑计算方法。

$\mathrm{\τ}=y_i\⊗c_0$ 或

τ为生成的掩膜，y_i为去除“尖锐”信息的数据，x_i，x_j为原始数据，

为关系运算符（包括<、≤、≯、>、≥、≮、≠、＝等），⊕为数学运算符（包括±、×、÷等），c₁，c₂，c₃为常数。

以上八个步骤可以看作是遥感数据的预处理过程。

第九步骤，特征指标计算。

这一步骤与图1中特征指标计算步骤内容相同，不再赘述。

第十步骤，建立回归方程。

这一步骤与图1中建立回归方程步骤内容相同，不再赘述。

第十一步骤，回归检验。

在本步骤对前一步骤建立的回归方程进行回归检验。可以采用如下三种回归检验方法之一：

1）总变异

$=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{(y_j-{\hat{y}}_j)}^2+\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{({\hat{y}}_j-\stackrel{\&OverBar;}{y})}^2=\mathrm{SS}+\mathrm{MS}$

SS为回归平方和，是回归值与均值差值的平方和，反映了自变量X的变化引起Y的波动，自由度df_回=n（n为自变量个数）。MS为残差平方和，是实测值y_j与回归值

差值的平方和，是由试验误差及其它因素引起的，自由度df_残=m-n-1。总变异自由度为m-1。如果观测值给定，总变异确定，可以用SS和MS来衡量回归效果，SS越大回归效果越显著，MS越大回归效果不好。

2）指标决定系数

为了检验总的回归效果，定义了无量纲指标-决定系数R²来表示，

R²反映回归离差对总变异的贡献比例。R=R^1/2称为复相关系数，反映全部自变量与因变量的相关程度。R²和R值越大，回归效果越好。

上述位总体回归效果检验，不能说明每个自变量x₁，x₂，…，x_n对因变量y都重要，有些自变量对因变量可能不起作用，或者作用为其他自变量所代替，这就需要把这些自变量从回归方程中剔除，建议每个自变量x_i是否显著，假设H0：a_i=0，i=1，2，…n。

3）F值检验

在H0：a_i=0假设下，

$F(n,m-n-1)=\frac{\mathrm{SS}/n}{\mathrm{MS}/(m-n-1)}$

对给定置信度α，从F值分布表中查与β对应的临界值F_β，如果|F_i|〉F_β，拒绝假设H0，认为n个自变量总体回归效果显著，反之，总体回归效果不显著。

第十二步骤，反推回归方程系数。

这一步骤与图1中反推回归方程系数步骤内容相同，不再赘述。

第十三步骤，得到矿物成分含量。

这一步骤与图1中得到矿物成分含步骤内容相同，不再赘述。

第十四步骤，生成彩色图像。

在获得了矿物成分含量数据之后，将这一数据通过假色彩方式叠加到遥感图像中，使得人们能更直观地观察到矿物成分含量分布情况。

采用本发明的矿物成分遥感信息提取方法，利用波段相关斜率、波段斜率指数、波段相关吸收深度指数和波段相关多项式指数这四个特征指标，以指数的方式将波段反射率或DN值的区别进行放大，从而得到了准确性更高的矿物成分含量反演结果。图4至图5通过实例说明了相对于现有直接利用波段反射率、DN值来反演矿物成分含量和本发明矿物成分遥感信息提取方法的不同效果。图4显示的是一个地区实测的铜元素的地球化学，图5是同一地区遥感地球化学（采用现有方法反演结果），图6是同一地区遥感地球化学（采用本发明方法反演结果）。可以看到，与图4的实测结果比对，图5显示的利用现有方法反演结果中，两个圆圈部分的信息是缺失的。而利用本发明的方法获得的结果（图6）则不存在这样的缺失。

图7显示了本发明基于特征指标的矿物相关遥感信息提取装置的一个实例。包括顺序连接的如下部件：预处理部件、特征指标计算部件、回归方程部件、系数推演部件和结果获取部件。

预处理部件用于执行图1中“对遥感图像预处理”步骤。

特征指标计算部件用于执行图1中“特征指标计算”步骤。

回归方程部件用于执行图1中“建立回归方程”步骤。

系数推演部件用于执行图1中“反推回归方程系数”步骤。

结果获取部件用于执行图1中“得到矿物成分含量”步骤。

图8显示了本发明基于特征指标的矿物相关遥感信息提取装置的另一个实例。包括顺序连接的如下部件：预处理部件、特征指标计算部件、回归方程部件、系数推演部件、结果获取部件和图像形成部件。

预处理部件用于执行图2中第一至第八步骤。

特征指标计算部件用于执行图2中第九步骤。

回归方程部件用于执行图2中第十和第十一步骤。

系数推演部件用于执行图2中第十二步骤。

结果获取部件用于执行图2中第十三步骤。

图像形成部件用于执行图2中第十四步骤。

值得注意的是，以上所述仅为本发明的较佳实施例，并非因此限定本发明的专利保护范围，本发明还可以采用等同技术进行替换。故凡运用本发明的说明书及图示内容所作的等效变化，或直接或间接运用于其他相关技术领域均同理皆包含于本发明所涵盖的范围内。

Claims (17)

1.基于特征指标的矿物成分遥感信息提取方法，其特征在于包括如下步骤：

A、对遥感图像进行预处理；

B、对经过步骤A处理的遥感图像的像元的特征指标进行计算，所述特征指标包括：

波段相关斜率SL：SL=（B_i-B_j）/（B_j-B_m）；

波段斜率指数SI：SI=e^((Bi-Bm)/Bi)；

波段相关吸收深度指数SD：SD=e^{((Bi-Bk)/(Bi-Bm))}；

波段相关多项式指数SP：SP=e^{((Bi-Bm)/Bi+(Bj-Bk)/Bm)}；

其中，B_i为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最大值对应的波段；B_j为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第二大值对应的波段；B_k为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第三大值对应的波段；B_m为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最小值对应的波段；

C、建立回归方程：

Con_ij=a×SL_ij+b×SI_ij+c×SD_ij+d×SP_ij+e

其中，Con_ij为像元处的矿物成分含量，SL_ij为所述像元的波段相关斜率，SI_ij为所述像元的波段斜率指数，SD_ij为所述像元的波段相关吸收深度指数，SP_ij为所述像元的波段相关多项式指数，a、b、c、d为系数，e为回归分析截距；

D、根据实际样品反推出所述回归方程中a、b、c、d和e的值；

E、根据步骤C得到的回归方程计算像元处的矿物成分含量。

2.根据权利要求1所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：还包括步骤F：基于步骤E得到的所述像元处的矿物成分含量形成彩色图像。

3.根据权利要求2所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：在步骤C中还包括对所述回归方程进行回归检验的步骤。

4.根据权利要求1至3之一所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：步骤A所述预处理包括如下步骤：

选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合形成第一优化数据，或对属于ETM数据的所述原始数据选择7、4和3波段组合形成第一优化数据。

5.根据权利要求4所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：步骤A所述预处理还包括如下步骤：在所述第一优化数据中选取研究区域：分割第一优化数据，判断分割区内每个波段的直方图是否为正态分布，如果是正态分布，则将该分割区作为研究区域；如果不是正态分布，则将所述分割区划分为多个次级区域，使得所述次级区域内的每个光谱波段的直方图为正态分布，并将所述次级区域作为所述研究区域。

6.根据权利要求5所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：步骤A所述预处理还包括如下步骤：将所述第一优化数据中不包含全部波段数据的边界区域的去除步骤。

7.根据权利要求6所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：在步骤A所述预处理中，对所述第一优化数据还按如下公式进行处理：

公式为：

$R=\frac{\mathrm{pi}\&times;L\&times;d^2}{\mathrm{Esun}\&times;\cos \left(A\right)}$

其中：R为行星反射率；pi=3.14;d表示所述原始数据成像当天的日地距离；A是太阳高度角；Esun值为大气层外相应波长的太阳光谱辐照度；L是辐亮度，可通过下式求出：

L=gain×DN+bias

其中：gain为增益，bias为偏移；

对于热红外数据需要温度和比发射率的分离，公式如下：

$T=\frac{c_2}{\mathrm{\&lambda;}\ln (1+\frac{{\mathrm{\&epsiv;c}}_1}{{\mathrm{\&lambda;}}^{5}R})}$

$\mathrm{\&epsiv;}=\frac{R{\mathrm{\&lambda;}}^5(e^{\frac{c_2}{\mathrm{\&lambda;T}}}-1)}{c_1}$

其中：T为温度，λ为波长，ε为比发射率，c₁，c₂为常数，c₁=3.74818×10^-4Wμm²，c₂=1.43878×10⁴Kμm；R为光谱辐射亮度，可以通过下面公式计算出：

$R={\mathrm{LMIN}}_{\mathrm{\&lambda;}}+\left(\frac{{\mathrm{LMAX}}_{\mathrm{\&lambda;}}-{\mathrm{LMIN}}_{\mathrm{\&lambda;}}}{\mathrm{QCALMAX}}\right)\mathrm{QCAL}$

QCAL是数据的实际辐射，LMIN_λ是QCAL=0时的光谱辐射值，LMAX_λ是在QCAL=QCALMAX的光谱辐射值，QCALMAX是数据的图像辐射值。

8.根据权利要求7所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：步骤A所述预处理还包括如下步骤：即对所述波段直方图上存在的偏离正态分布的值或者在所述波段直方图上有突出的值进行去除。

9.根据权利要求8所述基于特征指标的矿物成分遥感信息提取方法，其特征在于：步骤A所述预处理还包括掩膜步骤：即进行数值设置，使得所述第一优化数据中不需要进行研究的数据被去除。

10.一种基于特征指标的矿物相关遥感信息提取装置，其特征在于包括如下部件：

预处理部件：用于对遥感图像进行预处理；

特征指标计算部件：用于对经过所述预处理部件处理的遥感图像的像元的特征指标进行计算，所述特征指标包括：

波段相关斜率SL：SL=（B_i-B_j）/（B_j-B_m）；

波段斜率指数SI：SI=e^((Bi-Bm)/Bi)；

波段相关吸收深度指数SD：SD=e^{((Bi-Bk)/(Bi-Bm))}；

波段相关多项式指数SP：SP=e^{((Bi-Bm)/Bi+(Bj-Bk)/Bm)}；

其中，B_i为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最大值对应的波段；B_j为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第二大值对应的波段；B_k为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的第三大值对应的波段；B_m为所述像元的反射率或DN值，对应的波段为样本波段光谱图中反射率或DN值的最小值对应的波段；

回归方程部件：用于建立如下回归方程，

Con_ij=a×SL_ij+b×SI_ij+c×SD_ij+d×SP_ij+e

其中，Con_ij为像元处的矿物成分含量，SL_ij为所述像元的波段相关斜率，SI_ij为所述像元的波段斜率指数，SD_ij为所述像元的波段相关吸收深度指数，SP_ij为所述像元的波段相关多项式指数，a、b、c、d为系数，e为回归分析截距；

系数推演部件：用于根据实际样品反推出所述回归方程中a、b、c、d和e的值；

结果获取部件：用于根据步骤C得到的回归方程计算像元处的矿物成分含量。

11.根据权利要求10所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：还包括图像形成部件：基于结果获取部件得到的所述像元处的矿物成分含量形成彩色图像。

12.根据权利要求11所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：所述回归方程部件还执行如下步骤：对所述回归方程进行回归检验。

13.根据权利要求10至12之一所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：所述预处理部件用于执行如下步骤：

选取成像时间为春季或秋季的遥感图像作为原始数据；对属于ASTER数据的所述原始数据选择6、3和1波段组合形成第一优化数据，或对属于ETM数据的所述原始数据选择7、4和3波段组合形成第一优化数据。

14.根据权利要求13所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：所述预处理部件还执行如下步骤：在所述第一优化数据中选取研究区域：分割第一优化数据，判断分割区内每个波段的直方图是否为正态分布，如果是正态分布，则将该分割区作为研究区域；如果不是正态分布，则将所述分割区划分为多个次级区域，使得所述次级区域内的每个光谱波段的直方图为正态分布，并将所述次级区域作为所述研究区域。

15.根据权利要14所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：所述预处理部件还执行如下步骤：将所述第一优化数据中不包含全部波段数据的边界区域的去除。

16.根据权利要求15所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：所述预处理部件还执行如下步骤：对所述第一优化数据按如下公式进行处理：

公式为：

$R=\frac{\mathrm{pi}\&times;L\&times;d^2}{\mathrm{Esun}\&times;\cos \left(A\right)}$

其中：R为行星反射率；pi=3.14;d表示所述原始数据成像当天的日地距离；A是太阳高度角；Esun值为大气层外相应波长的太阳光谱辐照度；L是辐亮度，可通过下式求出：

L=gain×DN+bias

其中：gain为增益，bias为偏移；

对于热红外数据需要温度和比发射率的分离，公式如下：

$T=\frac{c_2}{\mathrm{\&lambda;}\ln (1+\frac{{\mathrm{\&epsiv;c}}_1}{{\mathrm{\&lambda;}}^{5}R})}$

$\mathrm{\&epsiv;}=\frac{R{\mathrm{\&lambda;}}^5(e^{\frac{c_2}{\mathrm{\&lambda;T}}}-1)}{c_1}$

其中：T为温度，λ为波长，ε为比发射率，c₁，c₂为常数，c₁=3.74818×10^-4Wμm²，c₂=1.43878×10⁴Kμm；R为光谱辐射亮度，可以通过下面公式计算出：

$R={\mathrm{LMIN}}_{\mathrm{\&lambda;}}+\left(\frac{{\mathrm{LMAX}}_{\mathrm{\&lambda;}}-{\mathrm{LMIN}}_{\mathrm{\&lambda;}}}{\mathrm{QCALMAX}}\right)\mathrm{QCAL}$

QCAL是数据的实际辐射，LMIN_λ是QCAL=0时的光谱辐射值，LMAX_λ是在QCAL=QCALMAX的光谱辐射值，QCALMAX是数据的图像辐射值。

17.根据权利要求16所述基于特征指标的矿物相关遥感信息提取装置，其特征在于：所述预处理部件还执行如下步骤：对所述波段直方图上存在的偏离正态分布的值或者在所述波段直方图上有突出的值进行去除。

Images