面向大规模真正射影像生产的快速遮蔽检测方法

钟 成，李 卉，黄先锋

（1.中国地质大学 教育部长江三峡库区地质灾害研究中心，湖北 武汉 430074；2.中国地质大学 行星科学研究所，湖北 武汉 430074；

3.武汉大学 测绘遥感信息工程国家重点实验室，湖北 武汉 430079）

传统正射纠正未考虑建筑物、树木等地面目标的投影变形，在城市大比例尺正射影像上，建筑物倾斜遮挡的情况非常严重，影响了影像信息的有效解译和利用.20世纪90年代以来，随着世界范围内城市三维仿真项目的兴起和城市大比例尺地图的应用需求日益强烈，利用城市三维模型制作“真正射影像”的研究和生产逐渐展开.

真正射影像制作的核心问题是高效稳健的遮蔽检测过程.根据摄影测量的基本原理，非正射投影造成的高层建筑倾斜会导致对低矮目标的遮蔽，若直接利用数字微分纠正方法为DSM（数字表面模型）格网寻找对应影像信息，低矮目标将获得高层建筑的影像信息，其结果是造成高层建筑在正射影像上多次出现，学术界称之为“鬼影”现象.这样的正射影像无法谈及几何精度或光谱精度等问题.解决该问题的思路是，首先检测出高层建筑对低矮目标的遮蔽情况，在正射纠正过程中这些遮蔽地区不分配影像信息，最后利用影像序列的重叠信息或利用仿真方法修复遮蔽地区.

已有遮蔽检测方法大多基于计算机图形学的ZBuffer消隐技术，比较目标距光源点的距离判断其遮蔽情况［1－4］.这类方法在处理大量复杂的矢量建筑物模型时效率很低［5］；而在处理栅格数据时，因难以捕获相机畸变和投影成像造成的栅格变形，将导致错误的判断竞争关系，造成大量伪遮蔽和伪可见噪声［6－7］.此外，该方法无法处理被高窄墙体遮蔽目标的伪可见问题，即所谓“M-Portion”问题［3，8－9］，改进的方法需要将墙面纳入遮蔽检测计算过程，大大增加了计算工作量和复杂度.产生这些问题的根本原因在于复杂多样的现实世界难以用简单的三维体块模型表示，基于图形学的Z-Buffer技术也未顾及摄影测量的特殊情况［10］.

1 快速遮蔽检测方法

从计算效率的角度出发，优先考虑采用栅格DSM参与计算；从计算稳健性的角度出发，需克服各类伪可见和伪遮蔽噪声.为此可采用径向投影高度角的变化情况进行遮蔽区域分析，其原理是，同一方位射线上逐渐远离天底点时目标的投影高度角逐渐增大；在进入被遮蔽地区时投影高度角突然变小，直到离开被遮蔽地区时又回到之前的数值.据此，可以容易发现被遮蔽地区.该方法对遮蔽地区和遮蔽源的扫描都是连续的，对遮蔽地区的判断仅需高层建筑的边缘栅格，可避免因栅格变形造成的错误竞争，也可避免因缺乏墙面信息造成的“M-Portion”问题［11］.

为全面检测遮蔽情况，从天底点到DSM边缘所有栅格的射线都需要分析太阳高度角的变化情况.此时，邻近天底点的栅格将被多次访问，导致效率低下.可行的方法包括径向扫描和螺旋扫描两种方法（图1），前者首先以不同半径将DSM划分为多个圆环，在每个圆环内根据完整扫描的要求制定合适的扫描和方向间隔；后者从天底点出发以螺旋方式扫描DSM栅格，最终完成对DSM的完整扫描.本文力图通过分析两种方法的计算步骤和算法复杂度，推导最优的配置方案，最终提出具有最高效率的遮蔽检测方法.

