一种基于航空影像的水体与居民地变化检测方法

摘要

本发明公开了一种基于航空影像水体与居民地的变化检测方法，包括以下步骤：(1)近期航空影像预处理；(2)旧矢量数据栅格化；(3)影像和栅格化的矢量数据分别生成同级金字塔影像；(4)计算影像的边缘特征和LAW能量特征；(5)利用已有旧矢量数据采集航空影像居民地的训练样本，获得居民地的特征值，根据特征值确定阈值，在航空影像上提取居民地；(6)利用已有旧矢量数据采集航空影像水系的训练样本，结合居民地提取结果，获得水系得特征值，根据特征值确定阈值，在航空影像上提取水系；(7)金字塔还原，将提取结果还原为原始分辨率；(8)将提取结果与旧矢量数据叠加，获得变化区域和变化类型，提示给用户。采用该方法能有效地利用旧矢量数据，快速检测变化区域，减少利用航空影像进行地图修测的工作量。

说明书

技术领域

本发明涉及一种基于航空影像的水体与居民地变化检测方法，属于影像处理领域，特别是属于航空影像处理与目标提取技术领域。

背景技术

自然界的变化和人类的活动改变了地表景观和土地利用形式。及时、有效地监测变化区域和变化类型，是城市变化检测和地图自动修测的基础工作。现代对地观测技术进入了一个能够及时、快速地提供各种对地观测数据的阶段。利用近期航空影像，实现数据快速更新，成为地图修测的一项重要手段。自动、高效、准确地从某地区近期航空影像中提取水体和居民地并识别其变化区域和变化类型，可大大减少地图修测的人工作业量，对地图修测有重要的意义。

利用影像实现水体和居民地的变化检测首要解决的问题是水体和居民地的提取。Andrea采用区域增长的方法从中小比例尺影像提取居民地；Heipke将纹理分析方法和建筑物自动检测算法结合提取居民地，并对GIS数据库进行半自动修测；Paola深入研究居民地提取算法，比较BP神经网络和最大似然分类两种方法提取居民地，分析算法的局限性和适用性；杨存建从对居民地的遥感信息机理入手，分析目标在TM影像各个波段上与其他地物的可分性；易祎以人工选取种子点的方法获得训练样区，采用区域增长方法提取目标；佘琼基于LBP/C(LocalBinary Pattern/Contrast)法提取居民地；朱莉娜将水体分为主干河流、湖泊和支流水系网，分别对主干河流采用模式分类，对湖泊采用区域增长，对支流水系网采用多模板匹配，从遥感影像上提取水体。由于水体和居民地的复杂性，上述方法存在两个缺陷。一是这类提取方法均需要人工选择大量的训练样本；二是在目标的提取中，需要人工给出提取阈值。因此，大图幅海量数据的更新需要大量的人工干预，无法适应生产需求。

发明内容

为了克服现有技术的不足，本发明提供了一种基于航空影像的水体与居民地变化检测方法，充分利用已有旧矢量数据，在无人工干预的条件下实现样本训练和阈值确定，并根据不同的置信度和变化类型将变化检测结果提供给作业员，提高了地图修测的自动化程度，缩短了数据更新周期。

实现本发明目的采用的技术方案是：一种基于航空影像的水体与居民地变化检测方法，包括以下步骤：

(1)对近期航空影像进行预处理，预处理操作包括：(a)对航空影像的正射纠正，(b)对正射影像的拼接，(c)根据同区域旧矢量数据范围，对拼接的正射影像进行裁减，(d)将裁减结果作为变化检测的输入影像；

(2)将待测区域水体和居民地的旧矢量数据栅格化为二值图像，待测区域内部与外部赋不同的标识；

(3)创建正射影像和二值图像的金字塔影像：用低通滤波器对正射影像和二值图像进行低通滤波，从而得到金字塔影像；

(4)在正射影像的金字塔影像上计算纹理能量特征，包括以下步骤：

I.利用微窗口对原始图像进行卷积运算，获得卷积图像，通过式$F(i,j)=\underset{k=0}{\overset{N-1}{\Sigma }}\underset{l=0}{\overset{N-1}{\Sigma }}A(k,l)I(k+i,l+j)$进行卷积运算，其中，A(k，l)为能量模板，I(k+i，l+j)为原始图像，F(i，j)为卷积结果图像，

II.对于所得卷积图像，通过式$\mathrm{TE}(i,j)=\frac{i}{{(2n+1)}^2}\underset{k=i-n}{\overset{i+n}{\Sigma }}\underset{l=j-n}{\overset{j+n}{\Sigma }}F(k,l)$求出各像点在以该像点为中心的宏窗口中的均值，并将所求得的均值作为该点的能量值，其中TE(i，j)是能量值，F(k，l)为步骤I所得的卷积结果图像，n为宏窗口大小，

