CN105335966B - 基于局域同质性指标的多尺度遥感影像分割方法 - Google Patents

Abstract

本发明公开了基于局域同质性指标的多尺度遥感影像分割方法，依次包括基于二分K均值聚类的影像量化、基于地物上下文信息的多尺度分割以及基于SSIM与D‑S证据理论的区域合并等三个部分。本发明改进了JSEG中粗糙的影像量化，提出了基于地物上下文信息的多尺度分割策略及基于SSIM与D‑S证据理论的区域合并策略，通过对多组不同传感器类型的高分辨率遥感影像进行分割实验，证明了所提出算法能够准确定位对象的边界，同时具有更高的分割精度。

Description

基于局域同质性指标的多尺度遥感影像分割方法

技术领域

本发明属于遥感影像分析技术领域，特别涉及了一种遥感影像分割方法。

背景技术

图像分割技术实现了场景中具有地理意义对象的轮廓信息提取，是利用面向对象的图像分析(OBIA，Object-Based Image Analysis)技术进行遥感影像信息提取与目标识别的前提与基础。与普通图像不同，遥感影像具有多空间分辨率、多光谱分辨率、覆盖范围广泛、地物种类众多以及纹理特征丰富等特点。首先，遥感影像的多波段特性使传统针对单一波段的影像分割方法很难直接应用于多光谱或高光谱遥感影像分割中。另外，遥感影像分割通常会受到地物阴影、云层遮盖等干扰因素的影响，尤其在城市场景中，多样的地物种类及结构复杂的人造目标都给图像分割造成了困难。同时，遥感影像中丰富的纹理特征尽管能够有效描述地物的空间结构特征，同时也对有效抽取及表述对象的纹理特征提出了更高的要求。最后，遥感影像中丰富的纹理与空间语义信息必须在多个尺度下才能获得充分的表述，因此在分割过程中引入多尺度分析工具是提取对象多尺度下空间结构特征的重要手段。这些因素都使遥感影像分割领域在近三十年来始终充满了发展的动力。

目前，针对中、低分辨率遥感影像分割及其应用，学者们已经展开了广泛而深入的研究。例如，Laprade等人利用F检验对局部区域的光谱分布均质程度进行判断，提出了一种基于分裂与合并的分割算法，在航空遥感影像分割中取得了良好的效果；Dong等人针对合成孔径雷达SAR(Synthetic Aperture Radar)影像分割，提出了一种基于高斯—马尔科夫模型GMRF(Gauss-MRF)的分割算法，并与基于Gamma-MRF模型SAR影像分割算法进行了比较，该方法充分利用了相邻光谱间的相似性，显著提高了分割精度；Pan等人利用小波变换多尺度分析工具，采用分水岭分割与基于均值偏移Mean-Shift的聚类方法相结合，尤其适用于光学卫星遥感影像的分割。与此同时，尽管目前针对中、低遥感影像分辨率遥感影像的分割算法众多且算法可靠性、分割精度较高，但针对高分辨率遥感影像的分割研究依然在系统性与针对性中存在诸多不足。

随着遥感技术与计算机技术的不断发展，以SPOT 5、IKONOS、QuickBird等为代表的米级、亚米级高分辨率遥感影像已广泛应用于自然灾害监测与评估、土地资源规划等各个领域。空间分辨率的提高不仅带来了更加丰富的光谱、纹理特征以及形状、上下文等空间信息，同时也造成了更加突出的“同谱异物”及“同物异谱”现象，即不同种类地物的类间可分性降低，而相同种类地物的类内可分性降低。这些因素都给高分辨率遥感影像分割带来了新的挑战。Deng等人提出的JSEG是目前最为流行的彩色纹理分割算法之一，其采用的局域同质性指标J-value对局部区域光谱分布的同质性具有强大的检测能力，并已成功应用于遥感影像分割领域。尽管如此，JSEG算法也存在一些固有局限，如对对象间边界的细节特征不敏感、容易产生过分割等。

为克服这些局限，学者们已经提出了一些改进策略：Zheng等人根据隶属度函数将所有像素分类，利用模糊控制技术对影像量化进行了优化，从而使量化影像更好地保持了原始影像的光谱分布特征，并有效改善了过分割现象；Komati等人提出了三种改进JSEG算法，包括Fractal-JSEG、Fractal-only和Local Fractal Dimension，Fractal-JSEG和Fractal-only方法主要针对对象间边界的细节特征提取，而Local Fractal Dimension主要针对背景区域与前景目标具有相似光谱与纹理特征的应用场合，但这三种改进策略仅适用于普通图像分割。目前尽管JSEG在中、低分辨率遥感影像分割中已取得了良好的效果，但面对空间分辨率的提高给遥感影像分割带来的新挑战，直接采用JSEG算法很难取得满意的效果。

发明内容

为了解决上述背景技术提出的技术问题，本发明旨在提供基于局域同质性指标的多尺度遥感影像分割方法，克服传统JSEG算法在影像量化、多尺度分割及区域合并中存在的局限与不足。

