2017年8月8日，四川省九寨沟县发生7.0级地震，震源深度20 km，地震造成了大量的人员伤亡与财产损失，截止至8月13日，地震造成25人死亡，73000多间房屋不同程度的受损(http://news.cnhan.com[2017-08-24])。九寨沟熊猫海是此次地震受灾最严重的景区之一，震区诱发了大量的地震滑坡，该类灾害具有分布范围广、破坏性强的特点，尤其会引起道路生命线的损毁，造成实地调查与地震救援无法正常开展，严重影响地震应急调查与灾情评估。

地震发生之后，相关单位利用获取的卫星遥感影像数据开展了滑坡与道路损毁程度调查，解译获取了地震极震区漳扎镇地区周边的灾害分布情况。常规的基于遥感数据的滑坡解译主要以目视解译为主，该类方法精度高，但是面对震后覆盖范围广且数量众多的滑坡，其耗时耗力，不能满足地震应急与快速评估的时效性需求。针对此问题，国内外学者开展了遥感图像的滑坡自动识别与解译研究。付文杰和洪金益(2006)基于Landsat TM遥感影像，采用支持向量机的方法自动识别了影像中的滑坡体信息。花利忠等(2008)基于Landsat TM影像数据，结合不同的遥感指数，采用非监督与监督分类的方法实现了滑坡体的自动提取。当震后能获取震区不同时相的遥感影像时，可利用变化检测的方法识别滑坡体(Cheng 等，2004；Nichol和Wong，2005；Park和Chi，2008；李松 等，2010；Hervás等，2003)。以上方法主要基于中低分辨率遥感影像数据，处理的单元为像元，导致信息识别模式具有很大的局限性，阻碍不同尺度对象的关联，且提取的结果椒盐效应明显。

针对中低分辨率滑坡信息识别存在的问题，结合无人机遥感影像的高空间分辨率，面向对象分析方法应运而生并广泛地应用于滑坡信息的识别中。该方法以影像分割产生的同质影像对象作为分析的基本单元，综合利用图斑对象的光谱、纹理、形状、上下文等特征信息，对影像进行分类，实现信息的自动识别。Martin和Franklin(2005)基于Landsat ETM+影像数据，采用面向对象分析方法识别影像中的土滑滑坡与岩滑滑坡，总体识别精度达到65%。Trigila等(2010)采用面向对象的方法提取了意大利北部奥斯塔流域的滑坡体信息。牛全福等(2010)利用不同时相的高分辨率QuickBird影像，分别采用监督分类、密度分割和面向对象分类方法对玉树地震滑坡体进行了快速提取，并对比了提取结果。陈天博等(2017)基于无人机影像获取的DEM与DOM影像，结合地物特征，采用面向对象分类方法实现疑似滑坡区域的信息提取。一些学者还将面向对象与变化检测技术相结合提取震害目标信息。胡德勇等(2008)首先对不同时相的IKONOS和SPOT 5遥感影像进行多尺度分割获取影像对象，之后采用面向对象变化检测的方法提取滑坡体，与传统的基于像素的提取方法相比，提取结果精度更高，更符合现实。任玉环等(2009)基于同一地区不同时相的福卫影像与IKONOS影像，采用面向对象分析方法建立特征规则集分别对影像进行分类，之后采用分类后变化检测的方法提取了地震损毁道路信息。Lu等(2011)利用不同时相的遥感图像，基于面向对象变化检测方法，提取了意大利南部的滑坡点分布信息，提取精度达到81.8%。国内外相关的研究证明了面向对象分类方法针对高分辨率遥感影像能保持图斑的完整性，分类精度得到明显的改善和提高(Lobo 等，1996；Mauro和Eufemia，2001；Giada 等，2003；Huang 等，2003)。该类方法不仅依靠影像的光谱信息，还综合考虑影像对象的空间结构及纹理信息，提高影像信息识别准确度。但是其在分割参数与分类规则设置过程中，需要人为经验与参数不断调整，浪费大量的时间，不能满足实际应用需求。在震后快速应急阶段，需要兼顾识别效率和精度，实现滑坡信息的自动化提取。

鉴于灾情的紧迫性，为提高遥感地震应急时效性，本文旨在快速、准确的提取滑坡灾情分布信息。无人机具有机动、获取影像分辨率高的特征优势，在震后灾情详查中可精确地确定滑坡的空间分布。结合面向对象分析方法，研究影像特征规则集自动构建算法，建立震后短时间内滑坡灾害高精度判识模型，实现九寨沟景区滑坡空间分布的自动识别，提高滑坡识别的准确度与效率，突破获取数据后数人目视解译的工作模式。同时分析滑坡的空间分布特征，为灾后滑坡预测与危险性调查提供信息支撑。

根据中国地震局发布数据，8月8日九寨沟地震震中位于岷江断裂、塔藏断裂和虎牙断裂附近，中国地震局地质研究所徐锡伟等(2017a)根据区域地质构造环境推测，发震构造为塔藏断裂南侧分支和虎牙断裂北段(图1)，图1中活动构造的分布参考了徐锡伟等(2017b)和戴岚欣等(2017)对九寨沟地震震中附近构造的描述及推测。李渝生等(2017)基于野外调查的地表变形，推测可能发震断层位置(图1)。

九寨沟地震震后获取的卫星遥感影像数据受云影响较大，不能完整地解译震区的受损程度，尤其在九寨沟景区的熊猫海附近。基于以上情况，中国地震局地壳应力研究所携带无人机获取了景区的无人机影像。采用无人机对震区影像采集具有机动性强，影像分辨率高的优势，可弥补卫星数据区域不能有效成像的不足，同时解决人工调查不能深入的难题，已被广泛地应用于地质环境与灾害调查、震后应急与灾情评估等领域。本文以九寨沟地震熊猫海附近景区为实验区，以获取的震后无人机航拍正射影像(获取时间：2017-08-12)为数据源，采用的无人机型号为华奥X-UAV。影像获取面积为8.5 km²，空间分辨率为0.14 m。从地震前后影像对比图可见震后研究区内滑坡数量众多且分布范围较广(图2)。

图 1 九寨沟及其周围区域构造格架图 Figure 1 Tectonics frame of Jiuzhaigou and its adjacent area

图 2 研究区地震前后对比图 Figure 2 Comparison images of the study area pre- and post-earthquake

受强烈地震动以及地形影响，坡体松弛甚至解体，从而产生崩滑现象，这是大多数强震次生灾害发生的主要过程，快速滑动是地震次生滑坡的重要特征。黄润秋(2009)在调查了汶川地震滑坡之后，根据斜坡在强震条件下的失稳机制将滑坡划分为溃滑型、溃崩型、抛射型、剥皮型和震裂型5类。强震动作用下，地震波在传播过程中，导致坡体内溃裂破坏，随后，在强地震动持续作用下，坡体沿特定的“面”产生整体的溃散式下滑，相对于常规重力条件形成的弧形特征滑坡，强地震动拉张作用形成的溃滑滑坡主要以“后缘拉裂面”为主。在持续的强地震动作用下，溃裂破坏的坡体在滑床上很难稳定，整体“倾泻”而下形成溃滑，如果前缘没有阻挡，会产生高速远程的滑动。因此，地震滑坡形成的堆积物通常远离滑床，堆积于坡脚前部。在遥感图像上可通过形态、色调、纹理以及阴影等特征对滑坡进行判断。地震滑坡通常对地表覆盖破坏性强，尤其在植被覆盖度较高的区域，植被遭到严重破坏，下方的岩土体裸露在外，从而在影像上表现为色调浅，纹理特征明显，滑坡体所处的地面凹凸起伏(图3(a))。在遥感影像上滑坡体一般呈浅色调或白色调，特别是尚在发展或刚发生不久的崩滑点，由于出露部分多为新鲜结构面，对光谱具有较强的反射能力，使得影像上呈浅色调。在形态特征方面，滑坡在高分辨率遥感影像中呈现舌形、椭圆形、倒梨形、长椅形、牛角形、树叶形、菱形、叠瓦形或不规则形状等，表现为后壁和侧壁组成圈椅状或“双沟同源”等微地貌形态(图3(b))。地震滑坡滑移之后，堆积物向前滑动，常掩盖道路及其他建构筑物，造成道路堵塞、堰塞湖及建(构)筑物的破坏(图3(c))。

图 3 地震滑坡遥感影像示意图 Figure 3 Remote sensing images of earthquake-induced landslide

为提高震后灾情获取的时效性，除提高数据获取效率外，需构建一种快速、准确判识灾情的方法，在获取数据后短时间内快速识别灾情信息。本文提出的方法是在遥感影像获取后，采用多尺度分割方法对影像进行分割构建影像对象图斑；随后，利用特征规则集自动构建方法实现滑坡信息的提取(图4)。

影像对象是由具有相同形状、光谱、纹理、上下文特征的多个像元组成，其对象内部具有最大的同质性，而对象与对象之间具有最大的异质性。地震获取的高分辨率遥感图像中涵盖的地物类型多样，地物尺度大小不一，单一的分割尺度不能完整地识别滑坡信息。多尺度分割充分考虑了地物实体目标的多层次性，采用不同尺度的层次网络结构来揭示地物目标的大小与等级关系结构，克服了数据源尺度固定及单一尺度目标识别局限性的问题。高分辨率遥感图像分割的目的是充分利用影像的多个影像特征，分割之后的影像对象更接近于人眼判别的结果。

多尺度分割采用基于异质性最小原则的区域增长与合并算法，对象的平均异质性可以用色彩(光谱)异质性标准和形状(空间)异质性标准来综合表达。影像对象的平均异质性以光谱异质性h_color参数和形状异质性h_shape参数来计算(Baatz和Schäpe，1999)，公式如(1)所示

式中，w_color是自定义色彩权重，其取值范围为 $0 \leqslant {w_{\rm{color}}} \leqslant 1$ 。形状权重w_shape与色彩权重w_color之和为1。形状异质性由光滑度和紧致度来计算，如式(2)所示

式中，w_cmpct为自定义的紧致度权重，其取值范围为 $0 \leqslant {w_{\rm{cmpct}}} \leqslant 1$ 。本实验中w_color取值为0.65，w_cmpct取值为0.5。

图 4 面向对象滑坡识别与空间特征分析技术流程图 Figure 4 Technology flow chart of landslide identification and spatial distribution analysis based on object oriented method

分割阈值的不同，原始像元聚合形成的对象大小也有差异，分割阈值尺度参数的设定影响分类效果，因此需要选取最优的分割尺度。对于影像中的地震滑坡体，理想的分割效果是在设定的分割尺度下，获取的一个或多个影像对象能够清晰地表现出地物的边界，分割的对象不能太破碎，边界也不能太模糊，不合适的分割阈值对影像分割会造成“分割不足”或“过度分割”的现象。

对影像进行目视解译，发现影像中的地物类型主要包括植被、水体及滑坡，因此设置3层分割尺度来对影像对象进行描述。采用人工设置分割阈值的方式来设置分割阈值，设定的最优分割阈值分别为420，375，230。

特征规则集是利用地震滑坡对象的光谱、纹理、形状等特征，建立不同地物类型之间的特征与阈值集合。规则集构建的效率与完整性决定了后续滑坡灾害信息提取的效率与精度。规则集构建方法主要有目视判读和基于训练样本的特征自动选取，目视判读方法依赖于人为经验，需多次尝试方能选择出合适的特征与阈值，选取过程耗时耗力，且通常选用的特征不完整。基于训练样本的特征规则集构建是通过选择若干个影像对象的特征进行训练，从而获取合适的特征与阈值。本文采用的方法为基于训练样本的特征自动选取方法，由于识别的目标类型主要为地震滑坡信息，识别地物类型较少，因此本实验中采用SEaTH(Seperability and Thresholds)分类规则自动构建的方法来构建特征规则集。该方法基于导出的样本特征值，实现不同类别之间最佳特征的选择及阈值的自动确定，自动化构建特征规则集，提高信息识别效率。

SEaTH方法最早由Nussbaum等(2006)提出，可用于解决在面向对象分析过程中海量特征最优选择及显著特征阈值自动判定的问题。SEaTH方法包含两部分，第一部分是最优特征选择，第二部分是显著特征阈值确定。

(1)最优特征选择。SEaTH方法采用分离度来评估不同类别在某个特征上的关联程度。分离度采用J-M(Jeffries-Matusita)距离来描述，其计算公式如式(3)和(4)所示

式中，B表示巴氏距离，m₁和m₂分别表示两个类别的某特征均值， $\sigma _1$ 和 $\sigma _2$ 分别表示两个类别的某特征标准差。分离度取值范围为[0，2]，0表示不同的地物类别在某一特征上几乎完全混淆，2表示不同的地物类别在某一特征上能够较好地区分开。

(2)显著特征阈值确定。基于高斯概率分布公式来计算类别在某特征上的最佳阈值，其计算公式如式(5)

当 $P(x|C_1) = P(x|C_2)$ 时，C₁与C₂两种类别的混分最少，如图5所示，当特征值选择为X₁时，C₁类别与C₂类别的分离性最好。X₁即可确定为最优特征阈值。

图 5 显著特征阈值确定示意图(Blaschke，2010) Figure 5 Diagrammatic sketch map of determining threshold of salient features(Blaschke，2010)

阈值的计算公式如式(6)和式(7)所示

式中，n₁和n₂分别为C₁和C₂两种类别选择的样本数目，阈值T的值介于m₁与m₂之间。

在第一层次分割尺度上，选择光谱特征与纹理特征为主要特征，应用SEaTH建立植被与非植被识别特征规则集。在第二层次分割尺度上，通过光谱特征与形状特征的组合构建水体与非水体的识别特征规则集。最后在第三层次分割尺度上，综合光谱特征、形状特征和纹理特征建立滑坡体的识别特征规则集，自动确定的特征规则集分布如表1所示。

表 1(Table 1)

表 1 不同层次地物类别特征规则集 Table 1 Feature rule sets of different objects at different levels 层次 地物类型 特征规则集 第一层 植被 Homogenity ≤0.32 and 30 <Band B <146 非植被 Homogenity >0.32 or 30 ≤ Band B or Band B ≥146 第二层 水体 110<Band G <186 and 661<Area <9207 非水体 Band G≤110 or Band G ≥186 or Area ≥661 or Area ≥9207 第三层 滑坡体 Brightness >173.5 and Entory<0.96 and Second Moment >0.76 and Elliptic fit <0.44 非滑坡体 Brightness ≤173.5 or Entory ≤0.96 or Second Moment ≤0.76 or Elliptic fit ≥ 0.44

表 1 不同层次地物类别特征规则集 Table 1 Feature rule sets of different objects at different levels

基于构建的特征规则集，以分割对象为单元提取影像中地震滑坡信息(图6)，提取滑坡总共157处，总面积为1.627 km²，最小滑坡面积为111 m²，溃裂的堆积物将道路整体掩盖，部分滑坡堆积物滑移至熊猫海湖中。

图 6 地震滑坡影像信息提取结果 Figure 6 Automatic extraction result of earthquake-induced landslides

图 7 地震滑坡人工目视解译结果 Figure 7 Visual interpretation result of earthquake-induced landslides

为形象的展示地震滑坡的分布与滑移情况，将获取的无人机影像进行三维场景显示(图8)，便于滑坡的对比与分析。

图 8 地震滑坡三维场景显示 Figure 8 Three dimensional scene of earthquake-induced landslides

为定量化描述滑坡识别的精度，利用相同航片对滑坡进行人工目视判读，获取研究区滑坡分布(图7)，在此基础上，利用混淆矩阵计算用户精度UA(User’s Accuracy)，制图精度PA(Producer’s Accuracy)，总体精度OA(Overall Accuracy)与Kappa系数来衡量识别的准确性(表2)。采用面向对象自动提取方法滑坡体的识别总体精度为94.8%，识别准确度较高。与目视解译结果相比较，由于道路表面存在滑坡体崩滑产生的碎屑，因此将部分道路错分为滑坡体；部分堆积物往前滑移，将部分树木废屑推至湖中，导致在识别过程中误将湖中部分错分为滑坡，从而导致错分与误分。

表 2(Table 2)

表 2 滑坡识别精度分析 Table 2 Accuracy analysis of landslide identification results 类别 PA/% UA/% OA/% Kappa 滑坡体 98.5 95.4 94.8 0.827

表 2 滑坡识别精度分析 Table 2 Accuracy analysis of landslide identification results

地震应急期间，灾情信息的快速识别与评估对于地震应急救援至关重要。为检验本方法的效率，在利用自动分类方法进行识别的同时，利用配置相同的电脑，采用人工目视解译的方法开展滑坡灾害解译，解译者具有丰富解译经验与软件操作基础，两者对比发现，基于面向对象自动分类方法的提取效率更高(表3)。在电脑配置相同的情况下，面向对象的方法仅耗费人工解译一半的时间即可实现整幅图像地震滑坡信息的识别。随着计算机技术的发展，硬件设备能力的提升，自动化信息提取将逐步发展成为地震应急与灾情评估的主要方法。

表 3(Table 3)

表 3 不同方法滑坡识别效率对比 Table 3 The efficiency comparison of different methods for landslide identification 方法 解译总面积/km2 电脑配置 消耗时间/min 人工目视解译 6.2 4 G内存， i5-4570 CPU 90 面向对象 自动分类 6.2 4 G内存， i5-4570 CPU 42

表 3 不同方法滑坡识别效率对比 Table 3 The efficiency comparison of different methods for landslide identification

强震诱发滑坡的空间分布特征对于地震次生灾害风险评估、城镇区划、防灾减灾等具有重要的参考与指导意义。以九寨沟地震熊猫海附近157处滑坡分布为例，从地形、地质构造等方面分析地震滑坡的分布规律。

地震滑坡的形成是由地震动引起的，地震动的作用造成地表变形、山体堆积物松动，斜坡体承受的惯性力发生改变，破坏了滑坡体原有的平衡，促使滑坡体向下滑动及流动。地震波在传播过程中超过了岩体的强度，直接引起边坡岩体的崩塌，大规模的崩塌附着在下方的松散物上，推动了松散物向下滑动，从而形成滑坡。地形因素是孕育滑坡的重要因素，地形因子的统计分析可获取滑坡发生的地形特征，进而挖掘其分布规律。本文选择高程、坡度、坡度变率、坡向、坡向变率、地表粗糙度6个常用的因子作为地形因子，开展地震滑坡在地形空间分布规律的统计研究。坡度变率是地面坡度在微分空间的变化率，是在所提取坡度值的基础上对地面每一点再求算一次坡度，坡度变率表征了地表面高程相对于水平面变化的二次导数。坡向变率可反映地表所有的山脊线、山谷线，是一个能反映等高线弯曲程度的指标，是指在地表坡向提取的基础上，对坡向变化率的二次提取。地表粗糙度是指在一个特定的区域内，地球表面积与其投影面积之比，是反映地表形态的宏观指标。地形因子的详细提取过程可参考汤国安和杨昕(2012)。基于地理信息系统技术(GIS)，根据1:50000比例尺地形图中的等高线分布，建立研究区数字高程模型(DEM)，空间分辨率为20 m×20 m。基于提取的DEM，利用GIS空间分析功能，采用系统缺省参数分别计算获取坡度、坡度变率、坡向、坡向变率和地表粗糙度因子。

为统计地形因子与地震滑坡空间分布关系，采用相等间隔的分类方法将各个地形因子分为5类。然后将滑坡点与分类后的地形因子叠置分析，采用GIS的空间统计分析功能分析滑坡点在地形因子中的分布规律(图9)，发现地震滑坡主要分布在坡度为39.8°—69.8°，高程为2623—2935 m，坡向变率为14.8—44.3，地形起伏度为131—227 m，地表粗糙度为1.16—2.09的空间范围内。地震滑坡的空间分布与斜坡坡度、地形起伏度存在明显的正相关关系，与地表粗糙度存在负相关关系，与坡向、坡度变率等地形因子的关系不明显。许强和李为乐(2010)对汶川地震的统计结果显示，滑坡在坡向上的分布情况受地震波传播和断裂错动作用的影响，并将其归纳为“背坡面效应”与“断层错动方向效应”，灾害体坡向上的空间分布情况在一定程度上可为发震断层的位置提供依据。而在本实验中，并未发现滑坡体空间分布与坡向间的关系，推测是研究区距离发震断裂距离较远的原因，同时也可能受滑坡体数量较少的影响。

图 9 不同地形因子地震滑坡空间分布 Figure 9 Spatial distribution of earthquake-induced landslides under different topographic factors

地震滑坡受控于活动断裂，滑坡的空间分布与距离活动断裂的距离存在一定的关系，统计实验区地震滑坡到震区活动断裂的距离，发现熊猫海附近景区地震滑坡离塔藏断裂最近，距断层的最短距离为8.5 km，且距离断层越近，地震滑坡的分布密度越大，因此，推断研究区地震滑坡可能主要受塔藏断裂控制。

滑坡的空间分布在一定程度上也受地层岩性的影响，因研究区内植被覆盖密度高，无法通过遥感影像获取岩性，因此在本文中未分析与地层岩性的分布规律关系。地震滑坡空间分布规律的分析可确定影响滑坡的诱发因子，进而确定各因子诱发滑坡发生的最优阈值范围，为后续的滑坡危险性评价及滑坡预测提供信息支撑。后续可增大滑坡样本数量及影响因子的范围，进一步确定影响其空间分布的因子图层。

针对传统中低分辨率基于像元的滑坡信息提取方法存在的不同尺度对象关联受限，提取结果椒盐效应明显的问题，结合无人机机动、快速、获取影像分辨率高的特征优势，建立基于多层次分割模型，基于自动构建的多维特征集合，实现震区滑坡灾害空间分布的提取。信息识别精度和效率与目视解译结果进行实验对比分析，本文提出的方法结果均满足地震应急的需求，验证了方法的有效性。同时分析了地震滑坡的空间分布规律，为后续滑坡的危险性评价与灾害调查奠定基础。

利用多尺度分割方法构建影像的多层次分析模型的优势在于，避免了单一分割尺度下难以有效表达不同地物目标的最优分割效果；自动构建滑多维滑坡特征识别规则集合，充分利用了无人机高分辨率遥感影像中纹理与空间信息特征，避免了仅利用单一特征识别的信息遗漏，基于多特征约束实现信息的高精度识别。

通过深入挖掘了九寨沟地震中地震滑坡在孕灾环境、诱灾环境中的空间分布规律，确定了诱发滑坡的重要影响因素，可为后续滑坡危险性评价、灾害调查、预测防治等提供信息支撑，避免了灾害后的盲目治理与风险评估。

本文基于九寨沟地震无人机航拍影像建立的多层次滑坡识别方法在对象分割过程中，仍需要人为的介入，尚未实现真正的自动化处理，在后续的研究中需研究自动化分割算法，进一步提高信息识别精度，真正地在地震应急中发挥作用。同时，滑坡影响因子空间分布规律的分析需要构建震区完整的滑坡点目录，本次地震获取的无人机影像覆盖范围小，因此下一步需要利用卫星影像来要实现震区灾害点的普查，尤其是震后不受天气条件影响的雷达卫星图像。

志 谢 此次震前高分二号影像数据的获取得到了中国资源卫星应用中心的支持，在此表示衷心的感谢！