CN115638767B - 地面沉降监测方法 - Google Patents
地面沉降监测方法 Download PDFInfo
- Publication number
- CN115638767B CN115638767B CN202211390802.9A CN202211390802A CN115638767B CN 115638767 B CN115638767 B CN 115638767B CN 202211390802 A CN202211390802 A CN 202211390802A CN 115638767 B CN115638767 B CN 115638767B
- Authority
- CN
- China
- Prior art keywords
- data
- monitoring
- landslide
- deformation
- sequence
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000012544 monitoring process Methods 0.000 title claims abstract description 168
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000001914 filtration Methods 0.000 claims abstract description 41
- 238000005516 engineering process Methods 0.000 claims abstract description 26
- 230000004927 fusion Effects 0.000 claims abstract description 26
- 238000007637 random forest analysis Methods 0.000 claims abstract description 20
- 230000003044 adaptive effect Effects 0.000 claims abstract description 17
- 238000007781 pre-processing Methods 0.000 claims abstract description 6
- 238000006073 displacement reaction Methods 0.000 claims description 36
- 238000005259 measurement Methods 0.000 claims description 10
- 238000012806 monitoring device Methods 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 238000005316 response function Methods 0.000 claims description 8
- 239000002689 soil Substances 0.000 claims description 7
- 238000012935 Averaging Methods 0.000 claims description 6
- 230000002159 abnormal effect Effects 0.000 claims description 6
- 230000001186 cumulative effect Effects 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 5
- 238000009499 grossing Methods 0.000 claims description 5
- 230000007774 longterm Effects 0.000 claims description 5
- 238000013178 mathematical model Methods 0.000 claims description 5
- 230000006978 adaptation Effects 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 4
- 238000011160 research Methods 0.000 claims description 4
- 238000012360 testing method Methods 0.000 claims description 4
- 238000012876 topography Methods 0.000 claims description 4
- 238000009825 accumulation Methods 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims description 2
- 239000013598 vector Substances 0.000 claims description 2
- 238000004458 analytical method Methods 0.000 abstract description 4
- 238000004088 simulation Methods 0.000 description 11
- 230000001174 ascending effect Effects 0.000 description 9
- 238000010586 diagram Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 238000002360 preparation method Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000010219 correlation analysis Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000011439 discrete element method Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000011010 flushing procedure Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000008595 infiltration Effects 0.000 description 1
- 238000001764 infiltration Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000006641 stabilisation Effects 0.000 description 1
- 238000011105 stabilization Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Pit Excavations, Shoring, Fill Or Stabilisation Of Slopes (AREA)
- Testing Or Calibration Of Command Recording Devices (AREA)
Abstract
本发明的一种地面沉降监测方法，包括建立随机森林模型，获取数据并对获取到的数据信息进行预处理；对得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合；对得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测；对得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分，完成地面沉降监测。使用自适应Kalman滤波技术对包括RTK定位数据、无人机摄影测量数据、传感器数据在内的多源数据进行融合分析，将滑坡监测精度提高到了毫米级；使用RTK技术和土工带传感器，克服了天气、植被覆盖对监测的影响；使用灰色预测理论对山体滑坡监测点进行形变预测，结合蠕变切线角判据实现对山体滑坡预警等级的划分。
Description
技术领域
本发明涉及地质灾害防治技术领域，具体涉及一种地面沉降监测方法。
背景技术
山体滑坡是指山体上的土体或者岩体在河流冲刷、雨水浸润、重力、人类活动等因素的影响下，沿着某一滑面向下滑动的自然灾害现象，具有突发性和隐蔽性强、破坏力大的特点。传统的山体滑坡监测方法主要有宏观地质观测法、简易观测法、设站观测法、仪表观测法等，它们普遍遍存在的问题是数据的采集需要人工定期到现场进行，缺少实时性，并且监测精度较低。目前主流的山体滑坡监测方法是通过合成孔径雷达技术(SAR)进行监测，但是该方法受天气、植被覆盖的影响大，实时性差，难以实现长期监测。
全球卫星导航系统(GNSS)的高精度三维定位已经成为了地质灾害监测中的一种常用的技术手段。GNSS定位所需时间短、全天候使用、不受测站间障碍物影响、可提供三维定位信息、可实现监测自动化。GNSS实时动态差分(RTK)技术、精密单点定位技术(PPP)等新型的GNSS高精度定位技术在地质形变监测中的作用日益的凸显出来。RTK技术通过位置已知的固定的基准站和移动的形变监测站进行差分，可以基本实现实时定位。并且在静态情况下的定位精度可以保持在毫米级。
Kalman滤波的应用极其广泛，涵盖了飞行器的导引、导航及控制，机器人运动规划及控制，位移监测等诸多领域。作为位移监测中的重要数据处理算法，kalman滤波广泛地应用于建筑、桥梁、地质灾害的形变监测。
影响山体滑坡的因素复杂多变，并且具有很强的随机性。同时在诸多影响因素中有些是已知的，还有一些是暂时未被发现的，也就是这些复杂的影响因素对山体滑坡的作用具有典型的“灰色”特征。所谓“灰色”特征是指目标系统介于已知和未知之间。对于具有这样特征的系统，邓聚龙教授首创了一种灰色系统理论对其进行了深入的分析。灰色系统理论的内容包括灰色朦胧集、灰色关联分析、灰色序列生成、灰色模型，实现对目标系统的建模、评估、预测、优化等功能。灰色系统理论在解决信息不完全的系统问题时具有十分显著的优势。
发明内容
本发明提出的一种地面沉降监测方法，具体涉及到一种基于自适应kalman滤波和灰色预测理论的山体滑坡监测预警装置；可至少解决上述技术问题之一。
为实现上述目的，本发明采用了以下技术方案：
一种地面沉降监测方法，包括：
步骤一：建立随机森林模型，获取并分析气象水文、地质构造以及植被覆盖数据，对目标区域是否会发生山体滑坡进行预测，选出会发生山体滑坡的监测点；对判断得到的会发生滑坡的监测点进行重点监测，实时定位，采集地面沉降变化量，获取目标区域的数字高程数据及高精度水平数字信息，并对获取到的数据信息进行预处理；
步骤二：对步骤一中得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合；
步骤三：对步骤二中得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测；
步骤四：对步骤三中得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分，完成地面沉降监测。
进一步的，在步骤一中所述的建立随机森林模型选取滑坡重点监测点，采集这些监测点处的数据信息并进行预处理，其具体做法如下：
S11、对无人机摄影测量得到的DEM数据进行分析，并结合当地历史水文气象数据，得到目标监测区域的高程、坡度、坡向、粗糙度、地形起伏度、岩土体类型、断层缓冲区、降雨量、河流缓冲区、道路缓冲区、归一化植被等数据信息；
S12、使用S11得到的数据集建立随机森林模型，并测试建立的随机森林模型的预测准确率；
S13、将S11中得到的数据信息输入到S12建立起来的随机森林模型中，得到山体滑坡发生与否的初步预测结果；
S14、实时动态载波相位差分技术RTK进行高精度定位，精准获取监测装置的实时位置信息；
S15、使用土工带传感器采集地面沉降变化量；
S16、使用无人机拍摄获取数字高程数据、高精度水平数字信息及山体地质数据；
S17、对RTK定位接收到的经度、维度、高度数据进行中值滤波，剔除异常值；然后设置参考时刻，通过实时获得的经纬高信息计算任意时刻RTK监测装置与参考时刻相比移动的距离；
S18、土工带传感器按照网格状排列，对土工带传感器返回的数据，首先进行中值滤波来剔除异常值；然后再几何层面进行平滑，将所有行的土工带传感器取平均，将所有列的土工带传感器取平均，使用勾股定理得到针对该监测点的形变量；
S19、进行无人机飞行测量之后得到目标区域的数字高程地图DEM及高精度水平数字地图，对于不同期的摄影测量结果，通过对高精度水平数字地图作差，计算得到监测点处的水平位移量及位移方向；通过数字高程地图作差，计算得到监测点处的垂直方向的位移量，最终测量得到山体滑坡监测点的位移量。
进一步的，在步骤二中所述的对步骤一中得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合，其具体做法如下：
S21、对山体滑坡运动进行分析，建立监测点形变的运动数学模型；
S22、基于S21中建立的运动数学模型，使用自适应kalman滤波对RTK定位数据、土工带传感器数据、无人机摄影测量数据进行数据融合。
进一步的，所述在步骤三中所述的对步骤二中得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测，其具体做法如下：
S31、选取若干期的形变监测值组成长度为n的原始观测序列
S32、生成一次累加和序列按照灰色预测模型参数计算公式的需要，选后n-2、n-1期的数据分别生成一次累加和序列/>和/>
S33、依次生成和/>的一次累加和序列的紧邻均值序列/>和/>
S34、使用上述准备好的数据序列计算灰色方程x(0)(k)+ax(1)(k)＝b的系数a和b，分别求得序列和/>的灰色方程系数[a1,b1]、[a2,b2]、[a3,b3]；
S35、在上述步骤中求出了序列的灰色方程系数[a1,b1]，即已经得到了该序列的一次累加和序列即位置信息序列的灰色预测时间响应函数：
使用该预测时间响应函数对未来时刻的一次累加和序列即位置信息序列进行预测；对位置信息序列进行累减生成得到位移序列；
S36、分别计算序列和/>的均值和方差；
S37、通过位移预测数据前8期的数据序列与数据序列/>作差求出残差序列；然后计算该残差序列的均值和方差；
S38根据上述准备数据计算该灰色预测模型的均方差比值C和小概率误差P，根据提前设定的模型精度等级参考对该灰色理论模型的精度进行评估分析，进而判断山体滑坡灾害预测预报结果的可信度。
进一步的，在步骤四中所述的对步骤三中得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分；具体做法如下：
S41、建立蠕变切线角预警判据；
S42根据行业内的长期观测研究得到山体滑坡过程的预警模型，实现了对山体滑坡的合理等级预警。
由上述技术方案可知，本发明的地面沉降监测方法针对当前的山体滑坡监测技术精度低、实时性差的问题，提出了一种基于灰色理论的自适应Kalman滤波滑坡监测技术。该技术通过自适应Kalman滤波融合包括RTK定位数据、无人机摄影测量数据、土工带传感器数据在内的多源数据，极大程度上提高了山体滑坡监测精度。经仿真实验分析形变监测精度能够达到毫米级，误差在3％以内。使用灰色预测理论模型之后，该技术能够准确地预测未来一段时间内监测点的形变量，最后通过蠕变切线角判据可以准确地给出在未来一周内山体滑坡危险等级，从而做到对山体滑坡灾害的提前预防。RTK定位技术、土工带传感器的使用，使得山体滑坡监测功能受天气状况和植被覆盖等因素的影响大幅减弱，同时保证了监测系统的自动化，免去了人工定期到现场测算的麻烦。与已有的一些滑坡监测方法相比，该技术在提高滑坡监测精度、增强监测实时性、灾害及时预测等方面具有显著的优势，有很高的实用价值。
附图说明
图1是本发明实施例的方法流程图；
图2是土工带传感器的排布方式；
图3是蠕变切线角判据下的T-t曲线；
图4是无人机摄影测量得到的目标区域的数字高程地图；
图5是无人机摄影测量得到的目标区域的数字水平地图；
图6是本发明实施例的多源数据融合效果展示图；
图7是本发明实施例的多源数据融合曲线标注图；
图8是本发明实施例的警示级别数据仿真结果图；
图9是本发明实施例的警戒级别数据仿真结果图；
图10是本发明实施例的警报级别数据仿真结果图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。
如图1所示，本实施例所述的地面沉降监测方法，包括以下步骤，
步骤一：建立随机森林模型，获取并分析气象水文、地质构造以及植被覆盖等数据，对目标区域是否会发生山体滑坡进行预测，选出会发生山体滑坡的监测点；对判断得到的会发生滑坡的监测点进行重点监测，高精度实时定位，采集地面沉降变化量，获取目标区域的数字高程数据及高精度水平数字信息，并对获取到的数据信息进行预处理；
步骤二：对步骤一中得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合；
步骤三：对步骤二中得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测；
步骤四：对步骤三中得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分。最终实现地面沉降监测装置的功能。
其中，在步骤一中所述的建立随机森林模型选取滑坡重点监测点，采集这些监测点处的数据信息并进行预处理，其具体做法如下：
S11、对无人机摄影测量得到的DEM数据进行分析，并结合当地历史水文气象数据，得到目标监测区域的高程、坡度、坡向、粗糙度、地形起伏度、岩土体类型、断层缓冲区、降雨量、河流缓冲区、道路缓冲区、归一化植被等数据信息；
S12、使用S11得到的数据集建立随机森林模型，并测试建立的随机森林模型的预测准确率；
S13、将S11中得到的数据信息输入到S12建立起来的随机森林模型中，得到山体滑坡发生与否的初步预测结果；
S14、实时动态载波相位差分技术(RTK)进行高精度定位，精准获取监测装置的实时位置信息；
S15、使用土工带传感器采集地面沉降变化量；
S16、使用无人机拍摄获取数字高程数据、高精度水平数字信息及山体地质数据；
S17、对RTK定位接收到的经度、维度、高度数据进行中值滤波，剔除异常值。然后设置参考时刻，通过实时获得的经纬高信息计算任意时刻RTK监测装置与参考时刻相比移动的距离；
S18、土工带传感器按照网格状排列(5行5列)，对土工带传感器返回的数据，首先进行中值滤波来剔除异常值。然后再几何层面进行平滑，将所有行的土工带传感器取平均，将所有列的土工带传感器取平均，使用勾股定理得到针对该监测点的形变量；
S19、进行无人机飞行测量之后可以得到目标区域的数字高程地图(DEM)及高精度水平数字地图，对于不同期的摄影测量结果，通过对高精度水平数字地图作差，可以计算得到监测点处的水平位移量及位移方向；通过数字高程地图作差，可以计算得到监测点处的垂直方向的位移量，最终测量得到山体滑坡监测点的位移量。
其中，在步骤二中所述的对步骤一中得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合，其具体做法如下：
S21、对山体滑坡运动进行分析，建立监测点形变的运动模型；
S22、基于S21中建立的运动数学模型，使用自适应kalman滤波对RTK定位数据、土工带传感器数据、无人机摄影测量数据进行数据融合。
其中，在步骤三中所述的对步骤二中得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测，其具体做法如下：
S31、选取若干期的形变监测值组成长度为n的原始观测序列
S32、生成一次累加和序列按照灰色预测模型参数计算公式的需要，选后n-2、n-1期的数据分别生成一次累加和序列/>和/>
S33、依次生成和/>的一次累加和序列的紧邻均值序列/>和/>
S34、使用上述准备好的数据序列通过下列4式：
Y＝[x(0)(2),x(0)(3),…,x(0)(n)]T (17)
计算灰色方程x(0)(k)+ax(1)(k)＝b的系数a和b。分别求得序列和/>的灰色方程系数[a1,b1]、[a2,b2]、[a3,b3]。
S35、在上述步骤中求出了序列的灰色方程系数[a1,b1]，即已经得到了该序列的一次累加和序列即位置信息序列的灰色预测时间响应函数：
使用该预测时间响应函数可以对未来时刻的一次累加和序列即位置信息序列进行预测。对位置信息序列进行累减生成得到位移序列。
S36、分别计算序列和/>的均值和方差。
S37、通过位移预测数据前8期的数据序列与数据序列/>作差求出残差序列。然后计算该残差序列的均值和方差。
S38根据上述准备数据利用式
计算该灰色预测模型的均方差比值C。其中X为原始目标序列，为模型预测序列。为原始目标序列X的方差，/>为残差序列ε的方差。
根据上述准备数据利用式
计算该灰色预测模型的小概率误差P。式中为残差序列ε的均值。根据表1所示的模型精度等级参考对该灰色理论模型的精度进行评估分析，进而判断山体滑坡灾害预测预报结果的可信度。
| 模型精度等级 | 优 | 良 | 中 | 差 |
| 均值方差比A | ≤0.35 | 0.35～0.5 | 0.5～0.65 | >0.65 |
| 小概率误差P | ≥0.95 | 0.8～0.95 | 0.7～0.8 | <0.7 |
表1
其中，在步骤四中所述的对步骤三中得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分。最终实现地面沉降监测装置的功能，具体做法如下：
S41、建立蠕变切线角预警判据；
S42根据行业内的长期观测研究得到山体滑坡过程的预警模型，实现了对山体滑坡的合理等级预警。
以下举例说明：
第一步：使用已有山体滑坡各影响因子数据集建立随机森林模型，并测试建立的随机森林模型的预测准确率。对使用无人机摄影测量得到DEM数据进行分析，并结合当地历史水文气象数据，得到目标监测区域的高程、坡度、坡向、粗糙度、地形起伏度、岩土体类型、断层缓冲区、降雨量、河流缓冲区、道路缓冲区、归一化植被等数据信息，建立目标区域的山体滑坡各影响因子数据集。将目标区域的山体滑坡各影响因子数据集输入到已经建立的随机森林模型中，得到山体滑坡发生与否的初步预测结果。对判断得到的会发生滑坡的监测点进行重点监测，高精度实时定位，采集地面沉降变化量，获取目标区域的数字高程数据及高精度水平数字信息，并对获取到的数据信息进行预处理。
通过RTK定位得到监测点的经度、维度和天线高度，则监测点位移量dL与竖直位移量dL1和水平位移量dL2满足：
dL1＝H1-H2 (3)
其中θ1、φ1分别为在起始时刻监测点的经纬度；θ2、φ2为在观测时刻监测点的经纬度；R为地球半径，取平均值6371km。
基于土工带电阻值计算其形变量dL的计算公式如下：
式中，R(mΩ)为土工带在观测时刻的电阻测量值；R0(mΩ)是土工带在完成布置时的初始电阻值；L0(mm)是土工带在完成布置时的初始长度；k是比例系数；T(℃)为土壤温度。
因为土工带排布网络相对尺寸较小，所以在对一套监测设备进行分析时，忽略山体凹凸性，将土工带排布网络监测点处的切平面。土工带传感器排布方式如图2所示。10条土工带分成纵横两组，各有5条土工带。每一条土工带上有10个土工带传感器。10条土工带按照网格状排列，空心点是土工带传感器，实心点是两个土工带传感器重合点。RTK接收机安装在网络的中心位置。
首先对土工带数据进行中值滤波和平滑处理，得到监测点的纵向为位移信息和横向为位移信息y，则监测点的形变量:
通过无人机近景摄影测量技术获得监测区域的三维点云数据、倾斜摄影处理结果。使用ArcMap对多期测量数据中的影像控制点进行匹配之后，在某一期的三维模型中加入另一起的相对三维模型，计算出前后两期摄影测量中滑坡监测点处的绝对位移量，并获得形变区域的位置和大小。目标区域的数字高程地图(DEM，Digital Elevation Map)及高精度水平数字地图分别如图4和图5所示。
第二步：对重点监测点处的多源形变监测数据使用自适应kalman滤波技术进行数据融合。
对山体滑坡运动进行分析，建立监测点形变的运动模型，其描述方程如下：
其中，x(i)和分别表示形变监测点在时刻i的位置和速度；a(i)表示形变监测点从i时刻到i+1时刻形变的加速度。通过上述滑坡监测运动模型建立状态方程和观测方程如下：
其中， V(i)是均值为0，方差为/>的高斯白噪声序列，W(i)是均值为0，方差为/>的高斯白噪声序列。Z(i)为滑坡形变监测点位移、速度、加速度的观测向量。其中T为监测数据的采集周期，本文中T为2天。
基于上述建立的山体滑坡监测点运动模型，使用自适应kalman滤波对RTK定位数据、土工带传感器数据、无人机摄影测量数据进行数据融合。自适应Kalman滤波的原理可以通过如下方程加以描述：
时间更新方程：
Pi/i-1＝APi-1/i-1AT+Q (11)
状态更新方程：
其中，D表示观测噪声的等价协方差矩阵，观测噪声的等价权矩阵为D-1；ai表示自适应因子(0<ai≤1)。自适应kalman滤波的关键在于自适应因子ai的构造。
自适应因子的构造采用两段函数法，两段函数模型为：
式中，c为常量，一般取c＝1～2.5，0<ak≤1。
通过对山体滑坡进行运动建模，使用自适应kalman滤波可以实现对多源滑坡形变监测数据的融合功能，从而得到对于山体滑坡形变监测点的极为精确地形变量，在本装置中监测精度可以达到毫米级。本次步骤是山体滑坡监测的核心步骤，为之后的灰色预测理论模型预测做数据准备。
在本装置中，多源数据融合效果如图6所示。图中拟合实线为经过自适应kalman滤波之后的滑坡监测点的形变量。需要说明的是，图中的十字点既有土工带传感器的形变监测数据也有RTK定位监测的形变数据。此外，由于自适应kalman滤波结果同时还收到了无人机近景摄影测量计算得到的滑坡监测点的形变加速度影响，拟合实线并不处于土工带传感器监测数据和RTK定位监测形变数据的中间位置，这个效果图是很合理的。为了进一步展示自适应kalman滤波算法进行数据处理的效果，将滤波后的形变量数据单独显示，找到滤波收敛之后的较为稳定的区域，标记出数据波动的最大最小值，如图7所示。从图中可以发现，滤波收敛并稳定之后数据的最大值为3.864mm，最小值为3.686mm，二者相差0.178mm，满足监测精度要求。
第三步：通过上述多源监测数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测。
选取若干期的形变监测值组成长度为n的原始观测序列生成一次累加和序列/>按照灰色预测模型参数计算公式的需要，选后n-2、n-1期的数据分别生成一次累加和序列/>和/>依次生成/>和/>的一次累加和序列的紧邻均值序列和/>使用上述准备好的数据序列通过式(19)计算灰色方程x(0)(k)+ax(1)(k)＝b的系数a和b，可分别求得序列/>和/>的灰色方程系数[a1,b1]、[a2,b2]、[a3,b3]。
在上述步骤中求出了序列的灰色方程系数[a1,b1]，即已经得到了该序列的一次累加和序列即位置信息序列的灰色预测时间响应函数：
使用该预测时间响应函数可以对未来时刻的一次累加和序列即位置信息序列进行预测。对位置信息序列进行累减生成得到位移序列。分别计算序列和/>的均值和方差。通过位移预测数据/>前8期的数据序列与数据序列/>作差求出残差序列。然后计算该残差序列的均值和方差。通过位移预测数据/>前8期的数据序列与数据序列作差求出残差序列。然后计算该残差序列的均值和方差。根据上述准备数据利用式(20)、(21)和(22)计算该灰色预测模型的均方差比值C和小概率误差P，根据表1所示的模型精度等级参考对该灰色理论模型的精度进行评估分析，进而判断山体滑坡灾害预测预报结果的可信度。
在山体滑坡预测的问题中，影响系统预测结果的因素纷繁复杂，很难考虑到影响系统预测结果的所有因素。在此情形之下，灰色关联分析的优势便能很好的体现出来，它可以在已知信息不完全的情况下对系统各因素进行分析研究，从而得到较为精确山体滑坡预测结果。
第四步：对得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分。最终实现地面沉降监测装置的功能。
建立蠕变切线角预警判据，具体内容如下。通过一定时间的观测得到稳定变化时期监测点形变的速度v，将稳定变化时期的山体变形运动看作是匀速运动，可知v是一个定值，定义：
其中S(k)表示从起始监测时刻开始到当前监测时刻监测点处的累计位移量。T(k)为具有时间量纲的纵坐标值。由T(k)和t(k)绘制T-t曲线，如图3所示，通过T-t曲线可以得到切线角的表达式：
其中αk表示切线角，t(k)表示第k个监测时刻,Δt为采样时间间隔。
根据行业内的长期观测研究得到山体滑坡过程的预警模型，如表2所示。该模型综合考虑了山体滑坡监测点的变形速率V和切线角α包含的信息内容，对山体滑坡预警等级进行了5层划分，每一层中山体监测点的稳定性不同，即代表不同的山体滑坡可能性和危险性，实现了对山体滑坡的合理等级预警。
表2
为了更好的检验本文关于山体滑坡形变监测及预测预报数据处理系统的功能，使用综合考虑数据内容格式、山体滑坡监测一般形变指标的模拟数据进行数据处理系统仿真。与模拟数据相对应，下面分别对警示级别数据、警戒级别数据和警报级别数据的仿真结果进行分析讨论。
警示级别数据仿真结果如图8所示。均值方差比C1为0.2612，明显小于0.35，即均值方差比精度等级为“优”；小概率误差P1为1，明显大于0.95，即小概率误差精度等级为“优”，即该灰色理论预测模型精度等级为“优”，预测结果比较精确。经过自适应kalman滤波的监测点的精确形变数据通过灰色系统理论模型预测之后得到的多期数据dL_kalman_Grey_predict呈现上升趋势，即每隔两天滑坡监测点的形变量越来越大，预计第12期(6天后)采集到的形变量能达到23.21mm；形变速度呈现上升趋势，即滑坡监测点做加速形变，预计第12期(6天后)的形变速率能达到11.6mm/天；蠕变切线角angle呈现上升趋势，预计到第12期(6天后)能够达到79.0°；预警等级warning_results也呈现上升趋势，预计第12期(6天后)的预警等级为2级，即为“警示”等级，短期内有一定的概率出现山体滑坡灾害，与模拟的警示等级数据相符合。
警戒级别数据仿真结果如图9所示。均值方差比C1为0.0104，明显小于0.35，即均值方差比精度等级为“优”；小概率误差P1为1，明显大于0.95，即小概率误差精度等级为“优”，即该灰色理论预测模型精度等级为“优”，预测结果比较精确。经过自适应kalman滤波的监测点的精确形变数据通过灰色系统理论模型预测之后得到的多期数据dL_kalman_Grey_predict呈现上升趋势，即每隔两天滑坡监测点的形变量越来越大，预计第12期(6天后)采集到的形变量能达到52.77mm；形变速度呈现上升趋势，即滑坡监测点做加速形变，预计第12期(6天后)的形变速率能达到26.4mm/天；蠕变切线角angle呈现上升趋势，预计到第12期(6天后)能够达到84.9°；预警等级warning_results也呈现上升趋势，第9期(最近数据采集时刻)的预警等级为3级，预计第12期(6天后)的预警等级为3级，即为“警戒”等级，短期内有较大的概率出现山体滑坡灾害，与模拟的警戒等级数据相符合。
警报级别数据仿真结果如图10所示。均值方差比C1为0.0127，明显小于0.35，即均值方差比精度等级为“优”；小概率误差P1为1，明显大于0.95，即小概率误差精度等级为“优”，即该灰色理论预测模型精度等级为“优”，预测结果比较精确。经过自适应kalman滤波的监测点的精确形变数据通过灰色系统理论模型预测之后得到的多期数据dL_kalman_Grey_predict呈现上升趋势，即每隔两天滑坡监测点的形变量越来越大，预计第12期(6天后)采集到的形变量能达到60.49mm；形变速度呈现上升趋势，即滑坡监测点做加速形变，预计第12期(6天后)的形变速率能达到30.24mm/天；蠕变切线角angle呈现上升趋势，预计到第12期(6天后)能够达到85.56°；预警等级warning_results也呈现上升趋势，第9期(最近数据采集时刻)的预警等级为3级，预计第12期(6天后)的预警等级为4级，即为“警报”等级，短期内有非常大的概率出现山体滑坡灾害，与模拟的警报等级数据相符合。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
Claims (6)
1.一种地面沉降监测方法，其特征在于，包括以下步骤：
步骤一：建立随机森林模型，获取并分析气象水文、地质构造以及植被覆盖数据，对目标区域是否会发生山体滑坡进行预测，选出会发生山体滑坡的监测点；对判断得到的会发生滑坡的监测点进行重点监测，实时定位，采集地面沉降变化量，获取目标区域的数字高程数据及高精度水平数字信息，并对获取到的数据信息进行预处理；
步骤二：对步骤一中得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合；
步骤三：对步骤二中得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测；
步骤四：对步骤三中得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分，完成地面沉降监测；
在步骤一中所述的建立随机森林模型选取滑坡重点监测点，采集这些监测点处的数据信息并进行预处理，其具体做法如下：
S11、对无人机摄影测量得到的DEM数据进行分析，并结合当地历史水文气象数据，得到目标监测区域的高程、坡度、坡向、粗糙度、地形起伏度、岩土体类型、断层缓冲区、降雨量、河流缓冲区、道路缓冲区、归一化植被数据信息；
S12、使用S11得到的数据集建立随机森林模型，并测试建立的随机森林模型的预测准确率；
S13、将S11中得到的数据信息输入到S12建立起来的随机森林模型中，得到山体滑坡发生与否的初步预测结果；
S14、使用实时动态载波相位差分技术RTK进行高精度定位，精准获取监测装置的实时位置信息；
S15、使用土工带传感器采集地面沉降变化量；
S16、使用无人机拍摄获取数字高程数据、高精度水平数字信息及山体地质数据；
S17、对RTK定位接收到的经度、维度、高度数据进行中值滤波，剔除异常值；然后设置参考时刻，通过实时获得的经纬高信息计算任意时刻RTK监测装置与参考时刻相比移动的距离；
S18、土工带传感器按照网格状排列，对土工带传感器返回的数据，首先进行中值滤波来剔除异常值；然后在几何层面进行平滑，将所有行的土工带传感器取平均，将所有列的土工带传感器取平均，使用勾股定理得到针对该监测点的形变量；
S19、进行无人机飞行测量之后得到目标区域的数字高程地图DEM及高精度水平数字地图，对于不同期的摄影测量结果，通过对高精度水平数字地图作差，计算得到监测点处的水平位移量及位移方向；通过数字高程地图作差，计算得到监测点处的垂直方向的位移量，最终测量得到山体滑坡监测点的位移量。
2.根据权利要求1所述的地面沉降监测方法，其特征在于：在步骤二中所述的对步骤一中得到的多源形变监测数据使用自适应kalman滤波技术进行数据融合，其具体做法如下：
S21、对山体滑坡运动进行分析，建立监测点形变的运动数学模型；
S22、基于S21中建立的运动数学模型，使用自适应kalman滤波对RTK定位数据、土工带传感器数据、无人机摄影测量数据进行数据融合。
3.根据权利要求2所述的地面沉降监测方法，其特征在于：所述在步骤三中所述的对步骤二中得到的数据融合结果，使用灰色预测理论模型对目标监测点的形变量进行预测，其具体做法如下：
S31、选取若干期的形变监测值组成长度为n的原始观测序列；
S32、生成一次累加和序列，按照灰色预测模型参数计算公式的需要，利用选取的n-2、n-1期的数据分别生成一次累加和序列/>和/>；
S33、依次生成、/>和/>的一次累加和序列的紧邻均值序列/>、/>和/>；
S34、使用准备好的数据序列计算灰色方程的系数a和b，分别求得序列/>、/>和/>的灰色方程系数/>、/>、/>；
S35、在上述步骤中求出了序列的灰色方程系数/>，即已经得到了该序列的一次累加和序列即位置信息序列的灰色预测时间响应函数：
(1)
使用该预测时间响应函数对未来时刻的一次累加和序列即位置信息序列进行预测；对位置信息序列进行累减生成得到位移序列；
S36、分别计算序列、/>和/>的均值和方差；
S37、通过位移预测数据前8期的数据序列与数据序列/>作差求出残差序列；然后计算该残差序列的均值和方差；
S38、根据准备好的数据计算该灰色预测模型的均方差比值C和小概率误差P，根据提前设定的模型精度等级参考对该灰色理论模型的精度进行评估分析，进而判断山体滑坡灾害预测预报结果的可信度。
4.根据权利要求3所述的地面沉降监测方法，其特征在于：
在步骤四中所述的对步骤三中得到的目标监测点的形变预测值，使用蠕变切线角判据进行山体滑坡预警等级的划分；具体做法如下：
S41、建立蠕变切线角预警判据；
S42、根据行业内的长期观测研究得到山体滑坡过程的预警模型，实现了对山体滑坡的合理等级预警。
5.根据权利要求1所述的地面沉降监测方法，其特征在于：
所述步骤一包括：
通过RTK定位得到监测点的经度、维度和天线高度，则监测点位移量与竖直位移量d和水平位移量/>满足：
(2)
(3)
(4)
其中分别为在起始时刻监测点的经纬度；/>为在观测时刻监测点的经纬度；R为地球半径，取平均值/>；
基于土工带电阻值计算其形变量的计算公式如下：
(5)
*100(6)
式中，R为土工带在观测时刻的电阻测量值，单位：mΩ；是土工带在完成布置时的初始电阻值，单位：mΩ；/>是土工带在完成布置时的初始长度，单位：mm；k是比例系数；T为土壤温度，单位：/>；
首先对土工带数据进行中值滤波和平滑处理，得到监测点的纵向为位移信息y和横向为位移信息x，则监测点的形变量:
(7)
通过无人机近景摄影测量技术获得监测区域的三维点云数据、倾斜摄影处理结果；使用ArcMap对多期测量数据中的影像控制点进行匹配之后，在某一期的三维模型中加入另一期的相对三维模型，计算出前后两期摄影测量中滑坡监测点处的绝对位移量，并获得形变区域的位置和大小。
6.根据权利要求5所述的地面沉降监测方法，其特征在于，
所述步骤二具体包括：
对山体滑坡运动进行分析，建立监测点形变的运动模型，其描述方程如下：
(8)
其中，和/>分别表示形变监测点在时刻i的位置和速度；/>表示形变监测点从i时刻到i+1时刻的形变加速度；通过上述滑坡监测运动模型建立状态方程和观测方程如下：
(9)
其中，，/>，/>=/>，/>，/>是均值为0，方差为/>的高斯白噪声序列，/>是均值为0，方差为/>的高斯白噪声序列；/>为滑坡形变监测点位移、速度、加速度的观测向量；其中T为监测数据的采集周期；
基于上述建立的山体滑坡监测点运动模型，使用自适应kalman滤波对RTK定位数据、土工带传感器数据、无人机摄影测量数据进行数据融合；自适应Kalman滤波的原理通过如下方程加以描述：
时间更新方程：
(10)
(11)
状态更新方程：
(12)
=/>(13)
(14)
其中，表示观测噪声的等价协方差矩阵，观测噪声的等价权矩阵为/>；/>表示自适应因子，/>；自适应kalman滤波的关键在于自适应因子/>的构造；
自适应因子的构造采用两段函数法，两段函数模型为：
(15)
式中，c为常量，取，/>。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211390802.9A CN115638767B (zh) | 2022-11-07 | 2022-11-07 | 地面沉降监测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211390802.9A CN115638767B (zh) | 2022-11-07 | 2022-11-07 | 地面沉降监测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115638767A CN115638767A (zh) | 2023-01-24 |
| CN115638767B true CN115638767B (zh) | 2023-10-03 |
Family
ID=84949486
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211390802.9A Active CN115638767B (zh) | 2022-11-07 | 2022-11-07 | 地面沉降监测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115638767B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117870626B (zh) * | 2024-03-11 | 2024-05-07 | 中铁七局集团广州工程有限公司 | 一种用于熔岩地区的打桩浇筑过程的监测方法 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108332649A (zh) * | 2018-02-07 | 2018-07-27 | 桂林电子科技大学 | 一种滑坡形变综合预警方法及系统 |
| CN109389807A (zh) * | 2017-08-06 | 2019-02-26 | 中铁二院工程集团有限责任公司 | 蠕变型滑坡智能监测预警系统 |
| CN110568440A (zh) * | 2019-09-10 | 2019-12-13 | 四川省地质工程勘察院集团有限公司 | 一种基于DS-InSAR技术监测复杂山区形变的方法 |
| CN110986747A (zh) * | 2019-12-20 | 2020-04-10 | 桂林电子科技大学 | 一种滑坡位移组合预测方法及系统 |
| CN215006896U (zh) * | 2021-06-16 | 2021-12-03 | 深圳防灾减灾技术研究院 | 一种星地协同的边坡多风险因子联合实时监测与预警系统 |
| CN114818518A (zh) * | 2022-06-30 | 2022-07-29 | 深圳特科动力技术有限公司 | 一种陡坡防滑坡危险监测信息分析方法 |
| CN114898138A (zh) * | 2022-03-24 | 2022-08-12 | 杭州叙简科技股份有限公司 | 一种基于地形地貌及遥感数据进行森林火灾预测的方法 |
| CN115019476A (zh) * | 2022-06-09 | 2022-09-06 | 贵州大学 | 一种基于多源信息融合的滑坡时空信息监测预警方法 |
-
2022
- 2022-11-07 CN CN202211390802.9A patent/CN115638767B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109389807A (zh) * | 2017-08-06 | 2019-02-26 | 中铁二院工程集团有限责任公司 | 蠕变型滑坡智能监测预警系统 |
| CN108332649A (zh) * | 2018-02-07 | 2018-07-27 | 桂林电子科技大学 | 一种滑坡形变综合预警方法及系统 |
| CN110568440A (zh) * | 2019-09-10 | 2019-12-13 | 四川省地质工程勘察院集团有限公司 | 一种基于DS-InSAR技术监测复杂山区形变的方法 |
| CN110986747A (zh) * | 2019-12-20 | 2020-04-10 | 桂林电子科技大学 | 一种滑坡位移组合预测方法及系统 |
| CN215006896U (zh) * | 2021-06-16 | 2021-12-03 | 深圳防灾减灾技术研究院 | 一种星地协同的边坡多风险因子联合实时监测与预警系统 |
| CN114898138A (zh) * | 2022-03-24 | 2022-08-12 | 杭州叙简科技股份有限公司 | 一种基于地形地貌及遥感数据进行森林火灾预测的方法 |
| CN115019476A (zh) * | 2022-06-09 | 2022-09-06 | 贵州大学 | 一种基于多源信息融合的滑坡时空信息监测预警方法 |
| CN114818518A (zh) * | 2022-06-30 | 2022-07-29 | 深圳特科动力技术有限公司 | 一种陡坡防滑坡危险监测信息分析方法 |
Non-Patent Citations (2)
| Title |
|---|
| 几种聚类优化的机器学习方法在灵台县滑坡易发性评价中的应用;邱维蓉;吴帮玉;潘学树;唐亚明;;西北地质(01);第223-232页 * |
| 邱维蓉 ; 吴帮玉 ; 潘学树 ; 唐亚明 ; .几种聚类优化的机器学习方法在灵台县滑坡易发性评价中的应用.西北地质.2020,(01),第223-232页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115638767A (zh) | 2023-01-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108332649B (zh) | 一种滑坡形变综合预警方法及系统 | |
| Wasowski et al. | Remote sensing of landslide motion with emphasis on satellite multi-temporal interferometry applications: an overview | |
| Elkhrachy | Flash flood hazard mapping using satellite images and GIS tools: a case study of Najran City, Kingdom of Saudi Arabia (KSA) | |
| Wang et al. | GIS-based assessment of landslide susceptibility using certainty factor and index of entropy models for the Qianyang County of Baoji city, China | |
| CN110044327B (zh) | 一种基于sar数据与gnss数据的基础设施沉降监测方法及系统 | |
| CN107218923A (zh) | 基于PS‑InSAR技术的地铁沿线周边环境历史沉降风险评估方法 | |
| CN113096005B (zh) | 一种监测山体现今抬升速度的雷达时序差分干涉测量方法 | |
| Zhang et al. | Remote sensing monitoring of gullies on a regional scale: a case study of Kebai region in Heilongjiang Province, China | |
| KR102118802B1 (ko) | 무인 항공기를 이용한 하천 건천화 모니터링 방법 및 시스템 | |
| CN115638767B (zh) | 地面沉降监测方法 | |
| Halicki et al. | The accuracy of the Sentinel-3A altimetry over Polish rivers | |
| Bopche et al. | Landslide susceptibility mapping: an integrated approach using geographic information value, remote sensing, and weight of evidence method | |
| Saber et al. | Bias correction of satellite-based rainfall estimates for modeling flash floods in semi-arid regions: application to Karpuz River, Turkey | |
| Fei et al. | Dem Development and Precision Analysis for Antarctic Ice Sheet Using Cryosat‐2 Altimetry Data | |
| Herzog et al. | Capturing complex star dune dynamics—repeated highly accurate surveys combining multitemporal 3D topographic measurements and local wind data | |
| Liu et al. | Architecture planning and geo-disasters assessment mapping of landslide by using airborne LiDAR data and UAV images | |
| Mukhamedjanov et al. | The use of satellite data for monitoring rivers in the Amu Darya basin | |
| Saadatkhah et al. | Spatial patterns of precipitation, altitude and monsoon directions in Hulu Kelang area, Malaysia | |
| Hu et al. | Time-series spaceborne InSAR monitoring of permafrost deformation combined with backscattering characteristics: Case study from the Qinghai–Tibet engineering corridor | |
| Amanzio et al. | Integration of terrestrial laser scanning and GIS analysis for multi-temporal landslide monitoring: A case study of the Mont de La Saxe (Aosta Valley, NW Italy) | |
| Wang et al. | Spatial Downscaling of Remote Sensing Precipitation Data in the Beijing-Tianjin-Hebei Region | |
| Lal et al. | Assessment of yearly soil erosion of Ranchi District using RUSLE integrating geo informatics approach | |
| Peitzsch et al. | Detecting snow depth changes in avalanche path starting zones using uninhabited aerial systems and structure from motion photogrammetry | |
| Kishor et al. | Landslide hazard potential analysis using GIS, Patalganga valley, garhwal, western himalayan region of India | |
| Osman et al. | Determination of a new gravimetric geoid modelling for Sudan using the least-squares collocation technique |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |