城市地下管网是发挥城市功能和确保城市快速协调发展的重要基础设施。随着我国市政基础设施信息化建设的力度加大，依赖于施工设计图以及计算机辅助设计（computer aided design, CAD）图纸的传统城市地下管网已经不能满足设备监控查询的需要，因此，地下管线的逆向重建数据成为城镇现代化建设中不可或缺的基础资料，也成为城市决策的重要基础资源之一^[1]。由于城市地下管网中工井内部环境复杂，各类管线纵横交错，使用人工测量方法难以获取全部的井内信息；且井内存在积水淤泥等情况，人工测量存在一定的危险性。而激光雷达技术具有快速、非接触性、高精度、数字化、自动化等优点^[2]，用其获取基础数据可提高工井普查的安全性，获得更详尽的基础数据。利用激光雷达技术辅助地下电缆工井的逆向重建已成为目前的主流方法。

地下管网主要由管线与工井组成。在管线连接方面，卢丹丹等^[3]提出一套高精度自动三维建模的方法与思路，利用二维管线普查数据，按照管点和管线段的特点，生成地下管线三维模型，并可同步更新其拓扑关系；毕天平等^[4]以二维管线数据为基础，建立三维管网模型数据标准，进行三维管网整体自动建模及可视化；李清泉等^[5]将Polyvrt拓扑结构和空间实体模型结合对管网数据进行管理，生成管线模型并记录其拓扑关系。但以上研究中，三维空间数据组织一般将管网抽象为管点和管线段两种实体，采用两点一线的静态层次结构，空间拓扑关系也仅是点与线、线与线的关系，缺少管线与对应墙壁间的线面关系以及管点与管井间的点面关系。

电缆通过管孔与管线相互连接，管线交错关系复杂，工井可以清晰地反映出地下电缆的拓扑连接关系，对建立地下电缆工井网络模型有十分重要的意义。在工井建模方面，周京春^[6]将工井抽象化为井室主体以及井盖对象，进行CSG（constructive solid geometry）实体建模；刘军等^[7]基于Skyline提出建立城市地下管线三维系统的方法，将不同的工井抽象为不同的三维符号，采用3Ds Max建模后导入管线系统；王舒等^[1]将工井抽象为点模型，赋予其平面位置、高程、构筑物等信息；吴思等^[8]基于DirectX自动化建模，将地下管线井室分类并构建模型。钟远根等^[9]在传统管线数据标准基础上扩充了井室调查内容，以满足井室自动建模需要。以上建模方法在不同的场景具有不同优势，但是只采用整体方式建模，拓扑结构仅抽象为点线连接，缺少以面及管孔表达工井模型拓扑结构的方式。在检测、维护以及新增管线时可通过对应管孔找寻到对应管线。综合以上因素，对工井进行以面为单位的精细化建模及管孔表达，对地下管网施工、生产及维护都具有重要意义。

基于点云数据进行建筑建模方面，许多学者也进行了研究。Edelsbrunner ^[10]提出一种基于隐式曲面的算法Wrapping，在基于尺度空间和多尺度分析时重建效果较好，但大型曲面重建时效率较低^[11]。Vosselman^[12]利用三维霍夫变换（3D Hough transform）得到建筑物平面片，从而实现建筑物三维重建。贺彪等^[13]提出一种自动拓扑重建算法，可以有效重建多源异构建筑物，在建筑物内拓扑以及外拓扑方面重建效果良好。以上算法在模型曲面重建上都有非常重要的参考价值，但是由于扫描数据量大，直接进行曲面三角网重建效率非常低；同时，由于扫描设备以及扫描条件限制，在扫描物体特定区域或位置时有不同程度的数据缺失以及噪声，常规算法难以重建出曲面模型，并且无法对墙面孔洞进行建模。

