基于可变尺度Mean-Shift的农田高分遥感影像分割算法

引用本文

苏腾飞, 张圣微, 李洪玉. 基于可变尺度Mean-Shift的农田高分遥感影像分割算法[J]. 国土资源遥感, 2017,29(3): 41-50
SU Tengfei, ZHANG Shengwei, LI Hongyu. Variable scale Mean-Shift based method for cropland segmentation from high spatial resolution remote sensing images[J]. REMOTE SENSING FOR LAND & RESOURCES,2017,29(3): 41-50 复制到剪切板

doi: 10.6046/gtzyyg.2017.03.06
Permissions

《国土资源遥感》编辑部

基于可变尺度Mean-Shift的农田高分遥感影像分割算法

苏腾飞, 张圣微, 李洪玉

内蒙古农业大学水利与土木建筑工程学院,呼和浩特 010018

通信作者：张圣微(1979-),男,博士,教授,主要从事定量遥感、生态水文及气候变化等方面的研究。Email:zsw_imau@163.com。

第一作者: 苏腾飞(1987-),男,硕士,实验师,主要从事面向对象的遥感图像分析算法方面的研究。Email:stf1987@126.com。

收稿日期: 2015-08-17

基金项目: 国家自然科学基金项目“科尔沁沙地典型生态系统水热通量传输机理及其与植被耦合关系试验和模拟研究”(编号: 51569017)、“内蒙古典型草原水文过程及其扰动与触发草地退化的水文临界条件实验与模拟研究”(编号: 51269014)、内蒙古自然科学基金项目“半干旱区沙地典型生态系统水热通量传输机理研究”(编号: 2015MS0514)和中国博士后科学基金面上资助项目“西部地区博士后人才资助计划”(编号: 2015M572630XB)共同资助

摘要

为了提升农田高分遥感影像(high spatial resolution remote sensing image,HRI)的信息提取效果,提出了一种新的农田HRI分割算法。传统的Mean-Shift(MS)HRI分割算法仅利用全局或单一的尺度参数; 而常规可变尺度MS算法在尺度参数估算中也只考虑光谱信息。这些都导致其分割结果难以完整地展现不同尺度的农田区域。针对该问题,在MS算法的基础上进行了改进: 第一,提出了一种局部可变尺度参数的估计方法; 第二,提出了利用局部可变尺度进行MS滤波的模型函数。该改进算法主要包含3步: ①为了全面考虑不同波段的响应变化,在MS滤波核函数中采用了对角化的尺度参数矩阵,并将其与采样点密度估计模型相结合,导出了一种可变尺度MS滤波的迭代函数; ②为了提高算法的自动化程度,利用局部光谱变化与边界强度信息,提出了一种新的局部尺度参数估算方法; ③将MS滤波结果输入到基于分形网络演化方法(fractal net evolution approach,FNEA)的空间聚类算法中,得到最终的分割结果。利用RapidEye与OrbView3的2景HRI进行了算法验证。实验结果表明,所提出的改进算法能够优化农田HRI分割的精度。

关键词: 可变尺度; Mean-Shift; 农田分割; 高分遥感影像

文献标志码:A 文章编号:1001-070X(2017)03-0041-10

Variable scale Mean-Shift based method for cropland segmentation from high spatial resolution remote sensing images

SU Tengfei, ZHANG Shengwei, LI Hongyu

College of Water Conservancy and Civil Engineering, Inner Mongolian Agricultural University, Hohhot 010018, China

Abstract

In order to improve the effect of information extraction from high spatial resolution remote sensing images (HRI) of cropland, the authors put forward a new HRI segmentation algorithm. Due to the fact that the traditional Mean-Shift (MS) segmentation method only uses a global and single scale, and that some variable bandwidth MS only considers spectral information in their scale estimation process, and croplands with various sizes could be hardly extracted in one segmentation result, the authors improved a MS based approach to tackle this problem. The main consideration lies in two aspects: ① A local variable scale parameter estimation method is proposed; ② The model function for local variable scale is established for MS filtering. The proposed approach mainly consists of 3 parts: ① With the objective of comprehensively considering the response variation of different bands, the diagonal scale parameter matrix is adopted in the kernel function of MS filtering, and it is combined with sample point estimation model to derive the iterative function for variable scale MS filtering; ② For the purpose of increasing automation of the proposed method, local spectral variation and edge strength information are utilized to design a new local scale parameter estimation method; ③ For obtaining the final segmentation, the filtering result is used as input for the fractal net evolution approach (FNEA) which is a spatial clustering method. Two scenes of HRI acquired by RapidEye and OrbView3 were employed for experiment, and the results show that the proposed method can optimize the accuracy of cropland HRI segmentation.

Keyword: variable scale; Mean-Shift; cropland segmentation; high spatial resolution remote sensing image

Show Figures

0 引言

