CN110244298A - 一种InSAR数据升降轨联合滑坡分析方法 - Google Patents

Abstract

本发明公开了合成孔径雷达干涉技术领域的一种InSAR数据升降轨联合滑坡分析方法，包括以下步骤：S1、挑选目标区域最合适的时间序列SAR影像，基于PSP算法生产出该区域InSAR升降轨数据，选通过升降轨数据插值出经纬网格的形变速率；S2、读取对应范围的DEM，运用卷积计算的方法求出栅格点对应的坡度和坡向值，通过插值的方法将经纬网格与升降轨数据的经纬网格统一；S3、通过升降轨速率和坡度坡向，将升降轨形变速率分解为水平形变速率和垂直形变速率，得到滑坡的水平形变量和垂直形变量；S4、挑选滑坡上局部形变最大值作为初始化隐患点，获取沿坡面的剖面线，并对其中的三条剖面线1、2、3进行分析，得到沿剖面线的坡面水平和垂直方向的滑动情况。

Description

一种InSAR数据升降轨联合滑坡分析方法

技术领域

本发明涉及合成孔径雷达干涉技术领域，具体为一种InSAR数据升降轨联合滑坡分析方法。

背景技术

坡是指斜坡上的土体或者岩体，受河流冲刷、地下水活动、雨水浸泡、地震及人工切坡等因素影响，在重力作用下，沿着一定的软弱面或者软弱带，整体地或者分散地顺坡向下滑动的自然现象，引起建筑沉降的因素很多，包括自然因素和人为因素。自然因素包括地形、岩性、地质构造、以及地震，暴雨等诱发因素。人为原因包括修建公路，铁路，建筑和水库等。

滑坡常常给工农生产以及人民生命财产造成巨大损失。滑坡对乡村最主要的危害是摧毁农田、房舍、道路、农业机械、水利水电设施等；对城镇的危害常常是砸埋房屋，伤亡人畜，摧毁工厂、学校、机关单位、停水停电等；所以滑坡的监测非常重要。

常规沉降监测技术包括采用水准仪、测距仪、全站仪等常规测量仪器测定点的形变值，其优点是：

(1)能够提供变形体整体的变形状态；

(2)适用于不同的监测精度、不同形式的变形体和不同的监测环境；

(3)可以提供绝对变形信息。

滑坡面积大、范围广，利用水准测量和倾斜测量进行监测，需要耗费大量的人力物力；传统监测，一般测量的是当下绝对值，无法对未来做出预期。

InSAR技术能够高效率、高质量监测滑坡，InSAR技术通过雷达卫星向目标区域发射微波，然后接收目标反射的回波，得到同一目标区域成像的SAR 复图像对，若复图像对之间存在相干条件，SAR复图像对共轭相乘可以得到干涉图，根据干涉图的相位值，得出两次成像中微波的路程差，从而计算出目标地区的地形、地貌以及表面的微小变化，精度达到毫米级。而且由于卫星可以进行升降轨联合观测，使得InSAR可以同时观测到水平和垂直方向的形变。

InSAR监测滑坡有如下优点：

(1)精度高；

(2)监测范围广；

(3)自动化程度高。

基于此，本发明设计了一种InSAR数据升降轨联合滑坡分析方法，综合剖面线水平垂直形变结果，以及坡面的形变信息以及坡度坡向，对滑坡整体稳定性进行评价，以解决上述提到的问题。

发明内容

本发明的目的在于提供一种的InSAR数据升降轨联合滑坡分析方法，以解决上述背景技术中提出的滑坡面积大、范围广，利用水准测量和倾斜测量进行监测，需要耗费大量的人力物力；传统监测，一般测量的是当下绝对值，无法对未来做出预期的问题。

为实现上述目的，本发明提供如下技术方案：一种InSAR数据升降轨联合滑坡分析方法，包括以下步骤：

S1、挑选目标区域最合适的时间序列SAR影像，基于PSP算法生产出该区域InSAR升降轨数据，选出同一区域的升轨数据和降轨数据，根据区域范围生成经纬网格，通过升降轨数据插值出经纬网格的形变速率；

S2、下载该区域的STRM90的DEM数据，读取对应范围的DEM，运用卷积计算的方法求出栅格点对应的坡度和坡向值，通过插值的方法将经纬网格与升降轨数据的经纬网格统一；

S3、通过升降轨速率和坡度坡向，将升降轨形变速率分解为水平形变速率和垂直形变速率，得到滑坡的水平形变量和垂直形变量；

S4、挑选滑坡上局部形变最大值作为初始化隐患点，获取沿坡面的剖面线，并对其中的三条剖面线1、2、3进行分析，得到沿剖面线的坡面水平和垂直方向的滑动情况。

优选的，运用卷积计算的方法的具体步骤如下：

S2.1、通过自定义栅格单元计算窗口；

S2.2、根据计算窗口选择X方向的卷积模和Y方向的卷积模；