综上，本文提出一种针对具有孔洞的地下工井三维模型拓扑重构算法。该算法利用点云数据拟合出工井井室的边线，进行以面为单位的参数化建模，再针对每个面的孔洞情况对其进行拓扑重构，完成地下工井精细化三维重建，解决了直接采用点云进行构网时模型不完整及速度较慢的问题。

1 具有孔洞的地下工井模型拓扑重构算法 1.1 算法流程

本文设计的具有孔洞的地下工井模型拓扑重构算法流程如图 1所示，大致分为4步：①原始点云预处理，提取面点云并提取边界点；②通过边界点以面为单位构建无孔工井模型；③提取点云中孔洞中心点并对中心点进行投影；④在投影点处对墙面进行拓扑重构。

1.2 算法实现 1.2.1 点云数据预处理与点云边界点提取

直接获取的激光点云数据含有大量冗余信息，难以直接应用于实际重建工程，需要先对点云数据进行精简去噪等预处理。本文采取一种可移动网格划分的点云精简过滤算法，该算法主要通过空间网格对点云进行精简^[14-15]。

在提取点云边界点时，由于工井模型通常形状较为规则，通过采用Schnabel等^[16-18]提出的高效率的随机采样一致性算法（random sample consensus, RANSAC）能够有效准确地检测出平面、圆柱、规则型曲面等基础形状，可以提取出工井点云数据每个面的点云片。对单个点云面片使用主成分分析法（principal component analysis, PCA）获得点云模型的法向量；依据法向量可以将其旋转至与XOY面平行，并投影至该平面获得投影轮廓；采用扫描线方式^[19]提取轮廓的边线初始点集；对点集进行拟合，获得模型边界线关键点。利用边界关键点以及获取的点云参数来进行参数化建模。图 2展示了顶部点云视角及拟合后的边界关键点^[20]。

1.2.2 构建工井墙壁面模型

工井墙壁面模型中每个面均可以用参数方程表示。通过获得的顶面关键点数据，采用高度平移可以获得点云模型底面的关键点。单面墙壁可以通过提取出的上顶面及下底面关键点建立参数方程表示墙面信息，同时采用三角网构建墙面实体显示墙面。

在20组平均每组大约500万点云的实验数据中，将单面墙壁大致分为平面墙壁以及曲面墙壁两种。图 3展示了四曲面工井中3种不同的墙壁类型：多边形平面、曲面和矩形。

1）平面墙壁

对于平面墙壁，墙壁所在的无限平面参数方程由过平面中的任意一点P以及面法线N确定，并提取关键点。平面墙壁的关键点如图 4所示。将平面关键点以顺时针或者逆时针方向排列生成环类型的元素，然后以环的范围为限制界限对无限平面进行剪裁，生成可以进行拓扑运算的平面参数方程。

2）曲面墙壁

曲面墙壁在实际的点云数据中不一定是完美的部分圆柱面，所以无法规则地参数化建立曲面。本文则根据提取出的曲面墙壁上下边界线的关键点，建立两组等长数组。依据关键点生成B样条曲线，在数学上三维曲面被认为由两条曲线的笛卡尔积生成，所以通过B样条曲线来限定曲面的范围并决定曲面的形状^[21-22]。曲面墙壁关键点如图 5所示。

p次B样条曲线的定义为：

$ C\left( u \right) = \sum\limits_{i = 0}^n {N_{i, p}}\left( u \right){P_i} $

(1)

式中，P_i是控制顶点；N_{i, p}(u)是定义在非周期节点矢量上的第i个p次B样条基函数，其定义为：

$ {N_{i, 0}}\left( u \right) = \left\{ {\begin{array}{*{20}{l}} {1, {\rm{}}{u_i} \le u \le {u_{i + 1}}}\\ {0, 其他} \end{array}} \right. $

(2)

$ \begin{array}{l} {N_{i, p}}\left( u \right) = \frac{{u - {u_i}}}{{{u_{i + p}} - {u_i}}}{N_{i, p - 1}}\left( u \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{u_{i + p + 1}} - u}}{{{u_{i + p + 1}} - {u_{i + 1}}}}{N_{i + 1, p - 1}}\left( u \right) \end{array} $