近年来, 一种新兴的遥感影像解译方式— — 面向地理对象的影像分析(geographical object based image analysis, GEOBIA)^[1]在农田遥感分类研究中逐渐增多^{[2, 3, 4]}, 这主要归因于2个方面: ①农田遥感监测需要高空间分辨率遥感影像(high spatial resolution remote sensing image, HRI), 以准确提取作物信息。很多地区(如中国)的农田面积较小^[5], 导致中低空间分辨率的遥感影像易产生混合像素光谱, HRI可有效地解决这一问题, 很多研究都表明GEOBIA适用于HRI的解译^{[1, 2]}; ②在GEOBIA中, 除农田的光谱特征外, 形状和背景特征也能被有机地融合到作物种类判别算法中, 从而显著提高了作物识别精度。例如一些地区使用了圈形喷灌机, 使农田呈规则的圆形或扇形, 这可被GEOBIA用来区分农田和其他地物。可见, 利用GEOBIA和HRI开展农田分类研究具有重要意义。

GEOBIA首先就是影像分割, 其目的是获取地物的原始图斑; 然后通过特征选择和影像分类, 赋予原始图斑以真实地物的语义信息。不准确的分割结果会导致后续分类的严重误判, 因此影像分割是GEOBIA的关键步骤。虽然目前已有大量的影像分割模型与软件, 但如何提高农田HRI的分割精度依然是具有挑战性的课题。这是因为: ①HRI中的农田区域内部经常出现较大的光谱变化; ②在同一地区, 不同农田的形状、大小存在很大的差异; ③一些相邻农田的边界十分模糊。以上因素都给传统的影像分割带来了一定困难。

本文在Mean-Shift(MS)遥感影像分割算法的基础上进行了改进, 使之更好地适用于农田HRI的分割。之所以选用MS算法, 主要有2个原因: ①MS在GEOBIA中应用广泛^[6], 具有抗噪性好^[7]、能处理较大尺度影像^[8]等优点; ②基于MS的HRI分割算法可分为2步, 第一步是边界保持平滑滤波^[9]; 第二步是空间聚类。其中第一步可以在不使边界模糊的前提下有效地平滑农田内部的光谱值, 从而提高农田HRI的分割精度。众多学者对基于MS的遥感影像分割算法进行了较为深入的研究。Wang等^[10]将一种监督学习的边界提取算法与MS滤波相结合, 以优化较大尺度HRI的分割效果。Ponti^[11]利用MS算法处理从遥感影像中提取的植被指数, 实现了MS在农田遥感影像分割中的应用。Lang等^[7]提出了一种可变尺度参数的MS滤波算法, 并将其用于抑制极化SAR影像的斑点噪声。Banerjee等^[12]利用k最近邻(k-nearest neighbor, kNN)算法自动估计MS的尺度参数, 以提高MS分割算法的自动化程度。Huang等^[13]发展了一套自适应尺度参数选择的MS分割算法, 用以分割高光谱影像。Comaniciu等^[14]提出了一种数据驱动的可变尺度MS分割算法。Ming等^[6]则利用空间统计方法来估计可变尺度MS分割算法中的尺度参数。纵观上述研究可以发现, 大多MS分割算法仅采用了全局或单一的尺度参数; 一些可变尺度的MS算法在其尺度参数估计中也仅考虑了光谱变化信息, 而忽略了边界等空间结构信息。在实际农田的HRI中, 不同农田的面积和田间光谱变化具有较大差异, 这都会给常规MS分割策略带来困难。

为了提高MS分割算法在农田HRI分割中的性能, 本文提出了一种可变尺度的MS分割算法, 其过程主要分为3步: ①将对角化的尺度参数矩阵与采样点密度估计模型相结合, 导出可变尺度MS滤波的迭代函数; ②利用局部光谱变化信息与矢量梯度边界强度, 提出一种新的局部尺度参数估算方法; ③将MS滤波结果输入到基于分形网络演化方法(fractal net evolution approach, FNEA)^[15]的区域生长算法中, 得到农田HRI的分割结果。并利用中国和美国的2景HRI进行算法验证, 以证明该算法可显著优化农田HRI的分割精度。

1 算法原理

1.1 可变尺度MS滤波

MS最早由Fukunaga等^[16]提出, 用于多变量特征空间的密度中心估计。Comaniciu等^[9]将MS运用到影像特征空间分析中, 发展了基于MS的影像分割算法。该算法主要包含 MS滤波和基于空间聚类的影像分割2个步骤。其中第一步的目的是获取每个像素点的“ mode” ^[9], 即像素点所属类块的密度中心。这实际上是对原始影像进行了平滑滤波, 但又与普通的低通平滑滤波不同— — 由于采用了基于核函数的密度梯度估计, MS滤波可以有效地保持区域边界信息, 因此该方法也被称为“ 边界保持平滑滤波” , 该步骤对提高分割精度非常关键。

MS滤波的一般步骤是: ①对于影像中任意一个像素x初始化, 选代次数j=1并且y₁=x; ②当j< j_max且||y_j-y_j+₁||> t, 则第j+1次迭代滤波结果为

y_j+₁= $\begin{matrix} \frac{\sum_{x_{i} \in N (y_{j})} x_{i} K_{H} (x_{i} - y_{j})}{\sum_{x_{i} \in N (y_{j})} K_{H} (x_{i} - y_{j})} \end{matrix}$ , (1)

式中: t为收敛条件阈值(1个很小的正数); N()为滤波窗口内的像素集合; K_H()为核函数(多采用高斯函数形式); H为K_H()的尺度参数(亦称其为“ 带宽参数” ); ③得到滤波结果z_i=y_i_,_j+₁。