两种方法基本原理一致，区别在于对DSM扫描的方式，以及由此导致的不同的扫描面积、方位射线数和进行可见性判断的面积.因此，通过分析计算过程的复杂度，找到最优的分配方法，可找到复杂度最低的遮蔽检测方法.首先，基于算法原理利用严密的几何关系可以推出两种方法不同的扫描面积、射线数据和可见性判断面积；然后，利用计算数学理论推导出各个步骤包含的计算复杂度；最后利用微积分的极值理论求出具有最小复杂度的方法.

1.1 算法过程

首先，对于给定DSM 栅格（x，y，z），投影中心坐标为（Xc，Yc，Zc），计算并记录每个DSM 栅格的投影高度角a.

图1 径向扫描和螺旋扫描方法Fig.1 Illustration of radial sweep and spiral sweep method

然后扫描DSM，分析沿扫描射线的投影高度角变化情况.对于径向扫描，计算当前范围内由边缘栅格决定的方位射线，沿着射线方向逐个比较a变化，同时记录栅格的可见性，若发现a突然变小，则记录当前的amax值，投影高度角未超出amax的栅格都被认为是被遮蔽的地区，其可见性值设为零.

对于螺旋扫描方式，需要计算每个栅格与天底点确定的方位射线，然后比较其与射线上相邻且靠近天底点的栅格的a变化，同样记录当前的amax值，投影高度角未超出amax的栅格被认为是被遮蔽的地区，其可见性值设为零.

由上述分析可知，两种方法效率的影响因子包括不同的扫描面积、方位射线数和可见性判断面积，需要研究这些因子的优化方法，以提高算法效率.

1.2 计算复杂度分析

假定DSM为圆形且未进行分割，设定栅格尺寸为1，径向扫描的每条扫描射线的面积为

式中：r为DSM半径；d为DSM栅格的直径.而扫描射线的数量可以由下式获得：

可计算得到径向扫描的扫描面积为

从式（4）可以看到，径向扫描面积两倍于DSM实际面积，存在较大程度的重复扫描.考虑到若一个栅格的可见性已经判定，那么在重复扫描的过程中，不需要再重复判断，因此需要进行可见性判定的面积为πr2.

若DSM进行了n次分割，根据图1，其扫描面积、射线数量和可见性判断面积计算公式演变为

式中：Sr为径向扫描面积；Lr为径向扫描射线数；Vr为径向扫描需要进行可见性判断的面积.

对于螺旋扫描，容易判断其扫描面积和可见性面积等于DSM的面积.考虑到每个栅格都需要计算方位射线，可以得到

式中：Ss为螺旋扫描面积；Ls为螺旋扫描射线数；Vs为螺旋扫描需要进行可见性判断的面积.

通过计算数学分析，可知计算每一条射线参数需要5个基本操作（加减乘除），确定栅格的可见性需要4步操作.在径向扫描中需要3步操作来扫描一个栅格，而螺旋扫描中仅需要1步可实现栅格定位.考虑到这些基本操作的复杂度是常量，两种方法的计算复杂度可以估计为

式中：Cr代表径向扫描方法的复杂度；Cs代表螺旋扫描方法的复杂度.

式（7）表明径向扫描方法的复杂度主要依赖于DSM划分方式，而螺旋扫描的复杂度依赖于DSM的大小.在给定DSM的情况下，螺旋扫描的复杂度是一个常量，而径向扫描的复杂度是可以通过修改DSM划分方法而改善的.

1.3 求解最小复杂度

两种方法计算复杂度的差异可以通过下式计算：

若采用等差数列方法划分DSM，即ri－ri－1＝m（m＞0）可得

在给定DSM的情况下，r为常量，此时求最小复杂度演变为求解单变量函数的最小值问题.考虑到n和ΔC在［1，r］是连续的，且ΔC是可微的，对ΔC求极值可得到

则可获得最小ΔC的分割数为

带入式（9），计算得到最小的ΔC为

当r大于10时，ΔC的值小于0，也就是说通过选择最佳分割数，径向扫描的复杂度此时小于螺旋扫描.这证明螺旋扫描效率高于径向扫描的一般性结论不成立.

