基于InSAR多维形变信息的黄土滑坡类型及滑动模式分析方法与流程

本发明涉及黄土滑坡类型识别及滑动模式分析领域，特别涉及一种基于insar多维形变信息的黄土滑坡类型及滑动模式分析方法。

背景技术：

黄土高原作为中华名族古代文明的发祥地之一，是中国最重要的能源、化工基地。近年来随着经济的发展及人口的扩张，使得其成为世界上水土流失最严重和生态环境最脆弱的地区之一。受人类活动及工程建设的影响，地面沉降、地裂缝、崩塌、滑坡、泥流等地质灾害频发，形成一条复杂的灾害链，严重危及交通干线、重大工程及人们生命财产的安全。在频繁发生的地质灾害中，黄土滑坡分布广泛且集中，突发性很强，破坏力最为严重。因此对黄土滑坡灾害的防灾减灾工作刻不容缓。其中对黄土滑坡的类型及滑动模式分析，对黄土滑坡灾害的工程治理及防灾减灾工作尤为关键。

合成孔径雷达干涉测量(interferometricsyntheticapertureradar,insar)作为一项新型的空间对地观测技术，具有全天候、全天时、空间分辨率高、基本不受气象条件的影响等的优点，近些年来被广泛的应用于滑坡的识别与监测中。对于黄土滑坡的类型及滑动模式进行分析，现有技术主要是通过地质工作人员的野外调查或者常规的工程勘察方法获得，这些方法需要丰富的地质工作经验，不仅耗时费力，覆盖面有限，还需要较大的经济、人力、物力的投入；现有技术中，仅利用单一的sar数据集获得滑坡体的一维los向形变速率及形变时间序列，无法揭示滑坡的空间形变特征，且由于时间分辨率低、难以捕获滑坡体瞬间突变信号，对滑坡类型及滑动模式的研究存在很大的局限性，难以深入分析滑坡体形变机理的缺点。

技术实现要素：

本发明的目的在于针对黄土滑坡类型及滑动模式分析中耗时费力，覆盖面有限，滑坡类型及滑动模式研究的局限性，单一sar数据无法揭示滑坡的空间形变特征，难以深入分析滑坡体形变机理的问题，提供一种基于insar多维形变信息的黄土滑坡类型及滑动模式分析方法。

为了实现上述目标，本发明采用以下技术方案：

一种基于insar多维形变信息的黄土滑坡类型及滑动模式分析方法，所述的方法包括以下步骤：

步骤1，通过覆盖研究区域的t个sar传感器采集的sar数据，获取滤波insar差分干涉图；

步骤2，对滤波insar差分干涉图进行相位解缠，将解缠后的滤波insar差分干涉图编码至地理坐标系，并采样至相同的坐标格网和空间分辨率，得到地理坐标系下的解缠干涉图；

步骤3，对地理坐标系下的解缠干涉图，计算研究区域的多维地表形变速率，并计算研究区域的多维累积地表形变；

针对于研究区域，按时间顺序选取多维累积地表形变构成多维地表形变时间序列；

步骤4，对获得的研究区域的多维地表形变速率及多维地表形变时间序列，结合研究区域的遥感影像、地形图进行滑坡形变机理分析，确定黄土滑坡的类型及滑动模式。

进一步地，所述的获取滤波insar差分干涉图，包括：

通过覆盖研究区域的t个sar传感器，获得覆盖研究区域的t个升降轨sar影像数据集并利用无人机摄影测量方式获得外部dem数据，对获取的t个升降轨sar影像数据集中每一个sar影像数据集的影像分别进行两两做差形成差分干涉图，然后将每张差分干涉图减去外部dem数据，得到insar差分干涉图，对insar差分干涉图进行滤波，得到滤波insar差分干涉图。

进一步地，所述的计算研究区域的多维地表形变速率，包括：

采用下式3计算研究区域的多维地表形变速率：

上式中，a¹到a^t分别表示1到t个升降轨影像数据集的系数矩阵，所述的多维地表变形速率包括vn、ve、vu，其中vn表示南北向的地表形变速率，ve表示东西向的地表形变速率，vu表示垂直方向的地表形变速率；到分别表示1到t个数据集的观测相位值。

进一步地，计算研究区域的多维累积地表形变，包括：

采用下式4计算研究区域的多维累积地表形变：

