Sentinel-1 SAR数据在四川茂县山体滑坡灾害监测中的应用研究

2019-05-23高思远黄诗峰孙亚勇

中国水利水电科学研究院学报 2019年2期

高思远，黄诗峰，孙亚勇，王慧

（1.中国水利水电科学研究院，北京 100038；2.水利部防洪抗旱减灾工程技术研究中心，北京 100038）

1 研究背景

我国是世界上滑坡灾害最为严重的国家之一，2016年全国共发生地质灾害9710起，其中滑坡7403起，占地质灾害总数的76.2%，与上年同期相比，地质灾害发生数量、造成死亡失踪人数和直接经济损失均有所增加。及时快速获取滑坡灾害的相关信息，实现滑坡灾害预防、监测、治理及灾后评估，对灾中有效地组织救援工作、最大限度地挽救人民生命、降低经济损失具有重要的意义。在滑坡灾害中，快速高效地获取滑坡信息是一切工作的基础，对防灾、减灾和救灾发挥着至关重要作用。卫星遥感技术具有快捷高效、大范围对地观测的优势，已成为一种有效的自然灾害监测手段。根据传感器光谱范围不同，可分为光学遥感和微波遥感两类。由于滑坡灾害发生时往往伴随阴雨等天气条件，光学遥感应用受到很大限制。合成孔径雷达（Synthetic Aperture Radar，SAR）作为微波遥感的一种类型，不受云雨雾等自然条件影响，能够全天时、全天候的获取数据，在滑坡监测应用中得到了学者的高度重视。Weirich［1］利用多幅航拍的SAR数据以及Spot5数据进行滑坡敏感制图比较和研究。Noferini［2］利用SAR数据进行滑坡灾害的缓慢变形的监测。Chihira等［3］人对汶川地震发生严重地形形变的区域进行滑坡判断，并分析了滑坡的分布情况。VanWestern等［4］人以在雷达图像等多种卫星数据的基础上，实现了汶川地震区域内共60000处地震滑坡的的分布图。相比国外，我国对于利用SAR图像提取滑坡信息的研究起步较晚，成功开展应用研究的案例不多，数据处理过程复杂。因此，针对滑坡灾害的应急监测，如何利用SAR数据快速及时提取滑坡体信息成为目前广泛关注的焦点。本文根据滑坡体快速提取的需求，以滑坡体等地物的雷达后向散射原理为基础，提出雷达比值指数法，并结合灾前和灾中SAR图像统计的滑坡区及非滑坡区的实际雷达特征，分析雷达比值指数的合理性，最后，针对单景灾中SAR图像，运用雷达比值指数法开展滑坡体的快速提取。

2 研究区与数据

2.1 四川省茂县区概况茂县位于青藏高原向川西平原过度地带，属四川省阿坝藏族羌族自治州，处于岷江上游中段。地形以高山峡谷为主，地势西北高、东南低；水系为岷江水系，其中岷江由北向南贯穿全境。境内自然地貌特征呈现高山林立、河流深切。茂县地层以马尔康分区地层为主，主要有千枚岩、砂岩、板岩、白云岩、碳酸盐、泥页岩分布。境内茂汶断裂、龙门山断裂地质活动较频繁，岩石节理深度发育，地质条件脆弱。气候干燥多风，降雨年内分布不均，多集中在5-10月，最大日降水量可达75.2 mm，瞬时降雨大，导致茂县地区滑坡、崩塌和泥石流等地质灾变现象普遍发生。2017年6月24日5点45分，四川省阿坝州茂县叠溪镇新磨村新村组富贵山山体突发高位垮塌，该滑坡具体位于岷山南部的富贵山西坡，夹于松坪沟河谷和岷江河谷之间［5］，持续时间100s。茂县滑坡是2008年“5·12”汶川地震后发生的最大规模的滑坡，垮塌方量约1000万立方米，有62户100余人被掩埋，堵塞河道约2公里。

