雷达数据辅助下的滑坡易发性评价

赵金童,牛瑞卿,姚 琦,武雪玲

中国地质大学地球物理与空间信息学院，武汉 430074

0 引言

三峡工程作为人类重大工程活动对库区地质环境影响巨大，滑坡作为典型的地质灾害在库区屡屡发生，对当地居民生命财产安全和三峡大坝正常运营构成了严重危害。因此，对三峡库区范围进行滑坡灾害的易发性评价意义重大。

随着遥感、地理信息系统和计算机技术的发展，对于滑坡的易发性评价，国内外学者已经开展了多方面的研究。例如：Carrara[1]利用回归分析对意大利山区滑坡进行易发性评价；Lee[2]基于GIS技术通过逻辑回归模型进行滑坡易发性预测；Ballabio 等[3]采用支持向量机模型对意大利北部区域进行了滑坡易发性制图；晏同珍等[4]、殷坤龙[5]多次对三峡库区进行滑坡灾害空间预测；牛瑞卿等[6]利用粗糙集-支持向量机模型实现了三峡库区秭归到巴东段滑坡易发性分区；彭令等[7]、武雪玲等[8]分别在多源数据支持下得到了三峡库区滑坡稳定性分区图。虽然上述研究已经使滑坡空间预测达到了一定的精度，但其中选取的湿度指数因子都只是间接地反应了岩土体含水量，再通过岩土体含水量间接反应地下水信息，而缺少因子来直接表示地表水文状况。

随着微波遥感技术的持续发展，各种适用模型不断改进，大面积反演岩土体含水量已经成为了可能。Paloscia等[9]研究表明，基于神经网络算法利用Sentinel-1影像反演岩土体含水量(soil moisture content，CSM)是可行的；何连等[10]利用多时相Sentinel-1数据实现了农田地表岩土体水分的反演；Oh等[11]根据不同极化方式下粗糙度和介电常数与后向散射系数的相关关系提出了反演岩土体含水量的经验模型；Fung等[12]提出积分方程模型(integrated equation model，IEM)，该模型基于辐射传输方程建立地表散射模型；Fung等[13]对IEM进行了改进，得到一种高级积分方程模型(advanced integrated equation model，AIEM)； Ulaby等[14]提出了基于微波辐射传输方程一阶解的微波散射模型(michigan microwave canopy scattering model，MIMICSM)，主要针对植被覆盖的粗糙地面，是目前应用最多的植被覆盖地区散射模型，但其较为复杂；Attema等[15]提出的水-云模型将植被看作各项散射均质体，模型简单实用。

Hassaballa等[16]分别利用RADARSAT-1雷达数据和Landsat-7 ETM+光学数据反演吉隆坡Bukit Antarabangs滑坡地区地表岩土体含水量，研究其与该地区历史性滑坡的空间相关性，发现利用雷达数据反演得到的含水量与滑坡空间相关性的受试者工作特征曲线下的面积(area under curve，AUC)为0.7且高于光学数据，证明地表岩土体含水量高的地区更易发生滑坡且雷达数据相对光学数据能更好地提取岩土体含水量；Brocca等[17]引入散射计(ASCAT)数据衍生得到的岩土体含水量作为水文因子来预测意大利中部地区Torgiovannetto滑坡的张力裂缝，结果表明其相关系数高于采用降雨量作为水文因子的预测结果，也就意味着岩土体含水量与滑坡的发生具有较强的相关性。

本文在前人研究的基础上，选取三峡库区秭归段作为研究区，利用Sentinel-1雷达数据通过水-云模型去除植被影响得到裸露地表后向散射系数，采取同时相同范围不同极化方式的Sentinel-1雷达数据获取地表粗糙度参数，再利用前人在AIEM的基础上建立的后向散射模型反演得到岩土体含水量；以此作为反映地表水及地下水信息的因子替代其他湿度指数，以实现研究区滑坡易发性评价。

1 岩土体含水量反演和逻辑回归

1.1岩土体含水量反演

1.1.1 水-云模型

水-云模型(water-cloud model)由Attma等[15]于1978年提出，主要为了简化植被覆盖区的后向散射机制，提高模型适用性。该模型将植被层看作一个各项均质散射体，忽略了植被层及地表之间的相互多次散射，将植被覆盖地区总的后向散射简单描述为两部分：由植被层直接反射回的体散射项；经植被层二次衰减后地表的后向散射项。模型表述如下：

(1)

(2)