上式中，表示t个sar影像数据集总共的sar影像数量，分别为南北、东西以及垂直向的累积地表形变，表示第在第ti时刻地表的南北向形变速率；表示在第ti时刻地表的东西向形变速率；表示在第ti时刻地表的垂直向形变速率；δti为ti-1时刻与ti时刻的时间间隔。

本发明与现有技术相比具有以下技术效果：

1.本发明无需野外调查，对于作业人员的地质专业知识要求较低，对于作业人员难以到达的滑坡也能进行研究，适合大面积、高难度地区和危险地段的黄土滑坡类型及滑动模式的快速高效分析，与传统黄土滑坡类型及滑动模式分析时进行野外地质调查相比，没有局限性。

2.本发明仅需获取覆盖研究区域的升降轨sar影像数据集，经过数据处理获得研究区域的多维地表形变信息，结合研究区域的遥感影像、地形图进行黄土滑坡的类型及滑动模式分析，可以深入分析滑坡体形变机理；操作简单、自动化程度高、可靠性和效率高。

附图说明

图1为本发明提供的黄土滑坡类型及滑动模式分析流程图；

图2为本发明提供的利用insar技术获得的黄土-基岩接触面型滑坡一维、二维形变图；其中，(a)为升轨sar数据获得的滑坡los向形变速率，(b)为降轨sar数据获得的滑坡los向形变速率图，(c)为滑坡垂直向形变速率，(d)为滑坡东西向形变速率；

图3为本发明提供的试验区黄土-基岩接触面型滑坡的遥感影像及滑动模式；其中，(a)为试验区的遥感影像，(b)为基于传统野外地质调查方式获得的黄土-基岩接触面型滑坡的滑动过程，(c)为insar技术所获得的试验区二维形变时间序列；

图4为野外地质调查获取的试验区滑坡的现场照片；其中，(a)为滑坡的全貌，(b)为滑坡表面局部裂缝的放大图，(c)为滑坡后缘的裂缝；

图5为本发明提供的利用insar技术获得的浅层崩塌型滑坡一维、二维形变图；其中，(a)为升轨数据所获得的los向形变速率，(b)为降轨数据所获得的los向形变速率，(c)为垂直向形变速率，(d)为东西向形变速率；

图6为本发明提供的试验区浅层崩塌型滑坡的遥感影像及滑动模式；其中，(a)为试验区的遥感影像，(b)为基于传统野外地质调查方式获得的浅层崩塌型滑坡的滑动过程，(c)为insar技术所获得的试验区滑坡的二维形变时间序列；

图7为野外地质调查获取的试验区滑坡的现场照片；

图8为本发明提供的利用insar技术所获得的黄土滑坡渐进后退式滑动模式一维、二维形变图；其中，(a)为升轨数据所获得的los向形变速率，(b)为降轨数据所获得的los向形变速率，(c)为垂直向形变速率，(d)为东西向形变速率；

图9为本发明提供的试验区黄土滑坡渐进后退式滑动的遥感影像；

图10为基于野外地质调查方式获得的黄土滑坡渐进后退式滑动过程；其中，(a)为黄土滑坡的局部失稳示意图，(b)为黄土滑坡的第一次全局失稳示意图，(c)为黄土滑坡的第二次全局失稳示意图，(d)为黄土滑坡的第三次全局失稳示意图；

图11为利用insar技术获得的试验区滑坡的二维形变时间序列，其中，(a)表示图9中p2点的二维形变时间序列，(b)表示图9中p3点的二维形变时间序列，(c)表示图9中p4点的二维形变时间序列，(d)表示图9中p5点的二维形变时间序列；

图12为野外地质调查获取的试验区滑坡的现场照片；其中，(a)为图中iii区的放大图，(b)为滑坡的现场照片，(c)为图中白色矩形所示裂缝的放大图，(d)为图中i区的放大图，(e)为图中白色矩形所示的地洞的放大图。

具体实施方式

如图1至图12所示，本发明公开了一种基于insar多维形变信息的黄土滑坡类型及滑动模式分析方法，详细步骤如下：

步骤1，通过覆盖研究区域的t个sar(合成孔径雷达，简称sar)传感器，获取覆盖研究区域的t个升降轨sar影像数据集并利用无人机摄影测量方式获得外部dem(数字高程模型，digitalelevationmodel，简称dem)数据，对获取的t个升降轨sar影像数据集中每一个sar影像数据集中的影像分别进行两两做差形成差分干涉图，然后将每张差分干涉图减去外部dem数据，得到insar(合成孔径雷达干涉测量，syntheticapertureradarinterferometry,简称insar)差分干涉图，对insar差分干涉图进行滤波，得到滤波insar差分干涉图；

