高精度时空信息崩滑坡灾害监测技术

本文内容来源于《测绘通报》2022年第11期，审图号：GS京（2022）1346号

高精度时空信息崩滑坡灾害监测技术

何伟^1,3, 张伦宁^2,3, 曾露¹, 赵清¹, 余国新², 竺科英², 李柯欣²

1. 复旦大学工程与应用技术学院, 上海 200433;

2. 新疆农业大学经济管理学院, 新疆 乌鲁木齐 830052;

3. 上海华测导航技术股份有限公司, 上海 201702

基金项目：国家重点研发计划(2018YFC1505103);上海市“科技创新行动计划”(19511100200;20511103100);应急管理部国重开放基金(NDRCCPICC202005)

关键词：高精度时空信息, 形变监测, 中长基线, 精密定位, 地质灾害

引文格式：何伟, 张伦宁, 曾露, 等. 高精度时空信息崩滑坡灾害监测技术[J]. 测绘通报, 2022(11): 1-7. DOI: 10.13474/j.cnki.11-2246.2022.0316.

摘要

摘要 ：我国部分地区地质灾害爆发频率高,多年来崩滑坡、泥石流等高危险性地质灾害严重威胁我国人民财产与生命安全。因此,地质灾害及时预警一直是我国在应急方面亟待解决的难题。随着北斗、GPS等全球卫星导航系统的建设组网,GNSS技术逐渐成为当前提供时空位置信息的重要手段,具有数据实时性、设备低功耗、搭建便捷化等特性,成为解决地质灾害监测与预警的理想技术手段。本文利用BDS+GPS+GLONASS时空信息数据,构建了面向地质灾害监测的中长基线毫米级精度的崩滑坡监测解算方法,确定了双差观测的载波相位解算模型。最后在重庆新浦区域展开试验应用,实现了基于高精度时空信息的毫米级滑坡灾害的监测与预警。为预计该区域的地质灾害提供了数据支持,为人民的生产生活提供了技术保障。

正文

我国的地质构造与地震活动强烈、地质情况复杂且降水在不同时空分布不均, 导致地质灾害隐患多、分布广、灾情重、风险大，进而使我国成为全球范围内地质灾害严重且处于其威胁下人口较多的国家之一^[1]。截至目前, 在我国已查明的存在地质灾害隐患的地点约有32万处之多，威胁1300余万人的生命安全和5000多亿元的财产安全。受地质环境及气象条件的控制, 尤其是在川、滇、赣、湘、陇、渝、新、藏、晋、陕、桂、贵等12个中西部地区的省(自治区、直辖市), 地质灾害隐患较集中, 约超过全国隐患总数的80%。

GNSS时空信息系统在20世纪80年代投入民用，其代表为美国GPS, 现已广泛应用在导航、定位等相关领域^[2], 特别是GNSS控制测量技术在测量领域的重要性难以估量。随着GNSS高精度实时精密定位参考站等研究的应用, 基于GNSS高精度时空信息的自动化监测系统可探测出毫米级的位移变化, 逐渐成为地质形变实时监测的有效手段^[3-4]。

本文利用BDS+GPS+GLONASS时空信息数据，构建了面向地质灾害监测的中长基线毫米级精度的崩滑坡监测解算方法，确定了双差观测的载波相位解算模型。最后在重庆新浦区域展开试验应用，实现了基于高精度时空信息的毫米级滑坡灾害的监测与预警。

1 GNSS高精度实时精密定位系统设计1.1 总体设计

GNSS区域高精度增强服务目前以连续运行参考站系统为主^[5]。该系统是在数据通信网络的基础上，动态实时、连续快速、精密准确地在空间范畴中获得所需的地理特征及具体数据的信息现代化基础设施体系^[6]。

GNSS高精度实时精密定位系统的总体架构如图 1所示。包括：研究系统的总体架构与系统组成，系统的工作模式；各子系统之间的参数配置与接口、数据流程及格式、协议与标准；各子系统任务与指标分配；IGS跟踪站信息处理与区域增强网信息服务模式，IGS跟踪站信息融合方案。

