CN117784075A - 一种融合雷达遥感和空间分析的岩溶发育强度分级方法 - Google Patents

Abstract

本发明涉及地质勘探技术领域，具体为一种融合雷达遥感和空间分析的岩溶发育强度分级方法，包括：S1、获取工程区卫星合成孔径雷达影像和航空激光雷达点云数据；S2、开展时序形变解算处理和点云精细化分类；S3、开展基于精细数字高程模型的落水洞自动定位与体积参数提取；S4、获取落水洞中心50m缓冲区范围内5年以上的平均形变速率；S5、根据落水洞坐标，获取岩溶发育密度等级初步评价指标；S6、融合落水洞的体积和邻近范围平均形变速率，基于GIS空间分析获取岩溶发育等级综合评价结果。本发明有效解决了植被覆盖下岩溶落水洞隐蔽性高，几何参数、形变信息获取的难题，为喀斯特山区隧道工程勘察与水害防治工作提供了一种有效的定量研究技术手段。

Description

一种融合雷达遥感和空间分析的岩溶发育强度分级方法

技术领域

本发明涉及地质勘探技术领域，尤其涉及一种融合雷达遥感和空间分析的岩溶发育强度分级方法。

背景技术

突泥涌水是喀斯特山区隧道工程建设中的重大安全隐患之一，因此岩溶不良地质的勘察与分级评价对于保障工程安全具有重要意义。为了查明岩溶发育的程度，辅助隧道工程选址，传统方法多采用人工现场调查、钻探的技术手段，根据线岩溶率、遇洞隙率等参量加以评价。然而岩溶发育具有严重的不均匀性和隐蔽性，上述单点化或抽样调查的方法存在效率低且盲区多的问题，最终导致评价结果与实际不符。

落水洞是地下径流的重要补给源，其空间分布、体积规模、邻近区域沉降速率等信息直观反映了岩溶的发育强度。但落水洞表面常被茂密的喜水植被覆盖，极大增加了勘察难度。首先，人工调查依赖于资料收集和走访，主观性较强，导致大量隐患被遗漏。其次，光学遥感影像无法穿透植被，亦难以获取小尺度、高隐蔽性落水洞的痕迹。最后，即使综合上述两种方法，投入大量人力物力勘察到少量落水洞，依然缺乏有效手段探明其体积规模和沉降速率。

因此，亟需研发一种基于落水洞空间分布、体积和邻近区域沉降信息的岩溶发育强度分级方法，为隧道工程选址提供参考。

发明内容

本发明的目的在于，解决现有技术的不足之处，提出一种融合雷达遥感和空间分析的岩溶发育强度分级方法，该方法通过搜集工程区合成孔径雷达和激光雷达数据并进行处理，获取形变速率场和去植被影响的数字高程模型，进而提取落水洞空间位置、体积和邻近区域形变参数，融合GIS空间分析方法获得工程区岩溶发育强度分级结果。

一种融合雷达遥感和空间分析的岩溶发育强度分级方法，包括如下步骤：

S1，获取工程区卫星合成孔径雷达SAR影像和航空激光雷达点云数据；

S2，开展时序形变解算处理和点云精细化分类：对S1中获取的SAR影像进行干涉处理，获取工程区形变速率场Rate；同时开展点云精细化分类，获得去除植被影响的数字高程模型DEM；

S3，开展基于去植被影响数字高程模型的落水洞自动定位与体积参数提取：利用S2中获取的DEM开展水文分析，自动获取岩溶洼地，并提取落水洞空间位置和体积参数；

S4，获取落水洞中心50m缓冲区范围内5年以上的平均形变速率；

S5，根据落水洞坐标，获取岩溶发育密度等级初步评价结果；

S6，根据落水洞体积和邻近区域的平均形变速率，融合GIS空间分析技术，获取地表岩溶发育等级综合评价结果。

优选的，S1中获取的数据包括以下方面：

A，划定石灰岩，白云岩等可溶岩地层的工程区作为雷达遥感数据获取区域；

B，根据工程区地形起伏和坡体朝向，确定适宜的卫星雷达数据源，同时搜集相应的地形数据，精密轨道数据，大气校正数据；

C，根据植被郁闭度，设定仿地飞行正交航线，获取高密度航空激光雷达点云数据。

优选的，S3包括以下子步骤：

S3.1，利用DEM提取山脊，沟谷等地形特征线，并获取岩溶汇水洼地；

S3.2，利用DEM数据构建山体阴影和坡度图，设置透明度进行融合制图；

S3.3，利用S3.1中获取的洼地边界矢量叠对S3.2的融合图件进行裁剪，获得兴趣区域ROI；

S3.4，在ROI中开展近似圆形或椭圆形的多边形目标检索，将边界灰度值差异超过20，且地形深度超过2m的凹陷状目标初步判定为落水洞；

