CN101976437B - 基于自适应阈值分割的高分辨率遥感影像变化检测方法 - Google Patents

Abstract

基于自适应阈值分割的高分辨率遥感影像变化检测方法，首先将两个时期的高分辨遥感影像进行几何配准，然后利用差值法获得二者的差值图像，再通过最大类间方差算法获得变化目标提取的灰度值阈值，并利用灰度值阈值提取出变化目标像素。在检测后处理阶段，针对较杂乱变化图斑，利用形态学变换消除大部分检测出的噪声像素；再利用区域增长算法结合所需提取的变化区域包含像素个数的阈值设定，提取出需要大小的变化区域；最后用边缘检测的方法获得其变化区域的边缘轮廓，再叠加到原始影像上，描绘出完整的变化区域。本发明方法无需过多人工参与，速度较快，可达到良好的检测精度和目视效果。

Description

基于自适应阈值分割的高分辨率遥感影像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及影像土地覆盖变化检测的方法，适用于两幅不同时期具备相同土地覆盖范围的单波段或多波段高分辨遥感影像的变化检测。

背景技术

土地覆盖作为人类社会和自然界交互作用的界面具有重要的社会经济意义和生态环境价值。近年来我国城市发展迅速，及时、周期性获取土地覆盖变化的信息有利于管理者掌握城市土地覆盖的变化趋势，在城市建设中统筹考虑生态环境和经济发展的因素，以实现城市的健康可持续发展。随着我国卫星遥感技术的迅速发展，卫星遥感影像由于周期性、连续性、覆盖范围广等特点为我国土地覆盖变化调查提供了有力的数据支持，提高了土地调查的效率。然而受分辨率的限制，往往不能得到较为理想的土地覆盖变化检测结果。随着越来越多的民用高分辨卫星遥感影像的出现，利用高分辨遥感影像进行土地覆盖变化检测将成为一种趋势。与中低分辨率的遥感影像相比，高分辨影像一方面能够反映较为丰富的地物信息，另一方面由于其较多的地物纹理细节给变化检测带来一定难度。

传统的遥感影像变化检测方法一般可分为两类：一类是对不同时间的数据在像素级进行分析，例如差值法、比值法、回归分析法、相关分析法等；另一类是对不同时间的数据相互独立的分类，然后对分类结果进行比较性研究，一般常用的有分类后比较法等。前一种方法根据像素值，提取有关变化性质的信息以找出感兴趣的变化，但要求精确的时域标准化和图像配准。后一种变化检测方法直接给出了有关地物性质的信息，因而受配准及时域标准化结果的影响小，但是所使用的分类器的准确度对它们的影响很大，而且要想达到较好的分类效果，往往需要过多的人工参与。

尽管上述两类变化检测方法已广泛应用到诸多领域，但是仍然存在由于高分辨遥感影像具有较多的地物纹理细节而使得检测结果存在较多的噪声；或是变化检测过程较复杂，需要较多的人工参与，从而降低了变化检测效率的问题。

发明内容

本发明的技术解决问题是：克服现有技术检测噪声较多，检测结果目视效果不好的不足，提供一种速度快、精度好、人工参与较少的基于自适应阈值分割的高分辨率遥感影像变化检测方法。

本发明的技术解决方案是：基于自适应阈值分割的高分辨率遥感影像变化检测方法，步骤如下：

(1)对于同一空间覆盖范围的两个不同时期的遥感影像，以其中一幅作为参考影像，另一幅与之进行几何配准；

(2)利用差值法获取各个波段的差值图像，差值图像上的像素灰度值为经步骤(1)处理后的两个不同时期遥感影像上各个波段相对应图像位置像素灰度值之差的绝对值；

(3)对步骤(2)得到的各个波段的差值图像应用最大类间方差算法求得各个波段的差值图像上变化区域像素的灰度提取阈值，将各个波段差值图像上小于对应波段灰度提取阈值的所有像素的灰度值置为0，各个波段差值图像上不小于对应波段灰度提取阈值的所有像素的灰度值置为相同的非0值，然后将经过阈值分割的各个波段二值化差值图像合并仍得到一个二值化的差异图像；

(4)利用形态学变换去掉步骤(3)得到的差异图像上的噪声像素；

