城市环境是由建筑物、道路及其附属设施、民众等共同构成的一个复杂的动态系统。基础设施的科学安置对城市的高效运作至关重要，故科学高效的城市空间分析必不可少，而可视域计算是空间分析的一个重要组成部分^[1]。道路是城市中人、车流量最大的场景，其主要由行人、机动车、道路地物要素(交通标志牌、交通信号灯、行道树、路灯、道路等)构成，道路系统中的任何一个要素出现异常都可能导致交通事故。交通事故原因主要分为驾驶员违规驾驶的主观原因和道路基础设施安置不当的客观原因^[2]，由于基础设施安置不当导致事故时，成因之一则是交通标志牌被遮挡造成道路信息获取不足。当机动车位于交通标志的通视盲区时，驾驶员或视觉传感器无法获取完整的交通标志信息，因而不能及时作出正确合理的决策，有可能导致交通事故的发生^[3-7]。还有部分交通事故是由夜晚道路照明不良所导致，夜间树木对路灯光线的遮挡导致道路照明状况不佳，严重影响夜间行车安全^[8]。综上所述，对道路环境中相关地物的可视域计算与分析至关重要。

1 空间可视域分析

通视分析首先出现在军事领域中，随着经济与技术的发展，其在无线电通讯及空间规划领域也得到了广泛应用。由于其原理简单、应用范围广，几乎所有的GIS软件都有可视域计算功能，但对于大范围、遮挡严重的城区环境，可视域计算仍然面临很大挑战。

GIS可视域计算基本分为3类：(1)两点之间的通视计算；(2)可视域计算；(3)视点能见度计算。总体来说，传统的可视域计算方法首先对三维模型如数字高程模型(digital elevation model，DEM)、数字地形模型(digital terrain model，DTM)等插值，构建满足精度要求的格网；然后在模型中添加建筑物、植被等图层；最后计算给定视点的可视域^[9]。然而，目前三维模型构建的自动化水平、精度和完整度等很难满足城区环境可视域分析的要求。车载激光扫描系统可以高速度、高密度、高精度地获取道路及两侧地物的三维坐标和属性信息(如反射强度、回波波形等)，在无需建模的情况下直接对地物目标进行可视域分析，为大规模城市场景可视域计算与分析提供了一种全新的技术手段^[10]。目前，基于激光点云的城市场景可视域计算方法主要包括基于面元^[11-12]、基于体素与基于HPR(hidden point removal)的方法^[13-14]。基于DEM面元的可视域计算一般需将已经生成的DEM规则格网化，而直接对格网遍历计算，效率非常低下。为加速计算，文献[11]通过构建金字塔四叉树空间索引实现格网的快速查找索引，在提高查找、计算效率的同时，也带来了数据冗余的问题；同时，可视域计算结果的精度会受到格网分辨率的影响。文献[12]基于三角网面片计算可视域，三角网对地物表面的描述精度高于基于DEM插值形成的格网，其通过切片算法计算地形可视域。基于面元的可视域计算一般需要根据点云构建目标表面的三角网或规则DEM面片，而构建组织三角网、DEM是一项非常耗时的工作。基于HPR的可视域计算首先将视场范围内的点云根据到视点距离的反比关系映射到反转空间，并在映射空间中构建凸包，落在凸包表面的点即为原始可视点^[13]。基于HPR的方法不需要构建地物表面三角网，可以直接计算出视点视场内的可视域，但需要求解凸包。道路环境点云数据量巨大，求解凸包非常耗时，而且只能计算出点云中的直接通视点，不适合应用于通视空间的计算。

基于道路环境通视空间的计算需求以及车载激光点云的具体特点，本文采用基于体素的深度缓存方法计算视点在视场内的通视区域(通视空间)，避免了模型构建以及凸包计算的过程；同时，本文通过对道路场景交通标志牌的遮挡分析、路灯有效照明区域的计算，验证了方法的有效性。

2 可视域计算方法

道路环境中地物目标可视域计算的主要流程包括视场快速估计和可视域估计分析等两个步骤，方法的整体流程见图 1。

2.1 视场快速估计

视场估计是可视域计算的关键步骤，包括视场基本空间参数估计、视场精化估计。空间参数包括视点位置、视场范围。视点一般为研究对象所处的空间位置，视场一般根据特定应用场景中的视角、视距确定。

由于点云数据海量，需要在可视域估计之前确定有效视场范围，从而降低可视域估计的计算复杂度。本文采用由粗到精的动态视场确定策略，具体方法如下：在可视域估计中，噪声点的存在可能会导致局部可视域计算结果不准确，故对测量过程中的各种系统误差和随机误差导致的噪声点，本文通过PCL(point cloud library)^[15]中的统计离群点移除算法将其剔除。

