一种基于奇异性理论的高光谱矿物异常信息提取方法与流程

本发明涉及高光谱遥感矿物识别技术领域，尤其涉及一种基于奇异性理论的高光谱矿物异常信息提取方法。

背景技术：

随着遥感技术的发展，以航空高光谱为代表的定量化遥感技术已广泛应用于各个领域，成为地球观测的重要手段。与多光谱不同，高光谱具有丰富的光谱信息可用于进行矿物信息的识别。

目前，高光谱矿物识别的主要方法可归纳为两类：基于光谱相似性匹配的提取方法和基于光谱特征参量的提取方法。光谱相似性匹配法对于光谱特征较相似的矿物很难区分，基于光谱特征参量法需制定十分详细的提取规则，且实现过程复杂。上述两种方法均不可缺少专家知识的参与和标准矿物光谱的分析，所以，单纯使用某一种方法在应用过程中往往难以实现对矿物信息的精细化提取。

从地质应用角度来讲，矿物信息提取的目的在于提供与成矿相关异常信息和线索。然而，与成矿相关的矿物信息往往占少数，大部分矿物信息反映的是区域上的背景矿物信息。如何通过光谱分析反应矿物在空间尺度上的异常信息是目前高光谱地质应用的难点。

技术实现要素：

本发明的目的在于提供一种基于奇异性理论的高光谱矿物异常信息提取方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明采用的技术方案如下：

一种基于奇异性理论的高光谱矿物异常信息提取方法，包括如下步骤：

s1，获取高光谱影像数据和标准光谱数据，依据矿物光谱特征范围，分别截取高光谱影像及标准光谱的特征区间反射率数据，得到高光谱影像反射率数据和标准光谱反射率数据；

s2，对s1中得到的反射率数据进行连续统去除，得到去除连续统后的光谱数据；

s3，对s2获取的去除连续统后的影像光谱数据进行高斯拟合，得到拟合曲线；

s4，根据拟合曲线的最低点获得影像光谱的特征吸收位置；

s5，计算s3中高斯拟合后的影像光谱数据和s2去除连续统后的标准光谱数据的拟合度；

s6，如果s4中获得的影像光谱的特征吸收位置与标准光谱中的特征吸收位置一致，且s5中计算获得的拟合度大于设定的阈值，则确定该像元含目标矿物；

s7，在s6中得到的包含目标矿物信息的整幅光谱影像上逐像元计算奇异性指数，得到矿物异常信息分布图。

优选地，s2中采用的公式为：

rcr＝r/rc

其中，rcr为去除连续统后的光谱数据，r为原始光谱数据，rc为原始光谱的连续统。

优选地，s5中，所述匹配度按照如下公式进行计算：

其中，r²表示拟合度，取值范围在[0，1]之间，r²越接近于1表示拟合度越高，yi表示去除连续统后的标准光谱数据，是yi的平均值，为高斯拟合后的影像光谱数据。

优选地，s4中还包括步骤，利用如下公式计算所述特征吸收位置的吸收深度：

depthmin＝1-rcrmin

其中，depthmin表示吸收深度，rcrmin表示拟合曲线最低点的数值。

优选地，s6中还包括步骤，如果s4中获得的影像光谱的特征吸收位置与标准光谱中的特征吸收位置一致，且s5中计算获得的拟合度大于设定的阈值，则确定该像元含目标矿物，吸收位置对应的吸收深度则代表目标矿物的相对丰度。

优选地，s6还包括步骤：基于专家知识对影像光谱像元进行二次核对，如果通过，则完成目标矿物信息的提取，否则，继续调节阈值直至提取结果满足应用需求。

优选地，s7具体为：设定一系列滑动窗口，窗口大小l依次递增，lmin＝l1＜l2…＜ln＝lmax，在矿物信息的整幅影像上某一点滑动窗口，并计算每个窗口内的平均矿物密度，通过拟合矿物密度与窗口大小l之间的函数关系，按照如下公式计算奇异性指数：

<ρ(a)>＝ca^-δα

其中，c为分形密度；δα为奇异性指数；<ρ(a)>为一片大小为l，面积为a的区域中某种物矿密度。

