约束下考虑坐标分量误差相关性的直线拟合

宋占峰 ，郭捷佳 ，李 军

（中南大学土木工程学院，湖南 长沙 410075）

直线拟合不仅是曲线拟合研究的热点，并且在工程实践中被应用广泛[1-2]. 直线拟合是由n个测量点 (xi,yi),i=1,2,···,n，基于最小二乘准则找到一条最佳拟合直线.x、y分别为自变量、因变量，a、b分别为直线的斜率和截距，拟合直线的方程为y=ax+b.

普通最小二乘认为变量x或y无误差去估计参数a和b，该方法简单明了，但选择自变量和因变量不同，拟合的直线不同. 正交最小二乘认为x和y具有相同的精度，几何意义是测量点到拟合直线的垂直距离平方和最小，通常被认为是最佳拟合[3]，在铁路整正工程中被广泛采用[4-5].

为顾及自变量x存在误差，通常采用EIV （errorsin-variables）模型进行整体最小二乘估计直线参数[6-7].整体最小二乘通过建立变量的随机模型可以实现普通及正交最小二乘拟合直线，同时还可以实现加权整体最小二乘拟合直线[8]. 但此时还不能考虑自变量x和因变量y间的误差相关性. 当已知观测误差的随机特性时，参数的最优估计应符合这种随机特性[9]. Amiri-Simkooei 等[10]进一步提出了考虑观测值误差充分相关性的加权整体最小二乘法拟合直线.

直线是轨道线形的组成部分，为保证线路平顺性，需要估计直线参数，进行既有轨道线形整正. 采用全站仪采集的数据点，其纵、横坐标不仅精度不同，并且误差具有相关性. 因此，考虑观测值误差相关性的加权整体最小二乘才能处理这种情况，实现直线的最佳拟合. 整体最小二乘法是基于拟合直线得出的，要处理系数矩阵含有误差及与因变量y误差相关等情况，引入了EIV 模型及Kronecker 积等复杂矩阵运算[11]. 这在一定程度上造成了理解及工程应用上的困难[12]. 同时，线路中直线的拟合还受到相邻线元的约束.

因此，本文基于极大似然估计及拉格朗日条件极值原理建立直线拟合模型，引入表征测点位置的附加参数，推导出了顾及约束和观测值误差相关性直线拟合的通用方法. 实验验证了该方法能在任何误差分布情况下顾及约束估计直线参数及其精度.

1 纵、横坐标的误差相关性

由式（1）和误差传播定律可计算点P纵、横坐标的方差-协方差阵为

图1 全站仪测量原理Fig. 1 Measuring principle

图2 拟合原理Fig. 2 Fitting principle

2 约束下顾及误差相关性的直线拟合

2.1 直线拟合模型

直线上n个观测点坐标组成观测向量l，其联合密度函数为

线路中直线的拟合还要满足与相邻线元相切等约束条件，如已知相邻线元为圆曲线，其圆心位置（x0，y0）和半径R已知，则直线拟合受到的约束为

式中：B为系数矩阵，其秩为约束条件的个数；w为闭合差向量.

构造拉格朗日极值函数为

2.2 参数精度

求得参数的最佳估值，需进一步评定其精度. 依据协因数传播定律，由式（24）得到参数估值 Θˆ 的协因数阵为

2.3 迭代寻优过程

由于拟合模型舍去了高次项，参数初值为近似值. 因此，需要使用高斯-牛顿迭代寻优算法进行迭代计算，直至 δΘ 向量趋于0. 迭代终止条件设为‖δΘ‖≤ε ，ε为阈值，取为10−6. 算法流程见图3.

图3 高斯-牛顿迭代算法Fig. 3 Gauss-Newton iteration algorithm

3 试验算例

选取国内某既有专用线（设计速度30 km/h）直线段的复测数据，每20 m 测一个点，共20 个测点，坐标数据见表1. 使用的全站仪的测角精度为2″，测距精度为（2 + 2 × 10−6D） mm，D为距离，mm. 控制点A的坐标为（1000.012，1500.023），后视方位角为45°00′05″，由式（3）～（5）可计算出各测点纵、横坐标的方差及协方差，进而组成方差-协方差矩阵C.取先验单位权中误差 σ0=1 ，由式（10）可计算出实际测点的随机模型P1，其各项元素列于表1.表中：px、py和pxy分别为式（2）中c−1对应对角线元素和非对角线元素. 结果表明观测点的坐标分量不仅精度不同，并且误差具有相关性.

