本发明涉及视觉重建
技术领域：
，具体涉及一种破损圆柱排水管道内壁三维重建方法。
背景技术：
：城市地下排水管道主要用于输送生活污水并及时排除雨水，在城市中起到了至关重要的作用。管道通常埋于地下，常处于潮湿的环境中，很容易发生腐蚀、裂缝等缺陷，排水系统故障引发的路面塌陷、污水管道外溢、城市内涝等频频发生。为了有效预防此类灾害，获取管道病害数据，并及时监测评价变的尤其重要。目前，在工程中，普遍采用重建的方法评估病害指标，这些重建方法主要分为接触式检测方法和非接触式检测方法。接触式检测方法，需要人为控制测量仪，费时费力。非接触式方法主要包括是基于激光扫描法和机器视觉法，激光扫描法利用激光扫描管道内表面得到管道表面点云数据，但由于实际管道内结构复杂，获取的点云数据量庞大，计算时间长，噪声大，不利于分析、人机交互。机器视觉法对管道内壁拍摄的序列图像进行特征提取、特征点匹配以及进行基础矩阵估计算法实现管道侧壁的三维重建，该方法需要对每幅图像提取特征、匹配，因此计算量很大。为了解决上述问题，如何在获取管道内壁三维结构的同时，又能快速地三维重构实现多视角观察，改善人机交互体验，获取描述地下排水管道病害情况的关键指标。技术实现要素：本发明的目的在于提供一种破损圆柱排水管道内壁三维重建方法，该方法可用于地下排水管道内壁病害检测，其利用三维显微测量技术和摄像机标定技术获取并处理数据源，利用小孔成像模型、三维点云数据的插值拟合技术、OpenGL成像技术，将局部病害数据源重建为管道内壁纹理平面图，利用所有病害纹理重构图重建整个管道三维场景，解决了现有管道内壁三维重建数据采集复杂、三维模型重建速度慢，人机交互体验差等问题，重建后的内壁纹理平面图对病害评测指标分析提供了更有意义的参考价值，工作人员可以更加直观的多视角的观测病害。为了达到上述目的，本发明采用了以下技术方案：一种破损圆柱排水管道内壁三维重建方法，包括以下步骤：S110设备准备阶段：利用张正友标定方法对CCD摄像机标定，获取不同摄像机倍率下对应的摄像机内参矩阵M1、畸变参数，建立摄像机倍率与M1之间的映射表，并将摄像机装在管道机器人机身，机器人控制摄像机实现绕101轴的轴向旋转、绕102轴的水平旋转；S120数据采集阶段：利用管道机器人拍摄管道内壁序列图像，实时记录管道机器人拍摄图像时的管道半径R，行进位移w、摄像机与管道底部内壁的距离h、以及摄像机的轴向旋转角度β、水平旋转角度α和摄像机倍率m；S130数据处理阶段：从S120拍摄的排水管道序列图像中挑选每个管道内壁病害区域的序列图像，利用管道内壁有效病害区域占所在照片的百分比和聚焦评估函数为每个有效病害区域选取一张高质量的二维图像，并对筛选后的每幅图像利用S110得到的畸变参数校正，得到每个病害的校正图像；S140搭载数学模型阶段：为获取三维重构中管道内壁病害部分的纹理，针对S130获得每个病害的校正图像建立对应的小孔成像模型，每个模型包含3个参数：摄像机内参矩阵M1、摄像机外参矩阵M2、以及摄像机坐标系矩阵Zc，其中摄像机内参矩阵M1是通过S120中记录的拍摄图像时摄像机的倍率求取，外参矩阵M2是通过S120中记录的机器人行进位移、摄像机轴向旋转角度与水平旋转角度、摄像机相对管道底部内壁的距离h共同求取，摄像机坐标系中Zc是通过内参矩阵M1、摄像机轴向旋转角度、水平旋转角度、已知的管道半径R、摄像机相对管道底部内壁的距离h共同求取，求解建立的小孔成像模型，获取图像坐标系中(u,v)对应世界坐标系中空间坐标(Xw,Yw,Zw)，得到空间点数据；利用上述空间点数据求解需要重建的管道内壁病害区域的边界，将管道内壁病害区域投影至重建像素坐标系(o-u′,v′)，利用基于高斯核函数的加权插值技术，填充重建像素坐标系的各点像素，作为管道三维模型中病害区域纹理数据；S150：选取指定管道内壁图像作为管道三维模型中非病害区域的纹理数据；S160：利用病害区域纹理数据以及非病害区域的纹理数据，借助OpenGL重建三维管道场景。优选地，所述步骤S110中，机器人控制摄像机实现绕101轴的轴向旋转、绕102轴的水平旋转，包括：在摄像机组成的摄像机坐标系中绕z轴旋转；在摄像机坐标系中绕y轴旋转，分别沿着这两次旋转对应的坐标轴正半轴观察原点，水平旋转角度为(Xc,Yc,Zc)，轴向旋转角度为Zc，其中，α的旋转角度范围为[-90°,90°]，β的旋转角度为[0°,360°]。优选地，步骤S130数据处理阶段分为以下两个步骤：S131:获取管道内壁序列图像中病害区域占该图像的百分比，滤除百分比低于阈值0.6的管道内壁图像。S132：对S131获取的图像序列用离散改进型拉普拉斯算子对病害区域进行聚焦评估，选取评估值最大的一张图像，离散改进型拉普拉斯算子为：其中，I(x,y)表示像素点(x,y)处的亮度，step表示像素点(x,y)的邻域范围；S140所利用的数学模型为小孔成像模型，小孔成像模型利用了摄像机内参矩阵M1、摄像机外参矩阵M2、摄像机坐标系(Xc,Yc,Zc)下被摄物体的Zc值、摄像机旋转角度α和β、机器人行进位移w、摄像机与管道底部内壁的距离h，建立了二维图像坐标系(u,v)像素值与世界坐标系(Xw,Yw,Zw)坐标值之间的线性方程，所述线性映射方程如下：其中，摄像机内参矩阵M1是与摄像机倍率相关，摄像机外参矩阵M2＝g(w,h,α,β)，摄像机坐标系中Zc＝h(u,v,w,h,α,β,M1)，g(·)为外参估计函数，h(·)为深度估计函数；利用小孔成像获取的空间点数据，获取病害区域的边界，该边界可以描述为：由与管道母线平行的一个平面及与管道母线垂直的两个平面，这三个平面与管道内壁所围的封闭曲面。S140中基于高斯核函数的加权插值技术主要分为以下几个步骤完成：a)将S140描述的重建像素坐标系(o-u′,v′)分为两个区域，由(o-u,v)像素点描述的为A区域，否则为B区域；b)B区域像素用常量替换，管道图像A区域利用分水岭变换进行图像分割为若干个边界区域Ai；c)对每个Ai区域内的像素点，利用待插值点周围像素进行插值计算；d)若待插值点周围像素均属于同一区域Ai，则直接利用基于高斯核函数的加权插值计算待插值点像素；e)若待插值点周围像素不属于同一区域，则判断待插值点所属区域Ak，滤除周围非Ak区域像素，然后再利用基于高斯核函数的加权插值计算待插值点像素。附图说明下面结合附图和实施例对本发明进一步说明：图1为本发明实施例中基于局部病害管道内壁三维重建的详细流程图；图2为本发明实施例中实现摄像头两种独立旋转搭载管道机器人示意图；图3为本发明实施例中摄像机坐标系与世界坐标系关系；图4为本发明实施例中重建像素坐标系中各区域的描述示意图；具体实施方式为了能更清楚地理解本发明的上述目的、特征和优点，以下结合附图及实施例，对本发明进行进一步详细说明，以下实施例不构成对本发明的限定。图1描述了本发明的实施例中的一种破损圆柱排水管道内壁三维重建方法的具体流程，具体步骤如下：S110设备准备阶段：准备CCD摄像机及用于标定的黑白棋盘图，在每种视场下获取不同摄像机倍率的棋盘图，实现一次采样多次标定，将采集的图像按照摄像机倍率分组，利用MATLAB内张氏标定的函数和按摄像机倍率分组的图片集，求取摄像机不同倍率下摄像机5个内参数(fx,fy,u0,v0,s)以及5个畸变参数：三个径向畸变参数(k1,k2,k3)和两个切向畸变参数(k4,k5)。将内参数转化为内参矩阵，记为M1,建立摄像机倍率与内参矩阵之间的映射函数:M1＝L(m)。M1中s、u0、v0取各标定结果的均值，m为摄像机倍率，fx，fy分别表示图像水平和垂直的尺度因子，u0和v0表示在像素坐标系下图像坐标原点的坐标，s是描述两个坐标轴倾斜角的参数，M1结构如下：如图2所示，将标定好的摄像机装载，使机器人可以控制摄像头完成两种相互独立的旋转。首先将摄像机做第一层包装，如图2中110所示，使摄像机可以绕102轴旋转，即绕水平旋转，旋转角度范围为[-90°,90°]，再将110装载如机器人中，使摄像机绕101轴旋转，即绕轴向旋转，旋转角度范围为[0°,360°]。S120数据采集阶段：将管道机器人置于管道中，编码测距系统开始工作，记录当前管道内径半径R，操作机器人行走录制视频数据，行进路径与管道母线平行所采集的数据为有效视频数据，将有效视频数据区间转为对应视频的区间。在行进过程中同步记录行进距离编码器、轴向旋转和水平旋转的的角度编码器数值、摄像机相对管道底部内壁的距离h、以及摄像机倍率m，并将这五种记录值转化为关于视频时间的函数，分别为w＝C1(t),α＝C2(t)，β＝C3(t)，m＝C4(t)，h＝C5(t)。并记录沿途检测到每一个病害在视频中采样的时间区间，以病害为基准建立病害文档，在有效视频数据中选取每个病害的采样区间，并转为图片数据得到管道内壁序列图像，存贮在对应的病害文档，并在各病害文档中建立每幅图像对应w、α、β、m、h数值记录条。对每个病害进行S130至S132相同操作。S130数据处理阶段包括以下子步骤：S131：标记每幅管道内壁图像的有效病害区域，有效病害区域D用管道内壁图像中的矩形区域表示，并记录D中左上角p1(x1,y1)与右下角p2(x2,y2)像素坐标，计算该区域像素占整幅图像的百分比p。如果p＜0.6，则滤除该图像，如果p≥0.6，则保留该图像，建立p1、p2与管道内壁图像名称之间的对应关系。S132：对S131获取的图像用离散改进型拉普拉斯算子作为聚焦评估函数，选取评估值最大的一张图像作为最佳病害图像。上述离散改进型拉普拉斯算子为：其中，I(x,y)表示像素点(x,y)处的亮度，step表示像素点(x,y)的邻域范围。上述拉离散改进型拉普斯算子所使用的普拉斯模板为：0101-41010聚焦评估函数为：S133：将S132选取的最佳病害图像利用S110得到的畸变参数(k1,k2,k3,k4,k5)校正，得到每个病害的校正图像。S140：搭载数学模型阶段：将S133获取的校正图像进行如下操作：S141：获取校正图像的拍摄时间，利用S120创建的函数m＝C4(t)，获取校正图像对应的的摄像机倍率，利用S110获取的函数L(m),获取校正图像对应的内参矩阵M1。M1为：图3所示机器人在管道采样过程中，各个坐标系之间的关系，以及各测量值的几何关系。由于柱形管道切面为圆形具有对称性，以及S120中设定的病害图像在视频中采样时间区间，假定小车行进过程、行进方向始终与管道母线平行。211(Ow-Xw,Yw,Zw)构成世界坐标系(三维笛卡尔坐标系)，满足左手坐标系，其中起始位置摄像头所在管道的切面圆圆心为原点Ow,Xw轴为水平方向,Yw轴为垂直方向,Zw轴与管道母线平行方向为机器人行进方向。201(Oc-Xc,Yc,Zc)为没有发生旋转由摄像机构成的摄像机坐标系(三维笛卡尔坐标系)，各坐标轴与世界坐标系平行。设沿着摄像机坐标系的坐标轴正半轴观察原点，那么绕坐标轴的逆时针旋转为正方向旋转。202、203分别为摄像机发生水平旋转和轴向旋转摄像机坐标系变化示意图，虚线为变化后坐标轴位置，α为绕Zc旋转的角度，即水平旋转角度，β为绕Yc旋转的角度，即轴向旋转角度。h为摄像机相对管道底部内壁的距离，R为管道半径，w为机器人(摄像机)行进的位移量。S142：利用S120中的w＝C1(t),α＝C2(t)，β＝C3(t)，m＝C4(t)，h＝C5(t)，计算拍摄校正图像时，摄像机的行进位移w，摄像机相对管道底部内壁的距离h，摄像机绕图2中102旋转的水平旋转角度α，摄像机绕图2中101旋转的轴向旋转角度β。那么，外参矩阵M1和M2的平移向量记作(0,h-R,w)T，旋转矩阵C由水平旋转和轴向旋转生成，记作C＝B·A。其中，A为图2中绕102轴旋转的旋转矩阵，B为图2中绕101轴旋转的旋转矩阵,A、B、C的表示如下：机器人的状态可以描述为先进行平移，然后再进行旋转。由外参矩阵M2的平移向量确定四维平移矩阵，由旋转矩阵C确定四维旋转矩阵，那么，新的描述机器人当前状态的外参矩阵M2可以四维平移矩阵和四维旋转矩阵共同确定，表示如下：S143:利用S120中的α＝C2(t)，β＝C3(t)，m＝C4(t)，h＝C5(t)，计算摄像机相对于管道底部内壁的距离h，摄像机绕图2中102旋转的角度α，摄像机绕图2中101旋转角度β。由图3机器人状态分析可知，摄像机水平、轴向旋转过程相当于圆柱管道的等效旋转过程，具体如下：在摄像机坐标系(Xc,Yc,Zc)下分别绕Zc轴、Yc轴分别旋转-α、-β角度，那么在管道为圆柱曲面的约束条件下，管道的等效旋转矩阵为C′记作：在摄像机坐标系中，没有旋转的管道参数方程为：其中，x、y、z分别表示管道在圆柱坐标系下的参数，d表示管道厚度，t表示沿着管道z轴方向的长度，θ表示柱坐标系下自x轴逆时针旋转的角度。根据上述摄像机水平、轴向旋转过程相应的管道等效旋转矩阵C′，在计算摄像机坐标系中，管道按照旋转矩阵C′旋转后的方程为：(cosαcosβx′-sinαy′+cosαsinβz′)2+(sinαcosβx′+cosαy′+sinαsinβz′-d)2＝R2，(1)其中，x′、y′、z′分别表示管道等效旋转后在圆柱坐标系下的参数。利用透视投影知识以及S141中的M1可得到管道三维模型中重建像素坐标系与摄像机坐标系之间的关系如下：联立以上公式(1)和公式(2)，得出关于Zc的一元二次方程如下：其中，如果在公式(2)的等号两侧同时加负号，则拍摄视场就与当前摄像机视场关于(O-Xc,Yc)平面对称，这样会导致上述一元二次方程会出现双解，此时求出的视场并不是我们需要的，因此求得z′应为：S144：利用S141、S142、S143建立每张二维管道内壁图像的像素点与三维空间坐标值之间的线性映射关系，关系式如下：其中，S145:求解S144线性方程，由管道内壁有效病害区域内二维图像像素值，计算出对应的三维空间坐标值，得到有效病害区域的三维空间点数据。S146:根据S145求得的有效病害区域的三维空间点数据，计算病害区域边界，该步骤具体包括如下子步骤：1)提取管道内壁图像中矩形区域的四个顶点的像素坐标，查找S145中求解的点云数据，获取四个顶点像素坐标对应的世界坐标系中的空间坐标{(xi,yi,zi)|i＝1,2,3,4}；2)选取(xi,yi,zi)中zi的最小值zmin＝min(zi)和最大值zmax＝max(zi)，{(xi,yi,zi)|i＝1,2,3,4}在图2中世界坐标系(Ow-XW、YW)平面的投影点为{(x′i,y′i,0)|i＝1,2,3,4}；以原点Ow为极点，XW轴为极轴建立极坐标系，将投影点转为极坐标{(ρi,θi)|i＝1,2,3,4}，式中ρi＝R，θi＝atan2(y′i,x′i)，求取θi的最小值θmin＝min(θi)和最大值θmax＝max(θi)。若θmax-θmin≥180°则θ′max＝θmin，θ′min＝θmax-360°，由(zmin,zmax,θ′min,θ′max)确定最后需要拟合的管道内壁边界。S147：重建管道三维模型中病害区域纹理数据。图4展示了重建纹理的图像坐标系与像素坐标系之间的关系，并将重建像素坐标系(O-u′,v′)划分为A、B两个区域。A区域为圆弧封闭区域，由校正图像的采样像素决定，并将A区划分若干区域Ai，Ai区域像素利用基于高斯核函数的图像插值技术获取，对于插值区域为区域内部的插值点“×”，及插值区域穿越边界的插值点“+”，分别用不同的插值办法。B区域各位置像素用常量替换。具体步骤如下：1)由S146描述的管道内壁病害区域投影至二维平面(O-x′,y′)，投影后的平面图为长zmax-zmin，宽的矩形区域，将管道内壁病害区域点云数据投影至(O-x′,y′)平面。2)设重建图像分辨率为M×N。(M、N均大于CCD摄像机拍摄图像的横纵分辨率)3)将(O-x′,y′)分割为(M-1)×(N-1)的网格，构成重建像素坐标系(O-u′,v′)，网格内对应的点云数据的像素，作为网格的整体像素，这样将(O-x′,y′)坐标系中点云数据对应的像素填充至像素坐标系(O-u′,v′)。若点云数据为n个，则此时(O-u′,v′)坐标系中有n个像素点。A区域像素用RGB分别为142、142、142替代。4)利用分水岭变换将S133获得的该病害区域的校正图像分解为若干个互不重叠的区域(Di|i＝1,2,..,N′)，N′为分解的区域的个数，并利用Di所包含的像素，及1)中的投影关系，在(O-u′,v′)坐标系的A区域中找到对应区域Ai。5)在(O-u′,v′)中，对于任意点(i,j)，以(i,j)为中心进行扩大搜索，搜索步长为每次在横纵方向各增加1，直至搜索区域内像素数不少于2，设搜索到的像素坐标为{(u0,v0),(u1,v1)…,(un,vn)}，增加区域标记维度信息mi，则对应三维像素坐标为{(u0,v0,m0),(u1,v1,m1)…,(un,vn,mn)}。6)若m0＝m1＝…＝mn，则对(i,j)点的像素，利用基于高斯核函数的加权插值，具体步骤如下：步骤a:将搜索区域内的RGB像素空间转换为HIS像素空间，对于H分量和S分量均采用搜索区域内像素对应的分量均值替代，对于I分量利用步骤b和步骤c进行插值计算。步骤b:计算上述搜索到的横纵坐标方差，分别为二维高斯函数为：步骤c:(i,j)点的插值像素为I(i,j)＝α0I(u0,v0)+α1I(u1,v1)+L+αnI(un,vn),其中，7)若存在a、b使得ma≠mb，则对(i,j)点的像素，先利用众数和距离判别插值点所在区域，滤除区域外像素点，最后利用基于高斯核函数的加权插值，具体步骤如下：步骤a1:扩大搜索区域，直至搜索区域内像素数不少于4；步骤b1:若搜索区域内标记而且两个坐标轴倾斜角的参数与点云数目比满足s/n≥0.75,则认为该插值点所属区域也为k。若无满足此比率的标记，计算搜索区域内Ai中所有像素点与插值点距离的均值，均值最小的区域为该插值点所属区域，并标记该插值所属区域为k。最后得到插值点对应的三维像素坐标为(i,j,k)，并滤除搜索区域内所有标记非k的像素。步骤c1:将搜索区域内的RGB像素空间转换为HIS像素空间，对于H分量和S分量均采用搜索区域内像素对应的分量均值替代，对于I分量利用步骤d1和步骤e1进行插值计算。步骤d1:计算上述搜索到的横纵坐标的方差分别为二维高斯函数为步骤e1:则该点的插值像素为I(i,j)＝α0I(u0,v0)+α1I(u1,v1)+…+αnI(un,vn),其中，S150:以RGB分量分别取值为142、142、142，创建非病害区域的纹理数据。S160:按照obj模型中数据的描述结构，并利用S146得到病害区域边界的描述和S147获取的病害区域纹理数据以及S150创建得到的非病害区域的纹理数据，创建管道模型obj模型，借助OpenGL展现整个的obj管道的三维场景，从而实现在虚拟空间对病害区域多角度观察。当前第1页1 2 3