为了实现上述技术目的，本发明的技术方案为：

基于局域同质性指标的多尺度遥感影像分割方法，包括以下步骤：

(1)采用二分K均值聚类方法对影像进行量化；

(2)在量化影像中，计算各尺寸窗口下各个像素对应的局部同质性指标，并将该指标作为该像素的像素值，从而获得多尺度J-image影像序列，再利用地物的上下文信息对多尺度J-image影像序列进行分割；

(3)根据多尺度J-image影像序列中最小尺度影像的分割结果和D-S证据理论，将所有对象与其所有相邻的对象进行合并，得到最终的分割结果。

进一步地，步骤(1)的具体过程如下：

(a)设定量化影像中灰度级总数K＝256；

(b)初始化聚类表

(c)将所有像素作为一个类，加入到聚类表D_T中；

(d)从聚类表D_T中选择一个聚类进行分裂；

(e)使用K均值聚类方法将该聚类划分为两个聚类，遍历所有可能的二分情况，并计算对应的残差平方和SSE；

(f)选择SSE最小的两个聚类，更新聚类表D_T，当聚类个数达到256时，进入步骤(g)，否则返回步骤(c)；

(g)根据聚类表生成隶属关系T，获得调色板D＝[d₁,d₂,...,d₂₅₆]^T，其中d_y表示聚类中心，y∈[1,256]；根据隶属关系及聚类中心进行颜色替换，获得量化影像。

进一步地，步骤(2)的具体过程如下：

计算尺寸为M×M的特定尺寸窗口下所有像素对应的局部同质性指标J-value，并将各像素的J-value作为该像素的像素值，得到该尺度下的J-image，逐步减小M，从而获得一组多尺度J-image影像序列S_t＝(S₁,S₂...,S_T)；

对最大尺度J-image影像S_T进行分块处理，将S_T划分为由尺寸为M_T×M_T像素构成的子图像，其中M_T的值与计算尺度S_T的特定尺寸窗口的尺寸相同；

计算每个子图像的阈值T_J，将每个子图像内部小于其阈值的像素作为种子点，采用4-connectivity方法获得种子区域，其余像素按照J-value值从小到大逐个并入相邻的种子区域，从而获得最大尺度S_T下的分割结果；

利用S_T提取的对象边界将尺度S_T-1分割为由对象集合组成的J-image影像，其中G_T-1为对象总数；在S_T-1中，以每个对象为基本单元，采用与S_T中每个子图像相同的策略提取种子点，并进行区域增长，遍历所有对象，获得尺度S_T-1下的分割结果；利用S_T-1下获得的分割结果，对下一尺度S_T-2采用与S_T-1相同的分割策略，以此类推，直至对最小尺度S₁完成分割。

进一步地，局部同质性指标J-value的计算过程如下：

在量化影像中，将每个像素z对应的坐标z(x,y)作为该像素的像素值；令z(x,y)∈Z，Z为以z为中心，尺寸为M×M的窗口中所有像素的集合，且窗口中的角点被去除；

定义N为以z为中心的窗口中的像素总数，则均值m：

定义m_p为窗口中属于相同灰度级p的所有像素均值，Z_p为窗口中属于灰度级p的所有像素的集合，P为量化影像中的灰度级总数，则窗口中属于同一灰度级像素的方差和S_W定义为：

定义S_A为窗口中所有像素的总体方差：

则局域同质性指标J-value为：

J＝(S_A-S_W)/S_W (4)

进一步地，阈值T_J的计算公式：

T_J＝μ_J+aσ_J (5)

其中，μ_J和σ_J分别为各子图像内像素J-value的均值和标准差，a为预设常数。

进一步地，a＝0.2。

进一步地，步骤(3)的具体过程如下：

(A)将尺度S₁分割结果提取的对象边界映射到所有尺度的J-image中，每个尺度J-image都能获得统一的对象集合R_S＝(R₁,R₂...R_N1)，其中N₁为尺度S₁分割结果中的对象总数；

(B)对于每个尺度J-image的对象集合，确定其中任意一个对象的所有相邻对象，并计算该对象分别与其各个相邻对象之间的相似性SSIM；

(C)根据该对象与其相邻对象之间的相似性SSIM以及D-S证据理论，确定它们是否合并，并将需要合并的对象进行合并；

(D)遍历所有尺度J-image的所有对象，直到没有能够合并的区域为止，从而获得最终分割结果。

进一步地，对象与其相邻对象之间相似性SSIM的计算公式：

其中，A、B分别对象R_A与其相邻对象R_B对应的特征向量；μ_A，μ_B，σ_A，σ_B，σ_AB分别是A与B的均值、标准差、方差以及两者的协方差；C₁，C₂是为了防止当分母接近零时产生不稳定现象所添加的常数。

进一步地，C₁＝0.2，C₂＝0.8。

