0 引言

大部分油气输送管道是埋地管道，在长期运营之后，会由于地震、地层的移动、热输变形及第三方施工等原因形成凹陷、椭圆和凸起变形等管道缺陷^[1]。通过管道检测机器人定期对管道进行巡检并且对测得的数据进行处理分析，及时发现管道安全隐患并提前修复，对管道的安全运营显得非常重要。

管道的三维成像是处理变形检测数据最直观但也是最复杂的方式，而曲面的三维重建算法是实现管道三维成像的关键。国内外很多学者对管道的三维重建算法进行了研究。盛沙和R.L.HARDY等^[2-3]利用超声波内检测获取管道腐蚀曲面的离散数据点空间坐标，应用Multiquadric径向基函数插值法重构了管道的腐蚀曲面。闵令伟^[4]改进了迭代最近点(ICP)配准算法，并基于泊松重建算法对配准完的离散数据进行了曲面重建。闫勇刚等^[5]利用结构光扫描设备得到管道点云模型，采用主成分分析法(PCA)提取点云的数据特征，并最终拟合截面点数据，实现了管道的三维重建。该方法适用于估测大型管道的弯曲度，但较难解决管道局部变形的复杂问题。蒋必辞等根据测井信息，基于OpenGL对三维井眼轨迹进行了可视化研究^[6]。

Non-Uniform Rational B-Spline (简称NURBS)方法能用统一的数学形式表示规则曲面和自由曲面，且具有可影响曲线曲面的权因子，使形状更易于控制和实现^[7]。本文将NURBS曲面应用于管道的后期三维重建中，并借助OpenGL生成真实、高质量的管道三维模型，便于检测人员更加直观地了解管道变形情况。

1 管道检测数据获取

管道检测机器人主要有两种检测方式：接触式测量与非接触式测量^[8]。接触式测量主要通过探头测量管道内径；非接触式测量则具有许多检测方法，如激光检测法和超声波法等^[9]，它们都是通过测量管道截面的n个点来检测管道截面的变形情况。

管道检测机器人的定位有姿态定位与轴向定位。姿态定位一般通过对捷联惯导系统(SINS)采集的数据进行计算，得到三欧拉角(横滚角、俯仰角、偏航角)^[10]。轴向定位一般采用多里程轮的方式进行定位，然后通过地面检测标定盒或其他方式进行修正^[11]。

因此，本文所研究的管道检测数据为：管道截面n个点距管道机器人轴线的距离，截面三欧拉角(静态旋转)以及截面距前一个截面的距离。

2 基于NURBS曲面的三维重建算法 2.1 NURBS理论基础

NURBS是在B样条的基础上引入权因子，并用有理分式表示的一种曲线曲面表示形式。NURBS曲面通过节点矢量与控制点来控制曲面的形状，可通过少量的点拟合管道形状，方便灵活且适用范围广。NURBS曲面的方程定义为：

(1)

式中：ω_{i, j}为权因子，d_{i, j}为控制点，p为U方向次数，u为U方向节点值，q为V方向次数，v为V方向节点值，控制点有(n+1)(m+1)个，N_{i, p}(u)与N_{j, q}(v)为德布尔-考克斯递推公式。

(2)

本文采用的p与q均为3，所有的权因子都设为1。每个截面有n个型值点(首尾重合的点认为是两个点)，有m个截面，则U方向节点矢量有n+6个值，V方向节点矢量有m+6个值，控制点有(n+2)(m+2)个。由式(1)与式(2)可得：

(3)

式(3)实际为双三次B样条曲面方程，对该方程的计算主要有正算和反算两大类。正算是通过节点矢量与控制点计算曲面上各个型值点；反算是通过型值点计算节点矢量与控制点。首先通过检测数据得到若干型值点，由型值点反算曲面节点矢量与控制点，然后使用OpenGL通过节点矢量与控制点计算型值点，绘制管道的三维模型。

2.2 节点矢量的生成

节点矢量为一系列非递减值，一般由0到1，对应一条曲线的开始与结束，中间各值对应曲线上的一系列点。节点矢量分为U方向与V方向，定义U方向为管道截面周向、V方向为管道轴向。

