针对圆曲线拟合问题，许多学者进行了各种算法研究。文献[1]以圆曲线一般方程为基础，采用最小二乘原理求解出圆曲线参数；文献[2]同样以圆曲线一般方程为基础，以原半径与拟合半径的几何距离差值平方和最小为准则，使用泰勒公式展开成最小二乘形式迭代求解；文献[3]以圆曲线正交距离残差平方和最小为拟合准则，基于参数方程提出采用正交最小二乘法拟合圆曲线；文献[4]采用总体最小二乘(total least square，TLS)算法求解待定参数；文献[5]提出利用协方差传播定律定权的加权总体最小二乘(weight total least square, WTLS)算法解算圆拟合参数的思想，但缺少实例；文献[6]在参数方程的拟合模型中采用总体最小二乘算法求解拟合参数，但平差计算量较大；文献[7]以圆曲线参数方程为基础构建拟合模型，采用最小二乘(least square，LS)算法取得较好的效果，但中间参数初始值的确定较为复杂。

针对总体最小二乘算法，文献[8]根据系数矩阵既含非随机元素又含重复的随机元素的特点提出Partial EIV(partial errors-in-variables)模型，该模型将系数矩阵中非随机元素与重复的随机元素分离并进行平差计算，使平差模型更具一般性；文献[9]对Partial EIV模型进行移项变化，转换为间接平差形式，采用两步迭代法求解，简化了参数的表达形式且提高了计算效率；文献[10]概述了整体最小二乘估计的发展历史及研究进展，侧重分析了各种算法的本质特征。

本文首先分析现有的圆曲线模型的构建及其参数求解方法，然后以圆曲线的参数方程为基础，针对模型的系数矩阵特点，建立圆曲线拟合的Partial EIV模型，采用两步迭代法求解模型参数。结合实例数据和模拟数据，对比分析现有的以一般方程构建的圆曲线拟合模型对应的最小二乘算法、EIV模型的总体最小二乘算法。实验结果表明，本文提出的基于Partial EIV模型的总体最小二乘法求解由参数方程构建的拟合模型是可行的。

1 拟合模型构建 1.1 圆曲线一般方程形式

圆曲线一般方程形式为^[11]：

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

(1)

式中，(a, b)为圆心，半径$ r = \sqrt {{a^2} + {b^2} - c} $。

1.2 圆曲线参数方程的EIV模型

平面上确定圆心(a, b)和半径r，就可以唯一确定圆。假设圆上一点m(x, y), 采用圆的参数表达形式为：