通常, 尺度参数有2种表达形式, 其中最简单的形式 ^[9] 为

H₁=h²I , (2)

h=h_s=h_r , (3)

式中h_s和h_r分别为空间和光谱尺度参数。这一形式的优点是只需给出1个尺度参数h, 即可完成滤波计算。在实际计算中, h_s与h_r可以不相等, 此时的核函数为

$\begin{matrix} K_{h_{s}, h_{r}} \end{matrix}$ (x)=C₁k $\begin{matrix} ({||\frac{x_{s}}{h_{s}}||}^{2}) \end{matrix}$ k $\begin{matrix} ({||\frac{x_{r}}{h_{r}}||}^{2}) \end{matrix}$ , (4)

式中: k()为核函数K()的轮廓函数^[9]; C₁为归一化常数。

对于多光谱遥感影像, 固定的h_r很难适用于大多数实际情况, 因为不同波段对同一地物的响应可能会有较大差别。所以第二种尺度参数的形式^[9]更接近实际情况, 即

H₂=diag( $\begin{matrix} h_{s}^{2} \end{matrix}$ , $\begin{matrix} h_{s}^{2} \end{matrix}$ , $\begin{matrix} h_{r, 1}^{2} \end{matrix}$ , $\begin{matrix} h_{r, 2}^{2} \end{matrix}$ , …, $\begin{matrix} h_{r, p}^{2} \end{matrix}$ ) , (5)

式中: p为波段数; diag()表示对角矩阵。为了简便, 将式(5)写成向量形式, 即

h=( $\begin{matrix} h_{s}^{2} \end{matrix}$ , $\begin{matrix} h_{s}^{2} \end{matrix}$ , $\begin{matrix} h_{r, 1}^{2} \end{matrix}$ , $\begin{matrix} h_{r, 2}^{2} \end{matrix}$ , …, $\begin{matrix} h_{r, p}^{2} \end{matrix}$ ) 。 (6)

然后在式(6)的基础上, 进而提出改变了的核函数形式, 即

$\begin{matrix} K_{H_{2}} \end{matrix}$ =C₂k $\begin{matrix} ({||\frac{x_{s}}{h_{s}}||}^{2}) \end{matrix} \begin{matrix} \overset{p}{\prod_{b = 1}} \end{matrix}$ k $\begin{matrix} ({||\frac{x_{r, b}}{h_{r, b}}||}^{2}) \end{matrix}$ , (7)

式中: C₂为归一化常数; h_r,_b为波段b的尺度参数。将式(7)带入式(1), 可得基于H₂的MS滤波迭代函数, 即

y_j+₁= $\begin{matrix} \frac{\sum_{x_{i} \in N (y_{j})} x_{i} k ({||\frac{y_{j, s} - x_{i, s}}{h_{s}}||}^{2}) \overset{p}{\prod_{b = 1}} k ({||\frac{y_{j, r} - x_{i, r}}{h_{r, b}}||}^{2})}{\sum_{x_{i} \in N (y_{j})} k ({||\frac{y_{j, s} - x_{i, s}}{h_{s}}||}^{2}) \overset{p}{\prod_{b = 1}} k ({||\frac{y_{j, r} - x_{i, r}}{h_{r, b}}||}^{2})} \end{matrix}$ 。 (8)

虽然式(8)考虑了不同波段尺度参数的差异, 但仍难以满足实际情况的要求。在HRI成像过程中, 由于气候条件和人为干扰等因素, 同一波段对不同地点的同种地物会存在光谱响应上的差异。因此, 需要在H₂的基础上考虑局部尺度参数。1.2节将详细叙述提出的局部尺度参数估计方法, 此处只给出局部尺度参数(即可变尺度的滤波迭代函数)。

可变尺度的MS滤波算法有2种模型, 分别是气球和采样点估计函数, 在文献[14]中均有详细介绍, 此处不再赘述。因为文献[16]指出第一种模式的效果较差, 本文采用第二种模式。在文献[14]中式(17)的基础上, 将本文修改的式(6)代入其中。由于本文是按照式(7)的假设来进行推导的, 即认为空间信息与光谱信息是相互独立的, 因此在密度中心估计中像素的位置与各波段光谱信息是连乘的形式。仿照式(8), 将其中随窗口中心变化的尺度参数替换为随任意像素位置变化的尺度参数形式, 可得到本文最终的MS迭代函数, 即

y_j+₁= $\begin{matrix} \frac{\sum_{x_{i} \in N (y_{j})} \frac{x_{i}}{h_{s}^{2} \overset{p}{\prod_{b = 1}} h_{x_{i}, r, b}} k ({||\frac{y_{j, s} - x_{i, s}}{h_{s}}||}^{2}) \overset{p}{\prod_{b = 1}} k ({||\frac{y_{j, r} - x_{i, r}}{h_{x_{i}, r, b}}||}^{2})}{\sum_{x_{i} \in N (y_{j})} \frac{1}{h_{s}^{2} \overset{p}{\prod_{b = 1}} h_{x_{i}, r, b}} k ({||\frac{y_{j, s} - x_{i, s}}{h_{s}}||}^{2}) \overset{p}{\prod_{b = 1}} k ({||\frac{y_{j, r} - x_{i, r}}{h_{x_{i}, r, b}}||}^{2})} \end{matrix}$ 。 (9)

