2018年全球海上最大起重能力的万吨级重型盆座式起重机制造完成，该起重机的基座部分直径达30 m，采用3种厚度的钢板(60 mm、80 mm、90 mm)焊接制作而成。在尺寸测量控制方面，在焊接之前和焊接之后都需测量周长以指导生产，且完工后的实测周长将作为船厂下料切割钢板的参考，尽量确保起重机基座与船厂基座周长相等，减少两基座对接时上下接口半径方向上的错位量。传统方法是采用围尺法^[1]，用卷尺直接绕测周长，其误差有：垫片间直线替代圆弧误差的直圆误差；卷尺上下起伏引起的不同面误差；卷尺拉力松紧引起的拉力误差；卷尺标定误差；焊缝隆起误差；圆度误差。该方法精度较低，工作量大，不宜用于大尺寸工业测量。针对这种十分少见的超大同心异径圆圆周长的测量，本文提出了基于极坐标的测量方案。

1 极坐标法 1.1 概述

采用极坐标法测量工件圆的内壁点坐标，将圆心到内壁点的距离作为极径，通过板厚差改正将极径从内壁点延伸到理论圆上，计算内壁点的极角。首先，将极坐标系转换成笛卡尔直角坐标系；然后，采用最小二乘法拟合计算圆心，拟合圆心作为极坐标系的新原点；再次，迭代计算极径及圆心，直到相邻两次迭代圆心差值小于设定值。圆心确定后，首先，分别采用相邻极径及圆弧夹角计算同一分段圆弧长，若两计算值互差大于设定误差，则采用坐标放样的方法找出该段圆弧的两端点；然后，内插测点，并再次计算圆心、弧长等；最后，圆心不再或较小地变化，且各分段圆弧长误差符合要求，此时圆周长等于各分段圆弧长之和。

1.2 公式推导及实施步骤 1.2.1 极坐标法计算圆上各点

为简化以后计算，以概略圆心(或任意某点)为坐标系原点，第一点为X轴建立工业坐标系，概略等距地采集工件圆内壁测量点坐标(x_{w_i}, y_{w_i})并记录各测量点板厚改正值。当然也可根据现场情况任意加点。

将内壁点沿极径方向延长到圆，则圆上各点坐标如下：

