GNSS定位技术具有高精度、全天候、测站间无需通视、可连续实时获取三维变形信息等优势，被广泛应用于地质灾害监测领域^[1-2]。根据监测应用场景和需求的不同，常用的GNSS监测定位技术包括静态相对定位、实时动态相对定位(RTK)、精密单点定位(PPP)等^[3]。其中，RTK监测技术通过基准站与监测站的高精度载波相位差分解算，能够实时获取监测站高精度三维坐标，精度可达mm级，可实现对滑坡体的长期连续监测，因而成为目前应用最广泛的GNSS地质灾害监测技术^[4-5]。

滑坡体变形从发育到破坏一般会经历初始变形-等速变形-加速变形3个特征阶段，在不同阶段，监测体的状态变化呈现不同的规律^[6]。为有效抵抗观测粗差、抑制观测噪声、提取真实形变信号，GNSS监测中通常采用卡尔曼(Kalman)滤波估计方法，将监测体状态预测方程作为约束信息与当前历元观测信息进行融合，从而提升定位精度与可靠性^[7]。然而，当监测体处于加速变形阶段时，其三维坐标在前后历元间发生显著变化，滤波过程噪声设置不合理可能导致滤波发散^[8]。因此，有必要引入监测体位移量的实时获取方法，实现对监测体快速位移的实时准确识别。

GNSS载波相位历元间差分(TDCP)通过对单站载波相位观测值进行历元间求差，估计出该站在历元间的坐标相对变化量，精度可达cm至mm级，因而可用于探测识别监测体发生的历元间快速变形^[9-10]。然而，TDCP技术当前多用于动态定位中确定载体运动状态信息^[11]，尚无在地质灾害监测领域的应用研究。

本文提出一种基于TDCP位置变化量的地质灾害监测体快速变形识别及Kalman滤波过程噪声调节方法，通过TDCP解算监测体前后历元间三维坐标变化量，识别监测体发生的历元间显著位移，进而对Kalman滤波的过程噪声进行实时调节，可在抑制观测噪声、提高监测精度的同时，及时准确地反映监测体的快速变形。

1 基于Kalman滤波的GNSS实时相对定位模型

GNSS高精度地质灾害位移监测中，广泛采用RTK技术获取监测站相对于基准站的位移变化，精度可达cm级甚至mm级^[12]。短基线RTK定位中，通常使用如下双差非组合观测模型：

$ \begin{aligned} & \nabla \Delta P_f^{i j}=\nabla \Delta \rho^{i j}+\nabla \Delta \varepsilon_f^{i j} \\ & \lambda_f \nabla \Delta \varphi_f^{i j}=\nabla \Delta \rho^{i j}+\lambda_f \nabla \Delta N_f^{i j}+\nabla \Delta \xi_f^{i j} \end{aligned} $

(1)

式中，▽Δ为双差算子，i、j分别为参考星与非参考星，f为频率，P和φ分别为伪距和载波相位观测值，ρ为卫星与地面的几何距离，λ_f为频率f对应的波长，N为整周模糊度，ε、ξ分别为伪距和相位观测噪声。

GNSS实时定位会不可避免地受到观测误差的影响，导致定位结果出现异常波动。为抑制观测噪声的影响，提升监测精度和可靠性，通常使用Kalman滤波器进行参数估计，以便将监测体状态预测信息与观测信息进行有效融合。

Kalman滤波的状态方程和观测方程为：

$ \begin{aligned} & \boldsymbol{X}_k=\boldsymbol{\varPhi}_{k, k-1} \boldsymbol{X}_{k-1}+\boldsymbol{w}_k \\ & \boldsymbol{Z}_k=\boldsymbol{H}_k \boldsymbol{X}_k+\boldsymbol{\eta}_k \end{aligned} $

(2)

式中，X_k－1和 X_k分别为第k－1和第k历元的状态向量，即监测站三维坐标和双差模糊度参数，可表示为X_k=[x y z N^T]_k^T，其中N为双差模糊度列向量；Φ_{k, k－1}为从状态X_k－1到状态X_k的状态转移矩阵，w_k为系统噪声向量，Z_k为k历元的观测向量，H_k为设计矩阵，η_k为观测噪声。

