全球导航卫星系统GNSS和加速度计具有良好的技术互补特性，工程上经常将2种技术融合起来监测结构体的变形和高频振动特性，既能抵抗复杂环境中GNSS的电磁信号干扰，又能减少加速度计积分产生的误差累积，目前已广泛应用于对大桥、高层建筑物等结构的健康监测和同震形变反演等研究领域^[1-2]。但常规的GNSS/加速度计融合解算方法需要通过高精度GNSS量测信息约束来校正加速度计基线转换偏差^[3]，进而充分利用加速度计短时精度高的特点获取准确的形变位移推算值。GNSS在复杂的遮挡环境下极易受多路径误差的影响，出现低精度异常监测结果，导致常规融合算法中大量GNSS误差被基线转换偏差所吸收，无法获取高精度的融合变形监测序列。

针对含有异常值的观测数据，常用的解决方法有数据探测法^[4]和抗差估计法^[5]，在GNSS动态导航定位中均有较多应用^[6-7]。数据探测法从函数模型的角度通过假设检验等方法剔除粗差，完全消除异常观测值对参数估计的影响。但GNSS/加速度计耦合变形监测的卡尔曼滤波器(Kalman filter, KF)量测信息仅由GNSS形变位移结果构成，不存在多余的观测信息，若直接剔除GNSS异常量测信息会导致量测更新中断，造成滤波结果的发散。本文基于Huber选权迭代法，将方差膨胀思想引入GNSS/加速度计融合KF中，从随机模型的角度对GNSS异常值的量测噪声进行自适应调整，降低GNSS异常观测对加速度计推算位置量测误差的影响，进而准确估计基线转换偏差，形成高精度的加速度计推算形变约束先验信息，提高GNSS/加速度计融合变形监测的可靠性。本文在复杂遮挡环境下实施静态和动态2种场景的变形监测实验来验证该算法的有效性。

1 GNSS/加速度计融合变形监测滤波算法

GNSS/加速度计融合变形监测滤波算法主要通过KF模型来实现^[8]。由于加速度计原始观测值存在显著的基线转换偏差，并且会随时间逐渐发散，因此将该误差扩展为待估状态参数与位移和速度一起进行估计。同时，在利用加速度推算状态时间更新时，将基线转换偏差作为闭环修正参数对加速度计的原始观测值进行校正，最后利用GNSS位移解作为量测信息进行滤波校正，获得高精度形变位移序列。具体算法流程如图 1所示。

在GNSS/加速度计融合变形监测中，所有数据解算过程均在站心坐标系E、N、U下进行，状态参数可表示为9维向量：

$\boldsymbol{X}_k=\left[\begin{array}{lll}\boldsymbol{d}_{3 \times 1} & \boldsymbol{v}_{3 \times 1} & \boldsymbol{u}_{3 \times 1}\end{array}\right]^{\mathrm{T}}$

(1)

式中，d为三维方向位移向量，v为三维方向速度向量，u为加速度计3轴方向上的基线转换偏差向量。

状态方程采用顾及加速度计基线转换偏差的变加速模型，可表示为：

