0.   引　　言

煤炭是我国的主体能源。据统计，2023年我国原煤产量47.1亿t，同比增长3.4%，煤炭消费量约为31.6亿t，占能源消费总量的55.3%^[1]。随着煤炭长期开采，一些产煤地区逐渐形成采选、发电、化工等涉煤产业集聚区，其中各类工业场地和固体废弃物堆场如：煤矿开采场、发电厂、焦化厂、煤矸石与粉煤灰堆场等都可能对土壤与地下水造成严重污染^[2-4]。该类地区地下水污染具有污染源分布集中，影响范围相互叠加，污染机制复杂等特点^[5-7]，因此开展地下水污染风险评价对保护该类地区地下水环境具有重要意义^[8]。近年来，国内外相关学者对煤矿集聚区内部分典型场地如：闭坑矿井^[9]、尾矿库^[10]、露天煤矿^[11]等小范围或特定边界的地下水污染开展风险评价研究，然而对于涉煤产业集聚区全区域、多含水层、多类型污染源、多种类特征污染物的地下水污染风险评价尚缺乏研究。

地下水污染风险定义为含水层本质脆弱性和由于人类活动对其造成污染达到不可接受水平的可能性^[12]。当前国内外常用的地下水污染风险评价方法主要包括：迭置指数法、过程模拟法和统计方法，其中迭置指数法因方法简单、操作性强的特点，是应用最广泛的评价方法^[13]。迭置指数法是将风险用几类指标来表征，按照特定的评分原则得到对象的风险指数，最终根据风险指数进行分级。然而传统的迭置指数法在评价指标选取、权重计算与等级划分上往往具有较强的主观性^[14]。随着机器学习技术的不断发展，越来越多的学者采用机器学习模型预测地下水污染风险^[15]，该方法可有效避免主观作用的影响。随机森林（Random Forest, RF）^[16]作为集成机器学习模型的一种，对不确定性和非线性问题具有较好的预测效果，采用随机森林方法可替代传统的迭置指数法对地下水污染风险进行预测。

笔者以山西某涉煤产业集聚区为研究区，采用DRASTIC与PLEIK模型分别对研究区孔隙水与岩溶水脆弱性进行评价，采用层次分析法确定PLEIK模型各指标权重，综合应用污染源荷载与地下水脆弱性计算该地区地下水污染风险，并根据研究区取样点水质等级，使用随机森林（RF）分类预测算法，构建地下水污染风险分级预测方法。本研究以期为后续开展地下水污染防治重点区划定与地下水污染风险管控提供依据和参考。

1.   材料与方法

1.1   研究区概况

研究区位于山西省中部，总面积约为2200 km²，含煤面积1662.55 km²，煤炭资源储量202.2亿t，是我国重要的炼焦煤基地。研究区地理位置及概况如图1所示。研究区高程范围为749～2197 m，地形大势呈西北高东南低。多年平均降水量为400～600 mm，汾河及其支流为研究区内主要河流水系。研究区内地下水类型可划分为3类：第四系松散岩类孔隙水、煤系碎屑岩类裂隙水和奥陶系碳酸盐岩类岩溶水。其中，第四系松散岩类孔隙水与奥陶系碳酸盐岩类岩溶水是具有供水意义的含水层。碳酸盐岩类岩溶水在全区广泛分布，汾河以北地区分布有大面积的裸露灰岩，接受大气降水的直接补给，并在研究区东南部以泉水形式出露于地表。孔隙水主要分布在山间谷地、汾河及其支流河谷以及山前地带，除了山前地带的第四系含水层连续分布外，其它地区分布范围小。

图  1  研究区地理位置

注：中国地图审图号：GS(2019)1822号

Figure  1.  Geographic location of the study area

下载: 全尺寸图片 幻灯片

