类矩形引导的玉米遥感图像分割方法基本流程如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 方法流程示意图Fig.1 Framework of proposed method

原始GF-1数据经过Gram-Schmidt(GS)算法^[6](一种光谱锐化高保真的正交化图像融合算法)融合处理得到融合图像; 利用SUSAN算子对融合图像进行边缘提取, SUSAN模板的灰度差阈值通过迭代的方式自适应地确定; 利用类矩形引导的分割方法进行图像分割, 通过反馈确定分割规则; 最后对分割结果进行实验对比分析。

本文处理的图像是使用GS算法对原始GF-1的多光谱和全色数据融合得到的, 然后利用SUSAN算子对融合图像进行边缘提取。本文在SUSAN算子的基础上加入一种自适应选取阈值的方法, 以减少原算法中人为选取灰度差阈值对处理效果的影响。SUSAN算子能够在拐角处得到较为连续的边界, 而且不涉及梯度运算, 因此对光谱复杂图像的边缘检测效果优于传统的边缘检测算法; 另外, 处理结果与所用模板的尺寸和颜色、对比度参数选取的相关性不大, 能较好地保持图像的特征结构细节信息, 获得比较光滑的边缘图像, 这是在边缘提取时选用SUSAN方法的主要原因^[7]。具体步骤如下:

1)设置一个5像元× 5像元的方形检测模板, 以融合图像中的每一个像元为模板中心在图像中移动; 对于每个模板, 通过计算模板内每个像元与中心像元的灰度差得到该模板的灰度差直方图; 最后根据灰度差直方图通过迭代法确定该模板的阈值, 使不同对比度的图像都能够自适应地计算出每个模板内合适的t值。由于每个模板的阈值是根据模板内像元的灰度差值确定的, 因此能很好地检测到不同对比度下的灰度变化, 使对吸收核同值区(univalue segment assimilating nucleus, USAN)的判断更加准确。

首先计算模板中每个像元I(x, y)与中心像元I(x_c, y_c)的灰度差值, 取其均值作为迭代初始值t₀, 即

t₀= 125[I(x, y)-I(x_c, y_c)] , (1)

然后根据迭代初始值将灰度差直方图分成2部分, 计算下一个迭代值t_i+1, 即

t_i+1= 12∑m=0tim·h(m)∑m=0tih(m)+∑m=ti+1Cmaxm·h(m)∑m=ti+1Cmaxh(m), (2)

式中: m为模板中像元点和与中心像元点的灰度差值; h(m)为模板中具有该灰度差值的点的数量; C_max为灰度差值的最大值。每进行一次迭代后进行判断, 若|t_i+1-t_i|=0, 则迭代停止, 取t_i+1作为USAN模板的最后灰度差阈值。因为每个模板的阈值t是根据模板内的灰度差值确定的, 因此能很好地检测到不同灰度对比度下的灰度变化, 使对USAN区的判断更加准确。

2)比较模板内部每个像元的灰度值和模板中心像元的灰度值, 若模板内某个像元的灰度与模板中心像元(核)灰度的差值小于阈值t, 则认为该像元与核具有相同(或相近)的灰度。满足这一条件的像元组成的区域称为USAN区。根据亮度比较函数计算模板中的USAN区域, 即

c(x, y)= 1|I(x, y)-I(xc, yc)|≤t0其他, (3)

式中: I(x, y)和I(x_c, y_c)分别为模板中心像元与模板中其他像元的灰度值; c(x, y)为属于USAN区的像元; x_c和y_c分别为模板中心像元的横、纵坐标; x和y分别为模板中除中心像元外的其他像元的横、纵坐标; t为确定相似程度的阈值。

3)以最大USAN区域的像元个数的3/4为全局USAN固定阈值, 查找小于该固定阈值的像元, 得到融合图像在SUSAN算子中的边缘响应值, 形成边缘图像。

计算出以每个像元为核心的USAN区域的像元个数, 即

n(x₀, y₀)= ∑(x, y)≠(x0, y0)c(x, y) , (4)

式中n(x₀, y₀)为每个模板中属于USAN区的面积。

4)在求得每个像元的USAN区后, 通过n(x₀, y₀)与某个固定阈值g比较得到边缘响应R(x₀, y₀), 即

R(x₀, y₀)= g-n(x0, y0)n(x0, y0)< g0n(x0, y0)≥g, (5)

