遥感图像在民用和军事领域都有非常广泛的应用,产生了巨大的社会和经济效益。根据单幅遥感影像估计高程信息可以应用于滑坡、泥石流等自然灾害的检测。我国是一个自然灾害多发的国家,因此,单幅遥感影像高程值提取在遥感图像处理中具有重要的应用价值。
现今,在国内外并没有见到有关单幅山脉遥感影像高程估计研究的报道,但是,有不少学者致力于遥感图像的建筑物高度的检测[1]。现在通过高分辨率遥感影像获取建筑物高度信息的方法主要有两种:第一种是利用立体像对进行立体测量[2-4];第二种是从单幅遥感影像入手以算法推导为主,研究地面阴影与高度的关系,或研究单幅影像成像姿态与地面高程信息的关系[5-7]。利用高分辨率遥感影像立体像对专门完成建筑物高度的量算以此实现对建筑物高度检测的应用很少,此领域基本侧重于城市三维建模和数字测图应用。由于建筑物是比较规则的地物,运用它的这一特性许多学者开始研究如何从单幅航片、高分辨率遥感图像中提取建筑物高度的算法,这也是当前研究的热点。单幅影像的立体测量基本都是从成像姿态和共线方程原理入手建立几何关系[7],或是从建筑物阴影的长度来求解其高度,至今未有更精确、更简便的方法。主要包括两种方法,一是利用成像姿态和共线方程求解建筑物高度[5, 7]。单张影像的立体量测基本都涉及卫星和太阳的高度角、方位角以及与建筑物的几何关系、影像的摄影姿态等参数,一般都需借助共线方程建立几何关系,最后完成建筑物高度的解算。另一种方法是利用阴影长度检测建筑物高度[6, 8-9]。这种方法工作量小,成本相对航空摄影测量要低,其突出优点是数据获取方便,能得到实时性强的数据。首先需要检测建筑物阴影的长度和方向,然后根据卫星、太阳的高度角和方位角等与阴影的几何关系求解建筑物高度。在这种方法中,研究者都作了无地形干扰、建筑物外部轮廓简单且垂直地表的假设。由于实际情况往往与此不符,所以此方法存在局限性。本文给出了一种基于暗通道原理的单幅山脉遥感图像高程估计方法。
2009年,He等[10]提出了一种运用暗通道先验去雾的算法,将暗通道原理应用于单幅图像去雾,取得了较好的效果;同时,在去雾的过程中,产生了准确的景深图像,此图像能够反映景物与镜头之间的距离。He利用Levin等[11]提出的soft matting过程计算得到优化透射图。由于soft matting过程的时间和空间复杂度均较高,近来有不少学者提出了一些改进
方法,例如基于加权暗通道的方法[12],此方法对图像的边缘和非边缘部分使用不同的权值,从而达到优化的目的;基于透射梯
度优先和多分辨率处理的方法[13];基于联合双边滤波的方法[14]等。本文将该算法应用于山脉遥感图像高程估计,得到了一个大致符合事实的高程信息图,但由于太阳光照、山脉高度等因素,会造成遥感图像局部地区有阴影,而阴影的存在,掩盖了该地区原本的高程信息,很大程度上影响了所得高程信息图的准确性。目前,还没有算法来解决这个问题,本文的目标就是对阴影部分进行处理[15-18],最终得到准确的高程信息。
1 基于暗通道的高程信息提取 1.1 暗通道原理山脉遥感图像是卫星从太空向地面拍摄的,得到的图像是经过大气中悬浮颗粒吸收和散射之后的退化图像。这些图像清晰度较差、颜色失真,不利于图像特征提取及实际的应用;但是,可以利用图像中的雾霾来提取山脉的高程信息。图 1中,给出了本文算法的流程。
类似于计算机视觉和计算机图形学中的雾天图像退化模型[19],在遥感图像研究中,常用如下公式描述遥感图像退化模型:
$\mathit{\boldsymbol{I}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}})=\mathit{\boldsymbol{J}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}})\mathit{\boldsymbol{t}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}})+\mathit{\boldsymbol{A}}(1-t(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}}))$ | (1) |
其中:I(x, y)为拍摄到的遥感图像,J(x, y)为未退化的图像,t(x, y)为透射图即景深图像,它代表了未被散射从而到达相机的那部分光线,即离镜头较近的景物会保存更多的光线到达镜头,在景深图中显示为亮度较大,反之离镜头较远的景物会损失更多的光线,从而只有少量到达镜头,在景深图中显示为亮度较小。在遥感图像中,较高的地物可以看作是离镜头(卫星)较近,较低的地物就是离镜头(卫星)较远,这与自然图像的景深图完全吻合,所以t(x, y)可以被看作高程信息图。A为全局大气光向量。高程信息提取的目的就是通过I(x, y)求得中间量A,最终求得高程信息图t(x, y)。
结合上文分析,图 2以一幅雾天图像的实验来进一步形象地解释式(1) :
暗通道先验是基于无雾户外图像的观察:在绝大多数非天空的局部区域里,某一些像素总会有至少一个颜色通道具有很低的值,即这些小的图像块所对应的最小灰度值近似于0。本文给暗通道一个数学定义,对于任意的输入图像J,其暗通道可以用下式表示:
${{\bf{J}}^{{\rm{dark}}}}({\bf{x}},{\bf{y}}) = \mathop {\min }\limits_{({x_0},{y_0}) \in \Omega ({\bf{x}},{\bf{y}})} {\mkern 1mu} (\mathop {\min }\limits_{c \in \{ r,g,b\} } {\mkern 1mu} {{\bf{J}}^c}({{\bf{x}}_0},{{\bf{y}}_0})) \to 0$ | (2) |
其中:Jc是J的一个颜色通道,Ω(x, y)是一个以(x, y)为中心的小图像块。
按照暗通道优先理论,对于退化图像I,在雾较浓的区域,其暗通道图像的灰度值也会明显增加,因此退化图像的暗通道图像的灰度值可近似为雾的浓度,A的值近似等于雾最浓区域的值。假设已知空气光向量A,对式(1) 进行简单变换后可得:
${\bf{\tilde t}}({\bf{x}},{\bf{y}}) = 1 - \omega \mathop {\min }\limits_{({{\rm{x}}_{\rm{0}}},{{\rm{y}}_{\rm{0}}}) \in \Omega ({\bf{x}},{\bf{y}})} {\mkern 1mu} \left( {\mathop {\min }\limits_{{\rm{c}} \in \{ {\rm{r}},{\rm{g}},{\rm{b}}\} } {\mkern 1mu} \frac{{{{\bf{I}}^{\rm{c}}}({{\bf{x}}_{\rm{0}}},{{\bf{y}}_{\rm{0}}})}}{{{{\bf{A}}^{\rm{c}}}}}} \right)$ | (3) |
其中:ω为常量参数。根据式(3) 可得初始高程图。在自然环境下,远处的景物相比近处的景物会稍显模糊。如果完全去除雾气,图像会显得不自然,同时会丢失高程信息,因此添加ω来保存较远景物处的少许雾气。
1.1.3 估计空气光向量A求取A方法如下:
1) 根据式(2) 计算得到暗通道图像Jdark;
2) 从Jdark中提取前0.1%亮度最高的像素点;
3) 在原图像I中提取步骤2) 中得到的点中平均强度最大的点,将该点的R、G、B值赋值给A。
1.1.4 soft matting由于暗通道处理时使用了分块的方法,使得初始高程图存在比较明显的方块,导致不能较好地保持图像边缘,因此,He等[10]利用Levin等[11]提出的soft matting过程计算得到优化高程图。本文中,直接采用He的方法[10]来改进高程图。
图 3中,(a)为原有雾图像,(b)是求得r、g、b三通道最小值的图像,(c)是经以Ω(x, y)为模板的最小值滤波得到的暗通道图像,(d)是初始景深图,(e)为经soft matting优化的景深图。
图 4中,图 4(a)是山脉遥感图像,其中灰色虚线表示山脊,白色虚线表示山谷。可以注意到,在卫星拍摄的山脉遥感图像中,大多有阴影区域的存在。本图中,山脊线左侧存在大量阴影区域。而这些阴影区域,严重影响了暗通道原理在提取遥感图像高程信息时的准确性。图 4(b)是相应的高程图,颜色越红高程越大,颜色越蓝高程越小。可以明显看到,山脊线左右两侧的颜色分布并不对称,左侧阴影区域颜色偏红,即左侧高于右侧,这与事实不符。
在暗通道原理下得到的高程信息图,受阴影的影响。由于山脉遥感图像中的阴影与现实中其他一般图像的阴影有本质区别,一般图像中的阴影可以用边缘检测算法确定[20-21],但山脉遥感图像中的阴影边界受山脊、山谷的影响,并不能用边缘检测算法确定。
本文采用思想如下:
1) 经过对大量图像的观察,通过大量实验对阴影提出一个阈值T1,即r、g、b三通道的值均小于T1,则判定该点为阴影点;
2) 在判定为阴影的区域中,存在一些比较小而稀疏的部分,根据先验知识,知道它们并不是山脉阴影区域,因此,根据实验结果,对其面积S设置一个阈值T2,若S>T2,则保留此阴影区域;
3) 保留下的阴影区域还有一个缺陷,内部存在缝隙,本文使用形态学闭运算,先膨胀后腐蚀,尽可能去除阴影内部空隙,得到相对最佳的阴影模板,如图 5(a)所示。
1.2.2 去除阴影在已经确定了阴影区域的情况下,可以采用对阴影部分调整亮度的方法来获取原高程信息。本文方法建立在这样一个假设上,阴影只是掩盖了该地区原有的信息,一旦调整亮度,原有的信息将会重现。
值得注意的是,即使在同一幅图像上,山脉阴影的程度也会不同,因此,对不同的区域,要使用不同的参数。本文把图像中除阴影区域的平均亮度作为本图像的亮度标准,对于每块阴影区域,根据亮度差值与自身亮度的比值来确定参数,即亮度较大的区域,调整参数较小;亮度较小的区域,则调整参数较大。阴影参数根据式(4) 计算:
$\mathit{\boldsymbol{ }}\!\!\alpha\!\!\rm{ }(i)=\frac{{{l}_{u}}-{{\mathit{\boldsymbol{l}}}_{s}}(i)}{{{\mathit{\boldsymbol{l}}}_{s}}(i)}$ | (4) |
其中:i表示各阴影块的序号,ls(i)表示第i块阴影的平均亮度,lu表示图像中除阴影块的其他面积亮度,α(i)表示每块阴影区域的调整参数。
在得到调整参数后,对阴影区域进行调整,得到调整后图像,再对其应用暗通道原理,得到本文的结果,如图 5(b)所示。图中,山脊线两边的颜色分布趋于一致,与事实相符。
2 实验结果与分析 2.1 定性分析本实验的图像来自百度地图,截取山脉地区的卫星图,以及相对应的地形图作为参考。本文总共做了10组实验,取其中3组。
图 6是3组实验的结果,第一列是原遥感山脉图像,第二列是该地区的地形图,第三列是本文算法的结果(其中,右侧颜色条表示从红色到蓝色,高程越来越小), 同时,每一行表示一组实验。
实验1的结果如图 6第一行;对照第一列和第二列可知,图正中间为山脉的山脊线,其周围海拔应较高,图上部和下部各有一条山谷线,其周围海拔应较低;同时,本文截取的是北半球的山脉图,山脊线上方的区域会有阴影,经本文算法处理,去除阴影影响,则山脊线两边高程应相似。第三列图中,山脊线周围偏红,表示高程大,并且山脊线两边颜色分布相似,山谷线周围偏蓝,表示高程小,这与实际符合。
实验2的结果如图 6第二行;图像中间部分从左往右分别是山谷线、山脊线、山谷线;由于这是东半球的山脉图,所以山脊线左方的区域会有阴影。结果中,山脊线周围偏红,并且山脊线两边颜色分布相似,山谷线周围偏蓝,相对高程正确。
实验3的结果如图 6第三行;图正中间为山脉的山谷线,其周围海拔应较低,图上部和下部海拔应较高;又因为该山脉处于北半球,所以山谷线下方的区域会有阴影。第三列图中,山谷线周围偏蓝,表示高程小,并且山谷线两边颜色分布相似。
2.2 定量分析本实验的图像来自Google地球,高程数据来自谷地地理信息系统,截取山脉地区的卫星图,并对图像进行采样,获取对应的高程信息。本文总共做了20组实验,其中10组取采
样点,10组取采样面,并各取3组(实验4~6) 进行展示。
2.2.1 设置采样点分析对实验4,图 7(a)为山脉遥感图像;图 7(b)为本文算法的结果,亮度越大,高程越大;图 7(c)显示的是实验4的8个采样点,从山脊开始,均匀分布到山谷,再从山谷分布到另一山脊,则高程应符合这样一个趋势,高程从大到小,再从小到大。表 1实验4展示了这8个采样点的位置和高程,其中实际排序与本文排序这两列完全吻合。图 7(d)为8个山脉采样点的真实高程数据和本文结果的折线图,其中,横坐标对应的是各个采样点,纵坐标是采样点的高程(真实高程数据和本文的结果均已归一化),红色表示真实高程数据,蓝色表示本文的结果,从图中可以得到,两组折线的起伏趋势一致。
图 7(c)中:实验5的7个采样点,从山谷开始,均匀分布到山脊,再从山脊分布到另一山谷,那么它的高程信息应符合从小到大,再从大到小;实验6的8个采样点,点3、4、5、7在山谷部分,高程较小,点1、2、6、8在山脊部分,高程较大。表 1实验5和实验6分别展示了这些采样点的位置和高程。图 7(d)中,从两组折线易知,山脉采样点的真实高程数据和本文结果大小一致。
2.2.2 采样面定量分析在定量分析之前,先进行预处理,将采样面和本文算法的结果量化为10个等级(即亮度为10个等级,1表示最低,10表示最高)再进行比较,下面给出了定量评价的指标Ratio:
$Ratio=\frac{\mathit{\boldsymbol{GT}}\oplus \mathit{\boldsymbol{R}}}{m\times n}$ | (5) |
其中,GT表示该地区的真实相对高程,R表示本文算法的结果,⊕运算表示GT和R中高度相同的点(误差在一个亮度级),m×n表示采样面的总像素点。
图 8是三组实验的结果,(a)是原遥感山脉图像,(b)是采样面数据(其中,黑色表示地势低,白色表示地势高),(c)是本文算法的结果,(d)是GT和R的相似度图(其中,白色表示相同,黑色表示不同), 同时,每一行表示一组实验。
实验7见图 8第一行,由图(b)可以看到此区域中间为山谷,两边较高,同时图(c)也满足这样的地势,使用式(5) 求得Ratio=0.535 4;实验8见图 8第二行,观察地势可知图中左下为山谷,右上为山脊,图(c)的明暗也符合这一趋势,并求得Ratio=0.578 6;实验9见图 8第三行,地势恰与实验7相反,图中中间为山脊,四周较低,图(c)也是中间高,四周低,并求得Ratio=0.523 6。
3 结语本文提出了一种基于暗通道的单幅山脉遥感影像高程提取方法,并对山脉的阴影部分进行了改进。首先确定山脉遥感图像的阴影区域,并计算各阴影区域的调整参数,调整亮度达到去除阴影的目的;然后计算大气光向量,并使用soft matting,以生成优化后的高程信息图。实验结果表明,本文的方法能够有效地提取单幅山脉遥感图像的相对高程信息,同时克服山脉阴影的影响,但是,在处理阴影方面还有欠缺:第一,阴影检测可能会出现偏差,即不能准确定位阴影区域;第二,去除阴影的时候,只考虑亮度的调整效果一般,未来研究中还可以考虑饱和度、对比度等。
[1] | 周雨薇, 陈强, 孙权森, 等. 结合暗通道原理和双边滤波的遥感图像增强[J]. 中国图象图形学报, 2014, 19 (2) : 313-321. ( ZHOU Y W, CHEN Q, SUN Q S, et al. Remote sensing image enhancement based on dark channel prior and bilateral filtering[J]. Journal of Image and Graphics, 2014, 19 (2) : 313-321. ) |
[2] | 高翔, 赵冬玲, 张蔚. 利用高分辨率遥感影像获取建筑物高度信息方法的分析[J]. 测绘通报, 2008 (3) : 41-43. ( GAO X, ZHAO D L, ZHANG W. On the methods of obtaining the building height information from high-resolution remote sensing images[J]. Bulletin of Surveying and Mapping, 2008 (3) : 41-43. ) |
[3] | 徐玮, 高辉, 张茂军, 等. 折反射全景与遥感图像融合的建筑物高程自动提取方法[J]. 国防科技大学学报, 2011, 33 (2) : 81-88. ( XU W, GAO H, ZHANG M J, et al. The method of automatic building height extraction by fusing catadioptric panoramas and remote sensing images[J]. Journal of National University of Defense Technology, 2011, 33 (2) : 81-88. ) |
[4] | 王媛媛, 陈旺, 张茂军, 等. 折反射全向图像与遥感图像配准的建筑物高度提取算法[J]. 计算机应用, 2011, 31 (9) : 2477-2480. ( WANG Y Y, CHEN W, ZHANG M J, et al. Method for estimating building heights via registering catadioptric omnidirectional image and remote sensing image[J]. Journal of Computer Applications, 2011, 31 (9) : 2477-2480. ) |
[5] | COMBER A, UMEZAKI M, ZHOU R, et al. Using shadows in high-resolution imagery to determine building height[J]. Remote Sensing Letters, 2012, 3 (7) : 551-556. doi: 10.1080/01431161.2011.635161 |
[6] | CHENG F, THIEL K H. Delimiting the building heights in a city from the shadow in a panchromatic SPOT-image-Part 1:test of forty-two buildings[J]. Remote Sensing, 1995, 16 (3) : 409-415. doi: 10.1080/01431169508954409 |
[7] | GUIDA R, IODICE A, RICCIO D. Height retrieval of isolated buildings from single high-resolution SAR images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48 (7) : 2967-2979. doi: 10.1109/TGRS.2010.2041460 |
[8] | 王继周, 林宗坚, 李成名. GIS信息辅助的单影像立体量测[J]. 测绘科学, 2005, 30 (2) : 60-63. ( WANG J Z, LIN Z J, LI C M. GIS-supported 3D information acquisition of buildings from a single image[J]. Science of Surveying and Mapping, 2005, 30 (2) : 60-63. ) |
[9] | DURIEUX L, LAGABRIELLE E, NELSON A. A method for monitoring building construction in urban sprawl areas using object-based analysis of Spot 5 images and existing GIS data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2008, 63 (4) : 399-408. doi: 10.1016/j.isprsjprs.2008.01.005 |
[10] | HE K, SUN J, TANG X. Single image haze removal using dark channel prior[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33 (12) : 2341-2353. doi: 10.1109/TPAMI.2010.168 |
[11] | LEVIN A, LISCHINSKI D, WEISS Y. A closed-form solution to natural image matting[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30 (2) : 228-242. doi: 10.1109/TPAMI.2007.1177 |
[12] | 刘楠, 程咏梅, 赵永强. 基于加权暗通道的图像去雾方法[J]. 光子学报, 2012, 41 (3) : 320-325. ( LIU N, CHENG Y M, ZHAO Y Q. An image dehazing method based on weighted dark channel prior[J]. Acta Photonica Sinica, 2012, 41 (3) : 320-325. doi: 10.3788/gzxb ) |
[13] | 胡伟, 袁国栋, 董朝, 等. 基于暗通道优先的单幅图像去雾新方法[J]. 计算机研究与发展, 2010, 47 (12) : 2132-2140. ( HU W, YUAN G D, DONG Z, et al. Improved single image dehazing using dark channel prior[J]. Journal of Computer Research and Development, 2010, 47 (12) : 2132-2140. ) |
[14] | 方帅, 王峰, 占吉清, 等. 单幅雾天图像的同步去噪与复原[J]. 模式识别与人工智能, 2012, 25 (1) : 136-142. ( FANG S, WANG F, ZHAN J Q, et al. Simultaneous dehazing and denoising of single hazing image[J]. Pattern Recognition and Artificial Intelligence, 2012, 25 (1) : 136-142. ) |
[15] | LI Y, BROWN M S. Single image layer separation using relative smoothness[C]//Proceedings of the 2014 IEEE Conference on Computer Vision and Pattern Recognition. Washington, DC:IEEE Computer Society, 2014:2752-2759. |
[16] | LIU M, SALZMANN M, HE X. Discrete-continuous depth estimation from a single image[C]//Proceedings of the 2014 IEEE Conference on Computer Vision and Pattern Recognition. Washington, DC:IEEE Computer Society, 2014:716-723. |
[17] | 田建东, 王占鹏, 唐延东. 静态阴影检测的研究进展[J]. 信息与控制, 2015, 44 (2) : 215-222. ( TIAN J D, WANG Z P, TANG Y D. Static shadow detection:a survey[J]. Information and Control, 2015, 44 (2) : 215-222. ) |
[18] | 邓琳, 邓明镜, 张力树. 高分辨率遥感影像阴影检测与补偿方法优化[J]. 遥感技术与应用, 2015, 30 (2) : 277-284. ( DENG L, DENG M J, ZHANG L S. Optimization of shadow detection and compensation method for high-resolution remote sensing images[J]. Remote Sensing Technology and Application, 2015, 30 (2) : 277-284. ) |
[19] | NARASIMHAN S G, NAYAR S K. Vision and the atmosphere[J]. International Journal of Computer Vision, 2002, 48 (3) : 233-254. doi: 10.1023/A:1016328200723 |
[20] | LIOW Y T, PAVLIDIS T. Use of shadows for extracting buildings in aerial images[J]. Computer Vision, Graphics, and Image Processing, 1990, 49 (2) : 242-277. doi: 10.1016/0734-189X(90)90139-M |
[21] | LINDEBERG T. Edge detection and ridge detection with automatic scale selection[J]. International Journal of Computer Vision, 1998, 30 (2) : 117-156. doi: 10.1023/A:1008097225773 |