形变是普遍存在于自然界事物中的一种现象，它是指形变体在各种载荷作用下，其形状、尺寸与位置在空间域与时间域中的变化。所谓形变监测，是对形变体的形变现象采用精密测量仪器设备和测量方法进行监视观测的过程，其目的是确定形变体在各种载荷与外力作用下，形状、尺寸及位置变化的空间状态与时间特征。

形变监测技术包括常规大地测量技术、特殊形变测量技术和全球导航卫星系统(Global Navigation Satellite System，GNSS)技术。常规大地测量，是采用经纬仪、水准仪、测距仪、全站仪等常规测量仪器测定点的形变值，其优点是：(1) 能够提供形变体整体的形变状态；(2) 适用于不同的监测精度要求、不同形式的形变体和不同的监测环境；(3) 可以提供绝对形变信息^[1]。特殊形变测量技术是指随着光、电、机技术的发展，人们利用专门研制的一些特殊和专用的监测仪器对事物、建筑物进行形变监测，具有测量过程简单、可监测形变体内部的形变、容易实现自动化监测等优点,但通常只能提供局部的和相对的形变信息，如文献[2-3]中提到的两种特殊测量方法和仪器。常规大地测量技术与特殊形变测量技术受限于巨大的设备成本和精密的设备要求，难以长时间在恶劣环境中工作。GNSS在形变监测上的应用随着时代的发展和技术的进步，变得更加广泛。GNSS能够很好地应对前两种技术在形变监测中难以适应的场景，包括滑坡在内的自然灾害，用于桥梁、轨道等的形变监测。有助于专家更好地了解建筑物的形变状态与趋势，进而提前预防灾害，如文献[4-5]中提出的两种GNSS的应用情况。然而在运用GNSS对形变体进行定位监测时，需要对形变体监测位置的三维空间坐标进行一段时间的采样，并将这段观测时长内的所有采样数据进行解算得到一个相对更加精确的坐标值。而在缩短观测时长的情况下，采样数据的减少导致采样数据解算过程中，一些多路径误差难以消除，使定位精度下降。本文就此问题提出一种通过粒子滤波算法，建立形变体监测点的状态空间模型，结合GNSS观测值，得到监测点的最优状态估计。通过对高铁铁轨路基的形变监测验证方法的可行性。

1 BDS数据采集与处理 1.1 数据采集方案

为了开展基于BDS形变监测的高铁铁轨路基沉降形变的研究分析，本文利用了位于广东省内广州市与汕头市之间广汕高速铁路搭建的BDS形变监测系统。其中，基准站B001位于主控楼楼顶，主控楼的持力层为中等风化凝灰层，同时在施工桩基前做了强化夯实处理，基桩深度达30 m，成为整个广汕高铁的BDS形变监测系统项目中最稳固的地方。监测站M001、M002、M003、M004、M005、M006以及基准站B001均位于该铁路路段的一侧，6个监测站均位于地面，距离基准站最近的是监测站M001，距离为60 m，距离基准站最远的监测站是M006，距离为350 m，同时M006监测站也最靠近铁轨。6个监测站底座安装位置的路基接触点即为一一对应的M1、M2、M3、M4、M5、M6监测点，基准站对应的基准点为B1。基准站和监测站的位置分布如图1所示。基准站和监测站接收机都采用上海司南导航M300 Plus北斗/GNSS高精度形变接收机，同时各测站上均配置了AT340测量型天线，天线安装在观测墩上，观测墩严格水平。

1.2 采集周期与频率

只在特殊情形下，高速铁路路基上才会发生瞬间对路基的沉降形变量造成突变的情况(如山体滑坡、地震及突然施工等)。一般场景下，高速铁路路基形变体的形变极为缓慢，其在局部时间域内可以认为是稳定的，即时间域内形变体的形变量可忽略不计，高铁路基形变体的局部时间域为半年。故对于在局部时间域内的高铁路基的形变监测就可以根据监测出的数据分析达到预警的目的。本文中广汕高铁路基形变监测的监测周期为1个月(2018年5月15日-2018年6月15日)。该监测系统为24 h实时在线监测系统。监测站与基准站的采集BDS定位数据的采样间隔为15 s。

