一种基于NDVI和PanTex指数的高分辨率遥感影像新增建设用地图斑自动提取方法

摘要

本发明涉及一种基于归一化植被指数(NDVI)和建筑物存在指数(PanTex)的高分辨率遥感影像新增建设用地图斑自动提取方法，包括：输入前后时相高分辨率遥感影像，进行几何精校正和相对辐射校正；计算前后时相的NDVI和PanTex图像；对两期NDVI和PanTex图像分别进行非监督聚类；利用两期NDVI聚类图像，提取植被到建筑物、推填土的二值变化图像；利用两期PanTex聚类图像，提取植被、推填土到建筑物的二值变化图像；提取干扰地物区域；将提取的两幅变化图像并集运算，并将干扰地物掩膜去除，得到新增建设用地图像；分割后时相影像；计算每个分割图斑内的变化像素比例，提取的新增建设用地图斑。本发明能有效提取三类新增建设用地图斑，为土地利用变更调查提供辅助信息。

说明书

技术领域

本发明属于遥感影像处理技术领域，主要涉及土地利用变化检测中的新增建设用地检测方法，具体涉及一种基于NDVI和PanTex指数的高分辨率遥感影像新增建设用地图斑自动提取方法。>

背景技术

土地利用/覆盖情况是人类和土地相互作用的综合结果。同时，作为各种资源管理和地理信息服务所需要的最基本数据，土地利用/覆盖信息的获取、分析和更新，显得尤为重要。>

遥感影像数据以其宏观性、实时性，一直以来都是进行土地利用/覆盖及其变化检测的最重要手段。高分辨率遥感影像由于携带更多用于影像分析的细节信息，已成为土地利用/覆盖变化信息获取的重要数据源。然而一方面，影像数据量的增加和人工解译效率低成本高的矛盾更加突出；另一方面，由于高分辨率遥感影像本身具有的特点，使得传统的遥感影像分析和变化检测方法并不一定适用于高分辨率遥感数据的分析。可以说高分辨率遥感影像为土地利用/覆盖变化检测提供了更丰富信息的同时也带来了更多挑战。>

从技术层面来看，现有的土地利用/覆盖变化检测方法可以分为像素级和对象级两类。像素级变化检测方法通过对前后时相影像进行 逐像素分析来提取影像上发生变化的部分。这类方法对不同时相影像数据的几何校正和辐射校正精度要求较高，提取的变化像素易受噪声影响。对象级变化检测方法是指首先通过分割得到同质区域(对象)，然后提取对象特征，最后进行对象比较以提取变化对象。基于对象的变化检测方法以对象为基本单元，可以充分利用对象所固有的尺寸、形状和上下文信息。但基于对象的变化检测方法的难点在于对象的获取和寻找有利用变化检测的特征组合。常用的特征包括各种指数、形状和上下文信息、纹理，梯度等。在使用时，单独使用某一特征进行变化检测都可能造成漏检或误检，因此，应综合利用多特征实现变化检测。>

以影像实际发生的客观变化为依据，新增建设用地的检测目标可以分为以下三种类型：①前时相影像为植被覆盖，后时相有明显建筑痕迹；②前时相为植被覆盖，后时相有建设推填土痕迹；③前时相有推填土痕迹，后时相有明显建设痕迹。因此，应针对这三种类型选取合适的特征组合，以实现高分辨率遥感影像新增建设用地图斑的自动提取。>

发明内容

(一)发明目的>

本发明的目的是：针对高分辨率遥感影像新增建设用地的检测，提供一种基于光谱、纹理信息和分割的新增建设用地图斑自动提取方法，能够提取三种类型的新增建设用地图斑。>

(二)技术解决方案>

本发明提供了一种基于NDVI和PanTex指数的高分辨率遥感影像新增建设用地图斑自动提取方法，包含以下步骤：>

步骤1、对不同时相的高分辨率遥感影像进行几何精校正和相对辐射校正；>

步骤2、利用经过所述步骤1处理后的两期影像的多光谱波段，计算两期NDVI图像；利用经过所述步骤1处理后的两期影像的全色波段，计算两期PanTex图像；>

