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

宋占峰; 郭捷佳; 李军

doi:10.3969/j.issn.0258-2724.20200120

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

doi: 10.3969/j.issn.0258-2724.20200120

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

基金项目: 国家自然科学基金（51678574）

详细信息

作者简介:
宋占峰（1973—），男，副教授，博士，研究方向为道路与铁道线路优化设计方法，E-mail：songzhanfeng@csu.edu.cn

通讯作者:
李军（1973—），男，副教授，博士，研究方向为道路与铁道工程信息化及优化，E-mail：lijun_csu@csu.edu.cn

中图分类号: U212.3
计量
- 文章访问数: 498
- HTML全文浏览量: 275
- PDF下载量: 23
- 被引次数: 1
出版历程
- 收稿日期: 2020-03-24
- 修回日期: 2020-08-03
- 网络出版日期: 2020-08-24
- 刊出日期: 2020-08-24

Fitting a Straight-Line to Data Points with Correlated Noise Between Coordinate Components under Constraints

School of Civil Engineering，Central South University，Changsha 410075，China

摘要

摘要:
直线拟合在曲线拟合研究及工程实践中受到广泛关注，常用的普通最小二乘和正交最小二乘忽略了坐标分量误差相关性的存在. 基于此，首先论证了在铁路线路整正中全站仪测量坐标点的纵横坐标间存在误差相关性，同时线路中直线的拟合受到相邻线元的约束；然后，基于极大似然估计及拉格朗日条件极值原理，推导出了顾及约束和坐标分量误差相关性的直线拟合通用模型，并给出了高斯-牛顿迭代算法搜索最优解；最后，采用了实测的数据进行了验证及测试. 试验结果表明：该方法能在任何误差分布情况下考虑约束估计直线参数及其精度；考虑坐标相关误差时，参数估计精度在约束及无约束下分别提高了9.2%和2.7%；高斯-牛顿算法在约束及无约束情况下分别仅6次及3次迭代就搜索出最优直线.
- 直线 /
- 曲线拟合 /
- 参数估计 /
- 误差相关 /
- 条件极值 /
- 算法
Abstract:
Straight-line fitting has received extensive attention both in curve fitting research and engineering practice. The methods of ordinary least squares and orthogonal least squares fitting ignore the existence of the observation error correlation. The coordinate pairs of surveying points, obtained by a total station in railway realignment, not only have different levels of precision but also have correlated noise. Meanwhile, straight-line fitting is usually under constraints in the realignment. Thus, a straight-line fitting model was derived based on the maximum likelihood estimation and Lagrange conditional extremum theory, considering constraints and correlated noise between coordinate components, and a Gauss-Newton algorithm was presented to search for the optimum. The method was tested with the field surveying data. Experimental results show that the proposed fitting method is capable of estimating straight-line parameters and their precisions in all circumstances by specifying stochastic models. When considering correlated noise, the precision of estimated parameters improve 9.2% with a constraint and improve 2.7% without constraints, respectively. The Gauss-Newton algorithm takes only 6 and 3 iteration times with a constraint and without constraints respectively, for locating the optimum straight-line.
- straight-line /
- curve fitting /
- parameter estimation /
- correlated noise /
- conditional extremum /
- algorithm

HTML全文

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

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

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

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

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

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

在工程实际中，经常使用全站仪获取点位坐标，测量原理见. 图中：XOY坐标系中，A $({x_A},{y_A})$ 和B $({x_B},{y_B})$ 为控制点，P(x,y)为全站仪测量点， ${\alpha _{AB}}$ 和 ${\alpha _{AP}}$ 分别为AB边和AP边的方位角. 测得独立观测量水平距离s和水平角β，则P点坐标为

图 1 全站仪测量原理

Figure 1. Measuring principle

下载: 全尺寸图片幻灯片

$\left\{ {\begin{array}{*{20}{l}} {x = {x_A} + s\cos ({\alpha _{AB}} + \beta ) = {x_A} + s\cos \;{\alpha _{AP}}{\text{，}}}\\ {y = {y_A} + s\sin ({\alpha _{AB}} + \beta ) = {y_A} + s\sin \;{\alpha _{AP}}.} \end{array}} \right.$

(1)

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

${\boldsymbol{c}} = \left[ {\begin{array}{*{20}{c}} {\sigma _x^2}&{{\sigma _{xy}}} \\ {{\sigma _{yx}}}&{\sigma _y^2} \end{array}} \right] = {\boldsymbol{K}}\left[ {\begin{array}{*{20}{c}} {\sigma _s^2}&0 \\ 0&{\sigma _\beta ^2} \end{array}} \right]{{\boldsymbol{K}}^{\rm T}},$

(2)

式中： $\sigma _x^{\rm{2}}$ 和 $\sigma _y^{\rm{2}}$ 分别为x、y的方差； ${\sigma _{xy}}$ 、 ${\sigma _{yx}}$ 分别为x和y、y和x间的协方差； ${\sigma _s}$ 、 ${\sigma _\beta }$ （单位：″ ）分别为全站仪测距和测角的先验精度；矩阵 ${\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {\cos \;{\alpha _{AP}}}&{ - \dfrac{s}{\rho }\sin \;{\alpha _{AP}}} \\ {\sin \;{\alpha _{AP}}}&{\dfrac{s}{\rho }\cos \;{\alpha _{AP}}} \end{array}} \right]$ .