本方案中，通过sar传感器来获取覆盖研究区域的升降轨sar影像数据集并利用无人机摄影测量方式获得外部dem数据，本方案中有t个sar传感器，每一个sar传感器获取的升轨或者降轨影像构成一个数据集，t个sar传感器获取的sar数据集共同构成sar传感器获取的sar数据集；采用insar技术对sar传感器获取的sar数据集进行处理，具体过程为，首先，对sar传感器获取的sar数据集中的每一种sar数据集的影像进行两两做差形成差分干涉图，按照相同的处理方法对其他sar数据集进行处理，形成的差分干涉图是一系列图像；其次，用形成的每张差分干涉图减去外部dem数据，得到insar差分干涉图；最后，对insar差分干涉图做滤波处理，得到滤波insar差分干涉图；对图像进行滤波处理是为了使图像的清晰度更好，可以提高信噪比。

步骤2，对滤波insar差分干涉图进行相位解缠，将解缠后的滤波insar差分干涉图编码至地理坐标系，并采样至相同的坐标格网和空间分辨率，得到地理坐标系下的解缠干涉图；

本方案中，对每张滤波insar差分干涉图进行相位解缠，相位解缠的目的是使每张滤波insar差分干涉图的相位由主值或相位差值恢复至真实值，将解缠后的滤波insar差分干涉图编码至地理坐标系，通常情况下地理坐标系选择wgs84坐标系，也可以根据实际情况选择合适的坐标系，例如西安80坐标系等。将影像编码至地理坐标系，可在地理坐标系的含位置的数据与影像之间建立联系，把影像分配给地理坐标的含相应位置的数据记录，图像经过地理编码后，便可在空间中显示各影像的位置，这样方便可对信息执行进一步的分析；由于不同的sar传感器具有不同的空间分辨率，识别的地面目标大小不一致，因此难以将不同分辨率sar传感器形成的差分干涉图进行一起处理。此外，不同sar传感器获得的sar数据还存在地理基准不一致问题，即使编码到同一地理坐标系下，相同地物目标位置之间存在偏差。因此为了消除不同sar传感器获得干涉图空间分辨率不一致及存在地理偏差问题，需要将地理编码后的不同sar传感器数据集获得的解缠干涉图进行重采样，使其具有相同的坐标格网和空间分辨率。

步骤3，对地理坐标系下的解缠干涉图，计算研究区域的多维地表形变速率，并计算研究区域的多维累积地表形变；

针对于研究区域，按时间顺序选取多维累积地表形变构成多维地表形变时间序列；

具体地，本实施例中，计算研究区域的多维形变速率的方法为：

对于来自单一sar平台的数据集，仅有一个入射角θ和一个卫星飞行方位角α，根据小基线集(sbas)insar技术原理，一维的地表形变时间序列可以按照式1计算：

其中，a是一个m×n的系数矩阵，m为滤波insar差分干涉图的数量，n+1为t个sar传感器的数据集中的升降轨影像的数量，vlos为待求的卫星视线(los)向形变速率，为观测的相位值，a⁺是矩阵a的伪逆，为在tⁱ时刻los向的累积地表形变，为在tⁱ⁺¹时刻los向的累积地表形变，为在tⁱ⁺¹时刻los向的地表形变速率，δtⁱ⁺¹为tⁱ时刻与tⁱ⁺¹时刻的时间间隔。

对于t个来自不同sar传感器的升降轨影像数据集，则它们具有不同的入射角和卫星飞行方位角，假定：

vlos＝sv＝snvn+seve+suvu，

s＝{sn,se,su}＝{sinαsinθ,-cosαsinθ,cosθ}；

根据卫星los向形变与地面三维形变的投影关系，对于其中任一sar数据集t(t＝1,2,…t)，式1可以写成式2所示的形式：

其中，s是一个los向的单位向量，由南北、东西和垂直三个方向的分量sn、se、su构成；v是一个待求的地表形变速率向量，包含了北、东和垂直三个方向的形变分量vn、ve、vu，表示sar数据集t的相位观测值，表示sar数据集t的los向单位向量的南北分量；表示sar数据集t的los向单位向量的东西分量；表示sar数据集t的los向单位向量的垂直分量；对于t个来自不同sar传感器的升降轨影像数据集，式2可以写成式3所示的矩阵的形式：