本发明的有益效果是：本发明提供的基于奇异性理论的高光谱矿物异常信息提取方法，首先获得目标矿物在光谱影像上的特征吸收位置，然后，通过将影像上的特征吸收位置与标准光谱上的特征吸收位置进行比较，如果两者一致，而且，高光谱影像数据与标准光谱数据之间的匹配度达到设定的阈值，则认为该光谱影像像元含目标矿物，从而完成矿物信息提取；最后，根据提取出的矿物信息，在整幅的光谱影像上逐像元计算奇异性指数，得到矿物异常信息分布图。所以，采用本发明提供的方法，强化了对矿物弱信息的提取，而这类具有“突变”特征的矿物异常信息相对于大量背景矿物信息往往与成矿关系更密切，更加利于提取矿物信息。

附图说明

图1是本发明提供的基于奇异性理论的高光谱矿物异常信息提取方法流程示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例提供了一种基于奇异性理论的高光谱矿物异常信息提取方法，包括如下步骤：

s1，获取高光谱影像数据和标准光谱数据，依据矿物光谱特征范围，分别截取高光谱影像及标准光谱的特征区间反射率数据，得到高光谱影像反射率数据和标准光谱反射率数据；

s2，对s1中得到的反射率数据进行连续统去除，得到去除连续统后的光谱数据；

s3，对s2获取的去除连续统后的影像光谱数据进行高斯拟合，得到拟合曲线；

s4，根据拟合曲线的最低点获得影像光谱的特征吸收位置；

s5，计算s3中高斯拟合后的影像光谱数据和s2去除连续统后的标准光谱数据的拟合度；

s6，如果s4中获得的影像光谱的特征吸收位置与标准光谱中的特征吸收位置一致，且s5中计算获得的拟合度大于设定的阈值，则确定该像元含目标矿物；

s7，在s6中得到的包含目标矿物信息的整幅光谱影像上逐像元计算奇异性指数，得到矿物异常信息分布图。

从地质应用角度来讲，矿物信息提取的目的在于提供与成矿相关的异常信息和线索。然而，与成矿相关的矿物信息往往占少数，大部分矿物信息反映的是区域上的背景矿物信息。大量学者基于非线性理论和复杂性理论论证了成矿异常范围内与成矿相关的矿物和元素的分布规律服从幂律分布，而与成矿无关的背景区矿物和元素的分布服从正态分布。作为非线性理论的延伸，局部奇异性的根本原理是定量表征成矿物质富集和聚集的时空分布规律。本发明中，基于该思想，通过区域矿物奇异性指数来反映矿物在空间尺度上的异常信息，进而反映某个局部位置矿物信息相对整体背景的变化量，强化了对矿物弱信息的提取，而这类具有“突变”特征的矿物异常信息相对于大量背景矿物信息往往与成矿关系更密切。

本发明提供的上述方法中，首先通过对高光谱影像数据进行处理，获得目标矿物在光谱影像上的特征吸收位置，然后，通过将影像上的特征吸收位置与标准光谱上的特征吸收位置进行比较，如果两者一致，而且，高光谱影像数据与标准光谱数据之间的匹配度达到设定的阈值，则认为该光谱影像像元即为目标矿物，从而完成矿物信息提取；

最后，根据提取出的矿物信息，在整幅的光谱影像上逐像元计算奇异性指数，由于奇异性指数可以反映矿物在空间尺度上的异常信息，所以，通过计算奇异性指数，就可以得到矿物异常信息分布图。

所以，采用本发明提供的方法，强化了对矿物弱信息的提取，而这类具有“突变”特征的矿物异常信息相对于大量背景矿物信息往往与成矿关系更密切，更加利于提取矿物信息。

其中，s1中，获取的高光谱影像数据，可以进行辐射定标、大气校正等预处理，标准光谱数据可以从usgs光谱库或tsg光谱库中获取。

利用专家知识分析标准光谱(主要参考usgs及tsg光谱库)的诊断特征位置，在此基础上确定矿物光谱特征范围，该范围可以是一段也可以是多段，对于常见的矿物，三段光谱特征范围基本可以准确描述其光谱特征。表1中给出了几种矿物的光谱特征范围，如，白云母矿物的光谱特征范围为三段，2145-2285nm区间，2285-2405nm区间和2405-2528nm区间。

表1基于高光谱的常见矿物信息提取规则表