(3)

根据B样条曲线定义可知，给定控制顶点P_i、曲线次数P及节点向量U（向量中元素为u_i），B样条曲线即可确定。

依据已经获得的B样条曲线，通过笛卡尔积获得B样条曲面。曲面的数据包含u有理标志位r_u，v有理标志位r_v，曲面次数m_u、m_v和带权控制点。u、v的次数m_u、m_v都不能大于25，u，v的控制点个数n_u、n_v都需要大于等于2。

B样条曲面的参数方程如下：

$ \begin{array}{l} S\left( {u, v} \right) = \frac{{\sum\limits_{i = 1}^{{n_u}} \sum\limits_{j = 1}^{{n_v}} {B_{i, j}} \cdot {h_{i, j}} \cdot {N_{i, {m_u} + 1}}\left( u \right) \cdot {M_{j, {m_v} + 1}}\left( v \right)}}{{\sum\limits_{i = 1}^{{n_u}} \sum\limits_{j = 1}^{{n_v}} {h_{i, j}} \cdot {N_{i, {m_u} + 1}}\left( u \right) \cdot {M_{j, {m_v} + 1}}\left( v \right)}}, \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {u, v} \right) \in \left[ {{u_1}, {u_{{k_u}}}\left] \times \right[{v_1}, {v_{{k_v}}}} \right] \end{array} $

(4)

式中，B_{i, j}表示三维点；h_{i, j}表示权重因子。

1.2.3 墙面模型拓扑重构

工井墙壁模型上的管孔是电缆连接的基础，为了满足逆向重建的需求，需要对墙面模型进行布尔运算以完成拓扑重构。

三维模型拓扑重构需要墙面模型参数和投影点处的球模型参数S2。通过RANSAC算法提取得到点云边界线，再通过最小二乘法获得管孔的圆心点^[23]，对拟合出的管孔中心进行投影，获取管孔中心在墙面上的实际坐标，然后通过用户设置的半径或拟合半径生成S2球模型，再对S1模型求差运算，则获得带有孔洞的墙面模型。

模型拓扑重构流程如图 6所示，分为交叉检测与布尔运算两大步骤。首先以输入两个网格模型为基础，对输入的三维模型建立八叉树进行交叉检测，判断相交区域，以减少模型面的遍历次数；计算模型间交点，对交点进行分类；利用计算得到的交点进行网格划分。

两个模型之间的交点检测是布尔运算的基础，快速准确地搜索交叉线是布尔运算的关键。对于两个相交的模型，相交部分的空间不大，通过构建八叉树来划分相交空间，加快交叉检测的速度并降低内存利用率。

拓扑重构算法需要输入运算模型相交部分的三角面片集合A和B。进行布尔运算需要遍历A中的所有三角形边，同时遍历B中的所有三角面片（参数曲面表示）。将A中所有的三角边与B中所有的三角面片求交点，获得由交点分割出的新边。再按照相同步骤计算出B中三角形边在曲面A中的交点以及新连接得到的边线。

算法过程中的计算交点需要将所有的参数曲线统一表示为：

$ C\left( t \right) = \left( {x\left( t \right)} \right)\left( {y\left( t \right)} \right)\left( {z\left( t \right)} \right) $

(5)

所有的参数曲面统一表示为：

$ F\left( {u, v} \right) = \left( {x\left( {u, v} \right), y\left( {u, v} \right), z\left( {u, v} \right)} \right) $

(6)

则参数曲线和参数曲面的交点由如下算法确定：将两空间实体采用参数方程表示，参数方程求交需满足：

$ C\left( t \right) - F\left( {u, v} \right) = 0 $

(7)

令

