在平面直线y=ax+b拟合中，当测量点(x_i, y_i)均含有误差时，通常采用最小二乘平差、整体最小二乘平差。最小二乘平差忽略x方向误差，建立误差方程求解直线参数a、b，该方法简单便捷但理论不严密。整体最小二乘法将x方向误差纳入设计矩阵中，根据整体最小平差准则求解未知参数的最佳估值，理论严密，符合实际应用。

目前，对于空间三维直线拟合，大致分为最小二乘法^[1-3]、牛顿-梯度最优化法^[4]、整体最小二乘法等方法^[5-10]。最小二乘法，采用降维思想^{[1, 2]}忽略部分方向上的误差，将三维直线问题转换成二维平面直线拟合，最后还原到三维空间直线。郭际明等^[1]通过求解与测量点距离和最小的平面，并将空间点投影到该平面上，在该平面拟合一条平面直线，最后经过相应旋转变换得到目标空间直线。王伟锋等^[2]将空间直线投影到空间3个互相垂直的平面上，从而得到3条平面直线，分别拟合这3条直线，并分别过直线做垂直于当前平面的垂面，最后3个垂面相交即为所求的目标直线。该方法在建模时，没有同时考虑3个方向的误差，理论不严密且不符合实际^{[11, 12]}。

牛顿-梯度最优化法^{[4, 10]}以${\varepsilon ^2} = \sum\limits_{i = 1}^n {(\varepsilon _{i1}^2 + \varepsilon _{i2}^2\varepsilon _{i3}^2)} $最小的准则(其中ε_i1、ε_i2、ε_i3为x、y、z 3个方向的误差)，以两点为初始值，给定某一个精度，逐步迭代求取最佳的方向余弦(α β γ)^T。该方法在建模时同时考虑到3个方向误差，理论严密且符合实际。

整体最小二乘法^[6]将空间直线标准式转换成设计矩阵含有误差的模型，将z坐标纳入系数矩阵中，采用整体最小二乘迭代法，求解未知参数估值。但该方法未考虑观测值不等权，观测值有粗差等诸多情况。文献[7-9]提出了基于稳健总体最小二乘原理的三维空间拟合方法，即混合结构总体最小二乘算法的直线拟合算法，极大地丰富了整体最小二乘空间直线拟合解法。

最小二乘方法^[1]无法同时考虑到3方向的误差，而牛顿-梯度最优化法^[4]、整体最小二乘法^[6]涉及到大量数学公式推导，在常规工程测量拟合中适用性不高。本文提出了同时顾及3方向坐标误差的最小二乘解法，其原理是基于极坐标法建模将观测点的x、y、z 3个方向分离出来，建立最小二乘模型，解算出空间直线的天顶距、方位角和极径的最佳估值。该方法同时兼顾观测值的不等精度问题，且模型简单直观。

1 空间直线拟合极坐标方法

该空间直线拟合方法，基于极坐标方法，将观测值x、y、z分离开。在极坐标中，确定一空间点位坐标，需要一个方位角A、一个天顶距Z、一个极径ρ，如图 1所示。在拟合空间直线中，假设该基准点在该直线上，那么拟合空间直线的必要观测数与观测值数的关系为：当点数为N时，观测值数n=3N，必要观测数t=2+N。

在拟合空间直线中需要事先确定一个基准点，经实验分析该基准点的点位误差直接影响到拟合直线的精度，因此, 基准点最好是空间已知点数据，尽量避免选择精度不高的观测点作为直线的基准点。那么在缺失已知基准点情况下，如何确定空间直线？本文将该基准点的点位坐标纳入未知参数中，其基准点的近似值选取可以采用以下两种方法：①将所有观测值点位坐标取平均得到中心点，采用此方法在列观测方程时，中心点两边的观测点会存在方位角、天顶距相差180°的问题，可以通过设置极径正负区分中心点两边的坐标，从而解决基准点两侧观测值的方位角、天顶距相差180°问题；②选取该直线的端点观测值为近似值。为了方便起见，本文采用第②种方法。将基准点设置为未知参数进行参数估计，其未知参数个数为3+2+N(N为点个数; 2为空间直线的天顶距及方位角; 3为基准点的点位坐标个数)。其中, 含有2+N个独立参数，3个相关参数。由于未知参数具有相关性，导致法方程秩亏，本文采用伪逆解法求未知参数。

观测方程为：