在矿物光谱特征范围内，截取高光谱影像数据和标准光谱数据，得到对应的特征区间反射率数据。

s2中，对s1得到的反射率影像和标准光谱均进行连续统去除操作，可以消除地形、光照等对光谱的影响。

其中，连续统为反射光谱曲线中反射峰之间的连线。连续统去除即以反射光谱除以连续统光谱，处理后各个反射峰将变为1，反射峰之间的数值均小于1。

由于高光谱数据获取过程中影响因素较多，噪声存在会使得光谱数据往往呈锯齿状，因此，本发明中，在s3中，对s2获取的去连续统后的影像光谱进行高斯拟合，得到优化后的光影像谱曲线。

而拟合曲线的最低点即为特征吸收位置。

根据影像光谱数据和标准光谱数据的拟合度以及特征吸收位置，提取目标矿物信息，即判断是否含目标矿物。

其中，标准光谱的特征吸收位置是通过利用专家知识对标准光谱进行分析总结的信息提取规则，如表1中所示，白云母矿物在2145-2285nm区间内的2190/2210/2225nm位置均有吸收特征，在2285-2405nm区间内的2348nm附近及2405-2528nm区间的2440nm附近均有吸收特征。

如果影像光谱的吸收位置符合提取规则，且拟合度大于设定阈值的像元则包含目标矿物。

最后，s7中，通过计算各像元的奇异性指数，得到矿物异常信息分布图，完成矿物异常信息的提取。

本发明的一个实施例中，s2中可以采用的公式为：

rcr＝r/rc

其中，rcr为去除连续统后的光谱数据，r为原始光谱数据，rc为原始光谱的连续统。

本发明的一个优选的实施例中，s5中，所述拟合度可以按照如下公式进行计算：

其中，r²表示拟合度，取值范围在[0，1]之间，r²越接近于1表示拟合度越高，yi表示去除连续统后的标准光谱数据，是yi的平均值，为高斯拟合后的影像光谱数据。

本发明中，s4中还可以包括步骤，利用如下公式计算所述特征吸收位置的吸收深度：

depthmin＝1-rcrmin

其中，depthmin表示吸收深度，rcrmin表示拟合曲线最低点的数值。

根据吸收深度可以判断目标矿物的含量。

本发明中，如果s4中获得的影像光谱的特征吸收位置与标准光谱中的特征吸收位置一致，且s5中计算获得的拟合度大于设定的阈值，则确定该像元含目标矿物，吸收位置对应的吸收深度则代表目标矿物的相对丰度。。

本发明中，s6中还可以包括步骤：基于专家知识对影像光谱像元进行二次核对，如果通过，则完成目标矿物信息的提取，否则，继续调节阈值直至提取结果满足应用需求。

本发明中，采用标准光谱进行拟合度匹配的提取方法在一定程度上简化了单纯基于光谱特征进行矿物信息提取的步骤，而基于专家知识的矿物信息提取规则进一步加强了单纯利用光谱匹配方式进行矿物信息提取的精度。

本发明实施例中，s7具体可以为：设定一系列滑动窗口，窗口大小l依次递增，lmin＝l1＜l2…＜ln＝lmax，在矿物信息的整幅影像上某一点滑动窗口，并计算每个窗口内的平均矿物密度，通过拟合矿物密度与窗口大小l之间的函数关系，按照如下公式计算奇异性指数：

<ρ(a)>＝ca^-δα

其中，c为分形密度；δα为奇异性指数；<ρ(a)>为一片大小为l，面积为a的区域中某种物矿密度。

通过滑动窗口的方式进行奇异性指数计算，过程简单易用实现，通过分析不同尺度下矿物的空间变化规律，得到某个局部位置矿物信息的变化量，可以强化对矿物信息的提取，突出区域上具有“突变”特征的矿物异常信息。

具体实施例：

为了更好地说明本发明的方法和步骤，以某地区hymap机载高光谱数据为例，进行白云母矿物异常信息填图试验。

(1)试验所用的设备为图形工作站，规格型号为dellprecisiont7600，操作系统为windows7(64位)，cpu为2.66ghz，内容为32gb，硬盘为1t。

(2)具体步骤如下：

步骤一：反射率区间数据的获取。