本研究共收集了165个涉煤污染源资料，其中包括16个煤矿开采区、127个工业污染源（包括洗煤厂、焦化厂、发电厂等）以及22个危险废物堆场（包括粉煤灰与矸石堆场）。研究区是典型的煤炭产业集聚区，区内形成了“采煤−洗煤−发电”，“采煤−洗煤−炼焦−化产”等产业链，涉煤产品种类丰富，包括焦炭、冶金焦、铸造焦、焦油、粗苯、煤气、洗精煤和硫铰等。研究区内地下水主要的污染模式包括：在产煤矿矿井水串层污染、矸石与粉煤灰堆场淋滤液入渗、废弃煤矿矿井水自流外排、煤炭洗选场、焦化厂与燃煤发电厂等场地生产过程中的跑冒滴漏等，其污染模式如图2所示。

图  2  研究区地下水污染模式

Figure  2.  Patterns of groundwater contamination in the study area

下载: 全尺寸图片 幻灯片

1.2   地下水污染风险评价方法

地下水污染风险可表征为地下水污染发生概率与污染危害性的乘积^[17]。笔者通过在地下水脆弱性的基础上叠加地下水污染源荷载计算地下水污染风险指数，以表征地下水污染风险大小，计算公式如下：

式中：$I_{\rm{R}} $为地下水污染风险指数；$I_{\rm{D}} $为地下水脆弱性指数；$I_{\rm{P}} $为地下水污染源荷载综合指数。

1.2.1   地下水污染源荷载评价

地下水污染源荷载评价可用于评估人类活动对地下水潜在危害的程度^[18]。通过污染物毒性、污染物释放可能性和可能释放污染物的量，计算得出单一污染源荷载评价结果，计算公式如下：

式中：P为单一污染源荷载指数；T为污染物毒性；L为污染源释放可能性；Q为可能释放污染物的量。

根据研究区内污染源的分布和污染特征，评价综合污染源荷载等级，将各单一污染源的计算结果加权求和形成污染源荷载综合指数，计算公式如下：

式中：W_i为第i个污染源类型的权重；P_i为第i个污染源的荷载指数。PI值越大，表示地下水污染源荷载越大。

1.2.2   地下水脆弱性评价

根据研究区内地下水类型，本研究中脆弱性评价主要针对2个含水层：第四系松散岩类孔隙水与奥陶系碳酸盐岩类岩溶水，针对不同类型的地下水应选择不同的评价模型。

1)孔隙水脆弱性评价。美国EPA于1987年提出的DRASTIC模型来评价研究区孔隙水脆弱性，该模型具有操作简便、模型参数易获取以及适用性广泛的特点，在平原地区孔隙与裂隙水含水层的脆弱性评价中具有显著优势，成为国内外最常用的脆弱性评价模型之一^[15,19]，研究区内孔隙水主要分布在研究区东部山前平原地带，因此选用DRASTIC模型对研究区孔隙水脆弱性进行评价，其计算公式如下：

式中：D、R、A、S、T、I和C分别表示地下水埋深、净补给量、含水层厚度、土壤介质类型、地形坡度、包气带介质类型和含水层渗透系数。下标r表示因子评分，下标w表示因子权重，$I_{\rm{D,K}} $为孔隙水脆弱性指数，值越大，表明孔隙水脆弱性越高。

2)岩溶水脆弱性评价。DRASTIC模型在岩溶含水层脆弱性评价中的具有一定的局限性，其中的部分参数在岩溶含水层适用性较差，因此采用PLEIK模型来评价集聚区岩溶地下水的脆弱性。该模型在指标选取上包含了岩溶含水层的基本特征，是当前国内外常用的岩溶水脆弱性评价的模型之一^[20-21]。该模型由保护性盖层、土地利用类型、表层岩溶带发育强度、入渗补给强度、岩溶网络发育程度5个指标构成，其计算公式如下：

