1 引言(Introduction)

地下水水化学是地下水环境的重要组成部分, 人类活动对地下水环境的影响不仅表现在地下水系统的物理变化, 还在很大程度上通过水化学成分的变化体现出来, 进而致使地下水水质产生变化, 即“地下水水化学异常”(王焰新等, 2005; 李培月, 2014).用来评价地下水环境背景值的理想数据是未受人类活动影响的地下水数据(Shand et al., 2007), 但人类活动的广度、深度都在不断地影响着地下水的赋存环境, 因此, 识别出受人类活动影响而产生的水化学异常对客观地获取水化学背景值及确定人类活动对地下水影响程度都至关重要(Gorelick et al., 2010; 郭高轩等, 2010).

对比分析国内外对地下水环境背景值的研究方法, 欧盟分别于2000年和2006年颁布了《水框架指令》(EU WaterFramework Directive, 2000)和《地下水指令》(EU Groundwater Directive, 2006), 针对地下水阈值建立体系中天然背景值的问题做了一定的研究, 提出了异常数据剔除的5条基本准则, 广泛地在欧洲各地区的背景值研究中运用预筛选法(The pre-selection approach)、组分分离法(The component separation method)及运用空间分布的相关方法, 并获取了欧洲不同地区地下水组分的背景值(Müller et al., 2006; Preziosi et al., 2014; Ducci et al., 2016; Caro et al., 2017; Nguyen et al., 2016).一些学者利用累积频率曲线图、模糊聚类法及多元统计法等方法对地下水异常值剔除进行了当地的背景值计算和人类活动影响程度确定等研究(Güler et al., 2012; Cruz et al., 2015; Macdonald et al., 2005; Re et al., 2014; 张英, 2011).宇庆华等(1991)利用拉氏积分函数法和Thompson法对地下水水化学异常的识别进行了探讨研究.近年来部分学者(曾颖, 2015; 彭聪等, 2017; 张小文等, 2017; Zhang et al., 2017)分别通过利用水化学图及Durov图与数理统计相结合的方法, 对柳江盆地及滹沱河流域地下水进行异常识别均取得较好的效果.

英国地调局为响应欧盟《水框架指令》, 针对英格兰、威尔士及苏格兰的地下水开展专门的地下水水化学背景值的研究, 该项目掌握了各区域地下水质量, 其水化学演化进程的主要控制因素, 并为污染影响评价提供条件(Shand et al., 2007; Macdonald et al., 2005), 其地下水质监测井布设原则是尽量避开人类活动剧烈的区域或是利用保护得当专门用来监测地下水水质的观测井.相较之下, 我国地下水水质调查研究起步较晚, 早期的地下水水质资料缺乏, 同时, 我国没有开展专门的背景值研究工作.2005—2015年中国地调局组织实施了全国地下水污染调查评价.本次污染调查主要为1:25万区域调查及部分1:5万重点城市调查, 调查涵盖了中国陆域的主要平原盆地区, 获取了大量的地下水水质样本, 调查范围广, 采样精度小, 已有的水化学图法或Durov图法均充分考虑了水文地球化学的演化过程和水化学组分的内在联系, 但其运用过程复杂, 工作量大, 更适用于区域尺度小且采样精度高的地区.因此, 本研究尝试借鉴英国地调局的关于地下水异常数据剔除的原则, 通过调查采样信息识别出受人类活动影响的地下水异常, 并结合数理统计方法, 力求在大尺度范围内快速精准地对地下水无机常规组分进行异常识别, 并验证其异常识别效果.

2 基本思路与方法(Basic thoughts and methods) 2.1 基本思路

水化学图法利用可以反映地下水某种形成作用的水化学图来识别异常, 如溶滤、混合、离子交换等一系列过程, 利用其各组分间的关系转换为二维坐标, 并运用马氏距离进行异常值的识别, 充分考虑了各水化学组分间的内在联系及水文地球化学演化过程(曾颖, 2015; 彭聪等, 2017; Zhang et al., 2017).Durov图法则是在水化学图法的基础上有所继承和发展(张小文等, 2017).

无论是水化学图法还是Durov图法都对数据的要求严格, 其研究基础均是在1:5万水文地质调查所采集水样的基础上开展的异常识别研究, 采样精度高, 水文地质调查资料详尽, 水文地质单元尺度小, 地下水的补径排条件明确.本次研究基于全国地下水污染调查评价, 其主要调查方式为1:25万区域调查及部分1:5万重点区(中国地调局, 2008), 其调查范围广, 采样精度远远小于1:5万水文地质调查的精度, 大区域内水文地质条件复杂, 其采样点主要布设在地表存在明显污染源的周围.因此, 为了满足大尺度范围、采样精度较低、水文地质条件复杂的地下水水质调查而提出信息筛选法.