采用上述技术方案带来的有益效果：

本发明提供的影像分割方法克服了传统JSEG算法在影像量化、多尺度分割及区域合并中存在的局限与不足，提出了基于二分K均值聚类的影像量化策略、基于地物上下文信息的多尺度分割策略以及基于SSIM与D-S证据理论的区域合并策略。

通过对高分辨率多光谱卫星遥感影像与航空遥感影像的实验表明，本发明所采用的影像量化策略能够保持原始影像的光谱与纹理信息；多尺度分割策略能够准确定位对象的边缘，有效提取地物的细节特征，同时能够克服阴影等干扰因素的影响；利用局域同质性指标实现的多尺度区域合并能够有效避免单纯采用光谱特征造成的误合并现象。通过目视分析与精度评价表明，与传统JSEG算法及著名商业软件eCognition相比，所提出算法具有更好的分割效果及分割精度，同时避免了eCognition分割结果中狭长相邻对象间通常存在的虚假单元现象。

附图说明

图1为实施例中M＝9特定尺寸窗口示意图；

图2为实施例中上海地区SPOT 5遥感影像图；

图3为实施例中量化影像17灰度级图；

图4为本发明的流程图；

图5为实施例中实验一的量化影像图；

图6为实施例中实验一的区域合并前影像图；

图7为实施例中实验一的区域合并后影像图；

图8为实施例中实验一的JSEG算法分割结果影像图；

图9为实施例中实验一的eCognition分割结果影像图；

图10为实施例中河海大学航空遥感DOM影像图；

图11为实施例中实验二的量化影像；

图12为实施例中实验二的区域合并前；

图13为实施例中实验二的区域合并后；

图14为实施例中实验二的JSEG算法分割结果；

图15为实施例中实验二的eCognition分割结果；

图16为实验一的精度评价示意图；

图17为实验二精度评价示意图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

本实施例首先介绍了JSEG算法的基本原理与实现流程，并分析了JSEG算法在高分辨率遥感影像分割中存在的关键问题与不足，进而阐述了本发明所提出算法的改进策略与实现流程；并对实验结果进行了分析和比较。

JSEG算法的基本原理

传统JSEG算法的分割过程主要包括影像量化、空间分割及区域合并三个步骤。在影像量化中，首先将原始影像转换到LUV彩色空间，进而采用对等滤波器组PGF(Peer GroupFiltering)对影像进行平滑，然后采用经典的硬阈值HCM(Hard C-Means)方法进行颜色量化，从而获得压缩后的单波段量化影像。

在空间分割中，首先利用局部同质性指标J-value计算多尺度J-image影像序列。J-value的计算过程如下：

在量化影像中，将每个像素z对应的坐标z(x,y)作为该像素的像素值。令z(x,y)∈Z，Z为以z为中心，尺寸为M×M像素窗口中所有像素的集合。为保证各个方向的一致性，窗口中的角点被去除。以M＝9为例，像素z对应的窗口如图1所示：

在量化影像中，定义N为以z为中心窗口中的像素总数，则均值m：

定义m_p为窗口中属于相同灰度级p的所有像素均值，Z_p为窗口中属于灰度级p的所有像素的集合，P为量化影像中的灰度级总数，则窗口中属于同一灰度级像素的方差的和S_W可定义为：

定义S_A为窗口中所有像素的总体方差：

则局域同质性指标J-value为：

J＝(S_A-S_W)/S_W (4)

将像素z对应的J-value值作为该像素的像素值，遍历所有像素，可获得M×M像素尺度下的J-image。通过改变M的值，可获得多尺度下的J-image影像序列。M值越大，所获得的J-image越粗糙，尺度越大。通过公式(1)～(4)可以看出，J-value反映了以像素z为中心的特定尺寸区域内光谱分布的均质程度，J-value越大，则像素z越有可能位于对象的边缘。反之，z可能位于对象的中心。基于这一特性，在某一尺度J-image中，首先根据阈值T_J确定种子区域，分割出影像中对象的中心区域，T_J的定义如下：

T_J＝μ_J+aσ_J (5)

其中，μ_J和σ_J分别为像素z所在窗口内J-value的均值和标准差。a为预先设定的阈值，通常取a∈[-0.4,0.4]。将小于阈值T_J的点作为种子点，进而采用4-connectivity获得种子区域。计算非种子区域中像素的J-value均值，所有小于均值的像素构成增长区域。若该增长区域仅与一个种子区域相邻，则并入该种子区域。利用下一尺度J-image对剩余像素重新计算增长区域，直到最小尺度。对剩余像素按照J-value值从小到大逐个并入相邻的种子区域。在最后的区域合并中，利用区域颜色相似性信息并基于欧氏距离最小准则对相邻区域进行合并，获得最终分割结果。

JSEG算法存在的关键问题

根据高分辨率遥感影像的特点，通过分析JSEG算法的基本原理及相关文献，总结JSEG算法在高分辨率遥感影像分割中存在的关键问题主要包括：

