滑坡易发性研究可较准确地预测特定区域内潜在滑坡发生的空间概率. 目前滑坡易发性预测过程中广泛使用的机器学习被认为具有比数理统计模型更好的非线性预测能力^[1-2]. 按照是否利用已知样本数据作为模型输出变量，可将机器学习模型分为无监督和全监督两大类.

无监督机器学习可直接利用控制因素作为输入变量进行滑坡易发性建模，其主要包括K均值聚类^[3]、SOM神经网络^[4]. 全监督机器学习须同时具备模型输入和输出变量，并从给定的训练数据中学习出非线性函数，再根据该非线性函数计算输入样本的预测值. 由于建模过程充分利用已知滑坡样本的先验信息，全监督机器学习能预测出更准确的滑坡易发性^[5]. 常用的全监督机器学习包括BP神经网络 (back-propagation neural network, BPNN)、支持向量机^[6-7]、卡方自交互侦测决策树(Chi-squared automatic interaction detector, CHAID)^[8]、逻辑回归^[9]、随机森林^[10]等. CHAID模型是一种简单易用的非参数分类器，它不需要对数据有任何的先验假设就能得到较高的预测精度且其计算效率较高^[11]. BPNN是一种成熟的机器学习模型，能精确拟合出各种复杂非线性问题，被较早地用于滑坡易发性预测领域^[12].

虽然全监督和无监督机器学习预测滑坡易发性已取得一系列成果，但仍存在一些不足. 一方面，尽管无监督机器学习在训练和测试过程中不需要已知滑坡和非滑坡样本作为模型输出变量，但是缺乏滑坡和非滑坡先验知识的引导，无监督机器学习的建模准确性难以保证^[13]. 另一方面，基于全监督机器学习的滑坡易发性预测建模，存在野外获取滑坡样本数据难度较大且不完整，在整个研究区随机选取非滑坡样本造成大量误差的情况. 为了解决全监督和无监督机器学习存在的已知滑坡编录不完整、随机选取的非滑坡不准确的问题，本文拟利用半监督机器学习^[14]开展易发性预测建模. 半监督机器学习是一种能综合利用无监督和全监督机器学习优点的算法，其可在标记样本较少的情况下，充分利用大量的无标记样本提高学习性能. 自训练算法具有结构简单、预测精度高、应用范围广等优点，是典型的半监督机器学习模式. 本文拟利用自训练算法的半监督学习模式将传统机器学习从全监督模式扩展成半监督模式，实现半监督机器学习模型构建.

半监督机器学习在扩充已知滑坡编录数量时拟借鉴张为等^[15]提出的潜在滑坡识别思路，即解译的直接目标是发生滑坡概率高且灾害损失较大的特定地点. 研究思路如下：1）将机器学习预测的初始极高滑坡易发区与高分遥感影像叠加，达到大幅度缩小遥感解译范围和提高遥感解译针对性的目的。2）通过高分遥感或目视解译研究区极高概率的潜在滑坡. 3）从初始极低滑坡易发区选取非滑坡样本，确保认定为非滑坡的区域发生滑坡的概率极低. 4）有效扩充后的滑坡样本、准确的非滑坡样本再次导入全监督机器学习进行易发性预测，可实现更加准确可靠的滑坡易发性制图.

本文以江西省南康区为研究对象，基于全监督机器学习的CHAID和BPNN模型，拟采用自训练算法构建半监督CHAID (semi-supervised CHAID, SSCHAID)和半监督BPNN (semi-supervised BPNN, SSBPNN)预测滑坡易发性并做对比分析，探讨半监督机器学习的建模效果和预测精度.

本文提出的半监督机器学习^[15]以解决易发性预测中存在的已知滑坡编录较少和非滑坡样本随机选取的缺点. 如图1所示，基于自训练算法的半监督机器学习主要体现为全监督机器学习的二次迭代. 1）首先通过RS和ArcGIS平台管理并空间分析研究区内的滑坡编录及相关控制因素；2）基于频率比(frequency ratio, FR)和相关性分析，确定滑坡与其控制因素间的联动关系；3）基于控制因素的FR、已知滑坡栅格单元和随机选择的非滑坡栅格单元，对全监督机器学习进行训练测试并预测出初始滑坡易发性值；4）将高分辨率遥感影像与初始滑坡易发性图叠加，从初始极高滑坡易发区中通过遥感和目视解译的方式选择发生滑坡概率极高的栅格单元作为潜在滑坡，以扩充已知滑坡样本；5）从极低易发区中随机选取栅格单元作为非滑坡样本；6）将扩充后的滑坡样本和准确选取的非滑坡样本，再次导入同一机器学习中进行训练测试，成功构建半监督机器学习并进行最终的滑坡易发性预测；7）对上述监督和半监督机器学习的滑坡易发性预测结果及其精度进行对比分析.