数理统计法在背景值研究中是一种高效的工具, 其应用也很广泛(Dalla et al., 2017; Jeong et al., 2017; Cidu et al., 2017).其功能即在研究区进行水文地质分区的基础上, 分析所得水质数据并运用统计学方法进行处理, 最终得出研究区背景值.利用数理统计法识别地下水异常值操作简便, 能够较好地识别出观测序列中明显的离群数据.

2.2 异常识别方法

英国地调局对地下水水化学背景值展开过专门的研究, 并确定了通过一些易检测且能典型地代表某一类污染水化学的理想指示组分(如NO₃、Cl、SO₄等)(Wendland et al., 2008), 将受人类活动影响的数据与天然未受人类活动影响的数据区分出来的方法.信息筛选法借鉴英国地调局的研究思路, 从地下水环境背景值的概念着手, 以水文地质调查资料为起点, 将同一含水层系统中人类活动剧烈、存在明显污染源及采样调查信息不全的水样点剔除, 均不参与背景值的获取研究, 只保留明确说明其周边未发现污染源的水点.本方法主要遵循将受到人类活动影响的数据均识别为异常数据, 相较于英国地调局的方法更为保守.具体筛选方法如下:

1) 根据研究区地下水埋深、含水层富水性及包气带介质等资料将采样数据中同属于同一含水层系统的数据筛选出来;

2) 调查采样资料中表明其地表环境有明显污染源的, 包括化工厂、皮革厂、矿渣、垃圾场、加油站、养殖场等进行筛选并作为异常数据;

3) 调查采样资料中地表水体环境, 内容中包括污染严重的河流、污水沟、排水渠等信息进行筛选并作为异常数据.

本次研究的数据量较大, 所采用数理统计的方法为拉依达准则(张敏等, 1997), 即各水化学指标的范围应该在其均值加减3倍标准差的范围内, 超出该范围的数据则利用拉依达准则迭代识别出来, 并作为异常值剔除.拉依达准则公式, 见式(1):

(1)

式中, X表示数组质量-体积浓度值; X表示数组质量-体积浓度的算术平均值; σ表示数组质量-体积浓度的标准差.

拉依达准则相对于其他众多的数理统计方法更为简单易操作, 拉依达准则对于统计数据数量大时结果更加稳定且保有对观测序列中粗差识别的良好效果(熊艳艳等, 2010).信息筛选法可以快速而又简便地获得大尺度区域的地下水环境背景值, 运用过程中存在潜在污染风险的水点均被信息筛选法剔除, 这对前期的水文地质调查及采样资料有着一定的要求, 同时, 也存在一定的主观性, 因此, 结合拉依达准则从统计学上将离群数据进行识别及剔除是很有必要的.

2.3 异常识别结果

利用异常数据剔除前后的数据, 经过数理统计概括, 得出指标浓度拐点图、水样频率直方图及箱型图, 并以这些图件对比各异常值识别方法识别效果.

指标浓度拐点图是以水点个数为横坐标, 以水点指标浓度值为纵坐标的一系列连接起来的二维坐标图, 它能够清楚地展现出各水化学指标的大小分布概况, 通常会出现斜率较大的拐点位置将数据分为受不同来源影响的数据集群(Morgantini et al., 2009).通过将异常值识别前后数据的拐点图, 可以看出各水化学指标中被识别为异常值的浓度大小及数量情况, 以此来判断异常值识别的效果.

水样频率直方图能够简明扼要地识别出地下水数据的主要分布情况(Runnells et al., 1998), 组距为底边, 以频数为高度的一系列连接起来的直方型矩形图, 用来表示一组数据的分布情况, 在异常值的识别过程中通过数据的分布情况可大致判断出离群数据集与非离群数据集之间的界限.

箱型图提供了一种只用5个点对数据集做简单总结的方式.这5个点包括中位数、Q1、Q3、分部状态的高位和低位, 很形象地分为中心、延伸及分布状态的全范围.它是研究地下水环境背景值的最有用的统计工具, 它能够通过图表展示出数据的集中趋势和突显出极端值的分布, 使其成为背景值研究中重要的指示图件(Tukey et al., 1977).

2.4 异常识别结果对比分析

信息筛选法的提出及其可靠性需要进一步验证, 通过5种水化学图法(piper三线图及各离子比例关系图)、信息筛选法、拉依达准则组合出5种异常值识别方法(水化学图法、信息筛选法、拉依达准则、水化学图法+拉依达准则、信息筛选法+拉依达准则)进行验证研究.