2.2 Sentinel-1 SAR影像及预处理Sentinel-1星座以欧洲及世界其他地区的自然灾害和环境事故快速监测为主要目标，由欧洲委员会（Council of Europe，UC）出资、欧洲航天局（European Space Agency，ESA）设计与开发，位于同一轨道平面内轨道相位差为180°的A和B两颗卫星组成，属于欧洲哥白尼计划（Copernicus Programme；又称Global Monitoring for Environment and Security，GMES，全球环境与安全监测计划）的中重要组成部分。Sentinel-1A于2014年4月3号发射，单颗星轨道周期为12 d。后续A、B两颗星常规重复周期为6 d，高纬度地区最快重访时间为1-3 d。Sentinel-1A/B搭载一台C波段（5.404GHz）合成孔径雷达（Synthetic Aperture Radar，SAR），包含4种成像模式：SM（Stripmap）、IW（Interferometric Wide swath）、EW（Extra-Wide swath）和 WV（Wave）模式，并且 SM，IW 和EW模式分为单极化（HH/VV）和双极化（HH+HV/VV+VH）成像，最大覆盖幅宽可达400km［6］。其中IW模式影像为Sentinel-1的主要模式拍摄数据，采用了TOPSAR（Terrain Observation with Prigressive Scanning SAR）成像技术，由3个条带拼接而成，地距分辨率为20m×22m，覆盖幅宽为375 km。

通过利用2017年6月24日滑坡灾害后单景IW模式、GRD类型的VH和VV双极化方式SAR影像，对茂县地区的滑坡体进行快速提取。针对Sentinel-1A SAR数据，利用欧空局发布的SNAP2.0软件进行数据预处理，主要包括定标、热噪声去除、轨道纠正、滤波、相互配准和正射校正等，最后生成空间分辨率为20m×20m的雷达后向散射系数数据。有关SAR数据介绍如表1所示。

表1 Sentinel-1A SAR数据

2.3 其他数据除了SAR数据外，为了对Sentinel-1数据进行正射校正，本研究收集了研究区30 m分辨率ASTER DEM数据。且为了对比分析滑坡体提取结果，收集了灾后无人机航拍影像。

3 滑坡体雷达遥感监测方法

3.1 雷达后向散射机理森林与裸土作为自然界的两种典型地物类型，与雷达入射波的相互作用具有不同的散射机理。由于入射雷达波性质、植被区地形、植被的覆盖度等多种因素影响，雷达波与植被相互作用的机理非常复杂。森林区雷达波散射既包含冠层顶部的面散射，也包含经冠层衰减的地面反射的面散射和穿透植被冠层后树干与地面之间的二次散射。此外，还有更加复杂的植被冠层体散射回波［7］。与裸土表面相比，森林冠层的体散射具有较强的多次散射和去极化过程。其中，去极化过程是指散射的雷达波的极化方向与入射波极化方向相比发生偏转，极化特性改变，呈现交叉极化回波强度的增强的现象。裸土区土壤的散射特性与雷达波参数（入射角、频率和极化方式）、土壤的含水量及土壤的粗糙度等因素相关，主要呈现为粗糙面散射。裸土区粗糙面散射以单次散射为主，多次散射为辅，雷达回波的去极化现象相对较弱。并且随着粗糙度的增加，后向的单次散射和多次散射增大，土壤粗糙度与后向散射系数间存在近似指数的依赖关系［8-9］，这种依赖关系随着入射角的增加更为显著［10-11］；同时裸土区的雷达后向回波强度与土壤含水量存在正比关系，随着土壤含水量的增加，裸土区后向散射系数增大［12］。