1.2 局部尺度参数估计

在农田HRI中, 要将各个农田完整地分割出来是较为困难的。一方面, 由于灌溉条件和耕作习惯等因素, 即便在同一片农田中作物之间的生长速率也存在差异, 这在HRI中反映为农田内部的光谱变化, 会导致同一片农田被错误地分为若干区域; 另一方面, 农田之间经常被较窄的小道或灌渠分隔, 在HRI中以模糊边界的形式显示出来, 这可能会导致一些光谱特征相近的农田被分割为同一片区域。综上可见, 局部光谱变化和边界强度信息对确定MS滤波参数是非常重要的。

前人在设计局部尺度参数的估计模型中, 基本上只考虑了局部光谱变化信息^{[6, 14]}。从定性上来说, 较大的尺度参数易引起过平滑(即边界信息的丢失); 反之, 较小的尺度参数难以对光谱变化较大的农田区域进行较好的平滑。鉴于此, 在提出的尺度估计模型中, 对边界强度较高的区域, 其尺度参数的估计结果应适当偏小; 而对于变化较大的农田内部, 其尺度估计结果应适当偏大。综合考虑以上因素, 局部尺度参数应与边界强度成反比, 而与农田内部的光谱变化成正比。为了避免因局部光谱变化较大而导致其被错误地当做边界信息, 在本文模型中, 采取了单调递减指数函数形式将边界强度信息与局部光谱信息结合到一起, 使其更适合农田HRI的分割。其计算公式为

h_r_{, edge}(x, y, b)=exp[-E_n(x, y)]h_r(x, y, b) , (10)

式中: E_n(x, y)为像素(x, y)的归一化边界强度值; h_r(x, y, b)为像素(x, y)在波段b的尺度参数, 即

h_r(x, y, b)= $\begin{matrix} \frac{c}{10} \end{matrix}$ exp $\begin{matrix} [- \frac{\ln 2}{c} (\underset{\frac{}{v \in W (x, y, b)}}{argmax v} - \underset{\frac{}{v \in W (x, y, b)}}{argmin v})] \end{matrix}$ , (11)

c= $\begin{matrix} \underset{v \in I (b)}{argmax v} \end{matrix}$ - $\begin{matrix} \underset{v \in I (b)}{argmin v} \end{matrix}$ 。 (12)

其中: v为像素的光谱强度值; W(x, y, b)为以像素(x, y)为中心的窗口所包含的波段b子影像; I(b)为波段b的整景影像; c为影像最大灰度差异常数; 10与ln 2为经验常数, 均通过大量实验得出。由式(11)可以看出, 对于变化较大的局部区域, 其尺度参数较小, 从而防止过度平滑。

HRI通常具有多个波段, 因此采用了一种基于矢量梯度的边界强度信息提取方法^[17]。该方法需要首先计算出2个正交方向的梯度影像, 本文采用高斯差分滤波来实现这一步骤, 得到的梯度矢量为

V(x, y)= $\begin{matrix} {(\begin{matrix} \partial g_{h, 1} & \partial g_{h, 2} & \dots & \partial g_{h, b} \\ \partial g_{v, 1} & \partial g_{v, 2} & \dots & \partial g_{v, b} \end{matrix})}^{T} \end{matrix}$ (13)

式中􀆟g_{h, 1}和􀆟g_{v, 1}分别为波段1水平和竖直方向的梯度值。边界强度值即为2× 2矩阵V^T· V的最大特征值λ 。在实际计算中需要采用归一化的边界强度, 即

E_n(x, y)= $\begin{matrix} \frac{λ (x, y)}{\underset{(p, q) \in I}{argmax} λ (p, q)} \end{matrix}$ 。 (14)

从式(10)可以看出, 对于影像中灰度变化较大或边界强度较高的部位, 所估计的局部尺度参数较小, 这可以避免过大尺度可能导致的过度平滑。另外, 对于内部灰度变化较大的区域, 由于其边界强度通常较小(农田内部的灰度变化通常较为平缓), 因此其估计的尺度参数不会特别小, 这有利于区域内部的平滑。对于边界较为模糊的部位, 式(10)第一项依然可以使其局部尺度较低, 从而降低过度平滑的概率。

1.3 基于FNEA的遥感影像分割

本文的主要目的是利用MS滤波算法提高农田HRI的分割效果, 因此在空间滤波阶段采用了一种主流的遥感影像分割算法— — FNEA来实现整体尺度的HRI分割。FNEA是一种基于区域生长的空间聚类算法, 在初始阶段, 各个像素被当作单独的斑块; 然后通过迭代合并相似的斑块, 使均一的区域被分割出来。FNEA利用均一性变化量来衡量2个图斑适于合并的程度, 而均一性变化量的计算则考虑了形状和光谱信息^[14]。在算法执行之前, FNEA需要设置形状异质性权重f_s、紧凑异质性权重f_c与斑块大小阈值f_size这3个参数; 前2个参数分别影响分割过程中光谱和形状信息在均一性变化量计算中的比重, f_size则决定分割结果中图斑的平均大小。一般而言, f_s与f_c分别被默认为0.1和0.5。由于本文主要考察MS滤波结果对农田HRI的分割效果, 因此在实验中FNEA的f_s与f_c均采用了默认设置。图1示出了本文算法的流程。

	Figure Option View Download New Window
	图1 本文算法流程Fig.1 Flowchart of proposed method