为了量化各方法的识别效果, 使用稳健对抗离群数据的距离值方法, 绝对离差中位数(MAD)是一种采用计算各观测值与平均值的距离总和的检测离群值的方法(Rousseeuw et al., 1998), 利用MAD进一步量化出各组数据相较于集中值的距离值(佘玉萍, 2016), 并利用距离值的集中值计算出各组数据的距离值偏离程度.该方法主要适用于正态分布的数据集, 分布类型为对数正态分布的数据需先进行对数化处理后运用该方法.

在此基础上进行改进, 本次量化研究由以下计算步骤完成:

1) 计算异常识别前后各项指标的中位数;

2) 计算异常识别前后各项指标与中位数的绝对偏差值;

3) 计算异常识别前后绝对偏差值的中位数;

4) 将异常识别前后绝对偏差值与绝对偏差值的中位数的比值作为基于MAD的偏离中心的距离值:

(2)

5) 异常识别前后各项指标最大值的最大距离值与均值差值的比值作为距离值的偏离程度:

(3)

式(2)、(3)中, dv_i为某指标偏离中心的距离值, X_i为某指标实测数据, (X_i)_med为某指标中位数, dv_max为异常识别剔除后距离最大值, dv′_max为异常识别前距离最大值, D_i为距离值偏离程度, 为偏离中心距离值的均值.

距离值偏离程度的相对大小来表示异常值的识别效果, 偏离程度越小代表识别效果越好, 通过上述量化分析结果, 5种水化学异常值识别方法均得到各水化学指标剩余数据的偏离程度, 综合所有指标的量化结果, 借此对比各水化学异常值识别方法的识别效果.

3 实例研究(Case study) 3.1 研究区概况

沙颍河流域是淮河流域最大的子流域, 流域跨河南及安徽两省, 总面积约39075.30 km², 其中, 河南省境内流域面积34467.00 km², 占流域总面积88.21%;安徽省境内流域面积为4608.30 km², 占流域总面积的11.79%.流域边界为东经111°56′44″~116°31′07″, 北纬32°29′24″~34°57′15″.沙颍河发源于河南省登封市嵩山, 流经平顶山、漯河、许昌、周口、阜阳等40余县市, 于颍上县沫河口汇入淮河.

沙颍河流域地形由西北向东南倾斜, 上游为伏牛山脉和外方山脉, 东南部为淮北沙颍河平原区.山区地面海拔高程一般在600~1500 m, 高度变化较大; 东南部平原地势坦荡开阔, 地面高程一般在30~100 m, 坡度较小.研究区地处暖温带向亚热带的过渡地带, 属暖温带半湿润大陆性季风气候区, 冬春干旱少雨, 夏秋闷热多雨, 年平均气温14~15 ℃.山区多年平均降水量达1401.7 mm, 平原区年降水量为600~900 mm.沙颍河流域土地利用类型以平地旱地和平地水浇地为主, 分别占流域平原区总面积的63.45%和23.40%.其他土地利用类型包括旱坡地、平地水田、防护林地、果园、菜地等, 共占流域平原区总面积的13.16%.

沙颍河流域内松散地层广泛分布, 在西部山前一带以洪积冲积为主, 平原区则为湖积冲积或冲洪积, 总厚度受地质构造控制.第四系松散地层内以卵石、卵砾石、砂砾石、粗砂、中粗砂、细砂、粉细砂及粉砂等构成全区地下水含水系统.卵石至中粗砂分布在山前一带, 其余则分布在平原区.粘土、亚粘土和亚砂土构成含水层之间的相对隔水层.含水层由山前至平原区颗粒由粗到细, 单层厚度由厚变薄, 层数由少变多(陈荦, 2013).以地层分层为基础, 水文地质要素为依据, 根据传统水文地质勘查, 认为埋深在50 m左右的粘土层为上下两个含水层组的界限(左正金等, 2007).以50 m为界限, 将整个沙颍河流域含水层系统分为浅层含水层和深层含水层, 深浅层含水层之间广泛分布一层14~20 m的黏性土层, 天然条件下水力联系微弱(陶建华, 2012; 叶念军等, 2009), 本次研究仅保留浅层含水层的数据.研究区区域地下水循环补径排分区如图 1所示, 补给区位于山前倾斜平原, 径流区和排泄区分布在平原区, 地下水主要由大气降水补给, 水力坡度在1/8000~1/12000之间, 地下水水平径流流速非常小, 排泄方式天然条件下主要以蒸发为主, 说明不同水文地质单元间的水力联系并不紧密.局部水循环是影响它的水质的主要因素.

3.2 数据预处理