若采用等比数列方法分割 DSM，即ri／ri－1＝m（m＞0），式（8）演变为

由于求解幂指数的极值比较复杂，反而会影响计算复杂度，本文不再讨论.同样的，也可以采取其他方法划分DSM，进而比较算法复杂度.

上述推导已经证明通过求解最优分割数，径向扫描方法的复杂度可以得到大幅度改善并获得比螺旋扫描更高的计算效率.因此，基于最优分割数的径向扫描方法可认为是一种高效稳健的遮蔽检测方法.

2 实验与讨论

考虑到LiDAR（激光扫描数据）可以快速有效地获取地表三维模型，已经日益成为城市三维建模的重要数据源，实验采用LiDAR点云构建高分辨率DSM；作为基础实验数据，引入Z-Buffer方法对比算法的稳健性和效率.

实验数据如图2所示，包含加拿大多伦多地区的LiDAR点云数据（图2a）和航空影像数据（图2b）.航空影像的分辨率为0.5m，LiDAR点云的分辨率为1.0m.航空影像记录了城市目标的灰度和纹理信息，LiDAR数据记录城市目标的三维坐标.在图2b中左下部分有重大偏移和遮蔽的地区，从图2c展示的局部点云可知，这些地区点云中不存在变形和遮蔽.因此，利用LiDAR构建的DSM实施航空影像的真正射纠正是可行的.

图2 实验数据Fig.2 Experimental data

图3展示了利用典型方法实施遮蔽检测的实验结果.图3a展示了Z-buffer方法产生的可见性矩阵，图3b展示其局部放大图，图3c展示了对其进行Z-buffer遮蔽检测的结果，图3d展示了径向扫描方法产生的可见性矩阵.从图3a～3c中，可以发现大量的离散噪声，以及由噪声组成的黑色曲线.这些是因为Z-buffer方法无法适应DSM起伏和分辨率差异引起的.如果采用滤波方法消除噪声，将降低影像的细节层次.而在图3d～3f中则没有这些离散噪声和曲线，表明基于径向投影高度角变化的检测方法可以有效地适应城市地区复杂的DSM起伏和分辨率差异.图3也表明，由于LiDAR数据分辨率低于遥感影像，由其构建的DSM生产的真正射影像中缺乏精确的断裂线，建筑物轮廓出现锯齿.

图3 遮蔽检测实验Fig.3 Detecting occlusions with typical methods

高窄的结构容易产生M－portion问题，在数据源中较难发现此类结构，为了验证算法对于 M－portion问题的处理能力，图4展示了两种方法对于高窄结构的检测结果.图4a展示了实验中人为添加了模拟的高窄结构.图4b是利用Z-buffer方法检测的结果.可见建筑顶面的遮蔽范围被检测出来，而部分被墙面遮蔽地区未被检测到.这是因为考虑到效率问题，一般未将墙面纳入DSM参与计算和分析，若墙面遮蔽地区大于建筑本身面积，则会导致伪可见的情况，即所谓“M-portion”问题.高且窄的建筑物最容易导致这类问题的出现.图4c是利用径向扫描方法检测的结果，可见建筑物顶面和墙面遮蔽的地区都被检测出来了.根据该方法的原理，仅需扫描建筑物边缘的栅格即可完成对顶面和墙面遮蔽的检测，无需额外改进和其他步骤.而绝大部分基于Z－buffer原理的方法及改进算法往往很难高效地发现并解决M－portion问题.从这一点来说，该方法极大地提高了遮蔽检测的稳健性.

图4 对高窄结构的检测Fig.4 Detecting occlusions with simulated DSM

本文对典型方法的计算效率做了对比实验，包括Z-Buffer方法，径向扫描方法、螺旋扫描方法和最优分割数径向扫描方法.实验设备是具有双核2.8G CPU和3.0GDDR内存的微机.图5展示了效率实验的结果.

图5 效率实验Fig.5 Performance of typical methods