图 1 GNSS高精度实时精密定位系统总体架构

基于新一代网络实时差分定位技术的区域增强信息反演是建立在非差处理模式的基础上，即在计算模型参数时必须对基准站的非差整周模糊度进行固定。这项操作的原理与双差模糊度不同, 由于非差模糊度接收机自身存在偏差且卫星相位极易出现偏差^[7-8], 导致其只具备非整数性质。因此, 拟对如下内容进行详细研究: ①分析模糊度组合特征(具体存在于诸如宽巷组合、几何无关组合、无电离层组合等双频非差模糊度线性组合中), 并结合参数化的电离层模型与对流层模型，确定L1/L2非差模糊度实数解解算策略；②剖析相位偏差中的时空层面特点，精准确立关于实时预报层面的模型，重新确立L1/L2非差模糊度整周特性^[9]并固定。

(1) 硬件延迟及初始相位偏差(UPD)确定技术。分析影响模糊度固定正确性、可靠性相关因素，重点研究未检验相位偏差(UPD)对精密单点定位模糊度解算的影响及其时空变化规律，找出符合实际情况的相位偏差求解模型，构建适合多种偏差影响的分离方法与提取模型，研究实时估计方法与预报模型，从而为模糊度快速、准确的收敛提供支持。

(2) 基于CORS的非差观测值残差OMC反演技术。分析电离层误差、对流层误差等大气因素对精密单点定位的影响，构建适当的数学模型反演CORS基站上空的大气误差，重点研究大气误差的区域性差异对单点定位结果精度的影响，构建大范围高精度非差观测值残差OMC改正模型。

1.2 BDS+GLONASS+GPS三系统联合解算研究

(1) 将BDS+GLONASS+GPS数据解算至统一的时空基准中，结合卫星导航发展成熟度，选择GPS时空基准为统一基准^[10]。

(2) 观测值线性化与误差模型改正，主要有卫星钟差、测站坐标潮汐、星地几何距离、天线相位中心偏差、对流层延迟、地球自转误差等层面的相应改正。最主要的误差分别是由于电离层和对流层的延迟误差^[11]对整个系统造成影响而导致的。只用模型改正的方法难以全面解决电离层延迟误差的问题^[12]，还应结合电离层无关组合。具体改正方法如下: ①短基线(一般小于10 km)，采用Klobuchar模型改正；②中长基线(大于10 km)，本文采用线性组合法消除其影响，公式为

(1)

式中，ϕ_c为无电离层组合观测值；ϕ₁为L1信号或B1信号的载波相位观测值；ϕ₂为L2信号或B2信号的载波相位观测值；f_B1、f_B2分别为L1/B1及L2/B2信号的频率。

针对对流层延迟误差的修正(干分量部分)，在参数估计的基础上进行模型改正，再利用Sasstamoinen模型理论, 可使精度达厘米级。但若基线长度较长，无法用模型改正, 则需要采用参数估计的方法, 即将对流层天顶湿延迟线性化, 确定为参数后再进行估计运作, 公式为

(2)

式中，Δ∇T_pq^ij是双差对流层延迟；Δ∇T_{pq, dry}^ij是双差对流层干延迟；Δ∇T_{pq, wet}^ij是双差对流层湿延迟；∇MF(θ_p^ij)、∇MF(θ_q^ij)分别是映射函数；ZD_{p, wet}、ZD_{q, wet}分别是p、q两站干湿分量(天顶方向)。

其中，Δ∇T_{pq, dry}^ij使用模型直接改正，由于天顶对流层湿延迟ZD_{p, wet}、ZD_{q, wet}演变缓慢, 可采用参数估计方法对对流层湿延迟误差进行改正，并将此参数结合其他参数一起进行估计。

(3) 对初始三差模型进行求解计算, 求得基本的基线向量, 反复对周跳进行探测与修复, 公式为

(3)