步骤3、对所述步骤2中得到的两期NDVI图像分别进行两类非监督聚类，聚类结果为植被和非植被；对所述步骤2中得到的两期PanTex图像分别进行两类非监督聚类，聚类结果为建筑物和非建筑物；>

步骤4、对所述步骤3中得到的NDVI聚类影像，提取耕地到推土、建筑物的变化图像；对所述步骤3中得到的PanTex聚类影像，提取耕地、推土到建筑物的变化图像；>

步骤5、利用经过所述步骤1处理后的两期影像，通过阈值法提取干扰地物，包括暗干扰地物，如水体、阴影，亮干扰地物，如云、噪声；>

步骤6、将所述步骤4中得到的两部分变化结果图像进行并集运算，并将由于干扰地物造成的伪变化区域掩膜去除，得到包含三种类 型的新增建设用地的二值结果图像；>

步骤7、对后时相影像进行分割，提取均值图斑对象；>

步骤8、将所述步骤7中得到的分割图斑和步骤6中得到的新增建设用地变化检测结果图像进行叠置，计算每个图斑内的变化像素比例，将比例超过设定阈值的图斑提取出来，作为新增建设用地变化图斑。>

所述步骤2中的NDVI能够很好的区分植被与非植被，其计算公式如下：>

$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{nir}}\text{-}{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{nir}}+{\rho }_{\mathrm{red}}}$

其中，ρ_nir和ρ_red分别代表某个像素近红波段和红波段的像素值。>

PanTex是一种基于结构信息的建筑物识别特征，能够提供简单、稳定的建筑物提取效果。PanTex定义了建筑物的识别模型，认为典型的建设用地是由建筑物与其周围的阴影组成的，它的对比度纹理测度在各个方向都高。纹理测度的定义基于灰度共生矩阵，灰度共生矩阵的计算如下：>

设f(x，y)为一幅二维数字图象，其大小为M×N，灰度级别为Z，则满足一定空间关系的灰度共生矩阵为：>

P(i，j)＝#{(x₁，y₁)，(x₂，y₂)∈M×N|f(x₁，y₁)＝i，f(x₂，y₂)＝j}>

式中#{x}表示集合x中的元素个数，P为Z×Z的矩阵，若(x₁，y₁)与(x₂，y₂)间的X方向偏移量为shiftX，Y方向偏移量为shiftY，则可以得到各个方向的灰度共生矩阵P(i，j|shiftX,shiftY)。>

进一步，在灰度共生矩阵的基础上计算对比度(Contrast)测度，计算公式如下：>

$\mathrm{Contrast}=\underset{i,j}{\Sigma }{(i-j)}^{2}P(i,j|\mathrm{shiftX},\mathrm{shiftY})$

PanTex的具体计算步骤如下：>

1.计算如下10个方向的对比度纹理测度，并对结果图像进行归一化处理；>

(shiftX，shiftY)＝{(1，-2)，(1，-1)，(2，-1)，(1，0)，(2，0)，(0，1)，(1，1)，(2，1)，(0，2)，(1，2)}>

2.将归一化处理后的10幅纹理测度图像进行逐像素取最小值，输出PanTex图像。>

所述步骤3中的两类非监督聚类采用K均值算法，得到的结果为一幅二值化图像。>

K均值聚类算法的步骤如下：>

1.选取K个类别的初始中心初始中心的选择对聚类结果有一定的影响，初始中心的选则一般有以下两种方法：>

a)任选K个类别的初始中心；>

b)将全部数据随机地分为K个类别，计算每类的重心，将这些重心作为K个类的初始中心。>

2.在第k次迭代中，对任一样本X按如下的方法把它划分到K个类别中的某一类。对于所有的i≠j，i=1，2，...，K，如果 $\left|\right|X-Z_j^{\left(k\right)}\left|\right|<\left|\right|X-Z_i^{\left(k\right)}\left|\right|,$则$X\in S_j^{\left(k\right)},$其中是以为中心的类。>

3.由第二步得到类新的类中心：>

$Z_j^{(k+1)}=\frac{1}{N_j}\underset{X\in S_j^{\left(k\right)}}{\Sigma }X$

式中，N_j为类中的样本数。>

