基于高分遥感影像的泥石流潜在形成区识别方法探讨

苏轶君

(福建省测绘院，福建 福州 350003)

0 引言

泥石流山地灾害会对人类生产和生活带来不同程度的损害，严重的泥石流会直接威胁人类生命和财产安全[1]。若能识别出泥石流灾害可能发生和造成影响的区域并提前布设防范措施，则能有效地减小灾害损失，但由于影响泥石流灾害发生的因素比较多，因此这方面的研究一直存在比较大的困难。

早期国内外学者对泥石流的探究主要在泥石流基本特征、形成机理和危害性方面[2-7]，20世纪70、80年代开始逐渐从理论到应用层面引入遥感技术进行泥石流灾害的探究，且引入监督分类、阈值法等遥感分类技术实现泥石流物源区自动化提取，也有通过多种影响因子对泥石流灾害进行危害性评价和风险性预测，且这大多数研究基于已发生泥石流区域灾害现状和变化情况展开，却很少有将GIS和RS技术结合对泥石流潜在形成区进行识别和预测的研究[8-18]，而本文将针对已发生泥石流形成区的影响因素，利用RS和GIS各自可实现的技术对潜在泥石流形成区进行识别和预测。面对质量良好的国产高分影像数据，若要充分发挥其高分辨率优势实现泥石流形成区识别和预测，融合方法的选择也是重要的环节，本文应用NNDiffuse[19-20]和Gram-Schmidt[21-22]两种主流融合方法展开研究，并引入土壤亮度指数特征，利用在统计学习中最有效、最普及的SVM支持向量机监督分类[25]和常用的ISODATA动态聚类[26]非监督分类方法，并融入GIS技术实现泥石流形成区的预测，旨在能找到一种适应于当前测绘行业现状，可操作性强的应用性方案，实现潜在的泥石流形成区的最佳预测，利用各部门有效手段减少灾害带来的损失。

1 研究区概况和分析

1978年至2017年间福建省已发生泥石流灾害164次，泥石流灾害涉及全省包含安溪县、蕉城区、城厢区、闽清县、延平区等46个县市，主要分布在福建北部和东部的山区丘陵地、农村地区和偏远地区(图1)。通过对164次泥石流发生情况调查，福建省泥石流灾害149次与沟谷发育有关，且发生在集水盆地、山坡陡峭区、低植被覆盖区、裸露地表、采矿用地和道路周边，图2是已发生泥石流灾害的典型区域分布情况。根据前人的研究[25-27]并结合本研究区内泥石流发生的上述特点，分析得出对泥石流形成区具有重要影响的因素，它包括地形因素、气候因素、土壤因素、沟谷发育情况、植被覆盖和人类活动。要引入这些因素对目标区域采用传统的全局性实地考察是不切实际的[28]，布设防控措施更是难上加难，而RS技术可以实现对植被覆盖情况和人类活动影响区域识别和快速提取，GIS技术可以引入地形因素和沟谷发育情况对泥石流形成区进行预测。

图1 福建省1978—2017年已发生泥石流分布Fig.1 distribution of debris flow happened from 1978 to 2017 in Fujian Province

图2 典型泥石流发生情况分布图Fig.2 Distribution map of typical debris flow

本文研究区为福建省宁德市周宁县最南端的咸村镇，平均海拔568 m，地形总体呈西高东低，且高差悬殊，现辖23个行政村和1个社区，其中三分之一的村落聚集在东部河谷平原地区，四周环山，其余村落零星分布在偏远山区，这些地形地貌也让其面临不小的泥石流灾害风险(图3)。咸村镇的气候条件较为湿润，雨季占据全年将近三分之一，相对湿度达到80%，满足泥石流灾害爆发的条件之一。咸村镇主要以茶叶和粮食农业生产为主，是周宁县主要产粮区和产茶区，泥石流灾害发生将会给当地经济和人民生活带来巨大的损失。

图3 咸村镇地形地貌和村落分布情况Fig.3 Landform and village distribution of Xiancun Town

2 数据准备和分析

2.1 数据预处理