$\left\{\begin{array}{l}\boldsymbol{d}_{k, k-1}=\boldsymbol{d}_{k-1}+\boldsymbol{v}_{k-1} \Delta t+0.5\left(\boldsymbol{a}_k-\boldsymbol{u}_k\right) \Delta t^2 \\ \boldsymbol{v}_{k, k-1}=\boldsymbol{v}_{k-1}+\left(\boldsymbol{a}_k-\boldsymbol{u}_k\right) \Delta t \\ \boldsymbol{u}_{k, k-1}=\boldsymbol{u}_{k-1}\end{array}\right.$

(2)

式中，下标k为时刻，Δt为加速度计采样时间间隔，a_k为k时刻加速度计的观测值，整理为矩阵形式：

$\boldsymbol{X}_{k, k-1}=\boldsymbol{\varphi}_{k, k-1} \boldsymbol{X}_{k-1}+\boldsymbol{B}_k \boldsymbol{U}_k$

(3)

式(3)可理解为KF的状态时间更新方程。由于Δt在加速度计工作期间保持不变，且该式不需要线性化，因此可认为该GNSS/加速度计融合变形监测系统是线性定常系统。式中，φ_{k, k-1}为状态转移矩阵，B_k为控制输入矩阵，U_k为控制输入向量，可分别表示为：

$\begin{array} {c} \boldsymbol{\varphi}_{k, k-1}=\left[\begin{array}{ccc}I & \Delta t & -0.5 \Delta t^2 \\ 0 & I & -\Delta t \\ 0 & 0 & I\end{array}\right], \\ \boldsymbol{B}_k=\left[\begin{array}{c}0.5 \Delta t^2 \\ \Delta t \\ 0\end{array}\right], \boldsymbol{U}_k=\boldsymbol{a}_k\end{array}$

(4)

由于加速度计是状态时间更新中的唯一噪声源，因此可以通过误差传播定律将加速度噪声q映射到位移、速度和基线转换偏差分量上来构造Kalman滤波的过程噪声矩阵Q_k：

$\boldsymbol{Q}_k=\boldsymbol{b} \boldsymbol{q} \boldsymbol{b}^{\mathrm{T}}=\left[\begin{array}{ccc}0.25 q \Delta t^4 & 0.5 q \Delta t^3 & 0.5 q \Delta t^2 \\ 0.5 q \Delta t^3 & q \Delta t^2 & q \Delta t \\ 0.5 q \Delta t^2 & q \Delta t & q\end{array}\right]$

(5)

式中，

$\boldsymbol{b}=\left[\begin{array}{ccc}-0.5 \Delta t^2 & -\Delta t & 1\end{array}\right]^{\mathrm{T}}$

(6)

状态协方差的时间更新方程可表示为：

$\boldsymbol{P}_{k, k-1}=\boldsymbol{\varphi}_{k, k-1} \boldsymbol{P}_{k-1} \boldsymbol{\varphi}_{k, k-1}^{\mathrm{T}}+\boldsymbol{Q}_k$

(7)

采用GNSS形变位移作为量测信息，构建量测方程：

$\boldsymbol{Z}_k=\boldsymbol{H}_k \boldsymbol{X}_k+\boldsymbol{V}_k$

(8)

式中，Z_k为量测向量，H_k为量测系数阵，V_k为量测噪声序列，可分别表示为：

$\boldsymbol{Z}_k=\boldsymbol{d}_{k_{-} \mathrm{GNSS}} \boldsymbol{H}_k=\left[\begin{array}{lll}\boldsymbol{I} & 0 & 0\end{array}\right]$

(9)

进一步通过计算KF增益矩阵K_k，即可解算待估状态及其协方差阵的高精度滤波解：

$\left\{\begin{array}{l}\boldsymbol{K}_k=\boldsymbol{P}_{k, k-1} \boldsymbol{H}_k^{\mathrm{T}}\left(\boldsymbol{H}_k \boldsymbol{P}_{k, k-1} \boldsymbol{H}_k^{\mathrm{T}}+\boldsymbol{R}_k\right)^{-1} \\ \boldsymbol{X}_k=\boldsymbol{X}_{k, k-1}+\boldsymbol{K}_k\left(\boldsymbol{Z}_k-\boldsymbol{H}_k \boldsymbol{X}_{k, k-1}\right) \\ \boldsymbol{P}_k=\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{H}_k\right) \boldsymbol{P}_{k, k-1}\end{array}\right.$

(10)

式中，R_k为常规的GNSS量测噪声矩阵。

2 异常GNSS观测的方差膨胀模型

在GNSS变形监测作业中，通常仅可获取形变位移序列，无法获取其监测误差。在开阔环境下，通常按照经验将GNSS量测噪声赋为定值，但在滑坡等复杂环境中，GNSS监测结果会不可避免地出现粗差，常值量测噪声无法反映GNSS位移的真实误差，进而污染Kalman滤波器导致解算结果发散。为解决上述问题，本文通过引入方差膨胀模型对GNSS异常监测结果的量测噪声矩阵进行自适应修正，以恢复GNSS/加速度计融合监测的真实形变位移序列。

方差膨胀模型可通过适当扩大异常观测值的方差，从优化随机模型的角度来降低异常观测对参数估计的影响^[9-10]。在GNSS/加速度计融合变形监测作业中，加速度计推算在短时间内的误差累积较少，变形监测误差主要来源于GNSS解算结果，因此可基于量测残差计算方差膨胀因子，进而自适应调整GNSS量测噪声。处理步骤如下：

1) 基于Kalman滤波量测方程构建标准化残差V_i：

