本文提出的基于VRS的GNSS滑坡变形监测方法的数据处理流程如图 1所示，可以分为两个部分：

图  1  数据处理流程

Figure 1.  Workflow of Data Processing

1）服务端VRS观测数据生成。服务端利用CORS基准站观测数据，首先固定双差宽巷（wide-lane，WL）和窄巷（narrow-lane，NL）整数模糊度，计算高精度站星几何距离；然后解算得到双差对流层和电离层延迟，建立区域对流层、电离层延迟模型；最后内插生成VRS观测值数据。

2）监测端基线解算及网平差处理。首先利用监测站坐标组成Delaunay三角网；然后逐基线进行解算，得到每条基线的向量；其次依据三角形闭合差为零的准则，对基线向量解算结果进行检验；最后将检验通过的基线解算结果进行网平差处理，得到最终的监测点坐标向量。

与网络RTK数据处理方法和流程类似，服务端先利用CORS基准站观测数据固定双差模糊度，考虑基准站间距达到20~70 km，一般采用无电离层组合观测值消除电离层延迟的影响^[15-16]，观测方程为：

式中，$\nabla \mathrm{\Delta }$表示双差算子；$i$、$j$为卫星编号，$a\mathrm{、}b$为测站编号；$\nabla \mathrm{\Delta }{P}_{a,b,\mathrm{I}\mathrm{F}}^{i,j}\mathrm{、}\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{I}\mathrm{F}}^{i,j}$分别为伪距和载波相位无电离层组合观测值；$\nabla \mathrm{\Delta }{\rho }_{a,b}^{i,j}$为双差站星距离；$m_{a,b}^{i,j}$为对流层投影系数；$\mathrm{\Delta }{T}_{a,b}^{}$为站间单差天顶对流层延迟；${\lambda }_{\mathrm{I}\mathrm{F}}$为无电离层观测值波长；$\nabla \mathrm{\Delta }{N}_{a,b,\mathrm{I}\mathrm{F}}^{i,j}$表示双差无电离层组合模糊度；$\epsilon (\nabla \mathrm{\Delta }{P}_{a,b,\mathrm{I}\mathrm{F}}^{i,j})$、$\epsilon (\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{I}\mathrm{F}}^{i,j})$分别为伪距和载波观测值的噪声。

由于无电离层组合观测值不具有整数特性，故将其分解为WL和NL模糊度进行固定^[17-19]，其中WL模糊度可采用MW（Melbourne-Wübbena）组合观测值平滑取整进行固定^[19]，固定方程为：

式中，$\nabla \mathrm{\Delta }{N}_{a,b,\mathrm{W}\mathrm{L}}^{i,j}\mathrm{、}\nabla \mathrm{\Delta }{N}_{a,b,\mathrm{N}\mathrm{L}}^{i,j}$分别表示双差WL和NL模糊度；$f_{\mathrm{1}}\mathrm{、}{f}_{\mathrm{2}}$分别表示L1和L2观测值的频率；int表示取整操作；$\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{M}\mathrm{W}}^{i,j}$表示双差MW组合观测值；$<\cdot >$表示历元间平滑滤波。

模糊度固定后，首先建立区域大气误差模型，然后根据待监测区域范围和监测站的概略坐标，得到需要虚拟的参考站数量及其坐标，最后利用线性内插模型（linear interpolation model，LIM）^[20-21]插值得到大气误差改正模型，同时虚拟生成对应的VRS观测值数据^[22]。

监测端数据处理时，将服务端数据处理生成的VRS观测值与监测站观测数据组成双差观测值。VRS虽然与待监测点相距较近，一般仅几十到几百米，但考虑VRS观测值实际由距离待监测点几十千米的CORS站计算得到，因此，双差观测值仍需要考虑对流层和电离层残余误差的影响^[23-24]，因此本文采用L1、L2双频非组合模型，基线解算观测方程为：