式中，Δ∇ϕ_AB^ir(t₁, t₂)为t₁与t₂间的三差观测值；Δ∇R_AB^ir(t₁, t₂)为t₁与t₂间的三差站间距；Δ∇N_AB^ir(t₁)为t₁时刻的双差模糊度；Δ∇N_AB^ir(t₂)为t₂时刻的双差模糊度；ε_AB^ir(t₁, t₂)为剩余残差项；ir为卫星对；A、B为基线站点。

令μ=[Δ∇N_AB^ir(t₁)-Δ∇N_AB^ir(t₂)]+ε_AB^ir(t₁, t₂)，若历元间不发生周跳，即Δ∇N_AB^ir(t₁)-Δ∇N_AB^ir(t₂)=0，则μ为本方程常数项, 此时，μ仅存在相邻历元间大气传播误差项与观测噪声项的影响。μ一般小于0.1 cycles。若μ大于1 cycles时，则t₁与t₂之间必存在周跳，Δ∇ϕ_AB^ir(t₁)与Δ∇ϕ_AB^ir(t₂)之间的周跳数为ROUND(μ)。其中，ROUND(μ)为常数项，并结合实际应用需求四舍五入。

为了提高模糊度固定的成功率，在三差探测的基础上采用以下策略进行修复：即当|μ-ROUND(μ)| < 0.25周期时，同时进行周跳修复工作; 否则将该卫星视为新升卫星, 重估其模糊度。此外，采用上述策略, 无论对于单频观测数据还是多频观测数据均适用。

2 面向监测系统的中长基线毫米级精度基线解算

基于中长基线毫米级精度基线解算策略，改进了北斗毫米级精度形变监测算法，为解决单历元北斗高精度形变监测技术提供了可能。

2.1 数据组合策略

针对不同研究目的，GNSS多种时空信息观测值与之对应有不同的组合策略。此外, 观测值于不同主体(如历元间、星间或站间)间求差。在不同的组合策略中，双差观测值是在高精度定位中常见的数据组合策略，它可使接收机钟差消失，卫星钟差参数误差完全清零, 而且对轨道、对流层、电离层等大多数误差具有减弱作用。结合本文要求, 确定采用具有高精度特性的双差组合。

2.2 误差改正

数据解算过程中的误差改正如图 2所示。

图 2 顾及的误差改正项

2.3 周跳探测与修复

周跳探测与修复是北斗数据处理中的重要步骤，只有经过周跳修复后的相位观测值才具有高精度距离观测值的意义。笔者所在的课题组根据前述的研究成果，通过以下3个步骤进行周跳探测与标记。

2.3.1 探测和标定周跳(M-W组合法)

利用GNSS观测值中载波、伪距观测值的表达式经推导可得

(4)

式中，φ₁、φ₂、P₁、P₂分别为GNSS系统的双频载波相位、伪距观测值；f₁、f₂为载波频率；N₁、N₂为双频上的模糊度参数。

由式(4)可知，利用相位、伪距观测值经过组合即可直接得到宽巷模糊度N₁-N₂，通过比较不同历元宽巷模糊度值即可实现周跳探测与标定。利用M-W组合法可消去电离层、测站、卫星等其他无关参数的影响, 使其仅处于观测噪声和多路径效应下, 而这两种影响完全通过多历元平滑法可以消除。因此，在大气折射、站坐标改正有误差出现时, 这一组合仍可以进行有效的周跳探测。

2.3.2 利用电离层残差组合(ionospheric residual)φ_IR探测与标定周跳

双频相位层面的不同观测值组合，可得到几何关系相应观测值(geometry free)，公式为

(5)

式中，φ_GF为无几何关系组合观测值；δ_I为以距离表示的L1频率上的电离层延迟改正值。利用前后历元上的无几何关系组合观测值做差即可得到电离层残差组合φ_IR，其仅与频率间、历元间的电离层残差之差相关，对于周跳在电离层变化较稳定时非常敏感，而且相位观测值的精度很高。因此，这一组合在周跳探测中具有重要作用，可以弥补M-W组合的不足。

