线性光谱解混的城市不透水面遥感提取方法与流程

本发明涉及遥感图像处理领域，更具体地说，涉及一种线性光谱解混的城市不透水面遥感提取方法。

背景技术：

采用遥感数据提取城市不透水面过程中，首先选取合适的数据源，经过几何校正和大气校正等预处理，对遥感影像上地物的分类和提取的方法有线性光谱解混方法、监督分类、非监督分类、决策树分类、人工神经网络等，常用的提取方法是线性光谱解混方法(LSMA)。在遥感提取城市不透水面过程中，传统的线性光谱解混方法存在两个问题：一、在森林区域，有少量不透水面如个别房屋等高反照率地物，在端元选取时容易被忽略，影响端元选取的准确度；二、有些区域像元反照率比较低，如植被、裸土等，对不透水面中低反照率部分的影响较大，容易出现大面积解混异常区域，一般线性光谱解混提取不透水面时，需要对提取结果采用阈值法去除异常区域。阈值的选择直接影响不透水面的提取精度，且具有经验性，人为影响因素较大。

技术实现要素：

有鉴于此，本发明目的是提供本发明提出一种线性光谱解混方法来提取城市不透水面。首先利用纯净像元指数结合手动端元选取的方式优化端元选取，确保森林区域少量高反照率地物的端元被合理选取；然后利用归一化建筑指数(Normalized Difference Building Index，NDBI)对低反照率结果进行掩膜，剔除不透水面中低反照率地物解混异常区域。

为了解决上述技术问题，本发明的技术方案是：一种线性光谱解混的城市不透水面遥感提取方法，

步骤一，获取遥感影像素材文件并进行数据预处理与水体掩摸；

步骤二，将遥感影像素材文件通过最小噪声分离变换得到数据文件一；

步骤三，将数据文件一通过纯净像元指数计算得出在数据文件一中每一像元的纯净度；

步骤四，根据纯净度预设值(阈值)，在数据文件一中筛选出大于该值的像元，即数据文件二；

步骤五，利用数据文件二并根据筛选结果在分量特征散点二维模型构建拟化分量三角形，根据三角形的三个顶角位置筛选端元；

步骤六，以步骤五中的端元为依据，结合线性光谱解混模型采用最小二乘法求解每个像元中各端元所占的比例，其数学模型为，错误！未找到引用源。；错误！未找到引用源。

式中i＝1,2…,M，M为光谱波段数；N为端元数目。R_i是波段i的反射率；f_k是端元k在像元中所占的面积比例；R_ik波段i上端元k的光谱反射率；ER_i是波段i的残差；

步骤七，对解混结果中的低反照率地物采用归一化建筑物指数进行掩膜处理。

进一步地：所述的端元类型设置为四个且为高反照率、低反照率、植被和土壤端元。

进一步地：步骤六，结合线性光谱解混模型采用最小二乘法求解每个像元中各端元所占的比例。

进一步地：在步骤七中，结合空间分辨率为2m的Google Earth高分影像，采用均方根误差RMSE，系统误差SE，平均绝对误差MAE以及精度验证趋势图对不透水面提取精度进行验证：错误！未找到引用源。错误！未找到引用源。；错误！未找到引用源。错误！未找到引用源。错误！未找到引用源。；错误！未找到引用源。错误！未找到引用源。错误！未找到引用源。；

其中，为本研究中从Landsat 8OLI数据中提取的样本区域内不透水面的比例；X_i为高分验证影像样本区域内不透水面的比例；N为样本数。

进一步地：步骤七中，根据Landsat 8第5波段和第6波段反射率计算遥感影像的归一化建筑物指数NDBI：其中，MIR为遥感影像中红外波段反射率，NIR为遥感影像近红外波段反射率。

进一步地，步骤四中，构建拟化分量三角形的方法如下，根据筛选结果的坐标值寻找极限位置的四个极点，分别对四个极点中的三个极点进行连线围成四个对照分量三角形，求取每一对照分量三角形中像元的数量，将像元数量最多的对照分量三角形作为拟化分量三角形。

进一步地，步骤四中，筛选端元的方法如下，以拟化分量三角形的三个顶点为中心分别做圆，圆的半径为R，当前纯净度为X，线性减小变量R或增加变量X的值，直至每个顶点剩余一个像元，该像元即为端元。

进一步地，步骤四中，定义线性率Q1和Q2，每次变化时，使R＝Q1*R或X＝Q2*X。

进一步地，步骤五中，结合高分影像，通过目视、比较分析的方法选取纯净的像元。

本发明技术效果主要体现在以下方面：本研究提出结合手动选取端元，即从MNF变换结果中将PPI阈值筛选出来的像元掩膜，作为手动选取端元的基础，像元数量减少，纯度提高，利用二维散点图，以高分影像为参考，继续选取影像空间顶点区域零散的点，确保植被和高反照率地物的端元具有各自代表性，使高反照率地物在大面积森林中可以被识别区分，提高端元选取的精度。