4.对于所有的j＝1，2，...，K，如果则迭代结束，否则转到第二步继续进行迭代。>

对NDVI图像聚类时，认为NDVI值较高的一类为植被，赋值为1，值较低的一类为非植被，赋值为0，输出为一幅二值化图像；>

对PanTex图像聚类时，认为PanTex值较高的一类为建筑物，赋值为1，值较低的一类为非建筑物，赋值为0，输出为一幅二值化图像。>

所述步骤4中的变化图像的生成，其具体步骤是：>

利用前后时相NDVI二值聚类影像，将前时相取值为1且后时相取值为0的像素视为从耕地到推土、从耕地到建筑物的变化像素，赋值为1，其他像素赋值为0，输出一幅二值化图像。>

利用前后时相PanTex二值聚类影像，将前时相取值为0且后时相取值为1的像素视为从耕地到建筑物、从推土到建筑物的变化像素，赋值为1，其他像素赋值为0，输出一幅二值化图像。>

所述步骤5中的阈值法提取干扰地物，其具体步骤是：>

利用前后时相多光谱影像的红外波段，设定以下阈值，认为低于阈值的像素是变化检测的暗干扰地物，包括水体、阴影等：>

MIN(band_NIR)+MEAN(band_NIR)×λ>

其中λ为调整参数。将低于阈值的像素赋值为1，否则赋值为0，输出一幅二值化图像；>

利用前后时相多光谱影像的所有波段，设定以下阈值，认为所有 波段都高于阈值的像素是变化检测的亮干扰地物，包括厚云、噪声等：>

MIN(band_i)+MEAN(band_i)×λ_i

其中下标i表示波段号，λ_i为第i波段的调整参数。将低于阈值的像素赋值为1，否则赋值为0，输出一幅二值化图像。>

将上面两步的干扰地物提取结果图像并集运算，即只要在任一幅二值化图像上某个像素的取值为1，则将该像素赋值为1，否则赋值为0，输出一幅二值干扰地物提取图像。>

所述步骤6中的将两部分变化结果图像进行并集运算，并将结果图像中由于干扰地物造成的伪变化区域掩膜去除，其具体步骤是：>

1.令所述步骤4中得到的两幅二值化变化地物提取图像分别表示为I₁、I₂，令所述步骤5中得到的二值化干扰地物提取图像表示为I_noise；>

2.如果在任一幅图像上某像素的取值为1，且该像素在所述步骤5中得到的二值干扰地物提取图像上的取值为0，则将该像素视为新增建设用地并赋值为1，否则赋值为0，即：EITHER 1IF((I₁＝＝1||I₂＝＝1)&&I_noise＝＝0)OR 0 OTHERWISE输出结果为一幅二值化新增建设用地提取图像，包含了耕地到推土、耕地到建筑物、推土到建筑物这三种新增建设用地类型。>

所述步骤7中的图像分割，采用均值漂移(Mean Shift)分割算法，其具体步骤是：>

1.设置核函数G(x)，收敛阈值为ε，空间参数为h，最小像素数为mp；>

2.对于图像上的任一像素，它的特征向量由2维空间坐标值和p维光谱值组合而成，即x＝(x，y，b₁，b₂，...，b_p)，计算均值偏移迭代公式：>

$m_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\Sigma }}G\left(\frac{x_i-x}{h}\right)x_i}{\underset{i=1}{\overset{n}{\Sigma }}G\left(\frac{x_i-x}{h}\right)}$

3.计算均值漂移向量：>

M_h(x)＝m_h(x)-x>

如果||m_h(x)-x||<ε，则结束循环，否则令x＝m_h(x)，继续执行步骤2；>

4.重复执行步骤2～3，直到所有像素点都经过均值漂移处理；>

5.合并均质区域，并将像素个数低于mp的零碎区域合并到与其公共边界最长的邻近区域中去；>

6.将均质区域矢量化输出为图斑文件。>

所述步骤8中的图斑内变化像素比例是指：>

1.将所述步骤5中得到的二值化新增建设用地提取图像与所述步骤6中得到的图斑文件叠合，用逐个图斑掩膜图像。>