1.3 数据传输

数据传输网络使用星形拓扑结构，当监控系统某一个网点发生故障时，不会影响其他网点的正常工作。监测站主机与各监测站的接受机进行实时数据交换，并转为RINEX文件储存^[6]，其主要流程如图2所示。

1.4 数据解算处理及路基监测点精度

形变监测集成软件(Dam Deformation Monitoring Software，DDMS)由总控、数据处理、数据分析和数据库管理４个模块组成，主要功能包括 BDS基线解算、坐标估计、形变量提取以及精度分析、基准稳定性分析、形变过程线显示、预警报警等^[7]。基于形变监测系统中基线短的特点，且需保证实时性和高精度的要求，广汕高铁路基BDS形变监测系统项目中关于BDS基线解算模型软件的开发主要借鉴参考了文献[7]中GPS基线解算软件的开发。具体模型以及处理策略受篇幅限制在这里不再详述。

广汕高铁路基监测的一个月周期处在局部时间域内，同时在此周期内的项目范围5公里未发生会导致路基产生突变式形变的特殊情况，因此认为在监测周期内项目的所有监测点与基准点的空间位置坐标不变。

本文分别收集了M006监测站与B001基准站于2018年5月30日(0：00~24：00)采集的BDS定位数据，M001监测站与B001基准站于2018年5月28日(0：00~24：00)采集的BDS定位数据，M003监测站与B001基准站于2018年5月29日(0：00~24：00)采集的BDS定位数据。选取M006监测站与B001基准站之间收集的BDS定位数据，结合开发的BDS基线解算模型软件将定位数据按每小时解算一次，解算出24个M6监测点的三维空间坐标值。将M6监测点的24个坐标值取平均，均值坐标即作为M6监测点的真值坐标。由此统计出M6监测点的水平方向(径向A、切向B)精度在3 mm内，高程H方向的精度在5 mm内。其误差序列如图3所示。其横轴 $ {K}_{1\;\mathrm{h}} $ 表示每小时解算一次的离散解算序列，纵轴 $ {\Delta }_{1\;\mathrm{h}} $ 表示与每小时解算一次的 $ {K}_{1\;\mathrm{h}} $ 时刻点一一对应的M6监测点的误差。而在将M6监测点坐标值的解算时间从每小时解算一次缩短为每15 min解算一次，即在采样数据减少为原来的1/4的情况下，因每次解算所用的采样数据大量减少，导致每15 min解算一次所统计的M6监测点的水平面方向 $ (A,B) $ 误差范围扩大到8 mm，高程 $ H $ 方向的误差范围扩大到11 mm。每15 min解算一次的采样数据来自2018年5月30日(0：00~24：00)，共解算出96个坐标值。其误差序列如图4所示，其横轴 $ {K}_{15\;\mathrm{m}\mathrm{i}\mathrm{n}} $ 表示每15 min解算一次的离散解算序列，纵轴 $ {\Delta }_{15\;\mathrm{m}\mathrm{i}\mathrm{n}} $ 表示与每15 min解算一次的 $ {K}_{15\;\mathrm{m}\mathrm{i}\mathrm{n}} $ 时刻点一一对应的M6监测点的误差。

由图3与图4比较，能明显看到由于解算时取用的采样数据减少，其在解算过程中消去的多路径误差变少，导致解算坐标值的精度明显下降。对于这种情形，本文尝试采用一种粒子滤波算法，结合对监测点解算坐标值(此处解算坐标值即粒子滤波算法中的观测数据)的分析，构建监测点坐标值数据的系统模型。在系统模型中引入观测数据，递推出观测数据同时刻的最优估计值，达到减少多路径误差并提升精度的目的。

2 研究方法 2.1 粒子滤波的基本原理与公式

粒子滤波是通过寻找一组在状态空间中传播的随机样本对后验概率密度函数进行近似模拟，然后以样本均值代替积分运算，从而获得状态最小方差估计的过程，而这些样本即被称为“粒子”^[8-9]。是一种通过非参数化蒙特卡洛模拟方法实现递推贝叶斯估计的算法，可以通过状态模型与观测模型由前一时刻的状态值递推得到当前时刻的状态值^[10-11]。