(5)应用区域增长的方式将步骤(4)所得图像上相互邻接的灰度值为非0值的像素先归并至同一区域，再根据所需提取的变化区域的大小确定区域所包括的像素个数阈值，将所包含像素个数小于像素个数阈值的区域中所有像素的灰度值置为0；

(6)通过Sobel边缘检测算子提取经步骤(5)处理后的差值图像上变化区域的边缘轮廓；

(7)将步骤(6)得到的各变化区域的边缘轮廓线叠加在步骤(1)中的参考影像上，得到遥感影像变化检测结果。

所述步骤(1)中对遥感影像进行几何配准的方法为：在两幅图像上选取至少6对同名控制点对，采用二次多项式建立两幅图像的函数关系，对待配准的图像进行空间坐标变换，再采用最近邻元法进行灰度内插，确定待配准图像各像素的灰度值。

所述步骤(3)中单波段差值图像提取阈值的确定方法为：

(31)将单波段差值图像按照灰度值划分为L个等级，统计各等级的像素个数N_i，i＝1，2，3...L，计算得到整幅差值图像的像素总个数N＝N₁+N₂+...+N_i+...+N_L；

(32)计算单波段差值图像上各等级像素出现的概率P_i＝N_i/N；

(33)通过公式σ²(K)＝a₀a₁(u₀-u₁)²计算类间方差σ²(K)，σ²(K)的最大值对应的K即为提取阈值，上式中a₁＝1-a₀，

所述步骤(4)中的形态学变化方法包括腐蚀和膨胀，结构元素采用3乘3的窗口，4邻域运算；对于腐蚀，当窗口中心元素本身为0或者窗口中心元素的4个邻域中有一个为0时，将窗口中心元素设为0；对于膨胀：当窗口中心元素本身不为0或者窗口中心元素的4个邻域中有一个不为0时，则将窗口中心元素设为相同的非零值；对差异图像进行形态学变换时，采用先腐蚀再膨胀的方法，去掉噪声像素。

所述步骤(5)中的区域增长方法为：

(51)设定最小区域阈值，最小区域阈值由所需要得到的变化检测区域的大小确定；

(52)对步骤(4)得到的二值化图像按照从左至右，从上至下的顺序进行逐像素的遍历，当遇到一个满足灰度值不为0且未进行区域标记的像素时，停止遍历，并将该像素作为此次区域增长运算的第一个种子点，对该种子点进行区域标记，并将该种子点存入堆栈，进行下一步；若遍历完整幅图像，未找到满足上述条件的像素时，则步骤(5)结束；

(53)顺序从存储种子点的堆栈中取出种子点，在步骤(4)得到的二值化图像上，遍历种子点的八个邻域内的像素，将灰度值不为0且未进行区域标记的像素与种子点合并至一个区域，同时将这些新合并的像素作为种子点顺序存入堆栈；

(54)重复步骤(53)，直到堆栈为空；

(55)统计步骤(54)得到的变化区域所包含的像素个数，若小于预先设定的最小区域阈值，则将该变化区域中所有像素的灰度值置0，并取消其区域标记；

(56)重复步骤(52)～(55)，进行下一个区域增长的运算，直至满足步骤(5)结束的条件。

所述步骤(6)中通过Sobel边缘检测算子提取边缘的方法为：对经步骤(5)处理后的差值图像进行如下运算：

fx′＝|f(x+1，y-1)+2*f(x+1，y)+f(x+1，y+1)-f(x-1，y-1)-2*f(x-1，y)-f(x-1，y+1)|，

fy′＝|f(x-1，y+1)+2*f(x，y+1)+f(x+1，y+1)-f(x-1，y-1)-2*f(x，y-1)-f(x+1，y-1)|，

上式中f(x，y)表示在图像上第x行，第y列处像素的灰度值，fx′为像素(x，y)在x方向上的梯度值，fy′为像素(x，y)在y方向上的梯度值；当fx′或fy′不为0时，将f(x，y)置为相同的非零值，当fx′，fy′都为0时，将f(x，y)置为0。

本发明与现有技术相比的优点在于：

(1)本发明方法采用最大类间方差算法获取差值图像上变化目标与非变化目标像素提取的灰度值阈值，具有自适应的特点，对于任何光学遥感图像都适用，在很大程度上减少了人工的参与，简化了变化区域像素提取的过程；