γ2(θ)=exp[-2BCvegsecθ]。

(3)

1.1.2 岩土体含水量反演模型

决定裸露地表后向散射系数的参数主要是CSM和地表粗糙度，而地表粗糙度通常由均方根高度(s)在垂直方向相关长度(l)在水平方向上的组合表示。为了简化模型参数， Zribi等[18]提出利用地表粗糙度参数Zs来表示粗糙度，Zs=s2/l。

李震等[19]通过高级积分方程模型的模拟数据，分别建立了同极化雷达后向散射模型和双极化雷达数据粗糙度计算模型，然后将两种模型结合即可求得CSM：

(4)

(5)

1.2 逻辑回归

逻辑回归模型(logistic regression model)是一种常见的机器学习模型，因其便捷高效，得以普遍应用。因为滑坡只有发生或不发生两种情况(一般1代表发生，0代表未发生)，本质上它是一个二分类问题，所以二元逻辑回归模型可以完美地解决滑坡易发性评价中的二分类问题[21]。

我们假设滑坡发生的概率为P，不发生的概率则为1-P。以P作为因变量，xi为自变量，建立回归方程：

(6)

式中：α为常数；βm为滑坡各影响因子的逻辑回归系数；xm即代表滑坡各影响因子的相对权重；m为影响因子数目。

2 研究区概况及数据源

2.1 研究区概况

研究区位于长江三峡库区秭归县境内，地理坐标为110°30′00″E—110°51′00″E， 30°55′00″N—31°04′00″N，面积约为293 km2，其中水系面积约41 km2。研究区位于湖北省西部、三峡库区库首位置(图1)，地势西南高东北低，海拔80～1 500 m，沿江两岸高山林立，为峡谷地貌；气候属亚热带大陆季风气候，雨量充沛。研究区内滑坡灾害频发。

2.2 数据源

本文主要采用的数据源有：1)图2a 为Sentinel-1A IW(interferometric wide swath)模式下 Level-1 地距影像(ground range detected，GRD)数据，编号为S1A_IW_GRDH_1SDV_20150826T103524_20150826T103557_007433_00A3CA_0028，用于提取地表岩土体含水量；2)图2b 为Landsat-8 OLI 卫星影像一景，编号为LC81250392015216LGN00，用来提取地表覆被信息；3)1∶5万比例尺地形图，用于提取地形、地貌等信息；4)1∶5万及1∶20万比例尺地质图，用于提取地层岩性、构造等信息；5)其他历史滑坡资料和野外调查资料，用于已发生滑坡的解译。

Sentinel-1A作为欧空局(European space agency ，ESA)哥白尼全球对地观测计划的首颗卫星于2014年4月发射，主要用于取代Envisat卫星提供全球范围的高分辨率C波段雷达数据。本文选取的数据为Level-1 GRD产品，它经过多视处理并经WGS84椭球投影，只有幅度信息而没有相位信息。Landsat-8 由美国国家航空航天局(national aeronautics and space administration, NASA)于2013年2月发射，其携带的OLI(operational land imager)陆地成像仪除了ETM+传感器包含的7个波段外还增加了蓝色波段和短波红外波段。本文所使用的遥感数据参数信息见表1。

3 岩土体含水量提取

3.1影像预处理

对Landsat-8影像处理首先进行辐射定标、大气校正等预处理，然后获取归一化差分水体指数(normalized difference water index，INDW)、归一化植被指数(normalized difference vegetation index,INDV)和土地覆被等信息，主要通过ENVI软件实现；对于Sentinel-1雷达数据的处理主要包括辐射校正、斑点滤波和几何校正，目的是获取植被覆盖地表总的雷达后向散射系数，主要通过ESA开源软件Sentinel application platform 实现。

图1 研究区位置示意图Fig.1 Location of the study area

a. Sentinel-1A雷达影像；b.Lansat8影像。图2 研究区遥感数据Fig.2 Remote sensing data of the study area

影像类型成像时间分辨率/m模式入射角/(°)频率/GHz极化方式Sentinel-12015-08-2620×20IW GRD32.95.405VH+VVLandsat-82015-08-0430×30OLI

注：IW GRD.干涉宽幅模式。

3.2 水-云模型应用

植被含水量即为单位面积植被中所含水的质量。作为水-云模型中一个重要的参数，其大面积提取比较困难，通常利用Landsat等陆地影像进行反演。2004年Jackson[22]发现植被含水量Cveg和Gao[23]提出的归一化差分水体指数INDW有很强的联系，即：

