CN106251310A - 一种多光谱遥感异常信息提取方法 - Google Patents

Abstract

本发明公开了一种多光谱遥感异常信息提取方法，该方法包括以下步骤：对获取的遥感影像进行预处理；对预处理影像进行掩膜处理；对掩膜处理后的影像进行比值‑主成分转换，得到若干主成分影像；对每种主成分影像进行分形分析，实现异常信息的提取。本发明具有一定的普适性和可操作性，能够最大限度的剔除各类非异常和假异常，突出与资源环境有关的真异常(如植被指数异常、蚀变异常等)，并实现异常分级，因此具有重要的理论和实践意义，可服务于地勘、矿山、环保等政府部门和相关企业。

Description

一种多光谱遥感异常信息提取方法

技术领域

本发明涉及遥感找矿与环境评价领域，具体涉及一种多光谱遥感异常信息提取方法。

背景技术

遥感技术是一项包含宏观、动态、快速、准确等特点以及相应优势的高科技技术方法。资源与环境遥感探测是资源环境的理论、方法与遥感的理论、技术、手段相结合的交叉学科，是全球变化研究中重要的理论基础和技术方法。异常识别与信息提取是资源环境遥感的主要研究课题之一。如具有异常光谱特征的地表围岩蚀变可作为热液矿床的主要找矿标志，且蚀变越强，成矿可能性越高；蚀变范围越宽，矿化规模越大；即或矿体隐伏，只要有一定面积的蚀变异常出露，遥感影像都有可能检测到。又如，在一些植被茂盛、风化层较厚、山陡谷深、原岩露头少的重金属矿区或污染区往往存在植物光谱异常，如叶面光谱反射率异常增高、红边位移、吸收深度降低、植被指数异常等。

遥感数据主要反映地表物质的光谱信息，由于土壤、植被、大气、冰川、水体、地形等干扰因素的存在，加之遥感影像波谱分辨率和空间分辨率的制约以及地质环境的复杂性和多变性，异常信息在遥感影像中主要表现为弱信息。多光谱遥感探测虽然有一些成功的典型案例存在，但它们不代表一般情况。它们的成功可能是由于这些异常信息在遥感影像中不表现为弱信息，如蚀变面积大、污染范围广、混合像元效应弱等。事实上，大部分出露于地表或埋深较浅的地质资源已被勘查，现行的异常信息提取技术，如掩膜、波段比值、主成分变换、门限法、分形等，有时很难适应资源环境遥感探测的需求。弱信息是多光谱遥感探测中的基本困难，前人在这方面亦开展了大量研究，但大多比较零散，仍然存在深化研究的空间。因此，开发出具有一定普适性的、可重复的、操作简便的多光谱遥感(弱)异常提取技术已成为当务之急。

发明内容

为解决上述问题，本发明公开了一种多光谱遥感异常信息提取方法，包括以下步骤：

步骤1，遥感影像选择及预处理,得到预处理影像；

步骤2，对预处理影像进行掩膜处理；

步骤3，对掩膜处理后的遥感影像进行遥感蚀变信息增强处理；

其特征在于，还包括如下步骤：

步骤4，对遥感影像增强处理后的波段进行分形求和计算，得到异常增强影像，提取异常信息。

进一步地，当遥感影像中存在植物异常，还包括：

步骤5：对植被异常增强影像进行混合筛分处理，剔除假异常。

进一步地，步骤1中所述的预处理是指对选取的遥感影像依次进行辐射定标、大气校正、几何校正和正射校正处理。

进一步地，步骤4所述的对主成分转化后的影像波段进行分形求和计算是指：

步骤41，通过下式构建分形模型，得到N(r)-r散点图：

N(r)＝Cr^-D(r＞0)

其中：r表示特征尺度，即影像DN值，(C＞0)为比例系数,D＞0为分维数，N(r)＝N(≥r)表示大于等于r的数目或和数；

步骤42，通过下式将N(r)-r散点图拟合为多段直线，每条直线对应于一个无标度区：