2.统计每个掩膜图像的总像素个数N和取值为1的像素的个数n_changcd(忽略背景值)，计算取值为1的像素占所有像素的比例即变化像素比例。>

设定阈值T，将变化像素比例高于阈值的图斑输出，作为检测到的新增建设用地图斑。>

(三)技术效果>

本发明与现有的技术方案相比具有如下的优点和有益效果：本发明将对新增建设用地提取效果较好的NDVI和PanTex特征进行结合采用，算法简单高效；利用Mean Shift分割算法产生均质图斑，并提出了变化像素比例的变化检测判别依据，消除了基于像素的变化检测算法的“椒盐”现象，且信息提取的结果为图斑形式，易于GIS数据库的更新；本发明不依赖于分类，无需样本学习，避免了因误分类造成的误差累积，不需要分类后处理。>

附图说明

图1是本发明实施例的基于NDVI和PanTex指数的高分辨率遥感影像新增建设用地图斑自动提取方法流程图；>

图2是前时相实验影像(北京市昌平区的ZY-1-02C影像，拍摄于2012年10月24日)；>

图3是后时相实验影像(北京市昌平区的ZY-1-02C影像，拍摄于2013年8月18日)；>

图4是后时相NDVI图像；>

图5是后时相PanTex图像；>

图6是二值化干扰地物提取结果图；>

图7是二值化新增建设用地提取结果图；>

图8是Mean Shift分割结果图；>

图9是提取的新增建设用地图斑。>

具体实施方式

本发明以两景北京市昌平区的ZY-1-02C影像为例，说明新增建设用地图斑提取的具体实施方式。前时相实验影像的拍摄时间为2012年10月24日，后时相实验影像的拍摄时间为2013年8月18日，实验影像如图2、3所示。下面结合附图对本发明进一步说明。>

如图1所示，是本发明实施例的基于NDVI和PanTex指数的高分辨率遥感影像新增建设用地图斑自动提取方法流程图，本实施例包括如下步骤：>

步骤1、对不同时相的高分辨率遥感影像进行几何精校正和相对辐射校正；>

步骤2、利用经过所述步骤1处理后的两期影像的多光谱波段，计算两期NDVI图像；利用经过所述步骤1处理后的两期影像的全色波段，计算两期PanTex图像；>

2.1：NDVI能够很好的区分植被与非植被，其计算公式如下：>

$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{nir}}\text{-}{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{nir}}+{\rho }_{\mathrm{red}}}$

其中，ρ_nir和ρ_red分别代表某个像素近红波段和红波段的像素值。图4是得到的后时相NDVI图像。>

2.2：PanTex是一种基于结构信息的建筑物识别特征，能够提供简单、稳定的建筑物提取效果。PanTex定义了建筑物的识别模型，认为典型的建设用地是由建筑物与其周围的阴影组成的，它的对比度纹理测度在各个方向都高。纹理测度的定义基于灰度共生矩阵，灰度共生矩阵的计算如下：>

设f(x，y)为一幅灰度级别为N的二维数字图象，设置窗口大小为M×M，则满足一定空间关系的灰度共生矩阵为：>

P(i，j)＝#{(x₁，y₁)，(x₂，y₂)∈M×M|f(x₁，y₁)＝i，f(x₂，y₂)＝j}>

式中#{x}表示集合x中的元素个数，P为N×N的矩阵，若(x₁，y₁)与(x₂，y₂)间的X方向偏移量为shiftX，Y方向偏移量为shiftY，则可以得到各个方向的灰度共生矩阵P(i，j|shiftX,shiftY)。PanTex建议的最优窗口大小与影像空间分辨率的乘积为50m，由于本发明实施例利用ZY-1-02C影像的HR波段进行PanTex计算，HR波段的分辨率为2.36m，因此推荐选取的窗口大小为M＝21。>

进一步，在灰度共生矩阵的基础上计算对比度(Contrast)测度，计算公式如下：>

$\mathrm{Contrast}=\underset{i,j}{\Sigma }{(i-j)}^{2}P(i,j|\mathrm{shiftX},\mathrm{shiftY})$

PanTex的具体计算步骤如下：>

1.计算如下10个方向的对比度纹理测度，并对结果图像进行归一化处理；>

