一种黄土滑坡灾害综合监测方法及系统

1.本发明涉及地质灾害监测技术领域，具体为一种黄土滑坡灾害综合监测方法及系统。


背景技术：

2.黄土高原是我国三大地质灾害最为严重的地区之一，水敏性是黄土最显著的特征，连续性强降雨量作用下雨水入渗到黄土内部致使土体饱和度迅速增高，土壤中水分的变化影响了土体的抗剪强度，从而诱发黄土浅层滑坡。有效开展此类滑坡的形变及诱因监测工作，对我国防灾减灾工作具有重要的意义。
3.目前，全球导航卫星系统gnss作为能够直接获取实时地表三维矢量变形的技术手段，被广泛应用于滑坡灾害监测领域。同时，全球导航卫星系统gnss还可以提供源源不断的l波段微波信号，以获取测站区域的大气水汽含量和土壤含水率等地表环境信息。
4.针对在以往滑坡灾害监测研究中，全球导航卫星系统gnss也仅仅只是被用于提供三维形变信息，而其能够提供监测站周围环境信息的能力却几乎被完全忽视的现象。并且现有技术中将导航卫星系统gnss所提供的监测站周围环境信息应用到滑坡灾害综合监测的研究还未涉及。


技术实现要素：

5.针对上述背景技术中所提出的问题，本发明提出了一种黄土滑坡灾害综合监测方法及系统，该方法为一种集“三维形变、大气水汽含量、土壤含水率”为一体的地基gnss遥感综合监测方法，验证了地基gnss在黄土浅层滑坡灾害综合监测方面具有出色的能力。
6.本发明提供了一种黄土滑坡灾害综合监测方法，包括以下步骤：
7.获取待监测区域的各项监测数据；
8.根据各项监测数据，利用gnss相对定位技术计算待监测区域的地表三维形变时间序列；利用gnss折射遥感技术计算待监测区域的大气水汽含量时间序列；利用gnss反射遥感技术计算待监测区域周围的地表土壤含水率时间序列；
9.获取的待监测区域的降雨量数据，分析其与大气水汽含量时间序列之间的响应关系，得到大气水汽含量与降雨量之间的关联性；
10.分析三维地表形变时间序列与地表土壤含水率时间序列之间的响应关系，得到待监测区域的滑坡变形与地表土壤含水率之间的关联性；
11.根据待监测区域的大气水汽含量与降雨量之间的关联性，以及解算的待监测区域大气水汽含量，对待监测区域进行降雨预警；根据待监测区域的滑坡变形与地表土壤含水率之间的关联性，以及解译的待监测区域的地表土壤含水率，对待监测区域进行滑坡综合监测。
12.进一步地，所述获取待监测区域的各项监测数据，包括：
13.在所述待监测区域设置gnss监测基准站、流动站，并获取所述待监测区域的gnss
测站观测数据；
14.获取所述待监测区域的广播星历数据；
15.获取所述待监测区域的精密星历数据；
16.获取所述待监测区域的气象数据。
17.进一步地，所述根据各项监测数据，利用gnss相对定位技术计算待监测区域的地表三维形变时间序列，包括以下步骤：
18.根据所述gnss测站观测数据，计算卫星的位置、速度、钟差；
19.根据所述gnss测站观测数据，利用伪距单点定位计算流动站、基准站的概略坐标；
20.根据所述gnss测站观测数据，计算所述基准站与流动站的非差残差项；
21.对所述gnss测站观测数据进行潮汐改正、对流层改正、天线相位中心改正；
22.对所述基准站与流动站的非差残差项进行卡尔曼滤波；
23.利用周跳探测中的m-w组合观测算法，并根据基准站与流动站的概略坐标，以及经过卡尔曼滤波处理后的非差残差项，确定整周模糊度；
24.提取流动站相对基准站的站心坐标neu；
25.根据流动站相对基准站的站心坐标neu，得到待监测区域的地表三维形变时间序列。
26.进一步地，所述根据各项监测数据，利用gnss折射遥感技术计算待监测区域的大气水汽含量时间序列，包括：
27.计算gnss测站观测数据中的天顶方向总延迟ztd；
28.从所述气象数据中获取待监测区域的气压、温度，并计算大气加权平均温度；
29.计算gnss测站观测数据中的天顶干延迟zhd；
30.计算湿延迟zwd，其计算公式为：
31.zwd＝ztd-zhd
ꢀꢀꢀꢀꢀ
(1)
32.根据湿延迟zwd，以及大气加权平均温度，计算不同时刻的大气水汽含量pwv，得到待监测区域的大气水汽含量时间序列。
33.进一步地，所述根据各项监测数据，利用gnss反射遥感技术计算待监测区域周围的地表土壤含水率时间序列，包括：
34.提取所述gnss测站观测数据中的信噪比snr数据；
35.选取观测卫星的高度角，并设置方位角，截取所述信噪比snr数据中发生高频振荡的信噪比snr数据；
36.利用低阶多项式将发生高频振荡的信噪比snr数据中的直、反射信号分离；
37.对信噪比snr数据中分离出的反射信号进行重采样；
38.对重采样后的所述反射信号进行ls频谱分析，得到gnss测站天线相位中心到地面的反射高h；
39.将所述反射高h、反射信号的波长λ、信噪比snr数据以及卫星高度角θ代入公式(2)，并利用非线性最小二乘法拟合得到相位参数ψ和幅值a；
[0040][0041]
根据所获取的相位参数ψ，建立线性反演模型获取待监测区域周围的地表土壤含水率时间序列。
[0042]
进一步地，所述分析大气水汽含量时间序列与降雨量数据之间的响应关系，得到大气水汽含量与降雨量之间的关联性，包括：
[0043]
绘制所述大气水汽含量时间序列和降雨量数据之间的双坐标轴关系图，从图中分析大气水汽含量时间序列和降雨量数据之间的关系，得到大气水汽含量与降雨量之间的关联性。
[0044]
进一步地，所述分析三维地表形变时间序列与地表土壤含水率时间序列之间的响应关系，得到待监测区域的滑坡变形与地表土壤含水率之间的关联性，包括：
[0045]
根据所述三维地表形变时间序列，得到滑坡累积位移；
[0046]
绘制所述滑坡累计位移与地表土壤含水率之间的双坐标轴关系图，得到所述滑坡累计位移与地表土壤含水率之间的关联性。
[0047]
进一步地，在所述绘制滑坡累计位移与地表土壤含水率之间的双坐标轴关系图时，通过时滞互相关分析方法定量分析滑坡位移速率对地表土壤含水率变化响应的滞后时间；
[0048]
依据所述滑坡位移速率对地表土壤含水率变化响应的滞后时间，绘制所述滑坡位移速率与地表土壤含水率之间的双坐标轴关系图，从图中分析所述滑坡位移速率与地表土壤含水率之间的相关性。
[0049]
本发明提供一种黄土滑坡灾害综合监测系统，包括：
[0050]
数据获取模块，用于获取待监测区域的各项监测数据；
[0051]
数据解算模块，用于根据各项监测数据，利用gnss相对定位技术计算待监测区域的地表三维形变时间序列；利用gnss折射遥感技术计算待监测区域的大气水汽含量时间序列；利用gnss反射遥感技术计算待监测区域周围的地表土壤含水率时间序列；
[0052]
第一响应关系获取模块，用于获取的待监测区域的降雨量数据，并分析其与大气水汽含量时间序列之间的响应关系，得到大气水汽含量与降雨量之间的关联性；
[0053]
第二响应关系获取模块，用于分析三维地表形变时间序列与地表土壤含水率时间序列之间的响应关系，得到待监测区域的滑坡变形与地表土壤含水率之间的关联性；
[0054]
综合监测模块，用于根据待监测区域的大气水汽含量与降雨量之间的关联性，以及解算的待监测区域的大气水汽含量，对待监测区域进行降雨预警；根据待监测区域的滑坡变形与地表土壤含水率之间的关联性，以及解译的待监测区域的地表土壤含水率，对待监测区域进行滑坡综合监测。
[0055]
与现有技术相比，本发明的有益效果：
[0056]
本发明提供了一种基于地基gnss遥感用于黄土滑坡灾害综合监测方法，该方法的优点在于基于一台gnss获取的数据集，分别采用gnss载波相位差分、gnss折射和gnss反射遥感技术，便可获取到滑坡测区的三维形变、大气水汽含量和土壤含水率信息，且gnss卫星所用的载波位于微波的l波段，受大雾和雨雪等恶劣天气的影响较小，且能够穿透云层，这有利于实现对滑坡形变及周围环境的长期稳定、全天候的观测。伴随着gnss技术的不断改进、全球cors站的不断完善和未来全球卫星导航系统卫星总数的不断增加，这将会为研究gnss解译滑坡地区环境与监测地表形变提供丰富的数据，进一步推进该技术在滑坡体区域环境监测方面的研究与应用，能够使gnss在滑坡灾害监测和预警等方面更好地发挥作用，此工作对于滑坡的综合监测具有重要意义。
附图说明
[0057]
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
[0058]
图1为本发明提供的一种黄土滑坡灾害综合监测方法的实施例中的综合监测体系图；
[0059]
图2为本发明实施例提供的待监测区域的三维形变时间序列。
[0060]
图3为本发明实施例提供的gnss解算后的大气水汽含量pwv和待监测区域的流动站相对于基准站的e方向形变；
[0061]
图4为本发明实施例提供的待监测区域的流动站相对于基准站的180
°‑
360
°
方位角内部分卫星菲涅尔反射区；
[0062]
图5为本发明实施例提供的各卫星解译土壤湿度与实测土壤湿度对比结果及相关性分析图；
[0063]
图6为本发明实施例提供的融合多卫星解译的土壤湿度结果与实测土壤湿度对比及相关性分析图；
[0064]
图7为本发明实施例提供的土壤湿度与位移关系图；
[0065]
图8为本发明实施例提供的位移速率与土壤湿度关系图；
[0066]
图9为本发明实施例提供的位移速率-土壤湿度时滞互相关序列；
[0067]
图10为本发明实施例提供的位移速率与土壤湿度变化速率关系图。
具体实施方式
[0068]
下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，但应当理解本发明的保护范围并不受具体实施方式的限制。
[0069]
如图1所示，本发明提供一种黄土滑坡灾害综合监测方法，包括以下步骤：
[0070]
步骤1：获取待监测区域的各项监测数据。
[0071]
步骤2：根据各项监测数据，利用gnss相对定位技术计算待监测区域的地表三维形变时间序列；利用gnss折射遥感技术计算待监测区域的大气水汽含量时间序列；利用gnss反射遥感技术计算待监测区域周围的地表土壤含水率时间序列；
[0072]
步骤3：获取的待监测区域的降雨量数据，分析其与大气水汽含量时间序列之间的响应关系，得到大气水汽含量与降雨量之间的关联性；；
[0073]
步骤4：分析三维地表形变时间序列与地表土壤含水率时间序列之间的响应关系，得到待监测区域的滑坡变形与地表土壤含水率之间的关联性；
[0074]
步骤5：根据待监测区域的大气水汽含量与降雨量之间的关联性，以及解算的待监测区域大气水汽含量，对待监测区域进行降雨预警；根据待监测区域的滑坡变形与地表土壤含水率之间的关联性，以及解译的待监测区域的地表土壤含水率，对待监测区域进行滑坡综合监测。
[0075]
在步骤1中，获取待监测区域的各项监测数据，具体包括：
[0076]
在待监测区域设置gnss监测基准站、流动站，并获取待监测区域的gnss测站观测数据；
[0077]
获取待监测区域的广播星历数据；
[0078]
获取待监测区域的精密星历数据；
[0079]
获取待监测区域的气象数据。
[0080]
使用teqc、rtklib等软件分析gnss测站观测数据中的利用率、多路径效应、周跳、信噪比、电离层延迟率和精度因子等数据，从而实现评估实时获取到的gnss测站观测数据质量的优劣。
[0081]
在步骤2中具体包括如下步骤：
[0082]
步骤2.1：根据各项监测数据，利用gnss相对定位技术计算待监测区域的地表三维形变时间序列，包括以下步骤：
[0083]
根据gnss测站观测数据，计算卫星的位置、速度、钟差；
[0084]
根据gnss测站观测数据，利用伪距单点定位计算流动站、基准站的概略坐标；
[0085]
根据gnss测站观测数据，计算基准站与流动站的非差残差项；
[0086]
对gnss测站观测数据进行潮汐改正、对流层改正、天线相位中心改正；
[0087]
对所述基准站与流动站的非差残差项进行卡尔曼滤波；
[0088]
利用周跳探测中的m-w组合观测算法，并根据基准站与流动站的概略坐标，以及经过卡尔曼滤波处理后的非差残差项，确定整周模糊度；；
[0089]
提取流动站相对基准站的站心坐标neu；
[0090]
根据流动站相对基准站的站心坐标neu，得到待监测区域的地表三维形变时间序列。
[0091]
步骤2.2：根据各项监测数据，利用gnss折射遥感技术计算待监测区域的大气水汽含量时间序列，包括：
[0092]
计算gnss测站观测数据中的天顶方向总延迟ztd；
[0093]
具体利用gamit/globk数据处理软件来gnss测站观测数据中的天顶方向总延迟ztd，在进行计算时选择全局投影函数模型以及gamita默认10
°
高度截止角。
[0094]
从气象数据中获取待监测区域的气压、温度，并计算大气加权平均温度；
[0095]
计算gnss测站观测数据中的天顶干延迟zhd，具体利用使用最常用的萨斯塔莫宁(saastamoinen)模型计算天顶干延迟zhd。
[0096]
计算湿延迟zwd，其计算公式为：
[0097]
zwd＝ztd-zhd
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0098]
根据湿延迟zwd，以及大气加权平均温度，计算不同时刻的大气水汽含量pwv，得到待监测区域的大气水汽含量时间序列。
[0099]
步骤2.3：根据各项监测数据，利用gnss反射遥感技术计算待监测区域周围的地表土壤含水率时间序列，包括：
[0100]
提取gnss测站观测数据中的信噪比snr数据；
[0101]
选取观测卫星的低高度角，并设置方位角，截取信噪比snr数据中发生高频振荡的信噪比snr数据；
[0102]
利用低阶多项式将发生高频振荡的信噪比snr数据中的直、反射信号分离；
[0103]
对信噪比snr数据中分离出的反射信号进行重采样；
[0104]
对重采样后的反射信号进行ls频谱分析，得到gnss测站天线相位中心到地面的距离，也称为反射高h；
[0105]
将反射高h、信号波长λ、信噪比snr以及卫星高度角θ代入公式(2)，利用非线性最小二乘法拟合得到相位参数ψ和幅值a；
[0106][0107]
根据所获取的相位参数ψ与土壤湿度具有显著的线性关系，建立线性反演模型获取待监测区域周围的地表土壤含水率时间序列。
[0108]
由于反射高的相位与土壤湿度具有显著的线性关系，故根据反射高的相位参数ψ，建立线性反演模型获取待监测区域周围的地表土壤含水率时间序列。
[0109]
在步骤3中，分析大气水汽含量时间序列与降雨量数据之间的响应关系，得到大气水汽含量与降雨量之间的关联性，包括：
[0110]
绘制大气水汽含量时间序列和降雨量数据之间的双坐标轴关系图，从图中分析大气水汽含量时间序列和降雨量数据之间的关系。
[0111]
在步骤4中，分析三维地表形变时间序列与地表土壤含水率时间序列之间的响应关系，得到待监测区域的滑坡变形与地表土壤含水率之间的关联性，包括：
[0112]
根据三维地表形变时间序列，得到滑坡累积位移；
[0113]
绘制滑坡累计位移与地表土壤含水率之间的双坐标轴关系图，得到滑坡累计位移与地表土壤含水率之间的关联性。
[0114]
步骤4：还包括在绘制滑坡累计位移与地表土壤含水率之间的双坐标轴关系图时，通过时滞互相关分析方法定量分析滑坡位移速率对地表土壤含水率变化响应的滞后时间；
[0115]
依据滑坡位移速率对地表土壤含水率变化响应的滞后时间，绘制滑坡位移速率与地表土壤含水率之间的双坐标轴关系图，从图中分析滑坡位移速率与地表土壤含水率之间的相关性。
[0116]
本发明提供一种黄土滑坡灾害综合监测系统，包括：
[0117]
数据获取模块，用于获取待监测区域的各项监测数据；
[0118]
数据解算模块，用于根据各项监测数据，利用gnss相对定位技术计算待监测区域的地表三维形变时间序列；利用gnss折射遥感技术计算待监测区域的大气水汽含量时间序列；利用gnss反射遥感技术计算待监测区域周围的地表土壤含水率时间序列；
[0119]
第一响应关系获取模块，用于获取的待监测区域的降雨量数据，并分析其与大气水汽含量时间序列之间的响应关系，得到大气水汽含量与降雨量之间的关联性；
[0120]
第二响应关系获取模块，用于分析三维地表形变时间序列与地表土壤含水率时间序列之间的响应关系，得到待监测区域的滑坡变形与地表土壤含水率之间的关联性；
[0121]
综合监测模块，用于根据待监测区域的大气水汽含量与降雨量之间的关联性，以及解算的待监测区域的大气水汽含量，对待监测区域进行降雨预警；根据待监测区域的滑坡变形与地表土壤含水率之间的关联性，以及解译的待监测区域的地表土壤含水率，对待监测区域进行滑坡综合监测。
[0122]
下面结合具体的实施例对本发明中的技术方案作进一步详细说明。
[0123]
1、获取待监测区域的地表三维形变时间序列。
[0124]
如图2为本发明的技术方案解算的实例区域的典型特征点上测站三个方向的时间序列形变图。
[0125]
图2(a)是gamit解算的2020年测站全年的形变时间序列，从中可以看出，2020年第
1天至245天时间段内，该测站在北n、东e、高u三个方向均无明显形变，形变主要发生在第245天至第313天，且东方向形变最为明显，此后该测站又无明显形变，三个方向均趋于稳定状态。
[0126]
为更加细致分析滑坡形变，本技术用track单历元解算了第245天至第335天的数据，结果如图2(b)所示，分析可知该时间段内测站点n方向仅变化了2cm左右，e方向变化了78cm左右，u方向变化了12cm左右。则可知滑坡主要滑动方向为东方向。并且，在形变时间序列的开始时段，测站的形变速率较慢且较为平稳，但是在第273天开始，测站形变迅速增大，在第287天至302天滑坡形变趋于稳定，第302天至第313天又出现一次较小的形变，到最后时段，测站的形变又变得平稳。
[0127]
2、获取待监测区域的大气水汽含量时间序列。
[0128]
图3为本发明的技术方案解算的实例区域在2020年第245天至第335天发生明显三次形变时的e方向形变(滑动的主要方向)、大气水汽含量pwv以及由传感器采集的降雨量。线型最粗的柱状表示降雨量，线型最细的线表示形变序列，线型次细的线表示大气水汽含量pwv。从图中可以看出，在第275天前后，存在较多的降雨量事件且降雨量较大，此时的形变速率达到最大，说明降雨量是诱导此次滑坡发生的一个重要因素。并且可以发现，在降雨量事件发生前大气水汽含量pwv总是存在累积的过程，而且降雨量后大气水汽含量pwv会迅速下降。由降雨量事件引起的滑坡通常具有一定的滞后性，则降雨量事件的监测对于滑坡灾害预警有着重要的意义，而大气水汽含量pwv和降雨量事件有较强的相关性。
[0129]
3、获取待监测区域周围的地表土壤含水率时间序列，并分析大气水汽含量时间序列与降雨量数据之间的响应关系。
[0130]
图4为本技术方案解算的实例区域测站在180
°
至360
°
方位角范围内不同低高度角下各卫星的第一菲涅尔反射区，即选择的反演土壤含水率的区域。不同的颜色代表了不同的卫星，该图结合google earth绘制而成。
[0131]
图5为本技术方案解算的实例区域四颗卫星解译的土壤湿度与实测土壤湿度的对比结果及相关性分析，由图5可以发现，该技术解算的各颗卫星结果曲线与土壤湿度曲线有较好的一致性，都随实测土壤湿度呈现增减的趋势，符合土壤湿度变化的情况，相关系数均是显著的，分别达到0.78、0.85、0.77和0.88。验证了gnss-ir解译结果的可靠性。
[0132]
图6为本技术方案解算的实例区域融合多卫星估计的测站的土壤湿度、传感器采集的降雨量和土壤湿度。线型最粗的柱状表示降雨量，线型最细的线表示gnss估计的土壤湿度，线型次细的线表示原位土壤湿度。从图中可以看出，伴随着第274天和第302天这两次降水事件，土壤湿度曲线明显上升，因后期依旧出现降雨量天气，土壤重量含水率在随后的第276和303天达到峰值，土壤重量含水率分别为33.34％和32.85％。在第276-301天与第303-328天时间段内，降雨量出现明显的递减趋势，土壤湿度发生缓慢回落，直到新的降水事件的出现，则再一次出现上浮波动。多颗卫星融合结果的吻合度比单颗卫星表现出了更大的一致性，对融合后结果与实测土壤重量含水率进行相关性分析(图6(右))，相关系数分别达到0.89，相比单颗卫星的相关性有了一定的提升。
[0133]
4、分析滑坡累计位移与地表土壤含水率之间的响应关系。
[0134]
图7为本技术方案解算的实例区域的监测点相应时段的累计位移与土壤含水率关系图，结果见图7所示，无论是土壤湿度的监测结果，还是gnss解译的土壤湿度值，在三次形
变发生时土壤湿度值都出现了明显的提升，实测土壤重量含水率在第257天、274天和302天分别从30.48％、30％和26.87％提升至31.92％、31.12％和28.43％，增幅依次为4％、4％和6％。解译的结果虽然由于反演误差的原因在第274天反而出现下降，但在随后的275天土壤重量含水率得到提升，在第257天、275天和302天的增幅依次为4％、5％和6％。结果表明，土壤中水分的变化影响了土体的抗剪强度，是黄土滑坡形变的最主要诱发因素之一。
[0135]
5、分析滑坡位移速率与地表土壤含水率之间的相关性
[0136]
图8为本技术方案解算的实例区域的监测站点的位移速率与土壤湿度关系图。由图8通过目视解译可以看出，在监测时间内，位移速率与土壤湿度有着明显的相同变化趋势，可以判断两者存在一定相关性。
[0137]
图9为本技术方案解算的实例区域的监测站点的位移速率与土壤湿度时滞互相关序列图。由图9可看出：土壤湿度的监测结果与gnss解译的土壤湿度值在滞后周期为0d时，它们对位移变化的影响达到峰值，实测土壤湿度的相关性极值为0.7，gnss解译的土壤湿度的相关性极值为0.65，均为高度相关。
[0138]
图10为本技术方案解算的实例区域的监测站点的位移速率与土壤含水率变化速率关系图。由图10可以看出土壤湿度变化速率达到峰值之后的第258、276与304天的形变速率也在三个时间段内达到各自的峰值，且形变速率随着实测及解译土壤含水率变化速率而变化，二者在整体的变换趋势上是一致的，说明土壤含水率变化速率的骤增同样是影响黄土滑坡形变的因素之一，二者具有良好的响应关系。
[0139]
最后说明的是：以上公开的仅为本发明的一个具体实施例，但是，本发明实施例并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。