III.采用劳(LAW)能量的三个能量模板分别计算每个像素的能量值，得到对应的三幅特征图像；

(5)利用待测区域的旧矢量数据采集航空影像居民地的训练样本，获得居民地样本的特征值，根据特征值确定阈值，在特征图像上提取居民地，提取居民地的步骤如下：

(5-a)将所有特征图像均处理，转入步骤(5-e)，处理一幅特征图像，则将待测区域的旧矢量数据二值化图中标记为居民地的像素作为训练样本，统计训练样本在该特征图像上的特征值；

(5-b)以特征值为横轴，特征值对应的像素个数为纵轴，建立训练样本特征值直方图；

(5-c)计算直方图参数，包括直方图的峰值、均值、方差、谷等，根据直方图参数确定不同置信度下的提取阈值；

(5-d)根据确定的阈值，提取居民地并标记，得到与当前特征图像对应的结果图像；

(5-e)逐点标记结果图像，若该点像素在所有特征图像上都被标记为居民地，则该点标记为居民地，否则该点标记为背景点；

(6)在正射影像上计算边缘密度，获得提取水系的特征图像，利用旧矢量数据采集航空影像水系的训练样本，结合居民地提取结果，获得水系样本的特征值，根据特征值确定阈值，在特征图像上提取水系，提取水系的步骤如下：

(6-a)根据待测区域的旧矢量数据和高置信度下提取居民地结果，分析水系变化率，若原水系区域变为居民地，则原水体消失，统计原水系变为居民地的区域，除以原水系的总体面积，得到水系变化率；

(6-b)根据水系变化率确定水系变化检测，若水系变化率大于给定阈值，则进入步骤(6-c)，否则进入步骤(6-f)；

(6-c)由高置信度居民地提取结果在水系矢量区域内采集训练样本，统计居民地样本的直方图，计算直方图参数，包括峰值、方差、均值、谷等；

(6-d)获得判断水系转变为居民地区域的阈值，即得到水系消失阈值；

(6-e)在边缘密度图像上逐像素判断，若该像素边缘密度大于水系消失阈值，则标定为消失，进入步骤(6-j)；

(6-f)利用水系矢量栅格化的二值图像在正射影像上采集水系训练样本；

(6-g)统计训练样本的特征直方图，并计算直方图参数，包括峰值、方差、均值、谷等；

(6-h)根据直方图参数计算不同置信度下水系提取的阈值；

(6-i)在不同置信度下，通过阈值和边缘密度图像提取水系；

(6-j)输出不同置信度下水系提取结果，或水系消失区域提取结果；

(7)将居民地和水系的提取结果与旧矢量数据叠加，获得变化区域和变化类型，步骤如下：

(7-a)在水系变化检测中，变化阈值大于给定阈值，则直接转入(7-e)；

(7-b)将所有提取结果图均处理，转入步骤(7-f)，否则，取某一置信度提取结果与旧矢量栅格图叠加，判断提取结果图中被标记为居民地或水系的像素，若该像素在旧矢量栅格图上被标记为背景点，则将该点标记为新增；否则，将该点标记为背景，依次判断各像素，得到二值图像；

(7-c)对二值图像进行形态学处理，消除矢量编辑中多边形趋近而引起的伪变化；

(7-d)逐点判断旧矢量栅格图上被标记为居民地或水系的像素，若该像素在提取结果图上被标记为背景点，则将该点标记为消失；否则标记为背景；

(7-e)采用步骤(7-c)操作消除伪变化区域；

(7-f)根据不同的标记值，赋予结果图像各像素不同的颜色值，将结果图提供给用户。

其中在步骤(6)中，通过以下步骤在正射影像上计算边缘密度特征：

I.用Sobel算子与图像进行卷积，生成Sobel边缘图像；

II.计算Sobel边缘图像的边缘密度，得到特征图像，即

$F(i,j)=\frac{1}{N\times N}\underset{k=-\frac{N}{2}}{\overset{\frac{N}{2}}{\Sigma }}\underset{l=-\frac{N}{2}}{\overset{\frac{N}{s}}{\Sigma }}S(k+i,l+j)$

其中，S(i，j)为Sobel边缘图像；F(i，j)为边缘密度图像；N为边缘密度窗口大小。

对边界的处理采用镜像对称的方法，即：

行向：S(-i，j)＝S(i，j)

S(i+k，j)＝S(i-k，j)，其中i＜＝N，k＝1，2，3，…，N为总行数