3.2 基于雷达比值指数的滑坡体提取方法根据植被与裸土区的散射原理分析可知，与植被区相比，裸露地表因无植被冠层覆盖，体散射较少，雷达回波的去极化现象较弱。一般裸土区的VH极化后向散射强度较小，VV极化散射强度一般较大，且不存在因植被冠层遮挡引起的地面散射衰减影响。但是，在植被覆盖的山坡发生滑坡时，滑坡区土层裸露，地表粗糙度和土壤含水量增加，VH极化和VV极化后向散射系数可能同时增大，降低了植被区与非植被区的单极化后向散射强度差异。因此，在山体滑坡区，相比单极化方式的VH或VV后向散射系数，VV极化与VH极化的后向散射系数比值应该更能突出裸露地表与植被覆盖区的雷达图像差异特征。为此，针对Sentinel-1雷达C波段的VV和VH双极化数据，提出基于雷达比值指数（Ratio Index，RI）的滑坡体快速监测法。雷达比值指数RI，计算公式如下：

式中：σVV为VV极化后向散射系数，linear形式；σVH为VH极化后向散射系数，linear形式。

4 茂县滑坡体散射特征分析

在茂县滑坡灾害发生前后，滑坡区的地表特征发生变化，由森林覆盖转为粗糙裸露地表。在滑坡前滑坡区雷达散射为森林冠层表层面散射、体散射及经冠层散射的地表面散射，并且森林冠层的面散射和体散射对后向散射强度贡献最大。在滑坡发生后，滑坡区的植被覆盖消失，呈现裸露地表，且地表粗糙度增大，对后向散射强度的贡献主要为裸土地表的粗糙面散射。除此之外，滑坡灾害发生时常伴随持续降水，土壤含水量的增加进一步增大滑坡区的雷达回波强度的变化。因此，为了说明基于雷达比值指数滑坡体监测法的合理性，针对Sentinel-1 SAR数据开展滑坡体的后向散射特征分析。

由于2017年6月12日与2017年6月24日前后都伴随降雨，因此选取与滑坡区植被和地形坡度类似的非滑坡区作为对比，分别统计6月12日与6月24日非滑坡区与滑坡区两类统计区的VH极化、VV极化后向散射系数的差异，如表2所示。

表2 灾前与灾中Sentenel-1 SAR图像后向散射特征统计

由表2可知，对于同类型统计区、同时相SAR数据，VV极化后向散射系数均大于VH极化后向散射系数，且以6月24日灾中滑坡区SAR数据不同极化的差异性最为显著。另外对于非滑坡区同极化不同时相SAR数据，6月24日VH极化、VV极化后向散射系数平均值仅分别略大于6月12日VH极化、VV极化后向散射系数，且差异较小；而对于滑坡区，6月24日VH极化后向散射系数平均值略大于6月12日VH极化后向散射系数，6月24日VV极化后向散射系数平均值明显大于6月12日VV极化后向散射系数。根据公式（2）可以计算不同时相SAR数据同极化后向散射系数变化幅度，计算结果见表3。

式中，pq为极化方式。

表3 针对灾前与灾中SAR数据的不同类型统计区的同极化后向散射系数变化幅度

由表3可知，对于灾前和灾中两时相的SAR同极化数据，滑坡区后向散射系数变化幅度明显大于非滑坡区。进一步分析可知，相对灾前，灾中非滑坡区两种不同极化方式后向散射系数增加幅度都不足10%，表明在森林冠层覆盖的非滑坡区，降雨引发的土壤水含量增大对非滑坡区总后向散射值影响很小。因此，相对灾前，灾中非滑坡森林区后向散射特征未发生显著变化；而相对灾前，灾中滑坡区VV极化后向散射系数增加幅度明显大于VH极化后向散射系数增加幅度。VH极化后向散射系数增加而幅度较小，表明尽管滑坡发生后森林冠层的减少甚至消失，会使VH极化总后向散射系数减小。但降雨引发的土壤水含量增大以及滑坡后土壤粗糙度增大，对后向散射系数复合贡献略大于森林冠层减少的贡献。综上分析，森林冠层的减少、土壤粗糙度和含水量增大均致使VV极化后向散射系数增大，而VH极化后向散射系数影响复杂部分情况下增大、部分情况下减少。因此，从地物的散射特征的角度，VV极化后向散射系数与VH极化后向散射系数比值更能突出滑坡区的图像特征。