考虑到管道截面周向n点均匀分布，且截面曲线为C²连续闭曲线，节点矢量的生成采用均匀参数化法，则U方向节点矢量如下：

(4)

在V方向，即管道轴向上，考虑到各截面间隔不均匀，节点矢量的生成采用积累弦长参数化法，即式(5)。前3个节点值为0，最后3个节点值为1，中间各节点值与该节点对应截面距初始截面累积弦长(截面间隔)成正比。

(5)

2.3 控制点的计算

控制点计算分为两步：①每条截面曲线控制点的计算；②轴向控制点的二次计算。

2.3.1 坐标系的定义

采用OpenGL的坐标系^[12]定义型值点的坐标系，如图 1所示。X轴自左向右，Y轴自下向上，Z轴自里向外，原点在屏幕中心。横滚角、俯仰角、偏航角分别为绕X轴、Z轴和Y轴旋转的角度。

2.3.2 截面控制点计算

由检测到的管道截面数据，可在Y-Z面得到截面各型值点坐标。考虑到管道机器人在检测过程中可能发生周向转动，因此测量了机器人的横滚角。将绕X轴旋转矩阵^[13]，即式(6)，与各个型值点坐标[x, y, z]^T相乘，即可得旋转后型值点的坐标。

(6)

每条截面曲线满足三次有理B样条曲线方程^[14]：

(7)

式中：d_i为控制点，N_{i, 3}(u)为德布尔-考克斯递推公式，即式(2)。

由于截面曲线首尾重合，所以有效型值点有n-1个。利用型值点与控制点关系可列n-1个方程，但有n+2个未知量(控制点)。因此，将截面曲线考虑为闭合曲线，首尾三控制点重合，使方程可解，如下式：

(8)

由于各截面周向旋转角度不同，所以轴向控制点线上各控制点不在周向同一角度，会在成像过程中发生扭转现象(见图 2)。因此，还需考虑对旋转后的曲线重新采样得到型值点。

对旋转后的曲线进行采样的方法是求曲线与若干条通过原点的直线的交点，对此国内外都有学者进行研究^[15-16]。这些算法具有一定的通用性，但是求解相对比较复杂。

本文研究的型值点比较特殊，u=0.0与u=0.5对应型值点连线经过原点，将闭曲线分成了两半，且经过原点的任意一条直线与曲线都仅有两个交点；在[0.0, 0.5)区间与[0.5, 1.0)区间，每半条曲线与Y轴和Z轴各有一个交点，求出曲线与坐标轴的交点，通过旋转曲线，即可求得曲线上任意角度处点的坐标，从而实现对曲线的均匀采样。周向旋转后的截面曲线如图 3所示。因此，本文采用一种新的采样算法，既简化了计算求解过程，又保证了较高的精度。

现以图 3中曲线在[0.0, 0.5)区间中与Y轴交点为例进行分析。以该交点为分界点，左半区间曲线Z < 0，右半区间曲线Z>0，特殊情况下，Z值等于0发生在边界处。因此可采用二分法不断逼近寻找该交点对应的节点值，继而可得该点坐标。同理，可得到曲线与坐标轴的4个交点。为了尽量保留曲线信息，应均匀地得到曲线上更多的点。考虑可旋转曲线(控制点)90°/n为n-1次，每次旋转得到在坐标轴上的4个点，根据旋转角度即可得到该点原坐标，如此可得到曲线上均匀分布的4n个点。

上述采样算法有一定的局限性，如检测点数量需为偶数、管道变形不能特别大(凹陷超过圆心)等。但对于大部分管道变形检测机器人而言，检测点数量一般为偶数，且管道变形一般为小变形，当发生特大变形时，管道实际已损坏，因此该算法适用于本文。得到曲线若干均匀分布型值点后，需再次反求出新的控制点为下一步计算做准备，通过得到的型值点信息也可得知该截面的异常信息，为生成检测报告做准备。

2.3.3 轴向控制点的计算

得到每条截面曲线在Y-Z面的控制点之后，需要将该曲线旋转移动到空间对应的位置。

