前言:　　 砷是自然界中广泛存在的一种类金属元素,是国际癌症研究机构确定的人类致癌物,可以引起皮肤癌、肺癌和膀胱癌等。人类可以通过多种途径接触到过量的砷,其中饮用受砷天然污染的地下水是最重要的途径。大量文献研究表明,地下水砷污染引起的慢性砷暴露可以严重影响人类健康,包括各种皮肤损伤、高血压、心血管疾病、糖尿病以及乌脚病等。目前全球约有2亿人正处于慢性饮水型高砷暴露的危险中,因此地下水砷污染已经成为全世界广泛关注的一个公共卫生问题。　　 我国是世界上受地下水砷污染影响最为严重的国家之一。据估计,我国大陆地区饮水型砷中毒病区中的高砷暴露人群已超过300万,随着调查的深入,慢性饮水型高砷暴露人群的数目正在不断扩大。自上世纪八十年代,在我国新疆奎屯地区发现第一例砷中毒患者后,又陆续在内蒙古河套平原和山西省大同盆地发现了新的饮水型砷中毒病区。鉴于地方性砷中毒的广泛性和严重危害,我国政府于1992年将其列为一种新的地方病,并投入大量的资金用于筛查和发现地下水砷污染区和受累及的人群。为了有效的筛查和识别潜在的地下水砷污染区,我国研究者研制开发了“10%抽样法”,并广泛应用于全国各地的地下水砷污染区调查中。最近的调查结果显示,在全国16个省292个县(约占全国总县数的12%)20517个村的445638口井中,有近5%的水井砷浓度超过了我国的饮用水砷浓度标准(50μg/L)。虽然“10%抽样法”在很大程度上加速了地下水砷污染区的发现,但是由于我国国土面积巨大,要想完成全国范围内的筛查,仍然需要几十年的时间。而一种可以预测某一地区地下水中是否出现砷污染的空间模型,可以加快水砷筛查速度,因为它可以指示出潜在的高砷地下水分布在哪些地区,指导相关部门或组织在这些地区进行优先采样。　　 最近的研究表明,基于水砷浓度和相关环境解释变量(如地质、气候、地形等)统计关系的地下水砷污染预测模型在东南亚地区成功建立。而Logistic回归也被成功的应用于预测地下水中砷浓度超过某一规定阈值的概率。　　 山西省是我国北方地区受饮水型地方性砷中毒影响最为严重的省份之一。与新疆和内蒙不同的是,山西省的人口密度大,且多集中生活在该省从北到南的几个盆地地区。自1994年在山阴县首次发现饮水型高砷暴露地区以来,经过多年的调查,初步确认大同盆地和太原盆地是该省两个最主要的饮水型地方性砷中毒病区。有大于90万人居住在这两个盆地,并且有3998个人被诊断为饮水型地方性砷中毒患者。这些患者主要分布在两个盆地的应县、山阴县、朔城区、文水县、汾阳市、平遥县和孝义市。除了这些已发现砷中毒患者的地区外,另外还有12个县存在高砷地下水。研究报道称在大同盆地采集的288份地下水水样中,砷浓度的范围是0.2-720μg/L;在太原盆地的采集的44个水样的砷浓度范围在0.1-115.5μg/L之间。　　 目前的研究普遍认为,控制地下水中砷迁移转化的机制主要有两种。第一种也是最重要的一个机制是:在还原性条件下,由微生物介导的铁氧化物还原溶解而导致吸附在矿物上的砷进入地下水中。这一过程主要发生在富含有机碳的年轻的冲积沉积物和三角洲沉积物中。在这个过程中,无机砷主要以As3+的形式进入地下水。第二种机制是:在高pH/氧化型含水层中,砷从铁氧化物表面解吸附,结果导致固体矿物中的砷释放,从而造成地下水的砷污染。这一过程主要发生在干旱和半干旱地区的内陆盆地或封闭式盆地,且无机砷主要以As5+的形式进入地下水。许多研究表明山西省盆地地区的地下水主要来自第四纪含水层,并且是这些地区最主要的饮用水来源。而除了一些富含碳酸氢钠的地下水外,这些地区地下水中高浓度的砷几乎都与还原性条件有关。　　 由于资金有限,当地政府只能将工作的重点放在已知的饮水型地方性砷中毒病村(地砷病病村)及其周边的村落。其他远离地砷病病村的地区也可能存在地下水砷污染,但是目前还没有被发现。在本研究中,我们的目的是建立一个基于Logistic回归的统计模型来预测山西省地下水砷污染地区的位置。首先,以现场采集的水砷数据作为因变量,以一系列环境解释变量作为自变量,分别建立以我国饮水砷浓度标准10μg/L和50μg/L为阈值的地下水砷污染预测模型;其次,依据建立的Logistic回归模型,结合相关的环境解释变量预测山西省地下水砷浓度大于10μg/L和50μg/L的概率,并建立地下水砷污染的概率地图和二元风险地图;最后,对模型的预测结果进行验证。　　 材料与方法:　　 1、研究地区　　 山西省位于我国北部,其地理坐标为北纬34°36’-40°44’,东经110°15’-114°32’之间。全省总面积约为15.6万平方千米,总人口大约为3374.60万。该省属于暖温带、温带大陆性气候,四季分明。年平均气温在-4-14℃,年平均降水量约为400-650mm。该省境内地形较为复杂,山地、丘陵、高原、盆地、平原等交错分布,且以山地和丘陵为主。山西省整个地貌是被黄土广泛覆盖的山地型高原,但从北到南分布着几个独立的盆地,分别是大同盆地、忻州盆地、太原盆地、长治盆地、临汾盆地和运城盆地。该省境内的平均海拔为1450m,最高峰是位于东北部的五台山(海拔为3058m)。　　 2、地下水砷浓度数据库　　 (1)现场采集的地下水砷浓度数据库:现场采集的地下水砷浓度数据来源于山西省2005-2007年中央转移支付项目饮水型地方性砷中毒调查。该调查采用“10%抽样法,,共筛查山西省959个村共5682口井。以村为单位,计算每个村的采样总井数、小于和大于10μg/L、50μg/L、100μg/L、200μg/L的井数、每个村的最高和最低水砷浓度。在GoogleEarth中对每个村庄进行地理空间定位,从而获取了712个村庄的经纬度数据。　　 (2)文献检索来源的地下水砷数据库:这些水砷数据来源于已发表的文献,其中删除了大同盆地的两个泉水采样点和运城盆地的一个水砷含量异常高的采样点(原文作者认为是由人为因素导致的污染),最终获得两个盆地共87个地下水采样点数据。这个数据集用于模型的验证使用。　　 3、饮水型地方性砷中毒病村数据收集与整理　　 在本研究中,地砷病病村按数据来源分为两部分。第一部分来源于山西省地方病防治所,其中包括应县的18个地砷病病村。第二部分来源于已发表的文献,其中包括山阴县的35个地砷病病村,朔城区的6个地砷病病村和平遥县的2个地砷病病村。这个数据库包含61个地砷病病村的村名和砷中毒患者检出率。我们以村为单位通过GoogleEarth对每个村庄进行地理空间定位,从而获取该村的经纬度数据。这个数据集用于分析地砷病病村与模型预测结果的空间分布关系。　　 4、环境解释变量　　 (1)地形相关地图:以空间分辨率为500m的数字高程模型为基础,在SAGA-GIS中将数字高程图的分辨率由500m增加到1km,并在地形分析模块计算生成以下与地形有关的地图:(1)地形湿度指数图,(2)集水区坡度图,(3)地形指数图,(4)修正的集水区图,(5)坡度图和(6)海拔图。　　 (2)遥感卫星图:从美国国家航空航天局下载覆盖整个研究区域连续两年的、空间分辨率为250m、8天合成的中分辨率成像光谱仪增强植被指数时间序列遥感影像。在ArcGIS中将影像分辨率从250m增加到1km,并应用主成分分析法对原始数据进行简化,取前8个主成分作为模型建立的辅助变量。　　 (3)全新世地质图:全新世地质与地下水高砷密切相关。在本研究中,我们根据山西省1:50万地质图,在ArcGIS中建立新的全新世地质图,并将其转化为空间分辨率为1km的栅格地图,从而建立二分类的山西省全新世地质图。　　 (4)土壤特征图:土壤相关信息来源于联合国粮农组织的土壤数据库。它包含了我国第二次全国土壤调查研究集成的1:100万土壤信息。我们选择一些可能与地下水砷污染有关的土壤特性作为模型建立的辅助变量,并根据以下规则建立空间分辨率为1km的二分类栅格地图:(1)土壤pH值是否高于7.0;(2)土壤是否显示盐碱土性质;(3)土壤中是否含有粘壤土,砂质粘土和壤土质土壤;(4)土壤电导率是否大于1dS/m;(5)土壤是否与冲积沉积物、湖相沉积物相关。　　 (5)水文学信息图:在中国的河流网数据库中,提取山西省的河流网数据。在ArcGIS中生成空间分辨率为1km的河流密度和河流距离栅格地图。这两个水文学特征可能会影响到地下水中砷的浓度。　　 (6)重力图:地球上的重力并不是均匀分布的。影响这一基本力量的因素主要有山脉和地表不同沉积物的密度差异。重力图常被用来估计沉积物厚度以及绘制地表以下的地图。由于在地质图中缺乏沉积物的厚度信息,因此本研究中,我们使用由丹麦科技大学Dr.OleBaltazar提供的空间分辨率为1km的重力图来提高和改进地质相关信息。　　 5、数据分析与统计　　 本研究中所有的数据统计分析都是在R软件中进行的。　　 (1)主成分分析:用来鉴别多维数据中各个变量的聚集模式。在本研究中,PCA用来探索辅助变量和采样点砷含量之间的关系。　　 (2)Logistic回归分析:Logistic回归被用来分析二元反应变量(高/低砷浓度)和环境解释变量之间的关系,从而建立一个以一定砷浓度为阈值的地下水砷污染风险概率预测模型。该模型如下:P[y=1|x]=[Exp(β0+β1X1+β2X2+...+βnXn]/[1+Exp(β0+β1X1+β2X2+...+βnXn)]　　 其中P是地下水砷浓度高于10μg/L或50μg/L发生的概率,y为因变量,β0是回归方程的截距,β1,β2,...,和βn是辅助变量X1,X2,...,和Xn的加权系数。我们使用后退逐步回归法,在将所有变量引入模型后,逐次把p值最大且p＞0.05的变量剔除出模型。每次剔除一个变量后重新拟合模型,按照上述剔除标准继续剔除变量,直至模型中所有变量的p＜0.05,从而建立最佳的回归模型。　　 (3)接受者操作特征曲线:接受者操作特征曲线是反映灵敏度和特异度连续变量的综合指标,是用构图法揭示灵敏度和特异度的相互关系。通过将连续变量设定出多个不同的分界点,从而计算出一系列灵敏度和特异度,再以灵敏度为纵坐标、1-特异度为横坐标绘制成曲线,曲线下面积越大,判断的准确性越高。　　 (4)地砷病病村位置与模型预测结果的空间分析:根据Logistic回归模型预测所得的二元地图,将地砷病病村与山西省地下水砷污染二元风险地图进行叠加分析,用以评价二者的空间分布关系及模型预测的准确程度。　　 结果:　　 1、水砷浓度数据描述　　 712个采样点主要分布在山西省从北到南的几个盆地地区,其中大同盆地(n=226),忻州盆地(n=31),太原盆地(n=320),临汾盆地(n=37),运城盆地(n=65),其他地区(n=33)。采样点水砷浓度范围变化较大,从＜10μg/L到500μg/L,平均值为24μg/L。其中22.75%(n=162)的采样点水砷含量超过了10μg/L且多分布在几个盆地;5.62%(n=40)的采样点水砷含量超过了50μg/L,水砷浓度最高的采样点出现在大同盆地。不同浓度采样点的空间分布表明水砷浓度较高的采样点主要位于盆地的中心位置且分布在河道的附近,而水砷浓度较低的采样点则远离盆地的中心位置和河道,多位于盆地边缘的山麓地带。大同盆地的桑干河、黄水河和南阳河流域,太原盆地的汾河流域是山西省境内受地下水砷污染影响最为严重的地区。　　 2、主成分分析　　 主成分分析结果显示有一组环境解释变量与第一主成分呈正相关关系,它们包括地形湿度指数、全新世沉积物、地形指数、修正的集水区、土壤结构、土壤盐分、以及部分遥感卫星图(EVI2,EVI4和EVI8)。该主成分主要代表了环境中的还原性条件。712个采样点的散点得分图显示地下水高砷与第一主成分呈正相关关系,表明山西省境内的高砷地下水与环境中的还原性条件密切相关。　　 3、Logistic回归模型　　 (1)10μg/L模型:23个环境解释变量在经过后退逐步回归法计算后,有7个变量(P＜0.05)进入了以10μg/L为阈值的Logistic回归模型中,这些变量分别是EVI1、EVI6、EVI7、地形湿度指数、重力、河流距离和全新世地质。其中EVI7、地形湿度指数、全新世地质与地下水高砷呈正相关(β值＞0),表示随着这些变量值的增加,地下水中砷浓度超过10μg/L的概率会增大;相反,EVI1、EVI6、重力、河流距离与地下水高砷呈负相关(β值＞0),表示随着这些变量值的增加,地下水中砷浓度超过10μg/L的概率会降低。　　 (2)50μg/L模型:50μg/L模型中包含6个环境解释变量,分别代表了地形特征、水文学条件、重力和土壤特性。其中地形相关的环境解释变量(地形湿度指数,β=22.93和地形指数,β=5.14)和土壤特性相关变量(土壤盐分,β=1.50)与高砷地下水呈正相关,表示随着这些变量值的增加,地下水中砷浓度大于50μg/L的概率会增大。而对于水文学变量(河流距离,β=-9.27)、重力(重力,β=-12.44)和遥感影像(EVI6,β=-4.78)与高砷地下水呈负相关,表示随着这些变量值的增加,地下水中砷浓度大于50μg/L的概率会降低。　　 4、地下水砷污染风险预测概率地图　　 10μg/L模型预测显示的概率地图表明临汾盆地北部的临汾市和洪洞县交界处,发生地下水砷浓度超过10μg/L的概率最高(＞0.8)。而50μg/L模型预测显示的概率地图表明忻州盆地的部分地区发生地下水砷浓度超过50μg/L的概率最高(＞0.8)。两个模型的预测结果都显示在大同盆地、忻州盆地、太原盆地和运城盆地发生地下水砷污染的概率较高(＞0.5),而东西两侧的山区和丘陵地带发生地下水砷污染的概率较低。　　 5、接受者操作特征曲线　　 (1)10μg/L模型:ROC曲线计算所得的10μg/L模型预测结果的最佳分界点概率值为0.22。　　 (2)50μg/L模型:ROC曲线计算所得的50μg/L模型预测结果的最佳分界点概率值为0.09。　　 6、地下水砷污染风险二元地图　　 (1)10μg/L模型:根据ROC曲线计算所得的0.22为最佳分界点概率值,将山西省地下水砷污染概率地图重新划分为一个二分类地图,从而建立山西省地下水砷污染发生的二元风险地图。结果显示山西省从北到南的若干盆地地区存在着大面积的地下水下砷污染高风险区域(水砷浓度大于10μg/L),主要分布在山西省境内的30多个县,共涉及8112km2的地区。　　 (2)50μg/L模型:根据ROC曲线计算所得的0.09为最佳分界点概率值,建立山西省地下水砷污染发生的二元风险地图。结果显示山西省有38的市县存在着面积大小不一的地下水砷污染高风险区域(水砷浓度大于50μg/L),共涉及3364km2的区域。　　 7、模型的验证结果　　 (1)10μg/L模型:20%数据集对模型进行内部验证的总正确率为68.3%,其中灵敏度为64.2%,特异度为81.8%。87个文献来源的采样点对模型进行外部验证的总正确率为63.2%,其中灵敏度为69.2%,特异度为60.7%。　　 (2)50μg/L模型:20%数据集对模型进行内部验证的总正确率为85.2%,其中灵敏度为72.7%,特异度为86.3%。87个文献来源的采样点对模型进行外部验证的总正确率为86.2%,其中灵敏度为66.7%,特异度为89.3%。　　 8、饮水型地方性砷中毒病村与模型预测结果的空间分布关系　　 地砷病病村与模型预测显示的概率地图叠加结果表明这些病村所在的位置其地下水砷污染发生的概率均较高。地砷病病村与山西省地下水砷污染风险二元地图的之间的空间分布关系表明大部分的病村都位于地下水砷污染发生的高风险区域。山西省61个地砷病病村与该省10μg/L模型预测显示的地下水砷污染风险二元地图的叠加结果显示,有90.2%(n=55)的地砷病病村落在预测显示的地下水砷污染高风险区域,9.8%(n=6)的地砷病病村落在预测显示的地下水砷污染低风险区域。61个地砷病病村与50μg/L模型预测显示的地下水砷污染风险二元地图的叠加结果显示,有91.8%(n=56)的地砷病病村落在预测显示的地下水砷污染高风险区域,8.2%(n=5)的地砷病病村落在预测显示的地下水砷污染低风险区域。　　 结论:　　 1、全新世沉积物图、地形湿度指数图、地形图、河流距离图、土壤盐分图、重力图和一些遥感卫星图在山西省地下水砷污染风险回归模型建立中起重要作用,有助于解释高砷地下水的分布。　　 2、模型显示山西省境内的高砷地下水与环境中的还原性条件密切相关,且干旱气候条件下的蒸发浓缩作用可能进一步增加了地下水中的砷含量。　　 3、该Logistic回归模型的建立,可以提供一种快速发现潜在地下水高砷区和饮水型砷中毒病区的方法,为饮水型地方性砷中毒的防治和未来的改水降砷工作提供理论依据和指导。　　 4、该Logistic回归模型建立和应用的方法可以被广泛应用于其他地区的地下水砷污染预测。同时该方法也可为其它地质成因污染物(如氟和硒)的模拟和预测提供借鉴。