5 茂县滑坡体监测

2017年6月24日5点45分，四川省阿坝藏族羌族自治州茂县叠溪镇新磨村富贵山山体顶部震裂山体突然滑动，高速撞击下方坡体，并沿程铲刮坡积块石土层，碎裂解体并转化为碎屑流，协同下部老滑坡堆积体，形成扩散型碎屑流引发高位顺层岩质滑坡。利用Sentinel-1数据，采用雷达比值指数监测方法，对四川省阿坝藏族羌族自治州茂县新磨村地区滑坡灾害进行了应急监测。在滑坡灾害发生前，滑坡体及周围地区植被完全覆盖；在灾害发生后，滑坡体地区高位垮塌，植被覆盖摧毁，伴随着降雨影响的裸露土壤、岩石碎屑物的后向散射特征体现出来。与灾前相比，灾后的滑坡体区雷达散射特征改变，滑坡体的雷达比值指数增大，雷达比值指数结果对比如图1（c）（d）所示，（c）图表示灾后雷达比值指数图，（d）图表示灾前雷达比值指数图。图1（a）（b）分别表示灾后VV、VH极化后向散射强度图，滑坡区均未有显著的体现。图1（a）（b）（d）分别与图1（c）对比可以看出，无论相对单极化的雷达图像或灾前雷达比值指数图，灾后雷达比值指数图表征的滑坡体特征更突出，高亮度显示。

为了从图像特征角度评价雷达比值指数对滑坡体与非滑坡体的识别效果，采用M.Shimada等［13］提出的分离参数R来分析雷达比值指数的优势和合理性，其中分离参数R计算公式如下：

式中：R为分离参数，R值越大说明分离效果越明显；μ1为6月24日（滑坡后）目标参数在森林区的平均值；μ2为6月24日（滑坡后）目标参数在滑坡区的平均值；σ1为6月24日（滑坡后）目标参数在森林区的标准差；σ2为6月24日（滑坡后）目标参数在滑坡区的标准差。

表5表明，与雷达单极化的数据相比，雷达比值指数RI更能突出滑坡体与周围森林区的差异，以此对滑坡体进行快速提取更为合理。最后，根据多次试验确定将RI≥8的部分作为滑坡体范围。另外，提取的茂县滑坡体范围与四川测绘信息局的茂县滑坡无人航拍图对比结果，如图2所示，两者几何形状相似且范围基本吻合，滑坡体长度为2378 m，最大滑坡宽度为1135 m，滑坡体面积为1.3 km2（无人机航拍数据，空间分辨率0.1 m，监测的滑坡体面积为1.62 km2，滑坡体长度约为2658 m，最大滑坡宽度为 1200 m）［14］。

图1 利用雷达比值指数的滑坡体区Sentinel-1 SAR图像对比（GCP1、GCP2、GCP3表示滑坡体所在位置）

表5 灾前与灾中滑坡区的Sentenel-1 SAR图像后向散射特征统计

6 结论

本文以森林和裸土等地物的雷达后向散射机理分析为基础，提出基于VV极化和VH极化后向散射系数的比值指数法，根据单景2017年6月24日灾中Sentinel-1 SAR图像，运用雷达比值指数法开展滑坡体的快速提取和灾害监测。监测和研究结果显示：（1）基于Sentinel-1 SAR数据快速提取的茂县滑坡体结果与四川测绘信息局的茂县滑坡无人航拍图监测结果基本一致，新磨村严重受灾。（2）与航拍图相比，通过精度分析可知，滑坡体面积相差19.7%，滑坡体长度相差10.53%，最大滑坡宽度相差5.42%，总体精度均在80%以上，满足快速提取精度要求。（3）研究结果表明：雷达比值指数能够显著的突出滑坡体的图像特征，该方法在针对Sentinel-1 SAR的滑坡提取与监测上具有较大的应用潜力和优势。（4）采用本文提出的滑坡体提取方法可快速获取滑坡体信息，为抢险救灾和灾后评估提供信息支持。