CN103578110B - 基于灰度共生矩阵的多波段高分辨率遥感影像分割方法 - Google Patents

Abstract

本发明公开一种基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，在基于灰度共生矩阵的纹理影像中采用降水分水岭变换单独分割每个波段影像，进而对各波段分割结果进行叠加。最后，提出了一种的基于多波段光谱信息的区域合并策略合并过分割结果中的碎片区域，最终实现影像分割。本发明分别对高分辨率ALOS和SPOT？5影像进行了实验，并与传统基于局域同质性梯度的分割方法进行了比较。实验结果表明，本发明能够准确定位对象的边缘，有效克服过分割及欠分割现象，具有更高的分割精度与稳定性。

Description

基于灰度共生矩阵的多波段高分辨率遥感影像分割方法

技术领域

本发明涉及一种基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，属于遥感影像分割技术领域。

背景技术

目前，高分辨率遥感影像在城市规划、环境评价及军事等很多领域都得到广泛应用，而面向对象的分析方法在高分辨率遥感影像研究中正越来越多的受到重视。有效的影像分割是面向对象的分析方法的基础和重要保证。目前，高分辨率遥感影像的分割方法主要分为五类：基于像元的分割方法，如阈值法、聚类法；基于边缘检测的分割方法；基于区域的分割方法；基于物理模型的分割方法和结合特定数学理论；工具的分割方法，如数学形态学、模糊数学、小波变换和人工神经网络等。

随着航天技术与传感器技术的不断发展，遥感影像空间分辨的不断提高，同时也为影像分割带来了新的挑战。与传统中低分辨率遥感影像相比，高分辨遥感影像中“同谱异物”现在依然存在，“同物异谱”现象更为突出。尤其在城市场景中，丰富地物的种类，多样的生态环境以及复杂的图像背景都给准确的影像分割造成了困难。在众多分割方法中，分水岭算法具有快速、有效、准确特点而受到广泛应用。但在高分辨率遥感影像中，在图像中通常存在大量的暗纹理细节，从而造成了大量的伪局部极小值，这些极小值产生的伪积水盆地会造成严重的过分割问题。针对这一问题，Vincent等提出的沉浸分水岭分割算法有效提高了分割精度，但是计算效率不高。Semt等提出了一种模拟降水的分水岭算法，该算法速度更快且没有保证了分割精度。在此基础上，陈秋晓等提出了一种新颖的基于局域同质性梯度的高分辨率遥感影像分割方法，该算法有效的解决了分水岭变换的过分割问题，但也存在如下问题：由于采用了过于粗糙的量化，使高分辨率遥感影像中大量细节信息的严重损失，从而容易产生欠分割现象，也很难准确定位对象的边缘；由于高分辨率遥感影像中存在大量具有相似梯度特征的对象，而此方法仅采用了梯度信息作为提取对象的依据，因此需要引入更多的特征信息对这些对象加以区分；最后，此方法所采用的区域合并策略仅对梯度影像中对象内部的灰度信息进行了比较，而忽略了高分辨率遥感影像中不同对象在不同波段影像中的光谱差异，从而可能造成误合并现象。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供一种新颖的多波段-多特征高分辨率遥感影像分割方法。首先利用灰度共生矩阵计算各个波段的纹理图像，进而采用降分水岭变换获得单一波段的分割结果。进而将各波段分割结果进行叠加获得过分割的碎片区域，最后采用基于多波段光谱信息的区域合并策略合并碎片区域。

技术方案：一种基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，主要分为三个步骤：波段影像纹理影像特征提取；降水分水岭分割；分割结果叠加与区域合并。

1，波段影像纹理影像特征提取

对每个波段影像单独进行分割，采用灰度共生矩阵方法提取遥感影像的多个纹理特征，计算各个波段的纹理影像：根据图像大小，选择尺寸为N×N像素的滑动窗口来计算灰度共生矩阵，计算灰度共生矩阵的熵、角二阶矩、对比度和一致性四个特征值，作为新的特征矩阵的窗口中心像元。移动纹理窗口，重复上述步骤，直到遍历整个影像，可以得到四个新的纹理矩阵。将得到的四个方向的纹理矩阵取平均值，获得最终的纹理特征矩阵，即纹理影像。

2，降水分水岭分割

基于所提取的纹理图像，采用降水分水岭变换方法对每个波段影像分割进行分割。

3，分割结果叠加与区域合并

将多波段影像中的分割结果中的边界进行叠加，从而将原始影像过分割为大量碎片区域，在此基础上对这些区域进行区域合并，获得最终的分割结果。

