基于半监督学习模型的地质环境承载力评价: 以府谷县为例

李宇新, 邓念东, 郭鹏程, 张富荣, 郭亚雷

(1.西安煤航遥感信息有限公司, 西安 710199; 2.中国煤炭地质总局航测遥感局, 西安 710199;3.西安科技大学地质与环境学院, 西安 710054)

地质环境系统是人类生产、生活赖以生存的基础载体。随着对地质环境研究的不断加深,人们认识到地质环境系统抵抗人为破坏具有一定的承载限度,所能承受人类改造作用的极限即为“地质环境承载力”[1]。

21世纪初,随着可持续发展理念深入人心,地质环境承载力内容与内涵不断延伸,逐渐演变为“人口-资源-环境-经济-社会”为一体的综合研究。彭立等[2]分析了汶川地震重灾区10个县域资源环境的变化趋势,确定了与资源环境承载力相适应的人口空间布局。姚治华等[3]从资源、环境、调节和社会经济4个方面建立评价指标体系,采用基于熵值—层级分析法的集对分析模型对大庆市地质环境承载力进行评价,对承载状态变化趋势提出改善策略。倪忠云[4]以西藏曲水-桑日地区为研究对象,基于生态地质环境空间分布和时间尺度两个角度,对气象、植被、土地利用等进行单因素分析,最终得出该地区生态地质环境承载力等级区划。李得立等[5]引入基于决策树的机器学习算法,对眉山市彭山区进行区域地质环境承载力评价,实现了评价方法的高效化和智能化。王念秦等[6]对骊山镇1980—2015年地质环境承载力进行动态评价,进一步分析总结地质环境承载力变化规律与趋势。谢正团[7]针对白龙江流域地质灾害高发区,选取地质灾害风险性和生态环境承载两个角度,采用模糊逻辑理论和层次分析法进行了地质环境承载力评价,创新性地利用时序干涉合成孔径雷达(Interferometric synthetic aperture radar, InSAR)的地表变形对评价结果验证。孙魁等[8]将保水采煤理念与地质环境承载力研究有机结合,提出了基于保水采煤的地质环境承载力的内涵与评价思路,为生态脆弱区煤炭开采和水资源保护提供理论支持。李仁伟[9]分别从滑坡易发性、水土质量评价和地质遗迹保护3个方面,基于GIS采用综合评价模型分析了神木市地质环境承载力。支泽民等[10]通过随机森林模型确定权重,对川藏铁路沿线县域开展地质环境承载力评价。

2019年9月,黄河流域生态环境保护与高质量发展上升为重大国家战略。府谷县地处黄河流域中上游、黄土高原为主的水土保持区,生态地质环境脆弱,矿产资源丰富。近20年来,随着府谷县城镇基础设施建设,能源化工基地建设以及矿产资源的深度开采,在带动社会经济高速发展的同时,致使原本脆弱的地质环境进一步恶化,滑坡、崩塌、泥石流、地面塌陷等地质灾害问题突出,诱发植被景观破坏、水土流失、土壤盐碱化等一系列环境问题[11]。如今,这些地质环境问题已经严重制约黄河流域生态环境保护和城镇化高质量发展,威胁当地人民生活和生产的健康安全[12],亟需开展区域性地质环境承载力专题与综合评价,服务于国土空间开发适宜性评价,以指导城镇防灾减灾和生态安全建设[13]。

其中,机器学习方法凭借高效、智能的数据处理能力,广受地质环境评价等空间预测领域研究学者的青睐。机器学习监督学习方法主要通过已知标签样本进行学习训练,从而达到较高准确率的预测效果。而地质环境承载力评价的目标策略是求解未知标签的学习过程,因此引入半监督学习模型,通过无监督学习和有监督学习模型组合形式进行地质环境承载力评价,并进一步分析影响因素与承载力潜在联系,为黄土高原生态环境脆弱区地质环境承载力评价提供参考。

1 研究方法

1.1 主成分分析

主成分分析(principal component analysis,PCA)是常见的降维算法之一,在数据噪音和冗余等领域具有广泛的应用,其本质是提炼原始数据中部分重要成分,通过特征值分解,有利于计算和反映原始数据的空间关系特征[14]。主要原理具体如下。

(1)构建评价指标相关系数矩阵R。

(1)

(2)

(2)计算方差百分比Pi与累积方差百分比Ai。

|λi-R|=0

(3)

(4)

(5)

(6)

式中:λi为R矩阵的特征值;uij为特征向量的分量;m为主成分数。

(3)计算主成分载荷Zij与各因子权重T。