$ \left\{ \begin{array}{l} x = a + r\cos \theta \\ y = b + r\sin \theta \end{array} \right. $

(2)

式中，θ∈[0 2π]是旋转角参数。式(2)采用平差值表示为：

$ \left\{ \begin{array}{l} \hat x = \hat a + \hat r\cos \hat \theta \\ \hat y = \hat b + \hat r\sin \hat \theta \end{array} \right. $

(3)

式中，$\left( {\hat a, \hat b} \right)$为圆心坐标平差值，$ {\hat r}$、${\hat \theta }$为半径和旋转角的平差值。令${\hat x}$=x+v_x, ${\hat y}$=y+v_y, ${\hat r}$=r⁰+δr, ${\hat \theta }$=θ⁰+δθ, ${\hat a}$ =a⁰+δa, ${\hat b}$ =b⁰+δb，将式(3)线性化后误差方程为：

$ \left\{ \begin{array}{l} {v_x} = \delta a + \cos {\theta ^0}\delta r - {r^0}\sin {\theta ^0}\frac{{\delta \theta }}{\rho } - {L_x}\\ {v_y} = \delta b + \sin {\theta ^0}\delta r + {r^0}\cos {\theta ^0}\frac{{\delta \theta }}{\rho } - {L_y} \end{array} \right. $

(4)

式中，$\left\{ \begin{array}{l} {L_x} = x - ({a^0} + {r^0}{\rm{cos}}{\theta ^0}) = x - {x^0}\\ {L_y} = y - ({b^0} + {r^0}{\rm{sin}}{\theta ^0}) = y - {y^0} \end{array} \right. $，θ⁰、x⁰、y⁰计算如下：

$ \left\{ \begin{array}{l} \cos {\theta ^0} = \frac{{x - a}}{{\sqrt {{{\left( {x - a} \right)}^2} + {{\left( {y - b} \right)}^2}} }}\\ \sin {\theta ^0} = \frac{{y - b}}{{\sqrt {{{\left( {x - a} \right)}^2} + {{\left( {y - b} \right)}^2}} }}\\ {x^0} = a + {r^0}\cos {\theta ^0}\\ {y^0} = b + {r^0}\sin {\theta ^0} \end{array} \right. $

(5)

式(4)写成矩阵形式为：

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{v_x}}\\ {{v_y}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0&{\cos {\theta ^0}}&{ - {r^0}\sin {\theta ^0}/\rho }\\ 0&1&{\sin {\theta ^0}}&{{r^0}\cos {\theta ^0}/\rho } \end{array}} \right] \cdot }\\ {\left[ {\begin{array}{*{20}{c}} {\delta a}\\ {\delta b}\\ {\delta r}\\ {\delta \theta } \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {x - {x^0}}\\ {y - {y^0}} \end{array}} \right]} \end{array} $

(6)

$ {\rm{令}}\;\;\;\;\;\;\;\mathit{\boldsymbol{l}} = \left[ {\begin{array}{*{20}{c}} {x - {x^0}}\\ {y - {y^0}} \end{array}} \right],\;\;\;\;{\mathit{\boldsymbol{e}}_l} = \left[ {\begin{array}{*{20}{c}} {{v_x}}\\ {{v_y}} \end{array}} \right] $

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 1&0&{\cos {\theta ^0}}&{ - {r^0}\sin {\theta ^0}/\rho }\\ 0&1&{\sin {\theta ^0}}&{{r^0}\cos {\theta ^0}/\rho } \end{array}} \right],\mathit{\boldsymbol{\xi }} = \left[ {\begin{array}{*{20}{c}} {\delta a}\\ {\delta b}\\ {\delta r}\\ {\delta \theta } \end{array}} \right] $

则式(6)简化为：

$ {\mathit{\boldsymbol{e}}_l} = \mathit{\boldsymbol{A\xi }} - \mathit{\boldsymbol{l}} $

(7)

式中,l为观测向量，e_l为对应的观测向量误差；A为系数矩阵；ξ为参数改正数；ρ= $\left( {\frac{{180}}{\pi }} \right)$×60×60；系数矩阵 A中元素求解时需采用带有误差的观测值(x, y)，故A含有随机误差，由此可以构建圆曲线拟合的EIV模型。

假设有圆曲线的一组数据为p_i(x_i, y_i)，需要由这组点拟合出圆曲线的参数，相应的误差方程可以根据式(5)表示出来。其中,

$ \begin{array}{l} \mathit{\boldsymbol{l = }}\left[ {\begin{array}{*{20}{c}} {{x_1} - x_1^0}\\ {{y_1} - y_1^0}\\ \vdots \\ {{x_i} - x_i^0}\\ {{y_i} - y_i^0} \end{array}} \right] \in {R^{2n \times 1}},\mathit{\boldsymbol{A}} = \\ \left[ {\begin{array}{*{20}{c}} 1&0&{\cos \theta _1^0}&{ - {r^0}\sin \theta _1^0/\rho }&0& \cdots &0\\ 0&1&{\sin \theta _1^0}&{ - {r^0}\cos \theta _1^0/\rho }&0& \cdots &0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1&0&{\cos \theta _n^0}&0& \cdots &0&{ - {r^0}\sin \theta _i^0/\rho }\\ 0&1&{\sin \theta _n^0}&0& \cdots &0&{ - {r^0}\cos \theta _i^0/\rho } \end{array}} \right] \in \end{array} $

R^2n×(n+3)，ξ =[δa δb δc δθ₁ … δθ_n]^T∈R^(n+3)×1。采用总体最小二乘解出 ξ，即确定了圆曲线方程。

