0 引　言

船体外板的成形质量极大影响着船舶的航行性能和安全性能。船体复杂曲面外板的成形需要经过冷弯加工和水火弯板。成形检测是水火弯板工艺中重要的工序，用于判断成形情况，并指导二次加工参数预测，直接影响外板的成形质量。目前这道工序通常由熟练的工人使用样板、样箱进行，存在经验性强、自动化程度低、检测效率及精度较低等问题^[1]。

在船体外板成形检测领域，国内外学者进行了全面详细的研究。Hiekata等^[2]开发了基于三维激光扫描仪的外板曲率偏差评估系统，通过获取目标外板点云数据与CAD模型进行比较，求得肋骨线曲率偏差，用于指导工人二次加工。实现了非接触式的检测，且摆脱了木制样板通用性不足的问题，但其求得的曲率偏差难以被工人应用于实际的工作场景。史一帆等^[3]采用摄影测量技术对成形船舶曲板进行图像获取，基于深度学习方法推断深度图得到曲板表面点云，并通过主成分分析法和迭代最近点法进行点云自动配准，但该方法得到的点云数据精度较低，配准结果难以被直接应用，并且未解决传统配准方法带来的点云穿透问题。赵亮等^[4]提出基于结构光辅助扫描的计算机双目视觉测量方法，通过提取图像中的激光中心线坐标，获取鞍形板表面的点云数据，测量系统的稳定性和精度较高，但受到激光扫描范围的限制，不适用于测量大型复杂曲板。时丹^[5]提出一种基于双目视觉的三维测量技术，该技术通过改进张氏标定法提高了检测的精度，并采用多阈值分割处理方法显著减少了噪点数量，但该技术的图像预处理时间较长，无法完成实时快速的检测工作。

当前的研究主要集中在使用三维激光扫描仪进行成形检测，该检测方式耗时较长，且对工作环境要求较高^[6]；并且现有的检测方法存在通用性不足、成形度评估效率不高的问题，难以适应船厂复杂的应用场景。

本文引入六自由度机械臂搭载线激光传感器进行测量，可实现高精度大范围的灵活扫描，测量精度最高可达0.04 mm；基于先进的特征点识别算法获取外板边界信息；将外板点云数据进行融合边界几何信息约束的ICP配准处理，并针对点云穿透问题，进行基于序列二次规划法的点云不穿透配准；基于Delaunay三角剖分曲面重建后，通过提取数字化样板从测量模型和设计模型中获取肋位点数据；基于样条插值的数学方法进行成形度评估，生成检测报告。该方法将检测精度提高到1 mm以内，相比传统方法显著改善了检测效率，可用于指导外板二次加工参数预测。

1 船体外板曲面点云数据获取研究 1.1 测量系统

本文的测量系统组成为ABB机器人、机器人控制柜、上位机及操作软件、交换机、线激光传感器及夹持器。在上位机建立服务器，以机械臂为客户端，传感器同为服务器，交换机作为通讯渠道。通过交换机将机械臂、上位机和传感器进行网口连接；基于C#语言编写传感器与机械臂之间通讯程序，建立上位机与传感器之间的端口连接，最终实现机械臂与传感器的信息交互功能。

1.2 TCP获取

线激光传感器检测的坐标信息需要转换成机器人末端执行器的坐标信息。在工业机器人领域，需要在机械臂末端法兰盘安装工具来进行作业，为了确定该工具的位姿，需要在工具上绑定工具坐标系TCS（Tool Coordinate System）,TCS的原点就是TCP（Tool Center Point）^[7]。为了实现传感器识别的坐标信息与机械臂TCP的坐标信息进行转换，并获取三维点云数据，需要在上位机端通过编程实现，实现流程如图1所示，分为原图获取、原图处理和TCP处理3个模块。编程采用C#语言，并结合ABB机器人专用编程软件Robot Studio实现数据信息交互功能。

2 船体外板边界信息获取及不穿透配准方法研究 2.1 边界信息获取研究

外板的边界信息能够帮助提高设计模型和测量模型点云数据的配准精度，并且能够实现对检测线的定位。相较于传统的三维检测方式，本文采用的测量方法具有更高的灵活性，可以实现对外板边界处点云信息的精准获取。