Cveg=1.44INDW2+1.36INDW+0.34；

(7)

(8)

其中，RNIR、RSWIR分别代表Landsat-8影像的5、6波段的灰度值。

图3反演得到的含水量只是浅层地表的岩土体含水量，但Uhlemann等[25]通过连续3年的地电监测数据发现滑坡重新激活前浅层岩土体含水量显著提升，同时邓孺孺等[26]的研究也发现深层岩土体含水量与浅层岩土体含水量正相关；两者研究表明浅层地表含水量也显著影响着滑坡的发生且与深层岩土体含水量成正比。虽然深层岩土体含水量对滑坡稳定性具有更重要的意义，但在现有技术手段下反演地表岩土体含水量可能是更现实的选择。

4 滑坡易发性评价

由于本文数据来源多样，坐标系、投影不尽相同，首先需要对数据统一至相同投影坐标系下。本文选取栅格单元作为评价单元，需将地质、地表覆被等矢量数据转为栅格数据，并将栅格单元大小统一至30 m×30 m。研究发现，土质滑坡较岩质滑坡更易受水的影响[27]，所以选取土质滑坡作为研究对象。

4.1 评价因子选取

滑坡的发生主要受控制因素和影响因素两大类因素的影响。控制因素指的是研究区的地形、地质构造等因素；影响因素指的是人类工程活动、地表覆被等因素。本文在前人研究的基础上结合研究区的特点选取地质(工程岩组指数(IPT)、断层距离指数(IDF)、斜坡结构指数(ISA))、地形地貌(数字高程(D)、地形粗糙指数(ITR)、地形表面纹理(ITST)、坡度(i))、地表覆被(归一化植被指数(INDV)、土地利用指数(ILU))和水文(水系距离指数(IDR)、地形湿度指数(ITW)、岩土体含水量(CSM))四大类共12个评价因子。为了确保各因子间的相对独立性，利用Pearson积差相关系数法对因子间相关性进行评价(表2)，发现相关系数的绝对值均小于0.75，可以作为因子评价滑坡易发性。

图3 岩土体含水量分布图Fig.3 Map of soil moisture content

4.2 因子归一化处理

因所选的评价因子处于不同类型、范围和数量级下，归一化处理可以提高数据处理效率，改善模型精度，所以对于连续型因子，我们采用极差标准化方法对其归一化；对于离散型因子(包括斜坡结构、土地利用和工程岩组)，我们根据已解译出的滑坡栅格单元中各因子不同属性占总滑坡栅格单元的比例作为各属性的值，再采用极差标准化方法对各属性值进行归一化，结果见表3。

表2 各因子Pearson相关系数

表3 离散型因子归一化处理

4.3 建立滑坡易发性评价模型

基于逻辑回归模型的滑坡预测主要思路就是选取合适的训练和测试样本建立模型并验证模型精度，进而将模型应用到整个研究区，得到研究区滑坡易发性分区图。本文分别使用岩土体含水量，以及前人常用的地形湿度指数为地表湿度指数因子进行易发性研究，并对不同湿度指数因子下易发性评价的精度进行分析，验证采用岩土体含水量因子进行评价的可行性。主要步骤如下：

1) 划分栅格单元。本文采用的数据基本上都是以栅格或矢量形式存在的影像数据，同时结合遥感影像的空间分辨率，将所有评价因子重采样为30 m×30 m大小的栅格单元。经统计可得，研究区共划分为280 670个栅格单元，其中历史土质滑坡单元6 830个，未发生滑坡单元273 840个。

2) 选择建模样本。我们从所有的滑坡栅格单元中随机抽取80%(5 464个栅格单元)和相同数量的未滑坡栅格单元作为训练样本，剩下20%(1 366个栅格单元)滑坡栅格单元作为测试样本。

4) 滑坡易发性评价。利用建立的模型对研究区全部栅格单元进行计算，结果输出即为滑坡易发性指数。

5) 结果分析。通过测试样本对滑坡预测结果进行评价，并对比分析采用岩土体含水量、地形湿度指数因子时的预测精度。

5 结果分析和验证

对模型输入整个研究区范围的因子，输出即为介于0.0～1.0之间的滑坡易发性指数，指数越大，滑坡越容易发生，反之滑坡越不容易发生。本文采用自然断点法对滑坡易发性指数进行分区，划分为不易发区(0.0～0.2)、低易发区(0.2～0.5)、中易发区(0.5～0.8)、高易发区(0.8～1.0)，得到岩土体含水量因子下滑坡易发性分区图(图4)，并对各分区所占面积比(表4)进行分析。