2 基于Partial EIV模型的参数估计

分析圆曲线参数方程拟合模型式(7)可发现，系数矩阵 A既含有随机元素又含有重复的非随机元素。为了保证重复元素的改正数一致，且减少平差计算量，采用基于Partial EIV模型的总体最小二乘求解拟合参数。Partial EIV模型是将系数矩阵中的随机元素提取出来，所对应的数学模型如下^[8]。

1) 函数模型：

$\mathit{\boldsymbol{l}}\mathit{ = }\left( {{\mathit{\boldsymbol{\xi }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{2n}}} \right)\left( {\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar t}}} \right) - {\mathit{\boldsymbol{e}}_l}\\ \mathit{\boldsymbol{t}} = \mathit{\boldsymbol{\bar t}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{t}}} $

(8)

2) 随机模型：

$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{l}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{t}}}} \end{array}} \right] \sim {N_{6n}}\left( {\left[ \begin{array}{l} 0\\ 0 \end{array} \right],\sigma _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{l}}}}&0\\ 0&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{t}}}} \end{array}} \right]} \right) $

(9)

式中，l 为2n×1的观测向量；h 是由系数矩阵 A按列拉直后的零元素和常数元素所组成的2n(n+3)×1的列向量；Q_l、Q_t分别是观测向量 y和列向量 a的协因数阵；B是构造出的2n(n+3)×4n的已知矩阵；I_2n为2n×2n的单位矩阵；σ₀²为单位权方差；⊗为Kronecker积；t 是由 A中不同随机元素所组成的4n×1向量，t表示 t 的真值，e_t为t 的随机误差。圆拟合的参数方程对应的向量t =[cosθ₁⁰ sinθ₁⁰ cosθ₂⁰ sinθ₂⁰ … cosθ_n⁰ sinθ_n⁰ -r⁰sinθ₁⁰/ρ r⁰cosθ₁⁰/ρ … -r⁰ sinθ_n⁰/ρ r⁰cosθ_n⁰/ρ]^T，采用Partial EIV模型进行参数计算时，所对应的已知向量 h 构造如下：

$ \mathit{\boldsymbol{h}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{h}}_1}}&{{\mathit{\boldsymbol{h}}_2}}&{{\mathit{\boldsymbol{h}}_3}} \end{array}} \right]^{\rm{T}}} $

(10)

式中，h₁为n个[1 0]^T组成的列向量，h₂为n个[0 1]^T组成的列向量，h₃为2n(n+1)×1的零向量。

以式(8)为基础，将第一个式子展开并移项^[9]为：

$ \left. \begin{array}{l} \mathit{\boldsymbol{l}} - \left( {{\mathit{\boldsymbol{\xi }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{h}} = \left( {{\mathit{\boldsymbol{\xi }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{B\bar t}} - {\mathit{\boldsymbol{e}}_l}\\ \mathit{\boldsymbol{t}} = {\mathit{\boldsymbol{I}}_{n + 3}}\mathit{\boldsymbol{\bar t}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{t}}} \end{array} \right\} $

(11)

令 $\mathit{\boldsymbol{L = }}\left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{l}}-\left( {{\mathit{\boldsymbol{\xi }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{2n}}} \right)\mathit{\boldsymbol{h}}}\\ \mathit{\boldsymbol{t}} \end{array}} \right]$, $\mathit{\boldsymbol{C = }}\left[{\begin{array}{*{20}{c}} {\left( {{\mathit{\boldsymbol{\xi }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{2n}}} \right)\mathit{\boldsymbol{B}}}\\ {{\mathit{\boldsymbol{I}}_{n + 3}}} \end{array}} \right]$，$\mathit{\boldsymbol{ \boldsymbol{\varDelta} = }}\left[{\begin{array}{*{20}{c}} {-{\mathit{\boldsymbol{e}}_l}}\\ {{\mathit{\boldsymbol{e}}_t}} \end{array}} \right]$，$\mathit{\boldsymbol{P}} = \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}_l^{-1}}&0\\ 0&{\mathit{\boldsymbol{Q}}_t^{-1}} \end{array}} \right]$，则式(11)等同于:

