本发明属于空间信息技术领域，特别涉及一种顾及三维拓扑的地下工井设施建模方法。
背景技术：
三维激光扫描技术早期兴起于工业领域的逆向工程中，主要用于工艺品设计、汽车、飞机以及船舶等工业制造业领域，后来逐渐被应用到测绘行业领域。而在城镇地下管线的发展中，随着城镇地下管线信息化建设水平的提高，地下管线信息表达形式已由传统二维符号逐步提升为三维立体。管线在地下的分布纵横交错，二维图形无法表现管线之间的空间关系。有些管线上下起伏，与地面垂直的一段管线在平面图上只能以一个点及相应注记来表示，视觉效果不直观。因此，在国内外竞相研制三维数码城市、三维数字小区的氛围下，研究地下三维管线及地下工井模型势在必行，因此利用地下工井高精度的点云数据信息进行地下工井及电缆模型的构建很好的符合现在的研究趋势。
地下管线结构复杂、信息量大，管线三维可视化表达对管道及工井三维建模的效率和精度要求较高。国外学者burger与ermes以近景摄影测量(close-range-photogrammetry)技术与结构实体(constructivesolidgeometry，csg)重建三维管线模型。他们的研究对于显示化学工厂等地上错综复杂的管线很有效，但是难以重建地下工井模型。朱良峰等针对城市地下空间信息三维建模、可视化分析的特点和要求，设计了面向空间结构信息和空间属性信息建模与可视化的三维空间数据模型；琚娟等研究了基于b.rep进行地下空间一体化建模，并给出了基于预分块技术的城市地下空间对象空间索引方法。但是地下空间建模研究目前大多还集中在地下管线三维建模、地层建模，对地下工井以及工井之间管线拓扑关系考虑较少，扈震等提出管件设施曲面动态套合和模型动态变换两种三维模拟方法进行地下管线模型构建，但是构建的模型与实际有一定误差。唐相桢通过交互实现二维模型的三维转化，但是由于需要二维数据，以及交互化操作，因此有一定局限性。李清泉提出了分段渐次推算的模型计算方法，但是该文章只介绍了圆形管线的生成，缺少对工井模型的研究。
技术实现要素：
本发明的目的利用三维激光点云数据特点在充分考虑其空间拓扑关系的情况下进行地下电缆设施的三维模型构建。
本发明提供了一种顾及三维拓扑的地下工井设施建模方法，步骤包括：
1)点云精简；
2)提取出工井点云主轴；
3)将原始点云旋转，使得原始点云沿着主轴方向正立与三维空间中；
4)构建三角面片实现工井模型的简易化建模
a)工井顶面点云提取并进行平面拟合；
b)工井底面点云提取并进行平面拟合；
c)依据顶面边界线、底面边界线和顶面和底面之间的高度构建三角面片以进行贴合点云数据的简易化建模；
5)工井模型部件构建。
优选方案是：所述点云精简的步骤包括：
1)将点云模型进行空间网格划分；
2)在每个网格内依据给定距离阈值对点云模型沿着空间三个方向划分成多个小网格；
3)所述小网格按照其对应编码值按升序或降序排列。
优选方案是：利用obb包围盒方法获取轴向，并将最长边的方向作为几何体的主方向，对于主轴选取错误的情况根据精度信息进行主轴调整。
优选方案是：按照四元数旋转方法将原始点云旋转，使得原始点云沿着主轴方向正立与三维空间中。
优选方案是：通过计算点云的最高点zmax，以及最低点zmin，设置分层个数n，则每层距离为(zmax-zmin)/n，依次从顶层往下提取点云直到满足一定个数的点云为止。
优选方案是：采取多次平面拟合的方式进行顶部点云提取，具体方法如下：
1)通过ransac算法进行平面拟合，在参数提取时通过重心化点云坐标，实现平面拟合的优化；
2)根据拟合参数再在整体点云中提取与该平面一定距离阈值的点云，并逐渐缩小阈值进行多次迭代处理得到最终没有缺失的顶面点云数据。
优选方案是：所述构建三角面片以进行贴合点云数据的精细化简易模型建模步骤包括：
1)首先在边集合中找到一条邻接三角形数量为1的边，该边的两个端点分别作为多边形边界的前一点与当前点；
2)在当前点的邻接三角形中寻找邻接三角形数量为1的边且该边的两个端点有且只有一点与当前点或前一点重合；
3)端点中非重合端点作为多边形边界的后续点，通过递归方法依次寻找边界其他后续点，直到与第一点重合时结束此次边界提取；
4)进入下一次循环提取下一条边界直至所有边界提取完毕。
优选方案是：采用ransac算法进行空间圆拟合，通过设置半径阈值提取井口参数。
优选方案是：所述工井模型部件构建包括工井管孔构建和工井内部电缆构建，所述工井管孔构建的步骤包括：
1)依据原始点云数据与生成模型参数进行点云与模型的平面分离，
2)对于分离出的点云平面通过旋转投影，提取边界线，从而拟合管孔中心及半径；
3)根据上述步骤提取出的圆心及半径构建球，使球与模型平面进行拓扑重构以及剖分构建带有管孔的模型。
优选方案是：所述工井内部电缆构建方法为采用两个节点的三次hermite插值实现电缆途经点的插值从而优化电缆模型。
本发明的有益效果如下：
1)本发明提出的算法充分考虑工井模型的拓扑关系，较一般建模有明显改进；
2)本发明能够适应多种点云类型，算法适应性强；
3)本发明模型贴合实际点云，构建的模型质量好；
4)本发明对比常用的曲面重建算法本发明构建模型速度快，实用性较强；
5)本发明与以往管线建模相比，本发明实现内部管线与外部点云的有机统一，很好的满足实际场景的重现。
附图说明
图1为点云精简消耗的时间折现图；
图2为提取工井点云主轴效果图；
图3为提取工井点云主轴效果图；
图4为点云分层效果图；
图5为顶部点云平面拟合图；
图6为点云边界及关键点提取；
图7为构建三角面片生成模型过程图；
图8为工井管孔构建过程图；
图9为模型管孔生成过程图；
图10为空间实体求交点流程图；
图11为直线电缆生成图；
图12为电缆模型三角面片构建图；
图13为电缆模型生成图；
图14为不同精简方法效果图；
图15为不同工井模型构建效果图。
具体实施方式
下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
应当理解，本发明所使用的诸如“具有”、“包含”以及“包括”术语并不配出一个或多个其它元件或其组合的存在或添加。
本发明提供了一种顾及三维拓扑的地下工井设施建模方法，步骤包括：
1)点云精简
由于通过扫描获取的点云数据具有冗余量大特点，直接将原始点云用于建模将会耗费大量的时间和资源，因此一般在进行点云数据配准之前都要进行点云数据精简过滤工作。点云数据精简过滤最初的一些方法只是简单基于点之间距离、曲率、法向等原则，而目前点云精简的方法主要集中在以下几种典型方法：包围盒法、几何图像精简法、曲率精简法、法向精度精简法。
基于以上已有算法的分析总结，本发明采用一种可移动网格划分的点云精简过滤算法，该算法主要通过空间网格对点云进行精简，精简速度优势较明显，能够达到预期效果，且该算法其中二次网格精简准则能够将上述算法的精简准则纳入其中。
针对大数据量点云精简而言，如果直接依据给定距离阈值对点云模型进行空间网格划分，即将点云整个包围盒沿着空间三个方向划分成无数个小网格，分别在小网格内进行点云的筛选精简，通常为方便小网格的二分查找，小网格整体应按照其对应编码值按升序或降序排列，因此在进行每次二分插入时，通过二分查找找到相应位置后，所有该位置后的小网格数据都要往后推，涉及到大量数据的拷贝问题，从而产生大数据点云精简的性能瓶颈。如图1所示。随着原始点云个数的逐渐增加，精简消耗的时间越来越多，大大降低了精简的效率。为解决这一性能瓶颈，本发明采用了在划分小网格之前先进行一次网格划分，对点云模型首先进行分块处理，使得小网格置于一次网格中，大大减少了小网格二分插入的时间消耗。
2)提取出工井点云主轴
利用obb包围盒技术可以得到三个方向，根据边长的大小，可以筛选出几何体的主方向初始值。本发明采用obb包围盒方法来获取轴向，并将最长边长的方向作为几何体的主方向，效果见图2、图3。
设原始输入点集为cloud＝{p1,p2,p3k,pk}，任取一点pi∈cloud，设pc为点云的质心，即：
cv＝∑(pi-pc)(pi-pc)t(1)
其中(pi-pc)是一个列向量，(pi-pc)t为对应的转置后的行向量，cv是一个对称的半正定的3×3矩阵。
其中x,y,z为点各个点对应的坐标值。设cv的特征值中最小的一个为λ，最小特征值对应的特征向量为v，将其单位化处理，也就是将它的长度压缩为1，则法向n＝v或者n＝-v。法矢的方向可以根据视点位置进行调整。
设视点为vp(xp,yp,zp)，对于任意采样点p(x,y,z)，若满足下列方程，则需要对法矢的方向进行调整。
n·(p-vp)＜0(3)
也就是取
n＝-v(4)
以上方法也称为pca(principalcomponentanalysis)主成分分析法。对上述方法进行了改进，对每一个点p赋予不同的高斯权重。这样不同距离的点对法线的影响程度是不一样的，距离点p越近的点影响越大，从而得到更加真实的估计结果。
cv＝∑ui(pi-pc)(pi-pc)t(5)
其中：
最终提取出工井点云主轴。
3)将原始点云旋转，使得原始点云沿着主轴方向正立与三维空间中
根据公式(4)中obb获取的轴向n，利用公式(7)计算n与坐标轴z的夹角值。
cosθ＝n·n1(7)
其中：n1＝(0,0,1)，符号·表示两个向量的数量积。
按照四元数旋转方法，只需要一个旋转轴和一个旋转角度，即可完成空间内的任意旋转(绕任意轴，旋转任意角度)。其中，旋转角度即为公式(7)中求出的θ，旋转轴可由公式(8)求解。
dir＝n×n1(8)
其中dir为旋转轴向，n，n1与公式(7)相同，符号×表示两个向量的向量积，所得结果为一个矢量，同样把dir单位化。按照四元数旋转公式(9)完成旋转过程。
p1＝upu-1(9)
其中p表示待旋转的原始点云中的任意点，p1表示按照四元数旋转后的对应点，并且
利用以上四元数旋转方法，将旋转体旋转轴旋转至z轴方向，计算出其对应的旋转矩阵，使点云数据分别乘以上述旋转矩阵，从而旋转体沿着主轴方向正立于三维空间中。
4)构建三角面片实现工井模型的简易化建模
a)工井顶面点云提取并进行平面拟合；
本发明通过自适应方法进行顶部点云提取，通过计算点云的最高点zmax，以及最低点zmin，设置分层个数n，则每层距离为(zmax-zmin)/n，依次从顶层往下提取点云直到满足一定个数的点云为止，如图4、图5所示。由于在分离过程中不同数据的分布不同，所以可能部分顶面分离不完全，所以本发明采取多次平面拟合的方式进行顶部点云提取，具体方法如下：
本发明通过ransac算法进行平面拟合，在参数提取时通过重心化点云坐标，实现平面拟合的优化。由于分离出顶部点云可能有缺失的情况，所以在通过平面拟合的基础上，根据拟合参数再在整体点云中提取与该平面一定距离阈值的点云，并逐渐缩小阈值进行多次迭代处理得到最终没有缺失的顶面点云数据，如图6所示。
将点云数据向指定的平面投影。基于上述所确定的平面，我们将所有数据投影到该平面上。平面方程设为
ax+by+cz+d＝0(10)
a,b,c,d为平面四个参数，先计算每一点到平面的距离d为：
d＝(a,b,c,d)·(x,y,z,1)(11)
其中的(x,y,z)为每一点的坐标，(a，b，c)已经经过了单位化处理，对应的投影点pro为：
pro＝(x,y,z,1)-(a,b,c,0)×d(12)
b)工井底面点云提取并进行平面拟合；
依据同样的方法提取点云数据的下底面并进行平面拟合，通过两个平面的参数信息计算两个平面之间的距离，以上底面边界线为基准实现上底面边界点往下底面投影
c)依据顶面边界线、底面边界线和顶面和底面之间的高度构建三角面片以进行贴合点云数据的简易化建模；
在提取边界之前，需要先建立三角形、边、顶点两两之间的拓扑关系，而后通过递归方法实现边界的提取。使用循环过程，首先在边集合中找到一条邻接三角形数量为1的边，该边的两个端点分别作为多边形边界的前一点与当前点。然后在当前点的邻接三角形中寻找邻接三角形数量为1的边且该边的两个端点有且只有一点与当前点或前一点重合。端点中非重合端点作为多边形边界的后续点，通过递归方法依次寻找边界其他后续点，直到与第一点重合时结束此次边界提取。之后进入下一次循环提取下一条边界直至所有边界提取完毕。通过两个平面的参数信息计算两个平面之间的距离，以上底面边界线为基准实现上底面边界点往下底面投影。依据边界信息以及模型高度构建三角面片以进行贴合点云数据的精细化简易模型建模。对于顶部井口边界点，本发明采用ransac空间圆拟合的方法提取工井井孔的圆心与半径，从而与整体边界点一起构建三角面片实现工井模型的构建，构建过程如图7所示。
5)工井模型部件构建。
a)工井管孔构建
工井管孔是电缆连接的基础，本发明依据原始点云数据与生成模型参数进行点云与模型的平面分离，对于分离出的点云平面通过旋转投影，提取边界线，从而拟合管孔中心，整体过程如图8所示；
由于实际点云分离出的平面为空间任意平面，所以利用罗德里格旋转矩阵将该平面进行旋转至xoy平面，设空间平面的初始法线为n＝(nx,ny,nz)，旋转后法线是n'＝(0,0,1)则旋转角度θ为：
则旋转轴ω(ωx，ωy，ωz)为：
ω＝n×n′(14)
则旋转矩阵为：
点云旋转之后进行xoy平面上的投影，投影后进行边界线提取本发明同样用ransac的方法对所有提取边界线进行空间圆拟合设定圆半径阈值，从而提取管孔中心，如图9所示。在实际获取的圆心点与模型实际表面具有一定的误差，采用射线投影的方法求取圆心在墙面三角网的投影点。以投影点为基准，采用三维模型剖切算法对墙面三角网进行拓扑剪切，完成简易模型的二次剪裁。目前已有的三维模型切割方法都存在着一定的不稳定性，面对复杂的目标实体，比如说曲面墙面进行剖切获取门窗或者管道孔的时候会出现剖切错误。结合opencascade以及cork在模型剖切上的应用，对算法加以修改，实现稳定的适用于建筑表达的拓扑切割算法。尝试融合toptop剖切三角网的方法，来解决未闭合面、非流形边、有自交的复杂三角网的剖切问题。最终生成带管孔模型。
工井模型上的管孔是电缆连接的基础，为了满足逆向重建的需求，需要对工井模型进行布尔运算以完成拓扑重构。
布尔运算于1847年在英国被发明提出，其目的在于处理二值之间的关系[24]。布尔运算由两个部分组成，第一部分是布尔运算算子，分为求交运算(common)，求并运算(fuse)，求差运算(cut)；第二部分是空间中需要进行布尔运算的两个体素参数s1目标物体以及s2布尔运算工具。布尔运算的实现依赖于参数体素中隐含的几何要素，比如体素参数的曲面参数方程等。如表1所示，参与到工井模型拓扑重构中的几种要素参数方程
前文中已经获取了布尔运算的第一个体素参数s1墙面模型以及s2布尔运算工具的投影位置，可以通过用户设置的管孔半径(边长)或ransac拟合出的管孔半径(边长)生成布尔运算的被剪裁元素s2，与工井墙壁的几何模型进行拓扑重构。
算法需要输入a、b两个目标空间体素。遍历a中的所有边(参数曲线表示)，同时遍历b中的所有面(参数曲面表示)。将a中所有的边与b中所有的面求交点，获得由交点分割出的新边。再将b中所有边与a中所有面进行求交运算，得到交点。
由交点分割得到的边进行分类处理，经过体素求交点后，a中与b中面相交的边被交点分割成若干段。此时判断a中这些分割以后的边是在b的内部，外部，还是重合三种情况。同理，b相对于a采取相同的操作。
求解新的交线。在求解新的交线之前需要先确定有效交点。即如果求解的交点同时满足在a的边界上和在b的边界内(含边界)，则该交点为有效点。同理，判断b相对于a的有效点。在判断出所有的有效点之后，将这些有效点进行重新连线，得到新的交线。
需要判断a、b中原有边界的去留。根据边分类处理的分类结果以及连接的新交线来完成布尔运算。当a和b做求交运算时，则a中的边相对于b是内部的和边界的，作为保留，舍去相对于b外部边，同时对b相对于a采取相同的操作。当a和b做求并运算或者差运算时，a中的边相对于b是外部的，作为保留，舍去内部。将保留下来的边界，重新按照拓扑重建思想由边构建环，再由环构建面，面构建壳和体的方式构建出a和b布尔运算之后的结果s。
通过以上过程，即计算交点、求解新的交线、边界的去留、拓扑重组之后，便实现了相交关系的空间模型边界拓扑重构，减少了拓扑数据冗余。算法的如图10所示。
算法过程中的计算交点，需要将所有的参数曲线统一表示为：
c(t)＝(x(t),y(t),z(t))(16)
所有的参数曲面统一表示为：
f(u,v)＝(x(u,v),y(u,v),z(u,v))(17)
则参数曲线和参数曲面的交点由下面的算法确定：
将两空间实体采用参数方程表示，参数方程求交需满足：
c(t)-f(u，v)＝0(18)
令
s(t，u，v)＝c(t)-f(u，v)(19)
求参数微分方程：
将微分方程各参数微分显示，并在两边点乘
在得到迭代方程之后，对曲线和曲面分别构建包围盒树，当曲线的包围盒树与曲面的包围盒树叶子相交时，取两个叶子内的参数值作为初始值t0，u0，v0来迭代，直到满足收敛条件，最终求取交点。
按照上述算法，通过模型之间的布尔运算，将模型之间的拓扑关系再次回归到模型内部拓扑关系的构建及拓扑结构的访问。
使用本算法只可以对相面进行单个孔洞的剪裁，为了实现墙面多孔洞的复杂墙体，采用迭代的方法，将上一步完成的墙面作为新的s1带入到布尔运算中，反复多次，完成墙面的拓扑重构。
b)工井内部电缆构建
电缆是地下电力工井的重要组成部分，但是通过管孔之间进行电缆直接连接既不美观也会不可避免出现电缆模型与工井模型产生交叉的情况，如图11所示。因此本发明采用两个节点的三次hermite插值[29]实现电缆途经点的插值从而优化电缆模型。
在电缆模型构建之前，我们通过三维交互方式拾取电缆连接起始点，由于电缆起始位置位于工井模型平面上，所以我们由此可以得到电缆起始点pa＝(xa,ya,za)，pb＝(xb,yb,zb)所在平面的法线na＝(nx,ny,nz)，nb＝(nx,ny,nz)，则曲面法线dir＝na×nb。符号×表示两个向量的向量积，所得结果为一个矢量，同样把dir单位化。dir即为曲线所在平面的法线，以该平面为基准构建高斯平面直角坐标系x00y0，以pa，pb中点为坐标原点。
因此在点pa，pb转换之后的坐标为p'a＝(-dis/2,0,0)，p'b＝(dis/2,0,0)，其dis为两点之间的距离。则设旋转矩阵为r平移向量为t(tx,ty,tz)则：
p′a＝r·pa+t(24)
p′b＝r·pb+t(25)
令：
p'c＝(0,0,1)(27)
则根据三个对应点的关系可以反解出旋转矩阵r[30]。
令点集p的中心为c，点集p′的中心为c'。
点云重心化
pci＝pi-c(28)
p'ci＝p'i-c'(29)
计算矩阵h：
最后对矩阵h进行svd分解得到
h＝u∑vt(31)
u，v为分解出的矩阵，则3*3的旋转矩阵r如下：
r＝vut(32)
则平移向量t(tx,ty,tz)如下：
t＝c-rc'(33)
令dirx＝(1,0,0)，diry＝(0,1,0)，n'a＝(nx,ny,nz)，n'b＝(nx,ny,nz)则点a在高斯平面直角坐标系的斜率为：
则点b在高斯平面直角坐标系的斜率为：
符号·表示两个向量的数量积。
如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶，二阶甚至更高阶的导数值，这就是hermite插值问题，本发明采用节点处插值函数与函数的值及一阶导数值均相等的hermite插值，从而实现弯曲管道途经点的插值，从而优化管道曲线。
两个节点的三次hermite插值多项式为：
令x0＝-dis/2，x1＝dis/2，y′0＝ka，y′1＝kb则在x0～x1之间进行插值，设置插值个数为m则插值步长δx＝dis/m则插入点为：
p'i(xi,h(x0+iδxi),0)(37)
则三维空间中曲线插入点pi：
pi＝r-1p'i-t(38)
设置电缆半径，在电缆的起始点上生成电缆的横截面圆，通过插值点与截面圆sweep方法生成光滑电缆，根据曲线插值点构建三角网，如图12所示。根据生成的三角网构建通过direct3d进行三角面片的绘制，最终生成电缆模型，如图13所示。
实验效果
首先通过几种精简算法的对比，得出移动网格二次划分算法精简算法效率高，精简效果好。
本发明在使用移动网格二次划分精简，同时使用按密度自适应精简，按百分比精简进行对比，各种方法的精简效果如图14所示。a为原始点云、b为移动网格二次划分精简、c为按密度自适应精简、d为按百分比精简。为多种算法精简结果对比结果如表2所示。
表2多种算法精简结果对比
由以上表格可得移动网格二次划分精简消耗时间少，精简效果较好。为验证地下电缆模型生成效果本发明采用多种数据进行工井模型构建以及内部电缆模型构建，通过对多种数据的建模验证本发明算法的适用性。
一般点云建模算法都是基于三维狄洛尼三角网进行曲面重建，本发明通过多种模型生成算法进行建模对比，得到在地下工井模型构建中本算法建模速度快，建模精度高，如图15所示。不同模型生成时间及误差如表3所示。
表3不同模型生成时间及误差
实验结果分析
通过上述工井电缆模型生成实验可知，本发明算法能够很好应用于不同类型的点云数据，具有较广的适用性，且每次建模的时间消耗很小。由模型与实际点云中误差中误差可知，本发明算法与实际工井之间拟合精度足够高，能很好的满足实际生产的需要。
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。