为验证直线的拟合方法的通用性，选用了对应正交最小二乘和普通最小二乘的两种随机模型对比计算，随机模型P2为单位矩阵，表示观测点的x坐标和y坐标具有相同的精度，即正交最小二乘拟合；模型P3表示x坐标的精度远高于y坐标，即普通最小二乘拟合.

由设计资料知该直线与圆曲线直接相接，圆心C的坐标为（1693.970，1443.300），半径R= 800 m. 为测试高斯-牛顿算法的效率和稳健性，选择的初始直线远离最优位置，初始直线参数列于表2 中的第0 次迭代，闭合差w为圆心到直线的距离减去半径，初始直线与圆曲线相离479.525 m. 算法经过6 次迭代收敛，w变为0，说明拟合直线满足约束，与圆曲线相切于直圆点（Z），图4（a）展示了该拟合过程，验证了方法的正确性.

表2 顾及约束和相关误差的直线拟合过程Tab. 2 Process of straight-line fitting with both a constraint and correlated noise

3 种随机模型的非零元素见表1. 约束下3 种随机模型的拟合直线参数和精度列于表3，算法均经过6 次迭代收敛，耗时均在0.6 s 以内. 基于P1得到的后验单位权中误差 σˆ0= 16.5 mm，截距的中误差σb= 60.0 mm，均为最小，相对于P2和P3得到的σb= 66.1 mm，精度提升了9.2%. 3 种随机模型下的直圆点坐标（xZ，yZ）列于表3，对应位置如图4（b），分别表示为Z1、Z2和Z3，其中Z2和Z3点重合.

表1 实地观测点坐标及采用的3 种随机模型Tab. 1 Coordinate pairs of field surveying data and three stochastic models for fitting

表3 约束下3 种随机模型拟合直线的参数估值及其精度Tab. 3 Parameter estimation of fitting line and their precisions of three stochastic models with constraints

令B为空矩阵，用同样的参数初值进行无约束的直线拟合，3 种随机模型的拟合直线参数和精度列于表4，算法均经过3 次迭代收敛，耗时均在0.5 s 以内，表明无约束拟合直线具有更高的效率.且无约束时的单位权中误差及参数精度均高于约束下的直线拟合，说明在无约束情况下，测点与拟合直线的贴合度更好. 但是，无约束拟合获得的直线不相切于已知圆曲线，如图4（b）所示，相离距离d= 0.18 m，不符合线路连续性的要求. 另外，在无约束拟合时，基于P1得到的后验单位权中误差σˆ0= 2.3 mm，截距的中误差 σb= 25.0 mm，相对P2和P3精度提升了2.7%. 说明不论是否有约束，考虑坐标相关误差时，获得的直线参数精度均为最高.

表4 无约束下3 种随机模型拟合直线的参数估值及其精度Tab. 4 Parameter estimation of line fitting and their precisions of three stochastic models without constraint

图4 约束及无约束的直线拟合Fig. 4 Straight-line fitting with a constraint and without constraint

4 结 论

在铁路维护中，全站仪获取的轨道坐标点具有误差相关性，普通最小二乘或正交最小二乘拟合直线不能考虑观测值之间的这种误差相关性. 提出的直线拟合通用模型，可以考虑坐标分量间的误差相关性实现直线拟合. 通过指定随机模型，可以实现普通最小二乘、正交最小二乘或加权整体最小二乘直线拟合，并揭示了其对应的几何意义.

通常的直线拟合方法未考虑约束条件，线路重构中直线的拟合受到相邻线元的约束. 提出的直线拟合通用模型，可以在约束下同时顾及误差相关性实现直线拟合. 顾及坐标相关误差时，可以提升估计参数的精度：约束及无约束下参数估计精度分别提高了9.2%和2.7%.

采用的高斯-牛顿算法运行效率高，能够快速得到直线参数的最佳估值及其精度. 在约束及无约束情况下分别仅6 次及3 次迭代就搜索出最优直线.

致谢：中南大学土木工程国家级实验教学示范中心项目（201905406）.