(shiftX，shiftY)＝{(1，-2)，(1，-1)，(2，-1)，(1，0)，(2，0)，(0，1)，(1，1)，(2，1)，(0，2)，(1，2)}>

2.将归一化处理后的10幅纹理测度图像进行逐像素取最小值，输出PanTex图像。>

图5是得到的后时相PanTex图像。>

步骤3、对所述步骤2中得到的两期NDVI图像分别进行两类非监督聚类，聚类结果为植被和非植被；对所述步骤2中得到的两期PanTex图像分别进行两类非监督聚类，聚类结果为建筑物和非建筑物：>

3.1：非监督聚类采用K均值算法，得到的结果为一幅二值化图像。K均值聚类算法的步骤如下：>

1.选取K个类别的初始中心初始中心的选择对聚类结果有一定的影响，初始中心的选则一般有以下两种方法：>

c)任选K个类别的初始聚类中心；>

d)将全部数据随机地分为K个类别，计算每类的重心，将这些重心作为K个类的初始中心。>

本发明实施例需要进行两类聚类，因此选取K＝2，并任选两个类别的初始聚类中心。>

2.在第k次迭代中，对任一样本X按如下的方法把它划分到K个类别中的某一类。对于所有的i≠j，i＝1，2，...，K，如果 $\left|\right|X-Z_j^{\left(k\right)}\left|\right|<\left|\right|X-Z_i^{\left(k\right)}\left|\right|,$则$X\in S_j^{\left(k\right)},$其中是以为中心的类。>

3.由第二步得到类新的类中心：>

$Z_j^{(k+1)}=\frac{1}{N_j}\underset{X\in S_j^{\left(k\right)}}{\Sigma }X$

式中，N_j为类中的样本数。>

4.对于所有的j＝1，2，...，K，如果则迭代结束，否则转到第二步继续进行迭代。>

3.2：对NDVI图像进行聚类时，认为NDVI值较高的一类为植被，赋值为1，值较低的一类为非植被，赋值为0，输出为一幅二值化图像；>

3.3：对PanTex图像进行聚类时，认为PanTex值较高的一类为建筑物，赋值为1，值较低的一类为非建筑物，赋值为0，输出为一幅二值化图像。>

步骤4、对所述步骤3中得到的NDVI聚类影像，提取耕地到推土、建筑物的变化图像；对所述步骤3中得到的PanTex聚类影像，提取耕地、推土到建筑物的变化图像；>

4.1：利用前后时相NDVI二值聚类影像，将前时相取值为1且后时相取值为0的像素，视为从耕地到推土、从耕地到建筑物的变化像素，赋值为1，否则赋值为0，输出一幅二值化图像。>

4.2：利用前后时相PanTex二值聚类影像，将前时相取值为0且后时相取值为1的像素视为从耕地到建筑物、从推土到建筑物的变化像素，赋值为1，否则赋值为0，输出一幅二值化图像。>

步骤5、利用经过所述步骤1处理后的两期影像，通过阈值法提取干扰地物，包括暗干扰地物，如水体、阴影，亮干扰地物，如云、噪声；>

5.1：利用前后时相多光谱影像的红外波段，设定以下闽值，认为低于阈值的像素是变化检测的暗干扰地物，包括水体、阴影等：>

MIN(band_NIR)+MEAN(band_NIR)×λ>

其中λ为调整参数，本发明实施例令λ＝0.2。将低于阈值的像素赋值为1，否则赋值为0，输出一幅二值化图像；>

5.2：利用前后时相多光谱影像的所有波段，设定以下阈值，认为所有波段都高于阈值的像素是变化检测的亮干扰地物，包括厚云、噪声等：>

MAX(band_i)-MEAN(band_i)×λ_i

其中下标i表示波段号，λ_i为第i波段的调整参数，本发明实施例令λ₁＝λ₂＝λ₃＝0.2。将低于阈值的像素赋值为1，否则赋值为0，输出一幅二值化图像。>