(2)本发明方法采用形态学变换，可消除大部分高分辨遥感影像因纹理细节较多，在通过灰度值阈值提取变化像素之后产生的众多噪声像素：边缘毛刺、区域内的岛、洞的现象及因几何配准误差导致的错检像素；

(3)本发明方法在区域增长的过程中加入区域大小的阈值设置，对由噪声像素所合并成的小区域再次进行了去除，提高了检测的精度；

(4)本发明方法最终的检测结果为一块块完整的区域，采用边缘检测方法检测出变化区域的边缘像素，获得变化区域边缘轮廓，将其叠加到原始影像上可清楚地勾绘出变化区域，使检测者一目了然；

(5)本发明方法对已配准的遥感影像，在检测之前只需对形态学变化次数以及变化区域大小的阈值进行简单的设置后即可自动完成，操作简单、快捷。

附图说明

图1为本发明方法的实现流程图；

图2为两个不同时期相同区域的单波段高分辨遥感影像；

图3为采用最大类间方差算法获取图2中两个不同时期影像的差异图像；

图4为采用本文发明的方法所获得图2中两个时期影像的变化区域边缘；

图5为两个不同时期相同区域的多段高分辨遥感影像；

图6为采用最大类间方差算法获取图5中两个不同时期影像的差异图像；

图7为采用本文发明方法所获得图5中两个个时期影像的变化区域边缘。

具体实施方式

如图1所示，本发明的具体实施步骤如下：

步骤1：对于同一空间覆盖范围的两个不同时期的遥感影像，以其中一幅作为参考影像，另一幅与之进行几何配准，配准原理如下：首先在两幅图像上选取至少6对同名控制点对，采用二次多项式建立两幅图像的函数关系，对待纠正影像进行空间坐标变换；然后再对待纠正的图像进行灰度内插，采用最近邻元法确定待配准图像各像素的灰度值，由于纠正后的影像各个像素的空间坐标非整数，所以需要将原图像上离此空间点最近的像素的灰度值赋给该空间点，形成重采样后的像素，采用此方法可保证改变影像原有的灰度值。

首先进行几何配准，是考虑到由于卫星轨道漂移、姿态、拍射角度的变化，使得两个时期的影像发生几何定位误差，从而导致两幅影像上相对位置所对应的实际像点不同，几何配准可消除这种误差，配准的精度控制在1个象元之内即可。

步骤2：利用差值法获取两幅影像各个波段的差值图像，差值图像上的像素灰度值为经配准后两个不同时期遥感影像上各个波段相对应图像位置像素灰度值之差的绝对值。

步骤3：将步骤2得到的各个波段的差值图像应用最大类间方差算法求得各个波段的差值图像上变化区域像素的灰度提取阈值，将各个波段差值图像上小于对应波段灰度提取阈值的所有像素的灰度值置为0，各个波段差值图像上不小于对应波段灰度提取阈值的所有像素的灰度值置为相同的非0值，然后将经过阈值分割的各个波段的二值化差值图像合并得到一个二值化的差异图像；具体操作方法如下：

(31)对于单波段差值图像的灰度值按从小到大的顺序划分为L个等级，等级范围设为[0，2，...L]，其中L对应影像的最大灰度值，统计各灰度级i的像素个素N_i，计算整幅影像的像素总个数N＝N₀+N₁+...+N_L；

(32)计算单波段差值图像各个灰度级像素出现的概率：P_i＝N_i/N；

(33)将单波段差值图像分为背景区域C₀与变化区域C₁两类，C₀类灰度级范围为[1，...K]，C₁类为[K+1，...L]，K为变化区域提取的灰度阈值，通过类间方差σ²最大求单波段差值图像灰度阈值K：由于类间方差σ²是K的函数：σ²(K)＝a₀a₁(u₀-u₁)²，其中a_i为类C_i的面积与图像总面积的比，则：

a₁＝1-a₀，u_i为类C_i的均值，则

分别将K_i代入σ²(K)式计算得出σ_i ²，其中K_i为单波段差值图像应灰度等级为i的像素的灰度值，通过求σ²最大得到变化区域像素灰度值提取阈值K；