或简写成

其中，

上式中，a¹到a^t分别表示1到t个升降轨影像数据集的系数矩阵，vn表示南北向的地表形变速率，ve表示东西向的地表形变速率，vu表示垂直方向的地表形变速率；到分别表示1到t个数据集的观测的相位值；表示方程的系数矩阵，表示待求的未知参数向量，即地表的三维形变速率；表示观测的相位值。

由于公式(3)中待求的未知参数的个数大于线性方程的个数，使得其系数矩阵秩亏，可采用奇异值分解(svd)或tikhonov正则化的方法获得方程的解，即三维地表形变速率。进一步采用式4可获得多维累积地表形变：

式中，表示t个sar影像数据集总共的sar影像数量，分别为南北、东西以及垂直向的累积地表形变，表示第在第ti时刻地表的南北向形变速率；表示在第ti时刻地表的东西向形变速率；表示在第ti时刻地表的垂直向形变速率；δti为ti-1时刻与ti时刻的时间间隔。

本实施例中，针对于研究区域，在一年当中，按时间顺序选取多个多维累积地表形变构成多维地表形变时间序列进行研究；以横坐标表示日期，纵坐标表示多维累积地表形变，如图3的(c)所示。

由于sar卫星对南北向的形变敏感性较低，如果仅当获得两个不同sar传感器的升降轨数据集时，式3、4中可以忽略南北向的形变，从而计算滑坡水平东西向与垂直向的二维地表形变。

步骤4，对获得的研究区域的多维地表形变速率及多维地表形变时间序列，结合研究区域的遥感影像、地形图进行滑坡形变机理分析，确定黄土滑坡的类型及滑动模式。

本方案中，滑坡的多维形变速率及时间序列，多维是指二维或者三维，本方案中实施例提供的是一维及二维形变图，仅当获得两个不同sar传感器的升降轨数据集时，可以忽略南北向的形变，从而计算滑坡水平东西向与垂直向的二维地表形变；根据获得的研究区域的二维形变速率，垂直向形变速率和东西向形变速率分析，判断在滑坡边界内，哪个方向的形变速率比较大，在滑坡边界内的覆盖范围广就可以得出滑坡以哪个方向的形变为主，例如，从图中得出滑坡在垂直方向上的滑坡边界中的形变速率比在东西向的形变速率大且覆盖范围广，就得出滑坡在垂直方向上存在形变；对于形变时间序列，将选出的日期作为横坐标，将根据对应的日期得到的垂直和东西方向的累积地表形变作为纵坐标，来构成坐标系，可以得到形变时间序列，本方案横坐标选用的日期为2016年1月到2016年12月，根据坐标系判断哪个方向的累积地表形变数据比较大，就是以哪个方向的形变为主，再以现有的滑坡类型作为对比，判断和黄土滑坡类型中的哪个滑坡类型的形变一样就属于哪个滑坡类型。为了验证判断的滑坡类型和滑动模式的可靠性，采用遥感影像和dem获得的地形图来进一步验证所判断的滑坡类型，可以和通过研究区域的二维形变速率及形变时间序列所判断的滑坡类型进行比对，确定分析的类型的正确性。

本发明的实验数据采用了真实的升降轨terrasar-x数据，升轨数据23景，降轨数据19景，入射角分别为41.2°和41.8°，像元分辨率距离向0.91米，方位向1.26米，总共覆盖的像元有5600×4200。覆盖地区为甘肃省黑方台黄土滑坡地区。

本发明首先对实验区域的升降轨terrasar-x数据进行单独处理，分别获得黑方台新塬滑坡群升降数据的形变速率图，进行实验区潜在滑坡的识别，共识别到三个潜在的黄土滑坡且升降轨数据所得到的结果一致。然后对这三个潜在的黄土滑坡的类型及滑动模式采用insar多维形变信息进行分析。高质量的升降轨解缠干涉图地理编码至wgs84坐标系下，并采样至3米空间分辨率的经纬度格网，采用公式3与公式4计算获得研究区域的二维形变速率及时间序列，结合试验区的googleearth影像及地形图进行滑坡类型及滑动模式的分析。

图2显示的是利用insar技术获得的黄土-基岩接触面型滑坡一维、二维形变图；(a)为升轨数据所获得的los向形变速率；(b)为降轨数据所获得的los向形变速率；(c)为垂直向形变速率；(d)为东西向形变速率。

