基于多数据源SBAS技术的滑坡监测和早期识别方法与流程

本发明涉及合成孔径雷达干涉测量技术，具体涉及一种基于多数据源sbas技术的滑坡监测和早期识别方法。
背景技术：
四川省阿坝藏族羌族自治州(简称阿坝州)位于四川省西北部，阿坝州地处青藏高原东南缘，横断山脉北端与川西北高山峡谷的结合部，其地貌以高原和高山峡谷为主。东南部为高山峡谷区，中部为山原区，西北部为高原区。阿坝州地震频发，如2008年“5.12汶川地震”和2017年“8.8九寨沟地震”都是其中影响较大的地震。阿坝州处于构造断裂带、地震带和高山峡谷地带，而受到断裂带的影响，当地的岩体比较破碎和松散，再加上地震所带来的次生灾害影响，在山峦重叠的阿坝州，滑坡地质灾害成为了影响当地居民生命财产安全的主要地质灾害之一。同时，受到强降雨的影响，滑坡地质灾害在阿坝州发生得更为频繁，滑坡地质灾害影响范围大、形变量大、突发性强，常规测量方法(如gps，精密水准测量等)又存在以下几点缺陷：1)外业工作量大；2)需消耗大量人力物力财力；3)布设点受地形影响；4)自动化程度不高；5)多用于点的监测，在大面积监测上没有优势。所以，在滑坡灾害频发的阿坝州，利用新的测量技术对其滑坡的早期识别和滑坡地质灾害的监测、预警成为保护当地人民生命财产安全的一项迫在眉睫的任务。合成孔径雷达干涉测量(interferometrysyntheticapertureradar，insar)是一门测量地形及地表形变的技术，相比较常规测量技术而言，insar技术具有全天候、宽覆盖、高精度、低成本的优势。但是，传统的差分合成孔径雷达干涉测量(differentialinterferometrysar，d-insar)技术容易受到大气扰动、时间失相关和空间失相关等因素的影响，其测量的精度受到了一定的限制，并且无法获得监测区域时间序列的地表形变。而时序insar技术如永久散射体干涉测量(permanentscatterinterferometry，psi)和小基线集(smallbaselinesubset，sbas)等技术及其衍生时序insar技术相继被提出。目前时序insar技术已广泛应用于城市地下水开采、地质灾害(滑坡、火山、地震等)和城市地表等形变测量当中。sbas技术是由berardina和lanari等研究人员提出的一种典型的时序insar技术，与传统的d-insar技术不同，sbas基于短时间基线、短空间基线获取干涉像对，其主影像有多个，在同样的观测时间、影像个数内，积累了尽可能多的干涉图，将时间、视角方面的差异最小化，最大限度地减小了时间去相干、空间去相干对观测结果的干扰。另一方面，相对于psi技术来说，sbas技术的优势在于能够通过奇异值分解(singularvaluedecomposition，svd)的技术，对形变模型进行求解，因此即便在干涉数据集中存在两个不相连的干涉数据子集，sbas方法仍然能够对模型中的形变项、大气影响项、轨道误差项进行估计并去除。sbas技术由于对干涉数据集中的干涉对进行时间基线以及空间基线进行限制，因此干涉数据集中的干涉图能够最大程度地减少时间失相关现象以及空间失相关现象的影响。sbas技术由于能够尽可能高地保证速率反演结果的空间覆盖率，而且sbas在进行影像中的大气相位屏(atmosphericphasescreen，aps)的估计以及去除是通过执行一个时间维上的滤波来实现的。不同于psi技术在对速率反演模型进行求解的时候完成相位解缠操作，sbas技术使用解缠后的干涉图进行模型的参数估计操作。因此，sbas技术因为在解的稳健性方面的优势，在包含较大面积的监测区域的形变提取工作中具有优势。目前，sbas技术虽然在滑坡灾害监测中有着广泛的应用，但其相比较于psi技术而言在高相干性的ps点上精度不高，且在滑坡灾害频发的阿坝州山区，gps和水准点的布设也出现了一定的困难，这也使得sbas技术在得到测量结果之后没有其他测量结果作为比较，其结果的准确性无法得到保证。技术实现要素：鉴于现有技术的不足，本发明所要解决的技术问题是提供一种基于多数据源sbas技术的滑坡监测和早期识别方法，以提高时序insar技术的准确性；该方法基于insar技术和多个(含两个)数据的sar影像，采用小基线集insar方法提取地表形变信息，并在多数据处理结果之间进行相互验证，从而提高时序insar技术在无其他检验条件下的准确性。为了解决上述技术问题，本发明采用如下技术方案：一种基于多数据源sbas技术的滑坡监测和早期识别方法，包括以下步骤：步骤一、获取覆盖监测范围的同一时段的多种sar影像数据；步骤二、公共主影像的选择；步骤三、sar影像的配准；步骤四、干涉对的选取；步骤五、差分干涉图的生成；步骤六、选取地面稳定点；步骤七、估算残余地形相位；步骤八、去除大气相位及残差；步骤九、建立模型反演沉降速率。进一步地，在步骤二中，基于时间序列的sar影像的分析技术，将不同数据源或同源不同图幅的sar影像分开处理，具体操作如下：首先从覆盖同一地区的sar影像中选取公共主影像，将其余影像与主影像配准生成干涉图，计算两幅图像干涉图的相干系数ωi,j，总体用下式表示：其中，tt、bt、dt分别为时间基线、空间垂直基线以及多普勒质心频率的阈值；通过计算每两幅图的相干系数ωi,j，得到一个系数矩阵并对矩阵中的每行求和，即：计算出第i幅图像作为公共影像时的相干系数ωi并比较ωi最大时的影像，此时选取的公共主影像其时间基线、空间垂直基线以及doppler质心频率基线达到最佳组合。进一步地，在步骤三中，在选取出公共主影像之后，将所有影像与主影像进行配准并重采样；当两幅干涉sar图像之间准确配准时，其干涉图中会出现干涉条纹，而insar又通过对干涉相位图的处理来提取地面高程及形变信息，因此影像配准是时间序列insar处理的关键步骤；具体配准操作如下：假设两幅sar复数影像分别为p1(x,y)和p2(x,y)，当两幅影像之间存在整数个数像素偏移量时，其干涉图的公式如下：i(x,y；m,n)＝p1*(x,y)·p2(x-m,y-n)(3)其中，(m,n)表示为整数个数的像素偏移量，p1*(x,y)表示为p1(x,y)的共轭；在图像的配准中，求信噪比snr的表达式为：其中，是干涉图频谱，可由公式(5)求得：其中，m和n分别表示干涉图频谱的行数和列数；当两幅图像精准配准时，其信噪比达到最大值，所以通过遍历所有可能偏移位置上的snr并找出最大snr所在位置及其对应的偏移量，确定干涉相对之间的偏移量，从而完成影像的配准。进一步地，在步骤四中，设定预定的时间基线和空间基线阈值，根据所设定的阈值，事先生成一些干涉条件较好的干涉对，从而初步剔除干涉条件差的干涉对。进一步地，在步骤五中，采用goldstein方法进行滤波以及最小费用流法进行解缠；在时间序列的干涉数据集中，对任意一景干涉图中的干涉相位表示为：φint＝φatmos+φtopo+φdefo+φobject+φorbit+φnoise(6)减去高程相位以后，获得差分干涉相位：φdiff＝φatmos+φtopo-error+φdefo+φobject+φorbit+φnoise(7)对于第j景差分干涉图中，方位向坐标a以及距离向坐标i的像元的干涉相位值用以下的方式表示：式中，j为影像的编号，范围为(1，……，n)，λ为景号信号的中心波长，d(tb,a,r),d(ta,a,r)为a、b两个时刻对于雷达视线方向的累计形变量；表示差分干涉图中残余的地形相位，为大气相位，则表示模型的总的噪声分量。进一步地，在步骤五中，在采用最小费用流的同时，还利用delauny三角网将被分成若干块的解缠结果并入整体当中进行优化，进一步提高结果的稳定性。进一步地，在步骤六中，采用在时间序列sar数据上提取稳健可靠稳定点的方式去除大气相位及残差，提取地面目标点采用以下方式进行：(1)提取时序高相干点：利用sbas方法中平均干涉图，选取时序高相干点，利用设定好的较短的时空基线保证相干点的数量；(2)提取时序高相干点后，在提取的时序高相干点中筛选出在两组或多组不同源干涉图中相同位置的高相干点；(3)选取的地面点满足在不同源干涉图中都为稳定点，同时尽量满足在不同源的sar数据中该点的回波特性相似。进一步地，在步骤七中，面对干涉图中可能会出现的大面积失相干以及相干性偏差的情况，将相干性阈值设置为大于0.18；对于第n景干涉图来说，干涉相位中的低通信号部分表示为如下形式：在上式中，tmn、tsn、a、r分别代表主影像时间、从影像时间、方位向坐标及距离向坐标，s(lp)以及δe(lp)代表形变信号以及地形误差的lp部分，记为噪声部分，而大气相位用φatmo(tmn,a,r)-φatmo(tsn,a,r)表示；此外，λ代表传递信号的波长，b⊥为垂直基线部分，θ为入射角；此后假定s(lp)(t0,a,r)≡0，d为主体；确定s(lp)(tn,a,r)，其中n＝0,…,n作为低通形变的时间序列关于像素在位置(a,r)以及与相对于参考时间t0的计算；在第m景干涉图每个像元的残余相位表示成如下的形式：在这里，v(hp)以及β(hp)分别代表平均速率以及残余形变的非线性部分，δe(hp)为高分辨率的地形误差，以及为噪声项。进一步地，在步骤八中，利用选取多个具有处于干涉相位相对稳定的区域的地面点，求取干涉数据集中每一景干涉图里的大气相位以及噪声项成分；在求取误差项以后，使用高次多项式插值的方法，对研究区域中的大气相位进行插值，达到恢复研究区域中每一景干涉图所包含的大气相位及噪声的目的。进一步地，在步骤九中，估计第n景单视的干涉图每个像元的残余相位(10)中的v(hp)以及δe(hp)最大化了时间相干因子γ(hp)(a,r)，从而得到：其中，δφm(a,r)是设定的相位残差模型，其表达式为：进一步的，对高通信号以及残余形变进行估计；在这里，假定从真实的相位信号高通部分模型误差是局限在(-π,π)之间，能够确定高通相位信号以及残余形变中的非线性部分是能够被简单减去的模为2π的部分；在模型(13)中，应用从第n景单视的干涉图每个像元的残余相位(10)的信号对模型(12)进行估计；模型(13)中，等号右边的已知项是对解缠相位的差分；β(hp)现在从模型(13)中进行反演获得；类比干涉相位中的低通信号部分(9)，模型(13)中的速率部分替换未知的部分β(hp)，即公式(14)中的速率矩阵；基于公式(14)，能够将模型(13)重写为模型(15)，便能够使用svd方法去计算模型(15)中的速率矢量；使用svd方法计算出法方程系数矩阵的广义逆，进而得到速度矢量的最小范数解，从而得到各个时段的形变量；此时，在遵循噪声影响以及解缠误差最小化的情况下解是稳定的，同时允许有效地结合从各个不同的子集中获得的信息。本发明所达到的有益效果是：本发明利用两组或两组以上具有相同覆盖区域、在相同时间段内的不同数据源的sar数据，采用sbas-insar技术分别对不同来源的sar数据进行处理，不同的信号对于地面的散射特性是不一致的，因此对于同一个信号、同一个误差项中误差，由于载波波长所引起的误差的量级也是不一致的。若在包含不同数据源反演结果的研究区域选取一定量的地面点，凭借地面点所反演的相位信息，就能够对这两种不同数据源干涉图相位中的误差项进行去除，因此在对于速率反演模型中的大气相位项、形变速率项进行估计的时候，模型中将会具有较高的信噪比特征，从而得到更为精确的形变速率反演结果。在选取地面点的同时，采用在不同数据源中选取相同地面点的方法，将两组或多组数据联系起来，在结果的时序分析中，可以大大减小多组数据对相同区域进行反演时时序性变结果的差异，达到了在缺省其他检核条件下的高精度、高准确性的反演结果，为滑坡地质灾害的发现和防治提供有效、可靠的技术支持。附图说明图1为本发明实施例的处理流程图。图2为一数据源的sar影像干涉相对图。图3为另一数据源的sar影像干涉相对图。图4为实施区域在一数据中的平均形变速率结果图。图5为实施区域在另一数据中的平均形变速率结果图。图6为实施区域待比较特征点的位置图。图7为特征点p1时序形变比较图。图8为特征点p2时序形变比较图。图9为特征点p3时序形变比较图。图10为发现滑坡发生地区位置图。图11为图10中滑坡位置的放大图。图12为滑坡发生地区时序形变分析图。具体实施方式为使本发明的目的、技术方案及优点更加清晰，以下结合附图及实施例，对本发明进行进一步阐释。应当理解，此处所描述的具体实施实例仅仅用以解释本发明，并不用于限定本发明。实施例：如图1所示，一种基于多数据源sbas技术的四川阿坝州滑坡监测和早期识别方法，包括以下步骤：步骤一、获取覆盖监测范围的同一时段的多种sar影像数据；步骤二、公共主影像的选择；步骤三、sar影像的配准；步骤四、干涉对的选取；步骤五、差分干涉图的生成；步骤六、选取地面稳定点；步骤七、估算残余地形相位；步骤八、去除大气相位及残差；步骤九、建立模型反演沉降速率。在本实施例的步骤一中，对于山区而言，在植被覆盖面积大以及滑坡、地震等地质灾害频发的条件下，为了监测地表量级较大的形变，同时，也为了雷达信号能穿透植被获得更多地表信息，本次实验采用l波段(alospalsar-1)的两组不同图幅的数据对阿坝州理县县城及周边地区。时间跨度为2007年1月4日至2011年3月2日，详细信息见表1、表2。表1：选用alospalsar-1雷达主要参数表中心频率127mhz波段l波段入射角34.3°方位向分辨率3.950099137472467m距离向分辨率6.5579600187500002m表2：alospalsar-1雷达数据日期表在本实施例的步骤二中，基于时间序列的sar影像的分析技术，将不同数据源或同源不同图幅的sar影像分开处理，具体操作如下：首先从覆盖同一地区的sar影像中选取公共主影像，将其余影像与主影像配准生成干涉图，计算两幅图像干涉图的相干系数ωi,j，总体用下式表示：其中，tt、bt、dt分别为时间基线、空间垂直基线以及多普勒质心频率的阈值；通过计算每两幅图的相干系数ωi,j，得到一个系数矩阵并对矩阵中的每行求和，即：计算出第i幅图像作为公共影像时的相干系数ωi并比较ωi最大时的影像，此时选取的公共主影像其时间基线、空间垂直基线以及doppler质心频率基线达到最佳组合。在本实施例的步骤三中，在选取出公共主影像之后，将所有影像与主影像进行配准并重采样；当两幅干涉sar图像之间准确配准时，其干涉图中会出现干涉条纹，而insar又通过对干涉相位图的处理来提取地面高程及形变信息，因此影像配准是时间序列insar处理的关键步骤；具体配准操作如下：假设两幅sar复数影像分别为p1(x,y)和p2(x,y)，当两幅影像之间存在整数个数像素偏移量时，其干涉图的公式如下：i(x,y；m,n)＝p1*(x,y).p2(x-m,y-n)(3)其中，(m,n)表示为整数个数的像素偏移量，p1*(x,y)表示为p1(x,y)的共轭；在图像的配准中，求信噪比(signaltonoiseratio，snr)的表达式为：其中，是干涉图频谱，可由公式(5)求得：其中，m和n分别表示干涉图频谱的行数和列数；当两幅图像精准配准时，其信噪比达到最大值，所以通过遍历所有可能偏移位置上的snr并找出最大snr所在位置及其对应的偏移量，确定干涉相对之间的偏移量，从而完成影像的配准。在本实施例的步骤四中，设定预定的时间基线和空间基线阈值，根据所设定的阈值，事先生成一些干涉条件较好的干涉对，从而初步剔除干涉条件差的干涉对；另外，对于干涉条件较好，干涉对较多的sar影像数据集，可以使用delauney三角网对构网方法选取干涉对；本次实验由于干涉对较少，且干涉条件相对差，所以并未使用delauney三角网构网，干涉相对如图2和图3所示。在本实施例的步骤五中，采用goldstein方法进行滤波以及最小费用流法进行解缠；其中，最小费用流法一种基于网络流的相位解缠方法，它在全局范围内进行搜索并根据流的大小和方向对相位矩阵进行积分，并以高质量区域为基础，向低质量区域积分，从而得到最小问题下的全局最优解。另外，考虑到阿坝州山区众多，干涉条件相对较差，可能出现干涉区域分散的情况，在采用最小费用流的同时，还需利用delauny三角网将被分成若干块的解缠结果并入整体当中进行优化，进一步提高结果的稳定性；与传统的d-insar技术的原理一致，在时间序列的干涉数据集中，对任意一景干涉图中的干涉相位表示为：φint＝φatmos+φtopo+φdefo+φobject+φorbit+φnoise(6)减去高程相位以后，获得差分干涉相位：φdiff＝φatmos+φtopo-error+φdefo+φobject+φorbit+φnoise(7)对于第j景差分干涉图中，方位向坐标a以及距离向坐标i的像元的干涉相位值用以下的方式表示：式中，j为影像的编号，范围为(1，……，n)，λ为景号信号的中心波长，d(tb,a,r),d(ta,a,r)为a、b两个时刻对于雷达视线方向的累计形变量；表示差分干涉图中残余的地形相位，为大气相位，则表示模型的总的噪声分量。在本实施例的步骤六中，采用在时间序列sar数据上提取稳健可靠稳定点的方式去除大气相位及残差，该步骤是形变监测的关键步骤。由于两组或多组数据之间需要进行结果的相互比对，在选取地面稳定点的同时，会将多种数据之间的干涉图联系起来。提取地面目标点采用以下方式进行：(1)提取时序高相干点：利用sbas方法中平均干涉图，选取时序高相干点，利用设定好的较短的时空基线保证相干点的数量；另外，根据步骤三中sar影像配准的像素偏移量计算方法选取时序高相干点，该方法可以提取部分时间跨度内保持高相干性的点，此类高相干点不需要保证在整个影像时域内均保持高相干性，从而增加选点的数量；(2)提取时序高相干点后，在提取的时序高相干点中筛选出在两组或多组不同源干涉图中相同位置的高相干点，此类高相干点在不同干涉图组中均能保持高相干性，可将多组数据联系在一起；(3)选取的地面点满足在不同源干涉图中都为稳定点，同时尽量满足在不同源的sar数据中该点的回波特性相似；此类稳定点不能位于残余地形相位上，需要远离已知的形变、滑坡区域，同时不能位于相位孤岛过着相位跃变的部分。在本实施例的步骤七中，估算形变速率时，考虑了解缠时的相干性阈值，在所有像对中相干性大于该阈值的像元才能参与反演，所以需要根据实际情况设置相干性阈值。在高山众多、植被覆盖区域面积大、滑坡灾害频发的阿坝州，所以面对干涉图中可能会出现的大面积失相干以及相干性偏差的情况，将相干性阈值设置为大于0.18；对于第n景干涉图来说，干涉相位中的低通信号部分表示为如下形式：在上式中，tmn、tsn、a、r分别代表主影像时间、从影像时间、方位向坐标及距离向坐标，s(lp)以及δe(lp)代表形变信号以及地形误差的lp部分，记为噪声部分，而大气相位用φatmo(tmn,a,r)-φatmo(tsn,a,r)表示；此外，λ代表传递信号的波长，b⊥为垂直基线部分，θ为入射角；此后假定s(lp)(t0,a,r)≡0，d为主体；确定s(lp)(tn,a,r)，其中n＝0,…,n作为低通形变的时间序列关于像素在位置(a,r)以及与相对于参考时间t0的计算；在第m景干涉图每个像元的残余相位表示成如下的形式：在这里，v(hp)以及β(hp)分别代表平均速率以及残余形变的非线性部分，δe(hp)为高分辨率的地形误差，以及为噪声项。残余地形相位的估算是利用线关系对像元进行逐个计算，需参考真实dem数据。参考的dem数据越精确，估算的结果越精确。在本实施例的步骤八中，利用选取多个具有处于干涉相位相对稳定的区域的地面点，求取干涉数据集中每一景干涉图里的大气相位以及噪声项成分；在求取误差项以后，使用高次多项式插值的方法，对研究区域中的大气相位进行插值，达到恢复研究区域中每一景干涉图所包含的大气相位及噪声的目的。在本实施例的步骤九中，估计第n景单视的干涉图每个像元的残余相位(10)中的v(hp)以及δe(hp)最大化了时间相干因子γ(hp)(a,r)，从而得到：其中，δφm(a,r)是设定的相位残差模型，其表达式为：进一步的，对高通信号以及残余形变进行估计；在这里，假定从真实的相位信号高通部分模型误差是局限在(-π,π)之间，能够确定高通相位信号以及残余形变中的非线性部分是能够被简单减去的模为2π的部分；在模型(13)中，应用从第n景单视的干涉图每个像元的残余相位(10)的信号对模型(12)进行估计；模型(13)中，等号右边的已知项是对解缠相位的差分；β(hp)现在从模型(13)中进行反演获得；类比干涉相位中的低通信号部分(9)，模型(13)中的速率部分替换未知的部分β(hp)，即公式(14)中的速率矩阵；基于公式(14)，能够将模型(13)重写为模型(15)，便能够使用svd方法去计算模型(15)中的速率矢量；然而，如果多于一个子集被列出，独立方程的数目小于未知的，即系统容易存在秩亏；因此，即便是最小二乘解也是不唯一的，因此需要增加约束才能求得唯一解。使用svd方法计算出法方程系数矩阵的广义逆，进而得到速度矢量的最小范数解，从而得到各个时段的形变量；此时，在遵循噪声影响以及解缠误差最小化的情况下解是稳定的，同时允许有效地结合从各个不同的子集中获得的信息。本次实验两组数据得到的同一区域的平均形变速率图如图4和图5所示，同时，为了探究两组数据结果在时序形变上的一致性，本次实验选取了3个点进行比较，得到了3点在不同的两组影像处理结果里视线方向(lookofsight，los)的形变速率，比较结果如图6至图9所示，从图6至图9中可以看出，两组结果并不存在明显差异，提高了再时序形变监测上的准确性。同时图10至图12为一区域中3点时序形变图，形变图表明该地区在2008年5月至2008年11月中有量级较大的下沉现象，可能存在滑坡地质灾害。、本实施例根据阿坝州地理位置及地形地貌特征，利用已有的同一区域的不同源sar数据，开展了基于小基线集技术的滑坡形变监测，同时增加地面点选点约束条件，加强不同源数据间的联系，缩小各结果之间的差异，达到了在缺省其它检核条件下的高精度、高准确性的反演结果，为滑坡地质灾害的发现和防治提供有效、可靠的技术支持。以上所述仅为本发明的优选实施例而已，并不用于限制本发明。应当指出，对于本
技术领域：
的普通技术人员来说，在不脱离本发明技术原理的前提下，还是可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。当前第1页12