(1)“影像量化”问题

JSEG算法首先进行影像量化的目的是实现影像平滑以及灰度级压缩，从而减少计算量。但量化后的影像通常仅包含几个或者十几个灰度级，而粗糙的量化会造成高分辨率遥感影像中大量细节信息的严重损失，也降低了不同种类地物间的类间可分性，从而难以准确提取对象的边缘信息。以2004年采集的中国上海地区SPOT 5全色—多光谱融合影像为例，尺寸为512×512像素，如图2所示。采用JSEG量化后影像仅包含17个灰度级，如图3所示。

对比图2、图3可以看出，量化影像中河流、道路及各种结构复杂的人造目标细节特征损失严重，且植被覆盖区域与周边道路等目标具有相似的光谱特征。因此，为了影像平滑及灰度级压缩而进行过于粗糙的影像量化对于高分辨率遥感影像分割已变得没有意义。

(2)“种子点提取及分布约束”问题

基于公式(5)，JSEG采用统一的阈值在最大尺度J-image影像中提取种子点，进而生成种子区域。种子区域一旦确定，在后续的多尺度分割中不再生成新的种子区域，非种子区域像素仅作为增长区域逐个并入现有种子区域。这样的种子点提取策略难以准确标记影像中的各种地理对象，主要体现在以下两个方面。首先，由于最大尺度J-image影像也最为粗糙，即每个像素的J-value反映了更大范围内的光谱分布同质性。尽管大尺度J-image能够有效平滑噪声等干扰因素的影响，却难以准确反映场景中的小目标及细节特征，因此容易造成欠分割现象。另一方面，高分辨率遥感影像中不同局部区域的纹理复杂程度通常差异较大，纹理特征复杂的区域应当提取更多的种子点以标记图像的细节特征，而均质程度较高的区域仅提取少量种子点即可。JSEG由于采用了统一的阈值进行种子点提取，较大的阈值尽管有利于提取图像的细节特征，但容易在均质程度较高的区域造成严重的过分割，反之则容易在纹理复杂的局部区域产生欠分割现象。

(3)“区域误合并”问题

高分辨率遥感影像中通常地物种类众多且纹理特征丰富，而普遍存在的“同谱异物”现象使相邻不同种类地物间可能具有相似的光谱特征。JSEG算法在区域合并时采用R.O.Duda等人提出的相邻对象灰度直方图间的欧氏距离作为相似性度量，仅考虑了对象内部的光谱特征，而忽略了纹理特征、尺度信息等可以进步一描述对象间差异的有用信息。因此，这样的区域合并策略难以克服“同谱异物”现象造成的干扰，从而容易产生错误的区域合并结果。

本发明的解决对策与算法实现流程

根据局域同质性行指标J-value的定义，J-value综合反映了原始影像所包含的光谱与纹理特征，同时还有具有多尺度特性。因此，利用J-value对对象的轮廓进行描述有利于更加深刻的分析对象内部的空间结构信息，同时有利于减少孤立点、噪声等干扰因素的影响。另一方面，上述诸多关键问题使JSEG算法在高分辨率遥感影像分割中难以取得理想的效果。针对这一问题，本发明提出了一种基于局域同质性指标的多尺度高分辨率遥感影像分割方法，主要包括基于二分K均值聚类的影像量化、基于地物上下文信息的多尺度分割及基于SSIM与D-S证据理论的区域合并三部分组成。算法实现流程如图4所示：

二分K均值聚类影像量化

为减少计算量及平滑影像，在区域分割之前首先对高分辨率遥感影像进行量化。采用二分K均值聚类方法进行影像量化，获得灰度级为256的8bit灰度影像，以尽量保持原始影像中的细节特征。所采用的具体量化策略如下：

假设原始影像的尺寸为P×Q像素，波段总数为C，用矩阵表示影像I为I＝[x₁,x₂,...,x_P·Q]^T。其中，x_i＝[c_1i,...,c_Ci]为彩色空间中的任意像素点，c代表原始影像的某一波段，设定量化影像中灰度级总数K＝256，输出调色板为D＝[d₁,d₂,...,d₂₅₆]^T，d_y表示聚类中心，y∈[1,256]。

Step1：初始化聚类中心

Step2：将所有像素作为一个类，加入到类表中：D_T＝{cluster(1)}；

Step3：从聚类表中选择一个聚类cluster(y)进行分裂；

Step4：开始循环。使用K均值聚类方法将该类划分为两个聚类，并计算归残差平方和SSE(Sum of Squares for Error)；

Step5：退出循环。选择SSE最小的两个聚类，更新聚类表；

Step6：当聚类个数达到256时，根据聚类表生成隶属关系T，获得调色板D，根据隶属关系及聚类中心进行颜色替换，获得量化影像。否则，返回Step3。

基于地物上下文信息的多尺度分割