CHAID是决策树中的回归预测方法，其特点是可分析大量的滑坡及滑坡控制因素，并对控制因素进行最优分类. CHAID建模步骤如下. 1）数据收集：离散或连续型数据均可作为模型输入变量，能接受缺失数据且无须标准化. 2）数据分析：每个分支代表一个判断结果的输出，每个叶节点代表一种分类结果，可以检查分类结果是否符合预期. 3）训练算法：在给出一批样本的情况下，每个样本都有自己的属性和分类结果，通过样本的全监督机器学习得到决策树以分类新的数据. 4）测试算法：使用经验数计算错误率，当错误率达到可接受的范围时表示CHAID模型可以被使用. CHAID中只要因变量与控制因素间没有显著卡方值则分类迭代的性能就会停止.

BPNN由输入层、隐含层和输出层组成，当输入节点为 $ n $且输出节点数为 $ m $时，BPNN表达从 $ n $个自变量到 $ m $个因变量的函数映射，对网络进行训练直到达到网络期望和实际输出值之间的最小目标误差. 在BPNN中通过大量的自适应学习可产生高度非线性的输入、输出对应关系. 输出节点表达式为

式中： $ {w_{ij}} $为输入层中节点 $ i $与隐含层中节点 $ j $的连接权值. $ {w_{ij}} $随学习进度变化的表达式为

式中： $ z $为迭代次数， $ {y_i} $为输入节点， $ {\delta _j} $为隐含层中节点 $ j $的局部梯度， $ \varphi $为学习率参数， $ \;\beta $为动量常数。本文使用的输出函数是S形函数.

根据研究区的基础地质资料，基于ArcGIS平台与遥感影像目视解译获取南康区的地形地貌、水文环境、地层岩性、地表覆被等控制因素. 基于控制因素的频率比、已知滑坡编录和随机选择的非滑坡栅格单元，预测研究区的初始滑坡易发性，并进一步提取南康区的极高滑坡易发区；通过对研究区历史滑坡形态和色调等特征的分析，建立区域滑坡遥感解译标志；根据滑坡遥感解译标志，并基于人机交互式目视，解译极高易发区内的滑坡隐患点. 具体解译流程如图2所示.

滑坡遥感解译标志分为直接标志和间接标志. 其中直接标志主要表现为滑坡平面形态多为圈椅形、舌形，剖面形态多为凹形、凸形、阶梯形等；滑坡整体呈浅色调，沿坡体两侧常形成沟谷，自然沟切割较深，有时会出现双沟同源现象；坡底下方由于滑体的挤压，影像有时会显现高低不平的地貌；滑体滑移停止以后滑坡裂缝逐渐发育成冲沟，在遥感影像上显现为带状影纹和明显的色调差异；滑坡舌延伸到平缓斜坡或河道，遥感影像上显现为较自然地面略高的舌状影纹. 间接标志主要表现为山坡沟谷出现沟槽突然改道、横断面显著变窄变浅；滑坡体上植被多为马刀树、醉汉林；不正常河流弯道、局部河道突然变窄，形成堰塞湖.

南康区位于江西省南部，地处北纬25°28′00″ ~ 26°14′24″，东经114°29′9″~144°55′24″. 全区海拔介于96~995 m，南北长约85.45 km，东西宽约42.6 km，总面积约1 844.96 km². 南康区地处南岭山脉东端北坡，地势西高东低、南北高、中部低. 主要地貌类型为丘陵、山地，沿章江、上犹江两岸有较广阔的河谷平原. 境内水资源丰富，属中亚热带季风性湿润气候，年降雨量约1 120 mm，降雨主要集中在4~6月的春雨、梅雨期和7~9月的台风雨期.

据统计，南康区1970年至2010年已发生或存在隐患的地灾点约338处，其中滑坡233处，崩塌105处. 滑坡体以第四纪堆积层为主，运动方式主要是牵引式整体滑动，单个滑坡几何形态一般呈圈椅形. 南康区滑坡按规模划分以中小型为主，滑坡平均面积约为1.0×10⁴ m²，堆积层厚度约为6 m；南康区稳定性较差的滑坡约占总滑坡的60%。滑坡的主要诱因是暴雨，其次是大规模工程建设. 全区地灾已造成上百人伤亡和几千万元的经济损失.