$ \left\{ \begin{array}{l} {{\tilde x}_i} = {{\tilde x}_0} + {{\tilde \rho }_{0i}} \times \sin \tilde Z \times \cos \tilde A\\ {{\tilde y}_i} = {{\tilde y}_0} + {{\tilde \rho }_{0i}} \times \sin \tilde Z \times \sin \tilde A\\ {{\tilde z}_i} = {{\tilde z}_0} + {{\tilde \rho }_{0i}} \times \cos \tilde Z \end{array} \right. $

(1)

将观测方程进行线性化为：

$ \left\{ \begin{array}{l} {V_{{x_i}}} + {x_i} = {x_0} + {\rho _{0i}}\sin Z\cos A + {\rm{d}}{x_0} + \\ \sin Z\cos A{\rm{d}}{\rho _{0i}} + {\rho _{0i}}\cos A\cos Z \times \frac{{{\rm{d}}Z}}{{\rho ''}} - \\ \sin A\sin Z \times {\rho _{0i}}\frac{{{\rm{d}}A}}{{\rho ''}}\\ {V_{{y_i}}} + {y_i} = {y_0} + {\rho _{0i}}\sin Z\sin A + {\rm{d}}{y_0} + \\ \sin Z\sin A \times {\rm{d}}{\rho _{0i}} + {\rho _{0i}}\sin A\cos Z \times \frac{{{\rm{d}}Z}}{{\rho ''}} + \\ \cos A\sin Z \times {\rho _{0i}}\frac{{{\rm{d}}A}}{{\rho ''}}\\ {V_{{z_i}}} + {z_i} = {z_0} + {\rho _{0i}}\cos Z + {\rm{d}}{z_0} + \\ \cos Z \times {\rm{d}}{\rho _{0i}} - {\rho _{0i}}\sin Z \times \frac{{{\rm{d}}Z}}{{\rho ''}} \end{array} \right. $

(2)

建立误差方程：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{CX}} - \mathit{\boldsymbol{l}} $

(3)

式中，

$ {\mathit{\boldsymbol{V}}_{\left( {3N,1} \right)}} = {\left[ {\begin{array}{*{20}{c}} {{x_1}}\\ {{y_1}}\\ {{z_1}}\\ \vdots \\ {{x_n}}\\ {{y_n}}\\ {{z_n}} \end{array}} \right]_{\left( {3N,1} \right)}};{\mathit{\boldsymbol{x}}_{\left( {3 + 2 + N,1} \right)}} = {\left[ {\begin{array}{*{20}{c}} {{{\hat x}_0}}\\ {{{\hat y}_0}}\\ {{{\hat z}_0}}\\ {\hat Z}\\ {\hat A}\\ {{\rho _0}}\\ \vdots \\ {{\rho _N}} \end{array}} \right]_{\left( {3 + 2 + N,1} \right)}};{\mathit{\boldsymbol{l}}_{3N,1}} = - {\left[ {\begin{array}{*{20}{c}} {{x_0} + {\rho _{0i}}\sin Z\cos A - {x_1}}\\ {{y_0} + {\rho _{0i}}\sin Z\sin A - {y_1}}\\ {{z_0} + {\rho _{0i}}\cos Z - {z_1}}\\ \vdots \\ {{x_0} + {\rho _{0i}}\sin Z\cos A - {x_N}}\\ {{y_0} + {\rho _{0i}}\sin Z\sin A - {y_N}}\\ {{z_0} + {\rho _{0i}}\cos Z - {z_N}} \end{array}} \right]_{\left( {3N,1} \right)}}; $

$ \mathit{\boldsymbol{B}} = {\left[ {\begin{array}{*{20}{c}} 1&0&0&{\frac{{{\rho _{01}}\cos A\cos Z}}{{\rho ''}}}&{ - \frac{{\sin A\sin Z{\rho _{01}}}}{{\rho ''}}}&{\sin Z\cos A}&0& \cdots \\ 0&1&0&{\frac{{{\rho _{01}}\sin A\cos Z}}{{\rho ''}}}&{\frac{{\cos A\sin Z{\rho _{01}}}}{{\rho ''}}}&{\sin Z\sin A}&0&{}\\ 0&0&1&{ - \frac{{{\rho _{01}}\sin Z}}{{\rho ''}}}&0&{\cos Z}&0& \cdots \\ \vdots &{}&{}&{}&{}&{}&{}& \vdots \\ 1&0&0&{\frac{{{\rho _{0N}}\cos A\cos Z}}{{\rho ''}}}&{ - \frac{{\sin A\sin Z \times {\rho _{0N}}}}{{\rho ''}}}&0& \cdots &{\sin Z\cos A}\\ 0&1&0&{\frac{{{\rho _{0N}}\sin A\cos Z}}{{\rho ''}}}&{\frac{{\cos A\sin Z \times {\rho _{0N}}}}{{\rho ''}}}&0&{}&{\sin Z\sin A}\\ 0&0&1&{ - \frac{{{\rho _{0N}}\sin Z}}{{\rho ''}}}&0&0& \cdots &{\cos Z} \end{array}} \right]_{\left( {3N,3 + 2 + N} \right)}} $

根据实际观测值的精度来定权。如在实际应用中，用GPS获得某条直线的三维坐标，将GPS的网平差结果及其方差作为直线拟合的输入。根据最小二乘准则，有：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{N}} = {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PB}}\\ {\rm{rank}}\left( \mathit{\boldsymbol{N}} \right) = 2 + \mathit{\boldsymbol{N}}\\ {\mathit{\boldsymbol{N}}^ + } = {\rm{pinv}}\left( \mathit{\boldsymbol{N}} \right)\\ \hat x = {\mathit{\boldsymbol{N}}^ + }{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_l}\\ X = {X_0} + \hat x \end{array} \right. $

