CN117805819B - 一种基于InSAR技术的道路选线地质稳定性评估方法 - Google Patents

Abstract

本发明公开了一种基于InSAR技术的道路选线地质稳定性评估方法，属于道路工程技术领域。该方法为：获取道路选线区域的广域InSAR总体形变速率情况；去除形变数据中的离群点并计算总体分位数；根据总体分位数识别活动断裂带网格并计算活动断裂带影响区域范围；根据总体分位数识别沉降网格并计算沉降影响区域范围；根据活动断裂带和/或沉降影响区域范围的识别结果，给出绕避或进一步调查分析的建议。本发明创造性地提出了一种基于InSAR数据的活动断裂带和沉降区识别方法，根据地表形变速率识别影响道路施工建设的区域，充分利用了InSAR技术的优势，实现活动断裂带和沉降区的自动快速识别，检测效率较高。

Description

一种基于InSAR技术的道路选线地质稳定性评估方法

技术领域

本发明涉及道路工程技术领域，尤其涉及一种基于InSAR技术的道路选线地质稳定性评估方法。

背景技术

合成孔径雷达（Interferometric Synthetic Aperture Radar，InSAR）技术作为一种天基遥感手段，能够快速获取大范围的空间形变态势，监测形变过程，具有快速、连续、范围广等特点，为地表形变监测研究提供一种有效的技术手段。

受地下空间快速开发与利用，以及自然矿产不断采掘、地下水侵蚀致使地下空洞等原因，地表沉降引发的地面设施破损事件时有发生。与此同时，断裂带频繁活动，容易致使地表错动，这不仅会对地面建筑物造成严重损毁，而且时刻威胁着线性交通基础设施的运行安全。

在实际工程应用中，地下采掘活动造成的沉降区空间分布和形变现状资料较为分散，多为历史调查资料现势性不强，不能反映当前形变发展情况；同时，断裂带的活动分布和变化资料也缺少直接观测数据，无法为线性交通基础设施（如公路、铁路等）的路线选择提供有效支撑。主要存在以下问题：

（1）传统地质调查现势性不高

传统地质调查采用资料收集的方式，例如收集历史采空区的分布数据，从何判断潜在的沉降区，而历史资料往往不够全面，无法反映现势性信息、加之可能存在盗挖乱采的情况，不能客观反映采空区的分布；此外，一些稳定的采空区，并未能表现出沉降特征，因此采空区的范围并不一定和沉降区重合，而InSAR技术能够较好的反映作业区形变状态现状，具有较好的现势性。

（2）传统地面调查寻找形变特征的方法容易受局部视角的影响

采用传统测绘手段进行水准、GNSS（Global Navigation Satellite System）测量等，虽然可以获取局部精确形变信息，但由于是逐点观测，只能反映局部形变态势，很难获取宏观形变态势；此外如果通过地面调查寻找形变特征的方式，如道路裂缝和房屋裂缝等方法，总体效率较低，且仅能发现已发生了较大规模不均匀沉降的区域，在整体均匀变形区域难以发现对应特征。

因此，如何采用有效方法获取长距离线性工程周围宏观地质背景及稳定性信息，有效标识采空区域及活动断裂带，将具有良好的实际工程应用价值。当前，虽然已有将InSAR方法应用于采空区、活动断裂带监测的案例，但并无采用计算机基于InSAR成果进行断裂带、沉降区范围自动识别、宏观评估的研究和专利。

发明内容

针对以上问题，本发明设计了一种基于InSAR形变监测成果的计算机活动断裂带和沉降区识别方法，为道路选线区域的广域地质稳定性评价提供技术支撑。

本发明所提出的技术方案如下：

一种基于InSAR技术的道路选线地质稳定性评估方法：

S1：获取道路选线区域的广域InSAR总体形变速率情况；

S2：去除形变数据中的离群点并计算总体分位数；

S3：根据总体分位数识别活动断裂带网格并计算活动断裂带影响区域范围；

S4：根据总体分位数识别沉降网格并计算沉降影响区域范围；

S5：根据活动断裂带和/或沉降影响区域范围的识别结果，给出绕避或进一步调查分析的建议。