Kalman滤波主要包含状态预测和测量更新2个环节。在状态预测环节，状态向量及其方差-协方差矩阵通过状态方程递推实现对待估状态的预测：

$ \begin{aligned} & \hat{\boldsymbol{X}}_{k, k-1}=\hat{\boldsymbol{X}}_{k-1} \\ & \boldsymbol{D}_{\hat{\boldsymbol{X}}_{k, k-1}}=\boldsymbol{D}_{\hat{\boldsymbol{X}}_{k-1}}+\boldsymbol{D}_{\boldsymbol{w}_k} \end{aligned} $

(3)

式中，$ \hat{\boldsymbol{X}}_{k, k-1}$和$\boldsymbol{D}_{\hat{\boldsymbol{X}}_{k, k-1}} $为由上一历元的状态向量预测得到的当前历元状态向量及其方差-协方差矩阵，D_wk为过程噪声方差-协方差矩阵。相位观测值未发生周跳时，前后历元间的双差模糊度参数保持不变；监测体处于稳定状态时，前后历元间的监测站三维坐标参数保持不变或变化较小，因此状态转移矩阵Φ_{k, k－1}为单位阵。

根据当前历元的观测值对状态预测值进行测量更新，可以得到：

$ \begin{aligned} & \boldsymbol{K}_k=\boldsymbol{D}_{\hat{\boldsymbol{X}}_{k, k-1}} \boldsymbol{H}_k^{\mathrm{T}}\left[\boldsymbol{H}_k \boldsymbol{D}_{\hat{\boldsymbol{X}}_{k, k-1}} \boldsymbol{H}_k^{\mathrm{T}}+\boldsymbol{D}_{\eta_k}\right]^{-1} \\ & \boldsymbol{V}_k=\boldsymbol{Z}_k-\boldsymbol{H}_k \hat{\boldsymbol{X}}_{k, k-1} \\ & \hat{\boldsymbol{X}}_k=\hat{\boldsymbol{X}}_{k, k-1}+\boldsymbol{K}_k \boldsymbol{V}_k \\ & \boldsymbol{D}_{\hat{\boldsymbol{X}}_k}=\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{H}_k\right) \boldsymbol{D}_{\hat{\boldsymbol{X}}_{k, k-1}}\left(\boldsymbol{I}-\boldsymbol{K}_k \boldsymbol{H}_k\right)^{\mathrm{T}}+ \\ & \quad \boldsymbol{K}_k \boldsymbol{D}_{\eta_k} \boldsymbol{K}_k^{\mathrm{T}} \end{aligned} $

(4)

式中，K_k和 V_k分别为增益矩阵和新息向量，D_ηk为观测向量对应的观测噪声方差-协方差矩阵。由此可见，滤波估计结果$\hat{\boldsymbol{X}}_k $与状态预测值$ \hat{\boldsymbol{X}}_{k, k-1}$和当前历元的观测值Z_k均相关，并受过程噪声D_wk和观测噪声D_ηk的共同调节：设定的过程噪声越小，当前历元的预报值对滤波估计值的贡献越大，滤波估计结果更依赖于预测值而非观测值。监测体未发生形变时，监测站三维位置在前后历元间保持不变或变化较小，此时可对位置参数设置一个较小的过程噪声，抑制观测噪声的影响，使滤波结果序列更加平滑。然而，当监测体变形发展到一定时期，变形速率呈现出不断加速增长的趋势，此时监测站的三维位置在前后历元间发生显著变化，对应的滤波位置状态预测值与实际状态将不再相符，甚至导致模糊度无法或错误固定，得到错误的解算结果，无法实现对监测体显著形变的实时准确探测。因此，如何获取监测体位移信息，从而合理设置过程噪声的大小，是基于Kalman滤波实现高精度地质灾害形变监测的关键。

2 基于TDCP历元间位置变化量的监测体快速变形识别方法

