引用格式：殷幼松，谢酬，李晓芳，等.贵州龙井村滑坡前兆形变的多源SAR影像分析[J].测绘科学,,45(10):55-63.贵州龙井村滑坡前兆形变的多源SAR影像分析殷幼松1，谢酬2，李晓芳1，康永辉3，刘广全1，李杰1(1.环球星云遥感科技有限公司，浙江湖州；2.中国科学院遥感与数字地球研究所，北京；3.中国测绘科学研究院，北京)摘要：为了获取贵州兴义市龙井村大型岩质滑坡的前兆形变信息，该文以RADARSAT-2和Sentinel-1A影像集为数据源，采用时序InSAR分析方法，解译得到该滑坡年初至临滑前的形变速率与时间序列。针对研究区单一影像源形变监测结果存在的离散点密度少、滑坡体时序演化规律刻画不足等问题，该文在充分融合相干分布式目标进行形变反演的基础上，结合两类影像源的空间分辨率和重访周期优势，获取到滑坡发生前周边的形变信息。分析结果表明，龙井村滑坡在监测时段内存在明显形变行为，且其发生前两个月内有加速下滑迹象，时间序列曲线切线角增大，滑动风险逐渐增加。关键词：时序InSAR；龙井村滑坡；RADARSAT-2；Sentinel-1A；形变0引言滑坡地质灾害是我国西南山区最为严重的地质灾害之一，具有发育程度高、危险性等级高、分布广等特点。近年来，针对传统滑坡地质灾害隐患点排查方法面临的诸如人员难到达、环境复杂难以判断、排查效率较低等突出问题，相关专家学者提出了滑坡地质灾害隐患点早期识别的空-天-地一体化“三查”体系，即基于高分辨率光学遥感、合成孔径雷达干涉测量技术(interferometricsyntheticapertureradar,InSAR)的全面普查；利用机载激光雷达、航拍技术开展的重点详查；专业地质人员的定点核查、定期巡查、地面监控设备采样信息的实时报送。就单次滑坡孕育的长期发展过程看，想要实现对其危险状态的成功预警，需要一个长链的监测，而将InSAR技术纳入到滑坡地质灾害监测的“三查”体系中去，能够充分发挥其优势特点，从更为长远的角度进行地质灾害隐患点的排查识别，为后续查核环节提供监测指导及判定依据，进而提升滑坡地质灾害多发区域的预警监测能力。随着InSAR技术体系不断的发展与完善，一些专家学者利用其开展了系列滑坡形变的监测工作。文献[1]利用永久散射体干涉测量（persistentscattererinterferometry，PSI）技术对滑坡是否处于活动状态进行了研究并将其用于新滑坡体的识别中去；文献[2]利用小基线集时序InSAR技术获取了中国三峡巴东地区的滑坡活动，并与库区水位变化进行了耦合分析；文献[3-5]先后利用融合分布式目标的时序InSAR技术对四川省丹巴县甲居滑坡和茂县新磨村滑坡开展了前兆形变研究及监测工作，指出在具有合适SAR影像源情况下，InSAR技术能够获取滑坡体前期形变特征并为其预警提供指导。年2月17日凌晨，贵州省黔西南州兴义市龙井村发生一次大规模顺层岩质滑坡，其体量约百万余方，是贵州省目前垮塌量最大的一例。本文以此滑坡为研究背景，借助于融合相干分布式目标和点目标的时序InSAR分析方法，使用RADARSAT-25m高分辨率和Sentinel-1A20米高重访周期时序影像集对该处滑坡前兆形变特征进行提取，发现在干涉条件满足情况下，综合利用高分辨率和高重访周期的SAR数据源不仅能够实现滑坡体隐患点的空间位置识别，同时也可以得到其临滑前的时序变形规律，这为相关核查和防范工作的开展提供了实施依据。同时，本文研究也为后续西南山区大范围滑坡地质灾害隐患点的排查识别奠定了基础，进一步体现出了InSAR技术在地质灾害“三查”监测体系中的普查应用价值。1研究区域概况研究区（图1）位于贵州省黔西南州兴义市境内，该区域地处南盘江造山构造区，断裂带充分发育，岩土体稳定性弱、完整性差，地灾隐患点分布广泛，尤以滑坡最为严重，是贵州省地灾重点防治区之一[6]。龙井村滑坡最初是在年由道路建设削减山体底部而引起，随后留下一处沿坡方向长约m的光滑圆面，并在滑坡光面东侧形成高约25m的垂直临空面，存在再次坍塌的危险。出于安全考虑，当地随即对该处区域开展周边岩体的稳定性监测。年6月在老滑坡区后缘出现0.9m宽的新裂缝，至年底相关裂缝已完全贯通且变形量不断增大，从年2月13日开始，该滑坡逐渐进入加速变形阶段，每日变形量达到4-5厘米，最终于年2月17日凌晨再次发生大型坍塌。图1龙井村滑坡位置及SAR数据范围Fig.1GeographicalLocationMapofLongjingVillageLandslideandSARDatasetsCoverage2监测数据与方法2.1SAR数据集本次研究收集到加拿大RADARSAT-2卫星获取的超精细宽模式升、降轨时序影像和欧空局Sentinel-1A卫星获取的干涉宽模式升、降轨时序影像。其中，RADARSAT-2升轨影像数为16期，均为灾前影像；降轨影像为3期，包含灾前和灾后。Sentinel-1A升轨和降轨影像数分别为35和9期，且均为灾前影像。SAR数据集的空间覆盖和基本参数如图1和表1所示。表1SAR数据集基本参数Tab.1BasicParametersofSARDatasets传感器名称成像模式轨道方向波长分辨率（Az×Rg）/m入射角/(°)获取日期RADARSAT-2条带模式升轨C-波段（5.6cm）2.5×.67-01-01—-02-13降轨28.13-01-07—-02-24Sentinel-1ATOPS模式升轨20×.32-01-02—-02-14降轨39.32-11-24—-02-.2相干目标时序InSAR分析方法相干目标时序InSAR分析方法通过对具有一定相干性的像元目标开展相位成分分析获取形变信息。而将高相干点目标与分布式目标进行融合处理是针对非城市区域复杂地表环境下反演地表形变的有效处理策略。它能够在充分提高离散点空间密度的基础上，增加解算过程中所构建的三角网络的稳健性以及相位解缠过程中对误差点的抑制。本文采用的时序InSAR分析方法基本流程如图2所示。图2数据处理流程图Fig.2FlowChartofDataProcessing从图2中可以看出，相干目标时序InSAR分析过程中关键处理步骤主要包括高质量稳定点目标和分布式目标的提取、分布式目标的相位优化以及时序相位分析。2.2.1相干点目标提取作为传统PSI技术用来获取形变信息的空间载体，相干点目标的提取方法较多，包括振幅离差法、子视谱相关法、平均振幅阈值法等，该类主要基于时序SAR影像的后向散射强度信息；基于相位信息提取的方法有差分干涉相位序列空间相关性分析法、形变时间模型假设的残余相位阈值剔除法及相干性阈值法等。当前的相干点目标提取通常将两种方法进行组合，以保证相干点目标的质量。本研究采用StaMPS软件处理方式对相干点目标进行候选点集提取[7-8]。1）振幅离差法。时序SAR影像的幅度稳定性在相干性较好的情况下，能够用来表征干涉相位的质量[9]。振幅离差可以用来对点目标质量做第一步筛选。2.2.2相干分布式目标处理在SAR影像上，相干点目标的像元数量与成像场景有很大关系。但通常情况下，这类目标数量远低于分布式目标单元。相干分布式目标的提取利用能够增加监测结果的点密度及空间分布，更加利于形变区域空间形态的判读和形变趋势的掌握。相干分布式目标与点目标散射特性上的差异决定了两者处理方法上的不同。分布式目标除需要考虑提取方法外，还需要对其相位进行优化。1）统计同质像元识别。该步是对其相位进行优化处理的保障。当前处理思想是将其识别过程转化为数学上的统计检验问题，从而利用参数或者非参数假设检验方法进行处理[10-11]。为充分提高运算效率，本研究采用FaSHPS算法[12-14]进行统计同质像元的提取。其思想是将影像集进行时间维平均，以平均后的参考像素值作为真值，在给定置信水平下，将周边待分类像元值与真值构筑区间进行比较，从而获取判别集。值得一提的是，对于分布式目标而言，只有一定相干性的散射单元才能够从相位中提取有用信息，当信噪比较低时，相位中反映出的是杂声相位，如水体及植被茂密区，其干涉相位为随机信号，不再保留任何形变成分。因此，提取用于后续分析的分布式目标时，对统计同质像元识别的结果还需要保留对相干性的考虑[15-16]。2）相位优化。与相干点目标不同，分布式目标通常由一个像元内的多个子散射体矢量叠加构成，就单个分布式目标散射单元而言，受信噪比水平限制，提取结果无法直接应用于时序差分相位分析过程中去，需要进行优化处理[8]。当前对分布式目标的相位优化可以分为两块，一部分是空间域自适应滤波，即对选取的统计同质像元在影像空间域上的滤波处理[17]。该方法主要是通过抑制各像素单元自身的噪声，加强相位中的共有分量，进而提高信噪比。时间域滤波方法不仅要重新组合生成干涉对，同时需要采用迭代方法进行相位的最优解算，因此具有很高的时空复杂度。在采用小基线集干涉构网方式的情况下，可认为小基线对干涉相位为最优相位，不再进行时间域滤波处理。2.2.3形变信息的分离提取利用小基线集网络之间的路径连通性，对未正确解缠的离散点识别剔除。对于大气影响相位的处理，可以利用其在空间一定范围内（1～3km）存在相关性，而在时间维上假定其为随机信号进行时空域的滤波分离。考虑到传统高斯窗口滤波存在的形变与大气分离不彻底、受解缠误差影响大等缺点，本文的时间维滤波针对每个相干目标选择合适的平滑系数，采用平滑样条分离大气效应，以提高大气影响相位估计的准确性[23]。最后将得到的大气相位分量从差分干涉相位中剥离后得到的便是感兴趣的形变信号。3结果与分析3.1多源影像幅度分析分别对两颗卫星的升、降轨数据集进行成像处理，得到图3和图4所示的滑坡体后向散射特征（不规则多边形范围）。从两幅图中可以看出，两颗卫星的降轨影像都受叠掩现象的影响，无法对滑坡体滑变前后形态进行反应和刻画，后续分析中将不再考虑两颗卫星降轨数据的使用。对比RADARSAT-2和Sentinel-1A升轨幅度图可以发现，虽然滑坡体在两种影像上呈现出相近的后向散射特征，但在5m分辨率的幅度影像上，老滑坡留下的悬立面形成的阴影和之前滑落遗留的光面可以清楚辨析，而20m分辨率的幅度影像则较为模糊，失去了对后向散射地物的细节表述能力。（a）升轨（b）降轨图3滑坡体在RADARSAT-2影像上的成像特征Fig.3CharacteristicsofLongjingVillageLandslidefromRADARSAT-2（a）升轨（b）降轨图4滑坡体在Sentinel-1A影像上的成像特征Fig.4CharacteristicsofLongjingVillageLandslidefromSentinel-1A3.2InSAR监测结果与分析1）处理过程。首先采用振幅离差法对RADARSAT-2和Sentinel-1A升轨时序影像集筛选相干点目标，初始阈值设为0.7。利用式（3）提取统计同质像元，搜索窗口根据两种影像分辨率的不同并结合滑坡体空间范围的大小分别设置为7×7和3×3，同时将同质像元个数平均相干系数大于0.35的同质像元选为分布式目标候选单元并对其进行空间域滤波。本文以小基线集组合（RADARSAT-2小基线集空间基线为m以内，时间基线为72d以内；Sentinel-1A小基线集空间基线为d以内，时间基线为36d以内，图5）作为干涉网络连通方式以代替利用不同连通路径对相位开展的加权优化处理，降低分布式目标时间域滤波的时空复杂度。（a）Sentinel-1A（b）RADARSAT-2图5时序影像差分干涉组合Fig.5InterferometricNetworksofTimeSeriesImages获取的相干目标在滑坡体周边的分布情况如图6所示，白色圆点代表分布式目标，黑色圆点代表点目标。从两种相干目标的空间分布看，点目标主要以高亮的后向散射单元为代表，体现出主体散射的特征，而分布式目标幅度背景较暗，散射特征不明显。受阴影现象的影响，靠近本次滑坡体一侧顶部和上部区域并未提取到相干目标（图6（a）中白色不规则多边形范围），但两种数据源在裸露岩面及滑坡体中下部都有很高的点密度。（a）RADARSAT-2（b）Sentinel-1A图6相干分布式目标和点目标在滑坡体周边的分布Fig.6DistributionofCoherentDistributedandPointTargetsamongLandslide利用式（2）对选取的相干目标候选点计算时间相干性，并依据相干性值直方图分布区间将该阈值设为0.55，以进一步筛选高质量散射单元。同时考虑到滑坡体周边的地形特征，将外部DEM非系统残差的搜索空间设置为30m，即构建三角网解算时，认定相邻目标点有效高差小于30m。对小基线集时序差分干涉相位进行3维迭代解缠后，将其转换到单主影像空间下，计算两者之间由闭合环解算不一致产生的相位残差。考虑到多视处理操作对干涉相位造成的非线性影响将会使得解缠相位闭合环与理论值有微小差异，本文将相位残差容忍值设为[21]，获取最终用于解译的相干目标，两种影像源的点数见表2。表2两种影像提取的高相干点个数Tab.2Numberofpointsselectedbytwodatesets监测方法相干目标个数/个RADARSAT-2Sentinel-1A相干目标）解译结果及分析。考虑到研究区内滑坡体周边局部高度差异不大，选择将解缠后的长短空间基线干涉对一起进行DEM系统残差解算，并重新进行迭代相位解缠处理，直至解算出的DEM系统残差完全收敛。最终获取的滑坡体周边线性形变速率如图7所示。从图7中可以看出，RADARSAR-2高分辨率时序影像结果中有3处存在明显形变迹象（方框范围），其中P1位置即为龙井村滑坡，可以从形变速率结果中进行很好的甄别。而对比Sentinel-1A时序影像结果，可以看出这3处区域在空间上表现出相同的形变态势，但无论量级上还是离散点的密度上都与RADARSAT-2影像结果存在一定差距，难以从形变速率上对隐患点加以判别。就本例而言，这里体现出高分辨率影像在滑坡隐患点的面上排查识别方面较于中等分辨率影像有更高的探测能力、更低的虚警率，从而可以提高解译结果的可靠性。a）RADARSAT-2（b）Sentinel-1A图7滑坡体周边形变速率（矩形框为滑坡体位置）Fig.7DeformationRateamongLandslide(RectangleisthePositionoftheLandslideBody)此外，对龙井村滑坡区域离散点形变信息进一步分析时发现（图8），两种数据源经过时序InSAR分析获取的明显形变区域位于滑坡体同侧中下部，而发生坍塌位置的顶部和东侧上部并不突出。该情况出现的原因可分为两种：一是顶部和东侧上部被植被充分覆盖；二是滑坡发生前顶部和东侧上部位移量较大。这两种情况导致了失相干现象的发生，进而使得InSAR技术未捕捉到形变信号。滑坡体中下部未受顶部悬立面遮挡影响且有部分裸岩，为相干目标的提取提供了环境条件。a）RADARSAT-2（b）Sentinel-1A图8滑坡体形变特征（椭圆形圈内为滑坡体位置）Fig.8DeformationRateofLandslide（theLandslideLocationisintheEllipticalCircle）为了充分体现两类数据源监测结果之间的互补性，绘制解译区域的离散点图，如图9所示。从图中可以看出，RADARSAT-2高空间分辨率监测结果存在3处明显的形变区域，分别对应上文提及的3处隐患点位置，且兴义滑坡的不稳定迹象可以从中清晰辨别。而在Sentienl-1A监测结果中，虽然滑坡体形变现象有所体现，但由于离散点数过于稀少，在尚未确定隐患点位置的情况下，很难判定其识别的可靠性，容易造成漏提和误选。（a）RADARSAT-2（b）Sentinel-1A图9滑坡体形变速率离散点图（长方形方框内为滑坡体位置）Fig.9ScatterPlotsofLandslide（thelandslideLocationisintheRectangularBox）为了准确获取滑坡体周边临滑前的时序形变演化特征，消除个别奇异点对监测结果的影响，对图8中的圆框范围内相干点位绘制过程曲线（图10）。从图10中可以看出，在滑坡隐患点时间序列演化规律的获取上，具有较高重访周期的Sentinel-1A监测结果优势突出。虽然两类数据源监测结果在整个研究时段内对滑坡体形变规律的刻画较为一致，但密集的时间维影像采样结果对滑坡体形变历程的反应更加完整、全面，有利于滑坡体发展趋势的预判，对提前预警起着关键作用。图10滑坡体形变点时序曲线Fig.10TimeSeriesCurvesofLandslide此外，图10的时序演化曲线同样反应出了滑坡体不同形变阶段的速率变化。在匀速形变阶段（年1月至年1月），滑坡体的形变速率约为28mm/a，而处于加速阶段（年01月至年2月中旬）时形变速率约为mm/a，增加了6倍之多。这与相关文献中对岩质滑坡的形变历程认知较为一致[24]。为了更好反应这种形变阶段的过度，引入滑坡体临滑预报判据中的位移切线角概念作为InSAR监测结果中滑坡隐患点所处阶段的判定标准，而位移切线角常指的是某一时刻形变曲线的切线与坐标横轴之间的夹角[25-26]。由于横纵坐标量纲的差异会导致计算出的位移切线角存在较大区别，本文计算位移切线角时纵坐标以“毫米”为单位，横坐标以“月度”为单位。与滑坡体临滑前的地面监测估计位移切线角时横坐标采用以“天”为单位不同，这里主要考虑的是InSAR监测结果形变量级会更小一些，且监测结果的时间维采样尺度会更长。对图10的两个阶段计算位移切线角，匀速形变阶段的位移切线角为-66°，而加速形变阶段的位移切线角则为-85°（符号的正负取决于形变量方向的选择），可见从均匀形变到加速形变位移切线角发生了明显变化。在依据地面设备监测的数据对滑坡临滑时刻进行判定时，大量滑坡实例的统计分析得出，位移切线角范围值通常为89°～89.5°，而对于InSAR监测结果计算的位移切线角而言，由于滑坡隐患点类型、所处地理环境、孕育成因等存在不同，各滑坡隐患点的位移切线角数值相差会很明显，但不变的是位移切线角都会随着形变的加速而逐渐变大。因此，滑坡体临滑预报判据中的位移切线角同样可以作为InSAR监测结果反应滑坡所处形变阶段的判据。从上述分析中可以看出，两种不同分辨率和重访周期的影像源在滑坡地质灾害的识别及形变发育过程获取方面充分表现出了各自的优势特点。所采用的RADARSAT-2影像源有着更高的分辨率，能够利用滑坡隐患点空间上存在的形变特征对其加以排查甄别；而Senetinel-1A影像源具有很快的重访周期，对识别的滑坡体能够追溯其形变演化历程，捕获异常阶段，进而判断滑坡可能发生的时段。4结束语本文对加拿大RADARSAT-2和欧空局Sentinel-1A两颗卫星多期影像集开展融合相干分布式目标和点目标的时序InSAR分析，充分借助于两种数据源各自的高分辨率和快重访周期优势，获取近期发生的兴义市龙井村滑坡前兆形变特征及时序演化规律。就该例而言，InSAR技术获取的形变信息与当前人类对于滑坡体形变演化规律的基本认识相一致，即大型滑坡体通常具有一定的前兆形变迹象，且随着临滑日期的接近，加速滑动的趋势表现得愈加明显。以本次研究为基础，后续将进一步结合RADARSAT-2和Sentinel-1A卫星的广域监测能力，对西南山区滑坡地灾隐患点开展高检出、低虚警的排查识别工作，充分凸显InSAR技术在地质灾害“三查”监测体系中的普查作用，为“人防+技防”的专业地质灾害监测体系的优良建设提供坚实的保障。致谢：河海大学蒋弥教授提供的FaSHPS软件包、加拿大MDA公司为本研究提供的RADARSAT-2卫星数据、欧洲太空局提供的Sentinel-1A数据和POD轨道数据，在此表示感谢。作者简介：殷幼松（—），男，安徽淮南人，工程师，硕士，主要研究方向为时序InSAR技术工程化应用。E-mail：.
转载请注明：http://www.longjingzx.com/ljly/8103.html