将式（2）展开得

$\sigma ^2_x = \sigma _s^2{\cos ^2}{\alpha _{AP}} + \frac{{\sigma _\beta ^2{s^2}}}{{{\rho ^2}}}{\sin ^2}{\alpha _{AP}}{\text{，}}$

(3)

$\sigma _y^2 = \sigma _s^2{\sin ^2}{\alpha _{AP}} + \frac{{\sigma _\beta ^2{s^2}}}{{{\rho ^2}}}{\cos ^2}{\alpha _{AP}}{\text{，}}$

(4)

${\sigma _{xy}} = {\sigma _{yx}} = 0.5\left( {\sigma _s^2 - \frac{{\sigma _\beta ^2{s^2}}}{{{\rho ^2}}}} \right)\sin (2{\alpha _{AP}}),$

(5)

式（3）～（5）中： $\rho {\rm{ \; = 206\;265''}}$ .

点P坐标观测值l = [x y] 的期望在拟合直线上，对应投影位置为 $\boldsymbol{\mu }(\boldsymbol{\hat \varTheta }) = [\hat x\;\;\hat a\hat x + \hat b]$ ，其中：参数估值 $\boldsymbol{\hat \varTheta }= {[\hat a\; \hat b\; \hat x]^{\rm{T}}}$ ， $\hat a$ 和 $\hat b$ 分别为a和b的估值， $\hat x$ 为点P位置x的估值. 坐标改正向量 $\boldsymbol{r}(\boldsymbol{\varTheta }) = {[{r_x}\;{r_y}]^{\rm{T}}} = {\left( { \boldsymbol{l} - \boldsymbol{\mu }(\boldsymbol{\varTheta })} \right)^{\rm{T}}}$ ，其中： $\boldsymbol{\varTheta }$ 为参数向量，r_x和r_y分别为x和y的坐标改正数. 点P坐标观测值正态分布密度函数为

$f(\boldsymbol{l};\boldsymbol{\varTheta }) = {\left| {2{\text{π}} \boldsymbol{c}} \right|^{ - 1/2}}\exp \left( { - 0.5{\boldsymbol{r}^{\rm{T}}}(\boldsymbol{\varTheta }){\boldsymbol{c}^{ - 1}}\boldsymbol{r}(\boldsymbol{\varTheta })} \right).$

(6)

最小二乘原理为 ${\rm{min }}\{{p_x}r_x^2 + {p_y}r_y^2\}$ ，其中 ${p}_{x}$ 和 ${p}_{y}$ 分别为 ${r_x}$ 和 ${r_y}$ 的权重. 当坐标分量同精度时，最小二乘原理为 ${\rm{min }}\{r_x^2 + r_y^2\}$ . 此时，中测量点P的坐标改正向量 ${[{r_{{x_1}}}\;{r_{{y_1}}}]^{\rm{T}}}$ 垂直拟合直线，投影点为P₁，即正交最小二乘. 当纵坐标x无误差时，其权 ${p}_{x}$ 为无穷大，则 ${r_x}$ 必然为0. 此时，测量点P的坐标改正向量 ${[0\;{r_{{y_2}}}]^{\rm{T}}}$ 垂直X轴，在拟合直线上的投影点为P₂，即普通最小二乘.

图 2 拟合原理

Figure 2. Fitting principle

下载: 全尺寸图片幻灯片

对于全站仪测得的采样点，坐标分量不仅精度不同（ $\sigma _x^2 \ne \sigma _y^2$ ），并且误差具有相关性（ ${\sigma _{xy}} \ne 0$ ）. 此时，普通最小二乘或正交最小二乘就不再适合.

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

2.1 直线拟合模型

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

$L(\boldsymbol{l};\boldsymbol{\varTheta }) = \prod\limits_{i = 1}^n {f({\boldsymbol{l}_i};\boldsymbol{\varTheta })} {\text{，}}$

(7)

式中：l_i = [x_i y_i].

$\boldsymbol{\varTheta }$ 的极大似然估计量为

$\boldsymbol{\hat \varTheta }{\rm{ = }}\arg \max \;L(\boldsymbol{l};\boldsymbol{\varTheta }).$

(8)

由式（6）～（8）可推出 $\boldsymbol{\varTheta }$ 的极大似然估计等同于