本文采用高分二号数据作为数据源(数据获取时间为2017年9月19日)，表1是原始数据4个多光谱波段和1个全色波段参数，覆盖研究区范围的4景影像，云量分别为12%、0%、2%和0%，图4是4景影像概视图和研究区分布情况。根据数据产品RPC函数，运用 ENVI平台完成高分二号数据预处理，单景影像依次完成辐射定标、大气校正、正射纠正、影像融合，四景影像的配准、镶嵌和研究区裁切，最后对遥感影像的质量和可行性分析。

图4 四景高分二号卫星数据情况Fig.4 GF-2 Four images satlite

表1 高分二号卫星数据参数Table 1 Data parameters of GF-2

2.2 主要融合方法介绍

为了挖掘卫星影像在泥石流形成区识别和提取中的潜力，方便定量统计和目视判读的定性研究，本文充分应用高分二号卫星影像空间分辨率的优势，应用多光谱和全色波段融合的影像展开研究。调研发现，无论是在高分一号、资源三号等国产卫星，还是Landsat 8、Worldview-2、SPOT5等国外卫星影像的融合上，对于光谱信息和局部细节保持度较好的NNDiffuse Pan Sharpening和Gram-Schmidt Pan Sharpening融合方法，在面对不同应用需求、不同传感器、同一传感器的融合表现上都有其表现效果最好的场景[19-24]，而基于这两种融合方法对泥石流形成区进行预测的研究少之又少。

(1)NNDiffuse融合算法[19-20]

NNDiffuse是由美国罗彻斯特理工学院(RIT)SUNW等人最先提出，其原理是首先将全色影像空间分辨率降采样到与多光谱影像一致，通过建立多光谱影像和全色影像线性响应向量，计算全色光谱波段的贡献向量，再运用线性混合模型得到融合影像。这种融合方法在可见光波段表现效果佳，且能很好地保持空间、光谱和色彩信息。

(2)Gram-Schmidt融合算法[21-22]

Gram-Schmidt是由Cralg A.Laben于1998年提出，该融合算法是将较低分辨率的多光谱影像转换到正交空间，再将高分辨率全色影像替换第一分量，最后利用反变换获得融合图像，这种融合算法可以克服波段数的限制，实现多波段融合，其在近红外波段表现效果佳，能有效保持原始光谱信息，且对真实地物的色彩保真度高。

图5是两种融合方法得到的影像和融合前影像中典型植被的光谱特征曲线图，融合前图像为原始光谱影像经过预处理后的影像，对比三幅影像中同一植被点的植被波谱曲线，三者光谱曲线都非常接近，融合后影像提高了图像空间分辨率但较好地保持了地物的光谱特征。图6是在1%线性拉伸的图像增强下NNDiffuse和Gram-Schmidt融合法得到的影像目视效果全局和局部的对比图，Gram-Schmidt融合影像中人工地物虽然亮度更亮，但不会因为亮度过大而削弱细节，适合实际工作中生产人员直接进行目视判读，NNDiffuse融合影像中地物更接近人眼真实看到的地物色彩，在实际工作中判读人员对地物的熟识度更高，可提高判读效率。

图5 融合前后植被光谱特征曲线Fig.5 Spectral characteristic curve of vegetation before and after image fusion

图6 NNDiffuse(左)和Gram-Schmidt(右)融合后研究区目视效果图Fig.6 NNDiffuse(Left) and Gram-Schmidt(Right) visual renderings of research area after fusion

3 泥石流潜在形成区识别和分析

3.1 识别方法介绍

RS技术可以对泥石流形成区植被覆盖情况和人类活动影响区域这两类影响因素进行一些判断和预测，具体包括无植被覆盖的自然地物和采矿用地、开垦地、居民地、道路等与人类活动密切相关地物的识别和提取。

3.1.1支持向量机法介绍

伴随信息化和现代化的不断发展，我国保险公司的财务风险管控也越来越受到社会关注，这不仅关系着国民经济的发展，也牵扯着社会关系的和谐与否。因此，建立完整的财务风险预警体系是当务之急。利用财务风险预警体系，可以在实际运行中为保险公司提供实时的风险动态，同时也能对潜在的风险进行预警通报以做好防控措施。可以有效降低财务风险发生的可能性，同时也能在风险真的发生之时以最有效的措施应对，减少因财务风险引起的企业财务损失。