$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{C\bar t}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} $

(12)

式中，C 为新构建的系数矩阵，并把 t作为参数平差。

令 Δ的估值为V，并在准则V^TPV =min下按间接平差原理解算平差参数 t的估值${\mathit{\boldsymbol{\hat t}}}$：

$ \mathit{\boldsymbol{\hat t}} = {\left( {{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{PC}}} \right)^{ - 1}}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{PL}} $

(13)

进而由${\mathit{\boldsymbol{\hat t}}}$重新构造成新的系数矩阵${\mathit{\boldsymbol{\hat A}}}$，再次按间接平差原理求解模型参数${\mathit{\boldsymbol{\hat \xi}}}$：

$ \mathit{\boldsymbol{\hat \xi }} = {\left( {{{\mathit{\boldsymbol{\hat A}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_l}\mathit{\boldsymbol{\hat A}}} \right)^{ - 1}}{{\mathit{\boldsymbol{\hat A}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_l}\mathit{\boldsymbol{l}} $

(14)

最后将计算出的${\mathit{\boldsymbol{\hat \xi }}}$前3项加上初始值(${{{\hat a}^0}, {{\hat b}^0}}$)、r⁰得到最终的圆心坐标和半径。

基于Partial EIV模型的总体最小二乘方法求解参数方程的具体计算步骤如下。

1) 初始值:采用一般方程形式的最小二乘解作为圆心初始值(${{{\hat a}^0}, {{\hat b}^0}}$)和半径的初始值r⁰，由(${{{\hat a}^0}, {{\hat b}^0}}$)、r⁰根据式(5)、(6)、(7)得到参数方程的系数矩阵初始值 A⁽⁰⁾和观测向量初始值 l⁽⁰⁾。再用式(7)求得的${\mathit{\boldsymbol{\hat \xi }}}$作为拟合参数的初值${{\mathit{\boldsymbol{\hat \xi }}}^{\left( 0 \right)}}$，并令${{\mathit{\boldsymbol{\hat t}}}^{\left( 0 \right)}} = \mathit{\boldsymbol{t}}$

2) 由 ${{\mathit{\boldsymbol{\hat \xi }}}^{\left( i \right)}}$、${{\mathit{\boldsymbol{\hat t }}}^{\left( i \right)}}$构造新的观测向量L⁽ⁱ⁾和系数矩阵C⁽ⁱ⁾。

3) 将步骤2)计算的 L⁽ⁱ⁾、C⁽ⁱ⁾代入式(13)计算${{\mathit{\boldsymbol{\hat t}}}^{\left( {i + 1} \right)}}$。

4) 把上一步求得的${{\mathit{\boldsymbol{\hat t}}}^{\left( {i + 1} \right)}}$重新构造成改正后的系数矩阵${{\mathit{\boldsymbol{\hat A}}}^{\left( {i + 1} \right)}}$。

5) 根据式(14)的迭代公式计算新的参数估值${{\mathit{\boldsymbol{\hat \xi}}}^{\left( {i + 1} \right)}}$。

6) 重复步骤2)~5)，直到前后两次计算的未知参数估值${{\mathit{\boldsymbol{\hat \xi}}}^{\left( {i } \right)}}$满足‖${{\mathit{\boldsymbol{\hat \xi}}}^{\left( {i + 1} \right)}}$-${{\mathit{\boldsymbol{\hat \xi}}}^{\left( {i } \right)}}$‖＜ε且‖${{\mathit{\boldsymbol{\hat t}}}^{\left( {i + 1} \right)}}$-${{\mathit{\boldsymbol{\hat t}}}^{\left( {i } \right)}}$‖＜δ时，迭代终止，输出所求参数的最终估值。

7) 将输出的${{\mathit{\boldsymbol{\hat \xi}}}^{\left( {i } \right)}}$前3个参数加上初始值$\left( {{{\hat a}^0}, {{\hat b}^0}} \right)$、r⁰，得到最终的圆心坐标和半径。

