基于gis的日遗化武含砷污染地统计学分析

基于gis的日遗化武含砷污染地统计学分析

日本侵略战争失败后，大量化学武器被掩埋并投入河流、湖泊和海洋。自2006年9月以来，在中国19个省市的近100个地点发现了日本遗产化学武器（日本文物武器）。日遗化武中装填的毒剂在日军秘密销毁和长期的埋藏过程中及“禁止化学武器公约”生效前简易处理过程中,由于泄漏和爆炸、焚烧销毁等原因,造成了部分土壤污染。至今,发现的日遗化武污染较集中的区域之一在吉林省某地。该地被污染的原因,据推测可能是日本投降时将化学武器就地销毁和掩埋,或其它时期对收集的日遗化武进行的简易销毁。由于当时技术条件有限,销毁时采取每隔一定的距离,挖掘一个深约1.0m、形状不定的销毁坑,放置数枚到数百枚的炮弹,利用常规弹药引爆,然后将销毁坑和未完全爆炸的炮弹掩埋。日遗化武大部分被引爆或炸毁,装填的含砷毒剂毒物也被分散到坑内或周边土壤中,并在长期的自然作用和人为扰动下,形成特殊的砷污染场地。自2004年该地日遗化武被发现后,通过对土壤中砷形态的分析,确认了该地的砷污染来源于日遗化武含砷毒剂的泄漏,并查明了该地砷污染的大致范围和污染程度,但由于销毁方式的特殊性和砷在土壤中的长期迁移及农耕、挖掘回收日遗化武等人为扰动,致使该地砷含量分布不同于自然形成的含砷土壤、含砷矿渣或工业源造成污染。为治理修复该地的砷污染,需了解和掌握该地砷分布的规律和特点。本文采用地统计学(Geostatisties),结合地理信息系统(GIS)技术,对该地砷污染的空间结构和分布进行解析与研究。地统计学现已被证明是分析土壤重金属空间分布特征及其变异规律最为有效的方法之一,是目前运用最为广泛的空间分析方法,有关这方面的国内外文献较多,但研究的对象多为污染程度相对较低的区域,母质对土壤重金属的贡献大于或接近于外界污染源的贡献,而销毁日遗化武含砷毒剂污染场地的污染状况、污染源分布等条件具有特殊性,在很多方面不同于上述研究对象的特征,导致土壤重金属空间变异特征呈现出不同于上述研究的规律,这方面的研究还是空白。另外,前人研究时,多采用简单随机采样方法,对于土壤污染物含量变化较大的污染区域,样点的代表性不高、采样密度不足,导致样点间的变异函数值的随机成分增加,小尺度结构特征被掩盖。本研究中的调查区域面积为26.19hm2,多由山坡和谷地组成,地形较为复杂,若按《HJ/T166-2004土壤环境监测技术规范》由变异系数和相对偏差计算的样品数为132个,但前人研究表明,地形复杂地区随取样密度增加,插值精度提高较快,增加取样数量可以较大幅度提高插值精度。为较好的反应该地砷污染的分布规律,本研究中,结合《HJ/T166-2004土壤环境监测技术规范》对污染事故监测土壤采样的要求,最终确定采取10m×10m的网格布点的方法取样,表层取样2668个,亚层取样637个样。1材料和方法1.1阳面坡面土壤条件该地在山区,纬度43.59°~43.62°,经度128.60°~128.62°。调查区域,在东南-西北走向的山麓阳面坡面和下方谷地,地势大体上东北高西南低;属于冷湿高山气候区,主要是花岗岩、玄武岩的风化物和坡积物形成的暗棕壤,呈弱酸性。在调查时,分为农田Ⅰ、Ⅱ区(17.78hm2)和林地(8.41hm2)三部分(图1),农田原为优质大豆种子基地,林地为人工种植杉木林。1.2样品的采集和测定1.2.1平面分布样品的采集与测定根据调查区域的污染形成原因、成土厚度及地下水位、当地农耕及农作物根系深度的事先调查结果,取样时按10m×10m的网格进行布点,网格交点即作为土壤调查取样点,取样深度取0.5m深度。采样时,先采集表层(-0.50~0m)混合土壤样品,四分法取分析样品约1.5kg,风干后,研碎,过100目筛,进行总砷测定;如果砷含量超过30mgkg-1(《GB19062-2003销毁日本遗弃在华化学武器土壤污染控制标准(试行)》),则在超标点对应位置,采集亚层(-1.00~-0.50m)混合土壤样品测定。表层采集样品2668个,亚层采集样品637个。因研究区域尺度较小,样点密集,不适合使用地理坐标,因此以纵向(北-南)0~800m,横向(西-东)0~840m建立相对坐标系,每个采样点位置采用相对坐标。1.2.2土壤剖面取样在具有代表性意义的农田地区,按HJ/T166-2004《土壤环境监测技术规范》的6.2.3.1部分采集剖面土壤样品。依据污染区范围和地势,以东西方向50m、南北方向20m的间隔,设置9个取样点(因空间不足,未在东北角设置取样点),并将每个土壤取样点从上至下将0~1.00m深度等距分层10层,各层内分别用土铲切取土壤然后集中起来混合均匀。1.3测定样品及方法所采集的土壤样品捏碎并铺于塑料盘上自然风干。土样风干后用玛瑙研钵研碎,过100目筛。土壤样品测试方法:利用硫酸-硝酸-高氯酸消煮,二乙基二硫代氨基甲酸银分光光度法测定,FI-HG-AAS法对测定结果进行验证。测定过程中采用GBW-07403土壤标准样品进行质量控制。1.4区间[x3以外的数据的消除由于特异值的存在会对变异函数具有显著的影响,因此计算变异函数前剔除这些特异值是十分有必要的。本研究在数据处理过程中采用域法识别特异值方法,即样本平均值加减3倍标准差,在此区间[x±3s以外的数据定为特异值,然后分别用正常的最大和最小值代替特异值,表层消除的特异值约占样本量的5.4%,表层消除的特异值约占样本量的1.8%。由于半方差函数的模型要求数据呈正态分布,否则会存在比例效应。消除特异值后,表层、亚层砷含量对数转换后数值经U检验都符合或近似符合正态分布。因此在半方差分析时,将表、亚层砷含量度进行对数转换,连同采样点的相对地理坐标输入地统计软件Arcgis10.0拟合半方差函数,并选择最佳拟合模型及其参数,用砷含量进行插值,得到研究区内土壤中各重金属含量的空间分布图。2结果与分析2.1土壤砷含量的测定对表层和亚层砷含量的测量值原始数据和原始数据消除特异值后经对数转化后的数据进行统计,结果见表1。由表1可知,从偏度系数和峰度系数来看,原始数据不符合正态分布;该地表、亚层砷含量原始数据变异系数分别为1.76.和1.52,均达到强变异程度;表层土壤中砷含量范围为0.09~758.44mgkg-1,亚层砷含量范围为0.23~1372.87mgkg-1,表层有57.38%、亚层有85.56%的样品超过该地背景值上限值;表、亚层砷含量的平均值分别为30.28mgkg-1、52.03mgkg-1,超过《土壤环境质量标准》(GB15618-1995)的国家二级标准和《销毁日本遗弃在华化学武器土壤污染控制标准(试行)》(GB19062-2003)规定的30mgkg-1。这些结果说明由于人为污染源的影响,使得该区某些网格中重金属元素的含量显著增加,且大量采样点砷含量超过国家二级标准,可能影响农业生产、人体健康。2.2空间相关因素2.2.1表层土壤砷含量空间分析特征统计分析只能粗略反映该地土壤砷含量情况,而不能反映该地砷含量的结构性和随机性,因此借助地统计学的半变异函数来进一步解释研究区土壤砷含量的空间结构性和随机性。半变异函数图通常可以被某些曲线方程所拟合,该拟合曲线方程就称为半变异函数的理论模型。半变异函数理论模型有球状(spherical)、高斯(Gaussian)、指数(Exponential)和线性(Linear)等。(1)半变异函数理论模型的确定以R2最大,RSS最小为原则选取的半变异函数模型为最佳模型[14~16]。在表4中,指数模型决定系数为0.913,在四种模型中最大;残差平方和为0.0528,在四种模型中最小,因此研究区域土壤砷含量半变异函数的最佳模型为指数模型。(2)指数模型参数分析半变异函数在原点处的数值称为块金值(C0),通常表示由于实验采样尺度大于空间变化尺度所引起的变异,属于随机变异。随机因素包括如施肥、农耕、种植、外在污染源等人为干扰。基台值通常表示系统内的总变异,它是结构性变异和随机性变异之和。结构性因素包括如气候、母质、地形、土壤类型等。理论上认为,在无人为扰动的自然条件下,土壤重金属半变异函数的块金效应为零。由表2可知,指数模型块金值为0.5957,具有明显的块金效应,而基台值(C0+C)为1.2983,表明研究区域内表层土壤除自身母质外还受到明显的人为干扰。另外,可能还存在采样误差、测定和短距离的变异引起的随机误差。结合研究区域实际,这种人为干扰是爆炸销毁时含砷毒剂散落到土壤中和当地在2004年被发现埋藏有日遗化武前,施肥、农耕、种植等活动所致。块金值与基台值之比(C0/(C0+C))是反映区域化变量空间异质性程度的重要指标,表示随机部分引起的空间异质性占系统总变异的比例,可作为衡量变量空间相关程度的尺度。按照区域化变量相关程度的分级标准,如果其值小于25%,表明变量的空间变异以结构性变异为主,具有较高的空间相关性;介于25%~75%之间为中等程度空间相关;大于75%时,以随机变量为主,则为空间弱相关。由表2可知,块金值/基台值(C0/(C0+C))为0.458,属于中等程度空间相关。研究区域内土壤砷的空间变异既有内在因子作用(或称为结构因素,如气候、地形、土壤类型等)外,一些外在因子(或称为随机因素,如农耕,施肥,种植等)也起了一定的作用,这些人为活动使空间相关性减弱。此外,该地区是日遗化武挖掘回收作业现场,挖掘回收活动也会使土壤中砷的空间相关性减弱。变程(A)表示随机变量在空间上的自相关尺度,它与观测以及取样尺度上影响土壤性状的各种过程的相互作用有关。在变程范围内,变量有空间自相关性,反之则不存在。由表2可知,有效变程为280.16m,表明在距离样点280.16m以内,砷含量具有空间相关性,反之则没有。(3)土壤砷含量最优空间插值模型分析对克里格插值模型进行选择优化能进一步准确而直观地剖析土壤砷含量的空间分布与变异特征。在确定最佳半变异函数的基础上,利用Arcmap中的Geostatistical扩展模块进行普通克里格、简单克利格、泛克利格三种空间插值模型插值,并得到各种插值后的交叉验证参数数值,结果见表3。根据克里格输出预测表面所生成的交叉验证参数评定标准:ME、MSE较接近于0,RMSE和ASE尽可能的接近,RMSSE较接近于1的插值方法精度相对较高,为最佳插值模型。比较表3的各种插值结果:简单克里格插值法的ME、MSE为1.773和-0.0541,比起普通克里格法(ME、MSE分别为2.864、0.0203)以及泛克里格法(ME、MSE分别为-2.172、-0.2187)更接近于0;RMSE、MSE分别为57.48、61.01,与普通克里格法(RMSE、MSE分别为54.06、94.94)以及泛克里格法(RMSE、MSE分别为54.95、68.89)相比,两者的值更加接近;RMSSE值为1.015比起0.6081和1.439也更加接近于1。综上所述,该地土壤砷的空间分布最佳插值效果是简单克里格法。2.2.2亚层土壤砷含量空间分析采取相同方法和步骤对亚层土壤中砷进行相关分析,并与表层土壤进行比较。(1)各最佳模型及参数在表1特征统计值计算的基础上,进行土壤亚层砷的最佳半变异函数及最佳克里格插值的确定,结果如表4。从表4看出,研究区域亚层土壤砷空间分析的最佳半变异函数为高斯模型,最佳克里格插值为简单克里格。2.3日遗化武从农田和林地地位分布的角度分析2.3.1表层土壤砷含量的空间插值利用前文分析的最优克里格插值模型对研究区域表层土壤砷的空间插值,并对插值图进行解释和进一步分析。在ARCmap中的Geostatistical扩展模块上进行简单克立格插值,进而得到研究区域土壤中砷含量的空间分布图(图2)。由图2看出,砷含量的空间分布具有明显的地域特征:(1)研究区域内高砷含量集中在农田与林地交汇的中部地带,并在此处形成了较大面积的高砷区域(砷含量为30~300mgkg-1),约占研究区域总面积的1/5。可推测该高含量区域为当年销毁日遗化武的中心地域,同时也是未被完全引爆的日遗化武集中埋藏地域。经爆炸散落和埋藏期间泄露的含砷毒剂,造成该区域的砷含量明显高于其它区域。(2)砷含量以农田林地交汇处为中心,向西北和东南方向出现递减趋势。这一规律,一方面与销毁时炮弹销毁坑的分布有关;另一方面可能也与研究区域的地理位置有关,研究区域是一片东南-西北走向的带状谷地,南北都是山丘,风向导致泄漏和销毁释放出来的含砷毒剂顺着东南-西北走向散落到土壤中并出现递减。另外,此处原为大豆种植区域,农耕时依据地势进行的东南-西北方向的耕地及打垄,将砷污染多次混合并扩大污染区域,使变化趋势更加明显。(3)该区域内的西北角和东南角均有小范围内的高含量斑块。原因可能是:1东南角和西北角靠近道路,方便设置零星的炮弹坑;2西北角靠近居民区,当地居民将散落的炮弹捡回,重新抛弃或埋藏。2.3.2亚层砷含量的空间插值在与表层相同的坐标上,对亚层土壤中的砷含量进行简单克里格插值,插值效果图见图3。从图2、图3整体上看,表、亚层的空间分布格局类似。高浓度区都集中在评价区域中部农田和林地交汇处。不同的是土壤表层砷含量绝大多数都高于亚层,部分样点存在亚层土壤砷含量高于表层的现象。在较深处形成的砷含量超标的原因:1销毁时,销毁坑所在位置;2销毁后,炮弹坑内未爆炸和未完全爆炸的炮弹,在埋藏过程中泄漏;3挖掘过程对污染土壤的扰动。土壤表层砷含量y与亚层x呈弱正相关(y=0.226x+77.416,r=0.231),说明该地表层和亚层无明显的线性关系。2.4污染土壤深度在具有代表性意义的农田地区,按1.2.2取样测定,结果见图4。从图4可以看出,该地农田内最大含砷量基本上都出现在地表(-20~0cm),大多数样点随着土壤深度的增加,土壤砷含量逐渐降低,并且降低幅度大,在垂直向下1.00m处的砷含量普遍很低,几乎都在5mgkg-1之下;砷含量超标的样品,几乎全部在地表-30cm以上,进一步表明除销毁坑外,该地污染土壤是由爆炸时含砷化合物散落造成的,污染的深度是爆炸散落、人为干扰和砷的自然迁移等原因共同影响的结果。砷含量的垂直分布结果,为修复该地农田时选择方法提供了基础的参考数据。3最佳插值模型及半变异函数本研究以吉林省某日遗化武污染场地为研究对象,在调查掌握土壤砷含量的基础上,利用地统计学方法探讨了土壤砷含量的半变异函数及最佳拟合模型,在地理信息系统(GIS)平台上通过对砷含量多种空间插值模型的比较分析,筛选出最佳插值模型,并分析得到土壤砷含量的空间分布特征:(1)研究区土壤表层(-0.50~0m)砷含量变异系数达1.76,具有明显的不均匀性。空间变异主要受销毁时的爆炸分散、埋藏的未爆含砷日遗化武的泄露、挖掘和农作物种植等因素影响,属中等程度的空间相关;其最佳半变异函数模型为指数模型,最佳插值模型为简单克里格插值;表层空间分布插值结果具有明显的地域特征。(2)土壤亚层