TDCP通过相邻历元间相位观测值差分，可避开整周模糊度求解问题，同时极大削弱或消除硬件延迟、电离层与对流层延迟误差等在时间域上变化缓慢的误差影响，通过单台接收机即可获取载体高精度的位置变化，精度可达到mm级^[9]。本文将通过TDCP方法获取监测体历元间位置变化量，识别监测体发生的历元间显著位移，进而对Kalman滤波的过程噪声进行调节，使滤波结果反映监测体实际运动状态，实现对监测体快速形变的实时准确探测。

2.1 TDCP位移解算

相位观测值未发生周跳时，历元间相位差分观测方程为：

$ \lambda_f \Delta \varphi_f^i=\Delta \rho^i-c \Delta \mathrm{d} t^i+c \Delta \mathrm{d} t_R+\Delta \xi_f^i $

(5)

式中，Δ为历元间差分算子，Δρⁱ为卫星i与接收机距离的历元间变化量，cΔdtⁱ和cΔdt_R分别为以距离表示的卫星i及接收机R钟差历元变化量，Δξ_fⁱ为观测噪声项。无周跳情形下，通过历元间差分可消除整周模糊度。同时，前后历元间电离层延迟、对流层延迟等误差变化较小，通过历元间差分也可基本消除其影响。

假设r_k^R、r_kⁱ分别为接收机R和卫星i在k历元的三维位置矢量，k历元及k－1历元卫星与接收机连线方向的方向余弦矢量为e_kⁱ、e_k－1ⁱ，则上述观测方程可表示为：

$ \begin{gathered} \lambda_f \Delta \varphi_f^i-\Delta D^i+\Delta g^i+c \Delta \mathrm{~d} t^i= \\ -\boldsymbol{e}_k^i \Delta \boldsymbol{r}^R+c \Delta \mathrm{~d} t_R+\Delta \xi_f^i \end{gathered} $

(6)

其中，

$ \begin{aligned} \Delta \boldsymbol{r}^R & =\boldsymbol{r}_k^R-\boldsymbol{r}_{k-1}^R \\ \Delta D^i & =\left(\boldsymbol{e}_k^i \boldsymbol{r}_k^i-\boldsymbol{e}_{k-1}^i \boldsymbol{r}_{k-1}^i\right) \\ \Delta g^i & =\left(\boldsymbol{e}_k^i \boldsymbol{r}_{k-1}^R-\boldsymbol{e}_{k-1}^i \boldsymbol{r}_{k-1}^R\right) \end{aligned} $

(7)

式中，TDCP待估参数为Δr^R=[dx dy dz]^T和cΔdt_R，分别为前后历元间监测体三维坐标变化量和接收机钟差变化量。当观测到同一GNSS系统的卫星数不小于4时，采用最小二乘估计方法即可解算得到监测体的高精度三维位置变化量Δr^R[5]。

2.2 基于TDCP历元间位移的快速变形识别及滤波过程噪声调节

本文通过引入TDCP历元间位移量反映监测体在前后历元间的位置变化，用于探测识别监测体是否发生显著变形，进而根据位移量的大小调整Kalman滤波状态预测过程噪声，使后者与位置状态的历元间变化保持一致，从而使Kalman滤波结果与监测站实际位置状态保持匹配，实现对监测体快速变形的准确反映。该方法主要包括TDCP形变探测和滤波过程噪声调节2个部分。

1) TDCP形变探测。首先基于TDCP模型求解监测体三维坐标历元间变化量Δr^R。为提高TDCP位移解算的精度和可靠性，采用验前观测值一致性检验和验后IGGⅢ抗差迭代相结合的质量控制方案。然后通过TDCP位移量判断监测体是否发生显著变形：位移量‖Δr^R‖>C时视为发生显著位移，否则视为TDCP结果噪声。阈值C主要与TDCP解算结果精度相关，本文取C=0.01 m。

2) 滤波过程噪声调节。在站心东北天(ENU)坐标系方向进行过程噪声调节。记当前历元TDCP探测到的E、N、U方向的形变量为Δr_ENU^R = RΔr^R=[Δr_E^R Δr_N^R Δr_U^R]，单位为m，其中R为地心地固坐标(XYZ)转换为ENU坐标的旋转矩阵。ENU方向的位置过程噪声谱密度为σ_wk，单位为m·s^-1/2。ENU方向滤波过程噪声调节因子α的计算方法为：