有益效果：与现有技术相比，本发明提供的基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，针对高分辨率影像的特点，算法对各波段影像单独进行了分割，且避免了影像量化造成的图像细节的损失；进而采用度共生矩阵提取图像的多个纹理特征增加的分割结果的可靠性，最后提出了一种改进的区域合并策略。通过对不同类型传感器，不同空间分辨率的两组遥感影像实验表明，本发明所提方法能够精确提取图像中的细节信息，准确定位对象的边缘，有效避免过分割及欠分割现象。与文献1所提出的方法相比，无论目视分析还是定量分析均表明，对于细节信息丰富的遥感影像，可以很好地解决常规算法欠分割的问题，而对于纹理特征单一的区域，也可以很好的避免过分割的情况，因此，本发明具有更高的分割精度及更好的分割效果。

附图说明

图1为本发明实施例的ALOS融合影像图；

图2为数据集1实验结果图，其中(a)为本发明实施例实验结果，(b)为文献1方法实验结果；

图3为2009年采集的中国上海地区的三波段全色-多光谱融合影像；

图4为数据集2实验结果图，其中(a)为基于局域同质性梯度分割算法实验结果，(b)为本发明实施例实验结果。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，包括如下步骤：

1多波段纹理影像生成

由于过于粗糙的量化会损失图像的细节信息，不对原始影像进行量化，而对每个波段影像单独进行分割。遥感影像中的纹理特征反映了对象的空间排列信息，而结合图像的纹理特征则能更加有效的提取遥感图像中的有用信息，从而提高图像分割的精度。灰度共生矩阵是图像纹理分析方法中的一种，在一定程度上反映了纹理图像中的各灰度级在空间上的分布特性，是纹理分析领域中最经常采用的方法之一。

灰度共生矩阵的数学基础是图像的二阶矩组合条件概率密度函数，其定义为：对于一个由N个离散灰度级组成的图像f(x,y)来说，灰度共生矩阵P_d,θ[i,j]定义为i，j联合出现的概率，其中，i是点(x,y)的灰度值，j是点(x+Δx,y+Δy)的灰度值，d是指灰度共生矩阵的生长步长；θ是指灰度共生矩阵的生成方向，是指满足d、θ条件的灰度值分别为i，j的点对数，是指满足d、θ条件的所有点对数；共生矩阵P(d,θ)是一个大小为N_q×N_q的方阵：

$P_{d,\mathrm{\θ}}\[i,j\]=\frac{P_{ij}^{d\mathrm{\θ}}}{\underset{i=1}{\overset{N_q}{\Σ}}\underset{j=1}{\overset{N_q}{\Σ}}{P}_{ij}^{d\mathrm{\θ}}}---\left(1\right)$

采用灰度共生矩阵方法提取遥感影像的多个纹理特征，计算各个波段的纹理影像：根据图像大小，选择尺寸为N×N像素的滑动窗口来计算灰度共生矩阵，确定公式(1)步长参数d。计算灰度共生矩阵的方向角θ取0°、45°、90°、135°，从而获得四个方向上的灰度共生矩阵。

Haralick等人在灰度共生矩阵的基础上提出了14种相关纹理量化的统计量，选取灰度共生矩阵的4个特征值，计算特征矩阵。根据公式(1)，选取的特征值如下：

①熵

$ENT=-\underset{i}{\Σ}\underset{j}{\Σ}{P}_{d,\mathrm{\θ}}\[i,j\]\lg P_{d,\mathrm{\θ}}\[i,j\]---\left(2\right)$

熵给出了一个图像内容随机性、无序性的度量，反映了图像中纹理的非均匀程度或复杂程度。当P_d,θ[i,j]中的数值均相等时，熵最大，反之，P_d,θ[i,j]的数值之间差别很大时，熵较小。

②角二阶矩

$ASM=\underset{i}{\Σ}\underset{j}{\Σ}{P}_{d,\mathrm{\θ}}^2\[i,j\]---\left(3\right)$

它是图像灰度分布均匀性或平滑性的度量，反映了图像灰度分布的均匀程度和纹理的粗细度。当P_d,θ[i,j]中的数值分布较集中于主对角线附近时，说明从局部区域范围内观察图像的灰度分布是角均匀的，即图像呈现较粗的纹理，相应的角二阶矩值较大，反之较小。

③对比度

$CON=\underset{i}{\Σ}\underset{j}{\Σ}{(i-j)}^2{P}_{d,\mathrm{\θ}}\[i,j\]---\left(4\right)$