式中: R(x₀, y₀)为边缘的响应, 最终目标边缘图像中USAN区的面积越小, 边缘的响应将越大; g为固定阈值, 本文设置g=3n_max/4, n_max为n(x₀, y₀)中的极大值。

利用改进的SUSAN算子对图像进行边缘提取。自适应地选取灰度差阈值的方法可以提高算法的鲁棒性, 使算法更好地保持图像的结构细节特征信息, 得到较为光滑的边缘图像, 同时具有一定的噪声抑制能力。

本文在Felzenszwalb等^[5]提出的基于图的图像分割理论基础上, 通过测量区域间的差异来计算边界的最小权值, 对区域之间有较低权值的2个区域进行合并。该方法基于区域间间距和区域内间距2个特征来判断2个区域是否应该合并: 若2个区域的区域间间距明显大于其中任意1个区域的区域内间距, 则认为这2个区域之间存在明显的界限(即不可以合并)^[5]。

定义区域间的间距为

Dif(C₁, C₂)= minυi∈C1, υj∈C2, (υi, υj)∈Eω (υ _i, υ _j) , (6)

式中: C₁和C₂分别为2个区域; υ _i, υ _j分别为不同分割区域映射在无向图中的节点; ω 为权重。在分别属于2个区域C₁和C₂且有边连接的所有点对中, 寻找权重最小的那对(若2个区域内的点没有边相连, 则定义间距为正无穷大), 即2个区域之间连接边的最小长度。

定义区域内间距为

Int(C)= maxe∈MST(C, E)ω (e) , (7)

式中MST(C, E)为区域C的最小生成树, 这里定义为区域内部的最大连接边长。当其中1个区域很小时, Int(C)并不能很好地反映其区域内间距(极端的情况是当C只含一个节点时, Int(C)=0)。可引入一个考虑区域面积的函数来补偿, 即

τ _补(C)=k/|C| , (8)

式中: k为常数, 当k越大时, 可以区分2个区域的界限就越明显; |C|为区域C的面积^[7]。

根据目标区域特点(农业玉米种植区土地多数为矩形的地块), 为了将图像分割成特定形状的区域, 将τ 定义为类矩形引导的相关函数, 即

τ (C)=k CMER|C|, (9)

式中: C_MER为区域C的最小外接类矩形(minimum enclosing rectangle, MER)面积; CMER|C|为区域C的类矩形度的倒数, 即C的外接类矩形面积与区域C的面积的比值。通过设置k值, τ 对类矩形形状的倾向性可强可弱。本文定义的类矩形, 其任意区域内的角度范围是[85° , 95° ]。

通过比较区域内间距和区域间间距之间的关系, 可以得到区域合并的判据, 即

D(C₁, C₂)= 真Dif(C1, C2)> MInt(C1, C2)假其他, (10)

其中 MInt(C₁, C₂)=min[Int(C₁)+τ (C₁),

Int(C₂)+τ (C₂)] 。 (11)

图像分割算法既不能单一使用像元强度值的变化作为分割依据, 也不能使用单一的分割阈值来决定分割与否。这里以图像中的每个区域表示图上的一个节点, 每一条连接节点的无向边都具有一个权重, 以衡量其连接的2个节点之间的不相似度。根据相邻区域在特征值上变化速度的大小, 动态调整分割阈值, 对经过边缘提取处理的图像进行分割。假设输入是1个有n个节点和m条边的图G, 那么输出是一系列区域。具体步骤如下:

1)将边按照权重值以非递减方式排序。

2)最初的分割记为S(0), 即每个节点属于一个区域。

3)按照以下方式由S(q-1)构造S(q)。记第q条边连接的2个节点分别为V_i和V_j, 若在S(q-1)中V_i和V_j分别属于2个区域并且第q条边的权重小于2个区域的区域内间距, 则合并这2个区域; 否则令S(q) = S(q-1)。

4)从q=1到q=m, 重复步骤2)。

5)返回S(m), 即为所求分割区域集合。

本文所使用的GF-1图像数据的获取时间为2013年7月30日, 覆盖区域为陕西省蓝田县。该图像包括4个多光谱波段: 蓝(0.45~0.52 μ m)、绿(0.52~0.59 μ m)、红(0.63~0.69 μ m)和近红外(0.77~0.89 μ m)波段, 分辨率均为8 m; 1个全色波段(0.45~0.90 μ m), 分辨率为2 m。该地区属于玉米种植大区, 地势较为平整, 对研究玉米与水稻、建筑物等背景区域的分割比较适宜。所使用的GF-1数据已经过正射纠正、辐射定标和大气校正等预处理。图2(c)为利用GS算法对GF-1全色图像(图2(a))和多光谱图像(图2(b))融合后的图像。考虑到实验区玉米田的形状具有较强的区域性, 有针对性地设计了类矩形的权重函数来引导分割。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 GF-1融合图像及边界提取结果Fig.2 GF-1 fusion image and edge detection result