在获得8bit量化影像后，为避免前文提出的JSEG算法“种子点提取及分布”问题，提出了一种基于地物上下文信息的多尺度分割策略。

首先通过逐渐减小特定窗口的尺寸M，可利用公式(4)计算量化影像单一尺度下的J-image影像S_t，从而获得一组多尺度J-image影像序列S_t＝(S₁,S₂...,S_T)，定义原始影像为S₀。在此基础上，对最大尺度J-image影像S_T进行分块处理，将S_T划分为由一个个尺寸为M_T×M_T像素构成的子图像。其中M_T的值与计算尺度S_T的特定尺寸窗口尺寸相同。由于J-image中每个像素z的像素值都是该像素对应的局域同质性指标J-value，因此计算每个子图像内的像素均值能够反映该子图像内光谱分布的均质程度。即子图像内部的灰度均值越大，该子图像内部的光谱特征越丰富，纹理特征越复杂，可采用更大阈值以提取更多的图像细节特征。反之，采用较小的阈值即可。因此，本发明改变了JSEG算法中对整幅影像采用统一阈值的种子点提取策略，而对每个子图像利用公式(5)单独计算阈值T_Jl(l＝1,2....L)，L为子图像总数，参数a根据经验统一设定为0.2。根据J-value的定义，纹理复杂的子图像通常具有更高的均值μ_Jl和方差σ_Jl，因此T_Jl更大。将每个子图像内部小于阈值T_Jl的点作为种子点，进而采用4-connectivity方法获得种子区域，其余像素按照J-value值从小到大逐个并入相邻的种子区域，从而获得尺度S_T下的分割结果。

由于J-image影像序列中单一尺度的J-image尺寸均与原始影像相同，因此可直接利用S_T提取的对象边界将尺度S_T-1分割为由对象集合组成的J-image影像，其中G_T-1为对象总数。在S_T-1中，以每个对象R_i为基本单元，采用与S_T中每个子图像相同的策略提取种子点，并进行区域增长。遍历所有对象R_i，获得尺度S_T-1下的分割结果。利用S_T-1下获得的分割结果，对下一尺度采用与S_T-1相同的分割策略，直到最小尺度S₁(即最精细尺度)计算完毕，从而利用地物上下文信息实现了由粗到精的多尺度分割结果。同时可以看出，分割过程中种子点也实现了由粗到精的提取过程，同时保证了种子点可根据局部区域内部的光谱分布均质程度不同而合理的分布。

基于SSIM与D-S证据理论的区域合并

在获得初步分割结果后，针对前文提出的JSEG“区域误合并”问题，提出了一种基于SSIM与D-S证据理论的区域合并策略。由于局域同质性指标J-value综合反映了原始影像的光谱、纹理特征，并具有尺度特性，因此在J-image中采用J-value描述对象间的差异较单纯采用光谱特征具有更高的可靠性。另一反面，需要选择合适的相似性度量进行相邻对象间的相似性描述。

首先将尺度S₁分割结果提取的对象边界映射到所有尺度的J-image中，因此在每个尺度J-image都能够获得同一的对象集合R_S＝(R₁,R₂...R_N)，其中N为尺度S₁分割结果中的对象总数。在某一尺度J-image中，计算任意两个相邻对象R_A和R_B的结构相似性如下：

其中，A、B分别R_A和R_B对应的特征向量；μ_A，μ_B，σ_A，σ_B， σ_AB分别是A与B的均值、标准差、方差和协方差；C₁，C₂是为了防止当分母接近零时产生不稳定现象所添加的常数，通常可设定为C₁＝0.2，C₂＝0.8。S(A,B)越大，则R_A和R_B相似度越高，即有可能属于同一对象。

为了进一步提高区域合并的可靠性，本文对度尺度下获得的相邻对象间相似性描述进行多尺度决策级融合，并提出了一种基于D-S证据理论的区域合并规则。D-S证据理论的相关概念如下：

定义U为识别框架，在U上的基本概率分配BPAF(Basic Probability AssignmentFormula)是一个函数m:2^U→[0,1]，且m满足

其中，若A满足m(A)＞0，A称为一个焦元(Focal Elements)，m(A)表示证据对A的一种信任度量。Dempster合成规则(Dempster’s Combinational Rule)定义如下：

对于U上的n个mass函数m_i(i＝1,2,...,n)的Dempster合成法则为：

其中，K为归一化常数，其反映了证据的冲突程度，定义如下：

定义U:{Y,N}。其中Y代表合并类，N代表非合并类。因此2^U的非空子集包括：{Y}，{N}，{Y,N}。定义尺度S_t中任意两个相邻R_A和R_B间的SSIM为SSIM_ABt，建立如下的BPAF：

m_ABt({Y})＝SSIM_ABt×α_t (11)

m_ABt({N})＝(1-SSIM_ABt)×α_t (12)

m_ABt({Y,N})＝1-α_t (13)

