CN102024258B - 具有边界保持特性的遥感图像多尺度分割方法 - Google Patents

Abstract

本发明提供一种针对高空间分辨率遥感图像具有边界保持特性的多尺度图像分割方法。该方法包括以下步骤：采用基于正交滤波器组对图像进行滤波后再求取图像局部能量图的预处理方法，准确定位和检测地物边界；对预处理后的图像进行分水岭分割，得到初始的分割结果；将初始分割结果进行矢量化，生成分割矢量图层，使每个分割区域用多边形矢量来表达，建立矢量多边形之间的邻接关系；对多边形的矢量图层，构建图；基于图的最小生成树算法结合阈值控制进行分割区域的合并；并将图合并过程迭代进行，直到形成所需的尺度分割层次及结构。

Description

具有边界保持特性的遥感图像多尺度分割方法

技术领域

本发明涉及遥感图像处理技术，具体的说，涉及一种具有边界保持特性的多尺度遥感图像分割方法。

背景技术

高空间分辨率遥感图像提供了丰富的关于地物几何结构、纹理细节、地物光谱等方面的信息，在大比例尺遥感制图、环境监测等方面有广阔的应用前景，使得可以在较小的空间尺度上观察地表的细节变化，监测人为活动对环境的影响。

高空间分辨率遥感图像应用的重要途径之一是利用图像分割技术将遥感图像分割成匀质区域，再利用地物的光谱信息、纹理信息、几何结构以及拓扑信息等识别地物。不同地物携带者不同的尺度信息，如湖泊、林地等的尺度不同于耕地、城区道路、针叶林、阔叶林等的尺度，因此如果可将图像进行多尺度的分割，使不同尺度地物的信息提取可选择在相应尺度的分割区域图层上提取，对空间尺度大的类别在分割尺度较大的图层中分析，对于空间尺度较小的类别在分割尺度较小的图层中进行处理，这种多尺度图层提取信息的方法，始终比在一个图层中进行操作更加合理与高效。

传统的遥感图像的分割方法主要包括两大类：先提取边界再获得区域的方法和先获得区域再扫描边界的方法。前者是通常所说的基于边缘的分割方法；后者是通常所说的基于区域的分割方法。基于区域的分割方法以直接寻找区域为基础，包括区域生长法、分裂合并法和基于图像聚类的分割方法等。基于边缘的分割方法计算简单、速度快，但对噪声干扰比较敏感；基于区域的分割方法将区域一致性作为准则，具有一定的抗噪能力，但计算量大，当图像尺寸较大时，分割速度较慢。

分水岭分割方法是目前使用较为广泛的一类新的分割方法，它是基于数学形态学的分割方法，可以和基于边缘的方法结合，也可以和基于区域的方法结合，它将图像梯度幅度看作地理意义上的地形表面，将最高的梯度值的象元看作是分水岭。如此，被分水岭围起的在任何像素位置的水将向下流向共同的局部极小处，因此形成一个“积水盆”。积水盆对应分割区域，分水岭形成闭合的区域边界。分水岭分割算法具有全局分割、边界闭合、实现效率高、准确率高、算法稳定和适用性强等优点，但存在严重的“过分割”问题。

引起分水岭“过分割”现象的主要原因是分水岭变换基于图像梯度幅度进行，易受噪声影响，而且梯度幅度图像中每个局部极小值将对应分割结果中一个单独的区域，在分割结束时会得到待分割对象的区域轮廓，而区域的个数是由局部极小值的数量决定的，图像中由于有纹理等特征的存在，局部极小值的个数会远大于实际的目标对象，这也就造成了大量的虚假轮廓，初始的分水岭分割结果会给出许多的小区域，其中许多区域并不是有意义的分割区域，这就是“过分割”。这些过分割的区域必须合并以形成一个有意义的最终的分割区域。初始的分割区域越多，形成最终分割区域所需要的合并步骤越多，时间的消耗也越长，因此，图像分割前需要预处理的工作。最常使用的预处理方法是低通滤波，如高斯滤波。在高斯滤波后使用梯度运算，继而进行分水岭变换。然而这样一类预处理方法对复合边界会导致一个系统性的定位误差，使用任何的线性滤波器都无法弥补。所谓复合边界指由于成像系统的作用，实际图像中的边界不是简单的阶梯函数，而是阶梯函数、尖峰函数和屋顶函数的组合。