为确保进行统计及异常识别研究数据的质量, 需要对沙颍河流域的每个地下水样点的室内测试数据进行明显错误检查, 电荷平衡检查和碳酸平衡检查, 将检查质量不合格的样品进行剔除.

电荷平衡检查根据阴阳离子平衡检查公式(式(4)), 将CBE>10%的样品剔除.

(4)

式中, CBE是指电荷偏离误差; z是指离子电荷数; m_c是指阳离子(Na⁺、K⁺、Ca²⁺、Mg²⁺、Fe²⁺、Fe³⁺、Mn²⁺)的摩尔浓度; m_a是指阴离子(HCO₃^-、CO₃^2-、SO₄^2-、Cl^-、NO₃^-、HSiO₃^-)的摩尔浓度.

碳酸平衡检查是根据水质分析表中的测试指标, 当pH < 8.34, CO₃^2-的浓度应小于目标检出限(3 mg·L^-1), 若值大于目标检出限值, 则认为碳酸不平衡, 并将碳酸不平衡样品剔除.

以沙颍河流域为例, 基于2006—2010年地下水水样采集工作, 经过上述数据质量检查后保留684组有效地下水水样数据, 其中包括601组浅层地下水数据, 83组深层地下水数据.本次地下水水化学异常值识别研究以浅层地下水为研究对象.

3.3 异常识别结果

为了明确地下水环境背景值的特征和演化规律, 根据地质结构、岩石性质、地形地貌条件、水文地质条件、包气带介质、含水层富水性等条件(叶念军等, 2009)将研究区域分为若干个水文地质亚区.沙颍河流域分为3个水文地质亚区, 分别为补给区、径流区及排泄区, 其中补给区共有180组数据, 径流区共有267组数据, 排泄区共有154组数据.运用信息筛选法、水化学图法分别与拉依达准则相结合不同组合的异常值识别方法对研究区内的有效浅层地下水数据进行异常识别研究, 结果见表 1.其中拉依达准则识别出异常数据个数最少, 而采用信息筛选法+拉依达准则识别出的异常数据个数最多.

以沙颍河流域补给区水文地质亚区为例, 利用Origin 8.5对该区水化学异常值识别结果进行分析.通过沙颍河流域补给区TDS浓度分布散点图(图 2a)及箱型图(图 2b)可以看出, 单独运用信息筛选法和拉依达准则的识别效果并没有将TDS浓度远远高于该区域的异常点识别出来, 识别异常后的剩余数据离散程度仍较高, 其余3种方法所得到的阈值差别不大, 识别效果较为一致.进而对比径流区和排泄区的相关图件, 可以发现3个水文地质分区中, 径流区的TDS识别效果, 仅单独运用信息筛选法的效果较差, 其余方法的识别效果较为一致; 排泄区的TDS识别效果, 水化学图法+拉依达准则和信息筛选法+拉依达准则的异常识别效果较之其余异常识别方法更好.综合所有无机常规组分的识别效果图件得出, 信息筛选法+拉依达准则及水化学图法+拉依达准则对异常的识别效果要优于其他方法.

对比运用异常值识别方法前后研究区补给区地下水水质中TDS的水样频率直方图, 见图 4, 可以看出大于1100 mg·L^-1的水点基本被识别为异常点, 说明此类水点所受人类活动影响大.单独使用信息筛选法或是拉依达准则所得结果中仍可见2000 mg·L^-1和2500 mg·L^-1附近有分布明显离群的数据, 因此信息筛选法或是拉依达准则的单独运用并不能对异常值的识别取得良好的效果, 而水化学图法和信息筛选法分别与拉依达准则的结合能够有效地识别出离群的异常值.

3.4 异常识别结果量化

要进行数据偏离程度的计算, 必须首先进行数据集的分布类型验证, 本次研究中拟推荐使用Kolmogorov-Smirnov(K-S)检验及Q-Q图法进行综合验证数据的分布类型, 首先通过SPSS 20选择数据集进行单样本K-S检验, 验证数据集是否满足正态分布假设, 而后再进一步利用Q-Q图更为直观地验证数据分布类型, Q-Q中各指标的点近似地在一条直线附近, 即说明整个数据符合此种分布类型, 通过Q-Q图上的数据情况进一步佐证沙颍河流域无机常规组分数据分布类型为正态分布或对数正态分布(图 4).