(34)对单波段差值图像根据阈值K进行二值化，获取二值化的差值图像，具体做法是对差值图像上灰度值小于K的像素赋值为0，灰度值大于或等于K的像素赋值为相同的非零值；

(35)对于多波段的影像，重复步骤(31)至(34)，获取各个波段的二值化的差值图像，然后对各个波段提取灰度值为非零值的像素进行合并，具体做法是：对各个波段的二值化差值图像上的像素灰度值进行判断，在合并后的二值图像上，将各个波段对应图像位置上灰度值为非零值的像素范围取并集，仍赋值相同的非零值，其余像素灰度值均赋为0。

步骤4：步骤3得到的差值图像存在很多噪声像素，例如：区域内较小的岛或者洞；变化区域边缘个别突出的像素所形成的毛刺；由邻接的变化像素所形成的狭长的区域；由配准误差引起的细小区域等。利用形态学变换去掉大部分的噪声像素：腐蚀则去掉边缘毛刺，及狭长、细小的区域以及变化区域内的小岛，膨胀则可去掉岛和洞。本实施方案所采用的形态学变换的结构元素为3乘3的窗口，为提高效率进行4邻域运算；对于腐蚀，当窗口中心元素本身灰度值为0或者窗口中心元素的4个邻域中有一个灰度值为0时，将窗口中心元素灰度值设为0；对于膨胀：当窗口中心元素本身灰度值不为0或者窗口中心元素的4个邻域中有一个灰度值不为0时，则将窗口中心元素设为相同的非零值。

在具体操作时，首先设定腐蚀和膨胀的次数，然后在步骤3所得二值化差值图像上按照从上到下，从左到右的顺序，依次逐像素执行腐蚀运算，完成一次腐蚀运算后，在得到的图像上再进行一次腐蚀运算，直到达到指定的运算次数为止；然后再在所得图像上执行指定次数的膨胀运算。此步骤中腐蚀和膨胀的次数可以不相同，但是相差不宜超过2次，否则会影响检测的精度。

步骤5：应用区域增长的方式将步骤4所得图像上相互邻接的灰度值为非0值的像素先归并至同一区域，再根据所需提取的变化区域的大小设定区域所包括的像素个数阈值，去掉小区域。具体操作过程如下：

(51)设定最小区域阈值，这个阈值的设定与检测者所要得到的变化检测区域的大小有关，例如检测者如果指定一个区域中若有a＊b个像素发生变化，则认为该区域是发生变化的，其中a，b分别为变化区域的行列数，则初始阈值可设定为a＊b；

(52)定义一个数组f[xsize][ysize]，其中xsize，yszie分别为图像的行数和列数，用f[x][y]为第x行，y列的像素所属区域做标记，首先将所有像素f[x][y]的值初始化0，表示图像上所有像素尚未进行区域的划分；定义一个堆栈d[xszie＊yszie]，用于存储种子点像素的位置；

(53)对步骤4得到的二值化图像按照从左至右，从上至下的顺序进行逐像素的遍历，当遇到一个满足灰度值不为0且未进行区域标记即f[x][y]＝0的像素时，停止遍历，并将该像素作为此次区域增长运算的第一个种子点，并对该点进行区域标记：令f[x][y]＝n(n依次取1，2，3......)，并将该点的位置：x＊yszie+ysize存放入堆栈d；若遍历完整幅图像，未找到满足上述条件的像素时，则步骤5结束；

(54)从堆栈d中顺序取出种子点，在步骤4得到的二值化图像上，遍历种子点的八个邻域内的像素，将灰度值不为0且未进行区域标记的像素合并至种子点所在的区域，并进行区域标记，即令各个像素的F[x][y]＝n；并将这些新合并的像素作为新的种子点，顺序放入堆栈d；

(55)重复进行步骤(54)的操作，直到堆栈d为空，此时完成了一个区域的生长过程。

(56)统计步骤(55)得到的变化区域所包含的像素个数，若小于预先设定的最小区域阈值，则将该变化区域中所有像素的灰度值置0，并取消其区域标记，令f[x][y]＝0；否则令区域标记n＝n+1；