式中，$\nabla \mathrm{\Delta }{P}_{a,b,\mathrm{1}}^{i,j}\mathrm{、}\nabla \mathrm{\Delta }{P}_{a,b,\mathrm{2}}^{i,j}$分别表示L1和L2双差伪距观测值；$\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{1}}^{i,j}\mathrm{、}\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{2}}^{i,j}$分别表示L1和L2双差载波相位观测值；$\mu$表示L2相对L1观测值的电离层比例系数；$\nabla \mathrm{\Delta }{I}_{a,b}^{i,j}$表示L1观测值上双差电离层延迟；${\lambda }_{\mathrm{1}}\mathrm{、}{\lambda }_{\mathrm{2}}$分别表示L1和L2观测值的波长；$\nabla \mathrm{\Delta }{N}_{a,b,\mathrm{1}}^{i,j}\mathrm{、}\nabla \mathrm{\Delta }{N}_{a,b,\mathrm{2}}^{i,j}$分别表示L1和L2观测值的双差模糊度；$\epsilon (\nabla \mathrm{\Delta }{P}_{a,b,\mathrm{1}}^{i,j})$、$\epsilon (\nabla \mathrm{\Delta }{P}_{a,b,\mathrm{2}}^{i,j})$、$\epsilon (\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{1}}^{i,j})\mathrm{、}\epsilon (\nabla \mathrm{\Delta }{L}_{a,b,\mathrm{2}}^{i,j})$表示观测值噪声。

为方便数据处理，本文采用静态逐次滤波解算式（3）中的待估计参数，包括基线向量、站间单差天顶对流层延迟、双差斜路径电离层延迟以及双差模糊度参数，其中基线向量在一次解算中认为保持不变，站间单差对流层延迟采用随机游走模型，而双差电离层斜路径延迟参数认为历元间不相关。站间单差对流层延迟和双差电离层斜路径延迟可利用服务端提供的对流层和电离层模型的精度作为先验信息进行约束。基线解算后，为进一步提高数据处理精度，本文将VRS与各监测站组成Delaunay三角网，添加三角网约束信息进行网平差，约束方程为：

式中，$\nabla \mathrm{\Delta }{\mathit{X}}_{\mathrm{V}\mathrm{R}\mathrm{S},a}$、$\nabla \mathrm{\Delta }{\mathit{X}}_{b,\mathrm{V}\mathrm{R}\mathrm{S}}$分别表示VRS与Delaunay三角网另外两个监测站a、b的双差基线向量；$\nabla \mathrm{\Delta }{\mathit{X}}_{a,b}$表示监测站a、b的双差基线向量。

屏边县地处云南省南部、红河哈尼族彝族自治州东南部，受东南海洋暖湿气流的影响，境内湿润多雨，地势北高南低，由北向南倾斜，由于河流的切割，地貌形成了“四河三山六面坡”的结构，高山横亘连绵，V型谷较多，地形极其复杂，是全州地质灾害隐患点较多的县（市）之一。2020年全县共有地质灾害隐患点139处，包括滑坡点94处，崩塌点30处，泥石流沟14条，地面塌陷1处，受地质灾害隐患威胁群众2 762户16 046人，潜在经济损失达5.2亿元。

本文选取其中30处滑坡灾害隐患点进行实验分析，如按常规的一点一基站方式，共需布设30个基站，若要实现全县所有隐患点监测则需要建设139个基站，将造成极大的成本压力。而红河哈尼族彝族自治州已建成了全域126个CORS基准站，平均站间距约20 km，部分CORS站点分布如图 2所示。充分利用已有的CORS基准站生成VRS代替实体基准站，进行滑坡形变监测应用，可大大降低监测成本。此外，VRS利用多基准站虚拟生成，消除了基准点不稳定导致基准变化的影响。

图  2  部分CORS站点分布图

Figure 2.  Distribution of Some CORS Stations

电离层延迟和对流层延迟模型精度是影响监测端数据处理精度的主要误差，本文选取基线长度分别为20 km、40 km和80 km以内的观测数据（分别对应图 2中红、黄、蓝三角网），分析电离层延迟和对流层延迟模型精度。利用多个基准站数据固定模糊度后，分离得到区域电离层和对流层延迟模型，采用LIM内插得到虚拟站与基准站之间双差电离层和对流层延迟，同时利用监测站实测数据计算某一卫星对的双差电离层和对流层延迟，比较两者的差异以分析电离层和对流层延迟模型的精度。