(7)

(8)

式(8)中:Sm为各因子主成分载荷与方差百分比乘积和。

(9)

1.2 反向传播神经网络

反向传播(back propagation,BP)神经网络是人工神经网络中的一种类型,该方法主要解决多层神经网络的学习问题,因此极大地促进神经网络的发展,被广泛应用于分类识别、回归预测等领域[15]。它属于有监督学习模型,基本思想为梯度下降法,主要可分为信息的前向传播和误差的反向传播两个步骤。BP神经网络由输入层、隐含层和输出层构成,层与层之间进行全连接,而同一层之间的神经元节点无连接(图1),其基本原理如下。

隐含层第i个神经元的输入和输出可分别表示为

(10)

(11)

输出层第j个神经元的输入和输出分别表示为

(12)

(13)

式(13)中:g(·)为线性函数。

输出层每个神经元的误差和总误差分别表示为

(14)

(15)

(2)接着,根据最速下降法计算误差为ωij的梯度[式(16)]与权值修正值[式(17)]

(16)

(17)

式中:g′(·)为修正后线性函数;Δωij为权值修正值;η为修正系数。

(3)同理,根据误差信号反向传播来调整输入层与隐含层之间的权值ωmi。

(18)

Δω=ηδv

(19)

2 研究区概况

府谷县位于陕西省最北部,地处秦、晋、蒙三省交界。地理坐标介于38°42′28″N～39°33′44″N,110°25′40″E～111°15′36″E,总面积约为3 229 km2。研究区属于中温带干旱大陆性季风气候,温差较大,四季分明,气候多变。降雨在空间分布上由西北向东南逐渐增大。月降雨量分布不均,主要集中在7—9月,占全年降雨量的69%。研究区水系发育,河流交错,主要有黄甫川、孤山川、清水川等,均为黄河支流。

研究区整体位于向西倾斜的单斜构造,新构造运动主要表现为震荡性上升,形成西北高、东南低的地势。西部地势相对平坦,主要为风蚀作用形成的剥蚀地貌;而东部沟壑纵横,以河流切割的侵蚀地貌为主。出露地层由老到新为奥陶系、石炭系、二叠系、三叠系、侏罗系、新近系和第四系。按照含水介质与富水条件不同,水文地质类型划分为松散岩类孔隙水、基岩裂隙水和碳酸盐岩类岩溶裂隙水三类。植被类型主要划分为灌丛类、草原类、针叶林和栽培植被。

研究区人类工程活动剧烈,主要表现为矿产资源开发利用、道路工程建设、城镇化建设等,对原本脆弱的生态环境加剧了扰动。根据现场调查,研究区共有地质灾害185处,其中崩塌140处、滑坡29处、地面塌陷8处以及泥石流8处,各行政区划地灾数量与密度统计如图2所示。

图2 各乡镇地质灾害密度统计Fig.2 Geological hazard density statistics of each town

3 评价体系的建立

3.1 评价单元的选取

评价单元指区域地质环境评价中信息构成的基础单元,它具有内部属性一致、不可分割的特点[16]。由于地质环境承载力评价尚属探索挖掘阶段,关于该领域的评价单元未进行详细对比与确定。然而,关于地质环境评价的相关研究值得借鉴,尤其是地质灾害易发性、风险性评价,常见的评价单元包括行政单元、栅格单元、斜坡单元和流域单元[17]。

斜坡单元是流域单元更细小的划分,而斜坡单元又由栅格基础单元组成。由于斜坡单元、流域单元具有操作繁琐、更适用于地形分析的局限性,而地质环境承载力内容涉及地形地貌、社会经济等多个方面,数据类型复杂,属性信息量庞大,因此选择栅格单元作为评价基础单元。结合数据源地质图1∶50 000比例尺,依据经验公式计算对应栅格单元大小,计算公式为

G=7.49+0.000 6S-2.01×10-9S2+

2.91×10-15S3

(20)

式(20)中:S为比例尺分母,即50 000;G约为30 m。

最终选择30 m×30 m栅格单元作为地质环境承载力评价单元。

3.2 评价指标预处理

结合研究区地质环境背景与地质环境承载力内涵范畴,从自然地理、生态环境、矿产资源、灾害地质和社会经济5个子系统共选取17类影响因素(图3),部分专题图如图4所示。

图3 地质环境承载力评价体系Fig.3 Evaluation system of geological environment carrying capacity

图4 主要影响因子专题图Fig.4 Thematic map of major impact factors