(57)重复步骤(53)～(56)，进行下一个区域增长的运算，直至满足则步骤5结束的条件时，完成骤4得到的二值化图像上所有变化区域的像素归并及小区域的去除工作。

步骤6：通过Sobel边缘检测算子提取经步骤5处理后的差值图像上变化区域的边缘轮廓；Solbel算子如下：

fx′＝|f(x+1，y-1)+2*f(x+1，y)+f(x+1，y+1)-f(x-1，y-1)-2*f(x-1，y)-f(x-1，y+1)|，

fy′＝|f(x-1，y+1)+2*f(x，y+1)+f(x+1，y+1)-f(x-1，y-1)-2*f(x，y-1)-f(x+1，y-1)|

上式中f(x，y)表示在图像上第x行，第y列处像素的灰度值，fx′为像素(x，y)在x方向上的梯度值，fy′为像素(x，y)在y方向上的梯度值；当fx′或fy′不为0时，将f(x，y)置为相同的非零值，当fx′，fy′都为0时，将f(x，y)置为0。

步骤7：将步骤6得到的各变化区域的边缘轮廓线叠加在步骤1的参考影像上，得到遥感影像变化检测结果。

实施例

1、实验数据

本实验数据分为两组，第一组为单波段不同时相覆盖相同区域的高分辨遥感影像，如图2所示：影像来源为分辨率2.36米的CBERS HR全色影像。为保证土地覆盖变化检测的时效性，本文选取时期较为相近的影像，分别于2008年1月1日、2008年12月30日拍摄，图像大小均为533＊589(单位：像素)；第二组为包含3个波段的不同时相覆盖相同区域的高分辨遥感影像，如图5所示：时相分别为2002年，2003年，图像大小均为1024＊1024(单位：像素)。两组数据在所采用的形态学变换过程均为先进行5次腐蚀运算，再进行6次膨胀运算，最终所提取的变化区域所包含的像素个数均大于400个。

2、实验结果分析

表1说明了两组数据在执行本算法的各个步骤阶段，所提取的变化区域个数的变化情况：在进行完本算法第三步，即采用最大类间方差法获取提取变化像素的灰度值阈值后，第一组数据提取的变化区域为2033个，如图3所示，第二组数据提取的变化区域为4373个，如图6所示，图中可以看出存在很多细小的噪声像素；在用形态学变换进行部分细小区域(岛、洞，以及因配准误差所引起的噪声)的去除后，所提取变换区域个数第一组为56个，第二组为182个，分别减少了97％，和96％，可见通过形态学变化可去掉绝大部分的细小区域；在运用区域增长算法将相互邻接的变化区域内在像素在逻辑上归并为一个区域后，通过设计最小区域像素个数的阈值(本实验数据均为400)，即提取出所需大小的变化区域，第一组数据为19个，如图4所示，第二组数据为51个，如图7所示。

表2说明了本发明方法所检测出来的结果与实际目测结果的对比：第一组数据共检测出19个变化区域，如图4所示，其中错检区域为1个，漏检的区域为3个；第二组数据共检测出51个变化区域，如图7所示，其中错检区域为10个，漏检的区域为8个；本算法所计算的变化检测精度的公式为：正确检测出变化区域的个数与实际目测检测出来的变化区域个数的比值，即100％＊(本算法检测出的变化区域的个数-错检的变化区域个数)/(本算法检测出的变化区域的个数-错检的变化区域个数+漏检变化区域的个数)，计算两组数据的变化检测精度分别为：第一组86％；第二组84％。

表1变化区域个数变化情况表

表2本发明方法所检测出来的结果与实际目测结果的对比表

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (1)

1.基于自适应阈值分割的高分辨率遥感影像变化检测方法，其特征在于步骤如下：

(1)对于同一空间覆盖范围的两个不同时期的遥感影像，以其中一幅作为参考影像，另一幅与之进行几何配准，配准的方法为：在两幅图像上选取至少6对同名控制点对，采用二次多项式建立两幅图像的函数关系，对待配准的图像进行空间坐标变换，再采用最近邻元法进行灰度内插，确定待配准图像各像素的灰度值；

(2)利用差值法获取各个波段的差值图像，差值图像上的像素灰度值为经步骤(1)处理后的两个不同时期遥感影像上各个波段相对应图像位置像素灰度值之差的绝对值；