纵向：S(i，-j)＝S(i，j)

S(i，j+k)＝S(i，j-k)，其中i＜＝N，k＝1，2，3，…，N为总列数

本发明充分利用待修测地区旧矢量数据对航空影像进行上述操作，可以快速有效地检测出水体与居民地的变化，并根据不同的置信度和不同的变化类型，将对应的不同检测结果提供给用户，减少了地图修测人员的工作量，提高了地图修测的自动化程度，缩短了数据更新周期。

附图说明

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

图1为本发明流程示意图。

图2是武汉地区1m空间分辨率的航空影像。

图3是武汉地区居民地矢量数据栅格化生成的1m空间分辨率的二值图像。

图4是武汉地区水系矢量数据栅格化生成的1m空间分辨率的二值图像。

图5是武汉地区两种置信度下居民地提取结果。

图6是武汉地区两种置信度下水系提取结果。

图7是武汉地区居民地变化检测结果。

图8是武汉地区水系变化检测结果。

图9是北京地区1m空间分辨率的航空影像。

图10是北京地区消失水系提取结果。

图11是北京地区水系变化检测结果。

具体实施方式

本实施例主要通过处理栅格化的旧矢量数据和经过预处理的航空影像来提取居民地，从而提取水系，最后居民地提取结果和水体提取结果分别与对应的旧矢量数据叠加，获得变化区域和变化类型，实现变化检测。其流程见图1。

本实施例选取武汉地区航空影像(1133×1079)如图2和北京地区航空影像(2049×1565)如图9，针对不同情况进行两组实验，具体实施步骤如下：

(1)分别输入武汉地区和北京地区航空影像，对航空影像预处理：

(a)航空影像的正射纠正；

(b)正射影像的拼接；

(c)根据同区域旧矢量数据范围，对拼接的正射影像进行裁减；

(2)将武汉地区和北京地区旧矢量数据分别栅格化为二值图像。武汉地区居民地矢量栅格化结果如图3，水系矢量栅格化结果如图4；

(3)居民地变化检测的具体实施步骤如下：

(3-a)采用高斯低通滤波器对栅格化的二值图像(图3)和预处理后的影像分别生成相同分辨率金字塔影像；

(3-b)利用栅格化的二值图像在航空影像的金字塔影像上采集居民地的训练样本；

(3-c)在航空影像的金字塔影像上计算纹理能量特征，具体实施步骤如下：

I.利用微窗口对原始图像进行卷积运算，通过式$F(i,j)=\underset{k=0}{\overset{N-1}{\Sigma }}\underset{l=0}{\overset{N-1}{\Sigma }}A(k,l)I(k+i,l+j)$获得卷积图像，其中，A(k，l)代表能量模板，I(k+i，l+j)代表原始图像即为图2和图9，F(i，j)代表卷积结果图像。

II.对于所得卷积图像，通过式$\mathrm{TE}(i,j)=\frac{1}{{(2n+1)}^2}\underset{k=i-n}{\overset{i+n}{\Sigma }}\underset{l=-j-n}{\overset{j+n}{\Sigma }}F(k,l)$求出各像点在以该像点为中心的宏窗口中的均值，并将所求得的均值作为该点的能量值，其中，TE(i，j)是能量值，F(k，l)代表步骤I所得的卷积结果图像，n为宏窗口大小。

III.采用劳(LAW)能量的三个能量模板分别计算每个像素的能量值，得到对应的三幅特征图像。

(3-d)统计训练样本的特征直方图，并计算直方图参数，包括峰值、方差、均值、谷，具体计算方法如下：以训练样本的能量特征值为横轴，能量特征值对应的像素个数为纵轴建立样本特征直方图H，依据下式计算直方图峰值、方差、均值、谷：

peak＝Max{X|(X，Y)∈H}

vally＝Min{X|(X，Y)∈H}

$N=\underset{i=1}{\overset{n}{\Sigma }}{Y}_i$

$\mathrm{mean}=\underset{i=1}{\overset{n}{\Sigma }}{X}_i\times \left(\frac{Y_i}{N}\right)$

$\mathrm{variance}=\sqrt{\underset{i=1}{\overset{n}{\Sigma }}{(X_i-\mathrm{mean})}^2\times \left(\frac{Y_i}{N}\right)}$

其中，peak代表样本特征直方图H的峰值，vally代表样本特征直方图的H谷值，N代表样本特征直方图像素总数，mean代表样本特征直方图H的均值，varia代表样本特征直方图H的方差值。

(3-e)根据直方图参数计算两种置信度下居民地提取的阈值，高置信度提取阈值为直方图峰值peak，低置信度阈值则根据所要求的置信度依据下式计算容差R：