5.3：将上面两步的干扰地物提取结果图像并集运算，即只要在任一幅二值化图像上某个像素的取值为1，则将该像素赋值为1，否则赋值为0，输出一幅二值干扰地物提取图像。图6是后时相多光谱影像与二值干扰地物提取图像叠置显示的效果，其中绿色像素为提取的干扰地物。>

步骤6、将所述步骤4中得到的两部分变化结果图像进行并集运算，并将由于干扰地物造成的伪变化区域掩膜去除，得到包含三种类型的新增建设用地的二值结果图像；>

6.1：输入所述步骤4中得到的两幅二值变化地物提取图像，将其分别表示为I₁、I₂；同时，输入所述步骤5中得到的二值干扰地物提取图像，将其表示为I_noise；>

6.2：如果在任一幅图像上某像素的取值为1，且该像素在所述步骤5中得到的二值干扰地物提取图像上的取值为0，则将该像素视为新增建设用地并赋值为1，否则赋值为0，即：>

EITHER 1 IF((I₁＝＝1||I₂＝＝1)&&I_noise＝＝0)OR 0 OTHERWISE>

输出结果为一幅二值化新增建设用地提取图像，包含了耕地到推土、耕地到建筑物、推土到建筑物这三种新增建设用地类型。图7是得到的二值化新增建设用地提取图像。>

步骤7、对后时相影像进行分割，提取均值图斑对象；>

7.1.1：将后时相多光谱影像作为待分割影像，采用均值漂移(Mean Shift)分割算法。>

7.1.2：设置核函数G(x)，收敛阈值为ε，最小像素数为mp；>

本发明实施例选取的核函数由两部分组成，表示如下：>

$G_{h_s,h_r}=\frac{C}{h_s^2{h}_r^p}g\left({\left|\right|\frac{x^s}{h_s}\left|\right|}^2\right)g\left({\left|\right|\frac{x^r}{h_r}\left|\right|}^2\right)$

其中选取的颜色参数为h_s＝6.5、空间参数为h_r＝8，收敛阈值为ε＝1e^-6，最小像素数mp＝100。>

7.1.3：对于图像上的任一像素，它的特征向量由2维空间坐标值和p维光谱值组合而成，即x＝(x，y，b₁，b₂，...，b_p)，计算均值偏移迭代公式：>

$m_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\Sigma }}G\left(\frac{x_i-x}{h}\right)x_i}{\underset{i=1}{\overset{n}{\Sigma }}G\left(\frac{x_i-x}{h}\right)}$

本发明实施例的分割影像为后时相的ZY-1-02C多光谱影像，波段数为3，即p＝3。>

7.1.3：计算均值漂移向量：>

M_h(x)＝m_h(x)-x>

如果||m_h(x)-x||<ε，则结束循环，否则令x＝m_h(x)，继续执行步骤2；>

7.1.4：重复执行步骤7.1.2～7.1.3，直到所有像素点都经过均值漂移处理；>

7.1.5：合并均质区域，并将像素个数低于mp的零碎区域合并到与其公共边界最长的邻近区域中去；>

7.1.6：将均质区域矢量化输出为图斑文件。>

图8为后时相多光谱影像经过Mean Shift分割得到的图斑与后时相影像叠置显示的效果。>

步骤8、将所述步骤7中得到的分割图斑和步骤6中得到的新增建设用地变化检测结果图像进行叠置，计算每个图斑内的变化像素比 例，将比例超过设定阈值的图斑提取出来，作为新增建设用地变化图斑。>

8.1：将所述步骤5中得到的二值化新增建设用地提取图像与所述步骤6中得到的图斑文件叠合，用逐个图斑掩膜图像。>

8.2：统计每个掩膜图像的总像素个数N和取值为1的像素的个数n_changed(忽略背景值)，计算取值为1的像素占所有像素的比例即变化像素比例。>

8.3：设定阈值T，将变化像素比例高于阈值的图斑输出，作为检测到的新增建设用地图斑。>

本发明实施例选取的变化像素比例阈值为T＝0.6。>

图9为本算法提取的新增建设用地图斑。>

实验结果表明，通过本技术方案，对新增建设用地图斑的提取能够取得较理想的结果。相比现有算法，本发明在保证检测率的同时，避免了人工选取训练样本的过程，实现了新增建设用地的自动提取。>

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