FNEA的实现参考了文献[18], 利用基于红黑搜索树的全局最优合并来实现快速、准确的HRI分割。本文算法基于C++语言完成, 实验平台为Windows7操作系统, 4Gb RAM, CPU为Intel Core I5-4200M。

2 算法实验

2.1 实验数据

本文实验共采用了2景农田HRI, 分别为RapidEye和OrbView3影像(以下简称为S1和S2), 其参数见表1。

表1 实验数据信息 Tab.1 Information about experiment data

上述2种数据均包含蓝光、绿光、红光和近红外波段, 而RapidEye还提供红边波段(690~740 nm)。为便于对比分割效果, 对S1和S2影像均采用了绿(Band2)、红(Band3)、近红外(Band4)波段组合。S1和S2影像Band4(R)Band3(G)Band2(B)假彩色合成结果如图2所示。

	Figure Option View Download New Window
	图2 实验用Band4(R)Band3(G)Band2(B)假彩色合成影像Fig.2 False color images composed of Band4(R)Band3(G)Band2(B)adopted for experiment

局部影像的位置分别见图2中的白色和蓝色方框, 下文中分别将S1中的白色和蓝色方框内的局部影像简称为S1a和S1b, 将S2中的白色和蓝色方框内的局部影像简称为S2a和S2b。在图2中, S1影像显示了中国山东地区的农田, 其田块的平均面积较小, 很多由田间小道和灌渠形成的边界比较模糊, 并且大多相邻的田块具有相似的光谱特征。相反, S2影像显示了美国加利福尼亚地区的农田, 其田块平均面积较大, 形状大多为规则的矩形, 而且许多田块之间的光谱变化较大。以上因素都会为准确分割影像带来一定的困难, 因此, 对S1和S2影像的分割实验可较为全面地验证本文算法对农田HRI分割的优化效果。

2.2 对比实验与结果分析

2.2.1 实验设置

为了验证本文算法的性能, 采用了与其他3种算法进行对比的实验, 包括基于可变尺度的MS分割算法^[14](简称“ VMS” )以及2种较新的MS分割算法。VMS采用的迭代滤波计算式与本文中的式(7)相同, 其尺度参数的估计仅考虑了光谱信息。另2种较新的算法中, 一种算法在文献[12]中被提出, 首先采用kNN进行全局最优光谱尺度参数h_r的估计, 然后利用估计的h_r进行MS滤波; 另一种算法是边界加权(edge weighted, EW)的MS滤波^[10]。

VMS, kNN和EW算法与本文算法均需要给出空间尺度参数h_s, 该参数也被当作滤波窗口的尺度。经过多次实验, 将4种算法的h_s均设置为5, 一些相关研究^{[7, 8, 10]}也都采用了相同的设置。除h_s外, EW还需设置h_r。在本文实验中发现, EW采用kNN估计的h_r可以得到满意的效果。另外, EW的输入还包括边界强度影像, 由1.2节中的矢量梯度方法提供。

为便于对比, 对4种MS滤波算法均采用FNEA完成最终的HRI分割。而对于S1和S2影像, 其f_size需要分别通过测试和调节给出最佳设置。经多次实验, 发现S1和S2影像的最优f_size分别为2 000和35 000。S2影像的参数更大的原因是因为S2影像中农田的平均尺度远大于S1影像。

为了定量地给出4种算法对农田HRI分割的优化效果, 采用了一种定量的遥感影像分割评价方法^{[18, 19, 20]}。该方法首先给出分割结果的精确率(P)和召回率(R), 分别反映过分割和亚分割错误, 并且数值范围均为[0, 1]。当P和R均接近1时, 说明分割结果是精确的; P接近1而R接近0时, 表明存在过分割错误; 反之, P接近0而R接近1时, 说明存在亚分割错误。再利用该方法根据P和R计算出F, 以全面地反映分割精度。同样, F的范围也是[0, 1], 其值越接近1说明越精确。另外, 为了更清晰地展现4种MS滤波算法对分割效果的提升, 还与原始影像作为FNEA输入时的分割结果进行了对比。

由于S1和S2影像的尺度较大, 观察全景影像的分割结果不便于进行细致、可靠的分割效果评价。因此, 在S1和S2影像中分别选取了2个具有代表性的局部影像S1a, S1b和S2a, S2b进行更为详细的分割结果评价。

2.2.2 实验结果

图3和图4为本文算法的中间结果, 其中图3为矢量梯度方法提取的边界强度。

	Figure Option View Download New Window
	图3 矢量梯度边界强度Fig.3 Vector gradient edge strength

由图3可以看出, 大部分明显的边界都被提取出来了。虽然一些较弱的边界没有被明显地显示出来, 但其边界强度值依然高于区域内部。图4是各个波段的局部尺度影像。

	Figure Option View Download New Window
	图4 局部尺度参数估计结果Fig.4 Results of local scale parameter estimation