对于高铁铁轨路基形变监测系统，采用粒子滤波算法递推的方式。其采用系统本身的状态转移方程和在缩短观测时长的观测资料之间，求得系统所处时刻的最优估计值。因而在一定程度上使用粒子滤波器可以解决由于BDS观测时长减少引起监测点定位精度下降的问题。

粒子滤波算法的数学模型包含状态转移方程和观测方程，其模型一般形式如下

状态转移方程

$ {\boldsymbol{X}}_{K}=f\left({\boldsymbol{X}}_{K-1}\right)+{\boldsymbol{Q}}_{K} $

(1)

观测方程

$ {\boldsymbol{Y}}_{K}=h\left({\boldsymbol{X}}_{K}\right)+{\boldsymbol{R}}_{K} $

(2)

式(1)~(2)中：正整数 $ K $ 为离散时间变量； $ \boldsymbol{X}、\boldsymbol{Y} $ 分别为系统状态值序列与系统观测值序列； $ f、h $ 分别为状态转移函数与观测函数； $ \boldsymbol{Q}、\boldsymbol{R} $ 分别为状态噪声和观测噪声。

假设状态噪声与观测噪声均服从高斯分布，可以得到粒子滤波递推公式如下：

(1) 给初值 $ {\boldsymbol{X}}_{0} \sim N\left({\boldsymbol{\mu }}_{0},{\boldsymbol{\sigma }}_{0}^{2}\right) $ ，即初值服从高斯分布。

(2) 粒子初始化：由先验概率分布 $ p\left({\boldsymbol{X}}_{0}\right) $ 生成采样粒子 $ {\boldsymbol{x}}_{0}^{\left(i\right)} $ ， $ i $ =1,2,···, $ {L} $ ； $ {L} $ 为采样粒子个数。

(3) 预测步。生成

$ {\boldsymbol{x}}_{1}^{\left(i\right)}=f({\boldsymbol{x}}_{0}^{\left(i\right)})+{\boldsymbol{q}}^{\left(i\right)} $

(3)

式中： $ {\boldsymbol{q}}^{\left(i\right)} $ 为~ $ N\left({\bf{0}},\boldsymbol{Q}\right) $ 的随机数，共 $ {L} $ 个。

(4) 更新步。

假设观测噪声为服从方差 $ \boldsymbol{R} $ 的高斯分布，设观测值为 $ {\boldsymbol{Y}}_{1} $ ，计算出 $ {\boldsymbol{x}}_{1}^{\left(i\right)} $ 粒子的对应权重 $ {\boldsymbol{w}}_{1}^{\left(i\right)} $ 。

$ {\boldsymbol{w}}_{1}^{\left(i\right)}={f}_{\boldsymbol{R}}\left[{\boldsymbol{Y}}_{1}-\boldsymbol{h}({\boldsymbol{x}}_{1}^{\left(i\right)})\right]\times {\boldsymbol{w}}_{0}^{\left(i\right)} =\frac{1}{\sqrt{2\text{π }\boldsymbol{R}}}\times {\mathrm{e}}^{-\frac{{\left[{\boldsymbol{Y}}_{1}-\boldsymbol{h}\left({\boldsymbol{x}}_{1}^{\left(i\right)}\right)\right]}^{2}}{2\boldsymbol{R}}}\times {\boldsymbol{w}}_{0}^{\left(i\right)} $

(4)

(5) 权重归一化。

$ {\boldsymbol{w}}_{1}^{\left(i\right)}=\frac{{\boldsymbol{w}}_{1}^{\left(i\right)}}{\displaystyle\sum {\boldsymbol{w}}_{1}^{\left(i\right)}} $

(5)

(6) 粒子重采样。