SVM的原理是根据线性函数计算训练样本之间的分离超平面，SVM在高光谱影像中有不错的表现，但随着特征维数的增加，训练样本固定且获取精准的训练样本困难，Hughes 现象越发明显，分类精度也有所下降[31-34]。本文试图基于固定的小样本学习，将SVM算法引入低维多光谱分类的泥石流易发区的提取，一是满足未来应用可操作性、继承SVM的推广性、多适用性，二是探究其潜在能力。Foody曾提出，采用混合像元的小样本训练集同样可以达到和使用纯净像元的精准度，且在多光谱分类中，混合像元比纯净像元更有效[21]，故本文以特征基元的形式进行样本采集。

3.1.2基于指数特征的ISODATA动态聚类法介绍

泥石流形成的条件是拥有大量松散固体物质，大气降水、地表水等水源条件，以及陡峭的沟谷等地形地貌条件，根据福建省已发生泥石流地区周边情况可见，泥石流绝大多数发生在植被条件覆盖差或地表覆盖为裸露松散的土质及岩石区域，土壤亮度指数(NDSI)对于影像中高亮的裸地、岩石等地物的光谱特征有增强的作用，NDSI指数特征的运算公式如下：

(1)

式中：R——多光谱影像的红波段；

G——多光谱影像的绿波段。

本文利用NDSI指数特征以增强泥石流潜在形成区的信息，并利用动态聚类法(ISODATA)将有效信息进行提取，ISODATA算法[30-34]是利用合并和分裂实现地类分类的聚类分析方法，其原理是随机选择初始聚类中心，根据最小距离将样本分配给聚类中心，从而得到初始聚类，两类地类聚类中心距离与阈值之间的关系决定地类是合并为一类还是分裂为两类，在不断的分开和合并中对聚类中心的值进行更新。单纯的ISODATA动态聚类法只是基于庞大复杂的影像信息进行统计分类，而由指数特征增强后的影像进行目标提取更具针对性。

3.2 地表覆盖提取结果及分析

分别基于NNDiffuse和Gram-Schmidt融合方法的结果数据展开泥石流潜在形成区的识别和提取，研究中的训练样本和验证样本依据1 m空间分辨率影像、Google Earth和外业调查数据进行选取，保证样本数据的精确性，以便后续的分类和精度评价。根据目标地物的特征大小，分类用样本以15×15像素大小进行面状矢量勾画，训练样本总数为177个，验证样本大小不定尺寸，验证样本总数为80个。SVM监督分类中影像所包含地类分为水域、云、林地、草地、耕地、裸露地表和建筑物七类，根据本文泥石流潜在形成区提取目的，将水域、云、林地、草地和耕地归为非目标地类，裸露地表和建筑物归为目标地类。ISODATA非监督分类基于NDSI指数特征结果进行二次迭代，分类结果为5类，这是通过1～5次迭代试验和2～10次类别数设置试验后得到的理想参数值，且最终得到的类别数依旧为5类，然后根据分类的结果提取出目标地类，图7 是目标地类提取的可视化图，可见SVM分类法得到的目标地类的覆盖面积小，但图斑个数比ISODATA分类法多，表2 是对目标地类和非目标地类进行混淆矩阵精度评价情况，表中的Kappa系数是用于衡量分类结果与验证样本一致性情况的精度评价方法，总体精度是被正确分类的像元总和除以总像元数，可见(1)利用SVM 法分类精度比ISODATA方法的高，主要原因在于前者有较为准确的训练样本，且对目标地类的提取更具针对性，后者虽然分类精度低，但分类后目标地类中涵盖了前者百分之九十以上的信息，且对目标地物的边界表现得更明显，将影像中更为细节的部分进行了增强，准确的边界有利于泥石流形成区的界定。(2)Gram-Schmidt融合法在SVM监督分类中的效果稍好于NNDiffuse融合法，但NNDiffuse融合法在非监督分类中的表现要比Gram-Schmidt好，NDSI指数是针对红绿可见光波段对目标地物表现进行的增强，而ISODATA能有效将裸地、岩石这些被NDSI指数增强的信息进行聚类，但同时又抑制泥石流形成其它有效信息的提取。