从图5可见，随着给定DSM尺寸增大，各类方法消耗的时间增加.其中，对于任何尺寸，Z-Buffer总是消耗最长时间，径向扫描方法次之，螺旋扫描方法又次之，最优分割径向扫描方法消耗的时间最短.这说明最优分割径向扫描方法在这些方法中具有最高效率.考虑到其他基于Z-Buffer原理的方法或改进算法往往比基于栅格数据的Z-Buffer消耗更长的时间［11］，可认为最优分割径向扫描方法是具有最高效率的遮蔽检测方法.此外，值得注意的是，最优分割径向扫描方法的效率曲线表现为线性函数，而ZBuffer方法的效率曲线表现为高次函数.随着给定DSM尺寸逐渐增大，两种方法消耗的时间差距会越来越大.

3 结语

近十几年来，国内外学者努力尝试寻找一种高效、稳健、精确的遮蔽检测方法.已有的基于ZBuffer技术的遮蔽检测方法并不适应摄影测量的特殊要求，已经被证明存在诸多难以克服的问题，而一些改进算法往往过于复杂，计算不稳定且效率低下.实验表明，基于径向投影高度角变化的方法可以有效解决因地形起伏、分辨率差异和高窄墙体引起的各类伪可见和伪遮挡问题，实现稳健、精确的遮蔽检测.对比实验表明，同样数据量下该方法消耗时间最少，且数据量增大时的时间增量也最少.本文的理论推导和实验表明，最优分割的径向扫描方法可以实现稳健、精确、快速、高效的遮蔽检测，有利于大规模高质量真正射影像制作.

［1］ Amhar F，Jansa J，Ries C.The generation of true orthophoto using a3D building model in conjunction with a conventional DTM ［J］.International Journal of Photogrammetry and Remote Sensing，1998，32（4）：16.

［2］ Skarlatos D.Orthophoto production in urban areas［J］.Photogrammetric Record，1999，16（94）：643.

［3］ Rau J Y，Chen N，Chen L.True orthophoto generation of built-up areas using multi-view photos［J］.Photogrammetric Engineering and Remote Sensing，2002，68（6）：581.

［4］ Bang K，Habib A F，Kim K，et al.Comprehensive analysis of alternative methodologies for true ortho-photo generation from high resolution satellite and aerial imagery［C］／／The 5th International Symposium on Mobile Mapping Technology MMT’07.Padua：International Society for Photogrammetry and Remote Sensing，2007：1-5.

［5］ 钟成，黄先锋，李德仁，等.真正射影像生成的多边形反演成像遮蔽检测算法［J］.测绘学报，2010，39（1）：59.

ZHONG Cheng，HUANG Xianfeng，LI Deren，et al.Polygon based inversion imaging for occlusion detection in true orthophoto generation［J］.Acta Geodaetica et Cartographica Sinica，2010，39（1）：59.

［6］ Braun J.Aspects on true orthophoto production［C］／／Proceedings of 49th Photogrammetric Week.Heidelberg：［s.n.］，2003：205-214.

［7］ Sheng Y W.Minimising algorithm-induced artefaccts in orthophoto generation：a direct method implemented in the vector domain［J］.The Photogrammetric Record，2007，22（118）：151.

［8］ Rau J Y，Chen N Y，Chen L C.Hidden compensation and shadow enhancement for true orthophoto generation［C／CD］∥Proceedings the of 21st Asian Conference on Remote Sensing.Taipei：［s.n.］，2000.

［9］ Sheng Y W，Gong P，Biging G.True orthoimage production for forested areas from large scale arial photograph［J］.Photogrammetric Engineering and Remote Sensing，2003，69（3）：259.

［10］ Zhong C，Li H，Li Z，et al.A vector-based backward projection method for robust detection of occlusions when generating true ortho photos［J］.GIScience and Remote Sensing，2010，47（3）：412.

［11］ Habib A F，Kim E，Kim C.New methodologies for true orthophoto generation［J］.Photogrammetric Engineering and Remote Sensing，2007，73（1）：699.