logN(r)＝-Dlog(r)+logC

步骤43，步骤42拟合的多段无标度区，应满足如下公式：

$\begin{array}{c}RSS=\underset{i=1}{\overset{i_0}{\Σ}}{\[\log N\left(r_i\right)+D_1\log \left(r_i\right)-\log C_1\]}^2+\underset{i=i_0+1}{\overset{i_1}{\Σ}}{\[\log N\left(r_i\right)+D_2\log \left(r_i\right)-\log C_2\]}^2\\ +\underset{i=i_1+1}{\overset{i_2}{\Σ}}{\[\log N\left(r_i\right)+D_3\log \left(r_i\right)-\log C_3\]}^2+\mathrm{...}\\ +\underset{i=i_n+1}{\overset{N}{\Σ}}\log N\left(r_i\right)+D_{n+2}\log \left(r_i\right)-\log C_{n+2}{\]}^2\→Min\end{array}$

其中，D₁、D₂、D₃…D_n+2分别表示多段无标度区中第1个、第2个、第3个…第n+2个无标度区的分维数，C₁、C₂、C₃…C_n+2为第1个、第2个、第3个…第n+2个无标度区所对应的比例系数，其中表示从小到大排列的灰度值,i＝1,2,3…,n＝0,1,2…，D_n所对应的像元数记为N_n，RSS为拟合残差的平方；

其中，D_n和D_n+1之间的拐点记为T_n，该拐点T_n对应的灰度值即为异常灰度值

进一步地，步骤5中所述的对植被异常增强影像进行混合筛分处理是指：

步骤51，获取植被异常增强影像的灰度直方图；

步骤52，采用混合筛分模型将影像灰度直方图分割为若干个正态分布概率密度曲线，取不同正态分布概率密度曲线的交点为阈值，找到均值最小的正态分布或均值最大的正态分布。

与现有技术相比，本发明具有以下技术效果：

本发明提出了一种快速有效提取多光谱遥感异常信息的技术方案，具有一定的普适性和可操作性，能够最大限度的剔除各类非异常和假异常，突出与资源环境有关的真异常(如植被指数异常、蚀变异常等)，并实现异常分级，因此具有重要的理论和实践意义，可服务于地勘、矿山、环保等政府部门和相关企业。

附图说明

图1为多光谱遥感异常信息提取总体流程图；

图2为实施例1中PC1和PC3影像logN(r)-log(r)分形模式图；

图3为实施例2中RVI影像logN(r)-log(r)分形模式图；

图4为实施例2中针对D₁(图4)对应像元的混合筛分模式图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

实施例1：德兴地区(28°50′-29°10′N,117°30′-118°00′E)遥感蚀变异常提取，具体包括以下步骤：

步骤1，遥感影像选择及预处理，预处理包括：辐射定标、大气校正、几何校正和正射校正处理；

在遥感影像选取时，尽量避开云、雾、植被等因素的影响；

本实例选取的遥感影像参数如下：传感器类型：Landsat 8OLI，空间分辨率：30m，轨道号120/40，获取时间：2016年3月28日02:37:51.75时，太阳方位角：134.45°，太阳高度角：55.78°，日地距离：0.9981577AU，云覆盖率：<0.02％，影像质量：9。

步骤2，对预处理影像进行掩膜处理；

掩膜对象为山地阴影、道路、水体等干扰地物，即2/7比值影像中大于1.549081的像素(或图斑)以及NDVI(归一化差分植被指数)图像中小于0.103273的图斑，掩膜区面积占总研究区的1.6268％。

步骤3，对掩膜处理后的遥感影像进行遥感蚀变信息增强处理，本实施例通过ENVI软件进行遥感蚀变信息增强处理，即比值-主成分转换，得到蚀变信息增强影像；