为了减少人为选择阈值对不透水面提取精度的影响，利用归一化建筑指数对不透水面的敏感性以及对植被土壤等地物的区分性，提出利用归一化建筑指数对低反照率结果进行处理，剔除不透水面的异常区域，大大提高不透水面的提取精度。

附图说明

图1：本发明的线性光谱解混流程图；

图2：传统的LSMA端元选取结果图；

图3：本发明的LSMA端元选取结果图；

图4：线性光谱解混结果图；

图5：传统LSMA结果与结合手动端元选取的LSMA结果的差异图；

图6：研究区域100个样本的空间分布；

图7：精度评价结果图；

图8：MNF变换结果图；

图9：PPI计算结果图；

图10：经过PPI阈值筛选的MNF变换结果图。

具体实施方式

以下结合附图，对本发明的具体实施方式作进一步详述，以使本发明技术方案更易于理解和掌握。

传统线性光谱解混模型提取不透水面流程包括：

(1)最小噪声分离变换(Minimum Noise Fraction Transform，MNF)——去除影像噪声，将遥感影像主要信息集中在前几个波段；

(2)纯净像元指数(Pixel Purity Index，PPI)计算——计算每个像元的纯净度，用于选择纯净度比较高的像元；

(3)端元选取——选择纯净度较高，并且光谱特征具有代表性的像元，用于线性光谱解混方程中的输入参数；

(4)线性光谱模型解混——假设影像中每个像元的反射率由该像元中所有地物端元的反射率以及其所占像元面积比例为权重系数的线性组合，其数学模型为：

错误！未找到引用源。，式中i＝1,2…,M，M为光谱波段数；N为端元数目。R_i是波段i的反射率；f_k是端元k在像元中所占的面积比例；R_ik波段i上端元k的光谱反射率；ER_i是波段i的残差；

(5)精度验证。在遥感影像上选择土壤、植被、高反照率和低反照率四个端元作为解混对象，采用最小二乘法求解各端元所占的比例。不透水面采用高反照率和低反照率组合表示。

本发明做出了以下改进，(1)在传统线性光谱解混模型基础上，在端元选取阶段，结合手动方式选取合适的端元，即从MNF变换结果中计算PPI，选取合适阈值筛选纯净度较高的像元，作为手动选取端元的基础，增加人为因素；分量两两之间的像元特征空间一般呈现三角形，端元分布在特征空间的三个顶点上。通常认为，纯度高的端元靠近边缘。因此，选取端元时尽量靠近三角形顶点区域，确保适用于研究区域。(2)在遥感数据OLI前四个波段，低反照率地物与植被、土壤等透水面的区分不明显，造成解混过程中的错分现象。在5、6波段，它们之间的光谱差异明显增大。计算遥感影像的NDBI值：错误！未找到引用源。，其中，MIR为遥感影像的中红外波段，即第6波段，NIR为遥感影像的近红外波段，即第5波段。然后通过与高分影像对比，确定NDBI阈值。将符合阈值条件的像元作为掩膜，对解混后不透水面低反照率结果进行裁剪，去除掉低于NDBI值的部分。此外，精度验证阶段，结合空间分辨率为2m的Google Earth高分影像，采用均方根误差RMSE，系统误差SE，平均绝对误差MAE以及精度验证趋势图对不透水面提取精度进行验证：错误！未找到引用源。；错误！未找到引用源。；错误！未找到引用源。；其中，为本研究中从Landsat 8OLI数据中提取的样本区域内不透水面的比例；X_i为高分验证影像样本区域内不透水面的比例；N为样本数。

从图2和图3可以看出经过筛选后的纯净像元数量大大减少，纯度普遍较高。通过结合手动选取端元，可以有目的性的将这种小面积不透水面选取出来，提高植被端元和不透水面高反照率端元的准确性，如图3(d)所示，高反照率端元的反射率高于植被端元的反射率。图4(b)为结合手动端元选取方式的LSMA模型提取的不透水面。由图4(a)和4(b)可以看出，结合手动选取端元的LSMA模型比传统的LSMA模型能更好地改善不透水面的提取精度，特别是在森林区域，如图5所示。图5为传统LSMA结果与手动端元选取优化LSMA结果的差异图，可以看出结合手动端元选取的方法能较好地减少大面积森林区植被的影响，而且对森林区域不透水面的提取也有较好的效果。这是因为经过手动选取端元后，植被与高反照率地物被区分开了，从图3(d)的光谱特性曲线也可以看出。