$ S\left( {t, u, v} \right) = C\left( t \right) - F\left( {u, v} \right) $

(8)

求参数微分方程：

$ {\rm{d}}S = \frac{{{\rm{d}}C}}{{{\rm{d}}t}}\left( t \right) - \frac{{\partial F\left( {u, v} \right)}}{{\partial u}}{\rm{d}}v - \frac{{\partial F\left( {u, v} \right)}}{{\partial v}}{\rm{d}}v $

(9)

将微分方程各参数微分显示，在两边点乘$\frac{{{\rm{d}}C}}{{{\rm{d}}t}} $，令$ M = \frac{{{\rm{d}}C}}{{{\rm{d}}t}}\left( {\frac{{\partial F\left( {u, v} \right)}}{{\partial u}} \times \frac{{\partial F\left( {u, v} \right)}}{{\partial v}}} \right)$，可以获得迭代方程：

$ {t_{i + 1}} = {t_i} - \frac{{\left[ {\frac{{\partial F\left( {u, v} \right)}}{{\partial u}}\left( {\frac{{\partial F\left( {u, v} \right)}}{{\partial v}} \times {\rm{d}}S} \right)} \right]}}{{M\left( {{t_i}, {u_i}, {v_i}} \right)}} $

(10)

$ {u_{i + 1}} = {u_i} - \frac{{\frac{{{\rm{d}}C}}{{{\rm{d}}t}}\left( {\frac{{\partial F\left( {u, v} \right)}}{{\partial u}} \times {\rm{d}}S} \right)}}{{M\left( {{t_i}, {u_i}, {v_i}} \right)}} $

(11)

$ {v_{i + 1}} = {v_i} - \frac{{\frac{{{\rm{d}}C}}{{{\rm{d}}t}}\left( {\frac{{\partial F\left( {u, v} \right)}}{{\partial v}} \times {\rm{d}}S} \right)}}{{M\left( {{t_i}, {u_i}, {v_i}} \right)}} $

(12)

当以上迭代方程满足收敛时，得到最终交点。新连接的交线可将所在面的切割情况等价视为两三角形相交的情况，需进一步讨论。当两个三角面片共面相交时，将问题划归至二维领域探讨。但当两个不共面的三角面片相交时，交点划分会出现3种情况：交于其中一个端点、交点在三角形边上以及交点位于三角形内部。

限于篇幅，本文仅探讨两个不共面三角面片相交时的3种情况。对于获得的交点进行交线连接后，相交部分的三角面片均有了不同程度的切割，需要将这些被切割开的部分三角形进行分类处理，判断其属于三角网A或者B中。如图 7所示，边ab与△def有一个交点n，顶点a(x_a，y_a，z_a)相对于面F_def(x_a，y_a，z_a) < 0。顶点a相对于包含△def的网格位于内部。采用相同的方法，顶点d在内部，c和e在外部。当然也有特殊情况，即两个三角形只有一个交点v(x，y，z)，并且对于顶点v(x，y，z)，F(x，y，z)=0并认为在网格上，对所有的交点进行位置属性判断，赋予外、内、上等属性，便可以判断被切割开的部分三角形所属的网格情况。

以上方法可以实现网格交点求解、交线连接、部分三角形判断属性这几个布尔运算中的重要部分。分割后得到的三角形与原有网格融合，能得到布尔运算后的网络模型。

按照上述算法，通过模型之间的布尔运算，将模型之间的拓扑关系再次回归到模型内部拓扑关系的构建及拓扑结构的访问。

使用本算法只可进行单个孔洞的剪裁，为了实现墙面多孔洞的复杂墙体，采用迭代法将完成的墙面作为新的S1代入布尔运算，完成墙面的拓扑重构。具体效果如图 8所示。

2 地下电缆工井拓扑关系构建

