一种融合多源SAR数据的潜在滑坡体积估算方法

一种融合多源sar数据的潜在滑坡体积估算方法
技术领域
1.本发明涉及地质灾害防治安全技术领域，尤其涉及一种融合多源sar数据的潜在滑坡体积估算方法。


背景技术：

2.近年来，在滑坡监测领域，“空-天-地”的监测手段不断发展，灾害的监测及预测得到了较大的提升。合成孔径雷达干涉(synthetic aperture radar interferometry,insar)技术作为一种主动遥感技术，因其具有监测范围广，不受天气限制，可以全天时、全天候获取影像的优势，在滑坡灾害的监测中得到了较好的应用。
3.在通常的灾害监测方面，目前较多的研究仍是基于地表变形进行监测的，无法具体地揭示滑坡的危险程度，且对于滑坡的危害及滑坡体积的研究，常用的单点接触式探测，其成本较高、耗时长并且受数据量限制，采用少量点进行大范围插值难免会影响滑坡面位置和形状的估计精度，不利于对滑坡体积、形变模式的精确估计，影响滑坡的危险性评价。
4.通过使用insar技术获取地表的面状形变监测结果，结合弹性位错模型寻找最佳滑动面，可以实现滑坡体积的估计，为灾害的危险性预测提供数据支撑。而现有的对潜在滑坡体积估计的研究方法多集中于钻孔及探地雷达等接触式的探测，其不仅成本高、耗时长，且不利于大型滑坡与深层滑坡的滑动面及深度探测。


技术实现要素：