$ \left\{ {\begin{array}{*{20}{l}} {{x_i} = {P_i} \times {\rm{cos}}{A_i}}\\ {{y_i} = {P_i} \times {\rm{sin}}{A_i}}\\ {{P_i} = \sqrt {{{({x_{{w_ - }i}} - {x_0})}^2} + {{({y_{{w_ - }i}} - {y_0})}^2}} + {w_i}}\\ {{A_i} = {\rm{ta}}{{\rm{n}}^{ - 1}}\frac{{{y_{{w_ - }i}} - {y_0}}}{{{x_{{w_ - }i}} - {x_0}}}} \end{array}} \right. $

(1)

式中，(x_i, y_i)为工件圆上任意一点坐标；(x₀, y₀)为工件圆心坐标；A_i为内壁测量点极角；P_i为极径；w_i为测量点板厚差改正(即工件内壁测量点到理论圆的距离)，可直接从设计图纸上量取。

1.2.2 圆心拟合

圆的表达式为：

$ {(x - {x_0})^2} + {(y - {y_0})^2} = {R^2} $

(2)

式中，R为圆半径。式(2)展开成以下形式：

$ {x^2} + {y^2} - 2x{x_0} - 2y{y_0} + x_0^2 + y_0^2 - {R^2} = 0 $

(3)

令

$ \left\{ {\begin{array}{*{20}{l}} {a = - 2{x_0}}\\ {b = - 2{y_0}}\\ {c = x_0^2 + y_0^2 - {R^2}} \end{array}} \right. $

(4)

则圆的方程式可等价转换为参数a、b、c有关的方程为：

$ {x^2} + {y^2} + ax + by + c = 0 $

(5)

只需求解出参数a、b、c，就可得圆心(x₀, y₀)及半径R。

令

$ f(a,b,c) = {x^2} + {y^2} + ax + by + c $

(6)

一般地，圆上某一点(x_i, y_i)到圆心(x₀, y₀)的距离d_i为：

$ {d_i} = \sqrt {{{({x_i} - {x_0})}^2} + {{({y_i} - {y_0})}^2}} $

(7)

由最小二乘法原理可知，当各点距离与半径差的平方和最小时，即

$ \sum {{{({d_i} - R)}^2}} = {\rm{min}} $

(8)

此时的x₀、y₀、R为圆心及半径的最或然解。

由于式(8)无解析解^[2]，在实际应用中，为减少运算量，将式(8)简化为：

$ \sum {{{(d_i^2 - {R^2})}^2}} = {\rm{min}} $

(9)

将式(7)代入式(9)有：

$ \sum {{{[{{({x_i} - {x_0})}^2} + {{({y_i} - {y_0})}^2} - {R^2}]}^2}} = {\rm{min}} $

(10)

结合式(2)、式(5)，将式(6)带入式(10)：

$ \sum {{{(x_i^2 + y_i^2 + a{x_i} + b{y_i} + c)}^2}} = {\rm{ Min }} $

(11)

即∑[f_i(a, b, c)]²=min，令：

$ F(a,b,c) = \sum {{{[{f_i}(a,b,c)]}^2}} $

(12)

对函F(a, b, c)分别求a、b、c的偏导数，并令其极值等于0：

$ \frac{{\partial F}}{{\partial a}} = \sum 2 (x_i^2 + y_i^2 + a{x_i} + b{y_i} + c){x_i} = 0 $

(13)

$ \frac{{\partial F}}{{\partial b}} = \sum 2 (x_i^2 + y_i^2 + a{x_i} + b{y_i} + c){y_i} = 0 $

(14)

$ \frac{{\partial F}}{{\partial c}} = \sum 2 (x_i^2 + y_i^2 + a{x_i} + b{y_i} + c) = 0 $

(15)

对式(13)~式(15)化简并展开为：

$ \sum {(x_i^2 + y_i^2)} {x_i} + a\sum {x_i^2} + b\sum {{x_i}} {y_i} + c\sum {{x_i}} = 0 $

(16)

$ \sum {(x_i^2 + y_i^2)} {y_i} + a\sum {{x_i}} {y_i} + b\sum {y_i^2} + c\sum {{y_i}} = 0 $

(17)

$ \sum {(x_i^2 + y_i^2)} + a\sum {{x_i}} + b\sum {{y_i}} + cN = 0 $

(18)

式中, N为测量点总数。

将式(16)~式(18)三元一次方程组等价转换成矩阵方程为：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\sum {{X^2}} }&{\sum X Y}&{\sum X }\\ {\sum X Y}&{\sum {{Y^2}} }&{\sum Y }\\ {\sum X }&{\sum Y }&N \end{array}} \right] \times \left[ {\begin{array}{*{20}{l}} a\\ b\\ c \end{array}} \right] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {\sum {({X^2} + {Y^2})} X}\\ {\sum {({X^2} + {Y^2})} Y}\\ {\sum {({X^2} + {Y^2})} } \end{array}} \right] = 0 \end{array} $

(19)

令$\mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {\sum {{X^2}} }&{\sum {XY} }&{\sum X }\\ {\sum {XY} }&{\sum {{Y^2}} }&{\sum Y }\\ {\sum X }&{\sum Y }&N \end{array}} \right], \mathit{\boldsymbol{W = }}\left[ {\begin{array}{*{20}{c}} a\\ b\\ c \end{array}} \right], \mathit{\boldsymbol{T = - }}\left[ {\begin{array}{*{20}{c}} {\sum {\left( {{X^2} + {Y^2}} \right)X} }\\ {\sum {\left( {{X^2} + {Y^2}} \right)Y} }\\ {\sum {\left( {{X^2} + {Y^2}} \right)} } \end{array}} \right] $，则式(19)可变成为：

$ \mathit{\boldsymbol{K}} \times \mathit{\boldsymbol{W}} - \mathit{\boldsymbol{T}} = 0 $

(20)

左乘矩阵K的逆矩阵K^－1，并移位解算可得：

