高 振，房 越

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

1 問題背景重述

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

2 模型的建立與求解

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

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

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

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

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

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

其中f( x, y, z )代表粒子的概率密度，σx, σy,σz分別為x, y, z上的傳播標(biāo)準(zhǔn)差，(x0,y0,z0)為概率密度最大的點(diǎn)，即濃度最大的點(diǎn)，即傳播源所在的位置。

Step3：污染濃度求解。已知從一定量的金屬污染物粒子從傳播點(diǎn)擴(kuò)散的概率密度方程。概率密度與在污染源點(diǎn)存在的污染物濃度相乘，得到在每一個(gè)點(diǎn)的金屬污染物的濃度，則測(cè)試點(diǎn)(x, y, z)的污染物濃度如下（Q為污染源處的濃度）：

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

我們可以得到在測(cè)試點(diǎn)(x, y, z)的最終的污染物濃度表達(dá)式：

Step4：簡(jiǎn)化方程。

然后對(duì)兩邊同時(shí)取對(duì)數(shù)，得到化簡(jiǎn)后的方程為：

2.2 污染源范圍的確定

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

2.3 非線性最小二乘法

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

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

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

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

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

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

Step2：運(yùn)用算法求解最小二乘問題。

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