图7 利用 RS技术的目标地类提取结果Fig.7 The results of target class extraction using RS technology

表2 目标地类的分类结果精度评价Table 2 Accuracy evaluation of classification results of target areas

在实际应用中，面对遥感影像的多源性和多时相性，庞大的数据量和人工处理巨大的工作量，应该充分发挥二者的优势获取有价值的潜在泥石流形成信息。

4 泥石流潜在形成区提取和分析

GIS技术可以实现地形因素和沟谷发育情况这两种影响因素对泥石流形成区的预测，且可以将前期RS技术可实现的自然地表覆盖和人类活动这两类因素融入到最终的判定和预测中。

利用ArcHydro水资源工具实现区域沟谷线、集水区、汇流点等要素属性信息的挂接和汇流分析，综合ArcGIS应用系统的空间分析和3D分析模块解决上一步的地表覆盖矢量图斑、集水区、汇流点的空间拓扑关系及属性查询筛选有效信息，实现泥石流形成区的最终预测，预测流程如图8所示。

图8 泥石流潜在形成区预测主要技术流程图Fig.8 Main technical flow chart for prediction of potential formationregion of debris flow

利用原始2 m DEM高程数据计算水流方向并填充洼地，填洼后重新计算水流方向，提取洼地，直到研究区无洼地，计算河流流量和沿水流方向到下游出水口的水流长度并生成河流网络，利用ArcHydro工具获取记录河网信息的河流链接数据，并生成包含河流属性信息的研究区所有集水区、沟谷线、汇流点和出水口的矢量数据。根据《T/CAGHP 006—2018泥石流灾害防治工程勘察规范(试行)》规范中极易发、中等易发、轻度易发和不易发生四种易发程度和坡度、流域相对高差和流域面积三种参数的量化评分表，并结合福建省已发生泥石流数据得到本文研究的泥石流易发生的三种参数量化值，即规范内轻度易发、中等易发和轻度易发统一归纳为研究区易发生泥石流灾害，具体参数见表3，筛选内容主要包括地表覆盖图斑120 m范围内的集水区，通过沟谷线平均坡度、相对高差和集水区的流域面积的属性和空间关系筛选出易发生泥石流的沟谷线，通过得到的沟谷线获得对应的集水区，最后根据地表覆盖与地形流水参数的空间关系筛选后得到筛选后地表覆盖，筛选后的隐患沟、集水区即是潜在的泥石流形成区，可视化预测情况见图9，图中的地表覆盖、集水区、隐患沟重叠部分比例大，当地表覆盖面积越大，最终的泥石流潜在形成区的涉及范围也越广，这可以减少在防治措施布设中的隐患区遗漏，而与泥石流形成区相关的地表覆盖面积越小，又可以减轻在布设经费上的支出，减轻财政负担。

表3 不同易发程度的泥石流灾害规范性量化参数和本文参数情况Table 3 Standardized quantitative parameters of debris flow disasters with different susceptibility and parameters in this paper

图9 四种方法和综合四种影响因素的泥石流形成区预测情况Fig.9 Prediction of debris flow formationregion by four methods and four factors

表4是四种方法对研究区泥石流潜在形成区预测结果中所有隐患沟的沟长、坡度和相对高差的平均值统计，以及集水区汇流面积平均值统计的数值分布情况，根据《T/CAGHP 006—2018泥石流灾害防治工程勘察规范(试行)》[35]数量化评分表划分已提取的泥石流潜在形成区的级别为极易发生泥石流区，这些数值可为未来的泥石流灾害预测工作提供经验性参考。

表4 研究区泥石流灾害隐患沟和集水区四种参数量化情况Table 4 Parameter quantification in potential region and water collection area

