论文部分内容阅读
摘要：针对当前地震后震区时空数据难以及时获取的问题，采用了一种改进均值漂移（Mean Shift）信息提取方法。首先，将原始影像划分为纹理区和均色区。均色区域直接利用Mean Shift算法获得；纹理区域则利用归一化分布密度值获取合适的带宽，再使用Mean Shift算法进行信息提取。通过构造代价函数判别相邻区域是否需要合并，以消除过分割区域。最后，提出了一种信息提取匹配指数对信息提取结果进行评价，并将该方法应用于汶川和芦山地震后获取的无人机高空间分辨率影像，进行信息提取实验。实验结果表明：所提出的改进Mean Shift算法提取精度优于传统Mean Shift算法，为后续的地震災情评估提供了基础数据。
关键词：地震灾区；均值漂移；信息提取；高分辨率影像
中图分类号：P315-39文献标识码：A文章编号：1000-0666（2016）03-0494-06
0引言
我国是一个自然灾害频发且种类繁多的国家，如何有效监测灾情、及时有效地处理综合信息并做出科学可行的救灾应急措施，这些都将成为政府日常管理的重要组成部分。对高分辨遥感影像处理的研究，需要把传统的灾害监测和防灾减灾技术与其他科技手段、特别是高新科技结合起来，综合应用遥感、测绘、地理和空间信息技术等多学科技术手段，为早日实现空、地、人的立体监测网和综合信息处理系统做准备（李德仁等，2012）。针对高分辨率遥感影像的信息提取算法研究已成为当前遥感应用研究领域的一个重要课题。
2013年4月20日四川省雅安市芦山县（303°N，10295°E）发生70级地震，震源深度13 km。震后灾区路断、桥断、通信中断，重灾区陷入信息隔绝状态，这是继2008年5月12日汶川大地震后四川省内发生的又一大地震。在地震灾区通信、交通被严重破坏的情况下，卫星遥感和无人机遥感技术成为快速获取灾情的最佳途径。在汶川地震、芦山地震这样的巨大自然灾害面前，如何快速有效地处理影像数据并提取震害信息，使其在灾后应急及灾情评估与灾后重建中发挥其应有的作用是对遥感技术应用提出的挑战，也是值得深入研究的问题。
随着计算机技术的发展，Mean Shift算法的应用越来越广泛，在空间目标聚类、动态目标跟踪、视觉图像处理等领域均取得了良好的效果（周芳芳等，2007）。由于高空间分辨率无人机影像具有光谱信息较弱、纹理信息丰富等特点，所以直接将Mean Shift算法应用于高空间分辨率遥感影像信息提取中很难得到预期的结果（鲁恒等，2015）。基于上述原因，本文提出了一种改进Mean Shift算法进行信息提取，以期能在地震发生后快速获取震后第一手空间数据。
1Mean Shift算法
当前大多数聚类算法普遍需要先验知识，但Mean Shift算法不需要先验条件，初始点可以选择数据集任意点。该算法在执行过程中，类别根据收敛点确定，面向各种采样点的分布和维度都能实现聚类。Mean shift算法是一个基于密度函数的梯度估计和模式搜索的迭代过程，从数据集合的某个点出发，通过计算均值向量，不断逼近数据集合分布的局部极值点。Fukunaga和Hostetler（1996）最早把Mean Shift算法引入到模式识别领域，Cheng（1995）改进了算法中的核函数，扩大了该算法的适用范围。Comaniciu和Meer（1999，2002）将该方法应用于图像特征空间的分析，证明了Mean Shift算法在满足一定条件下，可收敛到最近的一个概率密度函数的稳态点。密度函数梯度估计等于密度函数估计的梯度，即
21影像滤波
遥感影像主要包括两种噪音，即非高斯噪声和高斯噪声（李峰，彭嘉雄，2003；肖鹏峰等，2007），这可能导致影像信息提取出现错误。因此，在进行影像信息提取前需利用滤波技术进行过滤，笔者采用自适应平滑滤波。通过求解像元周围3×3矩形窗内8邻域像元3个色彩通道的方差和均值，相对于各个色彩通道，在像元减去均值绝对值大于12倍方差时，就认为该像元为噪声点。同时，在像元8邻域内像元值均采用均值代替。再通过色彩空间变化的方法处理滤波影像，并且采用（L，U，V）模型表示色度值、亮度值。
22影像区域分类
相同性质遥感影像区域表现平滑、灰度与色彩之间保持良好的一致性。因此，根据一致性就能将遥感影像划分成纹理区域和匀色区域。同时，还可通过简单分类处理方法，即中心直接采用x0点，在宽度W的矩形窗口范围，通过x0点色彩特征矢量（L，U，V）分布密度判断邻域色彩一致性，其中，所采用的计算公式为
σx0=kΣWgxi-x0h.（3）
核函数采用多元正态分布形式：
g（x）=（2π）-d/2exp-12‖x‖2.（4）
由式（4）可以看出，x0处平滑程度σx0数值成正比。在进行求解像元σxi值时，采用归一化处理得到σxi′。当σxi′趋于1时就说明平滑程度越高。同时，进行归纳所有σxi′，取定阈值T=09，将原影像转换为二值影像：当σxi′T，像元值为1，说明该区域为纹理区域。
23带宽可变的Mean Shift聚类
一般通过Mean Shift算法完成区域提取，匀色区域通常情况下倾向于使用较大带宽参数，能够进行连续和完整的提取。直接利用固定带宽为10的Mean Shift方法提取匀色区域。在提取变化剧烈区域时，根据地貌特征应选择较小的Mean Shift带宽，保证分析特征中充分考虑丰富的纹理信息。所求σxi′建立可变带宽的函数为
h=h01σxi′. （5）
其中，h0为基准Mean Shift带宽。经实验研究h0=12。
24区域合并
为了防止纹理区域出现的过分割问题，笔者采用区域合并方法解决该问题。首先构造代价函数，通过合并相邻区域，产生新区域，改变前后两者的变化程度通过特征向量来反映。若假设相邻区域ri和rj，合并后为rij，特征分量变化可采用以下计算公式： Δij=nijcij-（nici+njcj）.（6）
式中，nij=ni+nj表示合并区域rij包含的像元数。在经过加权平均方式计算出9维的向量空间（L，u，v，μ，σ，Sk，Ku，s，c）下的合并代价，具体计算公式为
C（ri，rj）=Σ9f=1kfΔij.（7）
其中，kf为各维向量的加权系数，且Σ9k=1=1。
然后进行邻接图合并，区域邻接图合并需要根据所有邻接图进行寻优，确定最佳的合并區域，保证合并代价最小。对于相邻区域ri和rj满足
rj=argminC（ri，rk），rk∈{R，eim∈E}.（8）
从式（8）可以看出， rj是ri合并代价最小的邻接区域，反过来ri也满足，即ri和rj分别是对方邻接区域中合并代价最小的，则ri和rj是满足局部最优合并条件的区域对。
[BT1-*8]3实验及结果分析
为验证笔者所提方法，在Matlab平台下编程实现基于改进的Mean Shift算法及传统的Mean Shift信息提取算法，将两种方法分别应用于所获取的芦山地震及汶川地震后震区无人机影像，通过两组实验来验证算法的有效性。所获取的芦山和汶川地震后原始无人机影像如图1所示。
传统Mean Shift提取算法使用带宽为04时的提取结果如图2a、c所示，使用改进Mean Shift提取算法的结果如图2b、d所示。对比图2a、b可发现，本文所提算法较之传统Mean Shift算法在影像中间建筑物聚集区提取边界更精细、准确，图中左上角的林地提取更为准确。对比图2c、d可发现，本文所提算法道路提取更为精细完整，建筑用地提取更加清晰、完整，碎片较少。总体来看，改进算法的提取结果与遥感影像中的实际地物有很好的对应关系。
为了定量验证方法信息提取的精度，需要建立评价指标。Neubert和Herold（2008）通过提取对象与实际地物的面积匹配度来判断分割尺度是否合理；于欢等（2010）通过计算分割对象边界与实际地物边界横向和纵向矢量距离来评判提取结果的好坏。上述匹配模型均忽略了分割对象与实际地物间的光谱匹配情况。笔者通过光谱和面积匹配情况进行综合考虑，建立了一种顾及面积和光谱的信息提取匹配指数EMI：
EMI=Hli-Hri·Ali-AriAri.（9）
式中，Ali代表最大提取对象的面积，对于实际地物A，提取出的对象A1、A2、A3都占有实际地物A的部分像元，但A1所切割A的面积最大，因此A1为最大提取对象，Ari表示作为参考的实际地物的面积，Hli为最大提取对象的灰度均值，Hri为作为参考的实际地物的灰度均值。以图2a、b提取出的建筑物为例说明计算方法，通过目视解译，随机在无人机影像选择150个实际地物数据作为样本参考数据，参考地物选择边界明显的建设用地、耕地、道路，如图3所示。
在获得参与评价对象的信息中，计算出所有
从表1可知，针对芦山、汶川两期震后无人机影像信息提取结果计算出来的MEMI值，本文所提算法均小于传统Mean Shift算法，验证了本文所提算法的可靠性。
4结论
本文从高分辨率遥感影像自身特点出发，对传统Mean Shift算法进行改进，取得了较为满意的信息提取结果。由于地震后获取的高分辨遥感影像纹理信息丰富，在使用传统Mean Shift算法进行信息提取时容易产生过分割和欠分割现象。改进的分割算法将影像划分为匀色区域和纹理区域，对不同性质的区域使用带宽可变的Mean Shift分割，改善了信息提取结果，保证震后无人机遥感影像分割的准确性，为灾情评估提供了基本的数据源。
参考文献：
李德仁， 眭海刚， 单杰.2012. 论地理国情监测的技术支撑[J]. 武汉大学学报（信息科学版）， 37（5）：505-512.
李峰， 彭嘉雄. 2003. 基于Wold特征的遥感图像的分割[J]. 测绘学报， 32（4）： 330-333.
鲁恒， 李龙国， 贺一楠. 2015. 考虑地形特征的无人机影像分区域加权平差拼接方法[J]. 农业机械学报，46（9）：296-301.
肖鹏峰，冯学智，赵书河等. 2007. 基于相位一致的高分辨率遥感图像分割方法[J].测绘学报，36（2）：146-151.
于欢，张树清，孔博等. 2010. 面向对象遥感影像分类的最优分割尺度选择研究[J].中国图象图形学报，15（2）：352-360.
周芳芳， 樊晓平， 叶榛. 2007. 均值漂移算法的研究与应用[J]. 控制与决策，22（8）：841-847.
关键词：地震灾区；均值漂移；信息提取；高分辨率影像
中图分类号：P315-39文献标识码：A文章编号：1000-0666（2016）03-0494-06
0引言
我国是一个自然灾害频发且种类繁多的国家，如何有效监测灾情、及时有效地处理综合信息并做出科学可行的救灾应急措施，这些都将成为政府日常管理的重要组成部分。对高分辨遥感影像处理的研究，需要把传统的灾害监测和防灾减灾技术与其他科技手段、特别是高新科技结合起来，综合应用遥感、测绘、地理和空间信息技术等多学科技术手段，为早日实现空、地、人的立体监测网和综合信息处理系统做准备（李德仁等，2012）。针对高分辨率遥感影像的信息提取算法研究已成为当前遥感应用研究领域的一个重要课题。
2013年4月20日四川省雅安市芦山县（303°N，10295°E）发生70级地震，震源深度13 km。震后灾区路断、桥断、通信中断，重灾区陷入信息隔绝状态，这是继2008年5月12日汶川大地震后四川省内发生的又一大地震。在地震灾区通信、交通被严重破坏的情况下，卫星遥感和无人机遥感技术成为快速获取灾情的最佳途径。在汶川地震、芦山地震这样的巨大自然灾害面前，如何快速有效地处理影像数据并提取震害信息，使其在灾后应急及灾情评估与灾后重建中发挥其应有的作用是对遥感技术应用提出的挑战，也是值得深入研究的问题。
随着计算机技术的发展，Mean Shift算法的应用越来越广泛，在空间目标聚类、动态目标跟踪、视觉图像处理等领域均取得了良好的效果（周芳芳等，2007）。由于高空间分辨率无人机影像具有光谱信息较弱、纹理信息丰富等特点，所以直接将Mean Shift算法应用于高空间分辨率遥感影像信息提取中很难得到预期的结果（鲁恒等，2015）。基于上述原因，本文提出了一种改进Mean Shift算法进行信息提取，以期能在地震发生后快速获取震后第一手空间数据。
1Mean Shift算法
当前大多数聚类算法普遍需要先验知识，但Mean Shift算法不需要先验条件，初始点可以选择数据集任意点。该算法在执行过程中，类别根据收敛点确定，面向各种采样点的分布和维度都能实现聚类。Mean shift算法是一个基于密度函数的梯度估计和模式搜索的迭代过程，从数据集合的某个点出发，通过计算均值向量，不断逼近数据集合分布的局部极值点。Fukunaga和Hostetler（1996）最早把Mean Shift算法引入到模式识别领域，Cheng（1995）改进了算法中的核函数，扩大了该算法的适用范围。Comaniciu和Meer（1999，2002）将该方法应用于图像特征空间的分析，证明了Mean Shift算法在满足一定条件下，可收敛到最近的一个概率密度函数的稳态点。密度函数梯度估计等于密度函数估计的梯度，即
21影像滤波
遥感影像主要包括两种噪音，即非高斯噪声和高斯噪声（李峰，彭嘉雄，2003；肖鹏峰等，2007），这可能导致影像信息提取出现错误。因此，在进行影像信息提取前需利用滤波技术进行过滤，笔者采用自适应平滑滤波。通过求解像元周围3×3矩形窗内8邻域像元3个色彩通道的方差和均值，相对于各个色彩通道，在像元减去均值绝对值大于12倍方差时，就认为该像元为噪声点。同时，在像元8邻域内像元值均采用均值代替。再通过色彩空间变化的方法处理滤波影像，并且采用（L，U，V）模型表示色度值、亮度值。
22影像区域分类
相同性质遥感影像区域表现平滑、灰度与色彩之间保持良好的一致性。因此，根据一致性就能将遥感影像划分成纹理区域和匀色区域。同时，还可通过简单分类处理方法，即中心直接采用x0点，在宽度W的矩形窗口范围，通过x0点色彩特征矢量（L，U，V）分布密度判断邻域色彩一致性，其中，所采用的计算公式为
σx0=kΣWgxi-x0h.（3）
核函数采用多元正态分布形式：
g（x）=（2π）-d/2exp-12‖x‖2.（4）
由式（4）可以看出，x0处平滑程度σx0数值成正比。在进行求解像元σxi值时，采用归一化处理得到σxi′。当σxi′趋于1时就说明平滑程度越高。同时，进行归纳所有σxi′，取定阈值T=09，将原影像转换为二值影像：当σxi′
23带宽可变的Mean Shift聚类
一般通过Mean Shift算法完成区域提取，匀色区域通常情况下倾向于使用较大带宽参数，能够进行连续和完整的提取。直接利用固定带宽为10的Mean Shift方法提取匀色区域。在提取变化剧烈区域时，根据地貌特征应选择较小的Mean Shift带宽，保证分析特征中充分考虑丰富的纹理信息。所求σxi′建立可变带宽的函数为
h=h01σxi′. （5）
其中，h0为基准Mean Shift带宽。经实验研究h0=12。
24区域合并
为了防止纹理区域出现的过分割问题，笔者采用区域合并方法解决该问题。首先构造代价函数，通过合并相邻区域，产生新区域，改变前后两者的变化程度通过特征向量来反映。若假设相邻区域ri和rj，合并后为rij，特征分量变化可采用以下计算公式： Δij=nijcij-（nici+njcj）.（6）
式中，nij=ni+nj表示合并区域rij包含的像元数。在经过加权平均方式计算出9维的向量空间（L，u，v，μ，σ，Sk，Ku，s，c）下的合并代价，具体计算公式为
C（ri，rj）=Σ9f=1kfΔij.（7）
其中，kf为各维向量的加权系数，且Σ9k=1=1。
然后进行邻接图合并，区域邻接图合并需要根据所有邻接图进行寻优，确定最佳的合并區域，保证合并代价最小。对于相邻区域ri和rj满足
rj=argminC（ri，rk），rk∈{R，eim∈E}.（8）
从式（8）可以看出， rj是ri合并代价最小的邻接区域，反过来ri也满足，即ri和rj分别是对方邻接区域中合并代价最小的，则ri和rj是满足局部最优合并条件的区域对。
[BT1-*8]3实验及结果分析
为验证笔者所提方法，在Matlab平台下编程实现基于改进的Mean Shift算法及传统的Mean Shift信息提取算法，将两种方法分别应用于所获取的芦山地震及汶川地震后震区无人机影像，通过两组实验来验证算法的有效性。所获取的芦山和汶川地震后原始无人机影像如图1所示。
传统Mean Shift提取算法使用带宽为04时的提取结果如图2a、c所示，使用改进Mean Shift提取算法的结果如图2b、d所示。对比图2a、b可发现，本文所提算法较之传统Mean Shift算法在影像中间建筑物聚集区提取边界更精细、准确，图中左上角的林地提取更为准确。对比图2c、d可发现，本文所提算法道路提取更为精细完整，建筑用地提取更加清晰、完整，碎片较少。总体来看，改进算法的提取结果与遥感影像中的实际地物有很好的对应关系。
为了定量验证方法信息提取的精度，需要建立评价指标。Neubert和Herold（2008）通过提取对象与实际地物的面积匹配度来判断分割尺度是否合理；于欢等（2010）通过计算分割对象边界与实际地物边界横向和纵向矢量距离来评判提取结果的好坏。上述匹配模型均忽略了分割对象与实际地物间的光谱匹配情况。笔者通过光谱和面积匹配情况进行综合考虑，建立了一种顾及面积和光谱的信息提取匹配指数EMI：
EMI=Hli-Hri·Ali-AriAri.（9）
式中，Ali代表最大提取对象的面积，对于实际地物A，提取出的对象A1、A2、A3都占有实际地物A的部分像元，但A1所切割A的面积最大，因此A1为最大提取对象，Ari表示作为参考的实际地物的面积，Hli为最大提取对象的灰度均值，Hri为作为参考的实际地物的灰度均值。以图2a、b提取出的建筑物为例说明计算方法，通过目视解译，随机在无人机影像选择150个实际地物数据作为样本参考数据，参考地物选择边界明显的建设用地、耕地、道路，如图3所示。
在获得参与评价对象的信息中，计算出所有
从表1可知，针对芦山、汶川两期震后无人机影像信息提取结果计算出来的MEMI值，本文所提算法均小于传统Mean Shift算法，验证了本文所提算法的可靠性。
4结论
本文从高分辨率遥感影像自身特点出发，对传统Mean Shift算法进行改进，取得了较为满意的信息提取结果。由于地震后获取的高分辨遥感影像纹理信息丰富，在使用传统Mean Shift算法进行信息提取时容易产生过分割和欠分割现象。改进的分割算法将影像划分为匀色区域和纹理区域，对不同性质的区域使用带宽可变的Mean Shift分割，改善了信息提取结果，保证震后无人机遥感影像分割的准确性，为灾情评估提供了基本的数据源。
参考文献：
李德仁， 眭海刚， 单杰.2012. 论地理国情监测的技术支撑[J]. 武汉大学学报（信息科学版）， 37（5）：505-512.
李峰， 彭嘉雄. 2003. 基于Wold特征的遥感图像的分割[J]. 测绘学报， 32（4）： 330-333.
鲁恒， 李龙国， 贺一楠. 2015. 考虑地形特征的无人机影像分区域加权平差拼接方法[J]. 农业机械学报，46（9）：296-301.
肖鹏峰，冯学智，赵书河等. 2007. 基于相位一致的高分辨率遥感图像分割方法[J].测绘学报，36（2）：146-151.
于欢，张树清，孔博等. 2010. 面向对象遥感影像分类的最优分割尺度选择研究[J].中国图象图形学报，15（2）：352-360.
周芳芳， 樊晓平， 叶榛. 2007. 均值漂移算法的研究与应用[J]. 控制与决策，22（8）：841-847.