本实例选择5/4(突出植被)、4/2(突出铁染异常)和6/7(突出羟基异常)三幅比值影像进行主成分转换(计算时使用协方差矩阵、进行正向旋转)，结果发现第一主成分(PC1)的方差贡献率为93.88％，第二(PC 2)、三(PC3)主成分的贡献率分别为4.67％和1.45％；第一主成分主要反映羟基蚀变，其特征向量为-0.987732，第二主成分主要反映植被信息，其特征向量为-0.960501，第三主成分主要反映铁染异常，其特征向量为-0.960835。

步骤4，分别对蚀变信息增强的影像波段进行分形求和计算，实现蚀变异常信息的提取；

其中，对蚀变信息增强的影像波段进行分形求和计算是指：

步骤41，通过下式构建主成分转化后影像的分形模型，得到N(r)-r散点图：

N(r)＝Cr^-D(r＞0)

其中：r表示特征尺度，即影像DN值，(C＞0)为比例系数,D＞0为分维数，N(r)＝N(≥r)表示大于等于r的数目或和数；

步骤42，通过下式将N(r)-r散点图拟合为多段直线，每条直线对应于不同地物信息：

logN(r)＝-Dlog(r)+logC

步骤43，一般而言，步骤42可拟合为两段或两段以上直线，其斜率的绝对值记为分维数D_n，n＝1,2,3…。若拟合为两段直线(即无标度区)，应满足如下公式：

$RSS=\underset{i=1}{\overset{i_0}{\&Sigma;}}{\&lsqb;logN\left(r_i\right)+D_{1}log\left(r_i\right)-log{C}_1\&rsqb;}^2+\underset{i=i_0+1}{\overset{N}{\&Sigma;}}{\&lsqb;logN\left(r_i\right)+D_{2}log\left(r_i\right)-log{C}_2\&rsqb;}^2\&RightArrow;Min$

其中D₁和D₂通常代表背景区和异常区的分维数，C₁和C₂为背景区和异常区比例系数，D₁和D₂之间的拐点记为T，拐点T对应的灰度值即异常灰度值(r_i表示从小到大排列的灰度值，i＝1,2,3…)。D₁和D₂所对应的像元数记为N₁和N₂。RSS为拟合残差的平方，两段直线拟合要求RSS取值最小。

步骤43，若步骤42可拟合为多个无标度区，则应满足如下公式：

$\begin{array}{c}RSS=\underset{i=1}{\overset{i_0}{\&Sigma;}}{\&lsqb;\log N\left(r_i\right)+D_1\log \left(r_i\right)-\log C_1\&rsqb;}^2+\underset{i=i_0+1}{\overset{i_1}{\&Sigma;}}{\&lsqb;\log N\left(r_i\right)+D_2\log \left(r_i\right)-\log C_2\&rsqb;}^2\\ +\underset{i=i_1+1}{\overset{i_2}{\&Sigma;}}{\&lsqb;\log N\left(r_i\right)+D_3\log \left(r_i\right)-\log C_3\&rsqb;}^2+\mathrm{...}\\ +\underset{i=i_n+1}{\overset{N}{\&Sigma;}}\log N\left(r_i\right)+D_{n+2}\log \left(r_i\right)-\log C_{n+2}{\&rsqb;}^2\&RightArrow;Min\end{array}$

其中，D₁、D₂、D₃…D_n+2分别表示多段无标度区中第1个、第2个、第3个…第n+2个无标度区的分维数，C₁、C₂、C₃…C_n+2为第1个、第2个、第3个…第n+2个无标度区所对应的比例系数，其中表示从小到大排列的灰度值,i＝1,2,3…,n＝0,1,2…，D_n所对应的像元数记为N_n，RSS为拟合残差的平方；

一般而言，D_1～2或D_n+1～n+2为代表异常的无标度区间。D_n和D_n+1之间的拐点记为T_n，该拐点T_n对应的灰度值即异常灰度值(n(<N)＝0,1,2…)。D_n所对应的像元数记为N_n，多段拟合同样要求RSS取值最小。