本发明所述的一种基于InSAR技术的道路选线地质稳定性评估方法，所述的步骤S1的具体步骤为：

S1.1：确定需要监测的区域范围，明确监测时段，进行InSAR卫星数据查档，获得所需InSAR数据；

S1.2：使用时间序列InSAR技术进行数据处理，获取监测区域的总体形变速率情况。

本发明所述的一种基于InSAR技术的道路选线地质稳定性评估方法，所述的步骤S2的具体步骤为：

S2.1：将监测点矢量数据转为栅格数据，若形变监测结果为栅格，则跳过该步骤；

S2.2：采用标准差法去除形变速率结果中的离群点；假设全部监测点位形变速率集合为X，平均值为μ，样本标准差为σ，使用以下公式来计算阈值从而去除离群值：

；

其中n是X中点的数量，T _s 为离群阈值下界，T _l 为离群阈值上界，形变速率小于T _s 或大于T _l 的值被定义为离群值，k为常数，k值设定为通过离群点值去除后，剩余点的数量≥99.9%的值；按上述方法进行离群点位的去除，形成去除离群点后的形变监测点，使用V _Raster 表示去除离群点后的栅格数据；

将V _Raster 的85%、20%、15%、10%分位数表示为V_0.85、V_0.2、V_0.15、V_0.1：

；

其中，Q _0.85为V _Raster 的85%分位数，Q _0.2为V _Raster 的20%分位数，Q _0.15为V _Raster 的15%分位数，Q _0.1为V _Raster 的10%分位数。

本发明所述的一种基于InSAR技术的道路选线地质稳定性评估方法，所述的步骤S3的具体操作为：

S3.1：定义覆盖监测区域的500m×500m网格为G ₅₀₀ ，按下述公式计算各个网格的形变速率范围、形变速率平均值：

；

其中，r1_k表示与第k个网格区域计算输出值，V _kpq 为形变速率，p为网格的x坐标，q为网格的y坐标，max(V _kpq )表示网格R _k 中的最大值，min(V _kpq )表示网格R _k 中的最小值，mean函数为取平均值函数，m1_k表示网格R _k 的平均速率值；

S3.2：定义r1的75%分位数为提取阈值：

；

其中，Q _0.75为r1的75%分位数；

S3.3：定义满足条件：①r1_k＞T_0.75且②存在距离网格k距离≤4的网格w，满足条件m1_w>V_0.85，则f为断裂带网格，用G _500fracture 表示；若G _500fracture 不为空则进入后续步骤；

S3.4：定义覆盖监测区域的200m×200m网格为G ₂₀₀ ，计算网格内形变速率的平均值，计算公式如下：

；

其中，mean函数为取平均值函数，V _kpq 为形变速率，p为网格的x坐标，q为网格的y坐标；

S3.5：定义与G _500fracture 空间重叠的G ₂₀₀ 网格为起算断裂带网格，使用G _200fracture 表示，采用扩散算法，以G _200fracture 为扩散中心，若相邻网格满足条件1：①m2 _k >V _0.85或②m2 _k <V _0.15，则标记为G _200fracture ，并继续向邻域扩散，直至邻域无网格满足条件1为止，得到G _200fracture 分布；

S3.6：按行遍历G _200fracture 网格所覆盖区域，定义每行最左侧、最右侧的网格为边界网格，如果仅有一个网格，则该网格为边界网格，依据空间关系填补边界内空网格，形成无内空洞的面，最终形成的连续区域为活动断裂带影响区域Z _fracture 。

本发明所述的一种基于InSAR技术的道路选线地质稳定性评估方法，所述的步骤S4的具体操作为：

S4.1：定义覆盖区域的100m×100m网格为G ₁₀₀ ，按空间范围排除已识别得到的Z _fracture 网格，计算网格内速率的平均值，得到速率值；

S4.2：标记速率值≤V_0.1的网格为G _sink ，以G _sink 为扩散中心，若相邻网格满足条件1：速率值≤V_0.2，则标记为G _sink ，并继续向邻域扩散，直至邻域无网格满足条件1为止，得到G _sink 分布；若沉降区面积过小，即单个G _sink 的网格总数≤4时，则取消标记；