(3)对步骤(2)得到的各个波段的差值图像应用最大类间方差算法求得各个波段的差值图像上变化区域像素的灰度提取阈值，将各个波段差值图像上小于对应波段灰度提取阈值的所有像素的灰度值置为0，各个波段差值图像上不小于对应波段灰度提取阈值的所有像素的灰度值置为相同的非0值，然后将经过阈值分割的各个波段二值化差值图像合并仍得到一个二值化的差异图像，差异图像上的非0值像素取各个波段非0值像素的并集；

单波段差值图像提取阈值的确定方法为：

(31)将单波段差值图像按照灰度值划分为L个等级，统计各等级的像素个数N_i，i＝1，2，3...L，计算得到整幅差值图像的像素总个数N＝N₁+N₂+...+N_i+...+N_L；

(32)计算单波段差值图像上各等级像素出现的概率P_i＝N_i/N；

(33)通过公式σ²(K)＝a₀a₁(u₀-u₁)2计算类间方差σ²(K)，σ²(K)的最大值对应的K即为提取阈值，上式中

a₁＝1-a₀，

(4)利用形态学变换去掉步骤(3)得到的差异图像上的噪声像素；形态学变换方法包括腐蚀和膨胀，结构元素采用3乘3的窗口，4邻域运算；对于腐蚀，当窗口中心元素本身为0或者窗口中心元素的4个邻域中有一个为0时，将窗口中心元素设为0；对于膨胀：当窗口中心元素本身不为0或者窗口中心元素的4个邻域中有一个不为0时，则将窗口中心元素设为相同的非零值；对差异图像进行形态学变换时，采用先腐蚀再膨胀的方法，去掉噪声像素；

(5)应用区域增长的方式将步骤(4)所得图像上相互邻接的灰度值为非0值的像素先归并至同一区域，再根据所需提取的变化区域的大小确定区域所包括的像素个数阈值，将所包含像素个数小于像素个数阈值的区域中所有像素的灰度值置为0；

区域增长方法为：

(51)设定最小区域阈值，最小区域阈值由所需要得到的变化检测区域的大小确定；

(52)对步骤(4)得到的图像按照从左至右，从上至下的顺序进行逐像素的遍历，当遇到一个满足灰度值不为0且未进行区域标记的像素时，停止遍历，并将该像素作为此次区域增长运算的第一个种子点，对该种子点进行区域标记，并将该种子点存入堆栈，进行下一步；若遍历完整幅图像，未找到满足上述条件的像素时，则步骤(5)结束；

(53)顺序从存储种子点的堆栈中取出种子点，在步骤(4)得到的图像上，遍历种子点的八个邻域内的像素，将灰度值不为0且未进行区域标记的像素与种子点合并至一个区域，同时将这些新合并的像素作为种子点顺序存入堆栈；

(54)重复步骤(53)，直到堆栈为空；

(55)统计步骤(54)得到的变化区域所包含的像素个数，若小于预先设定的最小区域阈值，则将该变化区域中所有像素的灰度值置0，并取消其区域标记；

(56)重复步骤(52)～(55)，进行下一个区域增长的运算，直至满足步骤(5)结束的条件；

(6)通过Sobel边缘检测算子提取经步骤(5)处理后的图像上变化区域的边缘轮廓，方法为：对经步骤(5)处理后的图像进行如下运算：

fx′＝|f(x+1，y-1)+2*f(x+1，y)+f(x+1，y+1)-f(x-1，y-1)-2*f(x-1，y)-f(x-1，y+1)|，

fy′＝|f(x-1，y+1)+2*f(x，y+1)+f(x+1，y+1)-f(x-1，y-1)-2*f(x，y-1)-f(x+1，y-1)|，

上式中f(x，y)表示在图像上第×行，第y列处像素的灰度值，fx′为像素(x，y)在x方向上的梯度值，fy′为像素(x，y)在y方向上的梯度值；当fx′或fy′不为0时，将f(x，y)置为相同的非零值，当fx′，fy′都为0时，将f(x，y)置为0；

(7)将步骤(6)得到的各变化区域的边缘轮廓线叠加在步骤(1)中的参考影像上，得到遥感影像变化检测结果。

Images