本实施例中，经主成分转换后，PC1和PC3图像各有1048572个像元值，大部分像元与区域矿化蚀变无关。虽然“门限法”所圈定的异常多为假异常或非异常，但有时在客观上起到了缩小靶区的作用，如对于PC1，本实例选取DN值大于0.2(约为均值+1倍方差)的像元进行分形分析；对于PC3，选取DN值小于-0.023(约等于均值-1倍方差)的像元进行分析，因此PC1影像中的待分析像元数降至119963个，PC3中的待分析像元数降至118666个。这里之所以选择一倍方差是为了避免遗漏有用信息。图1为Matlab平台下获得的PC1和PC3的分形模式图。

对于PC1，落入无标度区D₁中的像元数占总像元数的76.75％，但未见热液蚀变或矿床发育；如图2所示，D₃和D₄对应的异常图斑则与已知矿区完全吻合。对于PC3，落入D₁的像元数占总像元数的68.20％，但真正的异常为D₄和D₅，异常像元占0.80％。德兴矿床为亚洲最大的露天铜矿，蚀变异常显著而完整，因此分级异常提取效果较好，无需开展进一步统计分析。

其中所涉及的某些计算参数如下表。

表1所涉及的若干计算参数

实施例二：德兴地区(28°57′-29°04′N,117°34′-117°47′E)植被异常信息提取，具体包括以下步骤：

步骤1，遥感影像选取及预处理，预处理包括：辐射定标、大气校正、几何校正和正射校正；

为了突出植被异常，应尽量选取晚春或夏季影像。本实例选取的影像参数如下：传感器类型：Landsat 8 OLI，影像景号：LC81200402013223LGN00，轨道号：120/40，获取时间：2013年8月11日02:40:10.39时，空间分辨率：30m，太阳方位角：117.22°，太阳高度角：64.72°，日地距离：1.0135101AU，云覆盖率：<0.79％，影像质量：9。

步骤2，对预处理影像进行掩膜处理；

掩膜对象为居民地、道路、水体、裸地等非植被区，即2/7比值影像中大于1.493781的像素以及NDVI图像中小于0.18046的像素，掩膜区面积占总研究区的11.4928％。

步骤3，对掩膜处理后的遥感影像进行遥感蚀变信息增强处理，即通过ENVI软件进行波段比值计算；

本实例选用RVI(比值植被指数)进行异常提取，其中，RVI＝R/NIR，其中R对应OLI影像第4波段，NIR对应于OLI的5波段；

步骤4，对遥感影像增强处理后的波段进行分形求和计算，得到异常增强影像，提取异常信息

其中，对波段比值后的影像波段进行分形求和计算是指：

步骤41，通过下式构建分形模型，得到N(r)-r散点图：

N(r)＝Cr^-D(r＞0)

其中：r表示特征尺度，即影像DN值，(C＞0)为比例系数,D＞0为分维数，N(r)＝N(≥r)表示大于等于r的数目或和数；

步骤42，通过下式将N(r)-r散点图拟合为多段直线，每条直线对应于一个无标度区：

logN(r)＝-Dlog(r)+logC

步骤43，一般而言，步骤42可拟合为两段或两段以上直线，其斜率的绝对值记为分维数D_n，n＝1,2,3…。若拟合为两段直线(即无标度区)，应满足如下公式：

$RSS=\underset{i=1}{\overset{i_0}{\&Sigma;}}{\&lsqb;logN\left(r_i\right)+D_{1}log\left(r_i\right)-log{C}_1\&rsqb;}^2+\underset{i=i_0+1}{\overset{N}{\&Sigma;}}{\&lsqb;logN\left(r_i\right)+D_{2}log\left(r_i\right)-log{C}_2\&rsqb;}^2\&RightArrow;Min$

其中D₁和D₂通常代表背景区和异常区的分维数，C₁和C₂为背景区和异常区比例系数，D₁和D₂之间的拐点记为T，拐点T对应的灰度值即异常灰度值(r_i表示从小到大排列的灰度值，i＝1,2,3…)。D₁和D₂所对应的像元数记为N₁和N₂。RSS为拟合残差的平方，两段直线拟合要求RSS取值最小。