(4)

式中，N⁺为法方程N的伪逆。带入式(3)可解出V，则有：

$ \hat \sigma = \sqrt {\frac{{{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{PV}}}}{{n - {\rm{rank}}\left( \mathit{\boldsymbol{B}} \right)}}} $

(5)

2 空间直线精度评定方法

对于空间直线拟合精度评定，通常采用空间各点p_i(x_i, y_i, z_i)到拟合直线的距离之和以及直线度进行评定^[4]。空间直线的一般表达形式为：

$ \left\{ \begin{array}{l} {A_1}x + {B_1}y + {C_1}z + {D_1} = 0\\ {A_2}x + {B_2}y + {C_2}z + {D_2} = 0 \end{array} \right. $

(6)

那么，空间一点到该直线的距离为：

$ \left\{ \begin{array}{l} {\Delta _i} = \frac{{\left| {\left( {{A_1}{x_i} + {B_1}{y_i} + {C_1}{z_i} + {D_1}} \right){n_2} - \left( {{A_2}{x_i} + {B_2}{y_i} + {C_2}{z_i} + {D_2}} \right){n_1}} \right|}}{{\left| {{n_1} \times {n_2}} \right|}}\\ 直线度 = {\Delta _{\max }} - {\Delta _{\min }} \end{array} \right. $

(7)

式中，Δ_i表示空间一点到直线的距离；n_i=(A_i B_i C_i)。

3 实例分析

此次实验数据采用文献[4]的表 6.1空间直线仿真数据，没有物理意义，因此无单位，其模型误差为δ=±0.005。编写相应的Matlab程序，基准点的近似值选取空间直线样本数据的端点坐标(3.003 6, 2.996 0, 3.004 1)，其解算的出来的空间点位坐标的平差值以及点到直线距离Δ_i如表 1所示。其中，解算出的方位角最佳估值为63.446 195 5°，天顶距最佳估值为36.702 931 8°，基准点的最佳估值为(3.003 2, 2.998 1, 3.002 8)，单位权中误差为0.002 9，空间直线的表达式为：

$ \frac{{x - 3.003\;2}}{1} = \frac{{y = 2.998\;1}}{{2.000\;981\;94}} = \frac{{z - 3.002\;8}}{{3.000\;767\;34}} $

(8)

整体最小二乘法^[6]、最小二乘法(无迭代法)^[1]、牛顿-梯度最优法^[5]的解算结果与本文方法计算的空间直线方向矢量比较如表 2所示。

通过表 2可以看出，本文极坐标方法解算出来的直线方向矢量比最小二乘方法(无迭代法)解算结果更接近整体最小二乘的结果。原因是极坐标法跟整体最小二乘同时考虑到x、y、z误差，两者准则等价，因此，解算结果较为接近。

将确定的空间直线按照式(5)求解各观测点到拟合的空间直线的距离。并将本文方法的统计结果与整体最小二乘法^[6]、最小二乘法(无迭代法)^[1]、牛顿-梯度最优法^[5]结果比较，如表 3所示。通过直线度和误差总和ΣΔ_i可以看出，本文方法在空间直线拟合问题中具有可行性。关于整体最小二乘、无迭代法、极坐标法的残差和一致，邓小渊等^[8]指出距离残差和一致，因为它们都是根据Σ(Δx_i²+Δy_i²+Δz_i²)最小的平差准则或者是等价的平差准则下采用不同方法拟合的直线。极坐标法、整体最小二乘原理与其他最小二乘方法的区别是同时考虑x、y、z 3个方向误差，理论上更为严密。

4 平面直线、圆、球面等拟合问题

基于本文提出的空间直线拟合极坐标法的原理，可以通过扩展来解决诸多工程拟合问题，比如平面直线、圆^[11]、球面等。对于平面直线而言，可以将平面直线上的一点拆分成方位角A、极径ρ；平面圆拆分成方位角A、极径ρ；空间球面拆分成天顶距Z、方位角A、极径ρ。其极坐标原理，如图 2所示。

5 结束语

1) 基于极坐标的空间直线拟合算法克服了以往最小二乘在拟合空间直线无法同时考虑3个方向误差。通过建模将观测值的x、y、z 3个方向分解开，避免了系数矩阵的误差，模型简单直观。通过实例分析比较，该方法具有可行性，并且在精度与以前方法精度相当。

2) 本文提出的兼顾的不等权的加权极坐标法解决了以往文献不能考虑权的差异性问题，更贴近工程实际需要，适用性更高。

3) 该方法拆分观测值思路，同样适用于平面直线、平面圆、空间球面等工程测量拟合问题。且该方法变化灵活，可进一步加入抗差模型，提高其可靠性。