图2(d)是图2(c)经过灰度差阈值自适应选取的SUSAN算子提取的边界结果图, 该结果减少了人工设置参数所产生的干扰。从图2(d)中可以看出, 自适应选取灰度差阈值的SUSAN算子能较好地检测出不同对比度下的灰度变化, 没有出现边界大量丢失的现象; 检测出的边缘线细腻而光滑, 有较好的连续性, 有效地体现了玉米田地块的边缘特征, 且定位精度较高。

在基于图的分割算法中, 利用提取的边界区域构建无向图, 采用本文的类矩形引导分割规则进行分割处理。根据算法中的规定, 式(8)中的k是常数, k越大, 界定的可以区分2个区域的界限就越明显(表1)。

表1

Tab.1

表1(Tab.1)

表1 图像分割结果区域数目与k关系 Tab.1 Relationship between number of segmentation blocks and k k 30 50 100 200 250 300 区域数目 17 891 10 763 5 860 3 959 3 478 3 217

表1 图像分割结果区域数目与k关系 Tab.1 Relationship between number of segmentation blocks and k

从表1可以看出, 随着k的增大, 区域数目减少; 在k增大到 200以后, 区域数目的减少变缓。考虑到经过融合后的GF-1图像的像元尺寸为2.4 m× 2.4 m, 为了能更好地显示图像中的小目标, 同时避免因分割区域数量过大而导致的过分割现象, 通常在[200, 300]之间选取k值。本文选取k=250进行类矩形引导的分割。

图3―图6分别为本文分割方法结果图、基于图的分割结果图、分水岭分割结果图和人工解译玉米目标图。其中(a)都代表原图像整体分割结果, (c)和(b)分别为对应(a)中黑边框区域1和区域2的局部放大。

比较图3与图4分割结果图中局部放大区域可以看出, 加入类矩形阈值函数的本文分割方法在有效区分不同地物的同时, 减少了地块边缘的小块区域, 凸显了玉米地块的类矩形的形状特点。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 加入类矩形阈值函数的分割结果Fig.3 Segmentation results with near-rectangle guided threshold function

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 无类矩形阈值函数的分割结果Fig.4 Segmentation results without near-rectangle guided threshold function

图5中分水岭分割方法在玉米种植区分割结果显得更零碎些, 图像分割区域数目达到24 306个, 产生严重的过分割现象。这主要是由玉米田严重的同物异谱特征导致的, 特别是在地块边缘处产生了大量的小块区域。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 分水岭分割结果Fig.5 Segmentation results with watershed method

与图4比较, 本文的分割方法(图3)在有效分割不同光谱特征对象的同时, 减少了同物异谱带来的过分割现象, 保持了大面积玉米区的地块特征, 为进一步的玉米种植面积统计提供了有效的数据支持。

为了检验类矩形引导的分割结果与实际玉米作物分布的符合情况, 将本文方法的分割结果(图3)与人工解译的目标区域样本(图6)进行比较。人工解译样本是由专业人员结合影像目视解译与光谱特征分析获取的。通过解译获取本文的训练样本过程中, 首先选取局部区域进行地理坐标绘制, 投射到遥感图像上进行样本绘图; 再对整个图像区域进行人工解译; 最后对其进行野外调查和详细解译, 得到相对准确的目标区域参考分割样本。经比较可以看出, 本文提出的类矩形引导的分割结果比较符合玉米田的实际分布特征。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 人工解译玉米目标Fig.6 Artificial interpretation samples of corn fields

对不同分割结果的统计面积与人工解译获得的参考标准面积的比较(表2)也可看出, 本文方法的分割结果的统计面积更接近于玉米田的实际种植情况, 优于其他2种分割方法。

表2

Tab.2

表2(Tab.2)

表2 不同方法分割结果面积统计对比 Tab.2 Comparison of segmentation results with different methods(m2) 不同方法 分水岭分割 基于图的分割 本文方法 人工解译 参考标准 面积 277 469.82 357 462.16 429 578.94 455 363.71

表2 不同方法分割结果面积统计对比 Tab.2 Comparison of segmentation results with different methods(m2)