其中，α_t∈(0,1)代表了对尺度S_t相似性判别结果的信任度。较大尺度的J-image能有效减少噪声及孤立点的影响，而较小的尺度有利于反映对象内部的细节特征。因此，所涉及的阈值及参数可根据具体应用场合及经验人工设定。定义区域合并规则如下：

Step1：根据尺度S₁获得的分割结果生成区域邻接图RAG(Region AdjacencyGraphics)；

Step2：对任意一个对象R_A，根据RAG确定所有与R_A相邻的对象；

Step3：计算R_A与相邻的一个对象R_B，计算所有尺度下R_A和R_B的SSIM，即SSIM_AB＝{SSIM_AB1,SSIM_AB2...,SSIM_ABT}；

Step4：利用公式(9)计算m_ABt({Y})、m_ABt({N})。若m_ABt({Y})＞0.7且m_ABt({N})＜0.4，则合并区域R_A和R_B，进如下一步。否则，返回Step2；

Step5：更新RAG，重复Step2至Step4，遍历所有对象及其相邻对象，直到没有能够合并的区域为止，从而获得最终分割结果。

实验与分析

为验证算法的有效性与可靠性，选择两组不同空间分辨率、不同传感器类型的高分辨率遥感影像进行实验，并将实验结果分别与传统JSEG算法以及著名商业软件eCognition进行了比较。

德国Definiens Imaging公司开发的eCognition是目前公认最好的面向高分辨率遥感影像分类的商业软件，其通过模拟人类的感知系统，综合考虑了对象的光谱特征形状特征，实现了基于分型网络演化(Fractal Net Evolution Approach，FNEA)方法的多分辨率分割。eCognition的分割过程主要有三个参数需要进行人工设定，包括：尺度参数，尺度参数决定了分割结果中最小对象的尺寸，参数值越大，分割结果中对象的平均尺寸越大；形状参数，用于保持分割结果中对象的轮廓完整性；紧致度参数，用于提高具有相似光谱特征的对象的类间可分性。本文在实验中通过改变以上三个参数，选择分割效果最佳的参数组合用于比较实验。

实验一结果与目视分析

实验一采用2003年采集的中国上海地区SPOT 5全色—多光谱融合影像，空间分辨率为2.5m。多光谱影像包括红色波段、绿色波段及近红外波段，空间分辨率为10m。全色影像的空间分辨率2.5m，融合影像的尺寸为512×512像素示。实验中，设定JSEG所采用的J-image影像序列特定窗口尺寸为20×20像素、10×10像素、5×5像素，参数a＝0.2。本发明算法J-image影像序列特定窗口尺寸采用相同的设定。在公式(11)～(13)中，设定α₁＝0.95，α₂＝0.9，α₃＝0.8，即越精细尺度的J-image中检测结果的可信度越高。eCongition分割中设定尺度参数为75，形状参数为0.5，紧致度参数为0.5。本发明算法提取的量化影像如图5所示：

通过比较图3及图5可以看出，与JSEG算法的粗糙影像量化相比，本发明采用的影像量化策略能够有效保持图像的细节信息，量化影像能够充分反映场景中河流、道路以及其他结构负责人造建筑的轮廓细节特征，因此有更利于准确定位对象的边缘。三种方法实验结果如图6～图9所示。为了便于对场景中的典型地物进行目视分析，对其所在位置进行了标号，如图中位置所示。

对比图6、图7可以看出，本发明算法提出的区域合并策略使过分割现象得到了明显的改善，如位置且基本没有发生误合并现象。对比图7～9可以看出，本发明算法及eCognition分割精度明显高于传统JSEG算法，如对于位置中形状规则的人造建筑，JSEG算法仅能够提取对象的大致轮廓，而无法准确定位对象的边缘。进一步对比本发明算法及eCognition分割结果可以看出，对于场景中结构复杂的人造目标，本发明算法能够提取更加丰富的细节轮廓特征，例如位置但依然存在个别过分割区域，如位置而eCognition提取的对象轮廓更为完整，如位置的河流区域，位置的道路区域等，但eCognition分割结果中在狭长的相邻对象间存在明显的虚假单元，如位置等。

实验二结果与目视分析

第二组实验采用空间分辨率为0.5m的高分辨率航空遥感DOM(DigitalOrthophoto Map)影像，采集时间为2009年3月，所在地区为中国江苏省南京市河海大学江宁校区，影像尺寸为512×512像素。原始影像与量化影像如图10、图11所示。

通过图10可以看出，与第一组实验相比，更高的空间分辨率提供了更加丰富的对象细节特征，因此相邻对象间的边界也更加明显，影像中存在的干扰因素主要为大量地表高层建筑造成的阴影区域。实验中，设定本发明算法与JSEG算法计算J-image影像序列的特定窗口尺寸均为10×10像素、7×7像素、5×5像素，eCongition分割中尺度参数为115，其他参数设定均与第一组实验相同。三种算法实验结果如图12～图15所示。