对于粒子滤波算法应用过程中存在粒子退化的问题，即算法经过几次迭代后，绝大多数粒子权重变到很小甚至可忽略，而只有少部分粒子所占权重很大。随着无效采样粒子数目的增加，使得大量的计算浪费在对后验滤波概率分布几乎不起作用的粒子估计上，使算法估计性能下降^[12-14]。

重采样思路：在保持粒子数目不变的同时，将那些权重极小的不起作用的粒子舍去，用一些新的粒子来取代它们。通过将剩余不同权重的粒子按照它们自身权重在总权重中所占的比例复制。大权重粒子多复制，相对较小权重粒子少复制，最后粒子数补充恢复到原来粒子数。过程如下：

从均匀分布中 $ U({{0}},{{1}}] $ 中随机抽取 $ {L} $ 个独立分布随机样本 $ {\left\{{\boldsymbol{u}}_{1}^{\left(n\right)}\right\}}_{n=1}^{{L}} $ ,将抽取到的样本与归一化权值的累计和进行比较以选取粒子^[15]。从选取粒子的角度出发，在第 $ n $ 次选择中，若有权值累计和符合式(6)条件的粒子 $ {\boldsymbol{x}}_{1}^{\left(i\right)} $ ，则该抽取样本粒子选中^[16]； $ n $ =1,2,···, $ {L} $ ， $ j $ =1,2,···, $ {i} $ 。

$ \sum\nolimits _{j=1}^{i-1}{\boldsymbol{w}}_{1}^{\left(j\right)} < {\boldsymbol{u}}_{1}^{\left(n\right)}\leqslant \sum\nolimits _{j=1}^{i}{\boldsymbol{w}}_{1}^{\left(j\right)} $

(6)

最终可获得重采样之后的 $ K=1 $ 时刻的新粒子集 $ {\left\{{\boldsymbol{x}}_{1}^{\left(i\right)}\right\}}_{i=1}^{{L}} $ 。

(7) 状态最优估计值求取：由上述步骤(1)~(6)可得 $ K=1 $ 时刻的粒子 $ {\boldsymbol{x}}_{1}^{\left(i\right)} $ 与对应权重 $ {\boldsymbol{w}}_{1}^{\left(i\right)} $ ，则 $ K=1 $ 时刻的状态估计值为

$ \widehat{{\boldsymbol{X}}_{1}}=\sum\nolimits_{i=1}^{{L}}{\boldsymbol{w}}_{1}^{\left(i\right)}\times {\boldsymbol{x}}_{1}^{\left(i\right)} $

(7)

以上步骤通过粒子滤波求 $ K=1 $ 时刻的系统状态最优估计值，然后进入下一循环，将粒子 $ {\boldsymbol{x}}_{1}^{\left(i\right)} $ 引入状态转移方程求得 $ K=2 $ 时刻的粒子集 $ {\boldsymbol{x}}_{2}^{\left(i\right)} $ ；再由更新步求得权重 $ {\boldsymbol{w}}_{2}^{\left(i\right)} $ ，权值归一化等，从而得到K=1,2,3···的系统状态最优估计值 $ \widehat{{\boldsymbol{X}}_{1}} $ 、 $ \widehat{{\boldsymbol{X}}_{2}} $ 、 $ \widehat{{\boldsymbol{X}}_{3}} $ ···。其流程如图5所示。

2.2 高铁铁轨路基形变监测模型

在本文中，BDS是作为三维铁轨形变监测自动化系统的工具，使用BDS监测铁轨形变时，监测点的位置是工程本身独立坐标系中的坐标 $ (A,B,H) $ 。

根据粒子滤波 (Particle Filter，PF)模型，状态空间模型的建立是应用粒子滤波算法在实际问题中处理预测的基础。铁轨形变体的形变在局部时间域内铁轨的空间状态是稳定的情况下，本文观测系统BDS观测值的时间序列包含于这个局部时间域中。

本文采用一阶自回归模型作为局部时间域内铁轨形变的状态空间方程，采用一元线性回归模型作为铁轨形变的观测方程。具体的系统模型如下：

状态转移方程