$\overline{\boldsymbol{V}}_i=\frac{\boldsymbol{V}_i}{r_i}$

(11)

式中，下标i为三维方向E、N、U；r_i(单位m)为按照经验给出的GNSS形变位移误差；V_i为加速度计推算位移与GNSS位移的差值，即量测残差，可通过式(8)转换为：

$\boldsymbol{V}_k=\boldsymbol{Z}_k-\boldsymbol{H}_k \boldsymbol{X}_k$

(12)

2) 参考Huber法^[11]的等价权函数构建方差膨胀因子分段函数：

$\lambda_{i i}=\left\{\begin{array}{l}1, \left|\overline{\boldsymbol{V}}_i\right| \leqslant c \\ \frac{\left|\overline{\boldsymbol{V}}_i\right|}{c} \times d, \left|\overline{\boldsymbol{V}}_i\right|>c\end{array}\right.$

(13)

式中，c取1.0~1.5^[9]，c值越大说明算法能容忍的粗差量级越大；d取大于1.0的常数；标准化残差值V_i越大，说明GNSS量测信息精度越差，可相对应地扩大参数d的数值以增加方差膨胀系数，降低GNSS在融合解算中的比重。

3) 构建等价GNSS量测噪声矩阵。常规的量测噪声矩阵可表示为：

$\boldsymbol{R}_k=\left[\begin{array}{ccc}r_E^2 & r_{E N} & r_{E U} \\ r_{N E} & r_N^2 & r_{N U} \\ r_{U E} & r_{U N} & r_U^2\end{array}\right]$

(14)

量测噪声矩阵的各个元素可通过方差膨胀因子进行放大：