S2.3、对所述DEM进行傅里叶变换，之后再频率域里乘以得到的X方向的卷积模和Y方向的卷积模，再经过逆傅里叶变化生成X方向和Y方向的坡度；

S2.4、由X方向和Y方向的坡度计算坡度和坡向。

优选的，计算坡度和坡向的公式为：

Aspect＝Slope_sn/Slope_we (4)

Slope为坡度，aspect为坡向，Slope_we为X方向的坡度，Slope_sn为Y方向的坡度。

优选的，所述步骤S3中降轨形变速率分解公式为：

d_D＝d_V cos(θ_DV)+d_H cos(θ_DH) (7)

其中，d_D：降轨方向形变量；d_H：水平方向形变量；d_V：垂直方向形变量；θ_DV表示降轨视线方向和垂直方向夹角；θ_DV表示降轨视线方向和水平方向夹角。优选的，所述步骤S3升降轨形变速率分解公式为：

其中，d_A：升轨方向形变量；θ_AV表示升轨视线方向和垂直方向夹角；θ_AV表示升轨视线方向和水平方向夹角。

优选的，所述步骤S4中获取沿坡面的剖面线具体步骤如下：

S4.1、对于坡上某一点P0，可以求得该点对应的坡向a0，坡度b0，海拔 h0，并且初始化步长s0；

S4.2、按照步长s到达Pi+1获得i+1点的坡向ai+1，坡度bi+1，海拔hi+1；

S4.3、如果满足hi+1>hi，则进行判断ai+1-ai的绝对值是否小于某个阈值，如果是,则s＝s*0.9直到s<0.1s0终止循环，如果否，则s＝s*1.1，重新执行步骤 S4.2；

S4.4、如果不满足hi+1>hi，则判断bi+1是否连续小于阈值，如果是，则终止循环，否则i＝i+1，重新进入步骤S4.2；

S4.5、得到对于点P0沿坡向下的剖面线。

与现有技术相比，本发明的有益效果是：通过InSAR技术，获取滑坡升轨和降轨的PS点，PS点能够反映出该区域沿视线方向的形变信息，结合DEM 提取坡度坡向信息和将视线方向的形变信息转化为水平和垂直的形变信息，再根据重点区域自动提取沿滑坡方向的剖面线，根据剖面线对沿线进行水平和垂直形变分析，综合剖面线水平垂直形变结果，以及坡面的形变信息以及坡度坡向，对滑坡整体稳定性进行评价。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明分析流程图；

图2为本发明DEM求取坡度坡向流程图；

图3为本发明升降轨形变速率分解图；

图4为本发明沿坡面获取剖面线流程图；

图5为本发明PSP-InSAR降轨数据结果展示图；

图6为本发明PSP-InSAR升轨数据结果展示图；

图7为本发明初始隐患点分布展示图；

图8为本发明根据隐患点选择主滑面展示图；

图9为本发明主滑面1剖面分析图；

图10为本发明主滑面2剖面分析图；

图11为本发明主滑面3剖面分析图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

请参阅图1-11，本发明提供一种技术方案：一种InSAR数据升降轨联合滑坡分析方法，以金沙江流域某山区为例，包括以下步骤：

S1、挑选目标区域最合适的时间序列SAR影像，基于PSP算法生产出该区域InSAR升降轨数据，选出同一区域的升轨数据和降轨数据，根据区域范围生成经纬网格，通过升降轨数据插值出经纬网格的形变速率，如图5、6所示；

S2、下载该区域的STRM90的DEM数据，读取对应范围的DEM，如图 6所示，运用卷积计算的方法求出栅格点对应的坡度和坡向值，通过插值的方法将经纬网格与升降轨数据的经纬网格统一；

运用卷积计算的方法的具体步骤如下：

S2.1、通过自定义栅格单元计算窗口；

S2.2、根据计算窗口选择X方向的卷积模和Y方向的卷积模；

S2.3、对所述DEM进行傅里叶变换，之后再频率域里乘以得到的X方向的卷积模和Y方向的卷积模，再经过逆傅里叶变化生成X方向和Y方向的坡度；

S2.4、由X方向和Y方向的坡度计算坡度和坡向。

计算坡度和坡向的公式为：

Aspect＝Slope_sn/Slope_we (4)

Slope为坡度，aspect为坡向，Slope_we为X方向的坡度，Slope_sn为Y方向的坡度。

S3、通过升降轨速率和坡度坡向，将升降轨形变速率分解为水平形变速率和垂直形变速率，得到滑坡的水平形变量和垂直形变量，如图3所示；

降轨形变速率分解公式为：

d_D＝d_V cos(θ_DV)+d_H cos(θ_DH) (7)

其中，d_D：降轨方向形变量；d_H：水平方向形变量；d_V：垂直方向形变量；θ_DV表示降轨视线方向和垂直方向夹角；θ_DV表示降轨视线方向和水平方向夹角。

对于升轨同理，联立公式进行升降轨结合得到：