图4中,①地表切割深度、坡度、曲率因子数据通过30 m分辨率数字高程模型(digital elevation model, DEM)分析提取;②通过矢量化研究区1∶50 000地质图得到地形地貌、地质构造和地层岩性类因子图层;③通过Landsat 8 OLI_TIRS卫星数字产品提取归一化植被指数(normalized difference vegetation index,NDVI)数据;④通过1∶50 000地质灾害风险调查评价获得各地质灾害易发性数据;⑤资源储量通过搜集各矿区储量年报、采掘工程平面图资料来计算。

上述影响因素具有“三不同”的特点,即领域不同、来源不同、量纲不同,属于多源异构数据。因此,构建综合评价体系首先需要对原始数据进行预处理。原始数据可分为连续性数据和离散型数据。通常对连续性数据采用标准化来消除量纲、取值范围差异的影响[18]。一方面,由于各类影响因素中与地质环境承载力既有正相关的因素,如NDVI、距水系距离、资源可利用量占比等,也存在负相关的因素,如地表切割深度、坡度、地质灾害易发性等;另一方面,各项因素内部原始值分布较稳定,因此选择极差法公式式(21)、式(22)进行标准化,在将区间映射到[0,1]的同时,更能表达影响因素与地质环境承载力之间的正负关系。

(21)

(22)

式中:x′为影响因子标准化值;xmax和xmin分别为影响因子最大值和最小值。

将影响因素标准化后的结果进行整合,进行多重共线性检验。方差膨胀因子(variance inflation factor,VIF)和容差(tolerance,TOL)互为倒数,是多重共线性的常见检验指标[19]。一般认为,当VIF≥10或TOL≤0.1时,说明数据存在严重的共线性,不利用构建模型进行分析。结果如表1所示,17类影响因子间多重共线性较弱,无需再对因子属性进行剔除。最终,研究区数据属性库规格确定为3 553 256×17。

表1 影响因素多重共线性Table 1 Multicollinearity of influencing factors

4 地质环境承载力评价

一方面,评价模型的确定需要考虑研究对象与实际数据之间的适用性;另一方面,可以通过模型融合来弥补方法自身的缺陷,起到对结果优化的作用,常见的模型融合比,如主客观组合赋权、机器学习的集成学习模型等[20-21]。

由于样本数据量庞大,因子空间特征维度较高,而机器学习模型在处理这类数据时具有高效的辨识度,为避免主观赋权方法对评价结果的不利影响,引入数据挖掘半监督学习模型进行综合评价。所谓半监督学习模型,指的是无监督学习模型对样本数据首先进行初步评价,在其结果筛选新样本进行有监督学习模型的构建[22]。首先通过PCA模型进行客观赋权,为尽量避免由于忽略含有较重要信息、方差小的主成分造成的影响,再通过BP神经网络对综合指数叠加的结果进行回归预测。

4.1 无监督学习模型评价

通过IBM SPSS Statistics 25软件的因子分析功能进行主成分分析。样本数据检验结果显示,凯撒·梅尔·奥尔金(Kaiser-Meyer-Olkin-Measure of Sampling Adequacy, KMO)值为0.564(>0.5),Bartlett显著性P为0(<0.5),达到了显著水平,表明该数据适用于主成分分析。通过式(1)～式(6)计算得到特征值λi、方差百分比Pi和累计方差百分比Ai(表2)。按照累计方差百分比达到80%,提取X1～X11主成分代表原有17个因子变量,分别根据式(7)～式(9)计算得到主成分载荷矩阵Zij和各因子权重值T。通过在ArcGIS中综合指数法叠加,计算每一个评价单元的权重乘积和,得到初始地质环境承载力评价结果(图5)。

表2 总方差解释Table 2 Total variance interpretation

4.2 有监督学习模型评价

通过成分矩阵统计,绝对值大于0.5的影响因素作为该成分的主要贡献因子,而譬如地表切割深度和距道路距离2个因子并没有在任一主成分中具有较大贡献值,这就可能因特征值分解造成部分因子对结果贡献的缺失[23]。因此,基于PCA叠加的初步评价结果,通过随机选择样本将目标求解转化为回归预测。

依托Rapid Miner Studio软件来建立BP神经网络模型,该平台的数据挖掘任务涉及广泛,过程从嵌套的算子中产生,最终表示为流程设计,界面简洁,具有运行高效、可视化强的特点[24-25],模型构建结构可视化如图6所示。