$\min \left\{ {\sum\limits_{i = 1}^n {\boldsymbol{r}_i^{\rm{T}}(\boldsymbol{\varTheta }){\boldsymbol{c}_i}^{ - 1}{\boldsymbol{r}_i}(\boldsymbol{\varTheta })} } \right\}{\text{，}}$

(9)

式中： ${\boldsymbol{r}_i}(\boldsymbol{\varTheta }) = {[{r_{{x_i}}}\;{r_{{y_i}}}]^{\rm{T}}}$ ；c_i为测点i的方差-协方差矩阵.

令 ${\boldsymbol{r}^{\rm{T}}}\left( \boldsymbol{\varTheta }\right) = \left[ {\boldsymbol{r}_1^{\rm{T}}(\boldsymbol{\varTheta })\;\boldsymbol{r}_2^{\rm{T}}(\boldsymbol{\varTheta })\; \cdots\;\boldsymbol{r}_n^{\rm{T}}(\boldsymbol{\varTheta }]} \right)$ 及n个观测点的方差-协方差矩阵 $\boldsymbol{C} = {\rm{diag}} ({\boldsymbol{c}_1},{\boldsymbol{c}_2}, \cdots ,{\boldsymbol{c}_n})$ . 则 $\boldsymbol{C}$ 、协因数矩阵Q和权阵P即为随机模型，满足

$\boldsymbol{C}{\rm{ = }}\sigma _0^2\boldsymbol{Q} = \sigma _0^2{\boldsymbol{P}^{ - 1}}{\text{，}}$

(10)

式中： $\sigma _0^2$ 为先验单位权方差.

式（9）可写为

${\rm{min}}\left\{ {{\boldsymbol{r}^{\rm{T}}}\left( \boldsymbol{\varTheta }\right)\boldsymbol{Pr}\left( \boldsymbol{\varTheta }\right)} \right\}.$

(11)

即观测值正态分布时，极大似然估计和最小二乘估计等价. 易知，当P为单位矩阵I时，式（11）退化为 ${\rm{min}}\left\{ {{\boldsymbol{r}^{\rm{T}}}\left( \boldsymbol{\varTheta } \right)\boldsymbol{r}\left( \boldsymbol{\varTheta } \right)} \right\} = {\rm{min}}\left\{ {\displaystyle\sum\limits_{i = 1}^n {(r_{{x_i}}^2 + r_{{y_i}}^2} )} \right\}$ ，即正交最小二乘；当x权无穷大时，式（11）退化为 ${\rm{min}}\displaystyle\sum\limits_{i = 1}^n {r_{{y_i}}^2}$ ，即普通最小二乘.

第i个观测点的直线参数方程 ${\boldsymbol{\mu }_i}(\boldsymbol{\varTheta })$ 可表示为

$\left\{ {\begin{array}{l} {{{\hat x}_i} = {{\tilde x}_i} + {\rm{\delta }}{x_i},} \\ {{\hat y}_i} = (\tilde a + {\rm{\delta }}a) ({{\tilde x}_i} + {\rm{\delta }}{x_i}) + (\tilde b + {\rm{\delta }}b)\;{\rm{ =}}\\ \quad{{\tilde y}_i} + {{\tilde x}_i}{\rm{\delta }}a + \tilde a{\rm{\delta }}{x_i} + {\rm{\delta }}b + {\rm{\delta }}a{\rm{\delta }}b, \\ \end{array}} \right.$

(12)

式中： $(\tilde x,\tilde y)$ 和 $(\hat x,\hat y)$ 分别为点的近似坐标和估值坐标； $\tilde a$ 和 $\tilde b$ 分别为斜率和截距的近似值.

舍去高次微小量 ${\rm{\delta }}a{\rm{\delta }}b$ 后，式（12）线性化为

$\left\{ {\begin{array}{l} {{{\hat x}_i} = {{\tilde x}_i} + {\rm{\delta }}{x_i}{\text{，}}}\\ {{{\hat y}_i} = {{\tilde y}_i} + {{\tilde x}_i}{\rm{\delta }}a + \tilde a{\rm{\delta }}{x_i} + {\rm{\delta }}b.} \end{array}} \right.$

(13)

因此，线性化的直线方程可简化成矩阵形式如式（14）.

$\boldsymbol{\mu }({\boldsymbol{\hat \varTheta }}) = \boldsymbol{\mu }({\boldsymbol{\tilde \varTheta }}) + {\boldsymbol{A}}{\rm{\delta }}{\boldsymbol{\varTheta }}{\text{，}}$

(14)

式中： ${\boldsymbol{\tilde \varTheta }}$ 为 ${\boldsymbol{\varTheta}}$ 近似值； ${{\boldsymbol{A}}_{2n\times(n + 2)}} = \left[ {\begin{array}{*{20}{c}} {\bf{0}}&{\boldsymbol{I}} \\ {\boldsymbol{M}}&{\boldsymbol{N}} \end{array}} \right]$ ，其中 ${{\boldsymbol{M}}_{n \times 2}} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{{\tilde x}_1}}&{{{\tilde x}_{\rm{2}}}}& \cdots &{{{\tilde x}_n}} \end{array}} \\ {\begin{array}{*{20}{c}} 1&1& \cdots &1 \end{array}} \end{array}} \right]^{\rm{T}}}$ ， ${{\boldsymbol{N}}_{n \times n}} = {\rm{diag}}(\tilde a)$ .