发明内容

本发明的目的是提供一种针对高空间分辨率遥感图像具有边界保持特性的多尺度图像分割方法。

本发明的基本思路为：采用一种能够准确定位复合边界的图像预处理方法，使预处理后获得的图像类似于梯度图像，但图像中的边界位置更准确；对预处理后的图像进行分水岭分割，得到初始的分割结果；将初始分割结果进行矢量化，生成分割矢量图层，使每个分割区域用多边形矢量来表达，建立矢量多边形之间的邻接关系；对多边形的矢量图层，构建图；基于图进行分割区域的合并；将合并过程迭代进行，直到形成所需的尺度分割层次及结构。本发明采用的能够准确定位复合边界的图像预处理方法是基于正交滤波器组对图像进行滤波后再求取图像局部能量图的方法。该方法比先高斯滤波再进行梯度运算的常规预处理方法不仅能极大地减少分水岭过分割的区域，而且地物边界可以得到保持。

即使分水岭分割之前已经过预处理步骤，分水岭初始分割结果仍然包含大量的过分割区域，对过分割区域的合并是必不可少的。常规的基于阈值的合并准则只依赖阈值，随机地从某个区域开始合并操作，开始的区域不同，合并的次序不同，合同的结果也不同。本发明结合图论的最小生成树算法和阈值控制进行合并，使相似的邻接区域优先合并，快速有效地合并掉过大量“过分割”的区域，有效保持地物边界，避免单纯依据阈值不能有效获得匀质区域的问题。本发明还通过在矢量图层中进行过分割区域的合并处理，获得多尺度的分割结果并高效构建多尺度分割结果的层次结构。矢量化首先对初始分割结果进行，使每个分割区域用矢量多边形来表达，然后建立每个多边形之间的邻接关系，生成初始的底层小尺度区域结构，并在此基础上构建图，把初始分割中的区域看成图的顶点，把邻接区域看成是图的边，把邻接区域的特征差异看作是边的权值，结合最小生成树算法进行图合并，生成第二层较大尺度的矢量图层，同时构建该尺度的区域结构。该合并过程可以迭代进行，直到形成所需的尺度分割层次及结构。

本发明的技术方案提供的遥感图像具有边界保持特性的多尺度图像分割方法，其特征在于包括以下实施步骤：

A对输入的遥感图像提取分割波段；

B设计正交滤波器，构造多方向多尺度正交滤波器组，对分割波段滤波后计算图像局部能量图；

C简化图像局部能量图以此作为分水岭分割的输入图像；

D对图像的局部能量图进行分水岭分割，形成初始分割结果；

E对初始分割结果进行矢量化，形成初始分割结果矢量图层；

F扫描分割结果矢量图层，构建图；

G依据图论中的最小生成树算法结合阈值控制进行图合并，使相似的邻接区域优先合并，形成新尺度的分割结果的矢量图层，在此过程中，构建该尺度的区域结构；

H设置迭代次数，通过执行步骤F和步骤G进行迭代合并，产生多个尺度分割结果矢量图层，同时形成多尺度分割结果的区域层次结构。

上述实施步骤的特征在于：

步骤A中所述分割波段，对于遥感的全色图像，直接取该波段；对于多光谱图像，选取其中一个波段或波段组合结果作为分割波段，或在光谱变换后的分量上进行分割；

步骤B中一个正交滤波器包含一个均值为零的偶对称的滤波器和一个奇对称的滤波器，两个滤波器的L²范数相等而且正交，使用由该滤波器组构成的多方向多尺度的滤波器组分别对分割波段进行滤波，进而求得该波段的局部能量图，形成的局部能量图像具有很好的边界保持能力，能够准确定位及检测复合边界，这是与常规高斯滤波结合Sobel、Robert等梯度算子区别最大的地方；

