本发明属于区域环境污染及其公共健康风险的监测技术领域，特别是涉及一种面向基站分布稀疏区域精细定量飞灰有毒金属暴露人群的方法。
背景技术
飞灰是以煤为主的燃料经过高温燃烧后排放的一类工业副产品。2015年，全球飞灰年产量已达12亿吨。但全球平均实际回收利用率不足25％(yaoz.t.,jix.s.,sarκerp.k.,tangj.h.,gel.q.,xiam.s.,xiy.q.(2015)acomprehensivereviewontheapplicationsofcoalflyash.earth-sciencereviews,141:105–121.)。中国、印度、美国、俄罗斯、日本以及南非等国是这一工业副产品的主要排放国。飞灰平均直径小于20μm，涉及316种单一矿物和188种复合矿物，痕量金属含量为原煤的4-10倍，成分复杂。此类轻质颗粒物及其表面复杂的非均相反应在不利气象条件下会迅速恶化区域空气质量。目前，多国并未将其归于危险固废，疏于监管，由此而来的健康问题引人关注。迫切需要建立一种快速、高分辨率的方法有效监测区域飞灰有毒金属污染与其公共健康风险。然而，解答这一关键问题是目前极具挑战性的前沿领域之一。
近年，仪器在线测量颗粒物化学成分逐渐成为污染过程精细研究的重要手段。但受仪器设计原理的限制，只能测量部分粒径范围的颗粒物。这在反映颗粒物族群的代表性方面仍有较大局限，并对某些源类具有特征指示性的关键物种(例如，扬尘、飞灰、土尘等)难以实现在线监测(高健,李慧,史国良,丁爱军,游志强,张岳翀,王涵,柴发合,王淑兰.(2016)颗粒物动态源解析方法综述与应用展望.科学通报,61(27):3002-3021.)。同样，空气质量模拟可以获得不同行业排放来源时空贡献的多维分析。但受排放清单的限制，在很大程度上增加了扬尘、生物质燃烧、飞灰等关键物种源模拟结果的不确定性。颗粒物成份测量这种经典的方法是通过基于撞击原理、空气动力学原理和虚拟撞击原理有选择地进行颗粒物粒径分离，根据不同研究目的使用不同材质的采样膜收集所需粒径的颗粒物之后，在实验室进行离线分析。该类方法定量较为准确，但存在分析周期长、分析费用高、时间分辨率低以及样品化学成份极易被破坏等缺陷。
相对以上常用方法，样品磁性测定与化学成份相结合的方法为鉴别飞灰、烟尘、尘埃、大气悬浮粒、有毒金属来源以及标绘其传输、沉降等演变过程提供了一条新途径(thompsonr.,oldfieldf.environmentalmagnetism,allen&unwin,landon,1986.)。这种快速、经济、对样品无损的磁测技术在一定程度上也克服了因为添加工业催化剂而导致同位素示踪方法对痕量金属等污染物来源甄别不清的问题。
区域环境污染及其公共健康风险不仅是学术前沿，也是政府部门迫切需要解决的问题。与健康相关的行为和结果需经地理变量才能获得更为准确的解释。具有与健康相关的高时空精度信息为这一领域的研究与预测预警技术提供了可能，并在不断探索和扩展之中(shix.,kwanm-p.(2015)introduction:geospatialhealthresearchandgis.annalsofgis,21(2):1-3.)。目前，监测人群行为主要有两类数据源：社交性大数据和对地观测大数据。其中，个人移动电话等社交性大数据可以获得真实的人群短时空形迹。但网络覆盖密度的高低决定着该信息的准确性。由于涉及个人隐私等安全问题严重限制了数据的数量、范围以及时效，进而影响了该类信息的应用价值。
星载微光云图探测信息能够弥补此类局限。目前，星载微光影像有两类：美国国家航空航天局(nationalaeronauticsandspaceadministration,nasa)的国防气象卫星计划(defensemeteorologicalsatelliteprogram,dmsp)线性扫描业务系统(operationallinescansystem,ols)和2011年发射的新一代极轨运行环境卫星系统预备卫星计划(s-npp)可见光红外辐射仪(visibleinfraredimagingradiometersuite,viirs)上的白天/夜间波段(day/nightband,dnb)。dmsp-ols是近年广泛使用的人口数据产品的基础数据源。例如，thegriddedpopulationoftheworld、theglobalruralurbanmappingproject、landscan以及worldpop等产品。相对而言，s-npp/viirsdnb具有更高的时空分辨率(https://www.nasa.gov/mission_pages/npp/main/index.html)，并通过采用与viirs其他波段相一致的辐射校正，实现了以日为时间序列来定量探测地球夜间灯光、火光等辐射信号。同时，以更小的瞬时视场增强了空间图像的清晰度，多层灰度级削弱了像元的过饱和现象，辐射修订提高了定量观测精度，消除了像素间的跳变(caoc.,lucciaf.j.d.,xiongx.,wolfer.,wengf.(2013)earlyon-orbitperformanceofthevisibleinfraredimagingradiometersuiteonboardthesuominationalpolar-orbitingpartnership(s-npp)satellite.ieeetransactionsongeoscience&remotesensing,52(2):1142-1156.)。该数据源能够提供相对更为精细的人类活动的信息。其中所反映的综合信息也在城镇化、环境、灾害、渔业、能源、社会经济发展等领域显示了广阔的应用前景。
基于以上认识，本发明提出采用磁学信息示踪飞灰来源，通过多点同位协同模拟方法，提取人体健康风险下飞灰有毒金属分布区域，进一步采用降尺度算法将基于s-npp/viirsdnb辐射度反演的人口空间信息与超过人体健康风险阈值的飞灰源有毒金属含量集成于每一个栅格中，快速获得区域飞灰有毒金属暴露风险人群精细的空间表征，提高了基站分布稀疏/无覆盖区域人群的空间分析精度。
技术实现要素：
为达到快速获得区域飞灰有毒金属暴露风险人群精细的空间表征从而提高基站分布稀疏/无覆盖区域人群的空间分析精度的目的，本发明提供一种精细定量基站稀疏区域飞灰有毒金属暴露人群的方法，解决了现有的颗粒物监测技术中对某些类源具有特征指示性的关键物种监测的局限性问题，为基站分布稀疏/无覆盖区域提供一种快速、准确、精细、集成的飞灰有毒金属暴露人群格网空间计算方法。
本发明所采用的技术方案：本发明提供了一种精细定量基站稀疏区域飞灰有毒金属暴露人群的方法，包括以下步骤：
步骤1)：采集样品；任意选取一个研究区域，依照中华人民共和国国家标准城市用的分类与规划建设用地标准(gb50137-2011)中城乡用地14个小类布设采样点；打开地理信息系统软件arcgis10.2，导入美国陆地卫星landsat-8oli影像，将手持卫星定位系统(gps)与arcgis10.2中的手持卫星定位系统连接设置模块(gpsconnectionsetup)进行连接；平均采样网格间距为2.5km×2.5km，现场使用便携式磁化率仪(sm-30)测定每一采样点的体积磁化率κ，用作飞灰来源的示踪参数，每一元素样品由每一采样点1平方米范围内对角采集的4份样品混合而成；
步骤2)：测定样品元素含量，含量测定为有毒金属全量或有毒金属有效态含量；所述有毒金属为国家优先控制的铅、汞、镉、铬和类金属砷；
有毒金属全量的测定：铅、镉、铬和类金属砷的全量采用王水微波消解法进行样品处理，采用电感耦合等离子质谱仪(icp-ms)测定提取液中各有毒金属全量，单位为mg·kg-1；汞无需消解处理样品，直接采用ra-915m塞曼效应汞分析仪与配套pyro-915热解装置快速测定总汞含量，单位为ng·g-1；
有毒金属形态含量的测定：铅、镉、铬和类金属砷的形态含量利用连续萃取技术浸提获得，通过icp-ms测定，单位为mg·kg-1；甲基汞的含量通过萃取-乙基化结合气相色谱-冷原子荧光联用法分离测定，采用merx全自动烷基汞/总汞分析系统测定，测得结果由积分仪记录，通过峰面积计算其含量，单位为ng·g-1；甲基汞作为汞的有机形态，只要检出就有毒性风险，和其它金属的有效态的意义是一样的；
步骤3)：将步骤1)中获得的每一个目标采样点的体积磁化率κ和步骤2)中获得的有毒金属全量或有效态含量利用序贯指示同位随机协模拟算法(sics)进行同位协模拟，设定超过国家eqssgradeⅱ有毒金属限值的全量(gb15618-1995)或大于0mg·kg-1的有毒金属有效态含量作为协模拟过程中的硬指示变量，所述有效态是以bcr法的重金属形态分类为依据的碳酸盐结合态；体积磁化率κ是飞灰磁性示踪参数，作为协变量(软数据)，确定人体健康风险下飞灰有毒金属含量空间模式及其模拟预测的不确定性；本文确定人体健康风险下飞灰有毒金属含量空间模式涉及到了，空间模式的表达，即地理空间与分析空间、数据集与可视化、空间数据的标度及其存储类型；空间模式的转换，即分布中心计算与插值分析与抽样；和空间模式的测度，即空间随机过程：
序贯指示同位随机协模拟算法主要步骤：
1)指示编码；选取w个门槛值z1，…,zw,有毒金属全量或有效态含量z(sα)(α＝1,…,n)为离散化区域变量，设其局部硬数据i(sα；zw)，定义该数据的指示随机变量函数，
i(sα；zw)＝{1ifz(sα)≤zw；0otherwise}；(公式-1)
其中，α指一系列采样点s的位置；n指一系列采样点s中第n个采样点s的位置；sα指第n点的某一表土有毒金属全量或有效态；z(sα)指第n点的某一表土有毒金属全量或有效态的含量值；zw指位于z(sα)区间内的被设定的离散的阈值；即当第n点的某一表土重金属全量小于或等于设定的门槛值zw，通过指示随机变量函数i(sα；zw)将有毒金属全量转换为指示编码1，否则执行0；或者当有毒金属有效态含量小于等于设定的门槛值zw，即有毒金属有效态含量0mg·kg-1，通过指示随机变量函数i(sα；zw)将有效态含量z(sα)转换为指示编码1，否则执行0；
2)将飞灰磁性示踪参数κ作为协同变量，设其局部软数据u(sα；zw)，该函数为z(sα)提供条件先验概率函数：
u(sα；zw)＝prob{z(sα)≤zw|u(sα)}；(公式-2)
其中，u(sα)为第n点的飞灰磁性示踪参数κ的值；
3)估算全局条件概率分布值：
其中，f(zw)指表土有毒金属全量或有效态含量的条件先验概率函数；a指用于将协同变量飞灰磁性示踪参数κ进行离散化处理的界限值；z(s)指随机模拟中某一采样点s表土有毒金属全量或有效态含量值，z(s)是原来样点的z(sα)的实际值经过先验条件概率的计算过程后，变为的后验条件概率中的估算值，即利用该公式对条件先验概率函数prob{z(sα)≤zw|u(sα)}进行估算；
4)后验更新：
由相邻的多个局部先验累积分布函数更新至后验累积分布函数来实现全局样点s的不确定性建模，通过指示克里格法估计有毒金属全量或有效态含量的条件累积分布函数：
其中，prob是条件累积分布函数的英文缩写；z(s)指随机模拟中某一采样点s表土有毒金属全量或有效态含量值；zw指设定的门槛值，含义与步骤1)中zw相同；n指第n个采样点位置；n'指第n'个采样点位置；ik指指示克里格法；λ0指全局先验累积分布函数的加权因子；λα(s；zw)指与样点n相邻的硬数据的加权因子；vα'(s；zw)指与样点n'相邻的软数据的加权因子；u(s'α；zw)指样点n'硬数据的条件先验概率函数；i(sα；zw)指样点n有毒金属全量或有效态含量的指示随机变量函数；
5)空间相关性：
不同门槛值zw时，软数据与硬数据的空间相关性由马尔科夫-贝叶斯(marκov–bayes)算法可得：
其中，h为所有采样点位置对之间的距离向量；ci(h,zw)为硬数据的指示协方差；cy(h,zw)为软数据的协方差；ciy(h,zw)为硬数据与软数据的交互协方差；b(zw)为系数，由软数据对硬数据的校正而得；
序贯指示同位随机协模拟算法的实现方法：
将飞灰磁性示踪参数κ作为同位协同模拟变量，通过模拟数据集中α个采样点上的z(sα)和u(sα)，采用公式1-3，建立验证数据点位置z(sα)的先验概率；验证结果表明马尔可夫假设即公式-5成立后，基于公式-5计算不同门槛值zw对应的交互协方差；定义一条遍历所有待模拟位置的随机路径，通过以下步骤实现模拟：
(a)对于随机路径上某一采样点s，使用已经建立的指示协克里格估算该点有毒金属全量或有效态含量小于或等于既定该有毒金属变量门槛值的概率并获得相应的概率值；
(b)由已经获得的2个概率值外推建立该有毒金属变量在采样点s的先验条件累积分布函数，从累积分布函数随机获取一个模拟值并将该值赋予采样点s；
(c)将此模拟值的指示编码转入下一个未知的该有毒金属变量先验条件累积分布函数的建模中；
(d)沿设定访问每一个未采样点的随机路径重复(a)-(c)的步骤，遍历所有结点，每一个结点均参与一次模拟获得一个模拟值，完成一个序贯指示同位随机协模拟算法的实现；
(3)序贯指示同位随机协模拟算法模拟预测的空间不确定性评估：
重复m次序贯指示同位随机协模拟算法，每次模拟均使用不同的随机路径产生m个等概率的模拟实现；由贝叶斯公式更新模型参数概率分布估计得到的后验概率来表明其中的不确定性，计算公式为：
其中，m为随机模拟次数；zcl为人体健康风险下的有毒金属含量规定限值；n(i)为m次模拟中在点i的有毒金属变量模拟值z(i)大于阈值zcl的次数；
步骤4)：通过星载微光探测信息反演每一行政单元建成区内人群表征；基于s-npp/viirsdnb数据源，构建每个行政单元内微光探测信息与人口普查统计数之间的回归模型，估算每个行政单元上的人口数，通过每一个栅格上人口普查统计数与微光探测信息估算的总人口数的百分比进行校正并将校正后的人口数重新配置于街区尺度的空间格网上：所述s-npp/viirsdnb为新一代极轨运行环境卫星系统预备卫星计划可见光红外辐射仪的白天/夜间波段；
a)每个行政单元内s-npp/viirsdnb辐射度与人口普查统计数之间的回归模型如下：
其中，d指单一行政单元；
b)每一行政单元空间格网上人口的再配置，具体操作依照如下方程：
grid∑(d)pop(x,y)'＝grid∑(d)pop(x,y)×wd；(公式-9)
其中，grid指规则的空间格网单元；(x,y)指每一行政单元空间格网上任意一点的地理位置；wd指每个行政单元建成区内每一栅格单元中人口普查统计数与估算总人口数的百分比；
步骤5)：通过降尺度方法对步骤3)和步骤4)中所得结果进行计算，获得每一栅格处于飞灰有毒金属暴露风险中的人口密度，具体算法如下：
dpe(x,y)＝pop(x,y)'×e(x,y)；(公式-10)
其中，dpe(x,y)指每一栅格任一点超过人体健康风险门槛值的飞灰有毒金属暴露人口数；pop(x,y)'指每一栅格中任一点校正后的人口数；e(x,y)指每一栅格任一点人体健康风险下携带飞灰示踪信息的有毒金属含量。
本发明的有益效果是：
(1)本发明采用了新一代星载微光探测数据s-npp/viirsdnb，克服了社交性大数据的局限性，细化了基站分布稀疏/未覆盖区域城镇内部人口表征的空间粒度，实现了快速、集成、操作友好的人口密度精细格网定量表达。
(2)本发明结合遥感影像提取的人类活动建成区和街区人口普查统计数据，将s-npp/viirsdnb估算的人口数反演到空间分辨率为100×100米的空间格网中，数据质量好，容易获取，操作友好，格网单元的人口表征提高了城镇群内飞灰-有毒金属-人群空间准确集成分析精度；
(3)本发明采用便携式磁化率仪sm-30进行现场瞬时测量样点体积磁化率κ，可以清晰指印与燃煤相关工业排放来源的飞灰和有毒金属，克服了目前源解析(在线)技术和方法对此类有特征指示性的关键性物种来源诊断的局限性，能够作为一种燃煤飞灰原位快速、经济、对样品无损的微粒污染源示踪监测技术。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，以下描述中的附图仅为本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的实现流程图；
图2为本发明具体实例选取区内飞灰汞的交叉协变量图；
图3为本发明具体实例选取区内超过人体健康风险门槛值的飞灰汞序贯指示协同随机模拟图；
图4为本发明具体实例选取区内超过人体健康风险门槛值的飞灰汞序贯指示协同空间不确定性评估概率图；
图5为本发明具体实例选取区内每一行政单元建成区内s-npp/viirsdnb辐射度与人口普查统计数据之间的相关性图示；
图6为本发明具体实例选取区内校正后估算人口数的100m×100m格网反演栅格图；
图7为本发明具体实例选取区s-npp/viirsdnb辐射度反演人口表征准确性评估的蒙特卡洛模拟图；
图8为本发明具体实例选取区每一栅格超过人体健康风险门槛值的飞灰汞潜在暴露人口密度空间图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
图1为本发明实现面向基站稀疏区域飞灰有毒金属街区暴露人群定量方法的流程图。选取区为某亚洲中部干旱区绿洲城镇群，位于丝绸之路经济带核心区域，处于工业化中期初期阶段，以重工业为主导的”重型“产业结构，总面积约为1.42万平方公里(km2)，建成区面积约为500km2，常驻人口接近300万人，基站分布稀疏。境内石油、煤化工、钢铁、水泥、煤电、氯碱、冶金等重工业行业布局较为齐全，燃煤飞灰年产量超过1000万吨(t)，该选取区内仅氯碱工业每年排放含汞固废逾300t，释放大气汞约16-200t以及4t含汞废液。
第一部分：样品采集。实施中，依据gb50137-2011标准，将gps与arcgis10.2中导入的landsat8oli影像进行连接，实现2.5km×2.5km平均采样网格间距上的实时采样，采集选取区城乡用地14个小类表土(0-5cm)样品并记录目标样点经纬度坐标，1个目标样为1平方米范围内对角采集样混合而成，现场同时使用磁化率仪sm-30测定该目标样点的体积磁化率κ值，该选取区内共采集425个目标样样品，干样在实验室经冷冻干燥后，过200目筛，称取0.0010-0.2000g，采用lumex测汞仪对每个样品进行4次平行测定，快速获得其总汞(thg)含量；
第二部分：构建thg-飞灰暴露风险概率模型并评估其空间表征的不确定性。将第一部分的实验数据导入sgems或gslib软件包中执行sics模拟，每一次实现均为等概率模拟，在本实施例中，以国家eqssgradeⅱ限值(gb15618-1995)为标准，硬数据thg含量0.5mg·kg-1(ph＞6.5)作为人体健康风险门槛值zw，体积磁化率κ作为燃煤飞灰的示踪数据集(n＝351个)作为软数据执行同位协同模拟，351个采样点的thg值、κ值通过公式1-3建立验证数据点位置thg含量的先验概率，结果显示markov假设成立(公式-5),基于markov–bayes算法获得门槛值zw对应的κ的协方差模型以及thg与κ的交互协方差模型，建立二者的空间相关性(如图2所示)，对于随机路径上某一点z(sα)，使用上述基于markov–bayes算法已经建立的同位指示协克里格估算该点有毒金属全量或有效态含量小于等于既定该有毒金属变量门槛值的概率并获得相应的概率值，由已经获得的2个概率值外推建立该有毒金属变量在第n点的先验条件累积分布函数，从累积分布函数随机获取一个同位变量的模拟值并将该值赋予采样点n，将此模拟值的指示编码转入下一个位置n'先验条件累积分布函数的建模中，沿随机路径重复以上步骤遍历所有结点，使每个结点获得一个模拟值从而完成一个sics的实现，重复800次sics，每次模拟均使用不同的随机路径将产生800个等概率的模拟实现，以此提供超过人体健康风险门槛值的thg-飞灰空间协同模拟以及空间不确定性定量评估联合概率模型(如图3、图4所示)(具体步骤参见步骤3)；在概率值高于0.8的区域内，说明在800次模拟中超过其相应门槛值的次数超过了640次，而在概率值低于0.2的区域内，说明在600次模拟中超过其相应门槛值的次数不足160次，表明在这两类区域内的空间模拟结果的可靠性较高，其中概率值高于0.8的thg-飞灰分布区域可以作为人体健康风险优先管理热区；
第三部分：估算选取区每个行政单元建成区内人口数，在100m×100m网格上反演估算人口数并实现空间再配置。从美国海洋大气局/美国地学资料中心(noaa/ngdc)(http://ngdc.noaa.gov/eog/viirs/download_viirs_ntl.html)免费下载与外业采样同时段的s-npp/viirsdnb影像数据，并从美国地质勘探局(http://landsat.usgs.gov/)免费下载landsat8oli同期数据，每个行政单元内人口普查统计数据来自选取区统计年鉴。s-npp/viirsdnb为新一代极轨运行环境卫星系统预备卫星计划可见光红外辐射仪的白天/夜间波段。基于landsat8oli和选取区行政单元矢量图提取每一行政单元建成区，通过公式-7建立每个行政单元建成区内s-npp/viirsdnb辐射度与人口普查统计数据之间的相关性，回归模型显示二者存在高相关性，结果如图5所示，每个行政单元上估算的人口数由于线性回归的误差等原因，会使得每个行政单元内人口普查统计数据与空间格网上反演的估算人口数之间存在偏差，需对选取区内反演估算的总人口进行校正，依据公式-8、公式-9，通过对每一个栅格上人口普查统计数与反演估算总人口数的百分比进行校正，并将校正后的人口数重新配置在100m×100m网格上，实现反演估算人口数的空间再配置，结果如图6所示，蒙特卡洛模拟表明基于s-npp/viirsdnb辐射度反演人口分布在实施例350万总人口中的0.5-240万人之间的结果可靠，如图7所示；
第四部分：选取区人体健康风险门槛值约束下的thg-飞灰街区暴露人群定量表征。通过降尺度算法(公式-10)将校正后每一栅格的人口数图层与每一栅格任一点人体健康风险下携带燃煤飞灰磁信息的有毒金属含量图层集成并进行区域空间计算，获取每一栅格中超过人体健康风险门槛值的thg-飞灰暴露人群密度分布，如图8所示，该实施例经济、快速、精细的实现了这个丝绸之路经济带重要节点城镇群内thg-飞灰-人群准确集成的空间定量表征。
该实施例说明，通过多点同位协同随机模拟和降尺度的方法，集成磁测信息、新一代星载微光数据与有毒金属含量值提高了区域城镇内飞灰有毒金属暴露人群空间精度和可靠性。采用本发明的磁测方法示踪微粒污染源提高了甄别飞灰等某些源类中具有特征指示性关键性物种来源的准确性，现场、原位样品无损、快速的体积磁化率测定不仅提高了源解析的工作效率，也降低了调查费用。同时，新一代星载微光探测数据的应用对全球基站分布稀疏区域环境突发事件预警系统、环境流行病学及其社会经济等驱动要素的精细化时空分析与有效监测具有重要的理论意义及实用价值。
本说明书中的各个实施例均采用相关的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围。