由于点云具有数据量大、空间离散性、分布不均匀等特点，为了提高点云数据组织与管理效率，需要对点云构建空间索引^[16]。本文分如下步骤构建空间索引：首先基于$xoy $平面构建长、宽尺寸固定的规则体素为根结点，根结点的长度为$ {g}_{l}$，宽度为$ {g}_{w} $，同时计算根结点的相关几何属性，包括结点中点的数量、点号索引、最低高程值、最高高程值、根结点在$ xoy $平面的投影坐标等。若根结点在$ xoy$平面的投影坐标落入视场在$xoy $平面的投影区域，则被选为感兴趣根结点，如图 2(a)所示；然后对于感兴趣根结点，内部划分长、宽、高尺寸分别为$ {v}_{x}$、$ {v}_{y} $、$ {v}_{z} $的叶子结点体素，并统计叶子结点空间体素的相关属性，包括空间体素中点的数量、点号索引、体素几何中心。此种空间索引结构的深度为2层，而且每个根结点的叶子结点体素层数由根结点中最低高程与最高高程的差值自适应确定，故叶子结点的空结点数量较少，可以加快后续计算中的空间查询效率，详见图 2(b)；最后根据视场角度、距离关系进一步筛选感兴趣根结点下的叶子结点，快速准确确定参与可视域估计的有效视场点集。经过以上过程，位于视场内的叶子结点被精确筛选，减少了后续可视域计算的数据量，有助于提高可视域计算的效率。

2.2 可视域计算

根据§2.1确定的视点及视场范围内的叶子结点，本文通过基于体素的深度缓存算法计算点云中视点的可视域。将视场中的叶子结点由原始坐标转换到以视点为坐标原点的相机坐标系$ c-UVN $下，相机坐标系的定义如下：以视点为相机坐标系的原点$ {O}_{c}({x}_{c}, {y}_{c}, {z}_{c}) $；观察者位于相机后方，右手方向为$U $轴的方向；在相机平面上，按照右手法则，大拇指指向观察者，垂直于$U $轴的方向为坐标轴$ V$；按照右手法则，由$ U $轴旋转至$ V $轴的大拇指指向为$ N$轴。两坐标系的转换过程如下：

$ \left[\begin{array}{c}{X}_{c}\\ {Y}_{c}\\ {Z}_{c}\\ 1\end{array}\right]={\mathit{\boldsymbol{C}}}^{-1}·\left[\begin{array}{c}{X}_{u}\\ {Y}_{u}\\ {Z}_{u}\\ 1\end{array}\right] $

(1)

式中，$ ({X}_{c}, {Y}_{c}, {Z}_{c})$为叶子结点几何中心点$ P $在相机坐标系中的坐标；$ ({X}_{u}, {Y}_{u}, {Z}_{u}) $为其在原始世界坐标系中的坐标；$ {\mathit{\boldsymbol{C}}}^{-1}$为坐标转换参数矩阵的逆矩阵。$ {\mathit{\boldsymbol{C}}}^{-1} $由以下过程计算得到：已知世界坐标系中的视点坐标为(${x}_{c}, \mathrm{ }{y}_{c}, \mathrm{ }{z}_{c}) $，垂直于相机平面并指向目标的方向记为$ \mathit{\boldsymbol{r}} $，相机辅助竖直方向记为$ \mathit{\boldsymbol{u}} $，则有：

$ \mathit{\boldsymbol{N}}=(-1·\mathit{\boldsymbol{r}}).n\left(\right) $

(2)

$ \mathit{\boldsymbol{U}}=(\mathit{\boldsymbol{u}}\times \mathit{\boldsymbol{N}}).n\left(\right) $

(3)

$ \mathit{\boldsymbol{V}}=(\mathit{\boldsymbol{N}}\times \mathit{\boldsymbol{U}}).n\left(\right) $

(4)

$ \mathit{\boldsymbol{R}}=\left[\begin{array}{cc}\begin{array}{cc}{U}_{x}& {V}_{x}\\ {U}_{y}& {V}_{y}\end{array}& \begin{array}{cc}{N}_{x}& 0\\ {N}_{y}& 0\end{array}\\ \begin{array}{cc}{U}_{z}& {V}_{z}\\ 0& 0\end{array}& \begin{array}{cc}{N}_{z}& 0\\ 0& 1\end{array}\end{array}\right] $

(5)