步骤C中所述将图像局部能量图简化的方法是对局部能量图去除掉小的能量值，得到新的局部能量图；

步骤E中所述初始分割结果矢量化，是指初始分割结果中的区域边界以矢量多边形方式表示；

步骤F中所述的图是以初始分割的区域为图的顶点，邻接关系的区域形成图的边，邻接区域间的特征差异为图中连接两个顶点的边的权值；邻接区域间的特征差异通过邻接区域的光谱均值差异、纹理差异、邻接区域公共边界的平均边界强度定义；

步骤G中所述的依据图论中的最小生成树算法结合阈值控制的图合并，其特征是在最小树生成的过程中加入边的权值的阈值进行控制，使得最小树的生成会在中间被打断，成为若干个树，树的个数表示合并后区域的个数，树中的顶点表示合并的区域。所述的图论中的最小生成树算法主要包括Boruvka、Kruskal、Prim三种算法。

步骤H中所述的迭代合并过程是指每一迭代前都要对上一次区域合并的结果构建新的图，然后按照步骤G进行图合并，并将每一次迭代合并后的结果再次生成图，并作为下一次区域合并的输入。

步骤H中所述的区域层次结构包含着不同尺度间区域合并的父子关系，也称为层次关系。其中每个层的结构中包含该层所含的树队列，每个区域结构中包括该层所含的树的序号，每个树中的顶点(子区域)队列、树的父节点区域及树中的边。

本发明与现有技术相比有如下特点：本发明给出的基于局部能量图的用于分水岭分割的预处理方法能够准确定位和检测地物边界，同时有效地减轻分水岭分割的“过分割”；本发明中利用图论的最小生成树算法确定“过分割”区域合并的顺序，并构建了基于矢量图层多尺度分割结果层次化的结构表达，该结构支持不同尺度间分割区域的快速切换以及区域合并父子关系、邻接关系的快速访问、区域特征的快速计算和提取，并使层次间区域边界保持一致。

附图说明：

图1边界保持的遥感图像分割流程图

图2是一个尺度、六个方向奇偶滤波器示意图

图3是一个由分割区域生成的加权无向图示意

图4是区域合并过程的流程图

图5是区域邻接关系结构示意图

图6是区域层次结构示意图

图7是区域节点结构示意图

具体实施方式：

下面结合附图通过实施例对本发明作进一步的详细说明。

本发明所述的边界保持的遥感图像分割方法的实现流程如图1所示，图1包括12个单元。本实施例中的图合并将以最小生成树的Boruvka算法为例说明。

单元100为提取分割波段单元。遥感图像，尤其是多光谱图像，通常具有三个以上的波段。而图像分割往往在一个波段上进行，因此需要提取待分割的波段。对于全色波段，直接在该波段分割即可；对于多光谱波段，选取其中一个波段或在波段组合的图像上分割，也可在光谱变换后的分量上进行分割，如主成分变换后的第一分量或第二分量上进行分割。本实施例中选择在近红外波段上直接进行分割。

单元101为正交滤波器组构造单元。正交滤波器的设计在满足基本的奇、偶、正交性质的同时，还需要考虑这些滤波器有能力能够在一个比较宽的尺度范围内检测特征，特征响应具有很好的方向选择性，计算简单。本发明选择定向的奇对称滤波器为：

$f(x,y)={G^{\′}}_{{\mathrm{\σ}}_1}\left(y\right)G_{{\mathrm{\σ}}_2}\left(x\right)$

其中G_σ(x)表示标准偏差为σ的高斯函数，σ₁∶σ₂表示滤波器拉长的程度；对应的，选择定向的偶对称滤波器为：

$f_1(x,y)={G^{\′\′}}_{{\mathrm{\σ}}_1}\left(y\right)G_{{\mathrm{\σ}}_2}\left(x\right)$

将上述选择的正交滤波器用Hilbert变换表达，即：f₂(x，y)等于f₁(x，y)沿y轴的Hilbert变换，更具体地：