从图4可以看出, 边界位置基本具有较小的尺度参数; 而对于农田内部, 尤其是灰度变化较小的区域, 其尺度参数的数值较大。

对于滤波结果, kNN为S1和S2影像估计的最佳全局h_r分别为6.27和28.65, 前者小于后者主要是因为S1影像中田块的平均面积较小。

图5所示本文算法对S1和S2影像的最终分割结果。

	Figure Option View Download New Window
	图5 本文算法整体分割结果Fig.5 Whole segmentation results produced by proposed method

仔细观察图5可以看出, S1影像中主要的区域被较好地分割出来了; 而S2影像中较大的区域存在过分割现象。通过与S2影像原图对比可以发现, S2影像中大部分被过分割的农田, 其内部的灰度变化较大, 且其中存在树木和水塘等小尺度的地物, 这些都是导致过分割的主要原因。因传统MS, kNN与EW算法的最终分割结果与本文算法较为相似, 为节省篇幅, 未显示其整景分割结果。

表2给出了对4种算法分割结果的定量评价。

表2 4种算法分割实验结果定量评价 Tab.2 Quantitative evaluation of segmentation experiment results

从表2可见, 对于2景实验影像, 本文算法的F均最高, 这说明本文算法对农田HRI分割优化的效果较为显著。值得一提的是, 4种算法的P, R和F相差并不悬殊, 这在图6和图7局部影像的分割效果对比中也可以得到印证。

仔细观察图6第一行可发现, 本文算法的滤波结果对较大农田内部的平滑效果最佳(例如图6(a)右下部颜色较深的大块田地)。对于S1a和S1b影像, 本文算法都得到了较好的分割结果。例如, S1a影像右上部的几块农田之间有模糊的小道, 在图6第二行中, 只有本文算法较好地区分了这些田地(图6(e)); 其他算法在该区域普遍存在亚分割错误(即把相邻的农田合并到了一起)。对于S1b影像, 经目视对比可见, 4种算法的分割结果差异较小。kNN在S1b影像的分割中得到了较好的效果, 其中个别田块(例如图6(o)下部的矩形农田)被较为完整地分割出来。但在S1b影像右下部的细长农田分割中, 本文算法将其中1个农田错误地分为2个小区域(图6(m)); 除此以外, 本文算法在S1b影像其他区域的分割中都取得了较为不错的效果。

	Figure Option View Download New Window
	图6 S1局部影像的实验结果Fig.6 Test results of S1 subsets

与图6不同的是, 图7中的田块具有更大的面积, 但田块内部的光谱变化也更为显著。一些农田(特别是S2b影像中的田块)具有较大的光谱变化, 以至于4种算法都产生了过分割错误。在S2a影像的分割结果中, EW的过分割错误最为明显; 在S2b影像中, 4种算法均有明显的过分割错误。

	Figure Option View Download New Window
	图7 S2局部影像的实验结果Fig.7 Test results of S2 subsets

对图7进一步观察可以看出, 本文算法、VMS与kNN的分割结果相差较小。然而表2中的定量评价结果显示, 本文算法优于其他3种算法, 这主要是因为图7只显示了一部分面积较小农田的分割结果, 而表2中的定量评价结果是基于整景分割结果计算出来的。一些面积较大的田块被传统MS算法和kNN误分割为若干区域, 则是导致其F较低的主要原因。S2影像的分割实验可以更可靠地说明, 本文算法对农田HRI分割具有较好的优化效果。

3 结论

本文发展了一种基于可变尺度Mean-Shift(MS)的农田高分遥感影像(HRI)分割算法, 其主要贡献在于将局部灰度变化与边界强度信息加入到了局部尺度参数估计中, 从而较为显著地提高了MS滤波的边界保持能力与区域平滑效果。利用2景较大尺度的农田HRI进行了实验验证, 结果表明本文算法具有如下优势:

1)较高程度的自动化。由于本文提出的局部尺度参数估计方法不需要预先设置任何参数, 因而大大方便了用户的操作。

2)优越的边界保持平滑滤波。实验结果说明本文算法对农田HRI分割的优化效果要好于其他MS滤波算法。