由式（11）、（14），目标函数可写为

$\min \left\{ {{{\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right]}^{\rm{T}}}\boldsymbol{P}\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right]} \right\}.$

(15)

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

$\frac{{{y_0} - a{x_0} - b}}{{\sqrt {1 + {a^2}} }} - R = 0.$

(16)

线性化后为

$\begin{split} &{\frac{{ - {x_0}(1 + {{\tilde a}^2}) - \left( {{y_0} - \tilde a{x_0} - \tilde b} \right)\tilde a}}{{1{\rm{ + }}{{\tilde a}^2}}}{\rm{\delta }}a - }\\ &\quad{{\rm{\delta }}b + {y_0} - \tilde a{x_0} - \tilde b - R(1 + {{\tilde a}^2}) = 0}. \end{split}$

(17)

约束条件线性化后，矩阵表示为

$\boldsymbol{B}{\rm{\delta }}\boldsymbol{\varTheta } + \boldsymbol{w} = \bf{0}{\text{，}}$

(18)

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

构造拉格朗日极值函数为

$\begin{split} &\min \left\{ {{\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right]}^{\rm{T}}}\boldsymbol{P}\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right] +\right.\\ &\quad\bigg.{\boldsymbol{\gamma }^{\rm{T}}}\left( {\boldsymbol{B}{\rm{\delta }}\boldsymbol{\varTheta } + \boldsymbol{w}} \right) \bigg\}, \end{split}$

(19)

式中： $\boldsymbol{\gamma }$ 为拉格朗日乘向量.

为得到极小值，对式（19）求一阶导数并令其为0，可得

${\boldsymbol{A}^{\rm{T}}}\boldsymbol{PA}{\rm{\delta }}\boldsymbol{\varTheta } - {\boldsymbol{A}^{\rm{T}}}\boldsymbol{Pr}\left( {\boldsymbol{\tilde \varTheta }} \right) + {\boldsymbol{B}^{\rm{T}}}\boldsymbol{\gamma } = \bf{0}\bf{.}$

(20)

令 $\boldsymbol{N} = {\boldsymbol{A}^{\rm{T}}}\boldsymbol{PA}$ 及 $\boldsymbol{M} = {\boldsymbol{A}^{\rm{T}}}\boldsymbol{Pr}\left( {\boldsymbol{\tilde \varTheta }} \right)$ ，式（20）写为

$\boldsymbol{N}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{M} + {\boldsymbol{B}^{\rm{T}}}\boldsymbol{\gamma } = \bf{0}\bf{.}$

(21)

式（21）乘以 $\boldsymbol{B}{\boldsymbol{N}^{ - 1}}$ ，再减去式（18）可得

$\boldsymbol{B}{\boldsymbol{N}^{ - 1}}{\boldsymbol{B}^{\rm{T}}}\boldsymbol{\gamma } - \boldsymbol{B}{\boldsymbol{N}^{ - 1}}\boldsymbol{M} - \boldsymbol{w} = \bf{0}\bf{.}$

(22)

令 ${\boldsymbol{N}_{\rm{c}}} = \boldsymbol{B}{\boldsymbol{N}^{ - 1}}{\boldsymbol{B}^{\rm{T}}}$ ，由式（22）可求得

$\boldsymbol{\gamma } = \boldsymbol{N}_{\rm{c}}^{ - 1}\left( {\boldsymbol{B}{\boldsymbol{N}^{ - 1}}\boldsymbol{M}{\rm{ + }}\boldsymbol{w}} \right).$

(23)

将式（23）代入式（21）得

$\begin{split} &{\rm{\delta }}\boldsymbol{\varTheta } = {\boldsymbol{N}^{ - 1}}\left( {\boldsymbol{I} - {\boldsymbol{B}^{\rm{T}}}\boldsymbol{N}_{\rm{c}}^{ - 1}\boldsymbol{B}{\boldsymbol{N}^{ - 1}}} \right)\boldsymbol{M} \; - \\ &\quad{\boldsymbol{N}^{ - 1}}{\boldsymbol{B}^{\rm{T}}}\boldsymbol{N}_{\rm{c}}^{ - 1}\boldsymbol{w}. \end{split}$

(24)

当无约束时，B为空矩阵，则式（24）变为

${\rm{\delta }}\boldsymbol{\varTheta } = {\boldsymbol{N}^{ - 1}}\boldsymbol{M}.$

(25)

约束条件下参数的估值为

${\boldsymbol{\hat \varTheta }} = {\boldsymbol{\tilde \varTheta }} + {\rm{\delta }}{\boldsymbol{\varTheta }}.$