$ {{\boldsymbol{X}}_{K}={\boldsymbol{\varphi }}_{1}{\boldsymbol{X}}_{K-1}+{\boldsymbol{Q}}_{K}\text{，}\boldsymbol{Q}}_{K} \sim N\left({\bf{0}},{\boldsymbol{\sigma }}_{{\boldsymbol{Q}}_{K}}^{2}\right) $

(8)

观测方程

$ {\boldsymbol{Y}}_{K}={\boldsymbol{\varphi }}_{2}{\boldsymbol{X}}_{K}+{\boldsymbol{R}}_{K},{\boldsymbol{R}}_{K} \sim N({\bf{0}},{\boldsymbol{\sigma }}_{{\boldsymbol{R}}_{K}}^{2}) $

(9)

式(8)~(9)中： $ {\boldsymbol{X}}_{K} $ 为三维状态向量， $ {\boldsymbol{Y}}_{K} $ 为三维观测向量； $ {\boldsymbol{\varphi }}_{1} $ 为 $ 3\times 3 $ 状态转移系数矩阵， $ {\boldsymbol{\varphi }}_{2} $ 为 $ 3\times 3 $ 观测系数矩阵； $ {\boldsymbol{Q}}_{K} $ 为三维状态噪声向量，是服从于三维零期望向量、三维 $ {\boldsymbol{\sigma }}_{{\boldsymbol{Q}}_{K}}^{2} $ 方差向量的三维高斯噪声向量； $ {\boldsymbol{R}}_{K} $ 为三维观测噪声向量，是服从于三维零期望向量、三维 $ {\boldsymbol{\sigma }}_{{\boldsymbol{R}}_{K}}^{2} $ 方差向量的三维高斯噪声向量。

2.3 参数初始化

初始参数的设置对PF模型的预测具有重要作用，合适的模型参数可以快速捕捉系统状态，提高跟踪效率^[14]。其中主要的初始化参数包括 $ {\boldsymbol{X}}_{0}，{\boldsymbol{\mu }}_{0}，{\boldsymbol{\sigma }}_{0}^{2}，{\boldsymbol{\varphi }}_{1}，{\boldsymbol{\varphi }}_{2} $ 和 $ {L} $ 。本文设定粒子数 $ {L} $ 为200，离散时间变量 $ K=96 $ ；假设系统初始状态向量 $ {\boldsymbol{X}}_{0} $ 服从以初始观测值向量 $ {\boldsymbol{Y}}_{1} $ 为期望向量的正太分布，即 ${\boldsymbol{\mu }}_{0}={\boldsymbol{Y}}_{1}= \left[\begin{array}{c}\;\;\;0.511\;5\\ \;\;\;1.516\;3\\ -0.067\;7\end{array}\right]$ mm，根据系统观测值向量序列可求得方差向量 ${\boldsymbol{\sigma }}_{0}^{2}=\left[\begin{array}{c}1.3\;\;\;\;\;\\ 1.102\\ 2.32\;\;\;\end{array}\right]$ mm。根据铁轨在局部时间域内的形变特性，有 $ {\boldsymbol{\varphi }}_{1}={\boldsymbol{\varphi }}_{2}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right] $ 。模型参数初始化后，按照上述粒子滤波算法递推流程，将离散时间变量 $ K $ 对应的观测向量 $ {\boldsymbol{Y}}_{K} $ 引入建构的粒子滤波器中对系统状态进行递推估计，完成并输出对系统状态空间的最优估计 $ \widehat{{\boldsymbol{X}}_{K}} $ 。

3 实验仿真与结果

本文原始采样数据分别来自广汕高铁铁轨路基BDS形变监测系统于2018年5月30日(0：00~24：00)时间段对此段高铁的M6监测点持续采样收集，2018年5月28日(0：00~24：00)时间段M1监测点持续采样收集和2018年5月29日(0：00~24：00)时间段对M3监测点持续采样收集。