5.本发明提供了一种融合多源sar数据的潜在滑坡体积估算方法，目的是为了解决现有技术中存在的缺点。
6.为了实现上述目的，本发明提供如下技术方案：一种融合多源sar数据的潜在滑坡体积估算方法，包括如下步骤：
7.获取潜在滑坡区域的sar数据；
8.采用stacking insar技术对sar数据进行时序insar处理，获取升、降轨视线向形变速率图；
9.根据升、降轨视线向形变速率图、数字地形高程模型dem及光学影像确定滑坡表面变形的边界；
10.确定okada弹性位错模型反演滑坡滑动面时所需的几何参数及运动学参数；
11.将升、降轨视线向形变速率图转换成离散点，并添加对应的经纬度坐标，作为反演的表层约束条件；
12.输入表层约束条件于okada弹性位错模型对滑动面进行反演，并选择不同的平滑因子，进行多次反演计算；
13.通过反演得到多个模拟结果，将模拟结果与观测结果进行对比，选取最佳拟合度的模拟结果，从而获取滑动体距离地表的深度及滑动面的形变，确定出滑动面上的形变范围及滑动面形变方向；
14.通过滑动面形变面积与滑动体深度的积分计算滑坡体积，并根据滑动面上的滑动面形变方向确定滑坡的形变类型。
15.优选的，对所述sar数据进行形变监测，所述形变监测的过程包括：
16.数据获取：下载sentinel-1的升、降轨数据、30米的srtm dem数据及sentinel-1精密轨道数据，并通过对上述数据进行预处理后获取单视复数影像slc；
17.形变监测：采用stacking insar技术获取滑坡区域的形变速率图；
18.滑坡变形区域确定：根据获取的形变监测结果，对比升、降轨影像监测结果，结合dem地形图及光学影像，确定滑坡形变边界，并计算出该边界范围的面积。
19.优选的，获取所述形变速率图的具体过程包括：
20.通过设置一定空间和时间基线阈值，在满足条件的sar影像间进行干涉处理，生成多个干涉图集；
21.每个干涉图的相位主要包括dem误差造成的相位形变产生的相位大气延迟所产生的相位噪声所引起的相位
[0022][0023]
通过剔除相位中的获取形变产生的相位信息；
[0024]
采用stacking insar技术获取雷达视线向的年平均形变速率图。
[0025]
优选的，所述stacking insar技术在估计形变速率时，通过对若干个干涉解缠图进行相位叠加，获取形变速率图；
[0026]
其中，权重据主从影像的时间间隔计算，stacking insar的数学模型如下：
[0027][0028]
式中，ph_rate为相位形变速率，δti为第i个解缠图的两幅影像之间的时间间隔，为第i个干涉图的解缠相位值。
[0029]
优选的，反演所述滑动面时所需要的几何参数及运动学参数包括宽、深度、走向角、倾角、最大滑动量、滑动角，其中，参数确定的具体方法如下：
[0030]
初始深度设置：滑坡体的厚度，根据平均厚度进行迭代优化；
[0031]
初始宽度设置：沿主滑方向的的形变区域的宽度；
[0032]
走向角的设置：先确定滑坡的主滑方向，然后做主滑方向的垂线，垂线与北方向的夹角即为初始走向角；
[0033]
初始倾角设置：设置与边坡的坡度接近的区间进行搜索，通过不断调整确定最佳的倾角；
[0034]
最大滑动量：设置滑动面的最大滑动限制值；
[0035]
初始滑动角：设置一定的搜索范围，不断调整区间的最大最小的滑动角值，通过多次反演计算确定滑动角的大小。
[0036]
优选的，采用okada弹性位错模型对滑动面进行反演的过程包括如下步骤：
[0037]
结合insar技术获得的形变监测结果，确定滑坡的形变边界和区域；
[0038]
使用四叉树采样法对数据进行降采样；
[0039]
设定要反演滑动面区域的规模和方向，构建位错反演模型，并将滑动面沿着走向、倾向划分为多个子位错源，给定滑动面的长、宽、深度及子位错的分块情况。
[0040]
优选的，通过okada弹性位错模型计算每一个子块位错的形变，具体计算公式为：
[0041]
s＝[g
t
c-1
g]-1
gc-1
x
[0042]
其中，s为每个子位错源的滑动量，g为由地面形变与滑动面之间的联系建立的格林函数，c为观测值数据的协方差，x为地表的形变。
[0043]
优选的，所述滑坡体积的具体计算公式如下：
[0044][0045]
其中，a，b分别用a1，a2，b1，b2，h来表示，a1，b1分别为滑动面的外接椭圆的长短半轴，a2，b2分别为滑坡表面的外接椭圆的长短半轴，h表示滑坡体深度，则公式的具体变形公式为：
[0046]
a＝a1+z(a2-a1)/h
[0047]
b＝b2+z(b2-b1)/h
[0048][0049]
优选的，所述滑坡表面变形的参数包括变形范围、形变量级。
[0050]
本发明与现有技术相比具有以下有益效果：
[0051]
1、利用时序insar技术进行多源数据形变探测，获取不同数据的形变监测结果，通过比较，确定潜在滑坡在地表的面积和周长，并结合弹性位错模型，寻找一个最可能的滑坡基底滑动的位错平面，并获取滑动面的位置、形状和大小，并结合滑坡深度，利用该平面计算潜在滑坡的体积，精确判定潜在滑坡体的方量，为后续的数值模拟等滑坡灾害预测方法提供数据支持，具有较好的应用价值。
[0052]
2、本发明与质量守恒反演滑坡体厚度的方法不同，不需要使用流变参数等外部的参数，不仅获取了潜在滑坡体厚度，并反演了滑动面的形变信息，在此基础上，结合了insar的形变信息计算了滑坡体的体积。
附图说明
[0053]
图1为本发明提供的一种融合多源sar数据的潜在滑坡体积估算方法的具体流程图；
[0054]
图2为本发明提供的滑坡位错模型及参数示意图；
[0055]
图3为本发明提供的滑坡体体积估算示意图；
[0056]
图4为本发明提供的升、降轨形变速率图；
[0057]
图5为本发明提供的粗糙度和拟合残差曲线图；
[0058]
图6为本发明提供的升、降轨反演结果图；
[0059]
图7为本发明提供的最优滑动分布图；
[0060]
图8为本发明提供的滑坡体1各滑坡单元形变方向图；
[0061]
图9为本发明提供的滑坡体1的坡向图；
[0062]
图10为本发明提供的滑坡体2各滑坡单元形变方向图；
[0063]
图11为本发明提供的滑坡体2的坡向图。
具体实施方式
[0064]
下面结合附图，对本发明的具体实施方式作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
[0065]
一种多源sar数据融合的潜在滑坡体积估计方法结合insar技术获取的地面表层形变及弹性位错模型所确定的滑动面进行滑坡体积的估算。如图1所示，一种多源sar数据融合的潜在滑坡体积估计方法的技术方案如下：
[0066]
步骤1：采用sentinel-1a数据使用stacking insar技术进行sar数据处理，获取升、降轨视线向形变速率图。
[0067]
步骤2：根据形变速率图结果、数字地形高程模型dem及光学影像，确定滑坡表面变形的参数，包括变形范围、形变量级等。
[0068]
步骤3：确定okada弹性位错模型反演滑坡的滑动面时所需要的几何参数及运动学参数：宽、深度、走向角、倾角、最大滑动量、滑动角等参数；其中，参数确定的具体方法如下：
[0069]
(1)初始深度设置：滑坡体的厚度，根据平均厚度进行迭代优化。
[0070]
(2)初始宽度设置：沿主滑方向的的形变区域的宽度。
[0071]
(3)走向角的设置：先确定滑坡的主滑方向，然后做主滑方向的垂线，垂线与北方向的夹角即为初始走向角。
[0072]
(4)初始倾角设置：设置与边坡的坡度接近的区间进行搜索，通过不断调整确定最佳的倾角。
[0073]
(5)最大滑动量：设置滑动面的最大滑动限制值。
[0074]
(6)初始滑动角：设置一定的搜索范围，不断调整区间的最大最小的滑动角值，通过多次反演计算确定滑动角的大小。
[0075]
步骤4：将升降轨的形变速率图转成离散点，并添加对应的经纬度坐标，作为反演的表层约束条件。
[0076]
步骤5：对升、降轨视线向形变速率图进行降采样。通常采用四叉树与均匀降采样的方法，可以达到减少计算量，提高计算速度的目的。
[0077]
步骤6：输入表层约束条件于okada弹性位错模型进行滑动面的反演。选择不同的平滑因子，进行多次反演计算，根据拟合结果的l曲线选择了最优平滑因子，减少不确定性，提升反演结果的稳定性。通过对比反演得到的多个模拟结果与观测结果，选取最佳拟合度的模拟结果，从而获取滑动体的深度及滑动面的形变，确定出滑动面上的形变范围及滑动面形变方向。
[0078]
步骤7：通过积分计算滑坡体积，并根据滑动面上的滑动面形变方向，确定滑坡的形变类型。
[0079]
本发明主要分为三个关键的步骤：sar形变监测，滑动面确定以及滑坡体积估计。
[0080]
一、sar形变监测包括：
[0081]
1、sar数据获取：下载sentinel-1的升、降轨数据，通过预处理获取单视复数影像slc。
[0082]
2、形变监测：采用stacking insar技术获取滑坡区域的形变速率图。首先，通过设置一定空间和时间阈值，在满足条件的sar影像间进行干涉，生成多个干涉图集。对于每一个干涉图，相位主要包括dem误差造成的相位形变产生的相位大气延迟所产生的相位噪声所引起的相位
[0083][0084]
然后，通过剔除相位中的获取形变产生的相位信息。在数据处理过程中，sar坐标系下的地形相位可以通过外部的dem数据模拟获取，通过差分剔除模拟获得的地形相位从而获取差分干涉图。对于地形误差，通过构建高程误差与轨道基线之间的关系进行去除；于大气相位误差，通过构建dem与大气相位之间的二次多项式进行模拟并加以去除，获取形变相位图。最后，采用stacking insar技术获取雷达视线向的年平均形变速率图，stacking insar技术在估计形变速率时，权根据主从影像的时间间隔计算，stacking insar的数学模型如下：
[0085][0086]
式中，ph_rate为相位形变速率，δti为第i个解缠图的两幅影像之间的时间间隔，为第i个干涉图的解缠相位值。在stacking insar处理过程中，为了减少大气扰动产生的误差，要使用去除大气误差的解缠干涉图，并剔除误差较大的干涉图，在满足相干性条件下，选择时间基线较长的干涉图进行处理，从而获取较为准确的形变速率。
[0087]
3、滑坡变形区域确定：根据获取的形变监测结果，对比升、降轨影像监测结果，结合dem地形图及光学影像，确定滑坡形变边界，并计算出该边界范围的面积。
[0088]
二、如图2所示，基于okada弹性位错模型获取滑动面位置及形态的确定，其具体包括：
[0089]
根据遥感手段及gnss等监测方式只能获取地表的信息，对于地球内部的变化很难通过这些技术进行判定。但对于滑坡来说，滑动面对滑坡活动起着重要的作用。根据前人研究，一些大型滑坡的滑动面可以假想为构造断层。因此，为了获取滑坡的体积，我们可以通过采用okada位错模型，假设滑坡为各向同性介质，insar观测到的位移场是由各向同性弹性半空间内的平面位错模拟的，通过寻找观测值与模拟值的最佳拟合度对应的结果来确定滑动面的位置、形状和大小，进而为滑坡体积的估算提供数据。
[0090]
首先，结合insar技术获得的形变监测结果，确定滑坡的形变边界和区域，并使用四叉树采样法对数据进行降采样，提高计算效率。然后，设定滑动面的规模和方向，构建位错反演模型，并将滑动面沿着走向、倾向划分为多个子位错源，给定滑动面的长、宽、深度及子位错的分块情况，具体参数见图2。最后，通过okada弹性位错模型计算每一个子块位错的形变。具体计算公式为：
[0091]
s＝[g
t
c-1
g]-1
gc-1
x
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0092]
式中，s为每个子位错源的滑动量，g为由地面形变与滑动面之间的联系建立的格
林函数，c为观测值数据的协方差，x为地表的形变。在反演过程中，由于子位错源的滑动量有时存在跳变现象，需要进行一定程度的平滑，来确保反演结果的稳定性。
[0093]
三、如图3所示，对于滑坡体积估算的具体过程为：
[0094]
结合insar获取的表面形变监测结果，可以确定滑坡体的上表面面积，根据弹性位错模型反演的滑动面可以确定出滑坡体的下表面的面积，并且获取了滑坡体的平均深度h，一般情况下，假设滑坡的上表面与滑动面接近平行，此时，将滑坡体的上下表面选取一个最为接近的椭圆，其可以表示为图3。
[0095]
设滑动面的外接椭圆的长短半轴分别为：a1，b1；滑坡表面的外接椭圆的长短半轴分别为：a2，b2，滑坡体深度为：h。取积分单元dz。距离顶面为z高。通过积分可以计算出滑坡体的体积，具体计算公式如下：
[0096][0097]
其中a,b分别用a1,a2,b1,b2,h来表示：
[0098]
a＝a1+z(a2-a1)/h
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0099]
b＝b2+z(b2-b1)/h
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0100]
将a,b代入公式，可以求出滑坡体体积v:
[0101][0102]
实施例1：
[0103]
本发明以瓦堆村潜在滑坡为例进行实验，对潜在滑坡体的体积进行计算，并依据滑动面反演结果分析滑坡滑动类型。
[0104]
1、insar获取的形变速率图：根据形变速率结果，确定反演的滑动面的几何大小与位置。对于包含不同的形变方向滑坡体时，要选择多个滑动面进行反演。其形变速率图如图4所示，其中图4a为升轨监测结果图，图4b为降轨监测结果，白色框为滑动面的反演区域。
[0105]
2、okada模型的粗糙度和拟合残差曲线:如图5所示，在滑动分布反演过程中，通过权衡模型粗糙度和拟合残差之间的曲线来评定滑动分布结果，实验通过多次模拟反演，选定平滑因子为0.03时，模型的反演结果最优。
[0106]
3、okada模型的滑动面反演结果如表1所示：
[0107]
表1弹性位错模型反演结果
[0108]
滑动面长度宽度深度(km)走向角倾角平均滑动量11.60.70.06920350.2420.910.06025230.21
[0109]
滑坡形变观测结果与滑动分布的对比结果如图6所示。其中，图6中a、b、c分别为升轨的观测值、模拟值、残差值，e、f、g分别为降轨的观测值、模拟值、残差值。
[0110]
滑动面的最优滑动分布结果如图7所示。
[0111]
4、滑坡体积估算
[0112]
根据形变速率图，圈定滑坡表面的最大外接椭圆，并结合反演的滑动面，确定滑动面的外接椭圆，从而获取滑坡体上下表面的长半轴与短半轴的长度，并结合反演获得的深度值，计算出该区域存在的两个滑坡体的体积分别为：
[0113]
(1)滑坡体1
[0114]
图8与图9分别为滑坡体1各滑坡单元形变方向图与滑坡体1坡向图。根据图4、图8与图9内容计算得出，a1＝0.35km，b1＝0.19km，a2＝0.6km，b2＝0.26km，h＝0.069km。根据体积计算公式得出，滑坡体1的体积v＝1.69*107m3。结合走向与倾向的形变，可以大致判断出滑坡的形变方向，对于处于坡脚的滑坡体1，结合坡向联合分析，滑坡主要沿着坡向运动，走向与倾向均存在变形。滑坡的运动一般可以分为平移、旋转或复杂运动方式。旋转滑坡一般向下和向外移动，平移滑体运动方向与山坡平行，结构上滑动面较为平整。根据反演滑动面的结果，滑坡体1为沿坡向的平移式滑坡。
[0115]
(2)滑坡体2
[0116]
图10与图11分别为滑坡体2各滑坡单元形变方向图与滑坡体2坡向图。根据图4、图10与图11内容计算得出，a1＝0.40km，b1＝0.24km，a2＝0.5km，b2＝0.31km，h＝0.06km，带入体积计算公式可以得出，滑坡体2的体积v＝2.03*106m3.根据反演的结果，滑坡体2为平移式滑坡。
[0117]
以上所述实施例仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换，均属于本发明的保护范围。