本文算法的运算速度相对较慢, 对RapidEye影像进行分割用时近2 min。在未来可以考虑通过优化编程代码或简化收敛条件来提升计算效率。在下一步工作中, 也可以考虑本文算法在城市HRI或高分SAR影像分割中的应用。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Blaschke T, Hay G J, Kelly M, et al. Geographic object-based image analysis-towards a new paradigm[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 87: 180-191. [本文引用:2]
[2]	Kim H O, Yeom J M. Effect of red-edge and texture features for object-based paddy rice crop classification using Rapideye multi-spectral satellite image data[J]. International Journal of Remote Sensing, 2014, 35(19): 7046-7068. [本文引用:2]
[3]	Vieira A, Formaggio A R, Rennó C D, et al. Object based image analysis and data mining applied to a remotely sensed Land sat time-series to map sugarcane over large areas[J]. Remote Sensing of Environment, 2012, 123: 553-562. [本文引用:1]
[4]	Peña-Barragán J M, Ngugi M K, Plant R E, et al. Object-based crop identification using multiple vegetation indices, textural features and crop phenology[J]. Remote Sensing of Environment, 2011, 115(6): 1301-1316. [本文引用:1]
[5]	Liu M W, Ozdogan M, Zhu X J. Crop type classification by simultaneous use of satellite images of different resolutions[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(6): 3637-3649. [本文引用:1]
[6]	Ming D P, Li J, Wang J Y, et al. Scale parameter selection by spatial statistics for GEOBIA: Using mean-shift based multi-scale segmentation as an example[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 106: 28-41. [本文引用:3]
[7]	Lang F K, Yang J, Li D R, et al. Mean-shift-based speckle filtering of polarimetric SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(7): 4440-4454. [本文引用:3]
[8]	Michel J, Youssefi D, Grizonnet M. Stable mean-shift algorithm and its application to the segmentation of arbitrarily large remote sensing images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(2): 952-964. [本文引用:2]
[9]	Comaniciu D, Meer P. Mean shift: A robust approach toward feature space analysis[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 24(5): 603-619. [本文引用:6]
[10]	Wang L G, Liu G Y, Dai Q L. Optimization of segmentation algorithms through mean-shift filtering preprocessing[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(3): 622-626. [本文引用:3]
[11]	Ponti M P. Segmentation of low-cost remote sensing images combining vegetation indices and mean shift[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(1): 67-70. [本文引用:1]
[12]	Banerjee B, Varma S, Buddhiraju K M, et al. Unsupervised multi-spectral satellite image segmentation combining modified mean-shift and a new minimum spanning tree based clustering technique[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(3): 888-894. [本文引用:1]
[13]	Huang X, Zhang L P. An adaptive mean-shift analysis approach for object extraction and classification from urban hyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(12): 4173-4185. [本文引用:1]
[14]	Comaniciu D, Ramesh V, Meer P. The variable band width mean shift and data-driven scale selection[C]//Proceedings of the 8th International Conference of Computer Vision. Vancouver: IEEE, 2001: 438-445. [本文引用:4]
[15]	Baatz M, Schäpe A. Multiresolution segmentation: An optimization approach for high quality multi-scale image segmentation[M]//Strobl J, Baschke T, Griesebner G. Angewand te Geographische Informations Verarbeitung XII. Karlsruhe: Wichmann Verlag, 2000: 12-23. [本文引用:1]
[16]	Fukunaga K, Hostetler L. The estimation of the gradient of a density function, with applications in pattern recognition[J]. IEEE Transactions on Information Theory, 1975, 21(1): 32-40. [本文引用:1]
[17]	Qin A K, Clausi D A. Multivariate image segmentation using semantic region growing with adaptive edge penalty[J]. IEEE Transactions on Image Processing, 2010, 19(8): 2157-2170. [本文引用:1]
[18]	苏腾飞, 李洪玉. 一种两阶段区域生长的遥感图像分割算法[J]. 遥感技术与应用, 2015, 30(3): 476-485. Su T F, Li H Y. A two stage region growing method for remote sensing image segmentation[J]. Remote Sensing Technology and Application, 2015, 30(3): 476-485. [本文引用:1]
[19]	Jiang H B, Su Y Y, Jiao Q S, et al. Typical geologic disaster surveying in Wenchuan 8. 0 earthquake zone using high resolution ground LiDAR and UAV remote sensing[C]//Proceedings of SPIE 9262, Lidar Remote Sensing for Environmental Monitoring XIV. Beijing, China: SPIE, 2014: 926219. [本文引用:1]
[20]	Yi L N, Zhang G F, Wu Z C. A scale-synthesis method for high spatial resolution remote sensing image segmentation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(10): 4062-4070. [本文引用:1]

2014

0.0

... 面向地理对象的影像分析(geographical object based image analysis,GEOBIA)^[1]在农田遥感分类研究中逐渐增多^[2,3,4],这主要归因于2个方面: ①农田遥感监测需要高空间分辨率遥感影像(high spatial resolution remote sensing image,HRI),以准确提取作物信息 ...

... 很多地区(如中国)的农田面积较小^[5],导致中低空间分辨率的遥感影像易产生混合像素光谱,HRI可有效地解决这一问题,很多研究都表明GEOBIA适用于HRI的解译^[1,2] ...

2014

0.0

2012

0.0

2011

0.0

2014

0.0

2015

0.0

... 之所以选用MS算法,主要有2个原因: ①MS在GEOBIA中应用广泛^[6],具有抗噪性好^[7]、能处理较大尺度影像^[8]等优点 ...

... Ming等^[6]则利用空间统计方法来估计可变尺度MS分割算法中的尺度参数 ...

... 前人在设计局部尺度参数的估计模型中,基本上只考虑了局部光谱变化信息^[6,14] ...

2014

0.0

... 之所以选用MS算法,主要有2个原因: ①MS在GEOBIA中应用广泛^[6],具有抗噪性好^[7]、能处理较大尺度影像^[8]等优点 ...

... Lang等^[7]提出了一种可变尺度参数的MS滤波算法,并将其用于抑制极化SAR影像的斑点噪声 ...

... 经过多次实验,将4种算法的h_s均设置为5,一些相关研究^[7,8,10]也都采用了相同的设置 ...

2015

0.0

... 之所以选用MS算法,主要有2个原因: ①MS在GEOBIA中应用广泛^[6],具有抗噪性好^[7]、能处理较大尺度影像^[8]等优点 ...

... 经过多次实验,将4种算法的h_s均设置为5,一些相关研究^[7,8,10]也都采用了相同的设置 ...

2002

0.0

... ②基于MS的HRI分割算法可分为2步,第一步是边界保持平滑滤波^[9] ...

... Comaniciu等^[9]将MS运用到影像特征空间分析中,发展了基于MS的影像分割算法 ...

... ^[9],即像素点所属类块的密度中心 ...

... 通常,尺度参数有2种表达形式,其中最简单的形式 ^[9] 为 ...

... 式中: k()为核函数K()的轮廓函数^[9] ...

... 所以第二种尺度参数的形式^[9]更接近实际情况,即 ...

2014

0.0

... Wang等^[10]将一种监督学习的边界提取算法与MS滤波相结合,以优化较大尺度HRI的分割效果 ...

... 另一种算法是边界加权(edge weighted, EW)的MS滤波^[10] ...

... 经过多次实验,将4种算法的h_s均设置为5,一些相关研究^[7,8,10]也都采用了相同的设置 ...

2013

0.0

... Ponti^[11]利用MS算法处理从遥感影像中提取的植被指数,实现了MS在农田遥感影像分割中的应用 ...

2014

0.0

... Banerjee等^[12]利用k最近邻(k-nearest neighbor,kNN)算法自动估计MS的尺度参数,以提高MS分割算法的自动化程度 ...

2008

0.0

... Huang等^[13]发展了一套自适应尺度参数选择的MS分割算法,用以分割高光谱影像 ...

2001

0.0

... Comaniciu等^[14]提出了一种数据驱动的可变尺度MS分割算法 ...

... 前人在设计局部尺度参数的估计模型中,基本上只考虑了局部光谱变化信息^[6,14] ...

... FNEA利用均一性变化量来衡量2个图斑适于合并的程度,而均一性变化量的计算则考虑了形状和光谱信息^[14] ...

... 为了验证本文算法的性能,采用了与其他3种算法进行对比的实验,包括基于可变尺度的MS分割算法^[14](简称#cod#x0201c ...

2000

0.0

... ③将MS滤波结果输入到基于分形网络演化方法(fractal net evolution approach,FNEA)^[15]的区域生长算法中,得到农田HRI的分割结果 ...

1975

0.0

... 1 可变尺度MS滤波MS最早由Fukunaga等^[16]提出,用于多变量特征空间的密度中心估计 ...

2010

0.0

... HRI通常具有多个波段,因此采用了一种基于矢量梯度的边界强度信息提取方法^[17] ...

2015

0.0

苏腾飞, 李洪玉. 一种两阶段区域生长的遥感图像分割算法[J]. 遥感技术与应用, 2015, 30(3): 476-485.

Su T

, Li H

A two stage region growing method for remote sensing image segmentation[J]. Remote Sensing Technology and Application, 2015, 30(3): 476-485.

How to improve the accuracy and speed of image segmentation algorithm that has a significant effect upon the interpretation of remote sensing images.In this paper,a region\|growing\|based method is proposed for remote sensing image segmentation.This method is composed of two procedures which are local the best merge and global best merge.The first step focuses on performing high\|speed image segmentation.In implementation,a threshold is introduced into the first step to reduce erroneous region merges.The second step aims at increasing segmentation accuracy.In order to raise the running speed,an advanced data structure and red\|black tree are utilized to implement the second step.Finally,simulated remote sensing image and Orbview\|3 high resolution images are used to carry out segmentation experiment.A supervised image segmentation accuracy evaluation method is utilized to quantitively validate the performance of our algorithm.The experiment result indicates that the proposed method can achieve satisfactory segmentation in terms of segmentation accuracy and speed.

提高图像分割算法的精度与速度对遥感影像的解译工作具有重要意义。提出了一种基于区域生长的遥感图像分割算法,包括的两个步骤分别是局部最优合并和全局最优合并。第一步注重提高图像分割的速度,在其具体实现中引入了一个局部最优合并的阈值,以减少错误的合并;第二步侧重提高图像分割的精度,在其实现过程中利用了一种高级数据结构红黑树来提高搜索最优合并区域的速度。最后,利用模拟的遥感影像和Orbview\|3高分辨率影像开展了图像分割实验。利用一种监督的图像分割精度评价方法,定量评价了该方法的性能。实验结果表明:该方法在分割精度和分割速度方面都取得了令人满意的效果。

... 为了定量地给出4种算法对农田HRI分割的优化效果,采用了一种定量的遥感影像分割评价方法^[18,19,20] ...

2014

0.0

... 为了定量地给出4种算法对农田HRI分割的优化效果,采用了一种定量的遥感影像分割评价方法^[18,19,20] ...

2012

0.0

... 为了定量地给出4种算法对农田HRI分割的优化效果,采用了一种定量的遥感影像分割评价方法^[18,19,20] ...