图 3和图 4分别为某一卫星对双差电离层和对流层模型的内插值和计算值的时间序列。由图 3和图 4可以看出，对不同的基线长度，经过网络RTK内插模型计算的双差电离层和对流层延迟与实际计算的双差电离层和对流层延迟具有较好的一致性，可以较好地改正电离层和对流层延迟的影响。

图  3  不同基线长度双差电离层建模精度分析

Figure 3.  Accuracy Analysis of Double-Difference Ionospheric Modeling with Different Baseline Lengths

图  4  不同基线长度双差对流层建模精度分析

Figure 4.  Accuracy Analysis of Double-Difference Tropospheric Modeling with Different Baseline Lengths

图 5和图 6分别为所有卫星对双差电离层和对流层延迟的内插值和计算值的差异概率分布。由图 5和图 6可以看出，双差电离层和对流层延迟内插值和计算值的差异均呈现正态分布的特性，且随着基线长度的增加，内插值和计算值差异逐渐增大，说明随着站间距离的增加，电离层和对流层的内插模型的精度逐渐降低。当站间距小于40 km时，虽然双差电离层和对流层延迟模型的标准差（standard deviation，STD）仅能达到厘米级，但是差异的均值小于1 cm。而滑坡形变监测数据处理一般采用事后解算模式，通过多历元观测数据可以有效削弱内插模型厘米级噪声的影响，实现高精度数据处理。

图  5  不同基线长度双差电离层建模精度统计

Figure 5.  Accuracy Statistics of Double-Difference Ionospheric Modeling with Different Baseline Lengths

图  6  不同基线长度双差对流层建模精度统计

Figure 6.  Accuracy Statistics of Double-Difference Tropospheric Modeling with Different Baseline Lengths

考虑双差电离层和对流层延迟模型的内插精度影响，本文将CORS基准站站间距设置为小于40 km，服务端数据处理对应生成30个VRS观测数据，与90个监测站数据进行事后模拟，通过基线解算及网平差处理获得监测站的三维坐标变化。根据30个测站一周的观测数据，利用不同观测时长（2 h、4 h、6 h、8 h、10 h、12 h）数据的解算结果，统计N、E、U方向的均方根（root mean square，RMS）变化，如图 7所示。

图  7  不同观测时长数据的解算结果精度

Figure 7.  Accuracy of Data Processing Results with Different Observation Duration

由图 7可以看出，数据处理解算结果的精度与观测时段长短密切相关，2 h时长观测数据的解算结果仅能达到厘米级精度，在95%置信区间内，U方向的RMS可以达到25 mm；随着观测时长增加，解算精度逐渐提高，当观测时长达到8 h时，在95%置信区间内，N、E方向的RMS控制在5 mm以内、U方向的RMS控制在10 mm以内。

图 8为网平差前后所有基线的RMS均值。由图 8可以看出，经过网平差处理，基线解算结果在N、E、U方向的RMS均值分别从3.2 mm、3.8 mm和5.8 mm降至2.9 mm、3.0 mm和5.1 mm，精度分别提高了9%、14%和12%。选取某一监测站连续15天的8 h观测数据，采用本文方法和传统方法处理结果的时间序列如图 9所示。由图 9可以看出，传统形变监测方法由于基准站与监测站相距较近，其解算结果在N、E、U方向的精度分别达到0.51 mm、0.54 mm和1.26 mm；本文方法的解算结果在N、E、U方向的精度分别为3 mm、3 mm和6 mm，虽然解算精度不及传统方法，但也能够满足毫米级的形变监测应用需求。

图  8  网平差前后各基线精度统计值

Figure 8.  Statistics of Average RMS of All Baselines Before and After Network Adjustment

图  9  监测站数据处理结果时间序列

Figure 9.  Time Series of Data Processing Results of Monitoring Station

本文提出了一种基于VRS的GNSS变形监测方法，直接利用现有CORS观测数据虚拟生成VRS，代替实体参考站实现高精度变形监测。利用中国云南省屏边县约30处滑坡监测点进行实验分析，结果表明，当CORS基准站间距小于40 km时，双差对流层和电离层延迟模型内插精度优于10 mm，利用观测时长大于8 h的数据，可实现平面精度优于5 mm、高程精度优于10 mm的形变监测。