(26)

2.2 参数精度

求得参数的最佳估值，需进一步评定其精度. 依据协因数传播定律，由式（24）得到参数估值 ${\boldsymbol{\hat \varTheta }}$ 的协因数阵为

${\boldsymbol{Q}_{\boldsymbol{\hat \varTheta }}} = {\boldsymbol{N}^{ - 1}}(\boldsymbol{I} - {\boldsymbol{B}^{\rm{T}}}\boldsymbol{N}_{\rm{c}}^{ - 1}\boldsymbol{B}{\boldsymbol{N}^{ - 1}}).$

(27)

同理，当无约束时，式（27）变为

${\boldsymbol{Q}_{\boldsymbol{\hat \varTheta }}} = {\boldsymbol{N}^{ - 1}}{\rm{ = }}{({\boldsymbol{A}^{\rm{T}}}\boldsymbol{PA})^{ - 1}}.$

(28)

因此，附加约束后，参数的协因数将减小.

由式（10）可得 ${\boldsymbol{\hat \varTheta }}$ 的方差-协方差如式（29）所示.

${\boldsymbol{C}_{\boldsymbol{\hat \varTheta }}} = \hat \sigma _0^2{\boldsymbol{Q}_{\boldsymbol{\hat \varTheta }}}{\text{，}}$

(29)

$\hat \sigma _0^2 = \frac{{\boldsymbol{r}{{(\boldsymbol{\hat \varTheta })}^{\rm{T}}}\boldsymbol{Pr}(\boldsymbol{\hat \varTheta })}}{{2n - t + m}},$

(30)

式中： ${\hat \sigma _0}$ 为后验单位权中误差，反映了观测值与模型之间的拟合程度；t为 ${\boldsymbol{\hat \varTheta }}$ 中元素个数；m为约束个数.

2.3 迭代寻优过程

由于拟合模型舍去了高次项，参数初值为近似值. 因此，需要使用高斯-牛顿迭代寻优算法进行迭代计算，直至 ${\rm{\delta }}{\boldsymbol{\varTheta }}$ 向量趋于0. 迭代终止条件设为 $\left\| {{\rm{\delta }}{\boldsymbol{\varTheta }}} \right\| \leqslant \varepsilon$ ， $\varepsilon$ 为阈值，取为10⁻⁶. 算法流程见图3.

图 3 高斯-牛顿迭代算法

Figure 3. Gauss-Newton iteration algorithm

下载: 全尺寸图片幻灯片

3. 试验算例

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

表 1 实地观测点坐标及采用的3种随机模型

Table 1. Coordinate pairs of field surveying data and three stochastic models for fitting

点号	x/m	y/m	C 中非零元素			P₁中非零元素			P₂中非零元素	P₃中非零元素
点号	x/m	y/m	$\sigma _x^2$ /mm²	$\sigma _y^2$ /mm²	${\sigma _{xy}}$ /mm²	p_x	p_y	p_xy	${p_x}/{p_y}$	${p_x}/{p_y}$
1	688.639	1398.869	7.3371	9.7881	−0.8902	0.1378	0.1033	0.0125	1	10000
2	701.467	1383.525	7.3289	9.3014	−0.9080	0.1381	0.1088	0.0135	1	10000
3	714.294	1368.180	7.3340	8.8888	−0.9115	0.1381	0.1140	0.0142	1	10000
4	727.121	1352.835	7.3547	8.5485	−0.9080	0.1378	0.1185	0.0146	1	10000
5	739.953	1337.495	7.3948	8.2767	−0.9044	0.1371	0.1225	0.0150	1	10000
6	752.783	1322.152	7.4602	8.0684	−0.9071	0.1359	0.1257	0.0153	1	10000
7	765.609	1306.806	7.5584	7.9166	−0.9210	0.1342	0.1281	0.0156	1	10000
8	778.434	1291.460	7.6977	7.8126	−0.9487	0.1319	0.1299	0.0160	1	10000
9	791.262	1276.115	7.8867	7.7477	−0.9909	0.1289	0.1312	0.0165	1	10000
10	804.088	1260.770	8.1340	7.7132	−1.0462	0.1251	0.1320	0.0170	1	10000
11	816.915	1245.425	8.4471	7.7012	−1.1111	0.1207	0.1324	0.0174	1	10000
12	829.740	1230.078	8.8323	7.7054	−1.1805	0.1156	0.1325	0.0177	1	10000
13	842.564	1214.731	9.2939	7.7206	−1.2485	0.1100	0.1324	0.0178	1	10000
14	855.389	1199.384	9.8346	7.7437	−1.3087	0.1040	0.1321	0.0176	1	10000
15	868.213	1184.036	10.4556	7.7729	−1.3546	0.0979	0.1316	0.0171	1	10000
16	881.043	1168.694	11.1562	7.8077	−1.3803	0.0916	0.1309	0.0162	1	10000
17	893.875	1153.353	11.9355	7.8489	−1.3806	0.0855	0.1301	0.0150	1	10000
18	906.703	1138.009	12.7917	7.8981	−1.3511	0.0796	0.1289	0.0136	1	10000
19	919.537	1122.670	13.7218	7.9569	−1.2880	0.0740	0.1276	0.0120	1	10000
20	932.371	1107.331	14.7236	8.0277	−1.1886	0.0687	0.1261	0.0102	1	10000