inp为输入;out为输出;Optimize Parameters为优化参数;tra为训练;mod为模型;exa为示例集输入、输出;ori为默认输出值;ave为输出性能矢量;unl为预测数据输入端口;par为参数设置;thr:输出端口;lab为添加预测属性;fil为写入CSV文件;per为最优性能矢量;res为目标端口的结果

首先,对研究区随机提取的样本数据进行标签化,将基于PCA模型的初始评价结果作为目标值,按照7∶3将样本数据进行分割,采用十倍交叉验证得到评价性能指标。通过进化算法对BP神经网络进行参数寻优(图6步骤①～④),最优参数如表3所示。通过BP神经网络训练学习(图6步骤⑤～⑧),隐含层的权值计算结果如表4所示。输出层权值为{0.916,-0.369,0.409,-0.949,-0.348,0.445,0.021}。

表3 BP神经网络主要参数

表4 隐含层权值Table 4 Weight of hidden layer

其中,隐含层数通常依靠包含输入层和输出层的经验公式来确定;训练周期代表模型将误差反馈、调整权值过程的重复次数;学习率表示每一次权值改变的大小;动量是指将历史权值更新的一部分添加至当前权重,来有效避免局部极大值并平滑优化方向;终止误差代表模型达到该目标后停止训练;随机种子是指使用局部随机种子进行随机化,来确定初始权值。

最终对研究区数据进行回归预测,将结果输出和存储在格式为.csv文件中,导入ArcGIS Pro中生成地质环境承载力指数(geological environment carrying capacity index,GECCI),如图7所示。

图7 基于BP神经网络地质环境承载力指数Fig.7 Geological environment carrying capacity index based on BP neural network

本文半监督学习评价思路是在无监督学习模型PCA基础上,建立有监督学习模型BP神经网络进行的。换言之,是通过有监督学习模型对无监督学习模型初步评价结果的修正与调整,来弥补无监督学习在降维过程中对部分影响因素的“忽视”,来提高数据挖掘模型对地质环境承载力的识别性。通过性能度量(表5)和结果描述统计(表6),研究发现:①半监督学习模型7项误差指标检验值较小,表明该模型预测精度较高;②GECCI值较无监督学习模型相比,均值和中值均发生降低,而标准差提高,证明数据分布离散程度更高,并且GECCI降低的区域主要分布在小于均值的区域,说明半监督学习模型PCA-BP神经网络提高了地质环境承载力评价的敏感性,半监督学习模型对评价结果的预测、修正符合实际。

表5 模型性能度量指标

表6 评价结果描述统计Table 6 Descriptive statistics of evaluation results

4.3 地质环境承载力区划

目前,GECCI的定义主要为相对值而非绝对值,因此划分等级应侧重于体现组间性质的差异。自然间断法以方差作为衡量分类的标准,充分体现不同分组数据之间的差异性,更适合进行地质环境承载力区划[9,26]。在ArcGIS Pro中对上述GECCI结果进行区划,从低到高划分为5个承载力等级区域(图8),在空间分布上呈现出西北部高、中东部低的特点。

(1)地质环境承载力低区在全区呈条状分布,面积约为309.87 km2,占研究区总面积的9.69%。中部以孤山川、石马川、新城川为主的河谷地区,地势陡峭,岩土体破碎,地质灾害极易发育,该区域也是采矿活动密集地带;东部以府谷县城区为主,人口密度最为集中,地表破坏强烈,导致沿黄公路一线均为低承载力区;北部呈段状,零散分布在清水川、石峡沟一带;西部以大昌汗沟和大板兔川沟谷区为主,局部由采煤活动造成。

(2)地质环境承载力较低区呈现条块状,面积约708.32 km2,占总面积的22.15%。全区分布位置与承载力低区近似,为低区的外围延伸地带,在东北部分布集中,由于道路工程密集,在清水镇、黄甫镇呈片状分布。

(3)地质环境承载力中等区分布最广,面积约1 024.60 km2,占总面积的32.04%。主要分布在东北、西南两侧,东北部地质灾害发育较少,降雨条件相对充沛,植被覆盖程度较高;西南侧地势较为平缓,煤炭资源储量丰富,人口分布较少,受断层影响小。

(4)地质环境承载力较高区面积约796.02 km2,占总面积的24.90%。在中部、北部及西部承载力高区外围地带呈团状分布,多为境内河流上游支沟的沟顶区域,植被发育条件较好;在南部黄土梁顶亦集中分布,几乎无地质灾害分布。