$ \mathit{\boldsymbol{W}} = {\mathit{\boldsymbol{K}}^{ - 1}} \times \mathit{\boldsymbol{T}} $

(21)

由式(4)、式(21)可计算得到圆的拟合圆心及拟合半径为：

$ \left\{ {\begin{array}{*{20}{c}} {{x_0} = - a/2}\\ {{y_0} = - b/2}\\ {R = \sqrt {{a^2} + {b^2} + 4c} /2} \end{array}} \right. $

(22)

1.2.3 迭代计算圆心

将步骤§1.2.2转换后的拟合圆心坐标迭代代入步骤§1.2.1和§1.2.2，直至圆心坐标偏移小于设定值ε，为：

$ \sqrt {{{(x_0^\prime - {x_0})}^2} + {{(y_0^\prime - {y_0})}^2}} \le \varepsilon $

(23)

式中，(x′ ₀, y′ ₀)为迭代过程中的临时圆心坐标。

1.2.4 计算各分段弧长误差

拟合圆圆上各点半径r_i为：

$ {r_i} = \sqrt {{{({x_i} - {x_0})}^2} + {{({y_i} - {y_0})}^2}} $

(24)

分段圆弧圆心角θ_i为：

$ {\theta _i} = {A_i} - {A_{i - 1}} $

(25)

则相邻点圆度或半径差引起的分段弧长差Δ_i为：

$ {\varDelta _i} = \pm ({r_i} - {r_{i - 1}}) \times {\theta _i} $

(26)

根据经验可设置Δ≤0.2 mm，如果超差，可根据需要通过内插测量点使相邻点的Δ减小。由于测量点数较多，内插点宜采用坐标放样法，先找出第i点及第(i-1)点，然后根据Δ_i的大小决定内插测量点数。待一圈内插测量点完成，须再次按之前步骤计算，最终得到拟合圆圆心和半径。

1.2.5 圆周长计算及其中误差

若各分段弧长误差满足式(26)的要求，且圆心满足式(23)要求，则圆周长C为：

$ C = \sum\limits_{i = 1}^N {{r_i}} \times {\theta _i} $

(27)

由于N次弧长测量为等精度测量，不难得出周长中误差为：

$ m_c^2 = N(\theta _0^2m_r^2 + r_0^2m_\theta ^2) $

(28)

式中，m_c为周长中误差；θ₀为单位角度固定误差，由测量点数引起；m_θ为角度中误差，由仪器测角误差引起；r₀为圆理论半径；m_r为半径中误差，由仪器测距误差引起。由中误差方程可知，要想提高圆周长测量精度可从减少点数、提高半径精度和提高圆弧圆心角精度3方面入手。从式(24)、式(25)可知，提高半径精度和圆弧圆心角精度都和圆心相关，所以圆心的解算精度直接关系到圆周长精度，这是本文讨论的重点之一。

以本文工程并结合§1.3.1所列仪器标称精度，计算圆周长中误差。各参数赋值为：N=144，${\theta _0} = \frac{{360 \times 3{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 600}}{{N \times \rho }} $，m_r=0.5 mm，r₀=15 000 mm，m_θ=0.5″/ρ，常数ρ=206 265，计算圆周长中误差m_c= ±0.51 mm。由此可见，本文所述方案的可行性较好。

1.2.6 金属热膨胀系数对圆周长的影响

金属有热胀冷缩的现象，金属热膨胀系数为线胀系数。线胀系数是指固态物质当温度改变1 ℃时，其某一方向上的长度的变化和它在20 ℃(即标准实验室环境)时的长度的比值^[3]。根据线胀系数定义不难得到该圆的热膨胀改正公式为：

$ {C_2} = {C_1}/[1 + a \times ({t_2} - {t_1})] $

(29)

式中，C₁为改正前长度在实际温度下的实测值；C₂为改正后长度在设计温度下的改正值；圆热膨胀系数(由机械设计师给出)a=10.7×10^－6/℃；t₂为钢板实测温度；t₁为设计温度(取20 ℃)。

热膨胀改正公式的3种方法如下。

1) 假定工件温度为20 ℃，实测测量点坐标，计算圆周长C₁，然后将C₁和实测温度t₂代入改正式(29)得到改正后的圆周长C₂。

2) 令C₁=1，由式(29)计算出C₂，然后在全站仪测量过程中对测量点坐标实时改正，为：