下载: 导出CSV

| 显示表格

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

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

表 2 顾及约束和相关误差的直线拟合过程

Table 2. Process of straight-line fitting with both a constraint and correlated noise

迭代数/次	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	w/mm
0	−0.681654264	2210.153480	59093.318880	6.14777 × 10⁻²	93044.287650	−479524.603400
1	−0.935165106	2013.498778	5844.096808	6.07990 × 10⁻³	9201.714082	−59428.755230
2	−1.146887043	2183.687057	448.743795	9.66425 × 10⁻⁴	1148.244207	−9789.990571
3	−1.195130784	2221.779805	10.362946	3.19907 × 10⁻⁵	35.137584	−394.045005
4	−1.197201392	2223.403509	16.492003	5.52543 × 10⁻⁵	59.714483	−0.691388
5	−1.197204886	2223.406205	16.525002	5.55595 × 10⁻⁵	60.003229	−0.000002
6	−1.197204886	2223.406205	16.525002	5.55598 × 10⁻⁵	60.003515	0

下载: 导出CSV

| 显示表格

图 4 约束及无约束的直线拟合

Figure 4. Straight-line fitting with a constraint and without constraint

下载: 全尺寸图片幻灯片

3种随机模型的非零元素见表1. 约束下3种随机模型的拟合直线参数和精度列于，算法均经过6次迭代收敛，耗时均在0.6 s以内. 基于P₁得到的后验单位权中误差 ${\hat \sigma _0}$ = 16.5 mm，截距的中误差 ${\sigma _b}$ = 60.0 mm，均为最小，相对于P₂和P₃得到的 ${\sigma _b}$ = 66.1 mm，精度提升了9.2%. 3种随机模型下的直圆点坐标（x_Z，y_Z）列于表3，对应位置如图4（b），分别表示为Z₁、Z₂和Z₃，其中Z₂和Z₃点重合.

表 3 约束下3种随机模型拟合直线的参数估值及其精度

Table 3. Parameter estimation of fitting line and their precisions of three stochastic models with constraints

随机模型	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	x_Z/m	y_Z/m	迭代数/次	耗时/s
P₁	−1.19720489	2223.4062	16.525	5.55598 × 10⁻⁵	60.0	1079.9809	930.4478	6	0.494
P₂	−1.19722236	2223.4251	49.031	6.12012 × 10⁻⁵	66.1	1079.9772	930.4522	6	0.503
P₃	−1.19722233	2223.4250	76.478	6.12012 × 10⁻⁵	66.1	1079.9772	930.4522	6	0.496

下载: 导出CSV

| 显示表格

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

表 4 无约束下3种随机模型拟合直线的参数估值及其精度

Table 4. Parameter estimation of line fitting and their precisions of three stochastic models without constraint

随机模型	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	迭代数/次	耗时/s
P₁	−1.196269322	2222.672000	2.293655	3.10233 × 10⁻⁵	25.008852	3	0.461
P₂	−1.196259074	2222.663879	6.710739	3.16316 × 10⁻⁵	25.743943	3	0.414
P₃	−1.196259065	2222.663872	0.462488	3.16316 × 10⁻⁵	25.743942	3	0.438

下载: 导出CSV

| 显示表格

4. 结　论

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

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

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

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

图 1 全站仪测量原理

Figure 1. Measuring principle

下载: 全尺寸图片幻灯片

图 2 拟合原理

Figure 2. Fitting principle

下载: 全尺寸图片幻灯片

图 3 高斯-牛顿迭代算法

Figure 3. Gauss-Newton iteration algorithm

下载: 全尺寸图片幻灯片

图 4 约束及无约束的直线拟合

Figure 4. Straight-line fitting with a constraint and without constraint

下载: 全尺寸图片幻灯片

表 1 实地观测点坐标及采用的3种随机模型

Table 1. Coordinate pairs of field surveying data and three stochastic models for fitting