步骤43，步骤42拟合的多段无标度区，应满足如下公式：

$\begin{array}{c}RSS=\underset{i=1}{\overset{i_0}{\&Sigma;}}{\&lsqb;\log N\left(r_i\right)+D_1\log \left(r_i\right)-\log C_1\&rsqb;}^2+\underset{i=i_0+1}{\overset{i_1}{\&Sigma;}}{\&lsqb;\log N\left(r_i\right)+D_2\log \left(r_i\right)-\log C_2\&rsqb;}^2\\ +\underset{i=i_1+1}{\overset{i_2}{\&Sigma;}}{\&lsqb;\log N\left(r_i\right)+D_3\log \left(r_i\right)-\log C_3\&rsqb;}^2+\mathrm{...}\\ +\underset{i=i_n+1}{\overset{N}{\&Sigma;}}\log N\left(r_i\right)+D_{n+2}\log \left(r_i\right)-\log C_{n+2}{\&rsqb;}^2\&RightArrow;Min\end{array}$

其中，D₁、D₂、D₃…D_n+2分别表示多段无标度区中第1个、第2个、第3个…第n+2个无标度区的分维数，C₁、C₂、C₃…C_n+2为第1个、第2个、第3个…第n+2个无标度区所对应的比例系数，其中表示从小到大排列的灰度值,i＝1,2,3…,n＝0,1,2…，D_n所对应的像元数记为N_n，RSS为拟合残差的平方；

一般而言，D_1～2或D_n+1～n+2为代表异常的无标度区间。D_n和D_n+1之间的拐点记为T_n，该拐点T_n对应的灰度值即异常灰度值(n(<N)＝0,1,2…)。D_n所对应的像元数记为N_n，多段拟合同样要求RSS取值最小。

经证实，RVI＞0的像元一般为正常植被区或植被稀疏区，因此本实例计算只考虑RVI＜0的像元(共计50601个，占全区面积的16.70％)。图3为Matlab平台下获得的RVI的分形模式图，图中植被异常、假异常和非异常等对应于不同的分维数D_n(n＝1、2、3、4…)，亦即无标度区，不同分维数之间的拐点T_n(n＝1、2、3、4…)为区分不同地物的阈值，经证实，D₁为区域植被异常的反映，主要围绕着已知污染源(如露天采矿区、尾矿库、排污口等)分布，与实例一所圈定的铁染异常和羟基异常呈嵌套式接触关系，D₂主要为山地阴影，因此T₁＝-40.44730436即为圈定植被异常的分界点。

步骤5：对分形求和计算后得到的植被异常信息进行混合筛分处理，剔除假异常。

经查证，步骤4得到的异常信息仍然存在一些假异常(多为山地阴影区植被)，且不能进行异常分级，因此需要进行混合筛分处理。

所述的对异常影像进行混合筛分处理是指：

步骤51，获取异常信息的影像灰度直方图；

步骤52，采用混合筛分模型将影像灰度直方图分割为若干个正态分布，取不同成因总体的正态分布概率密度曲线的交点为阈值，一般而言，总体1(1代表均值最小的总体)或总体n(n代表均值最大的总体)即为异常。

不失一般性，可设样本的混合分布模型为：

$p(x;\mathrm{\&theta;})=\underset{i=1}{\overset{k}{\&Sigma;}}{a}_{i}f(x,{\mathrm{\&theta;}}_i)$

其中，参数集θ＝(k,θ₁,…,θ_k,α₁,…α_k)，k为分支的个数，f(x,θ_i)表示第i个分支的概率密度函数，θ_i为相应的参数，a_i为第i个分支的权重，并有 ${\&Sigma;}_{i=1}^k{a}_i=1,$ a_i＞0。

上述数学模型是MML-EM(最小信息长度准则与期望最大化算法)混合分布的一般公式，研究该模型的目的是精确地估计混合分布的分支和每个子分布的参数，如权重、均值和方差。详细的公式推导和算法实现可参阅Figueiredo和Jain发表于2002年的《Unsupervised Learning of Finite Mixture Models》一文。