一般总体最小二乘系数矩阵误差改正量为2n(n+3)，而本文Partial EIV模型的总体最小二乘系数矩阵误差改正量为4n，较前者减少了2n²+2n个误差改正量的计算。特别是当拟合点个数n较多时，模型的待估量大大减少。

3 圆曲线拟合的精度评定

拟合点与圆心之间的几何距离(或称为拟合前的原半径)与拟合后的半径之差d_i可以作为圆拟合精度评定的一个指标^{[2, 11]}。d_i的计算公式为：

$ d_i^2 = {\sum\limits_{i = 1}^n {\left[ {\sqrt {{{\left( {{x_i} - a} \right)}^2} + {{\left( {{y_i} - b} \right)}^2}} - r} \right]} ^2} $

(15)

对于模拟的圆曲线，可以把已知圆曲线的方程表示成一般方程形式或式(2)的参数表达形式。通过比较模拟点解算的参数平差值与真实值分别计算的残差值，可以评价拟合效果。

4 算例 4.1 实例数据

采用文献[7]提供的一组实例数据，坐标值如表 1，分别采用以一般方程构建的高斯-马尔科夫模型的LS算法，以参数方程构建的EIV模型、Partial EIV模型的TLS算法3种算法计算，结果如表 2所示。

从表 2可以看出，本文算法求解的拟合参数与EIV模型的TLS解基本一致，且两种算法的距离残差和均为1.227 599 1，这说明本文算法在圆拟合中的有效性。这主要是因为Partial EIV模型是EIV模型更为一般的表达形式，涵盖了一般总体最小二乘算法需要特殊处理的各种情况，例如部分元素为非随机量导致的 Q 奇异，A呈现结构性特征的结构总体最小二乘等。本文算法的系数矩阵改正量个数为24，较EIV模型的108个减少77.78%，可以看出本文算法大幅减少了平差计算量。

4.2 模拟数据

为进一步验证算法的可行性，采用仿真数据实验。假定拟合的圆方程为：

$ {\left( {x + 1} \right)^2} + {\left( {y + 2} \right)^2} = {3^2} $

(16)

利用Matlab在圆上取6个点，坐标如表 3所示。在x、y中分别加入中误差为0.01的随机误差进行圆拟合，同样采用上面3种算法求解拟合参数。模拟100次，取参数计算的平均值，结果如表 4所示。‖(Δa, Δb, Δr)‖₂为参数值与真值残差值的二范数。

从表 4可以看出，以参数方程建立的拟合模型的TLS解较一般方程的LS解更接近真值，说明顾及系数矩阵误差的TLS算法拟合效果好于LS算法，而以参数方程构建的Partial EIV拟合模型总体最小二乘解与EIV模型解基本一致，且均接近真值，这说明本文解法的有效性。本文方法的拟合距离残差和为0.242 239×10^-3，参数值与真值残差值的二范数为2.799 848×10^-3，在3种方法中均为最小，这说明本文方法在圆曲线拟合中有一定的优势。这是因为Partial EIV模型的总体最小二乘算法能保证系数矩阵相同元素的改正数一致，常数元素的改正数为零，较EIV模型的总体最小二乘算法更为合理严谨。同算例1，本算例中系数矩阵改正量数为24个，较EIV模型的108个减少77.78%。

5 结语

针对圆曲线拟合参数求解问题，提出一种基于Partial EIV模型的总体最小二乘算法。通过模拟数据和实例数据对算法的正确性、解算精度、适用性进行了验证，得到以下结论：

1) 不同形式构建的圆曲线拟合模型均可用于求解拟合参数。以一般方程构建的最小二乘法求解简单，但拟合精度较低。

2) 以圆曲线参数形式为基础构建的拟合模型，当观测坐标有误差导致系数矩阵有误差时，需采用总体最小二乘算法求解。

3) 本文给出的基于Partial EIV模型的总体最小二乘算法适用于圆曲线拟合，为圆拟合提供了一种新方法，且在拟合精度上具有一定优势。

4) 本文方法较EIV模型的总体最小二乘法，具有系数矩阵改正量少的优点。