本文研究主要数据源包括：1) 南康自然资源部门调查的滑坡资料；2) 30 m分辨率DEM 数据，用于提取地形地貌及水文因素；3) 利用江西省地质调查局测绘的1∶10万比例尺的地质图提取的地质信息；4) 30 m分辨率的Landsat TM8遥感影像1景 (2013年7月3日，轨道号119/041)，用于提取地表覆被因素；5) 从91卫图软件中下载的1.07 m分辨率航片 (2018年10月5日)用于识别潜在极高概率的滑坡. 由于30 m栅格既能有效表征研究区地形特征又能满足模型计算要求^[16]，本文的 DEM 和遥感影像原始分辨率均为30 m。

本研究依据江西相关县域滑坡易发性预测的文献资料、滑坡与其控制因素间关联值的计算结果，考虑相关控制因素获取的难易程度，从数据源中获取如表1所示的地形、地质、水文和地表覆被等4大类共11个控制因素^[17]. FR法是高效的定量分析法，可揭示各控制因素对滑坡发育的影响程度^[18].

式中： $ {l_i} $为某个控制因素的第 $ i $个属性区间内的滑坡面积， $ L $为研究区内滑坡的总面积， $ {s_i} $为控制因素第 $ i $个属性区间的面积； $ S $为研究区总面积.

由式（3）可知FR平均值为1. 由表1可知，FR>1表示控制因素某个区间内的属性有利于滑坡发育，FR<1表示该属性不利于滑坡发育. 在SPSS 23软件中通过相关性分析计算控制因素间的共线性，表明各控制因素的相关性不大，均可用作模型输入变量.

如图3所示，高程、坡度、坡向、剖面曲率、平面曲率和地形起伏度等地形因素通过GIS空间分析从DEM中获取. 南康区高程在163.450~360.967 m的FR>1，坡度在6.271°~21.211°时FR>1，表明滑坡主要发育在中等坡度的边坡. 变质岩、碎屑岩的FR分别为1.414、1.176.

通过获取南康区栅格单元距离水系距离和改进的归一化差异水体指数(modified normalized difference water index, MNDWI)表征水文环境对滑坡发育的影响. 如表1所示，当距离水系距离小于300 m时，FR>1. 本文主要采用归一化建筑物指数(normalized difference building index, NDBI)和归一化植被指数(normalized differential vegetation index, NDVI)表征地表覆被因素. 当NDVI为0.56~0.75时，FR>1；NDVI>0.75时，FI $\ll $1。说明一定强度的工程建设容易诱发滑坡，高植被覆盖区一般不容易发生滑坡.

将所选11个控制因素的FR归一至[0,1.0]，作为CHAID的输入变量和另外3个模型的输入变量. 在ArcGIS中将已发生的233处滑坡的面文件转为栅格单元，得到总计2 598个滑坡栅格单元. 将2 598个已标记的滑坡栅格和在南康非滑坡区中随机选择的2 598个非滑坡栅格单元组成模型训练测试数据集。将数据集随机划分为2个部分：70%的数据用于训练，其余的用于测试. 在CHAID模型训练测试过程中将已知滑坡栅格和随机选的非滑坡栅格的易发性标签分别设定为1和0. CHAID的输出变量反映滑坡易发性的各栅格单元位于[0, 1.0]的概率.

由交叉验证法得到CHAID模型主要参数：停止条件最大层数为5层，停止规则选择父节点、子节点的最小样本百分比分别为2%、1%，分割节点的显著性水平为0.05且收敛的最大迭代数为100次. 用训练测试后的CHAID预测得到整个南康区的滑坡易发性，采用自然间段点法^[8]并结合滑坡易发性分布规律，将南康区划分为5类滑坡易发性级别：极高 (9.9%)、高 (10.2%)、中等 (19.1%)、低 (36.9%)、极低 (23.9%)如图4(a)所示. 从极低到极高易发区的FR分别为0.024、0.286、0.953、2.218、4.851.

SSCHAID的建模步骤如下。1）在极高滑坡易发区通过高分遥感影像解译确定520个极高概率的潜在滑坡栅格单元，这些栅格占已知滑坡栅格的20%，可有效提高滑坡栅格数量，避免滑坡栅格数量太多导致可信度下降；2）用这520个潜在滑坡栅格扩充已知的2598个滑坡栅格，组成3 118个已知滑坡栅格并设标签设为1；3）在极低易发区中随机选取3 118个极高概率的非滑坡栅格并设标签为0；4）将扩充后的滑坡和非滑坡栅格数据再次用于CHAID的训练测试.