首先对M6的采集数据进行解算仿真；从使用观测时长为1 h的采样数据解算出一个解到使用观测时长为15 min的采样数据解算出一个解后，解算出96组有效数据。数据结果包含有监测点的三维空间位置信息：径向方向( $ A $ )、切向方向( $ B $ )与高程方向( $ H $ )。采用粒子滤波算法对96组数据进行滤波除误差处理，可以得到该数据序列的状态向量 $ {\boldsymbol{X}}_{K}= \left[\begin{array}{c}{A}_{K}\\ {B}_{K}\\ {H}_{K}\end{array}\right] $ 的状态最优估计向量 $ \widehat{{\boldsymbol{X}}_{K}}=\left[\begin{array}{c}\widehat{{A}_{K}}\\ \widehat{{B}_{K}}\\ \widehat{{H}_{K}}\end{array}\right] $ 。如图6~8分别为该监测点 $ A\mathrm{、}B\mathrm{、}H $ 方向的观测数据序列和PF滤波结果；其纵轴 $ {\varDelta }_{15\;\mathrm{m}\mathrm{i}\mathrm{n}\left(A\right)} $ 、 $ {\varDelta }_{15\;\mathrm{m}\mathrm{i}\mathrm{n}\left(B\right)} $ 、 $ {\varDelta }_{15\;\mathrm{m}\mathrm{i}\mathrm{n}\left(H\right)} $ 依次表示与每15 min解算一次的 $ {K}_{15\;\mathrm{m}\mathrm{i}\mathrm{n}} $ 时刻点一一对应的M6监测点 $ A\mathrm{、}B\mathrm{、}H $ 3个方向的误差。

由图6~8显示效果可知，建立的粒子滤波模型是合理可靠的，从图中能明显看到滤波估计值的精度比原始解算观测值精度有所提高。监测点每15 min解算一次的原始解算观测值 $ A\mathrm{、}B\mathrm{、}H $ 3个方向的误差范围经过滤波后分别从 $ (-8.5\;{\rm{mm}},5.5\;{\rm{mm}}) $ 、 $ (-3.2\;{\rm{mm}}, 6.8\;{\rm{mm}}) $ 、 $ (-\mathrm{7\;{\rm{mm}},11.8\;{\rm{mm}}}) $ 误差范围缩小至 $ (-2\;{\rm{mm}}, 1.5\;{\rm{mm}}) $ 、 $ (-1\;{\rm{mm}},2\;{\rm{mm}}) $ 、 $ (-2\;{\rm{mm}},4\;{\rm{mm}}) $ 的范围内。

为进一步验证方法的可靠性，本文在此仿真的基础上分别对2018年5月28日(0：00~24：00)时间段M1监测点的采样数据以及2018年5月29日(0：00~24：00)时间段M3监测点的采样数据以相同的观测时长15 min的采样数据解算一次的方式加入算法进行仿真验证，得到的监测点监测误差范围变化情况如表1所示。

表1中的 $ {\varDelta }_{A}、{\varDelta }_{B}、{\varDelta }_{H} $ 分别表示在加入粒子滤波算法前各监测点在不同时期以15 min的采样数据解算一次的坐标 $ A\mathrm{、}B\mathrm{、}H $ 3个方向误差范围，而 $ {\varDelta }_{A\left(\mathrm{P}\mathrm{F}\right)}、{\varDelta }_{B\left(\mathrm{P}\mathrm{F}\right)}、{\varDelta }_{H\left(\mathrm{P}\mathrm{F}\right)} $ 则表示加入粒子滤波算法(PF)后，各监测点在不同时期以15 min的采样数据解算一次的坐标 $ A\mathrm{、}B\mathrm{、}H $ 3个方向误差范围。可以从表中明显看到，加入粒子滤波算法后，各监测点的误差范围都在 $ \pm 5 $ mm以内，证实了方法的可靠性。

4 仿真结果总结

本文利用粒子滤波算法较好地解决了由于观测时长减少，监测点采样数据量减少而带来的监测点解算坐标值精度下降的问题。通过建立与监测点状态空间向量契合的状态转移方程同BDS采样数据的解算观测值相结合，去除了某些周期性误差的影响，求得监测点最优估计向量，从而提高了定位精度。通过对广汕高铁铁轨路基的实测数据进行分析，验证了该方法的可行性，为减少高铁空窗期维护时间以提升高铁运行效率提供有效方法。