(5)地质环境承载力高区以块状集中分布,面积约358.57 km2,占总面积的11.21%。主要分布在中部木瓜镇、孤山镇的平坦地区,受沟谷切割影响小,无地质灾害发生;侏罗纪、石炭——二叠纪煤炭资源埋藏区的承载力水平高;在北部哈镇、古城镇的边界地带,承载力高区亦集中分布,该区域对原始地质环境影响较轻。

5 影响因素作用分析

5.1 影响因素作用规律

地质环境承载力涉及环境、资源、社会等多学科系统,评价体系涵盖面广,内容庞杂。如何量化影响因素对地质环境承载力的作用程度,是识别地质环境承载力主控因素的重要步骤,亦是对同类地质环境背景地区的承载力评价体系结构的优化。

利用PCA-BP神经网络的地质环境承载力评价结果,结合上述影响因素与承载力的正负作用关系,将影响因素取值范围划分为5个等级,分别赋值1～5对应承载力由低到高。在Model Simulator界面动态提升某影响因素等级(图9),对应统计GECCI的增量百分比,来反映该等级下某影响因素对地质环境承载力的作用程度(图10)。

Input for Model为模型输入因子;Prediction为预测结果;Supports Prediction为支持预测程度;Contradicts Prediction为反对预测程度;CSI为地面塌陷易发性;DFSI为泥石流易发性;DT_ratio为采深采厚比;D_faults为距断层距离;D_rivers为距河流距离;D_roads为距道路距离;LSI为滑坡崩塌易发性;NDVI为植被归一化指数;P_density为人口密度;resource为资源可利用量占比;Important Fators for Prediction为影响因子预测的重要性

通过图10可以看出,部分影响因子在不同等级模拟中,对GECCI的影响程度差异较大;而另一部分因子对GECCI的作用相对不明显。主要规律如下。

(1)单因素承载力为低时,灾害地质条件是限制地质环境承载力的决定性因素,如崩塌滑坡易发性、地面塌陷易发性以及泥石流易发性的承载力提升时,分别引起GECCI提升3.89%、5.48%和2.69%。而随着单因素承载力的上升,灾害地质条件对GEECI的总体约束能力逐级降低,如在单因素承载力较高至高变化时,这3个影响因素分别导致GECCI提升2.98%、2.44%和2.98%。人口密度、距水系距离、土地利用类型以及坡度同样表现出上述规律,即承载力低等级对GECCI的作用程度大于承载力高等级。

(2)与之相反的是GDP密度、资源可利用量占比和距道路距离三类因子,在承载力等级为低时,对GEECI的影响分别为0.51%、0.17%和1.29%;而在高承载力下作用程度分别提升至8.55%、8.49%和3.52%,社会经济条件和矿产资源条件成为对地质环境承载力作用最显著的影响因子。

(3)其余7类影响因子,包括地表切割深度、采深采厚比、距断层距离、NDVI、地层岩性、曲率以及年均降雨量,对GECCI的作用程度相对较轻,介于0.17%～1.05%之间,并且在不同单因素承载力等级下,没有显示出较强的变化规律。

5.2 影响因素全局作用

在上述影响因子作用规律的基础上,最终通过Explain Predictions算子得出各影响因子的全局权重(图11)。

图11 地质环境承载力综合影响权重Fig.11 Comprehensive influence weight of geological environment carrying capacity

可以看出,研究区地质环境承载力主要影响因素为崩塌滑坡易发性、地面塌陷易发性、GDP密度、泥石流易发性、资源可利用量占比、距水系距离以及人口密度,相对重要性分别占12.59%、11.13%、10.11%、10.09%、9.71%、9.30%以及8.19%。这7类影响因素累积重要性占比超过70%,在黄土高原水土保持区的后续相关研究中,可直接作为影响因素构建评价体系。

6 结论

(1)半监督学习模型训练与验证结果表明,半监督学习模型PCA-BP神经网络进行地质环境承载力评价预测精度较高,较无监督学习模型相比,提高了对地质环境承载力评价的敏感性。

(2)府谷县地质环境承载力空间分布上呈现出西北部高、中东部低的特点。等级划分为高承载力区、较高承载力区、中等承载力区、较低承载力区和低承载力区,分别占总面积的11.21%、24.90%、32.04%、22.15%以及9.69%。

(3)影响因素作用分析结果表明,地质环境低承载力区主要受地质灾害控制;而地质环境高承载力区主要社会经济和矿产资源条件决定。

影响因素全局作用为评价体系的优化提供依据。结果表明崩塌滑坡易发性、地面塌陷易发性、GDP密度、泥石流易发性、资源可利用量占比、距水系距离以及人口密度共7项影响因素可优先用于黄土高原水土保持区地质环境承载力研究。