Segmentation of PolSAR Data Based on Mean-Shift and Spectral Graph Partitioning and Its Evaluation
-
摘要: 提出了一种基于均值漂移和谱图分割的极化SAR(PolSAR)影像分割方法。首先,通过均值漂移算法对PolSAR影像进行过分割处理,并基于Wishart统计分布和假设检验的方法构建边缘检测器,充分利用了PolSAR影像的全极化信息提取边缘信息;然后,在过分割和边缘信息的基础上构建相似性度量矩阵,并采用归一化割准则实现PolSAR影像的分割。该算法充分利用了均值漂移算法过分割的特点,降低了谱图分割算法的运算代价,并结合谱图分割算法全局优化的优点改善了PolSAR影像的分割结果;最后,利用Radar-sat-2全极化影像进行了实验,并采用改进的分割效果评价方法实现了精度评价。实验表明,该算法有效地实现了PolSAR影像的分割,显著提高了谱图分割算法的效率,分割结果优良,分割精度优于eCognition软件中的多尺度分割方法。
-
关键词:
- 极化合成孔径雷达 /
- 影像分割 /
- 均值漂移 /
- 谱图分割 /
- 归一化割中图法分类号:P237.3 /
- TP79 文献标志码:A
Abstract: A new segmentation approach of polarimetric synthetic aperture radar(PolSAR)data is pro-posed based on mean shift and spectral graph partitioning.First,Mean-shift algorithm is used to gen-erate the over-segmentation of PolSAR image.In order to extract edge information,we apply a set ofdetectors based on the Wishart distribution with a hypothesis testing method that has fully consideredthe polarization information in PolSAR images.Then,a similarity matrix is constructed based on theover-segmentation results and image edge information.The graph partitioning process is performed u-sing the normalized cut criterion.With this method,we improve the segmentation efficiency of spec-tral graph partitioning based on the over-segmentation results generated by Mean-shift.The quality ofsegmentation results is also improved as a result of the global optimization of spectral graph partitio-ning algorithm.We applied this method on Radarsat-2full polarization images and evaluated the seg-mentation results.The experiement showed that this scheme can realize PolSAR segmentation effec-tively,speed up the original algorithm,and also demonstrates a better result than eCognition’s Multi-resolution segmentation approach.-
Keywords:
- PolSAR /
- segmentation /
- Mean-Shift /
- spectral graph partition /
- normalized cut
-
土壤水分是影响生态、水文和气候等环境过程的重要参数,土壤水分的准确获取和长期监测对于作物生长、旱情监测与预报等具有积极作用。重量法和时域反射法等传统的土壤水分测量方法需要定时定点测量,无法进行大面积长时效的实时观测。遥感技术和理论的不断发展为大范围土壤水分的获取提供了有效途径[1]。
微波遥感具有全天时、全天候工作的特点,且具有一定的穿透能力,已成为土壤水分反演的重要手段。被动微波遥感起步较早、算法成熟,但空间分辨率较低,在较小尺度研究中具有一定的局限性。主动微波遥感中的极化合成孔径雷达（polarimetric synthetic aperture radar,PolSAR）空间分辨率较高,同时可提供对土壤水分反演有利的极化特征参数,国内外学者已对此进行了大量研究。文献[2]提出了一种双分量模型（表面散射和体散射）的极化分解方法,忽略了二面角散射分量,在农田土壤水分反演中取得了较好的效果。文献[3]基于Radarsat-2数据,利用多种分解方法提取多种极化特征参数,对美国黑云杉林的土壤水分进行反演,效果良好。上述研究表明,PolSAR极化特征参数在土壤水分反演中潜力巨大。然而,土壤水分和极化特征参数之间并非简单的线性关系,通常为了提高反演精度,往往考虑较多参数。简单的线性回归分析方法无法准确、全面地衡量大量特征参数在土壤水分反演中的作用,而机器学习方法对于参数的数量和类型没有限制,能够较好地弥补简单回归分析方法在土壤水分反演中的不足。因此,当考虑多极化特征参数时,常常与机器学习方法相结合。
随机森林（random forest,RF）[4]和支持向量回归机（support vector regression,SVR）[5]是目前较为常用的机器学习方法。文献[6]将后向散射系数作为输入变量、土壤水分作为输出变量,对RF和SVR等12种机器学习方法的性能进行了综合评估,认为RF性能最佳,适用于土壤水分反演。文献[7]基于Radarsat-2数据提取多种极化特征参数并进行特征筛选,评价了RF在土壤水分反演中的适用性。结果表明,RF能较好地模拟冬小麦各生长阶段的土壤水分含量,其性能优于人工神经网络（artifical neural network,ANN）和SVR。综上,机器学习与PolSAR数据相结合能够较好地实现土壤水分模拟。然而,在植被覆盖地区,土壤水分的反演精度常常受到植被覆盖和地表粗糙度的影响,仅考虑极化特征具有一定局限性。
因此,不少学者协同利用PolSAR数据和光学数据,以提升土壤水分反演精度。文献[8]基于Sentinel-1/2数据,结合水云模型对冬小麦区域的土壤水分进行反演,结果表明,融合植被指数（fusion vegetation index,FVI）的反演效果最优,且能在一定程度上去除冬小麦对后向散射系数的影响。文献[9]基于Sentinel-1/2数据提取后向散射系数和植被指数,利用SVR和广义回归神经网络（generalized regression neural network,GRNN）反演了生长初期冬小麦区的土壤水分,认为加入植被指数的SVR模型效果最优。文献[10]基于Sentinel-1和Landsat 8影像,建立SVR模型,将后向散射系数、归一化植被指数（normalized difference vegetation index,NDVI）、地形因子等作为特征变量,在洪泛平原的水分反演中取得了较好效果,并认为NDVI能够在一定程度上减少植被覆盖的影响。上述研究表明,协同利用PolSAR和光学数据能够在一定程度上减少反演的不确定性,使反演结果更接近真实情况。然而,已有研究多基于C波段SAR数据,对于L波段在土壤水分反演中的效果有待深入[11]。
本文以内蒙古额济纳绿洲为研究区,利用L波段ALOS-2 PALSAR-2影像和Landsat 8影像提取雷达、光学特征参数,通过平均精度下降（mean decrease accuracy,MDA）评估各参数重要性,建立并对比雷达、光学、雷达-光学协同的多参数组合方案的随机森林土壤水分反演模型,对额济纳绿洲区进行土壤水分反演。
1 研究区与数据
1.1 研究区概况
额济纳绿洲位于内蒙古自治区西部,气候干旱,降水稀少,年均气温8.3 ℃,年均降水量约34 mm,潜在蒸发量约3 700~4 000 mm,属温带大陆性气候。2019⁃09⁃24-2019⁃09⁃29,在研究区共布设土壤采样点111个,采集表层0~10 cm的土壤样品,在实验室通过烘干称重法[12]获取各采样点的土壤水分含量。土壤采样点覆盖林地、草地、耕地和裸地等地物类型,实测土壤水分含量为0.02%~15.34%。研究区位置和采样点分布如图1所示。图1采用2019⁃09⁃24的Landsat 8影像,波段组合为7、5、4。
![]()
图 1 研究区位置及采样点分布图Figure 1. Location of the Study Area and Distribution of the Sampling Points1.2 遥感数据及预处理
2014年日本宇航局发射ALOS-2卫星,搭载L波段全极化合成孔径雷达传感器PALSAR-2,卫星轨道高度为628 km,重访周期为14 d。本文采用2019-09-29的ALOS-2 PALSAR-2升轨影像,包含水平极化（horizontal polarization,HH）、垂直极化（vertical polarization,VV）、水平-垂直极化（horizontal-vertical polarization,HV）、垂直-水平极化（vertical-horizontal polarization,VH）这4种极化方式,视向为右视,产品等级为Level 1.1级,中心入射角为33.87°,空间分辨率为6 m。利用SARScape软件首先对ALOS-2影像进行多视、滤波和地理编码,提取后向散射系数σVV、σVH、σHV、σHH；然后利用PolSARpro软件对ALOS-2影像进行滤波处理,通过极化分解提取极化特征参数；最后通过距离多普勒地形校正和重采样,得到与光学影像相匹配的后向散射系数和极化特征参数图像。
Landsat 8卫星于2013年由美国成功发射,卫星上携带运行陆地成像仪（operational land imager,OLI）和热红外传感器（thermal infrared sensor,TIRS）2种载荷。其中,OLI包括9个波段,除全色波段的空间分辨率为15 m外,其余波段的空间分辨率均为30 m；TIRS包括2个波段,空间分辨率为100 m,卫星重访周期为16 d。本文采用2019-09-24的Landsat 8影像,云量为0.02%。为了消除光照和大气等因素的影响,首先对Landsat 8影像进行辐射定标和大气校正,然后提取地表温度（land surface temperature,LST）和绿叶指数（green leaf index,GLI）等光学特征参数。
2 研究方法
2.1 雷达特征参数提取
极化分解是PolSAR影像目标特征描述与参数提取的主要方法[13],其目的在于将复杂的散射机制表示为几种具有对应散射矩阵的单一散射分量。根据分解机制和针对目标的不同,已有学者提出了Freeman-Durden分解[14]（FD）、Cloude-Pottier分解[15]（CP）、Yamaguchi分解[16]（YM）、An &Yang分解[17]（AY）、H/A/α分解[18]、van Zyl分解[19]（VZ）和Bhattacharya & Frery分解[20]（BF）等多种极化分解方法。基于ALOS-2 PALSAR影像和上述极化分解方法,提取36个雷达特征参数,具体见表1,其中σVV、σHH为同极化后向散射系数,σHV、σVH为交叉极化后向散射系数,下标Vol、Odd、Surf、Dbl和Hlx分别代表体散射、奇次、表面、二面角和螺旋体散射分量。
表 1 雷达特征参数Table 1. Radar Characteristic Parameters参数类型 处理方法 参数个数 雷达特征参数 后向散射系数 标准强度和相位处理 4 σVV、σVH、σHV、σHH 极化特征参数 Freeman-Durden分解[14] 3 FDVol、FDOdd、FDDbl Cloude-Pottier分解[15] 3 CPDbl、CPSurf、CPVol Yamaguchi分解[16] 3 YM3Vol、YM3Odd、YM3Dbl 4 YM4Vol、YM4Odd、YM4Dbl、YM4Hlx An &Yang分解[17] 3 AY3Vol、AY3Odd、AY3Dbl 4 AY4Vol、AY4Odd、AY4Dbl、AY4Hlx H/A/α分解[18] 5 Alpha、Anisotropy、Entropy、Lambda、RVI van Zyl分解[19] 3 VZVol、VZOdd、VZDbl Bhattacharya & Frery分解[20] 4 BFVol、BFOdd、BFDbl、BFHlx 2.2 光学特征参数提取
基于Landsat 8 OLI/TIRS数据,提取11个光学特征参数,其中植被指数通过波段反射率计算得到,LST通过ENVI扩展工具采用大气校正法提取。各参数计算公式参见文献[21-31],包含NDVI、综合光谱响应指数（combined spectral response index,COSRI）、增强型植被指数（enhanced vegetation index,EVI）、GLI、全球植被水分指数（global vegetation moisture index,GVMI）、红外百分比植被指数（infrared percentage vegetation index,IPVI）、植被供水指数（vegetation supply water index,VSWI）、土壤改良植被指数
（modified soil adjusted vegetation index,MSAVI）、温度植被干旱指数（temperature vegetation
drought index,TVDI）、比值植被指数（ratio vegetation index,RVI）和LST。
2.3 随机森林建模与评价
随机森林是一种基于决策树和Bagging的集成学习方法。决策树是依据一定的划分规则对变量空间进行拆分,并根据相似实例对数据集分组的树形结构。Bagging是一种从原样本集中随机抽取训练样本的技术,这意味着某些样本可能在训练中多次使用,而一些样本可能永远都不会被使用[32],未被使用的样本称为袋外（out of bag,OOB）数据,可用来评估RF的预测效果。MDA是RF中评估特征变量重要性的参数,通过随机打乱OOB数据中某一特征变量值,计算OOB误差,根据其下降程度评估特征变量重要性。通常,OOB数据误差越大,参数越重要。
本文依据重要性评分结果选取最优特征参数方案建立RF模型,包括2个关键参数,即决策树数量N和节点分裂时输入的特征变量数M。采用N=｛400,800,1 200,1 600,2 000,2 400｝,M=｛1,2,3,4,5,6｝的网格搜索,获得最优参数为N=1 600,M=3。模型精度评定采用K折交叉验证,K=10,精度评定指标采用决定系数R2和均方根误差（root mean square error,RMSE）。
3 结果与分析
3.1 参数重要性评价
由于特征参数较多,因此本文采用MDA分别对雷达特征参数和光学特征参数的重要性进行评分,结果如图2所示。
从图2（a）可以看出,从散射方式看,奇次散射、表面散射和体散射在雷达特征参数中的重要性较高,而二面角散射和螺旋体散射的重要性相对较低。原因可能是由于研究区内部虽有胡杨、柽柳等植被覆盖,但大部分地物还是以低矮草地、戈壁、沙漠为主,因此主要散射机制仍为奇次散射、表面散射和体散射。（1）在奇次散射和表面散射分量中,FD、YM4、AY4、AY3和YM3分解的重要性程度相对较高,在重要性排名中占据前5；而BF、VZ、CP分解相对其他分解方式的重要性程度有所下降,分别占据7~9名。（2）在体散射分量中,YM3、AY3和FD分解的重要性程度相对较高,分别排名12~14位,CP分解的重要性偏低,对土壤水分含量变化不太敏感。（3）二面角散射分量的重要性程度大部分低于体散射分量,但也存在例外,如VZ、AY3分解的二面角散射分量分别优于YM4、AY4的体散射分量,重要性排名分别为15、17位。（4）螺旋体散射分量的重要性偏低,原因是由于该分量的引入是为了解决城区的人造地物不满足反射对称问题[33],研究区内人造地物相对较少,因而其重要性偏低。（5）在H/A/α分解中,Lambda对土壤水分反演的贡献率较高,重要性排名第6,仅次于表面散射分量；RVI相对Lambda重要性程度有所下降,主要原因是RVI能够量化区域植被影响差异[7],但研究区植被较为稀疏,因此重要性偏低；其他参数的重要性排序依次为Entropy>Anisotropy>Alpha。（6）在4个后向散射系数中,同极化后向散射系数的重要性高于交叉极化的后向散射系数,说明同极化后向散射系数相较交叉极化后向散射系数更适合于稀疏植被,与文献[7]在冬小麦返青期的研究较为一致。
从图2（b）可见,光学特征参数的重要性排序依次为VSWI>COSRI>IPVI>LST>RVI>TVDI>MSAVI>NDVI>GVMI>GLI>EVI。VSWI重要性最高,EVI最低,这说明VSWI对于0~10 cm的土壤水分敏感度较高,能够较好地反映土壤水分的变化情况；EVI是为改善高植被区大气校正和土壤背景的影响而提出的[34],研究区多为低矮植被和裸地,因此EVI重要性偏低。MSAVI、NDVI、GVMI和GLI相比EVI重要性程度有所提升,但总体依然较低,可能是GVMI不太适合小尺度区域的干旱监测,NDVI和GLI适用于中等植被覆盖区,在稀疏植被和较高植被覆盖区,其土壤水分敏感性有所降低[35]。相比NDVI,MSAVI能在一定程度上改善土壤背景影响,因此重要性略高于NDVI。其余光学特征参数均对土壤水分反演具有较高的重要性。
3.2 模型精度对比
3.2.1 基于雷达特征参数的模型
在雷达特征参数中,剔除重要性排名靠后的6个参数,即对土壤水分含量变化敏感度较低的参数,将剩余的30个参数按重要性排序,从VIMR4（VIM表示重要性,下标R表示雷达特征参数,4表示前4个参数）开始,依次放入2个特征变量,分别计算R2和RMSE,选出最优特征组合方案,建立RF模型,各模型精度如图3所示。
![]()
图 3 雷达特征参数组合模型精度对比Figure 3. Accuracy Comparison of Radar Characteristic Parameter Combination Models由图3可见,基于雷达特征参数组合的14种模型训练集R2普遍高于0.8,RMSE普遍低于1.5%,验证集R2普遍高于0.5,RMSE普遍低于2.5%。在所有模型中,VIMR30的R2和RMSE并非最优,说明模型精度与输入特征数量并非简单的正相关关系,特征选择是提高模型精度的重要环节。从VIMR4至VIMR12,随着输入特征增加,精度显著提升,R2从0.54提升至0.64,RMSE从2.48%降低至2.26%；从VIMR12至VIMR18,虽然精度有所提升,但提升速度相比之前有所减缓,R2由0.64提升至0.67,RMSE由2.26%降低至2.16%；从VIMR18之后,R2和RMSE均呈现不同程度的下降,推测是前后因子间存在强相关性造成了数据冗余,使模型在训练时可能存在轻微过拟合现象,尽管如此,VIMR18的R2和RMSE分别达到0.67和2.16%。C波段极化雷达数据在植被覆盖地表土壤水分研究中应用较多,笔者团队前期对研究区Radarsat-2 C波段数据的研究表明,总体样本中各参数组合方案的RF模型验证集R2均低于0.6,RMSE在3.5%~5%之间[36]。本文基于L波段的参数方案相比C波段在精度上均有提升,尤其是VIMR18效果更佳,说明相比C波段而言,L波段全极化SAR对土壤水分敏感性更高。文献[37-38]对比了不同NDVI状态下ALOS-2 PALSAR-2 L波段雷达数据和Sentinel-1 C波段雷达数据在法国南部、突尼斯中部地区土壤水分反演中的敏感性,认为不同NDVI状态下L波段对土壤水分的敏感程度优于C波段,与本文的研究结果相吻合。因此,选取VIMR18作为最优雷达特征参数组合模型反演土壤水分。
3.2.2 基于光学特征参数的模型
在光学特征参数中,按重要性排名从VIMO4（下标O代表光学特征参数）开始,每次放入1个特征参数,将其作为RF模型的输入变量,以土壤水分含量作为输出变量,建立RF模型,以R2和RMSE评估模型精度,结果如图4所示。
![]()
图 4 光学特征参数组合模型精度对比Figure 4. Accuracy Comparison of Optical Characteristic Parameter Combination Models由图4可以看出,基于光学特征参数组合的8种模型训练集R2均在0.89左右,RMSE在1.25%左右,验证集R2在0.5左右,RMSE在2.47%左右。从VIMO4至VIMO11,验证集R2呈现先缓慢增长、后逐渐降低、再增长、又降低的趋势；RMSE则呈现先降低、后增长、再降低、又增长的趋势；从VIMO6至VIMO9,模型精度的第一次降低可能是由于因子间的强相关性导致模型训练集的过度拟合所致；从VIMO9至VIMO11,验证集R2呈现先增长、后下降的趋势,但对比VIMO9和VIMO11两组方案可知,VIMO11虽然R2有所提升,但RMSE也呈现上升趋势,可能是由于采样点土壤水分低值数量较多、高值数量过少,导致出现低值点高估、高值点低估现象,最终导致验证集RMSE有所上升。由于在光学特征参数组合方案中综合使用R2和RMSE无法确定最优组合方案,因此,选取验证集R2位居前3位的VIMO6、VIMO7、VIMO10光学特征参数组合方案进行土壤水分反演。
3.2.3 基于雷达⁃光学特征参数的协同模型
将上述优选的雷达和光学特征参数组合方案进一步组合,得到3种协同反演方案,将3种方案的反演值分别与验证点的土壤水分含量实测值进行对比分析,计算R2和RMSE,结果见表2。
表 2 雷达-光学参数协同反演模型精度对比Table 2. Accuracy Comparison of Radar-Optical Parameter Integrated Inversion Model方案 组合类型 R2 RMSE/% 方案1 VIMR18+ VIMO6 0.72 1.99 方案2 VIMR18+VIMO7 0.71 2.00 方案3 VIMR18+VIMO10 0.70 2.02 由表2可知,3种协同反演模型R2均不低于0.7,RMSE均低于2.05%。其中方案1的R2为0.72,RMSE为1.99%,在3种协同模型中精度最高。与仅使用雷达特征参数的最佳模型VIMR18相比,R2提高7.46%,RMSE降低8.54%,表明光学特征参数能够在一定程度上消除植被覆盖的影响,弥补雷达数据在土壤水分反演中的不足。与仅使用光学特征参数的模型VIMO6相比,R2提升38.4%,RMSE降低22.6%,说明相较于单一数据源而言,多源数据融合能够在一定程度上提升土壤水分反演精度,使预测结果更接近真实情况。相比方案1,方案2和方案3总体精度虽然有所下降,但仍然较为理想,进一步说明融合多源数据的机器学习方法在土壤水分反演中的应用潜力,其适用范围相较单一数据源更为广泛。
3.3 多种模型土壤水分反演结果对比
图5为基于雷达、光学以及雷达-光学特征参数协同反演模型的土壤水分反演结果。由图5可以看出,仅使用光学特征参数的模型（图5（b））存在水分整体高估的现象,特别是在裸地和部分低植被覆盖区,这些区域易受土壤背景和天气的影响[35],从而导致预测精度较低。相比图5（b）,基于雷达特征参数（图5（a））和雷达-光学特征参数协同（图5（c））的模型对于裸地和低植被覆盖区的估算效果较好,与实测水分含量基本一致,表明雷达特征参数在裸地和低植被区土壤水分反演中的应用潜力。已有研究表明,对于植被覆盖较高地区而言,仅使用雷达特征参数难以完全消除植被覆盖和地表粗糙度的影响,由图5（a）可以看出,东部植被覆盖区土壤水分含量相比实测值有所低估,加入光学特征参数后（图5（c））,低估现象有所改善,表明多源数据融合能在一定程度上克服因影响因子考虑不足导致的精度受限问题,这与文献[39]的研究结果较为吻合。但由于研究区土壤水分含量整体较低,导致土壤水分高值采样点较少,高值低估现象仍然存在。总体来说,图5能反映研究区土壤水分含量的整体趋势,其中雷达-光学协同模型反演效果最佳,能够更好地反映研究区土壤水分的空间分布格局。
4 结语
本文基于L波段ALOS-2 PALSAR-2和Landsat 8影像提取雷达特征参数和光学特征参数,通过MDA优选特征参数,建立基于雷达、光学以及雷达-光学特征参数协同的RF土壤水分反演模型,定量评估模型精度,反演绿洲土壤水分。结论如下：（1）相比C波段而言,L波段对干旱荒漠绿洲地区的土壤水分含量敏感性更高；（2）雷达特征参数中,表面散射分量重要性最高,体散射分量次之,二面角和螺旋体散射分量重要性较低,4个后向散射系数中,同极化后向散射系数的重要性优于交叉极化后向散射系数,光学特征参数中,VSWI重要性最高,EVI最低；（3）雷达特征参数方案最优模型R2、RMSE分别为0.67、2.16%。光学特征参数方案模型精度普遍较低且精度相当,R2、RMSE分别为0.5、2.47%；（4）相比单一数据源而言,协同利用雷达和光学特征参数,模型适用范围更广,预测精度更高。最优协同反演模型的R2、RMSE分别为0.72、1.99%。相比仅使用雷达特征参数模型VIMR18,R2提升了7.46%,RMSE降低了8.54%,相比仅使用光学特征参数模型VIMO6,R2提升了38.4%,RMSE降低了22.6%。
本文研究表明,雷达-光学协同的随机森林模型在干旱荒漠绿洲区土壤水分反演中具有良好的适用性,能够较好地反映区域土壤水分空间分布格局。后期将采用更多的机器学习方法,探索协同利用L波段全极化雷达数据和光学数据在较深层次土壤水分含量反演中的效果和适用性。
-
[1] [11]。1.2 谱图分割算法 谱图分割(spectral graph partition,SGP)是建立在谱图理论基础上的一种聚类算法,基于该理论,可以将影像分割转化为图的划分问题[6]:将待分割影像映射为加权无向图G={V,E}(V:顶点,E:边)。顶点对应影像的每个像元,带权重的边对应像元间的相似性。这样,将影像分割转化为图的最优划分问题。具体过程可以分为三步:(1)选择合适的影像特征(边缘、灰度、纹理等),建立像元间的相似性度量矩阵,完成影像到图的映射;(2)选择划分准则,将图划分为一定数量的子图;(3)将图的划分结果映射到影像空间,完成分割。因此,首先要提取PolSAR影像的边缘信息,获取分割线索。1.2.1 边缘信息提取 本文进行边缘提取的边缘检测器如图2所示,共包含4个方向的边缘检测,检测窗口可根据影像实际情况指定(3×3,5×5,7×7,…)。图2 边缘提取的检测器 Fig.2 Edge Map Calculation Detectors若以检测器中心为边缘,则两边区域应存在一定的差异性。这种差异性可基于 Wishart统计分布构建假设检验衡量: H0:Σi=Σj,H1:Σi≠Σj(2)式中,Σi、Σj分别为i、j区域的期望协方差矩阵。设Θi、Θj分别为i、j区域的样本协方差矩阵数据集。假设样本间相互独立,i、j区域间是否相似可以通过似然比检验确定。检验统计量[1]如下:LH0(^ΣΘi,Θj) ^ΣinNi^ΣjnNjQ=LH1 ^Σi,^ΣjΘi,Θj)= ^Σn(Ni+Nj)( (3)式中,L(·)为不同假设条件下的似然函数;^Σ、^Σi、^Σj分别为Σ、Σi、Σj的最大似然估计量;Ni、Nj分别为i、j区域的样本数量;n为多视视数。如果Q的值较低,则原假设H0将被拒绝。基于这个检验统计量定义两个区域的差异性度量为:D(Si,Sj)=-1lnQ=n(4)(Ni+Nj)ln^Σ-Niln^Σi-Njln^Σj 式中,Si、Sj分别代表i、j区域;i、j两个区域差异性越大,D(Si,Sj)也越大,边缘强度值也越大。采用式(4)计算每个像元上4个方向的D(Si,Sj),并保留Dmax和它的检测方向。在此基础上进行边缘优化,对于任意像元p,其边缘强度值为Dmax(p),边缘方向为θ*(p)。比较垂直于边缘方向两边像元的边缘强度值,如果Dmax(p)大于等于两边的值,则保留该值,否则,设置为零。1.2.2 相似性度量矩阵构建边缘信息作为一种有效的影像分割线索,在PolSAR影像谱图分割算法中已经有所应用[1,8]。但这些算法只适用于从像元出发的情况。本文将试图找到能够代表过分割区域空间位置的某一像元,来解决这一问题,计算方法如图3所示。图3 过分割区域中心计算Fig.3 Positioning of the Center of aOver-segmentation Object 图3中,灰色多边形为某过分割区域S,x为其内部像元,L(α)为方位角为α时像元x到该区域边界的距离,因此,像元x在过分割区域内的可扩展程度可表示为:n E(x∈S)=i∏=1L(α×i),n=2π/α (5) 根据式(5),计算区域S内的所有像元的E值,max{E(x∈S)}所对应的像元即为能够代表区域S空间位置的像元。 然后即可利用影像的边缘信息来衡量区域间的相似性程度,如图4所示。图4(a)为影像的Pauli RGB显示,S1、S2、S3为3个过分割区域,以及其相应的空间代表性像元p1、p2、p3。两个区域间的相似性度量通过代表性像元间的连线上最大的边缘强度值来体现。如图4(b)所示,p1和p3之间有明显的边缘存在,其相似性就会较低,即图4(b)中顶点间边的权重W(p1,p3)较低;同理,W(p1,p2)则较高。因此,定义差异性度量:DC(x,y)=D*(z*),z* =arg maxz∈lD*(z) (6)式中,D*(·)代表边缘强度值;l是连接像元x和y的直线;z*为该连线上具有最大边缘强度值的像元。通过高斯核函数映射即可完成权重的计算:W(x,y)=exp{-DC2(x2,y)} (7)式中,σC为核函数映射的尺度参数(2σC)。由式(7)便可完成相似性度量矩阵的构建。图4 基于边缘信息的区域间相似性权重计算Fig.4 Illustration of Extracting the DissimilarityInformation from Edge Maps1.2.3 归一化割准则 谱图分割的划分准则种类繁多,常见的划分准则及其优缺点可[6]。本文选择应用较为广泛的归一化割准则完成影像分割。假设将一个图划分为两个子图A和B,原来连接A、B后来被删去的边的权重之和在图论中称为割(cut): cut(A,B)=x∈∑A,y∈BW(x,y) (8)Shi和Malik[13]将式(8)除以表现顶点集大小的度量,完成割的归一化,即为归一化割准则:Ncut(A,B)=assoccut(A(A,B,V))+assoccut(B(B,A,V))=cut(A,B) cut(B,A)+x∈∑A,v∈VW(x,v) y∈∑B,v∈VW(y,v) (9)式中,V=A∪B;assoc(A,V)代表A中的所有顶点与该图顶点间权重的和。这时,图的划分准则即最小化归一化割Ncut(A,B)。1.3 分割评价方法 鉴于PolSAR影像特殊的统计模型,分割评价中常用的区域间对比度、区域内部一致性等评价指标不再适用。本文将以在遥感影像分割研究中应用较多的正确分割的百分数作为指标,评价PolSAR影像分割结果。陈晓秋等[14]通过累计待评价分割图中的每个区域内与参考分割图重叠的最大面积计算正确分割的百分数。该方法计算出的精度与区域个数相关,分割区域数越多,精度越高。因此,只有当分割区域数相当的情况下,上述分割精度的比较才有意义。但这种比较的公平性难以保障,而且分割精度无法同时反映过分割和欠分割的程度。本文对上述方法进行改进,计算正确分割的百分数。如图5所示,虚线区域f为参考图分割 区域,与其有重叠的待评价分割区域分别是a、b、c、d这4个区域。灰色区域e是面积最大的重叠区,其面积为Se。改进方法考虑了区域e在区域d中的面积比例(Se/Sd)。若该比例很小,则e应属于区域f的欠分割部分,相应像元属于错误分割,不应参与正确分割百分数的计算。因此,本文定义一个反映欠分割程度的参数,欠分割比例(under segmentation ratio,USR):max(Sref∩Seva(1,2,3…N))USR=1-Seva(i) (10)式中,Sref为某参考分割区域面积,Seva为待评价分割区域面积,N为重叠的Seva个数,Seva(i)为与Sref重叠最大的Seva。如图5所示,USR=1-Se/Sd。依据USR,即可通过限制欠分割程度计算正确分割的百分数。例如,设置 USR=0.2时,只有重叠区域最大且欠分割比例不大于0.2的像元才会被计入正确分割的像元数。这样计算出的分割精度对于影像过分割和欠分割的程度均有所反映,可以更客观地评价分割方法的优劣。图5 分割评价原理示意图Fig.5 Illustration Chart of SegmentationEvaluation Principle2 实验与分析2.1 实验数据 本文实验采用的数据为2013年6月16日获取的Radarsat-2全极化数据(C-波段),升轨,模式为FQ18,入射角37.56°,方位向、距离向分辨率分别为4.96m和4.73m。覆盖区域为内蒙古海拉尔垦区上库力农场依根区域,截取部分农田区域作为实验区域,影像多视处理后大小为312像元×292像元,图6为影像的Pauli RGB显示。 评价分割结果的参考图如图7所示,通过专图6 实验区Pauli RGB显示Fig.6 Experimental PolSAR Data in Pauli RGB 家知识和地面真实数据获得。主要步骤包括:(1)首先获取该地区同季节的光学影像,利用ArcGIS9.3等软件勾绘出地块边界;(2)在Radarsat-2过境当天,开展同步的土地覆盖类别调查,采用差分GPS设备根据现场实际情况修正勾绘边界,定位记录地块属性;(3)利用ArcGIS 9.3等软件将修正后的边界、地块属性等数据数字化,并进行适当后处理,例如将属于同一类但没有明显边界的地块合并,属于同一类但有明显边界的不予合并;(4)与Radarsat-2影像数据匹配,获得最终的参考分割图。该区域的土地覆盖类别主要包括油菜、小麦、灌木、裸露耕地和其他共5类。其中,裸露耕地因粗糙度、垄向、含水量等因素的影响表现为三种不同的极化反映,因此细分为三类。图7 分割评价参考图Fig.7 Reference for Segmentation Evaluation2.2 分割实验 经过多视化、Lee滤波(3×3)等数据预处理步骤,首先提取影像的边缘信息,检测器窗口设置为7×7。图8为边缘粗提取结果;图9为边缘优化结果。图10是不同参数设置下的均值漂移的分割结果,分割区域的块数分别为349、173及96。可以看出,单纯尺度的调整无法使均值漂移算法实现全局较优的分割,即便在最大尺度的情况下,图10中的椭圆黑框区域范围仍处于过分割的状态,而矩形黑框内已经出现了欠分割的情况,由此可见,分割过于破碎是均值漂移难以克服的缺点。 选择中等尺度的均值漂移(N=173)结果作为谱图分割的输入。代表性像元的空间位置如图11中黄色的点所示(α=π/4),可以看到,这些点均处于区域相对中心的位置。图8 边缘粗提取结果 图9 边缘优化结果 Fig.8 Raw Result Fig.9 Optimizedof Edge Map Result of Edge Map在上述基础上,构建过分割区域间的相似性度量矩阵,并通过归一化割完成分割。这一步骤中,可以通过设置最终的分割区域个数控制影像分割的尺度。图12是不同分割区域个数设置的分割结果,其分割细节的表现不同。当区域个数设置为42时,矩形黄色框内的两块裸土没有分割开,设置为48时则被分割开。其主要原因在于,两块裸土间边缘强度较弱,因此在区域个数设置较小时,它们将被合并为一个区域。 为了进一步验证本文方法的有效性,使用eCognition图像分割软件中的多尺度分割方法,对该实验区Pauli RGB影像进行了不同尺度的分割,参数设置及分割结果如图13所示。随着分割尺度的增大,区域个数逐渐减少,过分割现象得到有效缓解(椭圆黑框区域),但是同时,欠分割现象也在逐渐加剧(方形黑框区域)。这也说明了eCognition软件的多尺度分割与均值漂移类似(图10),在单一尺度设置上很难获得全局最优的分割结果。综合比较图12及图13的分割结果,相比之下,本文方法(图12(b))的分割结果最为优良。不同大小的农田块基本都得到了较好的分割,成为独立的对象,区域一致性好,这为下一步的影像解译奠定了良好的基础。 算法的运算代价,实验影像为312行×292列,共91 104个像元,若采用传统的基于像元的谱图分割算法,则需通过C291 104次运算建立91 104×91 104大小的相似性矩阵,而本文方法只需C2173次运算建立173×173大小的矩阵,运算空间及时间需求均降低约28万倍。由此可见,均值漂移预分割明显降低了谱图分割算法的运算代价,而谱图算法全局优化的策略也使最终的分割结果优于目前比较流行的算法。2.3 定量评价与分析 为了更客观地评价分割算法,选用定量的方法评价分割结果的优劣。首先,分别选取本文和eCognition多尺度方法较好的分割结果(图12(b)和图13(b))进行精度评价;然后,在其基础上,固定尺度外的参数不变,分析分割精度随分割区域个数的变化;最后,分析了均值漂移分割区域个数对分割精度的影响。 经计算,本文方法的分割精度为83.6%,eCognition方法为75.1%(USR=0.3)。图14是不同方法正确分割像元的分布图。其中白色区域为正确分割的像元。可以看出,eCognition方法的分割精度很大程度上受到分割过于破碎的影响。 分割精度随USR的变化趋势如图15所示。分割精度随 USR的增大而增大,但本文方法的分割精度始终大于eCognition方法的分割精度。USR=0时,即表明不允许任何欠分割现象,实际分割时几乎不可能出现,因此精度为0;USR=1时,表明对于欠分割不做任何控制,只要与参考图重叠面积最大,则计入正确分割的像元数。 图16是不同方法分割精度随分割区域个数的变化趋势(USR=0.3)。可以看出,峰值两侧的分割精度均呈下降趋势,这验证了本文改进的分割精度评价指标对过分割和欠分割现象均有反映,克服了以往文献[14]中评价指标的缺陷。分割精度的峰值高于eCognition方法,且在42～65范围内时,精度明显高于eCognition方法;而在曲线两端,eCognition方法的精度要高于本文方法,但这两侧都趋向于不理想的分割结果,精度也较低。 图17是本文方法分割精度随均值漂移过分割区域个数的变化趋势(其他参数设置与图12(b)相同)。可以看出,随着过分割区域个数的增加,分割精度整体呈增长趋势,虽然有一定的波动,但在区域个数达到一定数目时趋于稳定。这说明,在利用均值漂移进行过分割时,区域个数不宜过少,在区域个数达到一定数目时(大于320),本文方法会表现出较好的稳健性。图10 不同尺度的均值漂移分割结果Fig.10 Segmentation Results of Different Scales Based on Mean-Shift图11 过分割区域中心图12 不同尺度的谱图分割结果 Fig.11 Over-segmentation Regional CenterFig.12 Segmentation Results of Different Scales Based on SGP图13 eCognition多尺度分割结果(Sh:形状C:紧致度Sc:尺度N:区域个数)Fig.13 Segmentation Results of Multi-resolution Segmentation Approach Based on eCognition图14 不同分割方法正确分割像元数分布Fig.14 Distribution of Correct Segments Pixels图15 分割精度随USR变化趋势Based on Different MethodFig.15 Segmentation Accuracies Accordingto Different SUR图16 分割精度随分割结果区域个数变化趋势Fig.16 Segmentation Accuracies According toDifferent Region Number图17 分割精度随均值漂移过分割区域个数变化趋势Fig.17 Segmentation Accuracies of SGP Accordingto Different Over-segmentation Region NumberBased on Mean-Shift3 结论与展望 在谱图分割框架下,本文提出了一种新的基于均值漂移和谱图分割的PolSAR影像分割方法。利用Radarsat-2全极化数据,验证了提出的方法,基于改进的精度评价方法和分割结果参考图,对于本文方法和eCognition的多尺度分割方法进行了定量的精度评价。相关实验表明:(1)本文提出的算法有效实现了PolSAR影像分割,组合算法能够有效降低谱图分割算法的运算代价,得到理想的分割结果;(2)本文分割算法的分割精度要优于eCognition软件的多尺度分割方法,分割结果更趋近于参考图的分割结果;(3)本文提出的改进的分割精度评价方法,可以有效地评价分割结果的优劣,对于过分割和欠分割都能有所反映,可以独立于分割结果的区域个数评价分割精度;(4)在均值漂移过分割区域个数达到一定数目时,本文提出的算法具有较好的稳健性。本文提出的算法仍存在一定的不足,例如整个分割算法中最优参数选择,算法的适用性,以及是否可引入其他特征等方面,尚需进一步的研究。参 考 文 献[1] Liu B,Hu H,Wang H,et al.Superpixel-basedClassification with an Adaptive Number of Classesfor Polarimetric SAR Images[J].IEEE Trans.Geosci.Remote Sens.,2013,51(2):907-924[2] Zhang Jie.Segmentation of Polarimetric SyntheticAperture Radar[D].Qingdao:Shandong Universityof Science and Technology,2012(张杰.极化SAR影像的分割[D].青岛:山东科技大学,2012)[3] Yu Jie,Liu Zhenyu,Yan Qin,et al.SemiautomaticObject-oriented Classification of SAR Images onMultiscale[J].Geomatics and Information Scienceof Wuhan University,2013,38(3):253-257(余洁,刘振宇,燕琴,等.多尺度下的半自动面向对象SAR影像分类[J].武汉大学学报·信息科学版,2013,38(3):253-257)[4] Deng Shaoping,Li Pingxiang,Zhang Jixian,et al.Segmentation of Multi-polarized SAR Imagery Basedon the Theory of Evidence[J].Science of Survey-ing and Mapping,2011,36(6):32-37(邓少平,李平湘,张继贤,等.一种基于证据理论的多极化SAR图像分割方法[J].测绘科学,2011,36(6):32-37)[5] Yang Xin.Study on Segmentation and ClassificationAlgorithm of Polarimetric SAR Images[D].Cheng-du:University of Electronic Science and Technologyof China,2008(杨新.极化SAR图像的分割和分类算法研究[D].成都:电子科技大学,2008)[6] You Li.Research on Image Segmentation MethodBased on Spectral Clustering[D].Changsha:Na-tional University of Defense Technology,2011(由里.基于谱聚类的图像分割方法研究[D].长沙:国防科学技术大学,2011)[7] Yan Chengxin,Sang Nong,Zhang Tianxu.Surveyon Graph Theory Based Image Segmentation Tech-nique[J].Computer Engineering and Applica-tions,2006,5:11-15(闫成新,桑农,张天序.基于图论的图像分割研究进展[J].计算机工程与应用,2006,5:11-15)[8] Ersahin K,Cumming I G,Ward R K.Segmenta-tion and Classification of Polarimetric SAR Data U-sing Spectral Graph Partitioning[J].IEEE Trans.Geosci.Remote Sens.,2010,48(1):164-174[9] Fowlkes C,Belongie S,Chung F,et al.SpectralGrouping Using the Nystrm Method[J].IEEETrans.Pattern Anal Math Intell,2004,36(2):214-225[10] Ma Xiuli,Jiao Licheng.SAR Image SegmentationBased on Watershed and Spectral Clustering[J].Journal of Infrared and Millimeter Waves,2008,27(6):452-457(马秀丽,焦李成.基于分水岭-谱聚类的SAR图像分割[J].红外与毫米波学报,2008,27(6):452-457)[11] He W,Jaeger M,Reigber A,et al.Building Ex-traction from Polarimetricsar Data Using Mean Shiftand Conditional Random Fields[C].The 7th Euro-pean Conference on Synthetic Aperture Radar(EU-SAR),Friedrichshafen,Germany,2008[12] Zou Tongyuan,Yang Wen,Dai Dengxin,et al.AnUnsupervised Classification Method of PolSAR Im-age[J].Geomaticsand Information Science of Wu-han University,2009,34(8):90-95(邹同元,杨文,代登信,等.一种新的极化SAR图像非监督分类算法研究[J].武汉大学学报·信息科学版,2009,34(8):90-95)[13] Shi J,Malik J.Normalized Cuts and Image Seg-mentation[J].IEEE Trans.Pattern Anal.Mach.Intell.,2000,22(8):888-905[14] Chen Qiuxiao,Chen Shupeng,Zhou Chenghu.Seg-mentation Approach for Remote Sensing ImagesBased on Local Homogeneity Gradient and Its Eval-uation[J].Journal of Remote Sensing,2006,10(3):357-366(陈秋晓,陈述彭,周成虎.基于局域同质性梯度的遥感图像分割方法及其评价[J].遥感学报,2006,10(3):357-366)
下载:
计量
- 文章访问数: 1372
- HTML全文浏览量: 77
- PDF下载量: 295