为便于目视分析，与第一组实验相同对典型地物所在位置进行了标记，如图中位置所示。通过对比图13、图14可以看出，本发明算法分割结果在区域合并后过分割现象得到了明显的改善，如位置的操场草坪及位置的道路区域。对比图14～16可以看出，尽管三种算法都有效识别了场景中的阴影区域，如位置b，但对于处于阴影区域中的对象提取，本发明明显优于其他两种算法，如位置JSEG算法提取的操场区域的轮廓最为完整，如位置但JSEG算法没有提取出操场中的球门区域这一细节信息，如位置a，而本发明算法和eCognition都有效提取了该区域。对于场景中结构复杂的人造目标如位置 本发明算法和eCognition的分割效果明显优于JSEG算法。与第一组实验相同，eCognition分割结果中在相邻的狭长对象间依然存在明显的虚假单元如，位置 同时对阴影区域中的对象分割效果不佳，如位置

精度评价

上文仅对两组实验结果进行了目视分析，下文进一步对不同算法的分割结果进行定量精度评价。选择Deng等人提出的遥感影像分割精度评价标准进行对两组实验进行定量分析^[10]：

首先利用公式(4)计算分割结果中每个对象对应的局域同质性指标J-value的值J_v，进而利用公式(14)计算精度评价指标

其中，Q为影像中的像素总数，V为分割结果中的对象总数，W_v为对象R_v内部的像素总数。当算法对应的越大时，说明分割结果中区域的平均内部均质程度越高，分割效果越好。基于这一准则，首先对不同算法中的J_v分布曲线进行分析，如图16、17所示。(曲线中每个区域对应的J_t按照就近原则在[0,1]区间内被量化为20个单位，纵坐标为J_t，横坐标为当前J_t值所占的比例，虚线为不同算法对应的精度指标)。

如图16、17所示，不同样式的曲线分别代表了不同算法中J_v值的分布情况。两组实验中本发明所提出算法提取的对象更集中于J_v较小的区域，因此平均均质程度最高即J最小，因此分割精度最高。JSEG算法分割精度则明显低于其他两种算法，与目视结果分析一致。两组实验中三种算法精度指标J评价结果如表1所示：

表1

通过表1可以看出，本发明提出的算法在区域合并后分割精度显著提高，与目视分析结果一致，因此，本发明所提出的区域合并策略是有效的。另一方面，与实验一相比，实验二中三种算法的分割精度都显著提高，由于实验二采用的航空遥感影像空间分辨率更高，因此对象的轮廓更加明显，同时场景中存在大面积的均质区域如操场、屋顶等所造成的。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (8)

1.基于局域同质性指标的多尺度遥感影像分割方法，其特征在于，包括以下步骤：

(1)采用二分K均值聚类方法对影像进行量化，具体过程如下：

(a)设定量化影像中灰度级总数K＝256；

(b)初始化聚类表

(c)将所有像素作为一个类，加入到聚类表D_T中；

(d)从聚类表D_T中选择一个聚类进行分裂；

(e)使用K均值聚类方法将该聚类划分为两个聚类，遍历所有可能的二分情况，并计算对应的残差平方和SSE；

(f)选择SSE最小的两个聚类，更新聚类表D_T，当聚类个数达到256时，进入步骤(g)，否则返回步骤(c)；

(g)根据聚类表生成隶属关系T，获得调色板D＝[d₁,d₂,...,d₂₅₆]^T，其中d_y表示聚类中心，y∈[1,256]；根据隶属关系及聚类中心进行颜色替换，获得量化影像；

(2)在量化影像中，计算各尺寸窗口下各个像素对应的局部同质性指标，并将该指标作为该像素的像素值，从而获得多尺度J-image影像序列，再利用地物的上下文信息对多尺度J-image影像序列进行分割；

(3)根据多尺度J-image影像序列中最小尺度影像的分割结果和D-S证据理论，将所有对象与其所有相邻的对象进行合并，得到最终的分割结果。

2.根据权利要求1所述基于局域同质性指标的多尺度遥感影像分割方法，其特征在于，步骤(2)的具体过程如下：

计算尺寸为M×M的特定尺寸窗口下所有像素对应的局部同质性指标J-value，并将各像素的J-value作为该像素的像素值，得到该尺度下的J-image，逐步减小M，从而获得一组多尺度J-image影像序列S_t＝(S₁,S₂...,S_T)；

对最大尺度J-image影像S_T进行分块处理，将S_T划分为由尺寸为M_T×M_T像素构成的子图像，其中M_T的值与计算尺度S_T的特定尺寸窗口的尺寸相同；

计算每个子图像的阈值T_J，将每个子图像内部小于其阈值的像素作为种子点，采用4-connectivity方法获得种子区域，其余像素按照J-value值从小到大逐个并入相邻的种子区域，从而获得最大尺度S_T下的分割结果；