S4.3：按行遍历G _sink 网格所覆盖区域，定义每行最左侧、最右侧的网格为边界网格；如果仅有一个网格，则该网格为边界网格；依据空间关系填补边界内空网格形成无内空洞的面，最终形成的连续区域为沉降影响区域Z _sink 。

本发明所述的一种基于InSAR技术的道路选线地质稳定性评估方法，所述的步骤S5的具体操作为：

S5.1：对于活动断裂带影响区域，公路路线绕避；对于沉降影响区域，公路路线绕避或结合专题资料开展更进一步的调查分析，确定其稳定性。

本发明基于InSAR形变监测成果，利用形变速率的空间分布信息，识别得到沉降影响区域范围和活动断裂带影响区域范围，实现了沉降区、活动断裂带的自动快速识别，并根据影响范围最终给出稳定性评估意见，确保道路选线合理。

附图说明

图1为本发明的技术流程图；

图2为形变速率分级设色图；

图3为500m网格的活动断裂带识别结果图；

图4为活动断裂带影响区域图；

图5为沉降影响区域图。

具体实施方式

下面结合试验例及具体实施方式对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。

如图1所示，本发明所述的一种基于InSAR技术的道路选线地质稳定性评估方法，包括以下步骤：

S1：获取道路选线区域的广域InSAR总体形变速率情况，具体包括如下步骤：

S1.1：根据公路走廊带确定需要监测的区域范围，明确监测时段，进行InSAR卫星数据查档，并结合实际项目预算及数据存档情况确定数据源，获得所需InSAR数据；

S1.2：使用时间序列InSAR技术进行数据处理，获取监测区域的总体形变速率情况。获取12期及以上的监测区域内的时间序列InSAR结果，可以是PSInSAR（PersistentScatterer InSAR）、SBASInSAR（Small Baseline Subset InSAR）、Stamps(StanfordMethod for Persistent Scatterers)等技术，技术的选择应当保证监测结果的空间点密度尽可能高。时间序列InSAR结果主要包含监测点空间位置、形变速率大小和逐期形变量大小，三个维度的信息。本技术主要利用监测点空间位置和形变速率大小。

S2：去除形变数据中的离群点并计算总体分位数，具体包括如下步骤：

S2.1：将监测点矢量数据转为栅格数据。栅格大小参考所使用的卫星雷达数据空间分辨率，一般定义为所用雷达数据的最高空间分辨率。若形变监测结果为栅格，则跳过该步骤。形变速率转换为栅格数据后，分级设色结果如图2所示；

S2.2：采用标准差法去除形变速率结果中的离群点。假设全部监测点位形变速率集合为X，平均值为μ，样本标准差为σ，使用以下公式来计算阈值来检测离群值：

；

其中n是X中点的数量，T _s 为离群阈值下界，T _l 为离群阈值上界，形变速率小于T _s 或大于T _l 的值被定义为离群值，k为常数，k值设定为通过离群点值去除后，剩余点的数量≥99.9%的值。按上述方法进行离群点位的去除，形成去除离群点后的形变监测点，使用V _Raster 表示去除离群点后的栅格数据；

S2.3：将V _Raster 的85%、20%、15%、10%分位数表示为V_0.85、V_0.2、V_0.15、V_0.1：

；

其中，Q _0.85为V _Raster 的85%分位数，Q _0.2为V _Raster 的20%分位数，Q _0.15为V _Raster 的15%分位数，Q _0.1为V _Raster 的10%分位数。

S3：根据总体分位数识别活动断裂带网格并计算活动断裂带影响区域范围，具体包括如下步骤：

S3.1：定义覆盖监测区域的500m×500m网格为G ₅₀₀ ，按下述公式计算各个网格的形变速率范围、形变速率平均值：

；

其中，r1_k表示与第k个网格区域计算输出值，V _kpq 为形变速率，p为网格的x坐标，q为网格的y坐标，max(V _kpq )表示网格R _k 中的最大值，min(V _kpq )表示网格R _k 中的最小值，mean函数为取平均值函数，m1_k表示网格R _k 的平均速率值；