将第2次训练测试后的CHAID模型用于南康区滑坡易发性预测，同样采用交叉验证法确定CHAID的参数，结果与第1次训练的参数差别不大. 为了方便模型对比，SSCHAID模型也将预测的滑坡易发性分为5个等级：极高(15.1%)、高(15.6%)、中等(14.8%)、低(16.3%)、极低(38.2%)如图4(b)所示. 从极低到极高滑坡易发区的FR逐渐增大，高、极高易发区的FR分别为1.38、4.27，极低、低易发区的FR分别为0.047、0.288，可见SSCHAID模型得到的极高、高滑坡易发区与已知滑坡分布特征非常吻合且其整体预测精度非常高.

再以所述11个控制因素的FR作为BPNN的输入变量. 本文研究发现随机划分出的多组不同训练测试集构建的BPNN模型预测效果差别很小，因此最终仅选择某次训练集和测试集划分结果开展BPNN建模^[8]. 使用R语言程序选择前述样本数据（2598个已知滑坡栅格和随机选取2598个非滑坡栅格）作为训练测试样本，载入NNET函数包建立单隐藏层的BPNN，经过循环迭代计算得到BPNN隐藏层节点数为12个，初始随机权重的范围[−0.1, 0.1]；当模型权重的衰减精度小于 $ 12\times10^{-2} $时不再迭代，最大迭代次数为5000次. 将得到的权重应用于BPNN训练测试并预测整个研究区的滑坡易发性，得到的滑坡易发性图如图4(c)所示.SSBPNN模型的建模步骤与SSCHAID模型一致. 最终的滑坡易发性图如图4(d)所示.

如图5所示，采用受试者工作特征曲线(receiver operation characteristic curves, ROC)下的面积 (area under ROC, AUC)^[8]分别评价4个模型预测精度. 图中，TPR为真阳率，FPR为假阳率. 由图可知，SSCHAID、CHAID、SSBPNN、BPNN的AUC分别为0.926、0.821、0.967、0.872，表明半监督机器学习大幅度提高了全监督机器学习的易发性预测精度. 进一步表明，通过扩充已知滑坡样本和准确有效地筛选出非滑坡样本，能大幅度提高机器学习的易发性预测性能. 由图还可看出，BPNN对于南康区的滑坡易发性预测精度略优于CHAID模型.

通过计算各模型测试样本与其易发性预测值的混淆矩阵可得到模型敏感度S(sensitivity)、真预测率(positive prediction rate, PPR)、总准确率(total accuracy, TA)共3个性能评价指标^[19]. S为在实际为滑坡的样本中被正确判断为滑坡的比例. PPR为在预测滑坡结果中，被正确判断为滑坡的比例. TA为正确分类的样本数与总样本数之比. 4个模型的统计学精度如表2所示，由表可知，SSCHAID、SSBPNN的S均显著高于BPNN、CHAID模型的，表明半监督机器学习预测滑坡易发性的准确率更高，且全监督机器学习模型存在的缺点得到大幅度改进.

SSCHAID、SSBPNN的易发性预测精度在CHAID、BPNN模型基础上大幅度提升. 原因是全监督机器学习随机选取非滑坡栅格单元作为模型训练和测试用的输出变量，这些不确定的非滑坡样本导致模型训练和测试过程存在大量误差，降低了全监督机器学习建模精度. 半监督机器学习在建模过程中从极低易发区选择可信度非常高的非滑坡样本，训练和测试数据集的误差减少、建模精度提高；通过筛选极高概率的潜在滑坡扩充已知滑坡样本，可使得半监督机器学习的训练测试样本具有更广泛的代表性，使得训练后的模型能更准确地反映滑坡与控制因素间的非线性函数关系. 综上分析可知半监督机器学习很好地利用并扩充现有的滑坡样本来指导建模过程，在全监督机器学习的基础上进一步提升滑坡易发性预测建模性能.

CHAID和BPNN均属于广泛使用的传统机器学习. CHAID能够从毫无规律的数据样本集合中推理出具体的表示形式及分类规则. BPNN通过自身的训练学习某种规则，在给定输入值时得到最接近期望的输出值，实现从输入到输出的非线性映射. 虽然CHAID建模过程比BPNN简单，但南康区易发性预测结果显示，BPNN的预测精度优于CHAID. 原因是BPNN模型可以更好地从复杂的滑坡及其控制因素大数据中表示滑坡易发性的属性特征，且具有更强的非线性泛化性能. 此外，CHAID模型也存在缺点，比如1）该模型可以创建复杂的树结构但是缺乏推广依据，可能导致过拟合现象；2）该模型是基于启发式算法寻求在每个节点上的局部最优决策返回，这种思路有时难以实现全局最优；3）该模型很可能在某些类占主导地位时创建出有偏异的树，导致建模性能受影响.