利用S_T提取的对象边界将尺度S_T-1分割为由对象集合组成的J-image影像，其中G_T-1为对象总数；在S_T-1中，以每个对象为基本单元，采用与S_T中每个子图像相同的策略提取种子点，并进行区域增长，遍历所有对象，获得尺度S_T-1下的分割结果；利用S_T-1下获得的分割结果，对下一尺度S_T-2采用与S_T-1相同的分割策略，以此类推，直至对最小尺度S₁完成分割。

3.根据权利要求2中基于局域同质性指标的多尺度遥感影像分割方法，其特征在于，局部同质性指标J-value的计算过程如下：

在量化影像中，将每个像素z对应的坐标z(x,y)作为该像素的像素值；令z(x,y)∈Z，Z为以z为中心，尺寸为M×M的窗口中所有像素的集合，且窗口中的角点被去除；

定义N为以z为中心的窗口中的像素总数，则均值m：

<mrow> <mi>m</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>z</mi> <mo>&amp;Element;</mo> <mi>Z</mi> </mrow> </munder> <mi>z</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

定义m_p为窗口中属于相同灰度级p的所有像素均值，Z_p为窗口中属于灰度级p的所有像素的集合，P为量化影像中的灰度级总数，则窗口中属于同一灰度级像素的方差和S_W定义为：

<mrow> <msub> <mi>S</mi> <mi>W</mi> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>P</mi> </munderover> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>z</mi> <mo>&amp;Element;</mo> <msub> <mi>Z</mi> <mi>p</mi> </msub> </mrow> </munder> <mo>|</mo> <mo>|</mo> <mi>z</mi> <mo>-</mo> <msub> <mi>m</mi> <mi>p</mi> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

定义S_A为窗口中所有像素的总体方差：

<mrow> <msub> <mi>S</mi> <mi>A</mi> </msub> <mo>=</mo> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>z</mi> <mo>&amp;Element;</mo> <mi>Z</mi> </mrow> </munder> <mo>|</mo> <mo>|</mo> <mi>z</mi> <mo>-</mo> <mi>m</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

则局域同质性指标J-value为：

J＝(S_A-S_W)/S_W (4)。

4.根据权利要求3所述基于局域同质性指标的多尺度遥感影像分割方法，其特征在于，阈值T_J的计算公式：

T_J＝μ_J+aσ_J (5)

其中，μ_J和σ_J分别为各子图像内像素J-value的均值和标准差，a为预设常数。

5.根据权利要求4所述基于局域同质性指标的多尺度遥感影像分割方法，其特征在于：a＝0.2。

6.根据权利要求2所述基于局域同质性指标的多尺度遥感影像分割方法，其特征在于，步骤(3)的具体过程如下：

(A)将尺度S₁分割结果提取的对象边界映射到所有尺度的J-image中，每个尺度J-image都能获得统一的对象集合其中N₁为尺度S₁分割结果中的对象总数；

(B)对于每个尺度J-image的对象集合，确定其中任意一个对象的所有相邻对象，并计算该对象分别与其各个相邻对象之间的相似性SSIM；

(C)根据该对象与其相邻对象之间的相似性SSIM以及D-S证据理论，确定它们是否合并，并将需要合并的对象进行合并；

(D)遍历所有尺度J-image的所有对象，直到没有能够合并的区域为止，从而获得最终分割结果。

7.根据权利要求6所述基于局域同质性指标的多尺度遥感影像分割方法，其特征在于，对象与其相邻对象之间相似性SSIM的计算公式：

<mrow> <mi>S</mi> <mi>S</mi> <mi>I</mi> <mi>M</mi> <mrow> <mo>(</mo> <mi>A</mi> <mo>,</mo> <mi>B</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>(</mo> <mn>2</mn> <msub> <mi>&amp;mu;</mi> <mi>A</mi> </msub> <msub> <mi>&amp;mu;</mi> <mi>B</mi> </msub> <mo>+</mo> <mo>)</mo> <mo>(</mo> <mn>2</mn> <msub> <mi>&amp;sigma;</mi> <mrow> <mi>A</mi> <mi>B</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>C</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msubsup> <mi>&amp;mu;</mi> <mi>A</mi> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;mu;</mi> <mi>B</mi> <mn>2</mn> </msubsup> <mo>+</mo> <msub> <mi>C</mi> <mn>1</mn> </msub> <mo>)</mo> <mo>(</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>A</mi> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>B</mi> <mn>2</mn> </msubsup> <mo>+</mo> <msub> <mi>C</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

其中，A、B分别对象R_A与其相邻对象R_B对应的特征向量；μ_A，μ_B，σ_A，σ_B，σ_AB分别是A与B的均值、标准差、方差以及两者的协方差；C₁，C₂是为了防止当分母接近零时产生不稳定现象所添加的常数。

8.根据权利要求7所述基于局域同质性指标的多尺度遥感影像分割方法，其特征在于：C₁＝0.2，C₂＝0.8。