$ \left\{ {\begin{array}{*{20}{l}} {{x_{{\rm{Ne}}{{\rm{w}}_ - }i}} = {C_2} \times {x_{{\rm{Ol}}{{\rm{d}}_ - }i}}}\\ {{y_{{\rm{Ne}}{{\rm{w}}_ - }i}} = {C_2} \times {y_{{\rm{Ol}}{{\rm{d}}_ - }i}}} \end{array}} \right. $

(30)

3) 金属热膨胀是线性变化，所以可在测量完所有测量点后，对坐标系进行比例变化，为：

$ \left[ {\begin{array}{*{20}{c}} {{X_{{\rm{New}}}}}\\ {{Y_{{\rm{New}}}}} \end{array}} \right] = {C_2} \times \left[ {\begin{array}{*{20}{c}} {{X_{{\rm{ Old }}}}}\\ {{Y_{{\rm{ Old }}}}} \end{array}} \right] $

(31)

方法2)和3)坐标改正后，计算所得的圆周长即为热膨胀系数改正后的周长。

1.3 全站仪性能要求及其精度影响 1.3.1 全站仪性能要求

考虑到工厂车间工作环境较恶劣，兼顾测量精度与仪器经济性，本文选择高精度工业全站仪。由于测量点数较多，采用可双向通讯控制的马达驱动全站仪，从而实现高精度自动测量^[4]。

1.3.2 仪器精度对弧长精度的影响

在极坐标系下，弧长s表述如下：

$ s = r \times \theta $

(32)

式中, 圆心角θ的单位为弧度。

以本文工程实例半径15 m的圆为例，采用索佳NET05测量，测角测距误差取仪器标称精度。计算弧长误差如下：

1) 仪器测角误差引起的弧长误差。令r=15 000 mm，θ=0.5″，带入式(32)可解算出弧长误差Δs为：

$ \Delta s = 15{\kern 1pt} {\kern 1pt} {\kern 1pt} 000 \times (0.5 \times \frac{\pi }{{3{\kern 1pt} {\kern 1pt} {\kern 1pt} 600 \times 180}}) = 0.036{\kern 1pt} {\kern 1pt} {\rm{mm}} $

(33)

即该仪器在15 m远处由测角误差引起的弧长误差为0.036 mm。此种情况下，很明显测角误差几乎可忽略不计。

2) 仪器测距误差Δr引起的弧长误差关系为：

$ \Delta r \times \theta \le \Delta s $

(34)

令Δs=0.1 mm，Δr=0.5 mm，则可解出θ=11°18′32″，即只要圆心角小于11°18′32″，该仪器测距误差引起的弧长误差将小于0.1 mm。

由此可知，全站仪的测角精度要求须考虑待测异径圆的半径大小，要求测角精度与圆半径成反比；全站仪的测距精度决定了细部测量点的最少点数和两点间的最大圆心角。

2 程序设计

为方便现场作业及现场数据处理，软件运行平台为PDA，采用微软公司的Visual Studio 2008 C++语言编辑器基于Windows CE 5.0目标平台进行了程序设计^[5]，数值计算方面的矩阵求逆采用了全选主元高斯-约当消元法^[6]，该消元法程序实现较容易，代码可移植性较强。程序流程图如图 1所示。

程序特点：①可指定初始圆心点；②可指定观测方向，如顺时针或逆时针；③可设定分段弧长限差，若超差将在列表框中以红色等颜色提醒(弧长限差为0.1 mm)；④可指定计算圆弧长或圆周长；⑤给出角度检验，可粗略判断该圆是否闭合或圆弧圆心角是否正确；⑥给出迭代拟合圆信息，可根据工程需要或经验，人工判断是否继续或中断程序。

3 工程实例