$ {\mathit{\boldsymbol{T}}}^{-1}=\left[\begin{array}{cc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& -{x}_{c}\\ 0& -{y}_{c}\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& -{z}_{c}\\ 0& 1\end{array}\end{array}\right] $

(6)

$ {\mathit{\boldsymbol{C}}}^{-1}={\mathit{\boldsymbol{R}}}^{\mathrm{T}}·{\mathit{\boldsymbol{T}}}^{-1} $

(7)

式中，$ .n\left(\right)$为对向量的归一化计算；$ \mathit{\boldsymbol{R}} $为坐标变换旋转齐次矩阵；$ \mathit{\boldsymbol{T}} $为平移齐次矩阵；N、U、V分别为相机坐标系三轴方向的单位向量；[U_x U_y U_z]、[V_x V_y V_z]、[N_x N_y N_z]分别为U、V、N 3个轴在原世界坐标系下的坐标分量。

经坐标转换后，视场范围内的叶子结点体素均被转换到以视点为原点的相机坐标系下，将视场中的所有叶子结点体素根据透视投影变换到距离视点距离为$ f$、栅格尺寸为$ {d}_{x}\times {d}_{y}$的像平面上，其中$ f$为视场中叶子结点到视点的最远距离，$ {d}_{x} $与$ {d}_{y} $为深度缓存平面的分辨率。

设$P $为相机坐标系中某一叶子结点的体素中心，坐标为$ ({U}_{P}, {V}_{P}, {N}_{P}) $，投影到平行于相机平面的像平面$ o-xy $的坐标为$ ({x}_{p}, {y}_{p}) $，则有：

$ \left[\begin{array}{c}{x}_{p}\\ {y}_{p}\end{array}\right]=\frac{f}{\left|{N}_{P}\right|}\bullet \left[\begin{array}{c}{U}_{P}\\ {V}_{P}\end{array}\right] $

(8)

待所有叶子结点体素均投影到像平面上后，统计落入深度缓存平面每个栅格中距离视点距离最近的非空叶子结点体素，将其作为直接可视体素，所有直接可视体素中的点的集合即为可视地物表面。如图 3所示，$ {P}_{1}$、$ {P}_{2} $、$ {P}_{3} $、$ {P}_{4} $为叶子结点体素的中心点，经透视投影到深度缓存平面的同一栅格中，其中$ {P}_{1} $距离视点距离最近，则$ {P}_{1} $所代表的叶子结点体素为可视体素，$ {P}_{1} $所代表的叶子结点体素中的点为可视点。

可视域计算根据输出分为通视空间与通视面的计算，若要求输出直接通视地物，则输出所有直接可视叶子结点体素中的点集；若要求输出通视空间，则输出像平面各栅格中到视点距离小于直接可视体素到视点距离的空叶子结点体素。

3 地物可视域分析实验

实验中数据处理参数取值见表 1。根结点的尺寸取值参考文献[17]；当叶子结点尺寸较小时，叶子结点数量较多，导致计算效率低，当叶子结点尺寸取值较大时，导致计算结果不满足可视域精度需求，实验发现叶子结点体素尺寸为$ 0.25\mathrm{ }\mathrm{m}\times 0.25\mathrm{ }\mathrm{m}\times 0.25\mathrm{ }\mathrm{m} $时，算法可以在计算效率和精度之间取得很好的平衡；深度缓存平面格网分辨率根据叶子结点体素尺寸确定，取相邻叶子结点体素的最小距离为缓存平面格网分辨率。

3.1 交通标志牌可视域分析

实验所用数据为由Optech Lynx SG1移动激光扫描系统获取的武汉大学部分路段点云数据(见图 4)。该设备有2个激光扫描仪，最大观测距离为250 m，测距精度可达5 mm，绝对精度可达5 cm。实验中视点位于距离地面高度1.5 m处(驾驶员眼睛距离地面的高度一般在1.0~2.5 m之间)，为得出此视点下对交通标志牌的通视性，视场为视点与交通标志牌的最小包围圆锥体。

实验先根据文献[17]中的方法提取出交通标志牌与道路，然后将道路栅格化，道路格网中心对于标志牌的可见度代表当前道路格网覆盖路面对目标交通标志牌的可视度。图 5(a)为视点A对于该交通标志牌的通视状况，图 5(b)为视点B对于该交通标志牌的通视状况。根据本文提出的可视域计算方法，对标志牌的可视度定义为视点对标志牌的直接可见面积占标志牌总面积的比例。图 5(c)为采用本文方法在道路上任意视点对该交通标志牌的可视度输出结果，图 5(d)为该道路对于该交通标志牌的可视度渲染输出结果的部分道路放大图。由此可见，本文方法可直接在点的基础上计算出交通标志牌的通视空间而不需要构建模型，可以及时发现道路旁侧的交通标志牌通视盲区，对安全驾驶有重要意义。

3.2 道路路面照明有效区域分析

实验所用数据为车载激光扫描系统获取的京承高速路段数据。根据参考文献[17]对非地面点分类并提取路灯，在地面点中提取道路。通过分类结果计算得到相关参数为：路灯灯高13 m，同侧相邻两灯距离为36.2 m，路灯仰角为15°，有效路宽为15.0 m。根据道路照明发光二极管(light emitting diode，LED)的光源特点及经过二次配光后的路面光照规律建立灯光在路面上的分布模型。目前，为满足道路照明要求及节约能源，道路照明一般将多个LED运用透镜再次配光，使其在路面上的照度呈沿路面分布的近似矩形光场。配光后单个路灯的有效照明区域沿道路延伸方向分布在±60°之间，垂直于道路延伸方向分布在±30°之间^[18]。经过再次配光后，单个路灯的路面照度根据LED路灯的光强配光曲线及路灯地面照度的分布规律得到路面照度的关系^[19]，从而推导出路面某点照度$ E(x, y, z)$与灯下点照度$ {I}_{o} $的近似计算公式：

$ E(x, y, z)=\frac{{I}_{o}s\left[{t}_{x}\right(x-{x}^{\text{'}}{)}^{2}+{t}_{y}(y-{y}^{\text{'}}{)}^{2}+1]}{\left[\right(x-{x}^{\text{'}}{)}^{2}+(y-{y}^{\text{'}}{)}^{2}+(h+z-{z}^{\text{'}}{)}^{2}{]}^{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}} $

(9)

式中，$s $为位于路面的灯下点与光源中心的连线距离；$ {I}_{o}$为灯下点$ {O}^{\text{'}}$处的光强；$ {t}_{x} $、$ {t}_{y} $为模拟配光参数；$ (x, \mathrm{ }y, \mathrm{ }z)$为路面点的坐标；$ ({x}^{\text{'}}, {y}^{\text{'}}, {z}^{\text{'}}) $为灯下点坐标；$ h$为路灯灯高。路灯安置简图如图 6所示。

通过对道路照明路面照度分布的模拟，最终得出路面照度分布结果，如图 7所示。图 7(a)为单灯照明路面照度渲染输出结果，单灯的有效照明区域为灯下矩形区域^[18]，即图 7(a)中矩形红框内区域；图 7(b)为同侧相邻两路灯的路面照度渲染输出结果，两灯之间路面上的照度为相邻两灯照度叠加的结果，有效照明区域为图 7(b)中红色矩形框内区域。表 2为矩形长、宽及矩形内的照度参数。根据表 2，假设单灯有效区域内最大照度为26.67 lx，则有最小照度为17.68 lx，平均照度为20.81 lx；在同侧相邻路灯的照明有效区域内，最大照度为27.51 lx，最小照度为17.68 lx，平均照度为21.57 lx。以上两个模型均满足《城市道路照明设计标准CJJ45—2006》中一级道路平均照度为20~30 lx的要求。当前实验数据中道路旁侧行道树及广告牌对路灯灯光均无遮挡，但对于一些行道树枝叶繁茂或者广告牌布设不合理的路段，由于树木枝叶及广告牌遮挡路灯光线，导致最终照到路面上的照度分布不满足道路照明要求，常见的情况是路面光照明暗相间、分布不均，车辆在如此照明条件下行驶时无法获得正确的路况信息，有可能发生车祸。后续可以根据具体研究路段中路灯的布设方式，建立路面照度模型并计算输出路面照度分布情况，根据路面照度分布情况判断路面照度是否符合路面照明标准，以便及时对照明不良路段开展故障排查工作，保障民众人身财产安全。

4 结语

本文提出了基于深度缓存的道路场景可视域快速计算方法，并将该方法用于分析城市环境中道路、驾驶员、交通标志牌、路灯、树木等之间的空间关系。通过估计道路两旁的交通标志牌的有效作用空间，分析交通标志牌的通视盲区，对于道路安全养护、无人驾驶环境的安全分析有决策支持作用。通过构建道路环境路灯的照明模型，分析路面照度分布，快速探测道路照明不良路段，对夜晚安全驾驶具有重要意义。总而言之，本文提出的方法可以直接快速地对道路环境中对象之间的通视性、可见度等进行计算分析，为基础设施科学安置、基础设施运行健康状况监测以及城市形态分析、城市规划等提供科学的辅助决策。

未来的工作包括根据该模型对城市建筑物窗口的直接可见绿地面积进行计算等^[20]。除此之外，为使可视域、可视度计算结果更加准确，未来的工作还包括路面动态目标探测及滤除，以消除动态目标对计算结果的影响。