2.3 污染物迁移数据处理与制图
2.3.1 离散数据网格化方法简介
一般实施的采样方法为不规则测网法。绘制污染等值线图时，要根据客观环境特征和数据本身的特点，选择合适的网格化方法。
（1）网格化方法的特征及应用条件 网格化是指通过一定的插值方法，将稀疏的、不规则分布的数据插值加密为规则分布的数据，以适合绘图的需要。原始数据的不规则分布，造成缺失数据的“空洞”。网格化则是用外推或者内插的算法填充了这些“空洞”。
（2）网格化方法 主要等值线图绘制软件有Surfer及Mapgis。
Mapgis中有4种网格化方法：距离幂函数反比加权网格化、Kring泛克里格法网格化、稠密数据中指选取网格化、稠密数据高斯距离权网格化。
Surfer中网格化方法有12种。Surfer网格化方法基本涵盖Mapgis网格化方法。Surfer网格化方法有：加权反距离插值法（inverse distance to a power）、克里格法（Kriging）、最小曲率法（minimum curvature）、谢别德法（Modified shepard’s method）、自然邻点插值法（natural neighbor）、最近邻点法（nearest neighbor）、多项式回归法（polynomial regression）、径向基本函数法（radial basis function）、带线性插值的三角剖分法（triangulation/liner interpolation）、移动平均法（moving average）、数据度量法（data metrics）和局部多项式法（local polynomial）。
Surfer应用更广泛，下面详细介绍Surfer中各种网格化方法的特征。
①加权反距离插值法 首先由气象学家和地质工作者提出的加权反距离法是一种加权平均内插，可以是准确插值或平滑插值，通常表现为准确插值。基本原理是设平面上分布一系列离散点，已知其位置坐标和属性值，P（x，y）为任一网格点，根据周围离散点的属性值，通过距离加权插值求P点属性值。实质是待插值点领域内已知散乱点属性值的加权平均，权的大小与待插值点领域内散乱点之间的距离有关，是距离n次方的倒数。
加权反距离插值法认为任何一个观测值都对邻近的区域有影响，且影响的大小随距离的增大而减小。在计算一个网格节点的Z值时，赋给数据点的权重是分数，所有数据点的权重和为1。权重与数据点到节点的距离成反比，愈靠近节点的原始数据点，其权重愈大。
该方法的优点是可以通过权重调整空间插值等值线的结构，但是其计算值容易受到数据点集群的影响，计算结果中常出现孤立点数据明显高于周围数据点的现象，表现为在网格区内围绕着某些数据点可能产生牛眼状（Bull’s eye）等值线。可以通过设置Smoothing参数平滑内插网格来消减牛眼效应。加权反距离插值法是一种非常快速的网格化方法，在小于500个数据点时，可以使用No Search（使用所有点）的搜索类型来快速生成网格。加权反距离插值法方程：
式中，hij为网格节点“j”与邻近点“i”之间的有效分离距离；Zj为网格节点“j”的内插值；Zi为邻近点；dij为网格节点“j”与邻近点“i”之间的距离；β为权重系数；δ为Smoothing参数。
加权反距离插值法高级选项：设置加权反距离插值高级选项Power和Smoothing参数。权重系数Power确定随着数据点到网格节点距离的增加，其权重降低的程度。随着Power逼近0，生成的表面逼近一个水平面，该平面通过数据文件中的所有观测点的平均值。随着权重系数的增加，生成的表面由最邻近点插值，导致表面变成多边形。多边形表现了最接近内插节点的观测表面。可接受的权重系数通常在1～3之间。
平滑参数Smoothing把“不确定性”因素与用户输入的数据联系起来，平滑参数愈大，计算相邻网格节点Z值时特异观察点的绝对影响就愈小。平滑参数大于0，则没有任何一个数据点对于某个节点的权重为1，即使该数据点正好位于网格节点上。
②克里格法 克里格法是一种在许多领域都很有用的地质统计网格化方法。最初是由南非金矿地质学家克里格根据南非金矿的具体情况提出的计算矿产储量的方法：按照样品与待估块段的相对空间位置和相关程度来计算块段品位及储量，并使估计误差为最小。后来，法国学者马特隆对克里格法进行了详细的研究，使之公式化和合理化。克里格法的基本原理是根据相邻变量的值（如若干样品元素含量值），利用变差函数所揭示的区域化变量的内在联系来估计空间变量数值。克里格法试图表示隐含在数据中的趋势，例如，高点会是沿一个脊连接，而不是被牛眼形等值线所孤立。
该方法总是尽可能地去描述原数据所隐含的趋势特征，以区域化变量理论为基础，以变差函数为主要工具，在保证研究对象的估计值满足无偏性条件和最小方差条件的前提下求得估计值。对于高值数据点会使之沿某一“脊”分布，而不围绕该点孤立插值，不形成牛眼状等值线。克里格法极为灵活，广泛地应用于各个科学领域，适于各种类型的离散数据，网格化精度高，是极佳的网格化方法。克里格法中包含了几个因子：变异图模型、漂移类型和矿块效应。其中变异图模型（Variogram Model）是用来确定插值每一个结点时所用数据点的邻域，以及在计算结点时给予数据点的权重。
Surfer提供了多种最常用的变异图模型，它们是指数、高斯模型、线性、对数、矿块效应、幂、二次模型、有理数二次模型、球面模型和波（空洞效应）。如果拿不准用哪一种变异图，可选用线性变异图，大多数情况下，效果较好。
③最小曲率法 这是一种在地学中广泛应用的网格化方法。由最小曲率法构成的插值表面像一个线性弹性薄板，试图在尽可能严格地尊重数据的同时，生成与原始数据点尽可能吻合的最平滑的曲面。最小曲率法不是准确插值，是典型的平滑插值。使用最小曲率法时要涉及两个参数：最大残差参数和最大循环次数参数，来控制最小曲率的收敛标准。
最大残差（max reciduals）：单位与数据相同，比较合适的值是数据精度的10%。缺省的最大残差为0.001（Zmax-Zmin）
最大循环次数参数（max iterations）：通常设为网格结点数的1～2倍。例如，对于50×50的网格，最大重复参数在2500～5000。
内部和边缘张性系数（internal and boundary tension）：设定弹性薄板内部和边缘弯曲度的参数。该值愈大，弯曲愈小。缺省值均为0。
松弛系数（relaxation factor）：算法参数，通常该值愈大，迭代算法会聚愈快。缺省值为1，一般不用另设定。
④谢别德法 谢别德法是一种加权反距离法的最小二乘法。与加权反距离插值法相似，但由于使用局部最小二乘法，消除或减少了绘制等值线时的“牛眼”效应。谢别德法可以是准确插值或者是平滑插值。可以设置网格化的平滑参数，允许进行平滑插值。随平滑参数值的增加，平滑效果愈明显。通常，该值在0～1最合适。
⑤最近邻点法 又称泰森多边形方法，泰森多边形（Thiesen，又叫Dirichlet或Voronoi多边形）分析法是荷兰气象学家A.H.Thiessen提出的一种分析方法。最初用于从离散分布气象站的降雨量数据中计算平均降雨量，现在GIS和地理分析中经常采用泰森多边形进行快速赋值。实际上，最近邻点插值的一个隐含的假设条件是任一网格点P（x，y）的属性值都使用距它最近的位置点的属性值，用每一个网格节点的最邻点值作为待定节点值。当数据已经是均匀间隔分布，要先将数据转换为Surfer的网格文件，可以应用最近邻点插值法；或者在一个文件中，数据紧密完整，只有少数点没有取值，可用最近邻点插值法来填充无值的数据点。有时需要排除网格文件中的无值区域，再搜索椭圆（search ellipse）设置一个值，对无值区域赋予该网格文件里的空白值。设置的搜索半径的大小要小于该网格文件数据值之间的距离，所有的无数据网格节点都被赋予空白值。在使用最近邻点插值网格化法将一个规则间隔的XYZ数据转换为一个网格文件时，可设置网格间隔和XYZ数据的数据点之间的间距相等。最近邻点插值网格化法没有选项，它是均质且无变化的，对均匀间隔的数据进行插值很有用，同时，它对填充无值数据的区域很有效。
最近邻点法用最邻近的数据点来计算每个网格结点的值。这种方法通常用于已有规则网格只需要转换为Surfer网格文件时，或数据点几乎构成网格，只有个别点缺失，该方法可以有效地填充“空洞”。通过设置搜寻椭圆半径的值小于数据点之间距离的方法，给缺少数据点的结点赋值为空白。最近邻点插值网格化法没有选项，它是均质且无变化的，对均匀间隔的数据进行插值很有用，同时，它对填充无值数据的区域很有效。
⑥多项式回归法 多元回归被用来确定数据的大规模的趋势和图案。可以用几个选项来确定需要的趋势面类型。多元回归实际上不是插值器，因为它并不试图预测未知的Z值。它实际上是一个趋势面分析作图程序。使用多元回归法时要涉及曲面定义和指定XY的最高方次设置，曲面定义是选择采用的数据的多项式类型，这些类型分别是简单平面、双线性鞍、二次曲面、三次曲面和用户定义的多项式。参数设置是指定多项式方程中X组元和Y组元的最高方次。多项式回归法用来确定用户数据整体的趋势或构造一种模型。多项式回归法实际上并不是一种插值，因为它并不试图预测未知的Z值。用户可以在表面选定（surface difination）框内选择4种曲面中的任一种。
⑦自然邻点插值法 自然邻点插值法（natural neighbor）是Surfer 7.0开始才有的网格化新方法，广泛应用于一些研究领域中。其基本原理是对于一组泰森（Thiessen）多边形，当在数据集中加入一个新的数据点（目标）时，就会修改这些泰森多边形，而使用邻点的权重平均值将决定待插点的权重，待插点的权重和目标泰森多边形成比例。实际上，在这些多边形中，有一些多边形的尺寸将缩小，并且没有一个多边形的大小会增加。同时，自然邻点插值法在数据点凸起的位置并不外推等值线（如泰森多边形的轮廓线）。
⑧径向基本函数法 径向基本函数法是多个数据插值方法的组合。根据适应的数据和生成一个圆滑曲面的能力，其中的复二次函数被许多人认为是最好的方法。所有径向基本函数法都是准确的插值器，它们都要为数据而努力。为了试图生成一个更圆滑的曲面，对所有这些方法都可以引入一个圆滑系数。可以指定的函数类似于克里格中的变化图。当对一个格网结点插值时，这些个函数给数据点规定了一套最佳权重。
⑨带线性插值的三角剖分法 是一种严密的准确插值，它的工作路线与手工绘制等值线相近。方法是在相邻点之间连线构成三角形，并且保持任一三角形的边都不与其他三角形的边相交。这样在网格范围内由一系列三角形平面构成拼接图。由于数据点平均分布，在通过地形变化显示断层线时，三角形法非常有效。因为每一个三角形都构成一个平面，所有的结点都在三角形中，其坐标被三角形平面方程唯一地确定。对于有200～1000个数据点，且平均地分配在网格区域里时，用带线性插值的三角剖分法最好。
（3）常用的网格化方法对比分析 不同的网格化方法可以得到不同的网格文件，用户应当选用最能代表自己数据特点的方法，选择网格化方法时应当考虑原始数据点数量的多寡。
①从原始数据点数量角度考虑 10个或10个以下的数据点，除了反映数据的一般趋势外，没有多大意义。这样少的点，带线性插值的三角剖分法无效，数据点<250个时，具线性变异图的克里格法、多重二次曲面法的径向基函数法都可以产生较好代表原始数据特点的网格。
中等数据量（250～1000数据点），带线性插值的三角剖分法网格化很快，并能生成很好代表原始数据特点的网格。克里格法和径向基函数法较慢，也可以产生高质量的网格。
大的数据量（>1000数据点），最小曲率法最快，网格足以代表原始数据特点。带线性插值的三角剖分法网格化较慢，网格有足够的代表性。
②常用网格化法特征 大部分情况下，具有线性变异图的克里格法是十分有效的，应首先予以推荐。其次是很接近的径向基函数法中的多重二次曲面法。这两种方法都能产生较好的代表原始数据的网格。但对于大量数据的网格化，克里格法比较慢。加权反距离法最快，但是围绕数据点，有产生“牛眼”效应的趋势。
谢别德法与加权的反距离法插值法相似，但没有产生等值线“牛眼”效应的缺点。带线性插值的三角剖分法对于中等数量的数据点，网格化很快。一个优点是，当有足够的数据点时，三角剖分法可以反映出数据文件所内含的不连续性。例如断层线。
③常用的8种网格化方法及其特征和应用条件 如表2.15所示。
表2.15 常用的8种网格化方法及其特征和应用条件
（4）网格化方法选择 选择网格化方法时应当考虑原始数据点数量的多寡，≤10个数据或<250个时，具线性变异图的克里格法、多重二次曲面法的径向基函数法都可以产生较好代表原始数据特点的网格。中等数据量（250～1000个数据点）时，带线性插值的三角剖分法网格化很快，并生成很好代表原始数据特点的网格。克里格法和径向基函数法也可以产生高质量的网格。大的数据量（>1000个数据点），最小曲率法最快，网格足以代表原始数据特点。带线性插值的三角剖分法网格化较慢，网格有足够的代表性。
推荐方法：大部分情况下，具有线性变异图的克里格法是十分有效的，应首先予以推荐。其次是径向基函数法。径向基函数法十分灵活，与克里格法产生的网格十分类似。改进的谢别德（Shepard’s）法与反距离加权插值法相似，但没有产生等值线“牛眼”效应的缺点。
（5）采样点值必须保持 必须保证采样点的值不被插值算法改变：反距离加权插值法、克里格插值法、径向基函数法，在插值点与取样点重合时，插值点的值就是样本点的值，其他方法不能保证如此。
（6）交叉验证 各种算法都可以检验插值质量，这就是“交叉验证”。即移去一个已知资料点的数据，用其他各点的数据来估计该点的数据，将插值数据和真实数据进行比较，以检验插值精度的方法。
一般这样的验证都是全交叉验证，即所有的资料点都要进行验证。对于验证的结果，运用绝对平均误差（MAE）、相对平均误差及均方根误差（RMSE）作为检验的质量标准。而检验结果的报告可以作为文件保存。
2.3.2 污染物扩散迁移推荐模型
进入土壤中的污染物可在土壤液相、气相和固相分配并达到平衡。表层、下层土壤及地下水中的挥发性污染物可扩散进入室外空气，下层土壤和地下水挥发性污染物可扩散进入室内空气，土壤中污染物可淋溶、迁移进入地下水。以下给出了土壤和地下水中污染物扩散迁移的相关模型。
2.3.2.1 气态污染物有效扩散系数计算模型
（1）土壤中气态污染物的有效扩散系数
式中，
为土壤中气态污染物的有效扩散系数，cm2/s；Da为空气中扩散系数，cm2/s；Dw为水中扩散系数，cm2/s；H'为无量纲亨利常数，cm3/cm3；θ为非饱和土层土壤中总孔隙体积比，无量纲；θws为非饱和土层土壤中孔隙水体积比，无量纲；θas为非饱和土层土壤中总孔隙空气体积比，无量纲。
其中
式中，ρb为土壤容重，kg/dm3，推荐值可查阅相关手册；ρs为土壤颗粒密度，kg/dm3，推荐值可查阅相关手册；Pws为土壤含水率，kg/dm3，土壤推荐值可查阅相关手册；ρw为水的密度，1kg/dm3。
（2）气态污染物在地基与墙体裂隙中的有效扩散系数
式中，
为气态污染物在地基与墙体裂隙中的有效扩散系数，cm2/s；θacrack为地基裂隙中空气体积比，无量纲；θwcrack为地基裂隙中水体积比，无量纲。
（3）毛细管层中气态污染物的有效扩散系数
式中，
为毛细管层中气态污染物的有效扩散系数，cm2/s；θacap为毛细管层土壤中孔隙空气体积比，无量纲；θwcap为毛细管层土壤中孔隙水体积比，无量纲。
（4）气态污染物从地下水到表层土壤的有效扩散系数
式中，
为地下水到表层土壤的有效扩散系数，cm2/s；hcap为地下水土壤交界处毛细管层厚度，cm，推荐值可查阅相关手册；hv为非饱和土层厚度，cm，优先根据场地调查数据确定，推荐值可查阅相关手册；Lgw为地下水埋深，cm，必须根据场地调查获得参数值。
（5）土壤-水中污染物分配系数
式中，Ksw为土壤-水中污染物分配系数，cm3/g；Kd为土壤固相-水中污染物分配系数，cm3/g。
其中
式中，Koc为土壤有机碳/土壤孔隙水分配系数，L/kg；foc为土壤有机碳质量分数，无量纲；fom为土壤有机质含量，g/kg，根据场地调查获得参数值。
（6）室外空气中气态污染物扩散因子
式中，DFoa为室外空气中气态污染物扩散因子，［g/（cm2·s）］/（g/cm3）；Uair为混合区风速，cm3/s；A为污染源区面积，cm2；W为污染源区宽度，cm2；δair为混合区高度，cm。
（7）室内空气中气态污染物扩散因子
式中，DFia为室内空气中气态污染物扩散因子，［g/（cm2·s）］/（g/cm3）；ER为室内空气交换速率，次/d；LB为室内空间体积与气态污染物入渗面积比，cm；86400为时间单位转换系数，s/d。
（8）流经地下室底板裂隙的对流空气流速
式中，Qs为流经地下室地板裂隙的对流空气流速，cm3/s；π为圆周率常数，3.14159；dP为室内和室外大气压力差，g/（cm·s）；Kv为土壤透性系数，cm2；Xcrack为地下室内地板（裂隙）周长，cm；μair为空气黏滞系数，1.81×10-4g/（cm·s）；Zcrack为地下室地面到地板底部厚度，cm；Rcrack为室内裂隙宽度，cm；Ab为地下室内地板面积，cm2；η为地基和墙体裂隙表面积占室内地表面积的比例，无量纲。
2.3.2.2 污染物扩散进入室外空气的挥发因子计算模型
（1）表层土壤中污染物扩散进入室外空气的挥发因子
VFsuroa1为表层土壤总污染物扩散进入室外空气的挥发因子（算法一），kg/m3；VFsuroa2为表层土壤总污染物扩散进入室外空气的挥发因子（算法二），kg/m3；VFsuroa为表层土壤总污染物扩散进入室外空气的挥发因子（算法一和算法二中的较小值），kg/m3；τ为气态污染物入侵持续时间，a；d为表层污染土壤层厚度，cm，必须根据场地调查获得参数值；31536000为时间单位转换系数，s/a。
（2）下层土壤中污染物扩散进入室外空气的挥发因子
如下层污染土壤厚度已知，污染物进入室外空气的挥发因子：
式中，VFsuboa1为下层土壤中污染物扩散进入室外空气的挥发因子（算法一），kg/m3；VFsuboa2为下层土壤中污染物扩散进入室外空气的挥发因子（算法二），kg/m3；VFsuboa为下层土壤中污染物扩散进入室外空气的挥发因子（算法一和算法二中的较小值），kg/m3；Ls为下层污染土壤上表面到地表距离，cm，必须根据场地调查获得参数值；ds为下层污染土壤厚度，cm。
（3）地下水中污染物扩散进入室外空气的挥发因子
式中，VFgwoa为地下水中污染物扩散进入室外空气的挥发因子，L/m3。
2.3.2.3 污染物扩散进入室内空气的挥发因子计算模型
（1）建筑物下方土壤中污染物进入室内空气的挥发因子
Qs=0时
Qs>0时
如下层污染土壤厚度已知，污染物进入室内空气的挥发因子：
式中，VFsubia1为下层土壤中污染物扩散进入室内空气的挥发因子（算法一），kg/m3；VFsubia2为下层土壤中污染物扩散进入室内空气的挥发因子（算法二），kg/m3；VFsubia为下层土壤中污染物扩散进入室内空气的挥发因子（算法一和算法二中的较小值），kg/m3；Lcrack为室内地基或墙体厚度，cm，推荐值可查阅相关手册；ξ为土壤污染物进入室内挥发因子计算过程参数；31536000为时间单位转换系数，s/a。
（2）地下水中污染物扩散进入室内空气的挥发因子
Qs=0时
Qs>0时
式中，VFgwia1为地下水中污染物扩散进入室内空气的挥发因子（算法一），kg/m3；VFgwia2为地下水中污染物扩散进入室内空气的挥发因子（算法二），kg/m3；VFgwia为地下水中污染物扩散进入室内空气的挥发因子（算法一和算法二中较小值），kg/m3。
2.3.2.4 污染物迁移进入地下水的淋溶因子计算模型
土壤中污染物迁移进入地下水的淋溶因子
如下层污染土壤厚度已知，污染物迁移进入地下水的淋溶因子：
式中，LFsgw1为土壤中污染物迁移进入地下水的淋溶因子（算法一），kg/m3；LFspw-gw为土壤孔隙水中污染物迁移进入地下水的淋溶因子（土壤孔隙水与地下水中污染物浓度的比值），无量纲；LFsgw2为土壤中污染物迁移进入地下水的淋溶因子（算法二），kg/m3；LFsgw为土壤中污染物迁移进入地下水的淋溶因子（算法一和算法二中的较小值），kg/m3；Ugw为地下水的达西速率，cm/a，推荐值可查阅相关手册；δgw为地下水混合区厚度，cm，推荐值可查阅相关手册；I为土壤中水的渗透速率，cm/a，推荐值可查阅相关手册。