基金项目:全时程灾情综合分析与决策技术研究(2018YFC1504503),国家重点研发计划专项(2018YFC1504403),国家重点研发计划(2018YFC1504503)及地震行业科研专项结余资金使用计划项目联合资助.
(中国地震局地质研究所,北京 100029)
(Institute of Geology,China Earthquake Administration,Beijing 100029,China)
pre-assessment of earthquake disaster losses; UAV; surface buildings extraction; morphological imageprocessing
备注
基金项目:全时程灾情综合分析与决策技术研究(2018YFC1504503),国家重点研发计划专项(2018YFC1504403),国家重点研发计划(2018YFC1504503)及地震行业科研专项结余资金使用计划项目联合资助.
选取新疆阿图什市下辖的琼哈拉峻村为研究区,将小型旋翼无人机拍摄的图像作为数据源,分别采用面向对象以及面向像元2种影像分析方法对研究区的房屋建筑进行提取,并对2种算法的提取结果进行对比,分析了各自的优势。结果 表明:面向对象方法可以有效地去除椒盐噪声对分类带来的影响,保证房屋形态的完整性,但影响内部相似的光谱、纹理信息若对应多种物体则会导致影像对象的错分。在面向像元的提取方法中加入了改进的数学形态学算法,可有效的抑制椒盐噪声,保持建筑物边缘的连续性与完整性,较好地解决了面向对象方法中部分农田与房屋出现错分的问题。
Taking Qiongharajun Village in Atushi City,Xinjiang as the research area,using the image taken by the small-scale rotorcraft UAV as the data source,we extract the housing buildings in the research area by object-oriented and pixel-oriented image analysis methods,compared the results of the two algorithms,and analyzed their respective advantages. The results show that the object-oriented method can effectively eliminate the impact of salt and pepper noise on classification and ensure the integrity of the housing morphology. However, if the similar spectral and texture information in the image corresponds to a variety of ground objects, the image object will be misclassified. The improved mathematical morphology algorithm is added to the pixel-oriented extraction method, which can effectively suppress salt and pepper noise, maintain the continuity and integrity of the building edge, and solve the problem of some farmland and houses in the object-oriented method.
引言
地震灾害损失预评估是震后及时、科学地进行地震应急指挥的关键之一(安基文等,2015; 聂高众等,2002)。地震导致人员伤亡的因素中有80%是房屋破坏造成的。在进行地震灾害损失预评估的过程中,关键要做好评估区域房屋建筑状况,以及人口规模的调查(聂高众等,2012)。而对于研究区域的房屋建筑提取工作则是上述调查的有力保证。
由于无人机遥感技术相较于星载遥感技术拥有诸多优势(Berra et al,2016; Dandois,Ellis,2013; Messinger et al,2016),在地震灾害损失预评估的实际工作中无人机逐渐开始代替传统的星载遥感技术对地表建筑物信息进行获取(付博,2018; 臧克等,2010; 李祖传等,2010; 彭大雷等,2017; 周洋等,2017; 陈晋等,2018)。在调查过程中每一个区域的调查点繁多,且每一个调查点的调查时间有限,所携带的大多数为飞行准备时间短、操作灵活、易于携带的小型旋翼无人机,而大多数小型无人机所搭载的传感器为数字相机,拍摄得到的影像虽然图像分辨率高,但是光谱分辨率很低,即所返回的影像数据波段数较少。本文利用面向对象以及面向像元2种影像分析方法对研究区的房屋建筑进行提取,并对2种算法的提取结果进行对比,分析了在该数据条件下2种方法各自的优势。
1 数据来源及影像拼接
本文选取的研究区为新疆阿图什市下辖的哈拉峻乡琼哈拉峻村(40.178°N,76.832°E)。哈拉峻乡90 m空间分辨率数字高程模型(Digital Elevation Model,DEM)叠加相应的山地阴影Hillshade图,如图1a所示; 无人机飞行路线规划如图1b所示。
本文使用小型四旋翼无人机,为大江精灵3专业版,相机有效像素1 240万,在新疆琼哈拉峻村由软件控制自动飞行获取影像数据。为了便于对拍摄得到的影片进行拼接,保证拍摄时一定的航向与旁向重叠度(程效军等,2005),设置旁向重叠度为75%,航向重叠度为50%,通过飞行拍摄得到69幅高分辨率航拍影片。
实际飞行时将无人机飞行高度(H)设置为相对地面100 m,由于无人机相机焦距(h=4 mm)和感光元件尺寸(d=0.001 58 mm)已知,参考公式:
R=(H×d)/h(1)
式中:R为地面空间分辨率。无人机拍摄影像的地面空间分辨率约为3.9 cm,优于航天平台可见光遥感图像的空间分辨率,达到高分辨率遥感影像的标准。
采用北京天地智绘科技有限公司开发的EasyUAV软件对拍摄得到的影像进行拼接,该软件采用的图像拼接技术主要分为3个步骤:
(1)通过图像预处理对图像中几何畸变进行矫正,使参考图像和待拼接图像不存在明显的几何畸变。本次飞行拍摄得到的影像为正射影像,其最常见的几何畸变是由像主点或中心点向外扩展时光轴倾斜度增大、视野变宽、比例尺缩小所带来的。由于本次飞行的飞行高度较低(100 m),单幅图像的像主点与边缘点相差的实际距离不大,在处理该几何畸变时利用畸变的图像提取阴影信息,沿着图像失真的逆过程恢复图像的原貌。实际的复原过程是设计一个反向滤波器,使其能从失真图像中计算得到失真图像的估值,根据规定的误差准则,最大程度地接近真实图像。预处理过程中还可以对图像噪声进行抑制,去除数据传输中带来的不可预测的随机信号,但这种随机信号极易形成图像上的噪声点。噪声点的去除通常由软件在图像输出时完成,本次输出的图像分辨率较高,人为观察没有发现明显的噪声点,如果此时再用传统的高斯滤波、均值滤波、中值滤波等低通滤波器对图像进行二次处理,反而会使图像的分辨率降低,降低提取精度。
(2)图像配准主要对参考图像和待拼接图像中的匹配信息进行提取,在提取出的信息中寻找最佳的匹配,完成图像之间的对齐。由于无人机拍摄影像分辨率较高,拍摄地区房屋数量较多,所以房屋的角点可以很好地作为配准点。使用逐一配准法对图像进行配准,其原理为:在搜索图S中以某点作为基点(搜索图S的长宽分别为M,N),截取一个与模板T(长宽分别为U,V)大小一样的分块图像,这样的分块图像有(M-U+1)×(N-V+1)个,配准的目标就是在(M-U+1)×(N-V+1)个分块图像中找一个与待配准图像最相似的图像,这样的基准点就是最佳配准点。
(3)在完成图像配准后对图像进行缝合,对缝合边界进行平滑处理,让缝合自然过渡,最终生成研究区DOM影像(图2)。从图2中可以看出,
图1 哈拉峻乡Hillshade图(a)和无人机飞行路线规划图(b)
Fig.1 Hillshade map of Halajun Township(a)and UAV flight route planning map(b)由于拍摄时间在冬季,新疆阿图什市下辖琼哈拉峻村内农田、草地等其他地物的光谱、纹理特征与房屋较为相似,光谱分辨率低极易出现同谱异物现象(田新光,2007; 陈济才等,2014),而这种情况可能会对接下来的分类产生一定影响。
2 研究方法
3 实验过程
3.1 面向对象的影像分析方法本文采用的图像分割方法如下:基于前文所列公式计算得到每一次分割时每个影像对象的异质性为F,利用基于异质性最小原则的区域合并算法(即多尺度分割的思想)进行图像的分割。首先将单个像元合并为一个个较小的对象,然后小的对象可以经过若干步骤逐渐合并成大的对象,合并过程必须保证影像对象的异质性小于给定的阈值。
在实际的分割过程中,由于房屋建筑的光谱特征差异性较小,设置光谱因子的权重为0.2,形状因子的权重为0.8,由于目视观察发现房屋屋顶纹理较为平滑,设置光滑度权重为0.7,紧密度权重0.3,经多次尝试阈值s设置在350 em为最佳。以任意一个像元为中心开始分割,第一次分割结束后由式(2)计算相邻像元的异质性F是否小于给定的阈值s,如果不满足阈值条件,用上一次分割出的影像对象为基础进行第二次分割,重复上述步骤直到相邻2个影像对象的异质性小于给定阈值s,完成最后的分割。利用多尺度分割的分割结果如图6a所示。
图6 图像分割结果(a)和基于面向对象的分类结果(b)
Fig.6 Image segmentation results(a)and object-oriented classification results(b)由图6a可以看出,在分割过程加入了紧致度以及表示屋顶的光滑度因子可以相对有效地对研究区不同地物进行分割,在计算异质性的过程中弱化光谱信息在分割中所占权重,提高形状因子所占权重,使得每一个影像对象的边缘较为完整。对分割后的图像提取每个影像对象3个波段的灰度均值以及标准差作为光谱特征,利用每个对象内部的灰度共生矩阵计算ASM与IDM值,以1:1加权求和作为纹理特征。ASM值反映了每一个对象内灰度分布均匀程度和纹理粗细度,计算公式如下:
ASM=∑ki=1∑kj=1[G(i, j)]2(18)
IDM反映图像纹理的同质性,度量图像纹理局部变化的多少,计算公式如下:
IDM=∑ki=1∑kj=1(G(i, j))/(1+(i-j)2)(19)
利用提取出的特征向量对分割结果进行分类,如图6b所示,其中绿色表示农田,蓝色表示房屋建筑,黄色表示道路,黑色表示阴影。
3.2 面向像元的形态学处理的算法实现首先对原始图像进行图像增强处理,提取原始图像每一点的灰度值做频数直方图,如图7a所示,在图中可以看出该图像的灰度分布集中分布在50~255,0~50除了0灰度存在多数像元(图像处理之后处于边缘的无效值),其他灰度值几乎无像元分布。针对这种情况采用分段线性拉伸的方法对图像进行增强,将图像中灰度50~255的区域进行拉伸。原始图像中房屋的灰度值主要集中在250~255区域且有一个明显的波峰。利用灰度拉伸后的原始影像作为底图,提取250~255灰度区域的图像作为前景图像,对前景图像二值化,结果如图7b所示。
在前景图像的提取过程中,容易提取出较多非房屋类型的像元,而用简单的去噪方法无法有效地去除这种大面积噪音,还可能误将有效的房屋像元一同去掉。鉴于这一问题,本文利用数学形态学图像处理的方法去除图像中复杂的背景,具体操作如下:①考虑到无关像元均为细碎斑点,较少数连续的无关像元面积也均远小于房屋面积这一事实,采用长水平线作为结构元B进行重建的开操作。在重建过程中,图像首先被腐蚀,
图7 原始影像的频数分布直方图(a)和前景图像(b)
Fig.7 Frequency distribution histogram of original image(a)and foreground image(b)如同标准的形态学开操作一样,然而这个被腐蚀的图像在图像重建过程中被用作标记图像F,原始图像作为模板G重新执行形态学重建过程:将标记图像F初始化为h1(标记F必须是G的子集,也就是FG)。建立结构元B=ones(1,11)定义连通性,本文使用8连通。重复hk+1=(hkB)I G,直到hk+1=hk,完成重建的开操作。②对重建后的图像执行顶帽重建,即从原始图像中减去重建的开操作,执行结果记为f-thr。③将结构元换成长竖直线B=ones(11,1),对原图像进行顶帽重建,重建结果记为f-thi。④以f-thr作为模板,以min(f-thi,f-thr)作为标记,最后执行形态学重建,作为最终的去噪结果。对重建之后的图像用canny算子进行边缘检测,其中阈值向量T设定为0.04~0.1,平滑滤波器的标准差sigma设定为0.1。
对图像中房屋的内部孔洞、裂隙进行处理:本文基于重建的膨胀方法,利用迭代的方式进行开操作,准确地恢复腐蚀之后的物体形状。在实际操作中,由于研究区内的房屋多为矩形,且房屋面积均大于5 m2,因此在执行开操作时选择矩形结构元且阈值设为5。对图像进行迭代的开操作,将房屋中孔洞及裂缝进行填补,对房屋内部的形态进行重建,之后同样利用长宽分别为5的矩形结构元对图像进行迭代的形态学闭操作,对房屋外部轮廓进行重建,去除细的突出,使房屋边界平整。最后对形态学重建后的图像进行填充。整体算法的操作流程如图8所示。加入形态学处理后的边缘检测结果与最终的提取结果如图9所示。
3.3 精度评价地表真实数据的获取方法如下:以数字正射影像图为底图对研究区进行矢量化构图操作,利用目视判读的方式提取出房屋所在位置以及面积信息,生成Shapefile文件,并以此文件作为模板在研究区域截取相关房屋数据,对截取出的图像进行二值化。为了对研究中提取的分类结果进行精度评价,将截取出的图像进行重分类,将分类后的图像(图 10)作为真实的房屋位置与面积数据。
利用混淆矩阵对分类结果进行评价:分别将面向对象的房屋提取结果与数学形态学方法进行房屋提取的结果作为分类数据,将矢量化后的房屋位置及面积信息(图 10)作为地表真实数据建立混淆矩阵,计算制图精度、用户精度、漏分误差以及错分误差,计算结果如表1所示。
4 讨论与结论
在新疆琼哈拉峻村内通过小型无人机获取影像数据,分别利用面向对象和面向像元的方法进行分类,从面向对象的影像分类结果可以看出:面向对象方法可以有效地去除椒盐噪声对分类带来的影响; 在分割过程中合理地加入地物的形态、纹理信息可以使分割后每一个对象的边缘保持连续性,这也保证了影像分类时房屋形态的完整性。但是,在结果中出现了部分区域错分的情况。观察错分区域的正射影像可以发现:错分区域多为农田,该区域不仅光谱特征与房屋相似,还拥有与房屋屋顶相似的纹理以及矩形形状。由于面向对象的分类是基于每一个影像对象的分类,而在影像的多尺度分割时只能考虑相邻影像对象之间的异质性差异,不能考虑非相邻对象之间可能会出现特征相似的现象。研究区农田与房屋存在多数特征相似的情况,若使用面向对象的方式提取影像对象内部特征进行分类,则容易导致影像对象的错分。
利用面向像元的方法分离房屋则会存在提取的地物完整性不高、易出现椒盐噪声的问题。在实验过程中首先需要考虑如何在不影响真实房屋像元的情况下剔除噪声,且在房屋提取过程中应尽量保持房屋的形态完整。本文在前人提出的数学形态学方法上做出改进,较好地实现了像元级房屋的提取,解决了面向对象方法中出现的由于影像对象内部特征相似导致的部分农田与房屋出现错分现象这一问题。
2.1 面向对象的影像分析方法面向对象影像分析方法将同质性相近的临近像元进行分割合并成影像对象,以影像对象作为分析处理的基本单元,对影像进行地物信息的提取(Baatz,2000)。本文利用图像中的光谱因子与形态因子计算该图像中影像对象之间的异质性(其中形态因子由光滑度和紧致度决定),计算公式为:
F=w·hcolor+(1-w)·hshape(2)
式中:F为异质性; w为光谱权重;(1-w)为形状权重。
光谱异质性hcolor的计算公式为:
hcolor=∑cwc[SMerge·σMergec-(Sobj1· σobj1c+Sobj2· σobj2c)](3)
式中:c为图层标识; wc为参与图像分割的图层权重; SMerge为合并后对象面积; σMergec为新生成的对象在图层c中的光谱标准差; Sobj1,Sobj2分别为待合并的2个相邻对象的面积; σobj1c,σobj2c分别为这2个对象在图层c中的光谱标准差。
形状异质性hshape主要由紧致度hcompact和光滑度hsmooth2个方面组成,hcompact和hsmooth取决于构成对象的像素个数:
hshape=wcompact·hcompact+(1-wcompact)·hcompact(4)
hcompact=nMerge·(lMerge)/((nMerge)1/2)-(nobj1·(lobj1)/((nobj1)1/2)+nobj2·(lobj2)/((nobj2)1/2)(5)
hsmooth=nMerge·(lMerge)/(nMerge)-(nobj1·(lMerge)/(nMerge)+nobj2·(lMerge)/(nMerge))(6)
式中:wcompact为紧致度权重;(1-wcompact)为光滑度权重; nMerge,nobj1,nobj2分别表示合并对象的最小外包正方形、待合并的2个相邻对象的最小外包正方形; lMerge,lobj1,lobj2分别表示这些对象的边长。
对利用以上公式计算得到的异质性与所给阈值进行对比,作为判断某区域是否合并的标识。
2.2 面向像元的形态学处理技术2.2.1 基于重建的膨胀与腐蚀方法膨胀是求局部最大值的操作,对于前景物体中一些细小的断裂处,如果小于结构元素的大小,这些断裂的地方就会被连接起来。而腐蚀与膨胀相反,是求局部最小值的操作。对于前景物体中一些细小的连接处,如果小于结构元素的大小,这些连接处就会被断开。膨胀与腐蚀的原理如图3所示。A被B膨胀的数学表达式可记为AB,作为集合操作,定义为:
AB={z|(B^)z∩A≠Ø}(7)
式中:Ø为空集; B为结构元; B^表示B的反射,z表示把B的反射用z平移,可以表示为z(zx,zy),x和y表示平移的横纵坐标,等式右边可以解释为B的反射被所有z平移后与A至少有一个非零公共元素。
A被B腐蚀的数学表达式可记为AB,作为集合操作,定义为:
AB={z|(B)zA}(8)
等式右边可以解释为用z平移的B包含在A中的所有的点z的集合。
本文首先定义一种基于重建的膨胀方法,为下述的开闭操作进行准备。定义F为标记图像,记录膨胀的起点,G为模板图像用来对膨胀进行约束且FG,则一次迭代的膨胀过程可以定义为:
D(1)G(F)=(FB)∩G(9)
则F关于G的n次迭代膨胀可定义为:
D(n)G(F)=D(1)G[D(n-1)G(F)](10)
在该递推公式中,集合求交在每一步都执行,由于每一次迭代模板G会限制标记F的生长,所以经过有限数量的迭代步骤后n总会收敛:
D(n)G(F)=D(n+1)G(F)(11)
基于重建的膨胀算法原理如图4所示。基于重建的腐蚀原理与上述过程类似。定义F为标记图像,记录腐蚀的起点,G为模板图像用来对膨胀进行约束且FG,则一次迭代的膨胀过程可以定义为:
D(1)G(F)=(FB)∪G(12)
则F关于G的n次迭代腐蚀可定义为:
D(n)G(F)=D(1)G[D(n-1)G(F)](13)
式中:并集操作在每一次迭代中必须被执行,保证每一次腐蚀过后的图像仍然大于等于模板图像,由于每一次迭代模板G会限制标记F的生长,所以经过有限数量的迭代步骤后n总会收敛:
D(n)G(F)=D(n+1)G(F)(14)
2.2.2 基于迭代的开操作与闭操作在传统的形态学开操作中,腐蚀典型地去除小的物体,且随后的膨胀趋向于恢复保留的物体形状。然而,这种恢复的精确度取决于形状和结构元之间的相似性。本文在基于重建的膨胀与腐蚀的基础上,利用迭代的方式进行开操作能够准确地恢复腐蚀之后的物体形状。
A被B形态学开操作表示为A B,定义为A被B腐蚀,然后再用B膨胀腐蚀结果:
A B=(AB)B(15)
几何学上,A B是B在A中完全匹配平移的并集。图5a显示了集合A和圆形结构元B,图5b显示了在A内完全匹配B的平移,所有这些平移的并集即图5c为形态学开操作的结果,剩下被割去的区域就是结构元不能完全在A中匹配的区域。形态学开操作除去了所有不能包含结构元的部分,平滑目标的轮廓,断开了细的连接部分,去除细的突出。A被B形态学闭操作表示为A·B,定义为A被B膨胀,然后再用B腐蚀膨胀结果:
A·B=(AB)B(16)
几何学上,A·B执行所有不与A重叠的B平移的补集。图5d显示了一些与A不重叠的B的平移,通过完成以上平移的并集操作,可以得到图5e所示的区域,该区域代表着完成闭操作后生成的区域。像开操作一样,形态学闭操作趋向于平滑物体的轮廓。然而,不同于开操作的是,闭操作一般连接窄的断裂并填满细长的缝隙,填满比结构元小的洞。
本文所使用的迭代开操作可以定义为:记F为标记图像,记录腐蚀的起点,G为模板图像用来对膨胀进行约束且FG,在每一次迭代过程中结构元B先对模板进行腐蚀,之后再对腐蚀后的图像进行膨胀重建。可以表示为:
O(n)G(F)=D(n)G[(G)nB](17)
而对于闭操作,对上述开操作的结果求补即为闭操作结果。
- 安基文,徐敬海,聂高众,等.2015.高精度承灾体数据支撑的地震灾情快速评估[J].地震地质,37(4):1225-1241.
- 陈济才,文学虎,李国明.2014.基于面向对象的高分影像地表覆盖典型要素快速提取对比研究[J].遥感信息,(4):37-40.
- 陈晋,习聪望,陈文凯,等.2018.基于无人机,高分卫星遥感影像的甘肃省陇南市建筑物空间化研究[J].地震研究,41(2):192-200.
- 程效军,朱鲤,刘俊领.2005.基于数字摄影测量技术的三维建模[J],同济大学学报(自然科学版),33(1):37-41.
- 付博.2018.基于无人机正射影像的建筑物震害识别研究[D].北京:中国地震局地质研究所.
- 李祖传,马建文,张睿,等.2010.利用融合纹理与形态特征进行地震倒塌房屋信息自动提取[J].武汉大学学报(信息科学版),35(4):446-450.
- 聂高众,安基文,邓砚.2012.地震应急灾情服务进展[J].地震地质,34(4):782-791.
- 聂高众,高建国,马宗晋,等.2002.中国未来10~15年地震灾害的风险评估[J].自然灾害学报,11(1):68-73.
- 彭大雷,许强,董秀军,等.2017.无人机低空摄影测量在黄土滑坡调查评估中的应用[J].地球科学进展,32(3):319-330.
- 田新光.2007.面向对象高分辨率遥感影像信息提取[D].北京:中国测绘科学研究院.
- 臧克,孙永华,李京,等.2010.微型无人机遥感系统在汶川地震中的应用[J].自然灾害学报,19(3):162-166.
- 周洋,明小娜,杨艳珠,等.2017.灾评新技术在云龙5.0级地震烈度调查中的应用[J].地震研究,40(1):161-166.
- Baatz M.2000.Multiresolution segmentation:an optimization approach for high quality multi-scale image segmentation[J].Angewandte geographische informations verarbeitung,XII,Karlsruhe,12-23.
- Berra E F,Gaulton R,Barr S.2016.Use of a digital camera onboard a UAV to monitor spring phenology at individual tree level[C].Geoscience and Remote Sensing Symposium(IGARSS),3496-3499.
- Dandois J P,Ellis E C.2013.High spatial resolution three-dimensional mapping of vegetation spectral dynamics using computer vision[J].Remote Sensing of Environment,136:259-276.
- Messinger M,Asner G P,Silman M.2016.Rapid assessments of Amazon forest structure and biomass using small unmanned aerial systems[J].Remote Sensing,8(8):615-615.