式中：P、L、E、I和K分别表示保护性盖层、土地利用类型、表层岩溶带发育强度、入渗补给强度、岩溶网络发育程度；下标r表示因子评分，下标w表示因子权重，$I_{\rm{D,Y}} $为岩溶水脆弱性指数，值越大，表明岩溶水脆弱性越高。

根据《环境影响评价技术导则 地下水环境》(HJ 610-2016)中的地下水环境敏感程度分级表，将研究区内3个泉域保护区设定为敏感区，并在计算脆弱性结果时分别赋予1.5倍权重^[22]，以提升评价结果的准确性。

1.2.3   层次分析法

本研究采用层次分析法确定PLEIK模型各指标权重，具体步骤如下：

1)构建评价矩阵。本研究借助专家打分法，通过两两比较的方式，确定i要素对j要素的重要程度（表1），从而构建评价矩阵X(n)，n为矩阵阶数；

表  1  层次分析法各因素打分表

Table  1.  Hierarchical analysis factors scoring scale

因素i比因素j	量化值
同等重要	1
稍微重要	3
较强重要	5
强烈重要	7
极端重要	9
两相邻判断的中间值	2、4、6、8

下载: 导出CSV 

| 显示表格

根据覆盖型岩溶区水文地质条件，PLEIK模型中5个指标的相对重要性排序为：保护性盖层>土地利用类型>表层岩溶带发育强度>入渗补给强度>岩溶网络发育程度；

2)权重计算。对评价矩阵的特征值与特征向量进行计算：

式中：λ、w分别为矩阵X的特征值、特征向量。本文使用Matlab软件计算评价矩阵的特征值与特征向量；

3)一致性检验。由于构建的判断矩阵具有主观性，因此需要对判断矩阵进行一致性检验，检验步骤如下：

式中：$I_{\rm{C}} $为一致性指标；λ_max为矩阵最大特征值；n为矩阵阶数；$R_{\rm{C}} $为一致性比率；$I_{\rm{R}} $为平均随机一致性指标，其取值见表2。当$R_{\rm{C}} $越小时，表明判断矩阵的一致性越好，当$R_{\rm{C}} $<0.1时，即可以认为判断矩阵具有一致性。

表  2  平均随机一致性指标标准值

Table  2.  Mean randomized consistency indicator standard value

阶数	1	2	3	4	5	6	7	8
IR	0	0	0.58	0.90	1.12	1.24	1.32	1.41

下载: 导出CSV 

| 显示表格

4)归一化处理

对特征向量进行归一化处理，计算方法如下：

式中：W为第i个指标对应的权重值，w为λ_max对应的特征向量，最终组成的向量(W₁，W₂，…，W_n)为评价矩阵各指标的权重向量。

经以上方法计算后得出PLEIK模型各指标权重，见表3。

表  3  PLEIK模型各指标权重

Table  3.  Weights of the indicators of the PLEIK model

指标	P	L	E	I	K
权重	0.28	0.25	0.19	0.16	0.12

下载: 导出CSV 

| 显示表格

1.3   数据来源

研究区污染源数据均来自资料收集与地下水采样测试；地下水埋深D数据来自于地下水监测井统计数据；地下水净补给量R基于降水入渗计算；含水层介质类型A与包气带岩性I来源于钻孔数据；土壤介质S数据来源于世界土壤数据库(HWSD)；地形坡度T数据来源于地理空间数据云的dem数据；含水层渗透系数C根据含水层类型经验值换算求得；保护性盖层P来自于文献^[23]；土地利用类型L数据来源于中国科学院资源环境科学与数据中心^[24]；其余数据来自山西省水文地质资料。

1.4   随机森林分类预测