$ \begin{aligned} & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\alpha_d= \\ & \left\{\begin{array}{l} \Delta \boldsymbol{r}_d^R /\left(\sigma_{w_k, d} \sqrt{t}\right), \left\|\Delta \boldsymbol{r}^R\right\|>C \\ 1, \left\|\Delta \boldsymbol{r}^R\right\| \leqslant C \end{array}(d=E, N, U)\right. \end{aligned} $

(8)

式中，t为历元时间间隔，单位为s。

该历元经过调节后的位置状态过程噪声方差-协方差矩阵为：

$ \begin{gathered} \boldsymbol{D}_{w_k}=\boldsymbol{R}^{-1} \boldsymbol{D}_{w_k, E N U}\left(\boldsymbol{R}^{-1}\right)^{\mathrm{T}}= \\ \boldsymbol{R}^{-1} \operatorname{diag}\left(\left(\alpha_d \sigma_{w_k, d}\right)^2 t\right)\left(\boldsymbol{R}^{-1}\right)^{\mathrm{T}}(d=E, N, U) \end{gathered} $

(9)

此时，Kalman滤波状态过程噪声量级与形变量一致，放松了对状态的约束，避免因状态约束不合理而导致滤波结果无法反映真实变形的问题。

基于TDCP历元间位移的快速变形识别及滤波过程噪声调节的主要流程如图 1所示。

3 实验验证与结果分析

采用精密微动螺旋位移平台和精密电控滑轨平台2组仿真实验，验证本文所提方法的可行性和有效性。位移平台实验在武汉大学测绘学院附2教学楼楼顶(3楼)进行，观测设备与环境如图 2(a)所示，电控滑轨实验在湖北高通空间技术有限公司办公楼楼顶(6楼)进行，观测设备和环境如图 2(b)所示。

位移平台实验采用安装有精密微动螺旋的三轴位移平台，平台设置3个移动方向相互正交的滑块，可通过调节其上的微动螺旋，带动GNSS天线在对应方向移动。位移平台的最小刻度为1 mm，最大范围为5 cm；监测站采用GNSS地质灾害监测设备JC-602，基准站使用Trimble Alloy测量型接收机及HX-CGX606A天线。实验中，将位移平台的移动轴向分别指向当地东(E)、北(N)和天(U)方向，从而使其位移在E、N、U方向上独立发生。

电控滑轨实验使用安装有精密电控滑块的线性位移滑轨(最大范围为55 cm)，通过设置滑块位移时间及速度，使GNSS天线在滑轨方向上运动，以模拟行程较大、速度较快的位移。实验中，监测站和基准站均采用GNSS表面位移监测设备DM-GNSS A300。

2组实验的详细信息如表 1所示，模拟位移信息分别如表 2和3所示，表中时间单位为GPS周秒。位移平台实验中，操作位移平台在6个预定时刻分别单独进行E、N、U方向的位移，位移时间小于1 s，位移距离大于1 cm，以验证本文方法探测快速位移的可行性与有效性。电控滑轨实验中，设置滑轨在预定时刻沿其安置方向进行快速位移，总位移量大于500 mm，以模拟在地质灾害监测中监测体处于临滑阶段时发生的持续快速变形，进而验证本文方法在模拟监测场景中的效果。

利用武汉大学测绘学院研制的GNSS高精度实时变形监测软件，结合本文基于TDCP的快速变形识别及过程噪声调节方案(以下简称方法A)，对上述快速变形模拟测试数据开展实时处理，并将处理结果与不含TDCP变形识别的固定过程噪声Kalman滤波定位方法(以下简称方法B)解算结果进行对比，2种方法使用的解算配置如表 4所示。为了保证分析结果的可靠性，2种解算方法在数据输入、基本数据处理策略上保持一致。实验中，将2种方法解算得到的监测体位移量与根据模拟位移信息得到的参考真值进行对比，以评估不同方法识别快速变形的性能。

3.1 位移平台实验结果分析

位移平台实验中，基于TDCP解算得到的监测站E、N、U方向的历元间位移量时间序列如图 3所示。由图可知，位移平台在某个方向发生位移时，对应方向上的TDCP位移值都显著大于其静止状态下的正常波动范围，表现为离群值。因此，TDCP可及时准确地探测监测站发生的快速变形。

方法A和方法B解算得到的位移平台E、N、U方向位移序列如图 4所示。由图可知，在模拟发生位移的历元时，方法A的滤波解算结果在短时间内即可识别到相应的位移量，大小与参考真值基本一致；方法B的滤波解算结果虽然在历元间波动较小，但未能准确识别发生的快速位移，位移序列经过较长时间的收敛后才与实际状态趋于一致。

2种方法的位移解算结果相较于参考位移量的误差序列如图 5所示。由图可知，整个测试过程中，方法A的E、N方向误差总体小于5 mm、U方向误差小于1 cm，表明该方法的位移解算结果始终与预设模拟位移保持一致。方法B的位移解算结果则在发生位移后的大部分时刻与真实值产生较大偏差。2种方法的位移序列RMS如表 5所示，其中方法A的E、N、U方向位移解算误差RMS可达1.7 mm、0.5 mm、3.0 mm，相较于方法B分别提升约72.2%、90.7%、78.5%。

3.2 电控滑轨实验结果分析

电控滑轨平台仿真实验中，基于TDCP解算得到的监测站E、N、U方向的历元间位移量时间序列如图 6所示。由图可知，TDCP可及时准确地探测监测站发生的持续多个历元的快速位移。

基于2种方法得到的沿滑轨方向的监测站位移序列及其误差序列分别如图 7和8所示。由图 7可知，对于模拟发生位移的连续8个历元，方法A的滤波解算结果均可准确识别滑轨上监测站的显著位移，最终准确反映出约50 cm的位移量，与参考值基本一致；方法B的滤波解算结果则出现较大偏差，没有准确反映出滑轨发生的真实位移。由图 8可知，整个解算过程中，方法A的位移序列与参考值的误差均小于1 cm，表明该方法位移解算结果与预设模拟位移基本一致；方法B的位移解算结果则在发生位移的时刻之后与模拟值产生较大偏差，且较长时间(图中约为600 s)无法收敛，未能反映出真实位移。

2种方法的位移序列误差RMS如表 6所示。由表可知，发生位移前后，方法A的位移解算误差RMS约为1.55~1.75 mm，与发生位移前方法B的位移解算误差RMS基本相当。

由上述结果可知，本文提出的基于TDCP的快速变形识别方法可以在保持位移监测精度的同时，准确识别和反映监测体的快速变形。

4 结语

本文提出一种基于TDCP历元间位置变化量的监测体快速变形识别方法。该方法通过TDCP估计监测体前后历元间的三维位置变化量，进而识别和判断监测体发生的显著位移。当监测体因自身形变导致原有的过程噪声难以匹配实际变形状态时，对Kalman滤波的系统过程噪声进行调节，使滤波状态估计结果与监测体实际变形状态保持一致，以反映监测体的快速变形。采用微动螺旋位移平台和电控滑轨平台仿真实验数据对本文方法的可行性与有效性进行分析，并与常见的固定过程噪声的Kalman滤波方案进行对比。结果表明，本文方法的位移序列相较于参考真值的误差在水平方向小于5mm、在高程方向小于1 cm，可实时、准确地反映监测体历元间大于1 cm的快速变形。同时，该方法可有效抑制观测噪声、实现mm级的监测精度。整体而言，本文方法可在有效抑制观测噪声、保证监测精度的同时，实时、准确地识别和反映监测体的显著变形，进而为地质灾害监测提供准确可靠的监测体变形信息，在GNSS地质灾害监测预警领域具有重要的应用潜力。

需要指出的是，由于监测体受到其内部结构和外部因素的共同作用，滑坡形变的特征较为复杂，与GNSS监测站所处的观测条件也不尽相同，本文方法仅能有效识别和反映监测体相邻历元间cm级以上的快速变形。对于地质灾害中另一类常见的慢速蠕动形变，则需要进一步研究更加有效、可靠的方法。此外，本文仅基于较为开阔场景的实验数据对所提方法的有效性进行验证，未来还需进一步研究适用于复杂监测场景、稳健可靠的变形识别和滤波过程噪声自适应调节方法。