图3显示的是试验区黄土-基岩接触面型滑坡的遥感影像及滑动模式；(a)为试验区的遥感影像；(b)为基于传统野外地质调查方式获得的黄土-基岩接触面型滑坡的滑动过程；(c)为insar技术获得的试验区滑坡的二维形变时间序列。

通过图2的(c)与(d)可以看出，该滑坡主要以东西向形变为主，且向西滑动，在垂直向上存在较小量级的形变。此形变特征与黄土-基岩接触面型黄土滑坡高度一致，由此我们可以判定此滑坡属于典型的黄土基岩接触面型滑坡，其滑动过程如图3中的(b)所示。在地下水的作用下，破裂面在黄土层中发育，并且向下层的基岩层传播。随着时间的发展，导致滑坡失稳，沿着基岩面滑动。为了验证insar多维形变信息所获得滑坡类型及滑动模式的可靠性，采用野外地质调查的方式来进行验证，此滑坡的现场照片如图4所示。从现场图可以看到，此滑坡的现场真实形变特征与insar技术所获得的二维形变特征高度一致，属于典型的黄土-基岩接触面型滑坡，证明了insar多维形变信息所分析的黄土滑坡类型的可靠性。

图5显示的是利用insar多维形变信息获得的浅层崩塌型黄土滑坡的一维、二维形变图；(a)为升轨数据所获得的los向形变速率；(b)为降轨数据所获得的los向形变速率；(c)为垂直向形变速率；(d)为东西向形变速率。

图6显示的是浅层崩塌型黄土滑坡的遥感影像及滑动模式；(a)为试验区滑坡的遥感影像；(b)为基于传统野外地质调查方式获得的浅层崩塌型黄土滑坡的滑动过程；(c)为试验区滑坡的二维形变时间序列。

通过图5的(c)与(d)可以看出，该滑坡主要以垂直向形变为主，在东西向上存在较小量级的形变。此形变特征与浅层崩塌型黄土滑坡高度一致，由此可以判定此滑坡属于典型的浅层崩塌型黄土滑坡，其滑动过程如图6中的(b)所示。破裂面在整个黄土层内发育，随着时间的发展，破裂面失稳发生滑动，该类型的滑坡通常具有较小的尺寸。为了验证insar技术所获得滑坡类型及滑动模式的可靠性，采用野外地质调查的方式来进行验证，此滑坡的现场照片如图7所示。从现场图可以看到，此滑坡的现场真实形变特征与insar技术所获得的二维形变特征高度一致，属于典型浅层崩塌型黄土滑坡，证明了insar多维形变信息所分析的黄土滑坡类型的可靠性。

图8显示的是利用insar技术获得的黄土滑坡渐进后退式滑动模式一维、二维形变图；(a)为升轨数据所获得的los向形变速率；(b)为降轨数据所获得的los向形变速率；(c)为垂直向形变速率；(d)为东西向形变速率。

图9显示的是试验区滑坡的遥感影像；图10显示的是基于传统野外地质调查方式获得的渐进后退式黄土滑坡的滑动过程；图11是insar技术所获得滑坡的二维形变时间序列。

从图8的(c)与(d)中可以看出，该滑坡整体向东滑动，垂直形变仅出现在滑坡的边缘。此形变特征与黄土滑坡的渐进后退式滑动模式高度一致，由此可以判定此滑坡属于典型的渐进后退式滑动模式，其滑动过程如图10所示。在地下水作用下，饱水黄土产生静态液化，使上部黄土以底部软弱基座为底滑面，顶部拉张裂缝为后侧边界，产生局部失稳，并激发超孔隙水压力引起黄土解体形成流态化堆积。在局部失稳发生后，台塬边形成弧形凹槽，形成新的局部临空面。滑坡后缘的应力继续发生改变，并产生新的近直立的拉张裂缝，产生第一次全局失稳，如此往复继续产生第二次、第三次全局失稳，在同一部位形成渐进后退式滑坡，故垂直向形变仅出现在滑坡体的边缘。为了验证insar技术所获得滑坡类型及滑动模式的可靠性，同样采用野外地质调查的方式来验证，此滑坡的现场照片如图12所示。从现场图可以看到，此滑坡的现场真实形变特征与insar技术所获得的二维形变特征高度一致，属于典型的渐进后退式滑动模式，证明了insar多维形变信息所分析的黄土滑坡类型的可靠性。

以上公开的仅为本发明的几个具体实施例，但是，本发明实施例并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。