2.3.3 利用LC观测值的单差(星间)残差观测值进行周跳探测与标定

将双频相位观测值组合得到无电离层延迟组合观测值LC, 通过LC在星间做差以消除接收机钟差的影响, 再利用已知的卫星、测站近似坐标及卫星钟差信息即可得到LC单差残差观测值。清除接收机钟差、电离层的影响，并减弱卫星-测站卫星钟差等无关参数的影响, 是周跳探测的有效手段。

2.4 参数估计方法

本文将L1或LC观测值利用法方程有效叠加，在此基础上对参数进行估计。

首先在已有法方程中加上参数约束, 然后逐个读取OMC值(改正过)和偏导数值(与之相对应), 利用双差映射矩阵法将其转变为双差观测值, 最后将相关参数重复输入。在此过程中需检查卫星星座变化状况, 找到映射矩阵再次形成的客观依据, 在法方程中输入模糊度数值, 再消除原有数值(参数消去恢复)。对所有观测值进行求和, 在基线长度和观测情况的实际结果上进行运算，以达到双差模糊度映射矩阵重新构成的目的。完成以上步骤即可求解法方程，可得出相应有效的参数估值，如图 3所示。

图 3 参数估计流程

2.5 模糊度固定策略

模糊度参数在理论层面需要的是整数, 但在模型中存在各种误差, 可能导致由法方程解算出的参数为非整数。为确保模型的精确性必须将相关参数值固定为整数值。

根据研究成果，采用Bootstrap模糊度固定方法，结合固定失败率的可靠性检验方法对模糊度进行固定。基于Bootstrap的模糊度固定流程如图 4所示。

图 4 模糊度固定流程

针对目前BDS快速发展的现状，本文研究了基于北斗卫星导航系统的高精度形变监测系统，并编制了相关软件。考虑目前北斗系统的精密星历产品仍处于试验阶段，因此软件采用了广播星历进行解算，可满足短基线的毫米级精度监测需求。

2.6 函数模型

双差模式消除了接收机、卫星钟差，削弱了其他各种误差，可以得到高精度的结果(尤其在短基线的情况下)，是形变监测的最佳选择。因此课题组选择双差观测值用于基于北斗卫星导航系统的高精度形变监测系统。BDS解算的双差观测方程为

(6)

式中，φ_jk^mn为经双差操作后的观测值(载波相位)；m、n为与双差一一对应的卫星号；j、k为两双差测站；φ_j^m、φ_jⁿ、φ_k^m、φ_kⁿ为相位观测值(单程载波层面)；X_j、X_k为信号两测站位置矢量(地固系中)；X^m、Xⁿ为信号发射时刻两卫星在地固系中的位置矢量；f_i为载波频率，其中f₁=1 561.098 MHz、f₂=1 207.14 MHz；c为光速；N_j^mN_jⁿ、N_k^m、N_kⁿ为相位观测值所对应的模糊度参数(单程载波层面)。在短基线情况下，对流层延迟、电离层延迟等改正项可认为全部被清除，在观测方程中不予考虑。

2.7 系统算法与实现

为实现使高精度卫星导航时空信息数据的形变达到可监测的程度, 应专门研发对其有针对性的形变监测软件, 具体处理流程如图 5所示。

图 5 数据处理流程

系统主要包括4部分：固定模糊度参数、分析时间序列、确定相关数据、形成与叠加法方程等，涵盖了数据输入到流出的全过程。

(1) 数据准备工作：一是确定卫星位置，具体方法为通过广播星数据转化后, 供后续计算内插使用; 二是初步观测周跳探测, 并使用TurboEdit法改进收集数据的准确性。

(2) 利用双差观测方程, 使模糊度与双差模糊度互为映射，并累加法方程。考虑系统的短基线、高精度要求, 解算时仅采用B1频率观测值。法方程的形成需考虑地球自转改正等误差。使用MGEX(multi-GNSS experiment)模型, 但尚未明晰接收机天线相位中心改正模型, 在系统中不予考虑。残差编辑解算法方程后得到浮点解, 并标记相应周跳和坏值, 再重构法方程, 重复此操作至残差编辑的周跳或坏值清零。

