一种区域尺度上滑坡对季节性降雨响应的分析方法，收集目标区域相关资料，如周围地质地层分布情况，收集与调查气温、降水等自然条件资料，收集目标区域内的历史SAR影像，对历史滑坡灾害的危害程度、诱因、成灾过程进行判断，初步划定分析的目标区域存在的大体范围；利用Stacking‑InSAR、时序InSAR技术对目标区域进行形变识别与监测，获取区域内滑坡位移数据；结合InSAR结果和现场调查数据，对综合选取的典型滑坡进行验证来保证滑坡位移数据的正确性；最后，利用小波分析工具，探究出降雨与滑坡位移量之间的定量关系。
一种区域尺度上滑坡对降雨响应的分析方法
技术领域
[0001] 本发明属于地质灾害技术领域，具体涉及了一种区域尺度上滑坡对降雨响应的分析方法。
背景技术
[0002] 降雨是大多数滑坡的重要触发因素，由降雨造成的滑坡在世界范围内达到三分之二。近几十年来，全球气候变化使降雨变得更加异常，特别是在高海拔和地形起伏地区，这显著增加了滑坡发生的频率。而降雨与滑坡的关系也一直是研究的重点和难点。
[0003] 研究滑坡形变对降雨的响应可以帮助我们了解滑坡的运动过程，从而减少灾害的发生。InSAR（Interferometic SAR）作为一种大范围监测小尺度变形的微波遥感技术，在滑坡识别与监测中发挥着越来越重要的作用，其监测出来的位移数据为滑坡的进一步研究提供基础数据。此外也有一些专家借助于模型来量化降雨对滑坡的影响，但模型一般具有极大的限制性和不确定性；而且如果只关注单一的时间维度，可能会忽略某些重要信息。
发明内容
[0004] 为了克服上述现有技术的不足，本发明的目的是提供一种区域尺度上滑坡对降雨响应的分析方法，其主要内容有四点，其一，收集目标区域相关资料，如周围地质地层分布情况，收集与调查气温、降水等自然条件资料，收集目标区域内的历史SAR影像，对历史滑坡灾害的危害程度、诱因、成灾过程进行判断，初步划定分析的目标区域存在的大体范围；其二，利用Stacking‑InSAR、时序InSAR技术对目标区域进行形变识别与监测，获取区域内滑坡位移数据；其三，结合InSAR结果和现场调查数据，对综合选取的典型滑坡进行验证来保证滑坡位移数据的正确性；最后，利用小波分析工具，探究出降雨与滑坡位移量之间的定量关系。通过将成熟的InSAR技术与时频分析工具小波分析的完美结合，从新的角度在大区域尺度上定量研究降雨与滑坡位移关系的不确定性。
[0005] 为了实现实现上述目的，本发明采用的技术方案是：
[0006] 一种区域尺度上滑坡对季节性降雨响应的分析方法，其包括以下步骤：
[0007] 步骤1、收集目标区域相关资料，对历史滑坡灾害的危害程度、诱因、成灾过程进行判断，初步划定分析的目标区域存在的大体范围；
[0008] 步骤2、根据目标滑坡区域的位置特点，基于InSAR技术处理雷达卫星影像来获取滑坡位移数据，包括以下三个方面：
[0009] 步骤2.1、对研究区SAR影像数据进行前期处理：
[0010] 获取目标区域SAR数据集后，需要对卫星雷达数据进行配准、去斜、滤波等预处理，通过引入精密轨道数据（POD）可有效去除因轨道误差引起的系统性误差，通过引入外部DEM、时空滤波去除部分误差形成相位解缠的差分干涉图；
[0011] 步骤2.2、大区域范围内滑坡灾害的识别：
[0012] 对前期处理过的SAR卫星数据集借助于Stacking技术进行大区域内滑坡灾害的识别，Stacking技术在获取M幅相位解缠的差分干涉图后，将解缠相位加权叠加，估计出平均相位变化速率，从而最大程度地减少了大气误差，提高了形变速率地解算精度，通过大范围的Stacking识别聚焦灾害严重的区域作为研究区，再对形变严重区进行长时间监测分析来获取更精确的位移结果；
[0013] 步骤2.3、区域滑坡的时序分析：
[0014] 运用联合PS点和DS点的时序InSAR分析方法，进行相位分析来获取时间序列形变，根据目标区域监测的形变结果，均匀地选取几个典型大滑坡的位移数据，通过分析典型滑坡位移量对降雨响应的共性来定量说明区域上滑坡与降雨的关系；
[0015] 步骤3、基于InSAR识别监测结果，结合现有位移资料，进行区域内典型滑坡的现场调查，对区域内滑坡监测精度进行评价，InSAR结果中剧烈形变的位置，发生了严重的破坏，也对应剧烈位移变化；
[0016] 步骤4、利用小波分析工具，探究出降雨与滑坡位移数据间的定量关系：
[0017] 步骤4.1、数据重采样：
[0018] 由于基于SAR卫星获取的滑坡位移数据以重访周期12天为一单元间隔，而对应时间段内降雨数据以一日为单元，为了减少数据误差，将日降雨数据重采样为以12天为间隔的数据，此外针对滑坡位移数据中由于卫星影像缺失可能存在的个别数据缺失，将根据InSAR原理采用线性插值来补充；
[0019] 步骤4.2、非线性位移数据提取：
[0020] 结合滑坡成因，分析监测的滑坡位移的时间序列，认为其位移的线性趋势主要是重力作用下的非弹性变形，那时序位移曲线上的小波动则应该与外部诱发因素密切相关，而对各诱因进行主成分、相关性分析后，则确定位移周期项是对主诱因降雨的体现，因此，采用时间序列加法模型将滑坡位移时间序列分解为线性趋势和非线性趋势两部分，线性分量通过线性最小二乘拟合计算，非线性分量定义为位移时间序列与之前计算的线性分量的差值，其表现出明显的周期性，因此，降雨对滑坡的影响就集中体现在降雨与非线性位移数据的定量关系上；
[0021] 步骤4.3、连续小波变换CWT：
[0022] 由于非线性时间序列位移、降雨数据都表现出明显周期性，将这两个数据借助于CWT转换到时频域，根据其二维图像中信号的时频分布来获取周期性、季节性信息，这即是对数据的时频处理；
[0023] 步骤4.4、交叉小波变换XWT和小波相干WTC：
[0024] 对连续小波变换后的双时间序列进行交叉小波变换CWT和小波相干分析WTC，可以得到周期性位移量与降雨量之间的交叉小波能量谱和交叉小波凝结谱，XWT的计算方法是将一个时间序列的CWT乘以第二个时间序列的CWT的复共轭，结果为二维图像，因此XWT高值的位置表示了两个时间序列中高值共振能量，WTC定义为两个时间序列的CWT在时频域的相干性，它是通过两者之间的标准化互相关来计算的，因此在两个CWT时频图局部相似，即相干的时频区域内，XWT的绝对值越大、WTC值会很高，即接近于1；
[0025] 步骤5、CWT结果表明：多年间降雨数据和形变数据都显著存在约1年的周期。此外，在每年的降水湿季，具有突出的动力信号，这是降雨影响下滑坡运动季节性表现。对CWT后的双时间序列进行XWT和WTC处理。图中的箭头表示相对相位关系，右箭头表示相同相位变化（正相关），左箭头表示相反相位变化（负相关）。向上的箭头表示滑坡运动先于降雨1/4个周期，向下的箭头表示其滞后于降雨。通过典型滑坡的一致性规律得到区域尺度上滑坡对降雨相应的相关性和滞后性。
[0026] 对CWT后的双时间序列进行XWT和WTC处理。箭头表示相对相位关系，右箭头表示相同相位变化（正相关），左箭头表示相反相位变化（负相关）。向上的箭头表示滑坡运动先于降雨1/4个周期，向下的箭头表示其滞后于降雨。通过典型滑坡的一致性规律得到区域尺度上滑坡对降雨相应的相关性和滞后性。
[0027] 所述的目标区域相关资料包括周围地质地层分布情况，气温、降水等自然条件资料，历史SAR影像资料。
[0028] 所述的进行相位分析来获取时间序列形变，主要包括三个步骤：PS点的选取、DS点的选取、PS点和DS点的结合分析；
[0029] PS点是指相干性较高的稳定点，是通过振幅离差指数（ADI）选取的。DS点在短时间内保持一定稳定性的目标体，在空间上分布广泛，但是，DS占主导地位的低信噪比像素不能直接用于时间序列分析，而分布散射体在空间分布上统计一致，因此将通过3倍多视滤波后提高的像素点作为DS点。最后将选取的PS和DS点合并在一起，与PS点共享相同雷达坐标的DS点被丢弃，将其在商业软件GAMMA中用最小成本流(MCF)算法进行相位解缠。然后根据解缠后的相位分量各自的特点，逐一计算并去除轨道误差、高程误差、大气扰动和噪声等相位分量。最后获得正确可靠的位移结果。
[0030] 本发明的有益效果是：
[0031] 通过使用不同的InSAR技术进行滑坡识别与监测，有效获取了区域内滑坡的高精度位移数据，并通过实地调查验证了其准确性和可靠性。作为可靠的位移数据源，将其与小波分析工具完美结合，可发现区域内滑坡运动与降雨的关系。即年际尺度上和年内尺度上降雨对区域内滑坡运动的周期性、季节性影响。且在时频分析中明确了主要外诱因降雨与区域滑坡变形的强相关关系及滑坡运动的滞后效应，并可以将相关性和滞后时间定量表示出。这些降雨对滑坡影响的定量解释，为研究滑坡运动提供了新的基础，在对提高今后的滑坡地质灾害防治和预测工作效率方面有重要理论意义和实际价值。
附图说明
[0032] 图1是本发明中某区域尺度上滑坡对降雨响应的分析方法示意图；
[0033] 图2是本发明中某区域内滑坡位移量与降水数据的CWT分析结果示例图。
[0034] 图3是本发明中某区域内滑坡位移量与降水数据的XWT、WTC分析结果示例图。
具体实施方式
[0035] 以下结合附图及实施例对本发明进一步叙述。
[0036] 如图1所示，一种区域尺度上滑坡对季节性降雨响应的分析方法，其包括以下步骤：
[0037] 步骤1、收集目标区域相关资料：
[0038] 收集区域有关的地质地层分布资料，气温、降水等资料，历史SAR影像，以及区域内历史滑坡灾害情况。根据上述资料对区域内滑坡的危害程度做出判断，初步确定目标区域初始范围。
[0039] 步骤2、根据收集的SAR影像数据，基于InSAR技术对目标区域进行识别与监测，获取区域内典型滑坡位移数据。包括以下三个方面：
[0040] 步骤2.1、对研究区SAR影像数据进行前期处理：
[0041] 获取目标区域SAR数据集后，需要对卫星雷达数据进行预处理，来保证后续相关处理方法的准确进行。通过引入精密轨道数据对轨道信息进行修正，可有效去除因轨道误差引起的系统性误差，在SAR影像配准、去斜中起到重要作用。通过引入外部DEM、时空滤波去除部分误差形成相位解缠的差分干涉图。
[0042] 步骤2.2、大区域范围内滑坡灾害的识别：
[0043] 对前期处理过的SAR数据借助于Stacking技术进行大区域内滑坡灾害的识别。Stacking技术在获取M幅相位解缠的差分干涉图后，将解缠相位加权叠加，估计出平均相位变化速率，从而最大程度地减少了大气误差，提高了形变速率地解算精度。通过大范围的Stacking识别聚焦灾害严重的区域作为研究区，这样不仅减少了工作量，提高了效率，而且可以很好的避开了解缠跃变和解缠跳变，然后再对形变严重区进行长时间监测分析来获取更精确的位移结果。这符合地质灾害研究的思路，是普查到详查的体现。
[0044] 步骤2.3、目标区域滑坡的时序分析：
[0045] 在典型的多时相雷达干涉(MTInSAR)算法：PSI和SBAS发展的基础上,本发明运用将两种算法结合的时序InSAR分析方法。PSI+DSI技术是联合永久性散射体（PS）和分布式散射体（DS）来改善形变信号空间采样率，进行相位分析来获取时间序列形变。主要包括三个步骤：PS点的选取、DS点的选取、PS点和DS点的结合分析。
[0046] PS点是指相干性较高的稳定点，是通过振幅离差指数（ADI）选取的。DS点在短时间内保持一定稳定性的目标体，在空间上分布广泛。但是，DS占主导地位的低信噪比像素不能直接用于时间序列分析，而分布散射体在空间分布上统计一致，因此将通过3倍多视滤波后提高的像素点作为DS点。最后将选取的PS和DS点合并在一起，与PS点共享相同雷达坐标的DS点被丢弃，将其在商业软件GAMMA中用最小成本流(MCF)算法进行相位解缠。然后根据解缠后的相位分量各自的特点，逐一计算并去除轨道误差、高程误差、大气扰动和噪声等相位分量。最后获得正确可靠的位移结果。
[0047] 根据目标区域监测的形变结果，均匀地选取几个典型大滑坡的位移数据，通过分析典型滑坡位移量对降雨响应的共性来定量说明区域上滑坡与降雨的关系。
[0048] 步骤3、基于InSAR识别监测结果，结合现有位移资料，进行区域内典型滑坡的现场调查，对区域内滑坡监测精度进行评价。InSAR结果中剧烈形变的位置，发生了严重的破坏，也对应剧烈位移变化。这保证了获取的滑坡位移数据的正确性，从而替代了传统的位移监测，降低了位移数据获取成本。
[0049] 步骤4、利用小波分析工具，探究出降雨与滑坡位移数据间的定量关系：
[0050] 步骤4.1、数据重采样：
[0051] 由于基于SAR卫星获取的滑坡位移数据以重访周期12天为一单元间隔，而对应时间段内降雨数据以一日为单元，为了减少数据误差，这就需要将日降雨数据重采样为以12天为间隔的数据。此外针对滑坡位移数据中由于卫星影像缺失可能存在的个别数据缺失，将根据InSAR原理采用线性插值来补充。
[0052] 步骤4.2、周期项位移数据提取：
[0053] 结合滑坡成因，分析监测的滑坡位移的时间序列，认为其位移的线性趋势主要是重力作用下的非弹性变形，那时序位移曲线上的小波动则应该与外部诱发因素密切相关，而对各诱因进行主成分、相关性分析后，则确定位移周期项是对主诱因降雨的体现。因此，采用时间序列加法模型将滑坡位移时间序列分解为线性趋势和非线性趋势两部分。线性分量通过线性最小二乘拟合计算，非线性分量定义为位移时间序列与之前计算的线性分量的差值，其表现出明显的周期性。因此，降雨对滑坡的影响就集中体现在降雨与非线性位移数据的定量关系上。
[0054] 步骤4.3、连续小波变换（CWT）：
[0055] 由于非线性时间序列位移、降雨数据都表现出明显周期性。将这两个数据借助于CWT转换到时频域，根据其二维图像中信号的时频分布来获取周期性、季节性信息，这即是对数据的时频处理，也是对此分析工具正确可用性的验证，提供了更好的可视化。如图2所示：在95%的置信度下，结果主要反应在锥内；其中高值频谱区域的位置表明了在特定时间存在显著的时间模式（即周期性、季节性）。
[0056] 步骤4.4、交叉小波变换（XWT）和小波相干（WTC）：
[0057] 对连续小波变换后的双时间序列进行交叉小波变换(CWT)和小波相干分析(WTC)，可以得到周期性位移量与降雨量之间的交叉小波能量谱和交叉小波凝结谱(图3)。XWT的计算方法是将一个时间序列的CWT乘以第二个时间序列的CWT的复共轭，结果为二维图像。因此XWT高值的位置表示了两个时间序列中高值共振能量。WTC定义为两个时间序列的CWT在时频域的相干性，它是通过两者之间的标准化互相关来计算的，因此在两个CWT时频图局部相似(即相干)的时频区域内，XWT的绝对值越大、WTC值会很高(即接近于1)；
[0058] 步骤5、CWT结果表明：多年间降雨数据和形变数据都显著存在约1年的周期。此外，在每年的降水湿季，具有突出的动力信号，这是降雨影响下滑坡运动季节性表现。对CWT后的双时间序列进行XWT和WTC处理。图中的箭头表示相对相位关系，右箭头表示相同相位变化（正相关），左箭头表示相反相位变化（负相关）。向上的箭头表示滑坡运动先于降雨1/4个周期，向下的箭头表示其滞后于降雨。通过典型滑坡的一致性规律得到区域尺度上滑坡对降雨相应的相关性和滞后性。
[0059] 最重要的是，图中箭头可表示滑坡位移与降雨量的相关性和滞后性。根据两者的相对相位关系，“→”表示相同的相位变化(正相关)，“←”表示相反的相位变化(负相关)，“↑”表示滑坡位移比降雨提前90天，“↓”表示滑坡位移滞后于降雨90天。在95%置信锥内，区域内典型滑坡一致性规律表示：区域内滑坡形变与降雨之间有很强的正相关关系,这更加清楚有力说明了降雨是区域上地表变形的主导因素。此外，根据观测到的频带箭头走向，发现区域上形变位移与降雨量之间存在的滞后周期，从而可以计算滑坡位移滞后降雨的影响约多少天。