地下工井建筑中，各个工井墙壁及附属管孔之间通过电缆的连接具备了物理和逻辑上的关系。为了实现地下管网三维数据的组织、空间查询和空间分析等目的，需要建立一个可以有效描述网络系统的三维模型，从而对空间关系进行更为贴切的表达。

本文采用一种混合型的拓扑模型结构CSG-BREP^[24]（constructive solid geometry-boundary representation）模型对地下工井电缆网络进行构建。该模型综合CSG^[25]模型的宏观组合性以及BREP^[26]模型的微观表达性。在宏观上，CSG模型将地下工井划分为小对象，对复杂的空间目标进行几何抽象，将复杂物体分解成CSG基本体素，如墙壁、电缆、井筒等，构建地下工井的CSG树；在微观上，采用具体的参数方程对某一对象进行表达，便于进行后续计算。CSG-BREP模型的拓扑结构是顶点组成边，边组成环，环组成面，面组成壳（shell）、体（solid）或者复合体（compound），由较低层级的拓扑对象逐级构建较高层级的拓扑对象，中间不可跨级构建^[27]。该模型可以全面细致地记录模型内部的拓扑关系。

地下电缆工井模型总体设计以建筑墙面为单位，面上包含电缆管孔、沟道孔等结构，如图 9所示。由多个面构成工井主体，管孔模型依附于墙壁表面并且承担井内的电缆连接，工井主体叠加工井井筒等体元素组成地下工井三维模型。对于带有管孔的墙面，本文基于参数曲线和参数曲面的求交运算来进行工井模型面与圆柱体之间的布尔运算。该模型在保证拓扑一致性的约束下，采用渐进式的构建方法，由边构建环，再由环构建面，用面构建壳和体，获得复合体三维模型。

本文算法构建的地下工井模型可以全部按单元分解，做到模型部件化，并且在墙面上建造出管孔模型，使地下电缆工井模型内与模型间的拓扑关系更加完善，除了点线关系、线线关系，还添加了点面关系以及线面关系。

为了实现地下管网三维数据的组织、空间查询和空间分析等目的，需要建立一个可以有效描述网络系统的三维模型，从而对空间关系进行完备和形式化的表达。地下管网各个组成对象之间的层次关系图如图 10所示。

3 实验与分析 3.1 实验环境

本文使用Visual Studio 2013软件，采用C++语言进行算法编写，搭载在北京建筑大学数字建筑研究所PC_D3D团队自主研发的软件——地下电缆大数据智库平台（UP_3D V2.0）中，以单独的拓扑模块存在。本文实验所用数据使用FaroS150三维激光扫描仪采集，扫描原始数据为.fls数据格式，通过Faro自带软件SCENE可将数据转换为PTS或PTX格式。

3.2 具有孔洞的地下工井模型三维建模实验

本文算法实现主要体现在软件的管孔功能界面内，算法建模效果如图 11所示。

图 11中展示了多曲面工井整体建模效果，包括工井内部的电缆连接。图 11中左侧工具栏为单体工井模型的层次结构，分为工井模型和电缆模型两大部分，其中Part1工井模型节点下是以面为单位的子节点，具有孔洞模型的面节点下会增加以管孔为单位的子节点（即Part1中被圈部分）；而Part2电缆节点下是由面上管孔作为连接点生成的电缆模型，选中某电缆节点时界面中高亮显示，同时右侧Part3会显示该电缆模型的属性信息。

本文算法能够很好地应用于带有孔洞的建筑点云模型，避免点云残缺所带来的墙壁表现不完整，并且完整地建立出墙面的管孔模型，提供能满足后续进行井内电缆连接的管孔拓扑结构，满足实际生产以及快速三维可视化操作浏览的需要。图 12展示了通过本算法进行建模的工井模型（部分）以及通过管孔进行连接的电缆模型。对于多个单体工井模型通过管孔的对应关系可以进行模型间连接，生成井间管道线路，构成完整的地下电缆网络。对井内任意一根电缆进行拓扑关系查询，查询结果如图 13所示。在图 13中，Part1表示当前电缆段所在工井的名称以及连接的起始点与终点位置，位置信息中包含了管孔所属面以及管孔序号；Part2表示与该条电缆在同条线路的其余电缆段名称及位置信息，这些电缆段的拓扑连接关系依靠管孔生成。