(3) 用Bootstrap+Decision Function方法固定模糊度。首先易知概率相关性，在实数模糊度及其方差-协方差矩阵法基础上进行求解计算；其次按概率大小将模糊度参数也如此排序, 确定最大概率值，并将与之相对应的模糊度参数固定取整；然后将其代入法方程中, 再将余下非整模糊度参数进行解算, 并把余下的模糊度参数重新排序，反复操作以上具体内容，直至无模糊度参数可以被固定；最后将以上整数模糊度参数代入法方程中求解运算，以得到固定解。

(4) 时间序列分析有两个运行层面，即频域与时域。在频域分析中，原始形变信号中以地球物理与气象因素为主要研究对象, 在使用高、低通滤波两方面从原始结果中提取真正形变数据。时域分析方法相对简明，对形变时间序列建模以达到预测未来的形变情况的目的。所有观测与解算信息均存储于数据库系统中。数据库系统通过ODBC技术与形变监测系统相连, 具有实时数据管理功能。

3 试验与应用3.1 监测站建设

本文选取了重庆奉节新浦滑坡段，建设了3个GNSS监测站，1个基准站。3个点位布设在一条监测剖面上，如图 6所示。

图 6 点位布置

3.2 设备监测情况

项目数据从2020年3月1日至2021年8月23日选取，设置数据实时采集，按每小时进行解算，整体监测情况如图 7所示。

图 7 整体监测

由图 7可以看出，该区域整体向北向滑动，滑动趋势范围内无密集村庄、农田等人类生产生活集聚区。

(1) X方向断面如图 8所示。

图 8 X方向断面

由图 8可看出，断面曲线X方向点位2变化最大，达443.7 mm。

(2) Y方向断面如图 9所示。

图 9 Y方向断面

由图 9可看出，断面曲线Y方向整体变化量不大，其中点位1变化最小，仅7.9 mm。

(3) H方向断面如图 10所示。

图 10 H方向断面

由图 10可以看出，断面曲线H方向点位1变化最明显达83.3 mm，点位2、3变化相对一致，分别为54 mm、43 mm。

(4) X-Y平面断面如图 11所示。

图 11 X-Y平面断面

由图 11可以看出，X-Y平面位移变化量点位2变化最为明显，达446.6 mm，点位3变化最小仅262.2 mm。

3.3 变化量汇总

本文设计的模型支持毫米级解算，并通过试验可以看出重庆新浦滑坡点2位置平面变化量最大，达446.6 mm, 点位1三维滑动量最大，达515.8 mm。其中，点位1由于建设地点较高，滑动H方向变化最大，下沉83.3 mm，而点位2平面方向上在X方向变化最大，为443.7 mm。目前在滑动趋势上未对较为集中的村庄、农田造成影响。

4 结语

随着高精度时空信息监测技术的不断发展与成熟，应用于地灾、地质形变监测将逐渐由准实时向实时性转变，对地质灾害预警的准确性也逐步提升。本文利用基于GNSS的可持续运行参考站数据，构建了面向监测领域的毫米级中长基线解算模型，并在重庆新浦展开了实地应用与验证。

经试验验证，模型具有的主要优点有：数据实时性优良、自动化程度高、监测结果输出的数据高度可靠，能够探测出滑坡毫米级的变化量，为管理部分提供了连续、稳定的监测数据，大大提高了作业效率，减少了人工劳动，为保证该区域村民的生产生活安全提供了数据保障。

作者简介

作者简介：何伟(1977-),男,博士生,主要从事高精度GNSS设备研发和应用研究的工作。E-mail:[e-mail]

通信作者：张伦宁。E-mail:[e-mail]

初审：纪银晓复审：宋启凡
终审：金   君

往期推荐

资讯