利用课题组现场采样获取的44组地下水水样预测污染发生概率（取样点位置见图1），将研究区的地下水分为：低污染风险（地下水水质I−III类）、中污染风险（地下水水质IV类）和高污染风险（地下水水质V类）3类^[15]。利用Matlab软件中的Machine Learning工具箱构建随机森林（RF）^[25-26]分类预测算法，将分类后的采样点水质数据按照3∶1的比例划分训练集与测试集，将地下水脆弱性评价模型与污染源荷载评价模型中的各指标作为输入层，将污染概率作为输出层。通过调试RF模型中决策树数量、分裂所需最小样本数以及最大深度3个超参数，多次训练后使模型预测准确性达到最高，并将训练好的RF模型用于预测研究区地下水污染风险等级。

2.   结果与讨论

2.1   污染源荷载评价结果

采用缓冲区分析、叠置分析、栅格转换、栅格代数计算方法等GIS空间分析方法，计算出研究区各类污染源荷载值，并利用图层叠置分析，识别重叠区域并计算累积污染荷载指数。通过栅格转换，将各类别污染荷载指数转换为栅格模型。根据各类别污染源权重，通过栅格代数计算加权求和得到综合污染荷载指数I_P。根据I_P计算结果，将研究区地下水污染源荷载等级分为低、中和高3个等级，划分标准见表4，评价结果如图3所示。

表  4  地下水污染荷载等级划分标准及占比

Table  4.  Pollution load classification criteria and specific gravity

荷载等级	IP	占比
低	≤40	41.08%
中	(40, 60]	32.19%
高	＞60	26.73%

下载: 导出CSV 

| 显示表格

图  3  研究区污染源荷载空间分布

Figure  3.  Spatial distribution of pollution loading

下载: 全尺寸图片 幻灯片

由图3可知，研究区地下水污染源低荷载区、中荷载区和高荷载区的占比分别为41.08%、32.19%和26.73%。高污染区主要分布在研究区北部煤矿开采区与东南部太原盆地等处，造成这些地区污染源荷载高的原因主要与涉煤产业集聚区污染源分布集中的特点有关，多个污染源的叠加效应明显^[27]，且涉煤污染源的特征污染物多以重金属与多环芳烃为主^[28]，毒性指标较其他类型污染源偏大，因此在量化研究区污染源荷载时其值偏高。研究区西南部污染源较少，污染荷载相对较低。

2.2   地下水脆弱性评价结果

2.2.1   孔隙水脆弱性评价结果

根据研究区水文地质条件，研究区第四系孔隙含水层分布如图4a所示。将相关数据进行处理得到DRASTIC模型各因子评分图（图4b—图4h）。根据式(4)，将DRASTIC模型各指标分值与其对应的权重相乘并叠加，生成孔隙水脆弱性指数，并划分为低、中、高3个等级，等级划分标准和面积占比见表5，分区分级结果如图4i所示。

图  4  DRASTIC模型评价结果

Figure  4.  DRASTIC model evaluation results

下载: 全尺寸图片 幻灯片

从图4i可以看出，对于研究区内的孔隙水含水层，低脆弱区主要分布在研究区西北部山区，该地区地下水埋深较大，且坡度较高，降水与地表径流不易汇集，地表污染物较难进入地下水，因而地下水脆弱性低^[29]。高脆弱区主要分布在研究区东部与东南部太原盆地地区，该地区地下水脆弱性偏高的主要原因是该地区地势平缓，地表污染容易在此处汇集，地下水埋深浅，污染物较容易进入地下水，且包气带及土壤介质渗透性强，对污染物的吸附能力较差。

2.2.2   岩溶含水层脆弱性评价结果

将研究区相关数据进行处理得到PLEIK模型各因子评分图（如图5a—图5e）。根据式(5)将PLEIK模型各指标分值与上述权重相乘并叠加，生成岩溶水脆弱性指数图，并划分为低、中等、高3个等级，等级划分标准与占比见表5，分区分级结果如图5f所示。

图  5  PLEIK模型评价结果

Figure  5.  PLEIK model evaluation results

下载: 全尺寸图片 幻灯片