S3.5，通过点云断面，高分影像对落水洞目标进行人工复核，查找遗漏并排除错误检测，并将最终的多边形结果记为落水洞向量Hole_i, i=1,…,n；

S3.6，记录向量Hole中各元素最小包闭矩形几何中心的空间坐标，记为向量（X_i,Y_i,Z_i）, i=1,…,n；利用边界矢量填方法计算出各落水洞的体积，记为向量V_i, i=1,…,n。

优选的，S4包括以下子步骤：

S4.1，利用S3中获取的落水洞中心坐标（X_i,Y_i,Z_i）, i=1,…,n生成半径50m的圆形缓冲区；

S4.2，利用缓冲区对S2中获取的形变速率场Rate进行裁剪，获得对应缓冲区元素的形变速率结果；

S4.3，计算裁剪后形变速率场各元素的形变算术平均值，并将其作为落水洞hi的平均形变速率，记为Rate_mean_i,i=1,…,n。

优选的，S5包括以下子步骤：

S5.1，将落水洞几何中心的空间坐标输入GIS系统，以工程区左上角为原点，1km×1km为搜索边长，通过点密度分析法确定岩溶发育密度等级初步评价指标；

S5.2，根据落水洞发育密度，确定岩溶发育等级，包括：岩溶强烈发育，落水洞密度＞5个km²；岩溶中等发育，1＜落水洞密度＜5个km²；岩溶弱发育，落水洞密度＜1个km²三个等级。

优选的，S6包括以下子步骤：

S6.1，在密度发育等级初步评价的基础上，进一步确定岩溶发育详细分级标准：Ⅰ级[0~0.3]：溶蚀影响较小区域；Ⅱ级（0.3~0.6]：溶蚀影响程度中等区域；Ⅲ级（0.6~1]：溶蚀影响严重区域；

S6.2，将落水洞体积和形变速率分别进行归一化处理，记为Norm_V_i, i=1,…,n和Norm_rate_i, i=1,…,n；

S6.3，为落水洞体积和平均形变速率设置相应的融合权重，记为P_V和P_Rate；

S6.4，基于归一化的落水洞体积，形变速率，以及融合权重，计算工程区落水洞溶蚀强度加权融合指标Index_i, i=1,…,n；

S6.5，利用S6.4获得的Index，结合GIS空间分析技术，计算获得工程区的岩溶发育等级综合评价结果；

S6.6，以综合评价发育等级为准，绘制工程区溶蚀影响空间分布图。

本发明的优点及技术效果在于：

本发明的一种融合雷达遥感和空间分析的岩溶发育强度分级方法，针对植被覆盖下岩溶落水洞隐蔽性高，调查低效的现状，通过航天、航空雷达遥感技术，获取了落水洞的空间位置、体积、邻近区域形变，融合GIS空间分析法获得工程区岩溶发育强度分级结果。现场复核及工程实践表明，该方法理论严谨，可靠性高，分级结果可以为隧道工程勘察和水害防治提供参考。

附图说明

图1为本发明的方法流程图；

图2为本发明的缓冲区形变速率裁剪示意图；

图3为本发明的落水洞自动识别与复核示意图；

图4为本发明的岩溶强烈发育区的综合分级结果示意图。

具体实施方式

下面详细描述本发明的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。

一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其获取方法包括如下步骤：

S1）获取工程区卫星合成孔径雷达（SAR）影像和航空激光雷达点云数据：

结合地质图、卫星影像图资料，划定石灰岩、白云岩等可溶岩地层的工程区作为雷达遥感数据获取区域；

②根据工程区地形起伏和坡体朝向，优选近5年或更长时序的C波段或L波段卫星合成孔径雷达影像，确保采样间隔至少在1景/月，即干涉处理至少60景数据。同时搜集相应区域的地形资料、精密轨道数据、大气校正数据；

③根据植被郁闭度，设定高重叠率的仿地飞行航线，获取高密度航空激光雷达点云数据，确保激光的植被穿透效果。

S2）开展时序形变解算处理和点云精细化分类：对S1）获取的SAR影像进行干涉处理，获取工程区形变速率场Rate。同时开展点云精细化分类，获得去除植被影响的数字高程模型（DEM）；

基于搜集的SAR影像和相关数据，利用小基线集或分布式散射体干涉技术，剔除低相干区域，获取工程区近5年或以上的年平均形变速率解算结果；

②基于高密度点云数据，通过种子点选取和逐级加密检索提取地面点，利用空间插值获得高分辨率DEM。

S3）开展基于去植被影响数字高程模型的落水洞自动定位与体积参数提取：利用S2）获取的DEM开展水文分析，自动获取岩溶洼地，并提取落水洞空间位置和体积参数：

利用DEM提取山脊、沟谷等地形特征线，并获取岩溶汇水洼地；

