
Citation: | LI Cun-jun, YANG Ru-gui, DENGHong-xia. Combination Prediction of Traffic-Volume in Intersections Based on Wavelet Analysis and Kalman Filter[J]. Journal of Southwest Jiaotong University, 2004, 17(4): 577-580. |
直线拟合不仅是曲线拟合研究的热点,并且在工程实践中被应用广泛[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点坐标为
{x=xA+scos(αAB+β)=xA+scosαAP,y=yA+ssin(αAB+β)=yA+ssinαAP. | (1) |
由式(1)和误差传播定律可计算点P纵、横坐标的方差-协方差阵为
c=[σ2xσxyσyxσ2y]=K[σ2s00σ2β]KT, | (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[ {cosαAP−sρsinαAPsinαAPsρcosαAP} \right]$.
将式(2)展开得
σ2x=σ2scos2αAP+σ2βs2ρ2sin2αAP, | (3) |
σ2y=σ2ssin2αAP+σ2βs2ρ2cos2αAP, | (4) |
σxy=σyx=0.5(σ2s−σ2βs2ρ2)sin(2α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 }$为参数向量,rx和ry分别为x和y的坐标改正数. 点P坐标观测值正态分布密度函数为
f(l;Θ)=|2πc|−1/2exp(−0.5rT(Θ)c−1r(Θ)). | (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\}$. 此时,图2中测量点P的坐标改正向量 ${[{r_{{x_1}}}\;{r_{{y_1}}}]^{\rm{T}}}$垂直拟合直线,投影点为P1,即正交最小二乘. 当纵坐标x无误差时,其权$ {p}_{x} $为无穷大,则 ${r_x}$ 必然为0. 此时,测量点P的坐标改正向量${[0\;{r_{{y_2}}}]^{\rm{T}}}$ 垂直X轴,在拟合直线上的投影点为P2,即普通最小二乘.
对于全站仪测得的采样点,坐标分量不仅精度不同($\sigma _x^2 \ne \sigma _y^2$),并且误差具有相关性(${\sigma _{xy}} \ne 0$). 此时,普通最小二乘或正交最小二乘就不再适合.
直线上n个观测点坐标组成观测向量l,其联合密度函数为
L(l;Θ)=n∏i=1f(li;Θ), | (7) |
式中:li = [xi yi].
$\boldsymbol{\varTheta }$的极大似然估计量为
ˆΘ=argmax | (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}}}$;ci为测点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) |
线路中直线的拟合还要满足与相邻线元相切等约束条件,如已知相邻线元为圆曲线,其圆心位置(x0,y0)和半径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) |
求得参数的最佳估值,需进一步评定其精度. 依据协因数传播定律,由式(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为约束个数.
由于拟合模型舍去了高次项,参数初值为近似值. 因此,需要使用高斯-牛顿迭代寻优算法进行迭代计算,直至${\rm{\delta }}{\boldsymbol{\varTheta }}$向量趋于0. 迭代终止条件设为$\left\| {{\rm{\delta }}{\boldsymbol{\varTheta }}} \right\| \leqslant \varepsilon $,$\varepsilon $为阈值,取为10−6. 算法流程见图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. 取先验单位权中误差 ${\sigma _0}{\rm{ = }}1$,由式(10)可计算出实际测点的随机模型P1,其各项元素列于表1. 表中:px、py和pxy分别为式(2)中c−1对应对角线元素和非对角线元素. 结果表明观测点的坐标分量不仅精度不同,并且误差具有相关性.
点号 | x/m | y/m | C 中非零元素 | P1 中非零元素 | P2 中非零元素 | P3 中非零元素 | |||||||
$\sigma _x^2$/mm2 | $\sigma _y^2$/mm2 | ${\sigma _{xy}}$/mm2 | px | py | pxy | ${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 |
为验证直线的拟合方法的通用性,选用了对应正交最小二乘和普通最小二乘的两种随机模型对比计算,随机模型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)展示了该拟合过程,验证了方法的正确性.
迭代数/次 | $\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−2 | 93044.287650 | −479524.603400 |
1 | −0.935165106 | 2013.498778 | 5844.096808 | 6.07990 × 10−3 | 9201.714082 | −59428.755230 |
2 | −1.146887043 | 2183.687057 | 448.743795 | 9.66425 × 10−4 | 1148.244207 | −9789.990571 |
3 | −1.195130784 | 2221.779805 | 10.362946 | 3.19907 × 10−5 | 35.137584 | −394.045005 |
4 | −1.197201392 | 2223.403509 | 16.492003 | 5.52543 × 10−5 | 59.714483 | −0.691388 |
5 | −1.197204886 | 2223.406205 | 16.525002 | 5.55595 × 10−5 | 60.003229 | −0.000002 |
6 | −1.197204886 | 2223.406205 | 16.525002 | 5.55598 × 10−5 | 60.003515 | 0 |
3种随机模型的非零元素见表1. 约束下3种随机模型的拟合直线参数和精度列于表3,算法均经过6次迭代收敛,耗时均在0.6 s以内. 基于P1得到的后验单位权中误差 ${\hat \sigma _0}$ = 16.5 mm,截距的中误差${\sigma _b}$ = 60.0 mm,均为最小,相对于P2和P3得到的${\sigma _b}$ = 66.1 mm,精度提升了9.2%. 3种随机模型下的直圆点坐标(xZ,yZ)列于表3,对应位置如图4(b),分别表示为Z1、Z2和Z3,其中Z2和Z3点重合.
随机模型 | $\hat a$ | $\hat b$/m | ${\hat \sigma _0}$/mm | ${\sigma _a}$ | ${\sigma _b}$/mm | xZ/m | yZ/m | 迭代数/次 | 耗时/s |
P1 | −1.19720489 | 2223.4062 | 16.525 | 5.55598 × 10−5 | 60.0 | 1079.9809 | 930.4478 | 6 | 0.494 |
P2 | −1.19722236 | 2223.4251 | 49.031 | 6.12012 × 10−5 | 66.1 | 1079.9772 | 930.4522 | 6 | 0.503 |
P3 | −1.19722233 | 2223.4250 | 76.478 | 6.12012 × 10−5 | 66.1 | 1079.9772 | 930.4522 | 6 | 0.496 |
令B为空矩阵,用同样的参数初值进行无约束的直线拟合,3种随机模型的拟合直线参数和精度列于表4,算法均经过3次迭代收敛,耗时均在0.5 s以内,表明无约束拟合直线具有更高的效率. 且无约束时的单位权中误差及参数精度均高于约束下的直线拟合,说明在无约束情况下,测点与拟合直线的贴合度更好. 但是,无约束拟合获得的直线不相切于已知圆曲线,如图4(b)所示,相离距离d = 0.18 m,不符合线路连续性的要求. 另外,在无约束拟合时,基于P1得到的后验单位权中误差${\hat \sigma _0}$ = 2.3 mm,截距的中误差${\sigma _b}$ = 25.0 mm,相对P2和P3精度提升了2.7%. 说明不论是否有约束,考虑坐标相关误差时,获得的直线参数精度均为最高.
随机模型 | $\hat a$ | $\hat b$/m | ${\hat \sigma _0}$/mm | ${\sigma _a}$ | ${\sigma _b}$/mm | 迭代数/次 | 耗时/s |
P1 | −1.196269322 | 2222.672000 | 2.293655 | 3.10233 × 10−5 | 25.008852 | 3 | 0.461 |
P2 | −1.196259074 | 2222.663879 | 6.710739 | 3.16316 × 10−5 | 25.743943 | 3 | 0.414 |
P3 | −1.196259065 | 2222.663872 | 0.462488 | 3.16316 × 10−5 | 25.743942 | 3 | 0.438 |
在铁路维护中,全站仪获取的轨道坐标点具有误差相关性,普通最小二乘或正交最小二乘拟合直线不能考虑观测值之间的这种误差相关性. 提出的直线拟合通用模型,可以考虑坐标分量间的误差相关性实现直线拟合. 通过指定随机模型,可以实现普通最小二乘、正交最小二乘或加权整体最小二乘直线拟合,并揭示了其对应的几何意义.
通常的直线拟合方法未考虑约束条件,线路重构中直线的拟合受到相邻线元的约束. 提出的直线拟合通用模型,可以在约束下同时顾及误差相关性实现直线拟合. 顾及坐标相关误差时,可以提升估计参数的精度:约束及无约束下参数估计精度分别提高了9.2%和2.7%.
采用的高斯-牛顿算法运行效率高,能够快速得到直线参数的最佳估值及其精度. 在约束及无约束情况下分别仅6次及3次迭代就搜索出最优直线.
致谢:中南大学土木工程国家级实验教学示范中心项目(201905406).
[1] | NONG Xingzhong, SHI Haiou, YUAN Quan, ZENG Wenqu, ZHENG Qing, DING Guofu. Review on BIM Technology Used in Urban Rail Transit Projects[J]. Journal of Southwest Jiaotong University, 2021, 56(3): 451-460. doi: 10.3969/j.issn.0258-2724.20200018 |
[2] | CAI Xuan, WANG Changlin. BeiDou Navigation Satellite System/Inertial Measurement Unit Integrated Train Positioning Method Based on Improved Unscented Kalman Filter Algorithm[J]. Journal of Southwest Jiaotong University, 2020, 55(2): 393-400. doi: 10.3969/j.issn.0258-2724.20170816 |
[3] | YUAN Weina, WANG Jiaxuan. Fast Time-Varying Sparse Channel Estimation Based on Kalman Filter[J]. Journal of Southwest Jiaotong University, 2018, 53(4): 835-841. doi: 10.3969/j.issn.0258-2724.2018.04.023 |
[4] | ZHAI Wanming, ZHAO Chunfa. Frontiers and Challenges of Sciences and Technologies in Modern Railway Engineering[J]. Journal of Southwest Jiaotong University, 2016, 29(2): 209-226. doi: 10.3969/j.issn.0258-2724.2016.02.001 |
[5] | WANG Zhifu, LIU Mingchun, ZHOU Yang. Estimation of Longitudinal Speed of In-wheel Motor Driven Vehicle Using Fuzzy Extended Kalman Filter[J]. Journal of Southwest Jiaotong University, 2015, 28(6): 1094-1099. doi: 10.3969/j.issn.0258-2724.2015.06.017 |
[6] | YAO Hongguang. Multi-resolution Wavelet Decomposition of Chinese Aviation Network[J]. Journal of Southwest Jiaotong University, 2013, 26(1): 141-146. doi: 10.3969/j.issn.0258-2724.2013.01.022 |
[7] | NIE Qinghui, XIA Jingxin, QIAN Zhendong. Short-Term Traffic Flow Forecasting and Reliability Analysis of Urban Road[J]. Journal of Southwest Jiaotong University, 2013, 26(5): 955-960. doi: 10.3969/j.issn.0258-2724.2013.05.027 |
[8] | XU Yantian, CHENG Pengfei, CAI Yanhui. Kalman Filter Algorithm for Medium-Range Real-Time Kinematic Positioning with One Reference Station[J]. Journal of Southwest Jiaotong University, 2013, 26(2): 317-322,356. doi: 10.3969/j.issn.0258-2724.2013.02.020 |
[9] | ZHOU Cong, XIAO Jian, WANG Song. Application of Multirate Unscented Kalman Filter to State Estimation in Vehicle's Active Front Steering System[J]. Journal of Southwest Jiaotong University, 2012, 25(5): 849-854,894. doi: 10.3969/j.issn.0258-2724.2012.05.019 |
[10] | WU Zhizhou, FAN Yujie, MA Wanjing. Spot Speed Prediction Model Based on Grey Neural Network[J]. Journal of Southwest Jiaotong University, 2012, 25(2): 285-290. doi: 10.3969/j.issn.0258-2724.2012.02.019 |
[11] | LI Yifan, LIN Jianhui, LIU Jianxin. Wheel-Rail Force Continuous Measurement Based on Combinational Forecast Model[J]. Journal of Southwest Jiaotong University, 2012, 25(4): 597-604. doi: 10.3969/j.issn.0258-2724.2012.04.010 |
[12] | MA Lei, SHI Xizhi. Recent Development on Consensus-Based Kalman Filtering in Multi-agent Systems[J]. Journal of Southwest Jiaotong University, 2011, 24(2): 287-293. doi: 10.3969/j.issn.0258-2724.2011.02.019 |
[13] | ZHANG Ming, HAN Songchen, HUANG Linyuan. Air Traffic Flow Combinational Forecast Based on Double Gravity Model and Artificial Neural Network[J]. Journal of Southwest Jiaotong University, 2009, 22(5): 764-770. |
[14] | HU Dan, XIAO Jian, CHE Chang. Support Vector Machine with Scaling Kernel and Its Application in Dynamic System Identification[J]. Journal of Southwest Jiaotong University, 2006, 19(4): 460-465. |
[15] | YUXiang, MAO Min, LIUJian-bing. Study on Combined Models of Demand Forecasting for Urban Traffic[J]. Journal of Southwest Jiaotong University, 2003, 16(1): 75-79. |
[16] | HE Zheng-you, DAIXiao-wen, QIANQing-quan. Prefilter Selection Methods in Fast Algorithm of Discrete Wavelet Transform[J]. Journal of Southwest Jiaotong University, 2000, 13(2): 183-187. |
[17] | LIShi-ling, LIHe-sheng, LI Zhi. Application and Design of the Wavelet Filter in Weak Signal Detection[J]. Journal of Southwest Jiaotong University, 2000, 13(1): 86-89. |
[18] | DAIXiao-wen, ZHONG Gui-ying, WUHao-zhong. Application of the Wavelet Transform to the Failure Diagnosis of Tilting Trains[J]. Journal of Southwest Jiaotong University, 2000, 13(6): 651-655. |
点号 | x/m | y/m | C 中非零元素 | P1 中非零元素 | P2 中非零元素 | P3 中非零元素 | |||||||
$\sigma _x^2$/mm2 | $\sigma _y^2$/mm2 | ${\sigma _{xy}}$/mm2 | px | py | pxy | ${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 |
迭代数/次 | $\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−2 | 93044.287650 | −479524.603400 |
1 | −0.935165106 | 2013.498778 | 5844.096808 | 6.07990 × 10−3 | 9201.714082 | −59428.755230 |
2 | −1.146887043 | 2183.687057 | 448.743795 | 9.66425 × 10−4 | 1148.244207 | −9789.990571 |
3 | −1.195130784 | 2221.779805 | 10.362946 | 3.19907 × 10−5 | 35.137584 | −394.045005 |
4 | −1.197201392 | 2223.403509 | 16.492003 | 5.52543 × 10−5 | 59.714483 | −0.691388 |
5 | −1.197204886 | 2223.406205 | 16.525002 | 5.55595 × 10−5 | 60.003229 | −0.000002 |
6 | −1.197204886 | 2223.406205 | 16.525002 | 5.55598 × 10−5 | 60.003515 | 0 |
随机模型 | $\hat a$ | $\hat b$/m | ${\hat \sigma _0}$/mm | ${\sigma _a}$ | ${\sigma _b}$/mm | xZ/m | yZ/m | 迭代数/次 | 耗时/s |
P1 | −1.19720489 | 2223.4062 | 16.525 | 5.55598 × 10−5 | 60.0 | 1079.9809 | 930.4478 | 6 | 0.494 |
P2 | −1.19722236 | 2223.4251 | 49.031 | 6.12012 × 10−5 | 66.1 | 1079.9772 | 930.4522 | 6 | 0.503 |
P3 | −1.19722233 | 2223.4250 | 76.478 | 6.12012 × 10−5 | 66.1 | 1079.9772 | 930.4522 | 6 | 0.496 |
随机模型 | $\hat a$ | $\hat b$/m | ${\hat \sigma _0}$/mm | ${\sigma _a}$ | ${\sigma _b}$/mm | 迭代数/次 | 耗时/s |
P1 | −1.196269322 | 2222.672000 | 2.293655 | 3.10233 × 10−5 | 25.008852 | 3 | 0.461 |
P2 | −1.196259074 | 2222.663879 | 6.710739 | 3.16316 × 10−5 | 25.743943 | 3 | 0.414 |
P3 | −1.196259065 | 2222.663872 | 0.462488 | 3.16316 × 10−5 | 25.743942 | 3 | 0.438 |
点号 | x/m | y/m | C 中非零元素 | P1 中非零元素 | P2 中非零元素 | P3 中非零元素 | |||||||
$\sigma _x^2$/mm2 | $\sigma _y^2$/mm2 | ${\sigma _{xy}}$/mm2 | px | py | pxy | ${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 |
迭代数/次 | $\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−2 | 93044.287650 | −479524.603400 |
1 | −0.935165106 | 2013.498778 | 5844.096808 | 6.07990 × 10−3 | 9201.714082 | −59428.755230 |
2 | −1.146887043 | 2183.687057 | 448.743795 | 9.66425 × 10−4 | 1148.244207 | −9789.990571 |
3 | −1.195130784 | 2221.779805 | 10.362946 | 3.19907 × 10−5 | 35.137584 | −394.045005 |
4 | −1.197201392 | 2223.403509 | 16.492003 | 5.52543 × 10−5 | 59.714483 | −0.691388 |
5 | −1.197204886 | 2223.406205 | 16.525002 | 5.55595 × 10−5 | 60.003229 | −0.000002 |
6 | −1.197204886 | 2223.406205 | 16.525002 | 5.55598 × 10−5 | 60.003515 | 0 |
随机模型 | $\hat a$ | $\hat b$/m | ${\hat \sigma _0}$/mm | ${\sigma _a}$ | ${\sigma _b}$/mm | xZ/m | yZ/m | 迭代数/次 | 耗时/s |
P1 | −1.19720489 | 2223.4062 | 16.525 | 5.55598 × 10−5 | 60.0 | 1079.9809 | 930.4478 | 6 | 0.494 |
P2 | −1.19722236 | 2223.4251 | 49.031 | 6.12012 × 10−5 | 66.1 | 1079.9772 | 930.4522 | 6 | 0.503 |
P3 | −1.19722233 | 2223.4250 | 76.478 | 6.12012 × 10−5 | 66.1 | 1079.9772 | 930.4522 | 6 | 0.496 |
随机模型 | $\hat a$ | $\hat b$/m | ${\hat \sigma _0}$/mm | ${\sigma _a}$ | ${\sigma _b}$/mm | 迭代数/次 | 耗时/s |
P1 | −1.196269322 | 2222.672000 | 2.293655 | 3.10233 × 10−5 | 25.008852 | 3 | 0.461 |
P2 | −1.196259074 | 2222.663879 | 6.710739 | 3.16316 × 10−5 | 25.743943 | 3 | 0.414 |
P3 | −1.196259065 | 2222.663872 | 0.462488 | 3.16316 × 10−5 | 25.743942 | 3 | 0.438 |