$\left\{\begin{array}{l}\bar{r}_i^2=\lambda_{i i} r_i^2 \\ \bar{r}_j^2=\lambda_{j j} r_j^2 \\ \bar{r}_{i j}=\lambda_{i j} r_{i j}\end{array}\right.$

(15)

式中，i和j分别为E、N、U三方向，$\lambda_{i j}=\sqrt{\lambda_{i i} \lambda_{j j}}$；r_i²、r_j²、r_ij分别为等价方差和协方差。等价GNSS量测噪声矩阵R_k可表示为：

$\begin{gathered}\overline{\boldsymbol{R}}_k=\left[\begin{array}{ccc}\bar{r}_E^2 & \bar{r}_{E N} & \bar{r}_{E U} \\ \bar{r}_{N E} & \bar{r}_N^2 & \bar{r}_{N U} \\ \bar{r}_{U E} & \bar{r}_{U N} & \bar{r}_U^2\end{array}\right]= \\ {\left[\begin{array}{ccc}\lambda_{11} r_E^2 & \lambda_{12} r_{E N} & \lambda_{13} r_{E U} \\ \lambda_{21} r_{N E} & \lambda_{22} r_N^2 & \lambda_{23} r_{N U} \\ \lambda_{31} r_{U E} & \lambda_{32} r_{U N} & \lambda_{33} r_U^2\end{array}\right]}\end{gathered}$

(16)

将等价GNSS量测噪声R_k代替式(13)中的R_k，即可降低异常GNSS形变结果的影响，提高KF系统的鲁棒性。

3 实验与分析

设计静态和动态2种实验场景对本文模型进行验证和分析。实验地点位于长安大学地学科技大厦东侧，周围存在较多楼宇和植被，易造成GNSS信号遮挡和多路径衍射。使用模拟形变平台将GNSS天线和加速度计固连，形成一个刚体结构(图 2)。GNSS监测站采用和芯星通UB4B0M型接收机，基准站使用华测RTK虚拟参考站，可以输出1 Hz的E、N、U三维形变位移，RTK的定位精度可以达到1 cm，其中高程精度略低于平面精度。基于上述结论，本文将GNSS三维方向上的量测噪声设为[0.01 0.01 0.03]，单位为m。加速度计设备使用中科院精密测量院研发的EWM1580型烈度计，可以输出100 Hz高频三轴加速度数据，通过Allan方差标定得到的噪声密度(加速度功率谱密度开方)为2.88 μg/$\sqrt{\mathrm{Hz}}$，加速度随机游走为4.65 μg·$\sqrt{\mathrm{Hz}}$，可参与构建KF过程噪声矩阵。

3.1 静态实验

本文采集到约700 s的GNSS/加速度计静态监测数据，在实验前将加速度计三轴与E、N、U方向大致对齐，尽可能降低基线转换偏差的影响。为验证本文所提算法在静态变形监测中的可靠性，实验设计GNSS、常规GNSS/加速度计融合算法(常规融合算法)和基于方差膨胀模型的GNSS/加速度计融合算法(改进融合算法)3种方案进行对比。

图 3为3种解算方法的监测序列结果，并以0为参考值分别统计均方根RMS和误差最大值，如表 1(单位m)所示。由图 3(b)可见，当GNSS监测结果无明显异常形变位移时，3种方法的变形监测序列基本一致。由图 3(c)可见，当GNSS受到多路径信号干扰时，其监测序列(蓝色曲线)出现较多的异常形变结果，E、N、U三维方向的RMS分别为0.311 m、0.250 m、0.657 m，最大值分别达1.468 m、1.086 m、0.741 m，其精度已无法满足高精度变形监测的要求。常规融合算法(绿色曲线)获取的监测序列略优于GNSS，其三维方向的RMS分别为0.297 m、0.227 m、0.612 m，相较于GNSS分别提升4.5%、9.2%、6.8%。但其误差最大值比GNSS更加突出，这是由于异常GNSS量测信息会引起加速度计基线转换偏差的错误校正，最终导致融合结果显著发散。改进融合算法(红色曲线)在三维方向的RMS分别为0.006 m、0.003 m、0.011 m，相较于GNSS分别提升98.1%、98.8%、98.3%，误差最大值分别为0.013 m、0.004 m、0.025 m。由此可见，改进融合算法显著抑制了GNSS的异常监测结果，提升了变形监测的可靠性。

图 4为基于改进融合算法获取的量测残差(红色曲线)和方差膨胀因子(蓝色柱形)时间序列，其中高程方向的量测残差量级显著大于水平方向。由图 4(b)可见，当GNSS无异常观测时，量测残差均小于1 m，此时方差膨胀因子为1，保持不变；由图 4(c)可见，当GNSS出现异常观测时，量测残差出现较大浮动，高程方向甚至出现大于2 m的结果，此时方差膨胀因子也随之扩大，需要通过构建等价量测噪声阵抑制GNSS变形监测中的异常结果。

图 5为加速度计原始观测值(黑色曲线)及常规融合算法(绿色曲线)和改进融合算法(红色曲线)获取的加速度计基线转换偏差序列。由图 5(b)可见，当GNSS无明显异常观测时，2种融合方法获取的基线转换偏差基本一致，且均无明显的错误估计现象；由图 5(c)可见，当GNSS出现异常观测时，常规融合算法估计的基线转换偏差的量级明显大于加速度原始观测值，且处于发散状态，进而导致了形变位移的错误估计。改进融合算法通过优化GNSS量测噪声矩阵，使解算得到的基线转换偏差序列变化趋势与加速度原始观测值一致，无明显的错误估计现象，这也保证了加速度推算能够获取准确的形变位移约束信息。

3.2 动态实验

本文采集了一段时长约500 s的GNSS/加速度计融合变形监测数据进行验证分析，在该时段中使用模拟形变平台分别在水平和高程方向实施动态形变测试(表 2)。需要注意的是，当加速度计处于动态状态时，其原始观测数据的噪声会比静态更加显著，在数据处理时应对加速度噪声q进行适当放大。

为验证本文算法在动态情况下的有效性，同样设计GNSS、常规融合算法和改进融合算法3种方案进行对比，图 6为3种解算方法的变形监测序列。当GNSS无明显异常观测时，3种解算方法获取的结果基本一致。由于GNSS异常观测主要出现在形变期间，因此本文仅统计2个形变时段的RMS和误差最大值，如表 3(单位m)所示。由图 6(b)可见，时段1内GNSS同样出现较多的监测误差，E、N、U三维方向的RMS分别为0.385 m、0.151 m、0.361 m，最大值达0.939 m、0.438 m、1.289 m，若这种异常结果发生在形变阶段，极易引起错误决策，导致严重的工程事故；常规融合算法的监测精度与GNSS相当，RMS仅有小幅度下降；改进融合算法获取的形变位移精度有明显改善，三维方向的RMS均在1.8 cm以内，误差最大值均在4.7 cm以内。由图 6(c)可见，时段2内GNSS监测序列在三维方向的RMS分别为0.072 m、0.014 m、0.307 m，该时段内GNSS异常观测结果较少；常规融合算法的变形监测精度相较于GNSS有一定程度的提高，但仍然维持在较低水平；改进融合算法在三维方向的RMS均在1.3 cm以内，误差最大值不超过3.5 cm，能够有效还原真实的形变位移序列。由于GNSS和常规融合算法的位移序列错误离群值较大，在一定程度上湮没了真实的位移变化情况，因此图 6(b)的E方向和图 6(c)的U方向分别给出改进融合算法解算得出的位移序列，以更清晰地显示2个时段的形变位移变化数值。

图 7为改进融合算法获取的量测误差(红色曲线)和方差膨胀因子(蓝色柱形)时间序列结果，由图可见，时段1的GNSS异常历元数显著大于时段2。量测误差与静态实验相比有一定程度的降低，但也同样触发了不同量级的方差膨胀，以获取更加合理的GNSS量测噪声矩阵。

图 8为加速度计原始观测值(黑色曲线)及常规融合算法(绿色曲线)和改进融合算法(红色曲线)解算得到的加速度计基线转换偏差序列。由图可见，在位移发生期间，加速度观测值出现显著变化，同时基线转换偏差也出现细微变化，符合实际情况。由图 8(b)可见，时段1内常规融合算法获取的基线转换偏差仍然存在明显的发散现象，甚至已经湮没原始加速度观测数据，这也会进一步导致加速度计推算位移结果的发散；而改进融合算法获取的序列更加合理，始终在加速度计观测值的变化范围内。由图 8(c)可见，由于时段2内GNSS异常观测历元较少，常规融合算法获取的序列结果仅有短时间的发散现象，但仍然差于改进融合算法获取的基线转换偏差序列结果。

4 结语

1) 当GNSS没有明显异常观测时，常规融合算法与改进融合算法获取的形变位移基本一致，均能得到可靠的变形监测结果。

2) 当GNSS出现异常观测时，常规融合算法的精度提升有限，无法满足高精度变形监测的要求；改进融合算法可根据量测残差自适应调整方差膨胀因子，进而优化GNSS量测噪声矩阵，获取高精度的形变位移序列和加速度计基线转换偏差结果。静态和动态实验结果表明，改进融合算法的形变位移RMS均在1.8 cm以内，相较于GNSS，其监测结果的提升幅度均在90%以上，可在复杂遮挡环境下提供可靠的变形监测结果。