点号	x/m	y/m	C 中非零元素			P₁中非零元素			P₂中非零元素	P₃中非零元素
点号	x/m	y/m	$\sigma _x^2$ /mm²	$\sigma _y^2$ /mm²	${\sigma _{xy}}$ /mm²	p_x	p_y	p_xy	${p_x}/{p_y}$	${p_x}/{p_y}$
1	688.639	1398.869	7.3371	9.7881	−0.8902	0.1378	0.1033	0.0125	1	10000
2	701.467	1383.525	7.3289	9.3014	−0.9080	0.1381	0.1088	0.0135	1	10000
3	714.294	1368.180	7.3340	8.8888	−0.9115	0.1381	0.1140	0.0142	1	10000
4	727.121	1352.835	7.3547	8.5485	−0.9080	0.1378	0.1185	0.0146	1	10000
5	739.953	1337.495	7.3948	8.2767	−0.9044	0.1371	0.1225	0.0150	1	10000
6	752.783	1322.152	7.4602	8.0684	−0.9071	0.1359	0.1257	0.0153	1	10000
7	765.609	1306.806	7.5584	7.9166	−0.9210	0.1342	0.1281	0.0156	1	10000
8	778.434	1291.460	7.6977	7.8126	−0.9487	0.1319	0.1299	0.0160	1	10000
9	791.262	1276.115	7.8867	7.7477	−0.9909	0.1289	0.1312	0.0165	1	10000
10	804.088	1260.770	8.1340	7.7132	−1.0462	0.1251	0.1320	0.0170	1	10000
11	816.915	1245.425	8.4471	7.7012	−1.1111	0.1207	0.1324	0.0174	1	10000
12	829.740	1230.078	8.8323	7.7054	−1.1805	0.1156	0.1325	0.0177	1	10000
13	842.564	1214.731	9.2939	7.7206	−1.2485	0.1100	0.1324	0.0178	1	10000
14	855.389	1199.384	9.8346	7.7437	−1.3087	0.1040	0.1321	0.0176	1	10000
15	868.213	1184.036	10.4556	7.7729	−1.3546	0.0979	0.1316	0.0171	1	10000
16	881.043	1168.694	11.1562	7.8077	−1.3803	0.0916	0.1309	0.0162	1	10000
17	893.875	1153.353	11.9355	7.8489	−1.3806	0.0855	0.1301	0.0150	1	10000
18	906.703	1138.009	12.7917	7.8981	−1.3511	0.0796	0.1289	0.0136	1	10000
19	919.537	1122.670	13.7218	7.9569	−1.2880	0.0740	0.1276	0.0120	1	10000
20	932.371	1107.331	14.7236	8.0277	−1.1886	0.0687	0.1261	0.0102	1	10000

下载: 导出CSV

表 2 顾及约束和相关误差的直线拟合过程

Table 2. Process of straight-line fitting with both a constraint and correlated noise

迭代数/次	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	w/mm
0	−0.681654264	2210.153480	59093.318880	6.14777 × 10⁻²	93044.287650	−479524.603400
1	−0.935165106	2013.498778	5844.096808	6.07990 × 10⁻³	9201.714082	−59428.755230
2	−1.146887043	2183.687057	448.743795	9.66425 × 10⁻⁴	1148.244207	−9789.990571
3	−1.195130784	2221.779805	10.362946	3.19907 × 10⁻⁵	35.137584	−394.045005
4	−1.197201392	2223.403509	16.492003	5.52543 × 10⁻⁵	59.714483	−0.691388
5	−1.197204886	2223.406205	16.525002	5.55595 × 10⁻⁵	60.003229	−0.000002
6	−1.197204886	2223.406205	16.525002	5.55598 × 10⁻⁵	60.003515	0

下载: 导出CSV

表 3 约束下3种随机模型拟合直线的参数估值及其精度

Table 3. Parameter estimation of fitting line and their precisions of three stochastic models with constraints

随机模型	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	x_Z/m	y_Z/m	迭代数/次	耗时/s
P₁	−1.19720489	2223.4062	16.525	5.55598 × 10⁻⁵	60.0	1079.9809	930.4478	6	0.494
P₂	−1.19722236	2223.4251	49.031	6.12012 × 10⁻⁵	66.1	1079.9772	930.4522	6	0.503
P₃	−1.19722233	2223.4250	76.478	6.12012 × 10⁻⁵	66.1	1079.9772	930.4522	6	0.496

下载: 导出CSV

表 4 无约束下3种随机模型拟合直线的参数估值及其精度

Table 4. Parameter estimation of line fitting and their precisions of three stochastic models without constraint

随机模型	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	迭代数/次	耗时/s
P₁	−1.196269322	2222.672000	2.293655	3.10233 × 10⁻⁵	25.008852	3	0.461
P₂	−1.196259074	2222.663879	6.710739	3.16316 × 10⁻⁵	25.743943	3	0.414
P₃	−1.196259065	2222.663872	0.462488	3.16316 × 10⁻⁵	25.743942	3	0.438

下载: 导出CSV

参考文献(12)