采用先进的边界特征点识别算法，对每条扫描线进行特征点计算并输出，可以得到外板边界上的点云数据集，见图2。由于线激光传感器与机械臂的采样运算频率不同，导致角点位置处的点云信息不完整，因此需要设计算法求取角点坐标。算法步骤如下：

1）对边界点云数据集进行降噪处理。对于偏离原边界的噪点，基于曲率值进行筛选，在曲率变化较大的位置处将该点去除。

2）样条插值处理。基于三次样条插值曲线将边界用数学形式表达。

3）基平面投影处理。将4条边界的函数方程投影到XY基平面，该平面为外板水平放置的平面。

4）求交点。将各边界投影的交点作为输入插值方程的参数，将求得函数值的平均值作为角点。

2.2 基于边界几何信息的点云配准方法研究

点云配准在步骤上分为粗配准和精配准^[8]。传统的粗配准方法将测量模型直接与设计模型进行匹配，由于数据量庞大导致耗时较长。本文将设计模型与边界点云模型进行初始配准，有效提高了粗配准效率，并且节约了精配准搜索匹配点对的部分时间。粗配准基于以下3个流程：重心配准、对角线向量叉乘配准及中点连线配准。

精配准采用点到点（Point–To–Point）的迭代最近点算法（Interative Closest Point，ICP）。该算法是一种通过迭代优化求出2组点云之间最佳变换关系的算法，算法根据2个点集P和Q，在P中寻找与Q点集中最接近的点进行配对，将所有配对点之间的欧式距离之和作为目标函数，求解目标函数最小时对应的旋转和平移矩阵，具有精度高、通用性强且简单直观的特点。目标函数可以表示为最小二乘问题^[9]，即：

$ E(R,T) = \mathop {\arg \min }\limits_{} \sum\limits_{i = 1}^n {{W_i}} {\left\| {R{q_i} + T - {p_i}} \right\|^2}。$

(1)

式中：R和T分别为旋转矩阵和平移向量；W_i为权重；q_i和p_i分别为2个点集中的点。求解该式子常用到奇异值分解法（Singular Value Decomposition，SVD）求解，即将R和T分别视为独立的变量进行计算。

计算平移向量时，先求出E关于T的导函数，并根据点集间的质心位置关系进行方程的简化。即：

$ \frac{{\partial E}}{{\partial T}} = \sum\limits_{i = 1}^n {2{W_i}} (R{q_i} + T - {p_i}) = 0 。$

(2)

令C_p和C_q分别为两点集的质心，根据${C_p} = \displaystyle\frac{1}{n}\sum\limits_{i = 1}^n {{p_i}} $，${C_q} = \displaystyle\frac{1}{n}\sum\limits_{i = 1}^n {{q_i}} $，代入式(2)，得：

$ E(R,T) = \mathop {\arg \min }\limits_{} \sum\limits_{i = 1}^n {{W_i}} {\left\| {R(q{}_i - {C_q}) - ({p_i} - {C_p})} \right\|^2}。$

(3)

令${k_{qi}} = ({q_i} - {C_q})$，${k_{pi}} = ({p_i} - {C_p})$，代入式(3)求目标函数最小时R的值。

$ E(R,T) = \mathop {\arg \min }\limits_{} \sum\limits_{i = 1}^n {{W_i}} {\left\| {R{k_q} - {k_p}} \right\|^2}。$

(4)

根据求得的旋转和平移矩阵对点集P进行变换，得到新的点集P'，即：

$ {{P}}' = R{{P}} + T。$

(5)

对新的点集进行最近点配对生成配对点集，并根据式(1)～式(5)进行迭代计算处理，直到生成的点集与基准点集的间距满足误差要求为止，误差阈值一般取0.1 mm，见图3。

2.3 点云不穿透配准方法研究

由于ICP算法本身并没有考虑物体表面的实际物理定义，仅对点云数据进行迭代计算转换矩阵，因此配准的结果出现设计模型和测量模型相互穿透的问题，不符合实际的工程情况。在实际工程中，设计模型即检测样板，需要放置在测量模型即加工后船体曲板的表面进行检测，因此需要对点云不穿透配准方法进行研究，使得检测方法符合实际工程应用情况。