②利用DEM数据构建山体阴影和坡度图。将坡度图设置30%透明度，叠加至山体阴影图层上方，进行融合制图；

③利用①获取的洼地边界矢量叠对②的融合图件进行裁剪，获得兴趣区域ROI；

④在ROI中开展近似圆形或椭圆形的多边形目标检索，将边界灰度值差异超过20，且地形深度超过2m的凹陷状目标初步定为落水洞；

⑤通过点云断面、高分影像对落水洞目标进行人工复核，查找遗漏并排除错误检测，并将最终的多边形结果记为落水洞向量Hole_i, i=1,…,n；

⑥记录向量Hole中各元素最小包闭矩形几何中心的空间坐标，记为向量（X_i,Y_i,Z_i）, i=1,…,n；利用边界矢量填方法计算出各落水洞的体积，记为向量V_i, i=1,…,n。

S4）获取落水洞中心50m缓冲区范围内5年以上的平均形变速率：

利用S3）获取的落水洞中心坐标（X_i,Y_i,Z_i）, i=1,…,n生成50m圆形缓冲区；

②利用缓冲区对S2）获取的形变速率场Rate进行裁剪，获得对应缓冲区元素的形变速率结果；

③计算裁剪后形变速率场各元素的形变算术平均值，并将其作为落水洞Hole_i, i=1,…,n的平均形变速率，记为Rate_mean_i,i=1,…,n。统计中忽略地形阴影和失相干导致的无效形变Nan值，计算公式如下，其中r表示缓冲区半径，j表示缓冲区内InSAR相干点个数，k表示InSAR矩形相干像元的边长：

j=⌈/⌉（1）

Rate_mean_i = (/j, i=1,…,n （2）

S5）根据落水洞坐标，获取岩溶发育密度等级初步评价结果：

①将落水洞几何中心的空间坐标输入GIS系统，以工程区左上角为原点，1km×1km为搜索边长，通过点密度分析法确定岩溶发育密度等级初步评价指标；

②根据落水洞发育密度，确定岩溶强烈发育（落水洞密度＞5个km²）、岩溶中等发育（1＜落水洞密度＜5个km²）和岩溶弱发育（落水洞密度＜1个km²）三个等级。

S6）根据落水洞体积和邻近区域的平均形变速率，融合GIS空间分析技术，获取地表岩溶发育等级综合评价结果：

在密度发育等级初步评价的基础上，进一步确定岩溶发育详细分级标准：Ⅰ级[0~0.3]：溶蚀影响较小区域；Ⅱ级（0.3~0.6]：溶蚀影响程度中等区域；Ⅲ级（0.6~1]：溶蚀影响严重区域；

②将落水洞体积和形变速率分别进行归一化处理，记为Norm_V_i, i=1,…,n和Norm_rate_i, i=1,…,n；

本方法采用了最值归一化来解决体积和形变量纲不一致的问题。归一化体积计算公式如下，其中V_i表示第i个落水洞的体积；V_i表示第i个落水洞的体积_；Norm_V_i为第i个落水洞的归一化体积；V_max和V_min对应落水洞体积的最大值和最小值：

Norm_V_i=（V_i-V_min）/（V_max-V_min） （3）

归一化形变速率计算公式如下，其中：

③为落水洞体积和平均形变速率设置相应的融合权重，记为P_V和P_Rate；

综合考虑落水洞体积和形变速率在岩溶汇水中的影响，本方法采用了如下定权方案，其中P_v表示体积权重，P_rate表示形变速率权重：

P_v=2/3,P_rate=1/3 （5）

④基于归一化的落水洞体积、形变速率，以及融合权重，计算工程区落水洞溶蚀强度加权融合指标Index_i, i=1,…,n，计算公式如下：

Index_i=Norm_V_i·P_v+Norm_rate_i·P_rate （6）

综上，得出岩溶发育加权融合指标如下：

表1 岩溶发育加权融合指标计算

⑤利用步骤④获得的加权融合指标Index_i, i=1,…,n，结合GIS空间分析技术，计算获得岩溶综合评价结果（9个级别），详见表2；

表2 岩溶发育等级综合评价

⑥以综合评价发育等级为准，绘制工程区溶蚀影响空间分布图。通过为溶蚀影响图设置相应的颜色带，从而为隧道工程选址提供参考，具体建议如下：

隧道工程宜优先选择穿越岩溶弱发育地区，其中以弱-Ⅰ级为最佳，弱-Ⅱ级和弱-Ⅲ级次之；当无法绕避，需穿越岩溶中等发育区时，宜优先穿越中等-Ⅰ级，中等-Ⅱ级和中等-Ⅲ级次之，并在设计阶段采取隧道衬砌加固、增强排导等措施；岩溶强烈发育区建议绕避。