$f_1(x,y)=\frac{d^2}{{\mathrm{dy}}^2}\left(\frac{1}{C}\exp \left(\frac{y^2}{{\mathrm{\σ}}^2}\right)\exp \left(\frac{x^2}{l^2{\mathrm{\σ}}^2}\right)\right)$

f₂(x，y)＝Hilbert(f₁(x，y))

其中σ是尺度，l是滤波器的长宽比，也是滤波器的形状参数，C是规范化参数。由Hilbert变换的性质可知上述滤波器对具有正交滤波器对的性质。

不同尺度、不同方向的滤波器组通过对上述定向的偶和奇对称的滤波器选择不同的σ，旋转指定的角度生成。如图2所示为选用的一个尺度、六个方向的偶对称滤波器和奇对称滤波器。

单元102是图像滤波和计算局部能量图单元。使用单元101得到的正交滤波器组对分割波段进行滤波，然后计算局部能量图。局部能量图的具体计算步骤为：

(1)给出方向个数、尺度数、窗口大小、滤波器长宽比，预先计算偶对称滤波器和奇对称滤波器；本实施中所给的方向个数为6，尺度数为3，窗口大小为15，滤波器长宽比为4。

(2)将图像进行对称扩展，通过快速傅里叶变换进行滤波，基本步骤是先把图像和滤波器分别进行快速傅里叶变换，在频域相乘，然后进行逆傅里叶变换，取实部得到结果；

(3)先计算各个方向上的定向的图像能量，再求和计算局部能量图。

图像在点(x，y)的局部能量为：

$E(x,y)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\sqrt{{O^2}_{i,\mathrm{odd}}(x,y)+{O^2}_{i,\mathrm{even}}(x,y)}$

其中O_i，even(x，y)和O_i，odd(x，y)分别是在点(x，y)图像在第i个方向对应的偶对称滤波器和奇对称滤波器的滤波响应，n为方向个数，本实例中n＝6。

单元103为处理单元102进行滤波后得到的局部能量图。

单元104进行局部能量图简化操作。简化方法为滤除掉小的能量，得到新的局部能量图。局部能量图的简化步骤为：

1)求取局部能量图像中的最大值tMax和最小值tMin；

2)选定滤除系数threshold，求取滤除的能量值tValue，

tValue＝tMin+(tMax-tMin)*threshold，threshold一般取0.02；

3)对局部能量图像进行后处理，使所有能量值小于tValue的能量置为0，其它不变得到新的局部能量图。

单元105对简化后的局部能量图进行分水岭分割，获得初始的分割结果。

单元107对单元106的初始分割结果进行矢量化操作，将初始分割区域以多边形矢量来表示，形成初始分割结果矢量图层单元108。

单元109为构建图的单元。将处理单元106获得的分水岭初始分割的结果表示为图结构，如图3所示，将所有初始分割的区域作为图的顶点，有邻接关系的区域构成图的边，邻接区域间的特征差定义为边的权值。如此构建的图为加权无向图，且图的顶点个数与初始分割区域的个数相同。

图中边的权值，即邻接区域的特征差通过区域之间的光谱距离、纹理距离和邻接区域公共边界的强度来定义：d(R_i，R_j)＝(d₁(R_i，R_j)，d₂(R_i，R_j))^T，其中d(R_i，R_j)为区域R_i和R_j的特征差，d₁(R_i，R_j)＝w×d_spectral(R_i，R_j)+(1-w)×d_texture(R_i，R_j)，

d_spectral(R_i，R_j)是两个区域之间的光谱距离，定义为两个区域均值间的欧式距离；d_texture(R_i，R_j)是两个区域之间的纹理距离，定义为两个区域方差间的欧式距离；d₂(R_i，R_j)是两个邻接区域在局部能量图上公共边界的强度，定义为公共边界的平均强度，这里w为光谱距离和d_spectral(R_i，R_j)和纹理距离d_texture(R_i，R_j)分配的权重，0≤w≤1，一般设置w＝0.5；n为邻接区域公共边界上点的个数，是区域R_i和R_j在局部能量图上公共边界上点(x_k，y_k)的强度。

在后续的处理单元中，约定：区域间特征差的阈值β指阈值向量β＝(β₁，β₂)^T；边的权值小于或大于阈值向量指边的权值的每个分量小于或大于阈值的对应分量，即d(R_i，R_j)≤β等价于d₁≤β₁，d₂≤β₂。

单元110为图合并单元。由于区域合并的次序不同，结果会不同。区域合并可依据图的最小生成树的生成顺序实现区域合并。图的最小生成树算法主要包括Kruskal算法、Prime算法、Boruvka算法。对于图像分割来说，合并成一棵树不是分割的目标，因此利用预先设定的区域之间特征差的阈值控制使分割区域不能最终合并成一个完整的树，这时所剩的树的个数表示合并后区域的个数，每个树中的节点表示合并的区域。

下面给出以Boruvka算法进行图合并的实施例，如图4所示，主要步骤如下：

A单元预先设置图合并的迭代次数和区域间特征差的阈值，以保证最后树的合并会在中间终止，并生成多尺度的合并结果；

B单元将单元109构建的加权无向图保存至一个列表结构中，其中列表的长度是边的数目，以边为序号，列表的每个元是一个结构，该结构中包括四个元素：边的序号、连接边的两个顶点及边的权值。图5所示为列表中元素的结构示意图；

C单元将B单元中的列表按照边的权值的第一或第二分量升序排序，值小的排在前。本发明的实施例中依据边权值的第一分量对边进行排序；

D单元用来生成若干树，每个树中的顶点认为是被合并了的区域。生成树的过程如下：首先初始化一个新树，树中有权值最小的边和连接该边的两个顶点；然后依据边的序号遍历所有的边，判断是否可将该边和连接边的两个顶点添加到树中。依次执行下列树生长的判断：如果边所连接的两个顶点已经在树中，就不添加这个边到树中；如果边的权值大于预先给定的阈值，就不添加这个边到树中；如果边的权值小于预先给定的阈值，则判断连接该边的两个顶点是否已经在树中，如果其中一个已经在树中，就将边也加入树中，并将另一个顶点也添加进树中；如果两个顶点都不在已有的树中，则初始化一个新树；

单元E是树的连接单元。整理D单元生成的树，将不在树中的边按照权值第一分量的大小重新以升序排列，遍历这些边，判断边的权值是否大于预先给定的阈值，如果边的权值大于预先给定的阈值，就不添加这个边到连接边的顶点所在的树中，防止产生闭环，否则添加进树中，使边的两个顶点连接起来，成为一个更大的树；

单元F是整理单元。整理经过E单元处理后生成的树，将每个树作为一个区域，该区域是树中所有顶点所代表的区域合并的结果，更新合并后区域的光谱特征、纹理特征、区域的邻接区域、邻接区域之间的公共边界及边界强度，并将合并后的区域进行层次结构表示，图6是单层(单尺度)合并区域的层次结构。该结构包括：层号(与尺度对应)、树(合并后的区域)的队列；图7是单个树的结构，该结构包括：树中包含的顶点队列(若存在)、树父节点(若存在)。最后将更新后的区域按照单元109的方法重新构建为一个加权无向图，

迭代执行单元A-单元F，直到满足迭代次数的要求；

单元G是迭代完成后的结果，是多层次的合并结果。合并结果具有下列特点：各层合并结果的区域边界与初始分割的结果的边界一致；层次越低，分割结果表示越精细。

单元111将初始分割结果的多层次区域合并的结果作为多尺度分割的结果，并将分割结果表示为区域层次结构。该区域层次结构包含着不同尺度间区域合并的父子关系，也称为层次关系。每个尺度或层次的区域结构包括该层所含的树(区域队列)，每个树(区域)中的顶点(子区域)队列、树的父节点区域及树中的边。该层次结构具有以下特点：不同层次间的分割区域构成了多尺度分割层次结构；上、下层分割区域的边界是一致的；一层的单个区域由它下层所有子区域的总和来表示，每一层又是以它的下一层为基础建立的；不同层的分割区域有合并前后的父子关系。

