技术领域
[0001] 本
发明涉及工业测量、钢结构检测领域，具体涉及一种钢构件圆柱体拟合算法。
背景技术
[0002] 在现代建造业中，人们对大型钢结构(尺寸在十米至数十米)的加工
精度要求越来越高，对建造技术提出了更高要求的同时，对建造
质量的准确评估和检验也有着迫切的需求。
[0003] 大型圆柱形刚钢构件广泛应用于大型
桥梁，涵洞，隧道，
船舶，航天，航空等领域，圆柱形钢构件的主要指标主要包括端面的平面度，柱体的圆度、挠度、直径等。相对于小型构件可用测量仪器直接测量，大型圆柱形钢构件通常都是固定的不可移动的物体，其直径往往非常巨大，且结构复杂，周围伴随着诸多遮挡物，测量空间也可能比较受限，往往很难在一个
位置完成测量，而多位置测量的数据如何进行同
坐标系转换又是一大难题。
[0004] 目前的大型圆柱形刚构件的测量一般采用撑杆法等几何测量的方式，这种方式除需要大量人
力外，检测慢，效率低，操作繁琐，主观性大，测量精度不高，自动化程度低，此外该方法需要整个构件处在开阔的场地，须保证其在某个位置能被整体测量，大大限制了测量条件。同时该方法存在的不确定因素多，即影响了建造效率，也无法对建造质量作出精确的评估。
[0005] 对于大型圆柱形刚构件传统方法也只能测量圆柱端面的形状，对于端面的平整度及整个柱体的圆度、挠度无法做出测量。因此提出一套完整的测量、分析、计算方法来满足大尺寸圆柱形钢构件高精度测量的要求。
发明内容
[0006] 有鉴于此，本发明提供了一种钢构件圆柱体拟合算法
[0007] 本发明提供了一种钢构件圆柱体拟合算法，包括以下步骤：
[0008] 步骤1：在待测一个大型钢结构的剖面圆上设置多个测量点，并布置高精度激光仪器采集测量点的三维坐标；
[0009] 步骤2：采用转站算法，把不同测量点测量的三维坐标转换到一个公共坐标系下；
[0010] 步骤3：利用步骤2转站后的坐标点计算拟合平面方程，得到平面方程参数，在拟合运算过程中，算法能自动剔除误差较大的数据，获取最佳拟合方程参数，以此平面作为剖面圆度的评定平面；
[0011] 步骤4：将测量点坐标投影到步骤3确定的评定平面上，并计算各点投影坐标以及到评定平面的距离；
[0012] 步骤5：利用平面圆拟合算法将步骤4投影到同一平面的坐标点作为圆周上的点进行平面圆拟合，得到圆的拟合半径和圆心；
[0013] 步骤6：重复步骤1-5，对多个剖面圆进行处理，在步骤5平面圆拟合的
基础上，将多个剖面圆采用最小二乘法进行柱体圆度计算即进行圆柱拟合。
[0014] 进一步地，步骤1中所述高精度激光仪器与计算机连接，实时传输所测数据给计算机进行处理计算。
[0015] 进一步地，步骤2中所述把不同测量位置测量的三维坐标转换到一个公共坐标系下的具体方法为：用两组对应的测量点数据计算出对应的转站关系参数，两组测量点数据中均有测量点编号进行对应，算法输入两组带测量点编号的空间三维数据后输出转站关系，其中包括有3个平移量和3个旋转量，采用奇异值分解算法，构造H矩阵，H矩阵为两组测量点数据去
重心化后的乘积，设目标坐标为PT(XT,YT,ZT)，待转换的坐标为PS(XS,YS,ZS)，H为3*3的矩阵，设P1(n*3)＝[XS,YS,ZS],(S＝0,1,2L n)，设P2(n*3)＝[XT,YT,ZT],(T＝0,1,2L n)，则H＝P1T·P2，求H矩阵奇异值，可直接求转换矩阵：R＝VUT，U和V为左右奇异矩阵，需要判断R的行列式值是否小于0，如果小于0则为反射矩阵，需要做转换，将V[0,3]的值乘以-1，再次计算R；n表示坐标点个数。
[0016] 设目标坐标点
云为Center2(x2,y2,z2)，待转换坐标为Center1(x1,y1,z1)[0017] 则平移向量为
[0018]
[0019] 有旋转矩阵后，可求得转换
角度，这里以欧拉角为例，旋转矩阵R如下所示：
[0020]
[0021] 可以求得旋转角为：
[0022] θx＝atan2(R32,R33)
[0023]
[0024] θz＝atan2(R21,R11)
[0025] 其中：θx为x的旋转角，θy为y的旋转角，θz为z的旋转角，Rij表示旋转矩阵第i行j列的值， 为自转角，θ为章动角，ψ为旋进角。
[0026] 进一步地，步骤3中所述平面方程参数的计算方法为：
[0027] 设平面方程为
[0028] ax+by+cz+d＝0
[0029] 由点云Pi(xi,yi,zi)、重心Center(x0,y0,z0)i表示点云pi中任意一点；
[0030] 得到去重心化的点云Pi(Xi,Yi,Zi)
[0031] 构造矩阵Dn*3[Xn Yn Zn],n＝0,1,2,3Kn
[0032] 求矩阵的奇异值USVT＝D.svd()
[0033] D的最小奇异值对应的奇异向量即为平面方程参数a、b、c
[0034] d＝-(ax0+by0+cz0)
[0035] 求坐标点到平面的距离，求中误差，设限差为3倍中误差，将剔除距离超过3倍中误差的点，然后重新拟合平面；
[0036] 其中，求点(x0,y0,z0)到平面的距离公式如下：
[0037]
[0038] 进一步地，步骤4中所述空间点投影到
指定的空间平面的方法为：
[0039] 由平行关系有下面投影方程：
[0040]
[0041] (x′，y′，z′)为各坐标点在投影平面上的投影坐标；
[0042] 由垂直关系求得测量点到投影平面的距离t
[0043]
[0044] 进一步地，步骤5中所述平面圆拟合方法如下：
[0045] 采用非线性最小二乘法，假设空间圆参数为观测值，初始值可以取上面的线性平差结果，或者圆心取重心，半径取点云到重心的距离平均值；
[0046] 根据距离公式构造误差方程
[0047]
[0048] 泰勒展开线性化，设
[0049] 展开后的法方程为V＝BX-l
[0050] 则B为[(x0-xi)/D,(y0-yi)/D,(z0-zi)/D,-1],i＝0,1,K,n
[0051] l为[R-D]
[0052] 求点云平面拟合方程，构造约束条件方程，其中
[0053] C＝[a b c]
[0054] Wx＝[-(d+ax0+by0+cz0)]
[0055] d为点(x0,y0,z0)到平面的距离
[0056] 得法方程
[0057]
[0058] X为方程解，Ks为联系数向量，BT、CT分别表示B、C的转置，解方程得[0059]
[0060] 解方程得圆心坐标的改正值(dx0,dy0,dz0,dR)
[0061] 用改正值修正初值，做
迭代计算，当X的改正值小于限差结束迭代，这里限差取10-6；
[0062] 半径R为各点到圆心的距离平均值，每一个坐标点到圆心的距离，再减去半径R就是它的挠度，平均挠度为所有挠度的平均值，最大挠度为所有挠度值中的最大值。
[0063] 进一步地，所述步骤6中采用最小二乘法进行柱体圆度计算的方法如下：
[0064] 取圆柱面上的点到轴线上的投影点的距离为圆柱半径，以此构造方程，具体方法如下：
[0066] x＝x0+at
[0067] y＝y0+bt
[0068] z＝z0+ct
[0069] 式中(a,b,c)为空间直线的单位方向矢量，(x0,y0,z0)为空间直线上离原点最近的点，t为直线上任意点到(x0,y0,z0)的距离，为了唯一表示直线，定义a>0，若a＝0则b>0；若a＝0且b＝0，则c>0；a，b，c不可能同时为0。
[0070] 过观测点Pi(xi,yi,zi)与中心轴线垂直的平面方程
[0071] ax+by+ch+Di＝0
[0072] 其中Di＝-(axi+byi+chi)
[0073] 将轴线方程带入得：
[0074] t＝-(ax0+by0+ch0+Di)＝a(xi-x0)+b(yi-y0)+c(hi-h0)
[0075] 中心轴线与该垂直平面的交点坐标：
[0076]
[0077] 观测点Pi(xi,yi,zi)与交点坐标的距离为：
[0078]
[0079] 得到各点误差方程为
[0080]
[0081] 其中，R圆柱为半径，7个参数分别为(a,b,c,x0,y0,z0,R)，在以下两个条件下求7个参数
[0082] a2+b2+c2＝1
[0083]
[0084] 令 误差方程线性化
[0085]
[0086]
[0087]
[0088]
[0089]
[0090]
[0091]
[0092] 其中
[0093]
[0094]
[0095]
[0096]
[0097] 初始坐标点可以取点云重心即所有坐标点的加权平均值，半径取坐标点到重心距离的平均值,解方程即得到7个参数的改正值，(a,b,c)取1，改正值即为初始值和迭代后值的差值，修正初值后迭代计算，迭代过程保证a、b、c不同时为0，当改正值小于限差10-6时，迭代结束。
[0098] 本发明提供的技术方案带来的有益效果是：通过本发明将多种算法整合为一体，对大型圆柱形钢构件尺寸测量具有很强的针对性，且适应各种实地建造环境，整套系统集成化高，便于携带，配套高精度激光测量手段，测量效率和测量精度都大大提高，解决大型圆柱形钢构件尺寸参数快速高精度测量的问题。
附图说明
[0099] 图1是本发明一种一种钢构件圆柱体拟合算法的功能结构图；
[0100] 图2是本发明一种钢构件圆柱体拟合算法的转站关系计算
流程图；
[0101] 图3是本发明一种钢构件圆柱体拟合算法的平面方程计算流程图；
[0102] 图4是本发明一种钢构件圆柱体拟合算法的投影计算流程图；
[0103] 图5是本发明一种钢构件圆柱体拟合算法的平面圆拟合流程图；
[0104] 图6是本发明一种钢构件圆柱体拟合算法的圆柱拟合流程图；
[0105] 图7是本发明一种钢构件圆柱体拟合算法的
软件主界面图；
[0106] 图8是本发明一种钢构件圆柱体拟合算法的计算转站关系界面图；
[0107] 图9是本发明一种钢构件圆柱体拟合算法的计算转站关系结果图；
[0108] 图10是本发明一种钢构件圆柱体拟合算法的平面圆度计算界面；
[0109] 图11是本发明一种钢构件圆柱体拟合算法的平面圆度计算结果图；
[0110] 图12是本发明一种钢构件圆柱体拟合算法的计算平面圆度限定半径的操作界面图；
[0111] 图13是本发明一种钢构件圆柱体拟合算法的圆柱拟合操作软件界面图；
[0112] 图14是本发明一种钢构件圆柱体拟合算法的软件圆柱拟合结果图。
具体实施方式
[0113] 为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。
[0114] 请参考图1，本发明提供了一种钢构件圆柱体拟合算法，包括以下步骤：
[0115] 步骤1：在待测剖面圆上设置多个测量点，并布置高精度激光仪器采集测量点的三维坐标；高精度激光仪器与计算机连接，实时传输所测数据给计算机进行处理计算。
[0116] 步骤2：采用转站算法，包括
正交矩阵法，奇异值分解法，单元四元数法，对偶四元数法4种转站算法，把不同测量位置测量的三维坐标，转换到一个公共坐标系下，如图2所示，具体方法为：用两组对应的测量点数据计算出对应的转站关系参数，两组测量点数据中均有测量点编号进行对应，算法输入两组带测量点编号的空间三维数据后输出转站关系，其中包括有3个平移量和3个旋转量，本
实施例采用奇异值分解算法，构造H矩阵，H矩阵为两组测量点数据去重心化后的乘积，设目标坐标为PT(XT,YT,ZT)，待转换的坐标为PS(XS,YS,ZS)，H为3*3的矩阵，设P1(n*3)＝[XS,YS,ZS],(S＝0,1,2L n)，设P2(n*3)＝[XT,YT,ZT],(T＝0,1,2L n)，则H＝P1T·P2，求H矩阵奇异值，可直接求转换矩阵：R＝VUT，U和V为左右奇异矩阵，需要判断R的行列式值是否小于0，如果小于0则为反射矩阵，需要做转换，将V[0,3]的值乘以-1，再次计算R；
[0117] 设目标坐标点云为Center2(x2,y2,z2)，待转换坐标为Center1(x1,y1,z1)[0118] 则平移向量为
[0119]
[0120] 有旋转矩阵后，可求得转换角度，这里以欧拉角为例，旋转矩阵R如下所示：
[0121]
[0122] 可以求得旋转角为：
[0123] θx＝atan2(R32,R33)
[0124]
[0125] θz＝atan2(R21,R11)
[0126] 其中：θx为x的旋转角，θy为y的旋转角，θz为z的旋转角，Rij表示旋转矩阵第i行j列的值， 为自转角，θ为章动角，ψ为旋进角。
[0127] 步骤3：利用步骤2转站后的坐标点计算拟合平面方程，得到平面方程参数，如图3所示，在拟合运算过程中，算法能自动剔除误差较大的数据，获取最佳拟合方程参数，以此平面作为剖面圆度的评定平面；
[0128] 设平面方程为
[0129] ax+by+cz+d＝0
[0130] 由点云Pi(xi,yi,zi)、重心Center(x0,y0,z0)，i表示点云pi中任意一点；
[0131] 得到去重心化的点云Pi(Xi,Yi,Zi)
[0132] 构造矩阵Dn*3[Xn Yn Zn],n＝0,1,2,3Kn
[0133] 求矩阵的奇异值USVT＝D.svd()
[0134] D的最小奇异值对应的奇异向量即为平面方程参数a、b、c
[0135] d＝-(ax0+by0+cz0)
[0136] 求坐标点到平面的距离，中误差取10-6，设限差为3倍中误差，将剔除距离超过3倍中误差的点，然后重新拟合平面；
[0137] 其中，求点(x0,y0,z0)到平面的距离公式如下：
[0138]
[0139] 步骤4：将测量点坐标投影到步骤3确定的评定平面上，并计算各点投影坐标以及到评定平面的距离，如图4所示，具体方法为：
[0140] 由平行关系有下面投影方程
[0141]
[0142] (x′，y′，z′)为各坐标点在投影平面上的投影坐标；
[0143] 由垂直关系求得测量点到投影平面的距离t
[0144]
[0145] 步骤5：采用非线性最小二乘法将步骤4投影到同一平面的坐标点进行平面圆拟合，得到圆的拟合半径和圆心，如图5所示，具体方法如下：
[0146] 假设空间圆参数为观测值，初始值可以取上面的线性平差结果，或者圆心取重心，半径取点云到重心的距离平均值；
[0147] 根据距离公式构造误差方程
[0148]
[0149] 泰勒展开线性化，设
[0150] 展开后的法方程为V＝BX-l
[0151] 则B为[(x0-xi)/D,(y0-yi)/D,(z0-zi)/D,-1],i＝0,1,K,n
[0152] l为[R-D]
[0153] 求点云平面拟合方程，构造约束条件方程，其中
[0154] C＝[a b c]
[0155] Wx＝[-(d+ax0+by0+cz0)]
[0156] d为点(x0,y0,z0)到平面的距离
[0157] 得法方程
[0158]
[0159] X为方程解，Ks为联系数向量，BT、CT分别表示B、C的转置，解方程得[0160]
[0161] 解方程得圆心坐标的改正值(dx0,dy0,dz0,dR)
[0162] 用改正值修正初值，做迭代计算，当X的改正值小于限差结束迭代，这里限差取10-6；
[0163] 半径R为各点到圆心的距离平均值，每一个坐标点到圆心的距离，再减去半径R就是它的挠度，平均挠度为所有挠度的平均值，最大挠度为所有挠度值中的最大值。
[0164] 步骤6：重复步骤1-5，对多个剖面圆进行处理，在步骤5平面圆拟合的基础上，将多个剖面圆采用最小二乘法进行柱体圆度计算即进行圆柱拟合，如图6所示，方法如下：
[0165] 取圆柱面上的点到轴线上的投影点的距离为圆柱半径，以此构造方程，具体方法如下：
[0166] 设中心轴线方程为
[0167] x＝x0+at
[0168] y＝y0+bt
[0169] z＝z0+ct
[0170] 式中(a,b,c)为空间直线的单位方向矢量，(x0,y0,z0)为空间直线上离原点最近的点，t为直线上任意点到(x0,y0,z0)的距离，为了唯一表示直线，定义a>0，若a＝0则b>0；若a＝0且b＝0，则c>0；a，b，c不可能同时为0；
[0171] 过观测点Pi(xi,yi,zi)与中心轴线垂直的平面方程
[0172] ax+by+cz+Di＝0
[0173] 其中Di＝-(axi+byi+czi)
[0174] 将轴线方程带入得：
[0175] t＝-(ax0+by0+cz0+Di)＝a(xi-x0)+b(yi-y0)+c(zi-z0)
[0176] 中心轴线与该垂直平面的交点坐标：
[0177]
[0178] 观测点Pi(xi,yi,zi)与交点坐标的距离为：
[0179]
[0180] 得到各点误差方程为
[0181]
[0182] 其中，R为圆柱半径，7个参数分别为(a,b,c,x0,y0,z0,R)，在以下两个条件下求7个参数
[0183] a2+b2+c2＝1
[0184]
[0185] 令 误差方程线性化
[0186]
[0187]
[0188]
[0189]
[0190]
[0191]
[0192]
[0193] 其中
[0194]
[0195]
[0196]
[0197]
[0198] 初始坐标点可以取点云重心即所有坐标点的加权平均值，半径取坐标点到重心距离的平均值,解方程即得到7个参数的改正值，(a,b,c)取1，改正值即为初始值和迭代后值的差值，修正初值后迭代计算，迭代过程保证a、b、c不同时为0，当改正值小于限差10-6时，迭代结束。
[0199] 考虑系统运行于微软Windows及Microsoft.NET Framework 4.0环境中，根据平台统一性和可重用性原则，使用C#语言进行开发，采用微软动态链接库(Dynamic Link Library)技术，实现在Windows
操作系统下共享函数库。同时实现模
块化，方便二次开发。最终以DLL(动态链接库)文件的方式形成成果，将各种算法封装于动态链接库中，供测量时使用。
[0200] 该程序用于调用测量辅助计算模块程序，并演示各种计算方法，同时输出各类计算结果，以便于使用者检验或者比较计算结果，主界面如图7所示，程序的文件数据格式为xyz坐标，以空格隔开每行一个点坐标。
[0201] (1)转站关系计算界面如图8所示
[0202] 读取数据，一次选中两个文件，打开，进行计算，会计算第一个文件的数据转到第二个文件的坐标系中参数。点击交换数据，再次计算可以反过来计算，计算结果如图9所示。
[0203] (2)平面圆度计算
[0204] 输入限差，会自动以该限差作为距离中误差的倍数进行计算，勾选求投影点和挠度，会自动调用计算方法。单击“读取数据”，打开文件，读取数据，操作界面如图10所示。
[0205] 勾选了显示绘图，会自动把结果绘制出来，勾选清空输出，每次计算会自动将上次输出的结果清除，关闭上次的绘图窗口。
[0206] 绘图中以挠度大小点的按从小到大
颜色以蓝色渐变至红色，如图11所示。
[0207] 限定半径的约束平差，需要输入约束半径如图12所示。
[0208] (3)圆柱拟合
[0209] 直接打开数据，进行计算，计算完后程序会自动绘制出图形，如图13所示，结果如图14所示。
[0210] 在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合，以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