在2.3.2节中已将截面曲线绕X轴旋转横滚角对应的角度，接下来要将截面曲线绕Z轴旋转俯仰角对应的角度、绕Y轴旋转偏航角对应的角度，即使用式(9)所示的绕Z轴旋转的旋转矩阵与式(10)所示的绕Y轴旋转的旋转矩阵乘以坐标即可。由截面与前一个截面的间隔与俯仰角、偏航角可得该截面与前一个截面的相对位置，迭代即可得当前截面(中心点)在空间中的坐标，将截面各控制点坐标进行平移处理即可将该截面移动到空间对应位置。

(9)

(10)

2.3.4 轴向控制点的二次计算

每条截面曲线在空间中的控制点坐标并非曲面的控制点坐标，若以此作为曲面控制点坐标，则会发现曲面形状并不能精确拟合各型值点。任意一个截面的控制点会对其他(临近)截面形状产生影响，所以需要在轴向上对控制点进行二次计算。

轴向控制点的二次计算过程为：在轴向上将截面的控制点当成型值点，再次计算控制点使得对应的曲线精确拟合各截面的控制点。在轴向上每条曲线有m个型值点(截面控制点)，通过其与m+2个控制点(曲面控制点)的关系可列m个方程，由于是开曲线，所以无法利用首尾重合的条件，这样就多出两个未知量。因此添加两个约束条件求解方程：开始两个控制点相同，尾部两个控制点相同。在进行轴向控制点的二次计算时，实际测量管道会比较长，无需对所有截面进行计算，只需计算屏幕上显示的截面即可，否则对计算机计算性能会有很高的要求。

3 基于OpenGL的管道三维成像

由前文反算所求得的曲面节点矢量与控制点需进行正算以生成管道曲面。这一过程无需编码实现，只需OpenGL提供相应的函数，利用节点矢量与控制点即可生成NURBS曲面。

管道三维成像流程大致可分为3步^[17]：首先通过gluNewNurbsRenderer()构建NURBS渲染器nurbs；然后通过gluNurbsProperty()设置渲染器模式等参数；最后以NURBS曲面节点矢量与控制点等为参数，通过函数gluNurbsSurface()完成曲面的显示。

为了取得较好的显示效果，分别使用线框模式与填充模式生成曲面，即将第二步与第三步执行两遍，并可以使用与计算控制点相同的方法计算颜色控制点，使用填充模式生成曲面后再次执行第三步，通过颜色控制点与求值器GL_MAP2_COLOR_4为生成的管道曲面设置颜色。

4 算法验证与结果分析

使用C++编程语言实现上述算法后，通过模拟数据生成了凹陷管道的三维模型(见图 4)，对算法进行了验证。具体验证与结果分析如下。

模拟数据每个截面采集8个点，周向旋转后二次采样得16个点。曲面上的点为要拟合的二次采样后的型值点。从图 4可以看出：各型值点都在生成的三维模型表面上，生成的弯曲管道三维模型上凹陷变形清晰可见，达到了预期效果。

常见管道变形主要有两种：凹陷变形与挤压变形。通过模拟数据生成了两种常见变形管道的三维模型，并对比了同种变形管道下，是否考虑管道机器人检测过程中可能发生的周向转动对图像结果产生的影响，结果如图 5~图 8所示。从图 5~图 8可以看出：基于NURBS曲面的三维重建算法对于管道局部变形的表面重构精度较高，其成像效果清晰、逼真，且当考虑检测机器人的周向转动时，生成的曲面过渡更加圆滑、显示效果好。

5 结论

(1) 针对管道变形检测机器人所测的数据，研究了基于双三次NURBS曲面的管道三维成像技术。

(2) 通过模拟常见变形管道的检测数据，利用C++编程语言和OpenGL生成了变形管道的三维模型，并对比了相同变形管道中，是否考虑检测机器人的周向转动对于成像效果产生的影响。研究结果表明：本文所采用的算法准确、有效，成像效果清晰、逼真；当考虑机器人的周向转动时，所生成的三维模型表面过渡更加圆滑。

(3) 本文研究的管道三维成像算法具有一定的通用性，稍加改变即可适用于各种管道变形检测机器人所测数据，但也有一定的局限性，如曲线采样算法要求机器人采得的截面点数量为偶数。