从图5f可以看出，对于研究区内岩溶水含水层，脆弱性等级主要为中，少量高脆弱区主要集中在研究区西北部的灰岩裸露区、汾河渗漏段以及煤层带压开采区，该地区岩溶发育程度较高，地层连通性好，降水与地表径流可直接补给，使得污染物与更易进入地下水中，这一结果与前人研究一致^[30-32]。研究区东部太原盆地上覆有第四系孔隙含水层，因而污染物不易进入下层的岩溶水，但该区域地形平坦、靠近城区，土地利用类型以工业用地为主，人为活动影响较大^[33]，因此拉高了该地区的脆弱性。研究区低脆弱性区占比较少，主要分布在研究区西侧。

2.2.3   地下水脆弱性综合评价结果

当评价范围内有多种地下水类型交互存在时，本着“应纳尽纳”的原则，将DRASTIC模型与PLEIK模型评价结果叠加取高值，得到研究区地下水脆弱性综合评价结果（如图6），研究区地下水综合脆弱性各等级占比见表5。研究区地下水低脆弱区、较低脆弱区、中等脆弱区、较高脆弱区和高脆弱区的占比分别为6.33%、82.59%和11.08%，中脆弱性区域占比约超过80%，表明研究区地下水脆弱性整体为中，高脆弱性地区主要分布在研究区东部与东南部的太原盆地以及西北部岩溶裸露区。

图  6  研究区地下水综合脆弱性评价结果

Figure  6.  Results of the groundwater vulnerability assessment for the study area

下载: 全尺寸图片 幻灯片

表  5  研究区地下水脆弱性评价分级划分标准及面积占比

Table  5.  Criteria for grading groundwater vulnerability assessment and specific gravity in the study area

脆弱性分级	孔隙水脆弱性			岩溶水脆弱性		综合脆弱性 评价结果占比
	ID,K孔隙	占比		ID,Y岩溶	占比
低	≤100	57.89%		≤4	4.97%	6.33%
较低	(100, 120]	12.73%		(4, 6]	89.09%	82.59%
高	＞120	29.38%		＞6	5.94%	11.08%

下载: 导出CSV 

| 显示表格

2.3   地下水污染风险评价结果

根据式(1)，将研究区污染源荷载结果与地下水脆弱性结果叠加取高值，将研究区地下水污染风险等级划分为低、中等和高3个等级。从地下水污染风险等级分区（图7）可知，研究区地下水污染低风险区、中风险区和高风险区的占比分别为3.55%、59.67%和36.77%，地下水污染风险以中风险为主。研究区地下水高风险区占比34.91%，而高污染荷载和高脆弱性区域占比分别为26.73%和13.29%，高污染风险区所占比例接近高污染荷载和高脆弱性区域占比之和，由于研究区地下水污染风险是由污染源荷载结果与地下水脆弱性结果叠加取高值所得，这表明研究区内高污染荷载与高脆弱性区重叠部分较少。然而有研究结果表明^[34-35]，对于其他类型的污染场地，污染源大多分布在地形相对平缓的区域，这就导致这类场地污染源荷载高的地区，其地下水脆弱性往往也相对较高，而本研究区评价结果与上述现象相反，这因为研究区内含煤地层大多埋藏在北部山区中^[36]，且煤矿开采区作为面源污染，占地面积较广、影响范围大，这导致研究区内脆弱性相对较低的山地地区污染源荷载相对较高。

图  7  研究区地下水污染风险评价结果

Figure  7.  Results of groundwater contamination risk evaluation

下载: 全尺寸图片 幻灯片

2.4   随机森林分类预测结果