S3.2：定义r1的75%分位数为提取阈值：

；

其中，Q _0.75为r1的75%分位数；

S3.3：定义满足条件：①r1_k＞T_0.75且②存在距离网格k距离≤4的网格w，满足条件m1_w>V_0.85(dis(w,k)≤4)为断裂带网格，用G _500fracture 表示，如图3所示。若G _500fracture 不为空则进入后续步骤；

S3.4：定义覆盖监测区域的200m×200m网格为G ₂₀₀ ，计算网格内形变速率的平均值，计算公式如下：

；

其中，mean函数为取平均值函数，V _kpq 为形变速率，p为网格的x坐标，q为网格的y坐标；

S3.5：定义与G _500fracture 空间重叠的G ₂₀₀ 网格为起算断裂带网格，使用G _200fracture 表示，采用扩散算法，以G _200fracture 为扩散中心，若相邻网格满足条件1：①m2 _k >V _0.85或②m2 _k <V _0.15，则标记为G _200fracture ，并继续向邻域扩散，直至邻域无网格满足条件1为止，得到G _200fracture 分布；

S3.6：按行遍历G _200fracture 网格所覆盖区域，定义每行最左侧、最右侧的网格为边界网格。如果仅有一个网格，则该网格为边界网格。依据空间关系填补边界内空网格，形成无内空洞的面，最终形成的连续区域为活动断裂带影响区域Z _fracture ，如图4所示。

S4：根据总体分位数识别沉降网格并计算沉降影响区域范围，具体包括如下步骤：

S4.1：定义覆盖区域的100m×100m网格为G ₁₀₀ ，按空间范围排除已识别得到的Z _fracture 网格，计算网格内速率的平均值，得到速率值；

S4.2：标记速率值≤V_0.1的网格为G _sink ，以G _sink 为扩散中心，若相邻网格满足条件1：速率值≤V_0.2，则标记为G _sink ，并继续向邻域扩散，直至邻域无网格满足条件1为止，得到G _sink 分布；若沉降区面积过小，即单个G _sink 的网格总数≤4时，则取消标记；

S4.3：按行遍历G _sink 网格所覆盖区域，定义每行最左侧、最右侧的网格为边界网格；如果仅有一个网格，则该网格为边界网格；依据空间关系填补边界内空网格形成无内空洞的面，最终形成的连续区域为沉降影响区域Z _sink ，如图5所示。

S5：根据活动断裂带和/或沉降影响区域范围的识别结果，给出绕避或进一步调查分析的建议，具体包括如下步骤：

S5.1：对于活动断裂带影响区域，公路路线绕避；对于沉降影响区域，公路路线绕避或结合专题资料开展更进一步的调查分析，确定其稳定性。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于InSAR技术的道路选线地质稳定性评估方法，其特征在于：

S1：获取道路选线区域的广域InSAR总体形变速率情况；

S2：去除形变数据中的离群点并计算总体分位数；

S3：根据总体分位数识别活动断裂带网格并计算活动断裂带影响区域范围；

S4：根据总体分位数识别沉降网格并计算沉降影响区域范围；

S5：根据活动断裂带和/或沉降影响区域范围的识别结果，给出绕避或进一步调查分析的建议；

所述的步骤S3的具体步骤如下：

S3.1：定义覆盖监测区域的500m×500m网格为G ₅₀₀ ，按下述公式计算各个网格的形变速率范围、形变速率平均值：

；

其中，r1_k表示与第k个网格区域计算输出值，V _kpq 为形变速率，p为网格的x坐标，q为网格的y坐标，max(V _kpq )表示网格R _k 中的最大值，min(V _kpq )表示网格R _k 中的最小值，mean函数为取平均值函数，m1_k表示网格R _k 的平均速率值；

S3.2：定义r1的75%分位数为提取阈值：

；

其中，Q _0.75为r1的75%分位数；

S3.3：定义满足条件：①r1_k＞T_0.75且②存在距离网格k距离≤4的网格w，满足条件m1_w>V_0.85，则f为断裂带网格，用G _500fracture 表示；若G _500fracture 不为空则进入后续步骤；