其中，d_A：升轨方向形变量；θ_AV表示升轨视线方向和垂直方向夹角；θ_AV表示升轨视线方向和水平方向夹角。

S4、挑选滑坡上局部形变最大值作为初始化隐患点，选用的初始化方法为选取局部形变最大值点，如图7所示，获取沿坡面的剖面线，初始化步长为100m，坡向之差的阈值为45°，

获取沿坡面的剖面线具体步骤如下：

S4.1、对于坡上某一点P0，可以求得该点对应的坡向a0，坡度b0，海拔 h0，并且初始化步长s0；

S4.2、按照步长s到达Pi+1获得i+1点的坡向ai+1，坡度bi+1，海拔hi+1；

S4.3、如果满足hi+1>hi，则进行判断ai+1-ai的绝对值是否小于某个阈值，如果是,则s＝s*0.9直到s<0.1s0终止循环，如果否，则s＝s*1.1，重新执行步骤 S4.2；

S4.4、如果不满足hi+1>hi，则判断bi+1是否连续小于阈值，如果是，则终止循环，否则i＝i+1，重新进入步骤S4.2；

S4.5、得到对于点P0沿坡向下的剖面线1、2、3。

得到剖面线如图8所示，并对其中的三条剖面线1、2、3进行分析，得到沿剖面线的坡面水平和垂直方向的滑动情况，如图9、10和11所示。

在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。

Claims (6)

1.一种InSAR数据升降轨联合滑坡分析方法，其特征在于，包括以下步骤：

S1、挑选目标区域最合适的时间序列SAR影像，基于PSP算法生产出该区域InSAR升降轨数据，选出同一区域的升轨数据和降轨数据，根据区域范围生成经纬网格，通过升降轨数据插值出经纬网格的形变速率；

S2、下载该区域的STRM90的DEM数据，读取对应范围的DEM，运用卷积计算的方法求出栅格点对应的坡度和坡向值，通过插值的方法将经纬网格与升降轨数据的经纬网格统一；

S3、通过升降轨速率和坡度坡向，将升降轨形变速率分解为水平形变速率和垂直形变速率，得到滑坡的水平形变量和垂直形变量；

S4、挑选滑坡上局部形变最大值作为初始化隐患点，获取沿坡面的剖面线，并对其中的三条剖面线1、2、3进行分析，得到沿剖面线的坡面水平和垂直方向的滑动情况。

2.根据权利要求1所述的一种InSAR数据升降轨联合滑坡分析方法，其特征在于：运用卷积计算的方法的具体步骤如下：

S2.1、通过自定义栅格单元计算窗口；

S2.2、根据计算窗口选择X方向的卷积模和Y方向的卷积模；

S2.3、对所述DEM进行傅里叶变换，之后再频率域里乘以得到的X方向的卷积模和Y方向的卷积模，再经过逆傅里叶变化生成X方向和Y方向的坡度；

S2.4、由X方向和Y方向的坡度计算坡度和坡向。

3.根据权利要求2所述的一种InSAR数据升降轨联合滑坡分析方法，其特征在于：计算坡度和坡向的公式为：

Aspect＝Slope_sn/Slope_we (4)

Slope为坡度，aspect为坡向，Slope_we为X方向的坡度，Slope_sn为Y方向的坡度。

4.根据权利要求1所述的一种InSAR数据升降轨联合滑坡分析方法，其特征在于：所述步骤S3中降轨形变速率分解公式为：

d_D＝d_Vcos(θ_DV)+d_Hcos(θ_DH) (7)

其中，d_D：降轨方向形变量；d_H：水平方向形变量；d_V：垂直方向形变量；θ_DV表示降轨视线方向和垂直方向夹角；θ_DV表示降轨视线方向和水平方向夹角。

5.根据权利要求4所述的一种InSAR数据升降轨联合滑坡分析方法，其特征在于：所述步骤S3升降轨形变速率分解公式为：

其中，d_A：升轨方向形变量；θ_AV表示升轨视线方向和垂直方向夹角；θ_AV表示升轨视线方向和水平方向夹角。

6.根据权利要求1所述的一种InSAR数据升降轨联合滑坡分析方法，其特征在于：所述步骤S4中获取沿坡面的剖面线具体步骤如下：

S4.1、对于坡上某一点P0，可以求得该点对应的坡向a0，坡度b0，海拔h0，并且初始化步长s0；

S4.2、按照步长s到达Pi+1获得i+1点的坡向ai+1，坡度bi+1，海拔hi+1；

S4.3、如果满足hi+1>hi，则进行判断ai+1-ai的绝对值是否小于某个阈值，如果是,则s＝s*0.9直到s<0.1s0终止循环，如果否，则s＝s*1.1，重新执行步骤S4.2；

S4.4、如果不满足hi+1>hi，则判断bi+1是否连续小于阈值，如果是，则终止循环，否则i＝i+1，重新进入步骤S4.2；

S4.5、得到对于点P0沿坡向下的剖面线1、2、3。