对比度反映图像的清晰度和纹理沟纹深浅的程度，对于粗纹理，P_d,θ[i,j]的数值较集中于主对角线附近。对比度能够与有效检测图像反差，提取图像边缘信息，增强线性构造等信息。

④一致性

$HOM=\underset{i}{\Σ}\underset{j}{\Σ}\frac{1}{1+{(i-j)}^2}{P}_{d,\mathrm{\θ}}\[i,j\]---\left(5\right)$

一致性反映图像纹理的同质性，度量图像纹理局部变化的多少，是区分不同目标的重要度量。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀。

根据得到的灰度共生矩阵计算相应的特征值，作为新的特征矩阵的窗口中心像元。移动纹理窗口，重复上述步骤，直到遍历整个影像，可以得到四个新的纹理矩阵。将得到的四个方向的纹理矩阵取平均值，获得最终的纹理特征矩阵，即纹理影像。

2降水分水岭分割

基于所提取的纹理图像，采用Smet和Pires等提出的降水分水岭变换方法对每个波段影像分割进行分割，降水分水岭变换的实质是寻找从图像的每个像素到图像表面高程的局部较小的下游路径。而集水盆地则定义为满足一下条件的所有像素的集合：这些像素的下游路径终止于同一个高程极小点。降水分水岭的具体步骤如下：

(1)为每一个像元寻找下游像元，并记录在数组中。对于每一个像元，寻找与其相比灰度最小的邻域像元。

(2)标记局部最小值的点。判断数组中的每一个元素是否为局部最小，若是的话则赋予一个新的标号，将该标号同样赋予与其连通并均为局部最小的区域。

(3)标识非局部最小像元。对于每一个非局部最小像元P，总是存在一个下游像元。若该下游像元已被标识，则将该标号赋予P，否则寻找下游像元的下游像元，直至找到已标识的下游像元，并将该标号赋予P。

经过降水分水岭变换图像被分割成标号从0开始的一系列标号区域，每一个局部极小值或者说一个具有极小值的连通区域都对应了一个集水盆地。可以看出分水岭变换的过分割问题很严重，为了解决过分割的问题，下面提出一种改进的区域合并策略。

3分割结果叠加与区域合并

将多波段影像中的分割结果中的边界进行叠加，从而将原始影像过分割为大量碎片区域，在此基础上对这些区域进行区域合并，获得最终的分割结果。鉴于文献1的合并策略中仅考虑了量化影像中对象内部的梯度信息，为进一步提高区域合并的准确性，本文提出了一种改进的区域合并策略：

第一步：把叠加后的边界映射到每一个波段影像中，定义过分割结果中区域为总数为n，影像中波段总数为C，某一波段中各区域的颜色均值向量为μ_n＝(μ_l1，μ_l2，…，μ_ln)，其中l＝1,2...C。

第二步：计算任意两个碎片区域R_A和R_B在每个波段中的欧式距离d_l(A,B)＝||μ_lA-μ_lB||，设定阈值T₁，若所有波段的欧式距离均满足d_l(A,B)≤T₁，则直接合并区域R_A和R_B。若没有一个波段满足d_l(A,B)≤T₁，则不进行合并。

第三步：对不满足第二步中两个判别条件的区域，进一步计算每个波段中R_A和R_B和方差σ_lAB，从而获得所有波段的方差均值σ_AB≤T₂，设定阈值T₂，若满足σ_AB≤T₂，则合并R_A和R_B。反之，不进行合并。

第四步：遍历所有区域，获得最终的分割结果。

实验结果与分析

为了验证方法的可行性与有效性，采用两组高分辨率全色多光谱融合影像进行了实验。为了进一步分析方法的性能，选择将实验结果与文献1(陈秋晓，陈述彭，周成虎.基于局域同质性梯度的遥感图像分割方法及其评价[J].遥感学报.2006.Vol.10,No.3,pp.357-365.)提出的方法进行了比较。

数据集1结果分析

实验数据集1采用ALOS卫星2007年采集的三波段全色-多光谱融合影像，空间分辨率为2.5m，所在地区为中国江苏省南京市江宁区河海大学，影像的空间分辨为2.5m。如图1所示。

实验中，由于波段数为3，因此C＝3，设定阈值T₁＝0.6，T₂＝0.3。文献1方法中，各阈值设定为：h-threshold＝39，合并阈值T＝100。实验结果如图2所示。