[1]	KRYSTEK M, ANTON M. A least-squares algorithm for fitting data points with mutually correlated coordinates to a straight line[J]. Measurement Science and Technology, 2011, 22(3): 035101.1-035101.9.
[2]	PETROLINI A. Linear least squares fit when both variables are affected by equal uncorrelated errors[J]. American Journal of Physics, 2014, 82(12): 1178-1185.
[3]	丁克良,沈云中,欧吉坤. 整体最小二乘法直线拟合[J]. 辽宁工程技术大学学报(自然科学版),2010,29(1): 44-47. DING Keliang, SHEN Yunzhong, OU Jikun. Methods of line-fitting based on total least-squares[J]. Journal of Liaoning Technical University (Natural Science), 2010, 29(1): 44-47.
[4]	宋占峰,彭欣,吴清华. 基于中线坐标的地铁调线优化算法[J]. 西南交通大学学报,2014,49(4): 656-661. doi: 10.3969/j.issn.0258-2724.2014.04.015 SONG Zhanfeng, PENG Xin, WU Qinghua. Optimization algorithm for horizontal realignment based on coordinate of metro centerline[J]. Journal of Southwest Jiaotong University, 2014, 49(4): 656-661. doi: 10.3969/j.issn.0258-2724.2014.04.015
[5]	宋占峰,王健,李军. 缓和曲线正交拟合的Levenberg-Marquardt算法[J]. 西南交通大学学报,2020,55(1): 144-149. doi: 10.3969/j.issn.0258-2724.20190130 SONG Zhanfeng, WANG Jian, LI Jun. Levenberg-Marquardt algorithm for orthogonal fitting of transition curves[J]. Journal of Southwest Jiaotong University, 2020, 55(1): 144-149. doi: 10.3969/j.issn.0258-2724.20190130
[6]	KARL P. On lines and planes of closest fit to systems of points in space[J]. Philosophical Magazine, 1901, 2(11): 559-572.
[7]	刘经南,曾文宪,徐培亮. 整体最小二乘估计的研究进展[J]. 武汉大学学报(信息科学版),2013,38(5): 505-512. LIU Jingnan, ZENG Wenxian, XU Peiliang. Overview of total least squares methods[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 505-512.
[8]	KRYSTEK M, ANTON M. A weighted total least-squares algorithm for fitting a straight line[J]. Measurement Science and Technology, 2007, 18(11): 3438-3442. doi: 10.1088/0957-0233/18/11/025
[9]	SONG Z F, DING H, LI J, et al. Circular curve-fitting method for field surveying data with correlated noise[J]. Journal of Surveying Engineering, 2018, 144(4): 04018010.1-04018010.9.
[10]	AMIRI-SIMKOOEI A R, ZANGENEH-NEJAD F, ASGARI J, et al. Estimation of straight line parameters with fully correlated coordinates[J]. Measurement, 2014, 48: 378-386. doi: 10.1016/j.measurement.2013.11.005
[11]	SHEN Y F, LI B F, CHEN Y. An iterative solution of weighted total least-squares adjustment[J]. Journal of Geodesy, 2011, 85(4): 229-238. doi: 10.1007/s00190-010-0431-1
[12]	鲁铁定,陶本藻,周世健. 基于整体最小二乘法的线性回归建模和解法[J]. 武汉大学学报(信息科学版),2008,33(5): 504-507. LU Tieding, TAO Benzao, ZHOU Shijian. Modeling and algorithm of linear regression based on total least squares[J]. Geomatics and Information Science of Wuhan University, 2008, 33(5): 504-507.

施引文献

期刊类型引用(0)
其他类型引用(1)

附加材料(0)

访问统计

点击查看大图

图(4) / 表(4)

计量

文章访问数: 498
HTML全文浏览量: 275
PDF下载量: 23
被引次数: 1

1. 纵、横坐标的误差相关性
2. 约束下顾及误差相关性的直线拟合
2.1 直线拟合模型
2.2 参数精度
2.3 迭代寻优过程
3. 试验算例
4. 结　论

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

doi: 10.3969/j.issn.0258-2724.20200120

作者简介: 宋占峰（1973—），男，副教授，博士，研究方向为道路与铁道线路优化设计方法，E-mail：songzhanfeng@csu.edu.cn

通讯作者: 李军（1973—），男，副教授，博士，研究方向为道路与铁道工程信息化及优化，E-mail：lijun_csu@csu.edu.cn

计量

出版历程

Fitting a Straight-Line to Data Points with Correlated Noise Between Coordinate Components under Constraints

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

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

2.1 直线拟合模型

2.2 参数精度

2.3 迭代寻优过程

3. 试验算例

4. 结 论

期刊类型引用(0)

其他类型引用(1)

计量

出版历程

目录

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

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

2.1 直线拟合模型

2.2 参数精度

2.3 迭代寻优过程

3. 试验算例

4. 结 论

作者简介:
宋占峰（1973—），男，副教授，博士，研究方向为道路与铁道线路优化设计方法，E-mail：songzhanfeng@csu.edu.cn

通讯作者:
李军（1973—），男，副教授，博士，研究方向为道路与铁道工程信息化及优化，E-mail：lijun_csu@csu.edu.cn

4. 结　论

4. 结　论