首先对hymap数据进行辐射定标、大气校正等预处理获得反射率数据，从usgs光谱库中获取白云母标准光谱数据。依据表1中白云母光谱特征范围，截取高光谱影像及标准光谱的2145-2285nm，2285-2405nm和2405-2528nm三段区间的反射率数据。

步骤二：光谱连续统去除。

连续统去除即以反射光谱除以连续统光谱，处理后各个反射峰将变为1，反射峰之间的数值均小于1，具体计算过程如下式所示：

rcr＝r/rc

其中rcr为去连续统处理后的光谱数据，r为原始光谱数据，rc为原始光谱的连续统。按照上述公式对步骤一获得的区间反射率影像和标准光谱进行去连续统处理。

步骤三：矿物信息提取。

首先对步骤二获取的去连续统后的影像光谱进行高斯拟合，得到优化后的影像光谱曲线。通过寻找拟合曲线的最低点，得到特征吸收位置，计算该位置的吸收深度，公式如下：

吸收深度计算公式：

depthmin＝1-rcrmin

depthmin其中表示吸收深度，rcrmin表示拟合曲线最低点的数值。

分别对优化后的影像光谱与和去连续统后的区间标准光谱进行匹配度计算，计算公式如下：

其中r²表示拟合度，取值范围在[0，1]之间，r²越接近于1表示拟合度越高，yi表示去连续统后的光谱数据，是yi的平均值，为拟合值。

当拟合度大于设定的阈值情况下，基于影像光谱和标准光谱的吸收位置的数值(影像光谱的吸收位置通过上述的计算得到，标准光谱的吸收位置通过对标准光谱进行分析得到，如表1中所示，如白云母矿物在2145-2285nm区间内的2190/2210/2225nm位置均有吸收特征在2285-2405nm区间内的2348nm附近及2405-2528nm区间的2440nm附近均有吸收特征)，可判断该光谱区间内吸收特征是否与目标矿物的特征相匹配，如果该光谱区间内吸收位置与目标矿物的标准光谱的吸收位置一致，则该光谱区间内吸收特征与目标矿物的特征匹配，可提取出满足提取条件的光谱信息，最终基于专家知识对提取的结果进行二次核对完成矿物信息的提取。

步骤四：奇异性指数计算。

假设在一片大小为l，面积为a的区域，其某种矿物密度可以表示为<ρ(a)>，根据奇异性理论，面积和密度具有如下关系：

<ρ(a)>＝ca^-δα

其中，c是分形密度；δα表示奇异性指数。

在二维空间中δα＝2-α，δα数值可反映当前区域矿物密度相对于周边富集或亏损的程度。

先设定一系列滑动窗口，窗口大小l依次递增，在本案例中窗口大小设置以像元为单位，分别是2pix，6pix，10pix，14pix，18pix，22pix，26pix。在影像上某一点滑动窗口，并计算每个窗口内的平均矿物密度，通过拟合矿物密度与窗口大小l之间的函数关系，可求得该点的奇异性指数δα和分形密度c。在步骤三中得到的矿物信息影像上逐像元计算可完成整幅影像的奇异性指数求解，得到矿物异常信息分布图。

δα＞0的区域空间分布范围是局限的，其分形维数往往小于2，然而，矿物密度相对不变或者贫化的地段(即δα≤0)空间分布是广泛的，其分形维数往往接近于2，矿物没有明显的富集的区域(δα≤0)即是背景矿物区。当δα＞0，说明随着面积的减小矿物密度急剧增加，这样的富集区域即是矿物信息异常区。则可以根据奇异性指数得到矿物异常信息分布图。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明提供的基于奇异性理论的高光谱矿物异常信息提取方法，首先获得目标矿物在光谱影像上的特征吸收位置，然后，通过将影像上的特征吸收位置与标准光谱上的特征吸收位置进行比较，如果两者一致，而且，高光谱影像数据与标准光谱数据之间的匹配度达到设定的阈值，则认为该光谱影像像元含目标矿物，从而完成矿物信息提取；最后，根据提取出的矿物信息，在整幅的光谱影像上逐像元计算奇异性指数，得到矿物异常信息分布图。所以，采用本发明提供的方法，强化了对矿物弱信息的提取，而这类具有“突变”特征的矿物异常信息相对于大量背景矿物信息往往与成矿关系更密切，更加利于提取矿物信息。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。