由模型与实际点云建模效果可知，本文算法能很好地满足实际生产的需要。除工井之外，本文算法对规则直通点云数据，如规则平面的房屋、建筑等，也有很强的适用性。

3.3 实验结果分析

为了验证本文算法对具有孔洞的地下工井点云数据的建模可行性及优选性，本文通过多种建模方法对3种不同类型的墙壁数据（曲面点云数据（图 14（a））、矩形面点云数据（图 14（b））以及不规则平面点云数据（图 14（c）））进行建模，对模型结果进行比较。图 14中第1列为顶部点云数据与无孔洞模型，第2列分别为使用本文算法、OpenCASCADE、CGAL狄洛尼曲面重建及Geomagic建模方法。

本文选定的对比方法有计算几何库(computational geometry algorithms library, CGAL)中三维狄洛尼三角网的曲面重建方法^[28]、常用点云建模软件Geomagic中的隐式曲面建模^[10]以及三维CAD开源库Open CASCADE。CGAL中三维狄洛尼三角网的曲面重建方法能够重建地形扫描、建筑物扫描和精细化扫描的点云模型；Geomagic采用的Wrapping表面重建算法是一种隐式曲面算法，基于Morse理论对狄洛尼三角形复杂度的定义对单形聚类进行溃减操作^{[10, 28]}；Open CASCADE是一个提供三维曲面造型和实体建模、CAD数据交换及可视化的软件开发平台，可用于开发除CAD外的计算机辅助制造(computer aided manufacturing, CAM)及计算机辅助工程(computer aided engineering, CAE)相关的软件。

定性分析方面，从图 14的实验结果样例可以看出，对于同一种数据，采用Geomagic以及狄洛尼曲面重建，模型均具有不完整性，缺失的点云在建模时会出现破损，无法还原出整体墙壁；墙面孔洞无法区分建模，会出现孔洞无法识别等情况；且采用Geomagic以及CGAL狄洛尼曲面重建进行建模得到的墙壁，墙壁的孔洞缺少三维空间中管孔的实际含义，即缺少与墙面之间的拓扑关系，无法承担井内电缆连接示意点的责任；采用Open CASCADE的三维建模算法进行建模时，在矩形平面墙壁与多边形平面墙壁有较好的效果，但是在对曲面墙壁进行拓扑运算时会出现曲面墙壁扭曲的情况。而本文提出的算法在墙面整体建模上没有缺失，并且可以完整建立墙面的管孔模型，能提供满足后续井内电缆连接要求的拓扑结构。

4 结语

本文提出一种针对具有孔洞的地下工井三维模型拓扑重构算法，通过采用一种混合拓扑模型结构CSG-BREP，可从宏观与微观上全面细致地记录模型内部拓扑关系，最终建立具有孔洞的地下工井三维模型，可解决直接采用点云进行构网情况下模型不完整及速度较慢的问题；在表达模型拓扑关系上突破了传统管网中单一的管孔与管线的点线关系及管线间的线线关系，增加管孔与墙壁的点面关系以及管线与墙壁的线面关系，在模型描述上更加准确。通过多组地下工井点云数据进行测试，可知本文算法能达到地下工井模型三维建模的要求。对比现有常用的点云建模方法与三维CAD建模方法，本文算法建模效果更加贴合实际，满足用户需求。

下一步的目标是结合多种方法，如CorkMaster、Toptop等^[29-30]来对目前的方法进行改善，使得对具有流形或自交的模型的布尔切割运算更加稳定，同时提高连续多次布尔运算下的相交检测和切割的精度与鲁棒性。