该海上重型盆座式起重机的基座主要是圆筒结构分3层，且每层直径达30 m，分别由不同板厚的钢板(如图 3)焊接成3个分段圆筒，各分段圆筒再上下组合焊接成长圆筒(如图 2所示)，为保证各分段圆筒顺利合拢，各分段接口的圆周长要求精度较高，尤其是最下面的接口，需要与新加坡船厂对接，周长实测误差要求小于等于15 mm。

3.1 准备工作

为减小工件由于重力引起的变形，以下端口为基准用水准仪将工件粗略地测设水平(高低差小于30 mm)。为保证各测量点概略地在同一横截面上，基于下端口用记号笔标记了一整圈圆。然后在圆上概略等距地标记了144个测量点，并逐点粘贴10 cm×10 cm的反射片。

3.2 现场测量

采用自由设站法，在工件内任意位置架设全站仪并采集特征点，观测顺序为逆时针从点1到点72.5共计144点。同时根据图 3所示，记录每个点对应的钢板厚度改正值(壁厚补偿分别为30.0 mm、40.0 mm、45.0 mm)。圆上每隔90°测量一次钢板温度，取4次平均值作为钢板温度。可得到测量点原始坐标如表 1和各点壁厚补偿如表 2对应列，钢板温度为15.7 ℃。

为尽量减少接头钢板厚度变化引起的弧长误差，应在接头两侧相近位置多布两点，如图 3所示。如果仪器在圆内可采用内测法，半径补偿为正值；如果仪器在圆外可采用外测法，半径补偿为负值。极端恶劣的现场环境下，也可以内外混测，但须注意半径补偿的正负号。

3.3 数据处理

程序自动读取测量文件各点原始测量坐标和对应的各点的壁厚补偿，然后迭代拟合计算圆心坐标、各分段弧长误差。得到中间过程数据如表 2的改正后的各测量点数据(表 2在迭代计算时将动态更新到程序界面)，根据各分段弧长误差及实际情况，决定是否需要内插测量点。内插时可采用坐标放样的方法放样出相邻两点的实际点位，然后按图 1所示的步骤内插测量点并再次进行数据处理，直到各分段弧长误差满足要求。

3.4 生成报告

通过上一步数据处理，可以得到如表 2所示各点极径、极角，两两相邻点之间的弧长误差，所有弧长的累计误差为12.4 mm，表 3中异径圆的圆心坐标、半径均值，及最重要的圆周长，其圆周长误差可用弧长累计误差表示。将表 3的圆周长和钢板实测温度带入式(29)可得到该圆的改正后周长为94 221.8±12.4 mm(温度为20 ℃)，由于半径误差ΔR=12.4/2π≈2 mm，小于错边设计值5 mm，可以认为该测量误差不会对以后的对接造成影响。新加坡船厂采用该实测周长指导钢板下料与船上圆筒制作，取得了很高的制作精度。

另外，可以根据所有测量点的极径(半径)，根据半径值变化评估圆度误差^[7]，即(最大值-最小值)=圆度，如图 4所示，半径最大最小相差值(圆度)为19.3 mm，可根据极径和极角计算任意位置起始的弧长。为便于剔除圆周长粗差及比较数据，程序同时给出了各点相连的多线段总长。

3.5 算法收敛性

从表 3可以看出，采用本文算法即使圆心初始值误差非常大甚至任意值，也只用了3次迭代即收敛(程序会给出迭代过程中的相邻两次圆心误差)。

4 结束语

对工业测量中，不同厚度材料制作的大型同心异径圆圆周长的测量，本文系统地进行了探讨。指出了传统测量方法的不足，同时基于全站仪给出了一整套解决方案。该方案具有可操作性强，现场测量便捷，仪器可任意设在圆内、圆外或内外混测，算法易于程序实现。同时具有良好的可扩展性的特点，如只需将板厚改正修改成同一数值，就可测量同一厚度材料制作的常规圆周长及其圆心，可测量并解算出类似非规则的同心异径圆弧的圆心及其弧长，可按半径差法评价圆度等。对用全站仪内测法标定立式金属罐容积时^[8]，采用本法可测得较高精度的周长，再反算出较高精度的内周圆半径。