()内为滑坡易发性指数。图4 岩土体含水量因子下滑坡易发性分区图Fig.4 Land slide susceptibility zoning map based on the CSM

易发性指数易发性分区所占面积比/%滑坡所占面积比/%0.8～1.00.5～0.80.2～0.50.0～0.2高易发区中易发区低易发区不易发区9.714.916.359.153.635.19.61.7

预测结果显示，土质滑坡多发生在河流沿岸，尤其在长江支流童庄河、青干河、吒溪河沿岸易发性指数较高，表明土质滑坡受水的影响较大，基本符合现实滑坡分布状况。

由表4可知，中高易发区面积仅占研究区总面积的24.6%，而区内滑坡面积却占到了总滑坡面积的88.7%；这也间接表明了我们采用岩土体含水量因子进行滑坡易发性评价是可行有效的。

选择研究区内典型土质滑坡新滩滑坡进行预测结果验证。根据已有历史资料和实际调查资料，新滩滑坡位于屈原镇西北2 km处长江北岸，与链子崖危岩体隔江相望，为古滑坡。最近一次滑动发生于1985年6月12日，滑坡体积约3 000万m3，南北长约2 km，滑动面积约为1.1km2；由于预报及时，新滩镇居民无一伤亡。在滑坡易发性分区图中，该地区被预测为高易发区，如图5所示。

为了检验2个不同湿度指数下滑坡的预测能力并进行对比分析，我们利用测试样本采用Chung等[28]提出的成功率验证法测试2种模型的预测精度(图6)。图6中：横轴代表滑坡易发性从大到小分成的100个区间，如横轴10.0%处即表示滑坡易发性指数为90.0%～100.0%的区间；纵轴即为对应易发性区间下实际滑坡累计发生的频率，曲线下面积越大，证明预测效果越好。

由图6可知，对于本文预测的滑坡易发性指数为80.0%～100.0%的区间(横轴20.0%处)，即对应易发性分区中的高易发区，采用岩土体含水量、地形湿度指数分别预测了64.6%、59.4%的实际滑坡落在高易发区中，采用岩土体含水量、地形湿度指数因子时曲线下面积即模型预测精度分别为80.2%、77.2%。按照滑坡易发性指数和滑坡累计发生频率将地形湿度指数由数字高程模型中提取的坡度、流域面积计算得到，它主要考虑了地形对地表水的再分配作用来描述岩土体水分的空间分布。但岩土体水分的空间分布还受植被、岩土体特性和人类活动等的干扰，这也决定了地形湿度指数因子只能间接反应岩土体含水量这一重要水文因子对研究区斜坡稳定性的影响。而通过雷达数据可以直接反演出去除植被影响后的裸露地表岩土体含水量，它也与区域范围内地下水的分布状况有很强的联系。同时滑坡预测结果也证明：利用岩土体含水量因子比地形湿度指数因子有着更高的预测准确率，意味着它作为直接反应岩土体含水量的水文因子比其他湿度指数能更好地反应地表及地下水信息，可以作为滑坡易发性评价的重要水文因子。

图5 典型滑坡验证Fig.5 Typical landslide validation

图6 滑坡累计发生频率与易发性指数曲线Fig.6 Relation between landslide accumulative frequency and prediction index

6 结论与展望

1)通过C波段Sentinel-1雷达数据反演获取岩土体含水量，将其作为一种新引入的水文因子代替地形湿度指数来表征地表及地下水信息。在此基础之上，再充分利用多源空间数据提取其他因子并对数据进行归一化处理，分别建立逻辑回归模型，定量获取滑坡易发性指数，得到研究区易发性分区图。

2)本文利用雷达数据获取岩土体含水量因子进行滑坡预测的精度要高于传统方法采用的地形湿度指数因子，证明岩土体含水量作为重要的水文因子引入到滑坡预测中是可行的并且精度较好。但由于岩土体含水量的反演受地形及地表覆被的影响较大，现有模型方法普适性较差，大面积反演水文因子仍存在较大困难。

3)由于条件限制，本文并未测定研究区实地岩土体含水量，如何相对准确的大面积获取岩土体含水量必然是后续研究的重点。随着遥感技术的发展，滑坡空间预测的突破点就在于如何结合利用多源多类型的遥感数据。