R＝F(vally，mean，variance)

在峰值peak的基础上加入一定容差得到低置信度提取阈值；

(3-f)根据不同的阈值分别在三幅特征图像上提取居民地，逐点标记像素。在当前阈值下，综合三幅特征图像提取结果，逐点标记像素，若该像素在三幅特征图像上均被标记为居民地，则该像素在当前阈值下标记为居民地，否则标记为背景。从而得到高置信度下居民地提取结果和低置信度下居民地提取结果两幅结果图；

(3-g)金字塔还原，将检测结果还原至原始分辨率；

(3-h)将两种置信度下居民地提取结果输出为一幅结果图，逐点标记像素，若该像素在高置信度提取结果中标记为居民地，则标记为255；若该像素在高置信度提取结果中标记为背景，而在低置信度提取结果中标志为居民地，则标记为128；若该像素在高置信度和低置信度提取结果中同时标记为背景，则标记为0。结果如图5；

(4)水系变化检测的具体实施步骤如下：

(4-a)根据旧矢量数据和高置信度下提取居民地结果，其中高置信度下居民地提取结果即图5中灰度值为255的区域，计算水系变化率。若水系旧矢量区域在高置信度居民地提取结果上为居民地，则原水体消失。统计原水体消失区域面积，除以原水体的总体面积，得到水系变化率。计算得到武汉地区水系变化率为2.5％，北京地区水系变化率为9.7％；

(4-b)根据水系变化率确定水系变化检测的策略。设定阈值为5％，若变化率大于给定阈值，采用策略一，进入步骤(4-c)，否则采用策略二，进入步骤(4-f)。北京地区水系变化检测进入步骤(4-c)，武汉地区水系变化检测进入步骤(4-f)；

(4-c)由高置信度居民地提取结果在原水系区域内采集训练样本，统计居民地区域样本的直方图，并计算直方图参数，包括峰值，方差，均值，谷等，具体计算方法同步骤3中(3-d)；

(4-d)根据直方图参数获得判断水体转变为居民地区域的阈值，即得到水体消失阈值，计算方法同步骤3中(3-e)；

(4-e)在边缘密度图像上逐像素判断，若该像素边缘密度大于消失阈值，则标定为消失，进入步骤(j)；

(4-f)利用栅格化的二值图像在正射影像上采集水系的训练样本；

(4-g)统计训练样本的特征直方图，并计算直方图参数，包括峰值，方差，均值，谷等，具体计算方法同步骤3中(3-d)；

(4-h)根据直方图参数计算两种置信度下水体提取的阈值，两种置信度下阈值的计算方法同步骤3中(3-e)；

(4-i)在两种不同的置信度下，通过阈值和边缘密度图像提取水系；

(4-j)输出结果。北京地区水体消失率大于给定阈值，根据策略一，输出消失水体提取结果图，如图10，武汉地区水体消失率小于给定阈值，根据策略二，将两种置信度下居民地提取结果输出为一幅结果图，逐点标记像素，若该像素在高置信度提取结果中标记为居民地，则标记为255；若该像素在高置信度提取结果中标记为背景，而在低置信度提取结果中标志为居民地，则标记为128；若该像素在高置信度和低置信度提取结果中同时标记为背景，则标记为0，输出水体提取结果图，如图6。

(5)变化检测及结果输出。将居民地提取结果如图5，消失水体提取结果如图10，或水体提取结果如图6，分别与对应的旧矢量数据叠加，获得变化区域和变化类型，具体方法如下：

(5-a)在水系检测中，变化阈值大于给定阈值，则直接转入步骤(5-e)；

(5-b)若所有提取结果图(图5和图6)均处理，转入(5-f)。否则，取某一提取结果图与旧矢量栅格图叠加。逐点判断提取图中像素，若该像素被标记为高置信度居民地或水系，即在提取结果中标记为255，并且该像素在旧矢量栅格图上被标记为背景点，则将该点标记为新增；否则，将该点标记为背景。依次判断各像素，得到二值图像；

(5-c)对二值图像进行形态学处理，消除矢量编辑中多边形趋近而引起的伪变化；

(5-d)逐点判断旧矢量栅格图上被标记为居民地或水系的像素，若该像素在提取结果图上被标记为背景，则将该点标记为消失；否则标记为背景；

(5-e)采用(5-c)操作消除伪变化区域；

(5-f)根据不同的标记值，赋予结果图像各像素不同的颜色值，并提供给用户，新增区域灰度设为255，消失区域灰度设为128，武汉地区居民地变化检测结果如图7，武汉地区水系变化检测结果如图8。北京地区水系消失检测结果如图11。