在沙颍河流域补径排3个水文地质亚区中, 确定pH的数据分布类型为正态分布类型, 其余指标数据分布类型为对数正态分布类型, 对pH以外的指标数据进行对数化处理, 再利用距离值偏离程度对各水文地质亚区异常识别后的剩余数据集群进行识别效果量化处理, 各方法的评分情况见表 2~表 4.可以看出, 5种异常值识别方法中, 单独运用信息筛选法和拉依达准则的量化结果表明其异常识别效果较于其他识别方法最差, 单独运用水化学图法所得到的异常识别结果较好, 水化学图法+拉依达准则和信息筛选法+拉依达准则所得的结果综合偏离程度最低, 即异常识别效果最好, 两种方法所得剩余数据的距离值偏离程度在不同水文地质单元内各有优劣.单独运用某种方法表现出略差于不同方法的结合运用, 表明不同方法的复合使用在识别异常值的研究中效果更好.

3.5 异常识别结果分析

将沙颍河流域3个水文地质单元所进行的异常识别结果整理得到异常识别结果空间分布状况, 见图 5.从整个沙颍河流域的异常值识别情况来看, 水化学图法+拉依达准则共识别出262组异常数据, 信息筛选法+拉依达准则共识别出276组异常数据, 两种方法共同识别出187组异常数据.其中水化学图法+拉依达准则单独识别出异常数据75组, 信息筛选法+拉依达准则单独识别出异常数据89组.共同识别出的异常数据187组数据, 其调查资料显示均为受人类活动影响较大或地表存在明显污染源; 水化学图法+拉依达准则所单独识别出的75组异常数据, 通过其水文地质调查采样资料可以看出, 这些采样点的周围环境为农业区, 周边无明显污染源.从统计数据上来看, 信息筛选法+拉依达准则单独识别出的异常数据与水化学图法+拉依达准则单独识别出的异常数据数量上均较多.为进一步解释此现象并验证异常识别结果, 以TDS为例绘制指标浓度拐点图.

通过图 6可以看出, 信息筛选法+拉依达准则的方法相较于水化学图法+拉依达准则的方法剔除的数据更多, 但两种方法的结果对比中可以看出, 浓度较高的离群数据均被这两种异常识别方法所识别及剔除, 以各区内TDS的异常识别结果为例, 其中补给区内大于1200 mg·L^-1的数据均被识别为水化学异常数据, 径流区内大于1500 mg·L^-1的数据均被识别为水化学异常数据, 排泄区内大于1100 mg·L^-1的数据均被识别为水化学异常数据, 异常识别结果也呼应了研究区水平径流流速非常小、不同水文地质单元间水力联系不紧密等水文地质条件的特点.两种方法所得剩余数据的值域范围相差不大, 信息筛选法所得剩余数据略高于水化学图法所得, 总体上对背景值的获取影响不大.但是, 信息筛选法+拉依达准则的方法筛选结果明显向Y轴偏移, 说明在数据量上得到了一定程度的抽稀, 其主要原因是根据信息会将地表环境存在污染源的, 但地下水不一定受到影响的点也作为异常剔除.

利用水化学图法侧重于水化学组分的内在联系, 其识别过程无论是人类活动影响所引起的人为因素导致的人为异常还是自然因素所引起的天然异常均被识别并剔除, 其应用范围更适用于水文地质条件相近、地域尺度小、采样精度高的研究区域; 信息筛选法侧重于人类活动影响可能引起的水化学异常, 地下水可能受到人类活动影响而并未表现出水化学特征异常的均被识别并剔除.信息筛选法+拉依达准则识别异常的速度快, 识别过程更加简便易操作, 识别异常数据的结果也取得不错的效果.通过不同异常识别方法的量化对比分析得出, 水化学图法+拉依达准则和信息筛选法+拉依达准则均能有效地识别出获取背景值过程中的异常数据, 运用信息筛选法+拉依达准则能够有效地替代水化学图法+拉依达准则, 既确保了异常识别效果的可靠性, 又简化加速了异常识别过程.

4 结论(Conclusions)

1) 实例研究表明运用信息筛选法+拉依达准则能够有效地识别出地下水水化学常规无机组分异常, 并且信息筛选法能够有效地替代水化学图法, 在异常识别的研究中有一定的现实可行性.

2) 信息筛选法结合拉依达准则的异常识别方法主要适用于1:25万等大尺度范围地下水污染的区域调查, 此类调查采样精度较低, 区域内水文地质条件复杂, 具备一定的水文地质调查资料、地表污染源调查资料及地表水体调查资料.

3) 异常识别效果量化对比分析结果表明, 运用水化学图法及信息筛选法等异常值识别方法时, 结合数理统计学的方法对于水化学异常的识别效果较之单独运用某一种异常识别方法的效果更佳.

致谢(Acknowledgement): 感谢中国地质科学院水文地质环境地质研究所孙继朝教授团队在本研究中给予的大力支持.