基于高斯扩散模型的重金属污染源求解

高 振，房 越

（南京邮电大学，江苏 南京 210003）

1 问题背景重述

近年来，人们对查证城市土壤地质环境异常、评价城市环境质量，研究人类活动影响下城市地质环境的演变模式等问题变得日益关注。根据文献[1]所给的附件采样点的位置海拔高度、主要重金属元素浓度分析重金属污染物的传播特征。建立模型确定污染源的位置。

2 模型的建立与求解

2.1 基于高斯扩散模型的偏微分方程的建立

Step1：正态分布的确定。采用统计学方法研究污染物质在土壤中的传播扩散模型[2]。如果从原点释放出一堆粒子，经过一段时间T之后，浓度区域一个稳定的状态。粒子呈现正态分布，如下图所示其满足的形式满足：

图1 正态分布图

其中f( y )代表粒子的概率密度，σ为此方向上的传播标准差，μ为概率密度最大的点，即浓度最大的点，即传播源所在的位置。

Step2：连续点源的扩散方程。假设污染物在传播的途中满足一下几个条件：

①土壤流场稳定，传播速度均匀；②金属污染物在x,y,z方向上满足正态分布；③金属污染物在传播的途中满足质量守恒；④污染源强连续均匀。得出如下的方程：

其中f( x, y, z )代表粒子的概率密度，σx, σy,σz分别为x, y, z上的传播标准差，(x0,y0,z0)为概率密度最大的点，即浓度最大的点，即传播源所在的位置。

Step3：污染浓度求解。已知从一定量的金属污染物粒子从传播点扩散的概率密度方程。概率密度与在污染源点存在的污染物浓度相乘，得到在每一个点的金属污染物的浓度，则测试点(x, y, z)的污染物浓度如下（Q为污染源处的浓度）：

c(x ,y,z )=Q ×f(x,y,z )

我们可以得到在测试点(x, y, z)的最终的污染物浓度表达式：

Step4：简化方程。

然后对两边同时取对数，得到化简后的方程为：

2.2 污染源范围的确定

根据附件中的数据得到的污染物三维立体图上颜色较深的区域，找出区域的的坐标及其对应的编号，即为污染较严重的区域。砷的污染源在样本点84，18附近。

2.3 非线性最小二乘法

结合我们已经求出的金属污染物传播方程，其为非线性多元函数。要求出其中的未知参数，我们选择用非线性最小二乘法。求解非线性最小二乘，其导数形式太过于复杂无法直接写出，我们通过不停的迭代计算寻找到函数的局部最小值，并认为该局部最小值能够使得我们的要求解的函数取得最优解。

Step1：局部最小值。非线性最小二乘问题形式，下面称为代价函数：

同时定义局部最小值为，对于代价函数F，有变量满足：

中，当x1, y1, z1分别代表x, y, z。求出所有满足该函数一阶导等于零的点，得驻点：

综上所述，局部最小值为函数的驻点，但是驻点也有可能是函数的局部最大值或者鞍点。为了区分到底是局部最小值、局部最大值或者是鞍点，我们保留泰勒展开式的二阶信息。通过判定展开式中的矩阵的相关性质进行判定：

当该矩阵为正定矩阵时，所求驻点是局部最小值。

Step2：运用算法求解最小二乘问题。

因高斯牛顿法中的泰勒展开只有在展开点附近才会有比较好的效果，因此为了确保近似的准确性我们需要设定一个具有一定半径的区域作为信赖区域。我们需要明确信赖区。在L-M算法中信赖区大小的确定也是运用增益比例来进行判定的。基于信赖区域能够重新构建更有效的优化框架。引入拉格朗日乘子将有约束优化问题转化为无约束优化问题。求出L-M算法增量表达式。通过MATLAB中“nlin fit”函数可以直接调用这种方法，求出的As金属的污染源坐标[18122，9315]、[12930，3400]。