图8为地下水污染风险分级和机器学习分类预测与采样点实际水质的对比图，从图中可知，基于叠置指数法所划分出的地下水污染风险与采样点实际水质一致率为75%，而经训练后的RF分类预测结果的一致率可达到98%，这表明RF分类预测的结果更能代表研究区地下水的污染风险，因此将训练后的RF模型保存并用于分类预测。在ArcGIS中利用渔网工具箱生成500 m× 500 m的网格，并提取网格中心点处的各评价指标参数，作为输入层，代入使用保存好的RF模型计算输出层，并与构建的网格进行连接，并划分地下水污染风险分区，如图8c所示。从图8b中可以看出，RF预测的地下水污染风险分级与基于叠置指数法计算出的风险分级结果具有较大差异，主要体现在低风险与中风险的划分上，基于叠置指数法计算出的中风险占比超过60%，而RF分类预测的低风险占比超过50%，而2种方法计算出的高风险区的占比大致相同，且都主要位于研究区东部与东南部太原盆地，以及北部岩溶裸露区。

图  8  RF预测地下水污染风险分级结果

Figure  8.  Risk classification results of groundwater pollution based on RF

下载: 全尺寸图片 幻灯片

研究区东部污染风险较高的原因与煤矿开采密切相关，相关研究表明，煤矿开采会导致含水层之间发生串层污染，使得具有高硫、酸性特征的矿井水进入含水层造成地下水的污染^[37]，因此需要采取相应治理手段，如采取封井回填、矿井水治理等^[38]避免污染的发生。研究区东南部为地下水脆弱性较高，且是焦化厂和洗煤厂等工业污染源的集中分布区，多个污染相互叠加影响，其特征污染物以重金属与多环芳烃^[39]为主，对人体健康风险较大，因此有必要对该地区加强管控，从源头上避免污染的发生。研究区北方的环境污染状况较为复杂，除工业污染源和矿山开采区之外，该地区还广泛分布着大量的粉煤灰与矸石堆场，未经处理的粉煤灰和煤矸石长期露天存放，在雨水的淋滤作用下，会向地下水含水层释放有害物质，从而导致地下水污染^[40]。面对这类污染问题，固体废物的资源化利用被认为是一种有效的处理策略，通过将粉煤灰与矸石转化为地质聚合物^[41]，不仅实现了煤矿固废的再利用，而且还能用于矿井水的治理。

综上所述，基于随机森林分类预测算法的涉煤炭产业集聚区地下水污染风险评价结果相较于叠置指数法其预测准确性显著提高，该方法可避免模型指标权重赋值时的主观性影响，使得评价结果更符合研究区实际情况，且算法实现较为简单。然而该方法仍具有一定的局限性，例如为提高模型预测的准确性，需要大量研究区地下水水质样本数据，这意味着现场取样工作会加重。此外，对不同场地的地下水污染风险评价，需重新训练机器学习模型。因此，将机器学习算法与多目标优化算法^[42]结合可有效降低训练时间，对该方法的广泛应用具有重要意义。

3.   结　　论

1）研究区污染源荷载值较高。高荷载区主要集中在5个区域，这些地区污染荷载偏高的原因主要与涉煤产业污染源普遍较为集中的特点有关，在量化污染源荷载时，多个污染源存在空间叠加效应，使得污染源荷载指数相对偏高。

2）研究区地下水综合脆弱性相对较低。对于孔隙水含水层，其高脆弱区主要集中在研究区东南部山前平原地带，该地区地下水埋深浅、坡度平缓、土壤与包气带渗透性高，地下水防污性能较差；对于岩溶水含水层，高脆弱区主要分布在汾河以北表层岩溶裸露区，该地区含水层连通性高，地表污染更容易进入含水层中。

3）研究区地下水污染风险整体相对偏低。低风险区、中风险区和高风险区占比分别为1.73%、63.36%和34.91%，高风险区主要位于研究区东南部化工园区以及汾河以北岩溶裸露区。

4）基于RF预测的地下水污染风险分区结果与实际水质更为接近，其预测准确率相较于叠置指数法提高了约23%，2种方法计算出来的高风险区均主要分布在研究区东部与东南部山前平原地带，以及北部汾河渗漏段与岩溶裸露区，因此需对上述地区开展有效的风险管控措施。