为便于目视分析分割结果，我们对图像中的典型地物用字母进行了标记，如果1、图2所示。其中，区域A、D和E为建筑物，区域B为操场，区域C为道路，区域F和G为绿化带。通过目视分析可以看出，文献1的分割算法能够基本的分割遥感影像的各个区域，但是图像的一些细节信息却没有很好地识别出来，如位置A、D和E，这些建筑物的纹理特征不是特别明显，因此基于局域同质性梯度的分割算法没有很好的将它们分割出来，但是本发明则对他们进行了精确的分割，解决了欠分割的问题，得到了很好的分割结果。而对于操场B，基于局域同质性梯度的分割算法没能区分出操场的跑道和草坪，而本发明则很好的识别出这些信息。我们还可以看出，本发明对于影像边缘定位比基于局域同质性梯度的分割算法更加准确。同时，对于植被覆盖部分即区域F和G，还有道路区域C，基于局域同质性梯度的分割算法还存在着过分割的现象，而本发明由于提出了改进的区域合并策略，很好的解决了分水岭变换常有的过分割情况。

数据集2结果分析

实验图像采用的是2009年采集的中国上海地区的三波段全色-多光谱融合影像，影像的空间分辨为2.5m；如图3所示。

实验中，由于波段数为3，因此C＝3，设定阈值T₁＝0.6，T₂＝0.3。文献1方法中，各阈值设定为：h-threshold＝39，合并阈值T＝100。实验结果如图4所示。

同样，我们对图像中的部分典型区域用字母进行了标记，如果图3、图4所示。由于该图像地物种类较少，且存在大量内部纹理特征单一且面积较大的区域，由于该图像地物种类较少，且存在大量内部纹理特征单一且面积较大的区域，因此通过两个算法结果的目视对比图我们不难看出，基于局域同质性梯度的分割算法得到的图像分割结果中欠分割的问题仍然存在，不能很好的识别出纹理特征不是很明显的细节区域，从而使得分割结果不尽理想。如图4的A、B、C和D区域，这些区域有一些细微的纹理特征，但是局域同质性梯度的分割算法没有将它们分割出来。且本发明更好的识别出区域A的边缘细节信息，因此对于这类细节信息丰富且纹理特征单一的遥感影像，本发明可以得到很好的分割结果。

定量分析

为进一步定量分析所提出算法的性能，采用Y.Deng.etal提出的评价标准比较两种算法的分割精度。定义分割结果中所有区域的J指标均值为

$\stackrel{\‾}{J}=\frac{1}{N}\underset{k}{\Σ}{M}_k{J}_k---\left(6\right)$

(6)式中，N为图像中像素总数，M_k为第k(k＝1,2...N)个区域的像素总数，J_k为第k个区域的J指标。J指标反映了分割结果中某一区域的颜色分布均质程度。某一算法的分割结果对应的越小，说明该算法的分割结果越好，即所提取的对象内部的平均均质程度越高。根据这一准则，分别对两组实验进行定量分析，结果图表1所示：

表1分割精度定量分析

如表1所示，本发明在两组实验中分割精度都明显高于文献1采用的方法，与目视分析结果一致，进一步验证了本发明的可行性与有效性。另外，第二组实验中两种算法的分割精度较第一组实验均有所提高，由于数据集2中地物种类较少，存在大量内部纹理特征单一且面积较大的区域所造成的。

Claims (2)

1.一种基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，其特征在于，主要分为三个步骤：波段影像纹理影像特征提取；降水分水岭分割；分割结果叠加与区域合并；

1)，波段影像纹理影像特征提取

对每个波段影像单独进行分割，采用灰度共生矩阵方法提取遥感影像的多个纹理特征，计算各个波段的纹理影像：根据图像大小，选择尺寸为N×N像素的滑动窗口来计算灰度共生矩阵，计算灰度共生矩阵的熵、角二阶矩、对比度和一致性四个特征值，作为新的特征矩阵的窗口中心像元；移动纹理窗口，重复上述步骤，直到遍历整个影像，可以得到四个新的纹理矩阵；将得到的四个方向的纹理矩阵取平均值，获得最终的纹理特征矩阵，即纹理影像；

2)，降水分水岭分割

基于所提取的纹理图像，采用降水分水岭变换方法对每个波段影像分割进行分割；

3)，分割结果叠加与区域合并

将多波段影像中的分割结果中的边界进行叠加，从而将原始影像过分割为大量碎片区域，在此基础上对这些区域进行区域合并，获得最终的分割结果；

降水分水岭的具体步骤如下：

(1)为每一个像元寻找下游像元，并记录在数组中，对于每一个像元，寻找与其相比灰度最小的邻域像元；

