高 振，房 越

（南京郵電大學(xué)，江蘇 南京 210003）

1 問題背景重述

近年來，人們對查證城市土壤地質(zhì)環(huán)境異常、評價城市環(huán)境質(zhì)量，研究人類活動影響下城市地質(zhì)環(huán)境的演變模式等問題變得日益關(guān)注。根據(jù)文獻[1]所給的附件采樣點的位置海拔高度、主要重金屬元素濃度分析重金屬污染物的傳播特征。建立模型確定污染源的位置。

2 模型的建立與求解

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

Step1：正態(tài)分布的確定。采用統(tǒng)計學(xué)方法研究污染物質(zhì)在土壤中的傳播擴散模型[2]。如果從原點釋放出一堆粒子，經(jīng)過一段時間T之后，濃度區(qū)域一個穩(wěn)定的狀態(tài)。粒子呈現(xiàn)正態(tài)分布，如下圖所示其滿足的形式滿足：

圖1 正態(tài)分布圖

其中f( y )代表粒子的概率密度，σ為此方向上的傳播標準差，μ為概率密度最大的點，即濃度最大的點，即傳播源所在的位置。

Step2：連續(xù)點源的擴散方程。假設(shè)污染物在傳播的途中滿足一下幾個條件：

①土壤流場穩(wěn)定，傳播速度均勻；②金屬污染物在x,y,z方向上滿足正態(tài)分布；③金屬污染物在傳播的途中滿足質(zhì)量守恒；④污染源強連續(xù)均勻。得出如下的方程：

其中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：簡化方程。

然后對兩邊同時取對數(shù)，得到化簡后的方程為：

2.2 污染源范圍的確定

根據(jù)附件中的數(shù)據(jù)得到的污染物三維立體圖上顏色較深的區(qū)域，找出區(qū)域的的坐標及其對應(yīng)的編號，即為污染較嚴重的區(qū)域。砷的污染源在樣本點84，18附近。

2.3 非線性最小二乘法

結(jié)合我們已經(jīng)求出的金屬污染物傳播方程，其為非線性多元函數(shù)。要求出其中的未知參數(shù)，我們選擇用非線性最小二乘法。求解非線性最小二乘，其導(dǎo)數(shù)形式太過于復(fù)雜無法直接寫出，我們通過不停的迭代計算尋找到函數(shù)的局部最小值，并認為該局部最小值能夠使得我們的要求解的函數(shù)取得最優(yōu)解。

Step1：局部最小值。非線性最小二乘問題形式，下面稱為代價函數(shù)：

同時定義局部最小值為，對于代價函數(shù)F，有變量滿足：

中，當(dāng)x1, y1, z1分別代表x, y, z。求出所有滿足該函數(shù)一階導(dǎo)等于零的點，得駐點：

綜上所述，局部最小值為函數(shù)的駐點，但是駐點也有可能是函數(shù)的局部最大值或者鞍點。為了區(qū)分到底是局部最小值、局部最大值或者是鞍點，我們保留泰勒展開式的二階信息。通過判定展開式中的矩陣的相關(guān)性質(zhì)進行判定：

當(dāng)該矩陣為正定矩陣時，所求駐點是局部最小值。

Step2：運用算法求解最小二乘問題。

因高斯牛頓法中的泰勒展開只有在展開點附近才會有比較好的效果，因此為了確保近似的準確性我們需要設(shè)定一個具有一定半徑的區(qū)域作為信賴區(qū)域。我們需要明確信賴區(qū)。在L-M算法中信賴區(qū)大小的確定也是運用增益比例來進行判定的?；谛刨噮^(qū)域能夠重新構(gòu)建更有效的優(yōu)化框架。引入拉格朗日乘子將有約束優(yōu)化問題轉(zhuǎn)化為無約束優(yōu)化問題。求出L-M算法增量表達式。通過MATLAB中“nlin fit”函數(shù)可以直接調(diào)用這種方法，求出的As金屬的污染源坐標[18122，9315]、[12930，3400]。