本文利用ARCGIS软件的空间分析功能，将初始极高滑坡易发区与高分遥感影像相叠加得到范围非常小的极高滑坡易发区，再从中依据滑坡影像和形态特征随机选择一部分斜坡单元作为新的滑坡样本^[20]. 在识别潜在滑坡时尽可能选择靠近农村居民区或公路网络且植被覆盖率较低的斜坡。原因是这类斜坡发生灾害的风险往往较大，这种方案可提高滑坡易发性预测结果的可靠度和实用性. 同时建议识别的潜在滑坡要广泛分布在整个研究区，以提高半监督机器学习建模的通用性，避免出现人为因素对建模的过度干扰.

这些扩充的滑坡样本不能保证是真实的滑坡，它们是高概率的滑坡隐患区. 滑坡样本扩充目的之一是增加模型训练测试样本，二是为了提高滑坡样本在研究区空间分布的均衡性. 为了避免降低已知滑坡编录样本的权重，通过反复试验本文将潜在滑坡样本扩充的比例确定为已知滑坡样本的20%. 在下一步研究中可定量考虑不同研究区、不同扩充比例下的半监督机器学习建模特征.

一部分初始滑坡易发性区域之所以被划分为极低易发区，就是因为这些区域在自身控制因素的非线性综合作用下发生滑坡的概率很低. 因此直接从极低易发区随机选择非滑坡栅格能有效保证这些非滑坡栅格的滑坡易发性是非常低的. 此外，不应该单独考虑某个滑坡控制因素(如坡度、水系、地层岩性等)选择非滑坡样本，这样选择出来的非滑坡样本难以能保证边坡自身发生滑坡易发性非常低.

从极低易发区随机选取非滑坡栅格作为负样本时，如果已知滑坡样本分布规律非常明确，且全监督机器学习受随机非滑坡样本的影响较小，就有可能导致半监督机器学习预测的滑坡易发性向极高和极低易发区集中分布。该现象与机器学习本身的过拟合也有一定关联. 若从极低易发区选择非滑坡样本出现该问题，建议同时从极低和低易发区选择非滑坡样本，避免极高和极低易发区比例过高及全监督机器学习出现过拟合.

本文缺少与半监督学习基准算法的对比，这些不同半监督机器学习模型导致的滑坡易发性预测不确定性问题将在下一步的研究中重点探讨. 为了进一步厘清半监督机器学习模型预测滑坡易发性的思路并实现建模自动化，下一步研究中可增加算法伪代码并编程实现本文的建模.

从4类模型预测的易发性结果可知南康区滑坡易发性分布规律与其控制因素相吻合. 1）南康区极高和高易发区主要分布在高程163~361 m、坡度9°~22°、距离水系距离小于250 m、变质岩和碎屑岩以及工程建设较频繁的区域. 较低的高程和中等坡度有利于边坡堆积层形成；距离水系较近的区域地下水渗流显著，且容易受水体的冲刷侵蚀作用；变质岩、碎屑岩等工程地质条件差、易风化等特性使坡体力学平衡容易被破坏；大量工程建设的区域边坡受工程破坏的概率更大. 上述自然环境的综合作用区域是南康区滑坡防治工作的重点区. 2）南康区极低和低滑坡易发区主要分布在海拔高于360 m、坡度非常小或者较大及工程建设较少的区域. 海拔较高的地区远离河流冲刷、蓄水条件较差、人类工程活动少且地表植被覆盖密集不易发生滑坡.

(1) SSCHAID、SSBPNN的易发性预测精度显著高于CHAID、BPNN模型. 半监督机器学习在滑坡易发性的滑坡易发性预测精度显著高于全监督机器学习模型的.

(2) 半监督机器学习为极高概率的潜在滑坡样本扩充提供新思路，并很好地解决了随机选取非滑坡准确性很低的问题.

(3) 本文尚未研究空间数据集的滑坡−非滑坡样本非对称特征问题，现有研究大部分采用1∶1的对称策略选择滑坡−非滑坡，与实际野外工程情况不相符. 下一步将在半监督机器学习的基础上考虑滑坡−非滑坡样本非对称特征问题，进行更符合实际野外工程情况的滑坡易发性建模，获取更合理的易发性图.