(2)标记局部最小值的点，判断数组中的每一个元素是否为局部最小，若是的话则赋予一个新的标号，将该标号同样赋予与其连通并均为局部最小的区域；

(3)标识非局部最小像元，对于每一个非局部最小像元P，总是存在一个下游像元；若该下游像元已被标识，则将该标号赋予P，否则寻找下游像元的下游像元，直至找到已标识的下游像元，并将该标号赋予P；

经过降水分水岭变换图像被分割成标号从0开始的一系列标号区域，每一个局部极小值或者说一个具有极小值的连通区域都对应了一个集水盆地；

区域合并策略为：

第一步：把叠加后的边界映射到每一个波段影像中，定义过分割结果中区域为总数为n，影像中波段总数为C，某一波段中各区域的颜色均值向量为μ_n＝(μ_l1，μ_l2，…，μ_ln)，其中l＝1,2...C；

第二步：计算任意两个碎片区域R_A和R_B在每个波段中的欧式距离d_l(A,B)＝||μ_lA-μ_lB||，设定阈值T₁，若所有波段的欧式距离均满足d_l(A,B)≤T₁，则直接合并区域R_A和R_B；若没有一个波段满足d_l(A,B)≤T₁，则不进行合并；

第三步：对不满足第二步中两个判别条件的区域，进一步计算每个波段中R_A和R_B的方差σ_lAB，从而获得所有波段的方差均值σ_AB≤T₂，设定阈值T₂，若满足σ_AB≤T₂，则合并R_A和R_B；反之，不进行合并；

第四步：遍历所有区域，获得最终的分割结果。

2.如权利要求1所述的基于灰度共生矩阵的多波段高分辨率遥感影像分割方法，其特征在于，

对于一个由N个离散灰度级组成的图像f(x,y)来说，灰度共生矩阵P_d,θ[i,j]定义为i，j联合出现的概率，其中，i是点(x,y)的灰度值，j是点(x+Δx,y+Δy)的灰度值，d是指灰度共生矩阵的生长步长；θ是指灰度共生矩阵的生成方向，是指满足d、θ条件的灰度值分别为i，j的点对数，是指满足d、θ条件的所有点对数；共生矩阵P(d,θ)是一个大小为N_q×N_q的方阵：

$P_{d,\mathrm{\θ}}\[i,j\]=\frac{P_{ij}^{d\mathrm{\θ}}}{\underset{i=1}{\overset{N_q}{\Σ}}\underset{j=1}{\overset{N_q}{\Σ}}{P}_{ij}^{d\mathrm{\θ}}}---\left(1\right)$

采用灰度共生矩阵方法提取遥感影像的多个纹理特征，计算各个波段的纹理影像：根据图像大小，选择尺寸为N×N像素的滑动窗口来计算灰度共生矩阵，确定公式(1)步长参数d；计算灰度共生矩阵的方向角θ取0°、45°、90°、135°，从而获得四个方向上的灰度共生矩阵；

选取灰度共生矩阵的4个特征值，计算特征矩阵；根据公式(1)，选取的特征值如下：

①熵

$ENT=-\underset{i}{\Σ}\underset{j}{\Σ}{P}_{d,\mathrm{\θ}}\[i,j\]\lg P_{d,\mathrm{\θ}}\[i,j\]---\left(2\right)$

当P_d,θ[i,j]中的数值均相等时，熵最大，反之，P_d,θ[i,j]的数值之间差别很大时，熵较小；

②角二阶矩

$ASM=\underset{i}{\Σ}\underset{j}{\Σ}{P}_{d,\mathrm{\θ}}^2\[i,j\]---\left(3\right)$

当P_d,θ[i,j]中的数值分布较集中于主对角线附近时，说明从局部区域范围内观察图像的灰度分布是角均匀的，即图像呈现较粗的纹理，相应的角二阶矩值较大，反之较小；

③对比度

$CON=\underset{i}{\Σ}\underset{j}{\Σ}{(i-j)}^2{P}_{d,\mathrm{\θ}}\[i,j\]---\left(4\right)$

④一致性

$HOM=\underset{i}{\Σ}\underset{j}{\Σ}\frac{1}{1+{(i-j)}^2}{P}_{d,\mathrm{\θ}}\[i,j\]---\left(5\right)$

根据得到的灰度共生矩阵计算相应的特征值，作为新的特征矩阵的窗口中心像元；移动纹理窗口，重复上述步骤，直到遍历整个影像，可以得到四个新的纹理矩阵；将得到的四个方向的纹理矩阵取平均值，获得最终的纹理特征矩阵，即纹理影像。