综合上述两个实施案例，现提出多光谱遥感异常信息提取方法如图1所示。该方法将复杂多变的地质环境简化与“异常显著与否”等判断条件，并给出了相应的解决方案，因此具有一定的普适性与可操作性。如该装置增设了野外植被光谱测试这一项工作内容，以期应对多光谱影像上不显著的植被异常。

Claims (5)

1.一种多光谱遥感异常信息提取方法，包括以下步骤：

步骤1，遥感影像选择及预处理,得到预处理影像；

步骤2，对预处理影像进行掩膜处理；

步骤3，对掩膜处理后的遥感影像进行遥感蚀变信息增强处理；

其特征在于，还包括如下步骤：

步骤4，对遥感影像增强处理后的波段进行分形求和计算，得到异常增强影像，提取异常信息。

2.如权利要求1所述的一种多光谱遥感异常信息提取方法，其特征在于，当遥感影像中存在植物异常，还包括：

步骤5：对植被异常增强影像进行混合筛分处理，剔除假异常。

3.如权利要求1所述的一种多光谱遥感异常信息提取方法，其特征在于，步骤1中所述的预处理是指对选取的遥感影像依次进行辐射定标、大气校正、几何校正和正射校正处理。

4.如权利要求1所述的一种多光谱遥感异常信息提取方法，其特征在于，步骤4所述的对主成分转化后的影像波段进行分形求和计算是指：

步骤41，通过下式构建分形模型，得到N(r)-r散点图：

N(r)＝Cr^-D(r>0)

其中：r表示特征尺度，即影像DN值，(C>0)为比例系数,D>0为分维数，N(r)＝N(≥r)表示大于等于r的数目或和数；

步骤42，通过下式将N(r)-r散点图拟合为多段直线，每条直线对应于一个无标度区：

logN(r)＝-Dlog(r)+logC

步骤43，步骤42拟合的多段无标度区，应满足如下公式：

$\begin{array}{c}RSS=\underset{i=1}{\overset{i_0}{\mathrm{\&Sigma;}}}{\&lsqb;\log N\left(r_i\right)+D_1\log \left(r_i\right)-\log C_1\&rsqb;}^2+\underset{i=i_0+1}{\overset{i_1}{\mathrm{\&Sigma;}}}{\&lsqb;\log N\left(r_i\right)+D_2\log \left(r_i\right)-\log C_2\&rsqb;}^2\\ +\underset{i=i_1+1}{\overset{i_2}{\mathrm{\&Sigma;}}}{\&lsqb;\log N\left(r_i\right)+D_3\log \left(r_i\right)-\log C_3\&rsqb;}^2+...\\ +\underset{i=i_n+1}{\overset{N}{\mathrm{\&Sigma;}}}{\&lsqb;\log N\left(r_i\right)+D_{n+2}\log \left(r_i\right)-\log C_{n+2}\&rsqb;}^2\&RightArrow;Min\end{array}$

其中，D₁、D₂、D₃…D_n+2分别表示多段无标度区中第1个、第2个、第3个…第n+2个无标度区的分维数，C₁、C₂、C₃…C_n+2为第1个、第2个、第3个…第n+2个无标度区所对应的比例系数，其中表示从小到大排列的灰度值,i＝1,2,3…,n＝0,1,2…，D_n所对应的像元数记为N_n，RSS为拟合残差的平方；

其中，D_n和D_n+1之间的拐点记为T_n，该拐点T_n对应的灰度值即为异常灰度值

5.如权利要求2所述的一种多光谱遥感异常信息提取方法，其特征在于，步骤5中所述的对植被异常增强影像进行混合筛分处理是指：

步骤51，获取植被异常增强影像的灰度直方图；

步骤52，采用混合筛分模型将影像灰度直方图分割为若干个正态分布概率密度曲线，取不同正态分布概率密度曲线的交点为阈值，找到均值最小的正态分布或均值最大的正态分布。