其中，V_0.85为所述步骤S2所得总体分位数中的85%分位数；

S3.4：定义覆盖监测区域的200m×200m网格为G ₂₀₀ ，计算网格内形变速率的平均值，计算公式如下：

；

其中，mean函数为取平均值函数，V _kpq 为形变速率，p为网格的x坐标，q为网格的y坐标；

S3.5：定义与G _500fracture 空间重叠的G ₂₀₀ 网格为起算断裂带网格，使用G _200fracture 表示，采用扩散算法，以G _200fracture 为扩散中心，若相邻网格满足条件1：①m2 _k >V _0.85或②m2 _k <V _0.15，则标记为G _200fracture ，并继续向邻域扩散，直至邻域无网格满足条件1为止，得到G _200fracture 分布；其中，V_0.15为所述步骤S2所得总体分位数中的15%分位数；

S3.6：按行遍历G _200fracture 网格所覆盖区域，定义每行最左侧、最右侧的网格为边界网格，如果仅有一个网格，则该网格为边界网格，依据空间关系填补边界内空网格，形成无内空洞的面，最终形成的连续区域为活动断裂带影响区域Z _fracture 。

2.根据权利要求1所述的一种基于InSAR技术的道路选线地质稳定性评估方法，其特征在于：所述的步骤S1的具体步骤如下：

S1.1：确定需要监测的区域范围，明确监测时段，进行InSAR卫星数据查档，获得所需InSAR数据；

S1.2：使用时间序列InSAR技术进行数据处理，获取监测区域的总体形变速率情况。

3.根据权利要求2所述的一种基于InSAR技术的道路选线地质稳定性评估方法，其特征在于：所述的步骤S2的具体步骤如下：

S2.1：将监测点矢量数据转为栅格数据，若形变监测结果为栅格，则跳过该步骤；

S2.2：采用标准差法去除形变速率结果中的离群点；假设全部监测点位形变速率集合为X，平均值为μ，样本标准差为σ，使用以下公式来计算阈值从而去除离群值：

；

其中n是X中点的数量，T _s 为离群阈值下界，T _l 为离群阈值上界，形变速率小于T _s 或大于T _l 的值被定义为离群值，k为常数，k值设定为通过离群点值去除后，剩余点的数量≥99.9%的值；按上述方法进行离群点位的去除，形成去除离群点后的形变监测点，使用V _Raster 表示去除离群点后的栅格数据；

S2.3：将V _Raster 的85%、20%、15%、10%分位数表示为V_0.85、V_0.2、V_0.15、V_0.1：

；

其中，Q _0.85为V _Raster 的85%分位数，Q _0.2为V _Raster 的20%分位数，Q _0.15为V _Raster 的15%分位数，Q _0.1为V _Raster 的10%分位数。

4.根据权利要求3所述的一种基于InSAR技术的道路选线地质稳定性评估方法，其特征在于：所述的步骤S4的具体步骤如下：

S4.1：定义覆盖区域的100m×100m网格为G ₁₀₀ ，按空间范围排除已识别得到的Z _fracture 网格，计算网格内速率的平均值，得到速率值；

S4.2：标记速率值≤V_0.1的网格为G _sink ，以G _sink 为扩散中心，若相邻网格满足条件1：速率值≤V_0.2，则标记为G _sink ，并继续向邻域扩散，直至邻域无网格满足条件1为止，得到G _sink 分布；若沉降区面积过小，即单个G _sink 的网格总数≤4时，则取消标记；

S4.3：按行遍历G _sink 网格所覆盖区域，定义每行最左侧、最右侧的网格为边界网格；如果仅有一个网格，则该网格为边界网格；依据空间关系填补边界内空网格形成无内空洞的面，最终形成的连续区域为沉降影响区域Z _sink 。

5.根据权利要求4所述的一种基于InSAR技术的道路选线地质稳定性评估方法，其特征在于：所述的步骤S5的具体步骤如下：

S5.1：对于活动断裂带影响区域，公路路线绕避；对于沉降影响区域，公路路线绕避或结合专题资料开展更进一步的调查分析，确定其稳定性。