在泥石流形成区预测中集水区的筛选分为第一次筛选和第二次筛选，第一次筛选具体指研究区所有 4 105 个集水区经过与地表覆盖图斑空间位置关系的筛选，第二次筛选是引入沟谷、集水区、汇流点和出水口后通过空间拓扑和属性关系的筛选；与潜在泥石流形成区相关的地表覆盖只有一次筛选，筛选前即应用RS技术分类提取的结果，筛选后即此分类提取的结果与集水区最终筛选结果进行的空间关系筛选，具体的筛选情况见表5，集水区第一次筛选率即集水区剔除的个数占总集水区个数比例，第二次筛选率即将不符合表3中参数的集水区剔除的个数占总集水区个数的比例，集水区正确个数是通过GoogleEarth Pro的3D图、DEM和高分辨率影像融合叠加的目视判读，及对目视无法确认的区域进行野外核实的结果。结果可见，Gram-Schmidt-SVM第一次筛选率最高，其次是Gram-Schmidt-ISODATA，在第二次筛选中NNDiffuse-SVM筛选率最高，其次是 NNDiffuse-ISODATA，但在最终提取准确率上NNDiffuse两种方法比基于Gram-Schmidt两种分类方法提取准确率高，可见，与流水有关的更全面的地形信息的引入有助于提高泥石流潜在形成区的正确提取；NNDiffuse融合方法比Gram-Schmidt融合方法对提取潜在形成区更有效；ISODATA方法比SVM方法更有效；地表覆盖提取精度高并不一定最终泥石流提取的效果就更好，泥石流潜在形成区提取还得依靠地表覆盖与沟谷线、集水区之间的空间和属性关系确定，因为地形要素是泥石流发生的重要因素之一。

表5 基于空间和属性关系筛选的集水区和地表覆盖情况Table 5 Catchment area and surface coverage based on spatial and attribute relation screening

综上可见，NNDiffuse-SVM方法在与影响泥石流形成的相关地表覆盖提取和泥石流形成区预测中效果总体最好，NNDiffuse-ISODATA方法在无先验知识情况下的提取效果总体最好，Gram-Schmidt-SVM的效果总体比Gram-Schmidt-ISODATA稍好，在实际应用中，生产人员专业水平的不同也会影响需要先验知识的SVM方法预测形成区的效果，综合专家知识、提取效率和适用性来看，建议采取ISODATA的方法进行泥石流潜在形成区的提取。从遥感数据快速提取影响泥石流形成的相关地表覆盖的方法是有效的，但更准确的泥石流形成区提取应该综合考虑植被覆盖情况、人类活动、地形因素和沟谷发育情况才能获得比较好的效果。基于当前测绘行业现状，既要发挥人工操作的灵活性、稳定性和准确性，又要充分发挥自动化提取的省时省力优势，泥石流形成区预测应该在发挥主观能动性基础上应用机器自动化，实现两种方法的优势互补。

5 结论与讨论

本文基于专家知识和机器自动化相结合的方式，实现对泥石流潜在形成区的半自动化提取，引入RS应用技术，实现与影响泥石流形成的地表覆盖快速提取，同时引入GIS应用技术实现空间信息、拓扑信息和属性信息查询和筛选，实现泥石流形成区的多角度综合提取。针对研究结果和研究过程中存在的问题，本文提出以下总结和展望：

(1)遥感影像的融合方法选择会影响泥石流形成区预测结果，基于NNDiffuse融合方法的泥石流形成区提取结果比Gram-Schmidt方法好。本文主要针对高分二号多光谱波段和全色波段融合影像进行研究，在对其它卫星影像泥石流形成区预测是否表现一致具有深入研究的必要。

(2)SVM的监督分类法在提取泥石流潜在形成区中效果最佳，而无需先验知识的ISODATA方法在泥石流形成区预测中有较好的表现，能基本满足当前测绘行业的应用需求。本文主要采用适合当下测绘行业推广性和适用性强的监督分类和非监督分类方法进行研究，分类方法的选择是否会影响提取效果，尤其是非监督分类对泥石流潜在形成区提取的潜力有待进一步挖掘。

(3)利用遥感应用技术提取地表覆盖和引入地理信息应用技术对沟谷信息筛选的应用方法有较大潜力实现泥石流潜在形成区提取，同时，在后续研究中将引入泥石流的流通区、堆积区等情况进行综合考量。