本发明的一个实施例在PC平台上实现，经实验验证，该分割流程能够得到较理想的分割结果，分割速度快，地物边界保持度高。分割生成的多尺度的区域结果具有良好的适用性，不同层次间分割区域的边界一致。同时，利用基于矢量图层的多尺度的分割层次结构，能在不同尺度间分割区域之间快速切换，能快速访问具有邻接关系和合并前后具有父子关系的区域，并快速计算和提取区域特征。这种多尺度的分割结果和分割结果的层次结构表达为一下步面向对象的遥感图像地物类型精细分类识别奠定了坚实的基础。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。

Claims (13)

1.提供一种针对高空间分辨率遥感图像具有边界保持特性的多尺度图像分割方法，其特征在于包括以下步骤：

A对输入的遥感图像提取分割波段；

B设计正交滤波器，构造多方向多尺度滤波器组，对分割波段滤波后计算图像局部能量图；

C简化图像局部能量图以此作为分水岭分割的输入图像；

D对图像的局部能量图进行分水岭分割，形成初始分割结果；

E对初始分割结果进行矢量化，形成初始分割结果矢量图层；

F扫描分割结果矢量图层，构建图；

G依据图论中的最小生成树算法图合并，形成新尺度分割结果的矢量图层，在此过程中，构建该尺度的区域结构；

H设置迭代次数，通过执行步骤F和步骤G进行迭代合并，产生多个尺度分割结果矢量图层，同时形成多尺度分割结果的区域层次结构。

2.根据权利要求1所述的方法，其特征在于步骤A中所述的分割波段的提取方法，对于遥感的全色图像，直接取该波段；对于多光谱图像，选取其中一个波段或波段组合结果作为分割波段，或在光谱变换后的分量上进行分割。

3.根据权利要求1所述方法，其特征在于步骤B中所述的正交滤波器包含一个均值为零的偶对称的滤波器和一个奇对称的滤波器，两个滤波器的L²范数相等而且正交。

4.根据权利要求1所述的方法，其特征在于步骤B中所述的图像局部能量图，是指利用多方向多尺度正交滤波器组进行图像滤波后求得的。

5.根据权利要求1所述的方法，其特征在于步骤C中所述的简化图像局部能量图的方法，是对局部能量图去除掉小的能量值得到的。

6.根据权利要求1所述的方法，其特征在于步骤E中矢量化指对初始分割结果中的区域边界以矢量多边形方式表示。

7.根据权利要求1所述的方法，其特征在于步骤F中的图的构建以初始分割的区域为图的顶点，有邻接关系的区域形成图的边，邻接区域间的特征差异表示为图中连接两个顶点的边的权值。

8.根据权利要求7所述的方法，其特征在于邻接区域间的特征差异是通过邻接区域的光谱均值差异、纹理差异、邻接区域公共边界的平均边界强度定义的。

9.根据权利要求1所述的步骤G中依据图论中的最小生成树算法的图合并，其特征是在最小树生成的过程中加入边权值的阈值控制，使得最小树的生成会在中间被打断，成为若干个树，树的个数表示合并后区域的个数，树中的顶点表示合并的区域。

10.根据权利要求1所述的步骤G中使用的最小生成树算法包括Boruvka、Kruskal、Prim三种算法。

11.根据权利要求1所述的方法，其特征在于步骤H中所述的区域层次结构包含着不同尺度间区域合并的父子关系，其中每个层的结构中包含该层所含的树队列，每个区域结构中包括该层所含的树的序号，每个树中的顶点队列、树的父节点区域及树中的边。

12.根据权利要求1所述的方法，步骤H中所述的迭代合并过程是指每一迭代前都要对上一次区域合并的结果构建新的图，然后按照步骤G进行图合并，并将每一次迭代合并后的结果再次生成图，并作为下一次区域合并的输入。

13.根据权利要求1所述的方法，形成多尺度分割结果的区域层次结构，其特征在于该层次结构具有以下特点：不同层次间的分割区域构成了多尺度分割层次结构；上、下层分割区域的边界是一致的；一层的单个区域由它下层所有子区域的总和来表示，每一层又是以它的下一层为基础建立的；不同层的分割区域有合并前后的父子关系。

Images