不穿透配准问题的目标是，在满足设计模型在测量模型上方的前提下，求解一个使得2组点云间距离最小的坐标转换矩阵。该问题本质上是一个非线性不等式约束的优化问题，可以表示为：

$ \left\{\begin{gathered}\min f(x)，\\ \mathrm{s.t.}g_i(x)\leqslant0,i=1,\cdots,n。\\ \end{gathered}\right. $

(6)

其中：f(x)为目标函数，即点云之间的距离差值的平方和；g(x)为约束条件，表示设计模型点云的高度值始终大于测量模型点云高度值。

本文采用序列二次规划法（Sequential Quadratic Programming，SQP）求解，该方法适用于平滑非线性约束的优化问题，在拟合和非线性回归模型构建方面表现更好，优化结果精确度高，适用于中等规模的优化问题。SQP的基本原理是将复杂的非线性约束优化问题转化为二次规划问题进行求解。二次规划问题即目标函数为二次函数，约束函数为线性函数的优化问题，其基本思想是利用泰勒级数展开非线性约束问题的目标函数，并将约束条件简化为线性函数，设当前的迭代点为x^k，对式(6)中的目标函数进行二阶泰勒展开，并简化约束条件可得：

$ \left\{\begin{aligned} & {\rm{min}} f(x^k)+\nabla f \left(x^k\right)^{\mathrm{T}}+\displaystyle\frac{1}{2}d^{\mathrm{T}}\nabla^2f(x^k)d，\\ & \mathrm{s.t.}g_i(x^k)+\nabla g_i\left(x^k\right)^{\mathrm{T}}d\leqslant0,\ \ i=1,\cdots,n。\end{aligned}\right. $

(7)

其中：d =x–x^k代表步长。配准流程如图4所示。

从经ICP配准的曲面中选取2组三维点云数据，设计模型Q（红色）包含158个数据点，测量模型P（蓝色）包含175个数据点，二者在YZ平面上的初始位置关系如图5所示，经检测存在44个穿透点。进行基于SQP的不穿透配准计算后，将得到的最优解代入坐标转换矩阵，如图6所示。穿透点数量为零，证明不穿透配准算法能有效处理三维点云不穿透配准问题。

3 数字化样板检测技术研究 3.1 数字化样板提取研究

模拟工人使用木制样板、样箱的检测方法^[10]，本文基于Delaunay三角剖分曲面重建方法处理配准后的点云模型，提取设计模型中的数字化样板对测量模型待检测位置进行检测。

Delaunay三角剖分是点云数据重建的关键技术手段，对于一个给定三维空间的离散点集，采用改进的分治算法，按照空外接圆性质、最小角最大化性质构建拓扑关系，见图7。

该方法相较于直接对离散点云数据处理的方法，更加符合外板曲面的数学特征，这使得切割得到的点云数据更加准确且平滑。根据输入的截面位置、截面法线向量，获取设计模型数字化样板的点云坐标。程序步骤如下：

1）输入截面位置信息。包含截面原点三维坐标、截面法线向量。

2）获取三角面片及顶点索引。从面片索引中获取顶点索引，进而获取顶点坐标索引。

3）根据截面信息筛选三角面片。首先将各三角面片在垂直于截面方向的投影求出，例如，根据截面位置x=x₀，则将各三角面片顶点在X方向的3个坐标x₁、x₂、x₃进行排序，根据关系x₀小于或等于三者的最大值，并且x₀大于或等于三者的最小值，得到经过截面的三角面片。

4）计算截面与三角面片交点。针对直角三角形，若截面与直角边相切，则输出直角边对应的2个顶点坐标，否则根据相似三角形法则计算交点。针对非直角三角形，根据分类结果，根据空间坐标计算方法计算交点坐标，见图8。

3.2 数据处理方法研究

对于获取到的设计模型和测量模型在待测位置处的点云数据，为了获取任意位置处的检测点偏差值，本文采用三次样条插值处理。三次样条插值是在由两相邻节点所构成的每一个小区间内用三次多项式来逼近，在节点处具有C2连续性，即相邻两段曲线在节点处具有相同的一阶导数和二阶导数，该方法对于大多数工程问题具有较好的光顺性。三次样条函数的基本形式为：

$ {S_{{i}}}({{x}}) = {a_i} + {b_i}x + {c_i}{x^2} + {d_i}{x^3}。$

(8)

n段子区间对应着4n个未知数，节点处的数值及导数条件为：

$ \left\{ \begin{gathered} {S_{{i}}}({{{x}}_{{i}}}) = {y_i}, \\ {S_i}({{{x}}_{{{i}} + {\text{1}}}}) = {S_{i + 1}}({{{x}}_{{{i}} + {\text{1}}}}) = {y_{i + 1}}, \\ S_i^{'}({{{x}}_{{{i}} + {\text{1}}}}) = S_{i + 1}^{'}({{{x}}_{{{i}} + {\text{1}}}}), \\ S_i^{''}({{{x}}_{{{i}} + {\text{1}}}}) = S_{i + 1}^{''}({{{x}}_{{{i}} + {\text{1}}}})。\\ \end{gathered} \right. $

(9)

此时已知条件包含（4n–2）个方程，还需要加入边界条件，本文采用非扭结（Not–A–Knot）边界条件，即两端点的三阶导数与近邻点的三阶导数相等，即：

$ \left\{ \begin{gathered} S_0^{'''}({{{x}}_0}) = S_1^{'''}({{{x}}_1})，\\ S_{n - 1}^{'''}({{{x}}_n}) = S_{n - 2}^{'''}({{{x}}_{n - 1}})。\\ \end{gathered} \right. $

(10)

求解方程组可得各段子区间上的对应的方程系数。求得插值函数后，输入区间范围、步长相同的坐标数组，求得坐标数组对应插值。

$ \begin{split} & \left[ {\begin{array}{*{20}{c}} { - {h_1}} & {({h_0} + {h_1})} & { - {h_0}} & 0 & \cdots & 0\\ {{h_0}} & {2({h_0} + {h_0})} & {{h_1}} & 0 & \cdots & \vdots \\ 0 & {{h_1}} & {2({h_1} + {h_2})} & {{h_2}} & 0 & \vdots \\ \vdots & 0 & \ddots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & {{h_{n - 2}}} & {2({h_{n - 2}} + {h_{n - 1}})}&{{h_{n - 1}}}\\ 0 & \cdots & \cdots & { - {h_{n - 1}}} & {({h_{n - 2}} + {h_{n - 1}})}&{ - {h_{n - 2}}} \end{array}} \right] \cdot \\ &\left[ {\begin{array}{*{20}{c}} {{m_0}}\\ {{m_1}}\\ {{m_2}}\\ {{m_3}}\\ \vdots \\ {{m_n}} \end{array}} \right] = 6\left[ {\begin{array}{*{20}{c}} {\rm{0}}\\ {\displaystyle\frac{{{y_3} - {y_2}}}{{{h_2}}} - \displaystyle\frac{{{y_2} - {y_1}}}{{{h_1}}}}\\ \vdots \\ {\displaystyle\frac{{{y_{i + 2}} - {y_{i + 1}}}}{{{h_{i + 1}}}} - \displaystyle\frac{{{y_{i + 1}} - {y_i}}}{{{h_i}}}}\\ 0 \end{array}} \right]。\end{split}$

(11)

式中：h_i=x_i+1–x_i，m_i=s_i"(x_i)。

为了评估设计模型与测量模型的差异，需要计算各插值点处的挠度差值和曲率差值。挠度差值根据设计模型和测量模型在插值点处的差值求得。为了准确反应检测线中各点的弯曲度变化，采用平均曲率公式进行曲率差值计算，计算公式为：

$ K = \frac{{\left| {{y^{''}}} \right|}}{{{{\left( {1 + {y^{'2}}} \right)}^{\frac{3}{2}}}}}。$

(12)

3.3 船体帆形板实测案例

本文选用经线加热工艺成型的帆形板，尺寸为1200 mm × 800 mm × 15 mm。采用40 mm/s的扫描速度以及100 Hz的采样频率进行扫描，扫描时长3.95 min，获取了1120705个数据点。

成形判别规则参考2016年发布的《中国造船质量标准》(CSQS)^[11]，对于双曲度板，规定了肋位方向与样箱的空隙允许极限为5 mm，标准范围≤4 mm；在长度方向与样箱的空隙允许极限为5 mm，标准范围≤3 mm。

3.3.1 横向成形检测

选取的横向检测线如图9所示，分别求出各条检测线处对应肋位点的平均曲率差值及挠度差值，基于挠度偏差求出成形度R_ω，计算公式为：

$ {R_\omega } = 1 - \frac{{\left| {{\omega _P} - {\omega _Q}} \right|}}{{{\omega _{std}}}}。$

(13)

式中：ω_P、ω_Q分别为P、Q点集检测点的挠度；ω_std为该位置处成形所需的挠度。

分别求出各条检测线上的200个等间距点处的挠度和曲率差值，见表1和表2。记录超出允许极限的点，根据各条检测线检测点的成形度求出对应检测线的成形度，得到各条检测线的检测报告，见表3。

统计数据显示，最大挠度偏差值为4.39027 mm，低于肋位方向与样箱的空隙允许极限，最大偏差发生在检测线11的检测点200位置处，分析数据发现该检测线的成形度为97.58%，仅在靠近板边位置存在较大偏差。14条检测线中有12条检测线的成形度在95%以上，平均横向成形度为96.77%，高于船厂经验规则，证明该板横向成形情况良好。

3.3.2 纵向成形检测

选取的纵向检测线如图10所示，检测方法同上。分别计算检测线上400个等间距点的挠度和曲率差值，见表4和表5，并记录超出允许极限的点，得到各检测线的检测报告，见表6。最大纵向挠度偏差值为3.53693 mm，低于长度方向与样箱的空隙允许极限；最大纵向挠度偏差发生在检测线7的点206位置处，检测线7平均成形度为94.01%，平均挠度偏差值为3.24683 mm，高于标准范围3 mm。分析数据可知其附近区域的检测线5、检测线6均存在挠度偏差较大的情况，因此需要对该区域进行二次加工处理。

3.3.3 大型板件成形检测可行性分析

受到实验条件的限制，本文选用的实验板件为小尺寸板。然而在船厂实际的生产过程中，板件的长度和宽度一般在几米到十几米不等，多为大尺寸板。为了分析本文采用的检测方法在大型板件检测中的可行性，通过小尺寸板的实验案例，对大型板件检测情况进行预测。

本文测量系统的机器人型号为IRB 4600–40/2.55，其最大工作半径为2.9 m，并且可以搭配长度为10 m的导轨使用，在装配夹持器及传感器后，测量系统工作半径可延伸至3 m以上。本文测量系统的检测范围，在延导轨方向为15.8～16.2 m，在垂直于导轨方向为3～3.4 m，测量覆盖面积为41.7～52.2 m²。因此本文的测量系统完全可以实现对船厂大尺寸板的检测工作。如板件尺寸超出检测范围，参考汽车零部件加工流水线，可以通过多个机器人协同作业来完成工作。

假设板件的尺寸为12000 mm × 3000 mm × 30 mm。为了减少测量时间，将扫描速度提高到100 mm/s，采样频率提高到120 Hz，基于实验案例预测的扫描时间约为59.25 min，获取大约2千万个数据点。对于实验案例中的小尺寸板件，用三维激光扫描仪扫描的时长约为24.5 min，可以预测对于假设的大尺寸板件，采用三维激光扫描仪的方法测量的时长约为367.5 min，相比之下本文的测量系统扫描效率高出83.88%。上述预测结果证明，本文的测量系统能够实现对大型板件的成形检测，并保证较高的效率。

4 结　语

本文针对船体外板成形检测中传统手工卡样板方法效率低、精度不足的问题，提出一种基于六自由度机械臂与线激光传感器的数字化检测方法，得到的结论如下：本文通过机械臂协同线激光传感器进行高效扫描，可快速获取外板曲面高密度点云数据，最高密度可达每平方米1250000点；提出的基于边界信息的粗配准方法和点云不穿透配准方法有效解决了复杂曲面匹配效率低及点云穿透问题，使得检测结果符合实际工程情况；基于改进Delaunay三角剖分的曲面重构技术结合数字化样板方法，可灵活提取任意数量的检测点信息，检测精度优于1 mm；结合成形挠度与曲率分析快速生成可视化检测报告，检测效率较传统方法提升75%以上，且检测结果可直接指导二次加工参数预测。研究成果为船体外板数字化检测提供了可靠的技术方案。