尽管结合手动端元选取方法可以较好地提高不透水面的提取精度，但是在森林区域仍存在较大误差，如图4(b)所示。为了进一步改善不透水面提取精度，本研究采用归一化建筑指数对低反照率地物进行优化处理。首先计算影像的NDBI值；统计发现计算的NDBI值在-0.66—0.45之间。通过与高分影像的对比发现，绝大部分不透水面的NDBI值在-0.15—0.45之间，因此设置其阈值为-0.15作为不透水面与透水面的临界值，将NDBI值小于-0.15的像元作为掩膜区域，对不透水面低反照率提取结果进行处理，去除异常部分，结果如图4(c)所示。通过对比图4(a)和4(c)可知，利用第5、6波段计算的NDBI值可以有效地突出不透水面，较好地区分较暗区域的透水面与不透水面，将森林区的透水面掩膜掉，更好地提取森林区的建筑物，这是因为NDBI对建筑区域具有一定的敏感性，较容易地从植被、土壤中将建筑区分出来。因此，将NDBI值对不透水面低反照率结果进行掩膜可以有效地提高不透水面的提取精度。

如图6，每个样本区域面积为230400平方米，对样本区域的高分影像进行人工解译以及矢量化不透水面并进行实地验证，并计算每个样本区域内不透水面的比例。

将不同方法提取的不透水面分别与验证数据进行比较。传统LSMA方法得到不透水面的RMS平均值为0.0093，RMSE为0.306，相关系数R为0.898，SE为0.21，MAE为0.24。而结合手动选取端元改进后的LSMA方法得到不透水面的RMS平均值为0.0079，从频率分布可以看出，绝大部分像元的RMS值都低于0.005，RMSE为0.293，相关系数R为0.932，SE为0.20，MAE为0.23。相对传统的LSMA方法，结合手动端元选取的LSMA方法可以较好地减少不透水面提取的异常值，如图7(a)和7(b)所示。然而，结合手动端元选取的LSMA方法得到的结果与验证数据相比仍存在较大的误差，很多透水面的地物(植被、土壤)被误分为不透水面了，如7(b)所示，这是不合理的。通过NDBI筛选去除异常部分后，得到不透水面的精度有了较大的提高，如图7(c)所示，RMSE减少到0.125，相关系数R高达0.943，SE减少为-0.035。总的来说，NDBI指数对可以较好地将解混结果中绝大部分异常区域剔除，如图4(c)所示，异常区域主要集中在大面积森林区域。这是因为NDBI指数可以较好地区分不透水面中低反照率地物和较暗的像元地物，对不透水面提取结果起到较好的优化作用。

选取合适的遥感数据源，采用2015年10月18日Landsat 8OLI数据，分别采用传统的线性光谱解混光谱解混方法和改进的线性光谱解混光谱解混方法提取不透水面，利用高分影像对不同方法提取的不透水面精度进行验证。整个过程可以在软件ENVI和Arcgis中实现：

在ENVI中，将遥感影像进行几何校正和大气校正。采用Bandmath计算影像的归一化水体指数，分析后选取0.4作为阈值，利用MASK功能去除水体。

MNF变换：选择Spectral/MNF Rotation/Forward MNF/Estimate Noise Statistics from Data。如下图8。

PPI计算：Spectral/Pixel Purity Index/Pixel Purity Index(PPI)New Output Band利用MNF变换结果图进行PPI计算。如下图9。

手动端元选取：使用Band Threshold to ROI Input Band工具，根据PPI计算结果选取合适的PPI阈值5，经过MASK筛选MNF变换结果图，如下图10。选取纯度较高的像元作为手动选取端元对象。通过目视、比较分析的方法，利用二维散点图，以高分影像为参考，继续选取影像空间顶点区域零散的点，选取植被、土壤、高反照率地物、低反照率地物四类。将选取的端元保存至ASCLL文件，将其光谱曲线图保存。

线性光谱解混：采用ENVI扩展工具——FCLS Spectral Unmixing工具，采用上述过程中保存的光谱特征曲线图解混遥感影像。得到的结果中包括植被、土壤、高反照率地物、低反照率地物波段，以及一个误差RMS波段：NDBI掩膜：计算全图NDBI，统计发现计算的NDBI值在-0.66—0.45之间。通过与高分影像的对比发现，绝大部分不透水面的NDBI值在-0.15—0.45之间，因此设置其阈值为-0.15作为不透水面与透水面的临界值，将NDBI值小于-0.15的像元作为掩膜区域，对不透水面低反照率提取结果进行处理。

精度验证：采用经过几何校正后的同时相Google Earth高分影像作为不透水面提取精度的验证影像。随机选取100个空间分布均匀的样本区域，每个样本区域面积为230400平方米，分别对NDBI掩膜结果和高分影像进行掩膜，得到100个验证区。在Arcgis中，对样本区域的高分影像进行人工解译以及矢量化不透水面并进行实地验证，并计算每个样本区域内不透水面的比例。在Excel表中分别列出解混结果验证区和高分矢量化验证区，求出SE、MAE、以及RMSE，并且绘制趋势图，同时添加趋势线，显示相关系数。

当然，以上只是本发明的典型实例，除此之外，本发明还可以有其它多种具体实施方式，凡采用等同替换或等效变换形成的技术方案，均落在本发明要求保护的范围之内。