最后，本发明的未尽述之处均采用现有技术中的成熟产品及成熟技术手段。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的实施例或示例中。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同限定。

Claims (6)

1.一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其特征在于，包括如下步骤：

S1，获取工程区卫星合成孔径雷达SAR影像和航空激光雷达点云数据；

S2，开展时序形变解算处理和点云精细化分类：对S1中获取的SAR影像进行干涉处理，获取工程区形变速率场Rate；同时开展点云精细化分类，获得去除植被影响的数字高程模型DEM；

S3，开展基于去植被影响数字高程模型的落水洞自动定位与体积参数提取：利用S2中获取的DEM开展水文分析，自动获取岩溶洼地，并提取落水洞空间位置和体积参数；

S4，获取落水洞中心50m缓冲区范围内5年以上的平均形变速率；

S5，根据落水洞坐标，获取岩溶发育密度等级初步评价结果；

S6，根据落水洞体积和邻近区域的平均形变速率，融合GIS空间分析技术，获取地表岩溶发育等级综合评价结果。

2.根据权利要求1所述的一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其特征在于：所述S1中获取的数据包括以下方面：

A，划定石灰岩，白云岩等可溶岩地层的工程区作为雷达遥感数据获取区域；

B，根据工程区地形起伏和坡体朝向，确定适宜的卫星雷达数据源，同时搜集相应的地形数据，精密轨道数据，大气校正数据；

C，根据植被郁闭度，设定仿地飞行正交航线，获取高密度航空激光雷达点云数据。

3.根据权利要求1所述的一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其特征在于：所述S3包括以下子步骤：

S3.1，利用DEM提取山脊，沟谷等地形特征线，并获取岩溶汇水洼地；

S3.2，利用DEM数据构建山体阴影和坡度图，设置透明度进行融合制图；

S3.3，利用S3.1中获取的洼地边界矢量叠对S3.2的融合图件进行裁剪，获得兴趣区域ROI；

S3.4，在ROI中开展近似圆形或椭圆形的多边形目标检索，将边界灰度值差异超过20，且地形深度超过2m的凹陷状目标初步判定为落水洞；

S3.5，通过点云断面，高分影像对落水洞目标进行人工复核，查找遗漏并排除错误检测，并将最终的多边形结果记为落水洞向量Hole_i, i=1,…,n；

S3.6，记录向量Hole中各元素最小包闭矩形几何中心的空间坐标，记为向量（X_i,Y_i,Z_i）, i=1,…,n；利用边界矢量填方法计算出各落水洞的体积，记为向量V_i, i=1,…,n。

4.根据权利要求1所述的一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其特征在于：所述S4包括以下子步骤：

S4.1，利用S3中获取的落水洞中心坐标（X_i,Y_i,Z_i）, i=1,…,n生成半径50m的圆形缓冲区；

S4.2，利用缓冲区对S2中获取的形变速率场Rate进行裁剪，获得对应缓冲区元素的形变速率结果；

S4.3，计算裁剪后形变速率场各元素的形变算术平均值，并将其作为落水洞hi的平均形变速率，记为Rate_mean_i,i=1,…,n。

5.根据权利要求1所述的一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其特征在于：所述S5包括以下子步骤：

S5.1，将落水洞几何中心的空间坐标输入GIS系统，以工程区左上角为原点，1km×1km为搜索边长，通过点密度分析法确定岩溶发育密度等级初步评价指标；

S5.2，根据落水洞发育密度，确定岩溶发育等级，包括：岩溶强烈发育，落水洞密度＞5个km²；岩溶中等发育，1＜落水洞密度＜5个km²；岩溶弱发育，落水洞密度＜1个km²三个等级。

6.根据权利要求1所述的一种融合雷达遥感和空间分析的岩溶发育强度分级方法，其特征在于：所述S6包括以下子步骤：

S6.1，在密度发育等级初步评价的基础上，进一步确定岩溶发育详细分级标准：Ⅰ级[0~0.3]：溶蚀影响较小区域；Ⅱ级（0.3~0.6]：溶蚀影响程度中等区域；Ⅲ级（0.6~1]：溶蚀影响严重区域；

S6.2，将落水洞体积和形变速率分别进行归一化处理，记为Norm_V_i, i=1,…,n和Norm_rate_i, i=1,…,n；

S6.3，为落水洞体积和平均形变速率设置相应的融合权重，记为P_V和P_Rate；

S6.4，基于归一化的落水洞体积，形变速率，以及融合权重，计算工程区落水洞溶蚀强度加权融合指标Index_i, i=1,…,n；

S6.5，利用S6.4获得的Index，结合GIS空间分析技术，计算获得工程区的岩溶发育等级综合评价结果；

S6.6，以综合评价发育等级为准，绘制工程区溶蚀影响空间分布图。