胡愷崢 袁海明 陳達寶 伊曉穎 周燕微

(紹興文理學院 數學系，浙江 紹興3312000)

隨著城市經濟的快速發(fā)展和城市人口的不斷增加，人類活動對城市環(huán)境質量的影響日顯突出.如何找到城市大氣和土壤環(huán)境污染的污染源并予以治理，成為人們日益關注的焦點.

對于城市污染問題的研究，一般是通過對城市中污染元素含量的測定，并根據一定的標準來判斷此城市的污染狀況.本文的研究方法是：在一定范圍內的具體位置進行數據測量，測定其污染元素及其含量，通過這些觀測點的污染元素含量來研究整個城市的污染狀況.本文以大氣和土壤污染為例，利用質量守恒定律來建立模型，解出污染物的濃度函數，從而確定污染源，以加快污染問題的治理進程.

1 大氣污染擴散模型

大氣中的污染氣體主要來源于工業(yè)廢氣、汽車尾氣的排放，而這些污染氣體主要是通過擴散進行傳播.所以，我們將根據流體動力學中的質量守恒定律來研究大氣擴散模型[1-2].

1.1 大氣污染擴散模型的建立

我們以污染源為原點建立坐標系，定義距離污染源l處的污染元素濃度為ρ(x,y,z,l)，因此，污染源處污染元素濃度為ρ(x,y,z,0).

圖1 空間域示意圖

空間域Ω內污染元素的增量為

從污染源釋放污染元素的總量為

其中p0為各個污染物濃度的平均值.

根據質量守恒定律和氣體泄漏的連續(xù)性原理[3]，即單位距離內通過所選球面S1向外擴散污染元素的質量與球面S1內污染元素增量之和，等于污染源在單位距離內向外釋放污染元素的質量，即Q0=Q+ΔQ，所以

又根據曲面積分的Gauss公式

所以

即

由此可得

這樣我們就得到了氣體擴散的模型為

1.2 模型的求解及污染源的確定方法

求解上述氣體擴散模型可得[5]，任意污染處P(x,y,z,l)的污染元素濃度為

(1)

其中k為常數.

通過樣本數據處理，污染物濃度的平均值為p0=176.4125ug/g.由于k和D只影響模型的峰度，不影響模型中心位置的坐標.經實驗說明k的取值范圍在[0.05,0.25]之間比較符合客觀事實，其中擴散系數D的取值范圍為[0,1].因此，在(1)式中，不妨令系數k=0.1，擴散系數D=0.1，距離l=1 000 000m.

圖2 氣體擴散機制圖

下面求污染源的坐標.先對大量測量數據進行篩選，將濃度偏差較大的異常點刪除，再把處理后的污染元素濃度值ρ(x,y,z,l)代入公式(1)，利用程序可得到該濃度下的一組解.如圖2所示，根據氣體擴散機制[5]，即污染源以球面的形式向外擴散，距污染源越遠污染元素的濃度越低，且離污染源相同距離處污染物元素的濃度相同，便可尋找到該組坐標下的圓心，由此可得該濃度下的污染源坐標.把測得的濃度代入(1)式，將計算得到的該濃度下的坐標和現實中污染元素的坐標進行比較，再對求得的污染源坐標進行坐標變換，可得現實坐標系下的污染源坐標.在上述的各系數的取值下，通過程序計算得到，其污染源位置為x=5515.6,y=10722.

我們將理論結果與實際情況進行對比，發(fā)現此氣體傳播模型能較好地確定污染源的位置，所以本模型能較好模擬出氣體擴散的規(guī)律，并找出污染源.這也為各種氣體泄漏、瓦斯爆炸等突發(fā)事件有效快速的解決提供了可能，從而節(jié)約了一定的社會成本.此外，也為治理城市污染提供了有效途徑.

2 土壤重金屬傳播模型

2.1 土壤重金屬傳播模型的建立

地表徑流沖刷是表層土壤中重金屬的主要傳播途徑[6]，因而對流是最主要的傳播特征，而擴散所起的作用相對來說較小.

記ρ=ρ(x,y,t)為t時刻位于點M(x,y)處重金屬元素的濃度.利用微元法，任取一小段時間[t,t+Δt]及任意區(qū)域ω，為方便起見，記x=(x,y)T，則在微小時間段[t,t+Δt]內，污染物通過區(qū)域ω的邊界?ω進入該區(qū)域,由流體的質量守恒定律可得

其中v為污染元素的流速，即污染物在單位時間內通過的距離.

由Green公式，上式可寫為

由t和ω的任意性，可得對流方程[7]

(2)

假設污染物在土壤中的流速較小，則由達西定律[8]可得

(3)

將(3)代入(2)，得到

即

給定初始時刻t=t0時污染元素的濃度

ρ(x,t0)=ρ0(x)，

由此得到重金屬的傳播模型為

(4)

2.2 模型的求解及污染源的確定方法

由于在模型的求解過程中需要用到高度函數h，因此需要先根據測得數據的高度信息確定高度函數h.又考慮到我們測量到的數據分布呈散亂狀況，因而采用散亂數據插值方法得到高度函數.利用Matlab軟件，運用“V4”插值法[9]得到高度地形分布圖，如圖3.

圖3 “V4”插值地形分布圖

由于模型(4)是一個一階雙曲型方程，可以利用特征法[10]求解.(4)的特征線滿足

(5)

記過(t0,x0)的特征線為x(t,x0)，重金屬元素濃度ρ滿足

(6)

從而沿特征線的解為

(7)

假設T時刻測量點處污染元素的濃度為ui(i=1,2,…,n)，則初始條件可寫為

ρ(xi,yi,T)=ρi,(i=1,2,...,n).

由于只有一個時間點上的截面數據，無法確定系數k的值，因此對時間t作伸縮變換，即令

t′=kt，

則方程(6)可改寫為(為方便起見，仍記t'為t)

由特征線法可得，沿特征線(7)污染元素的濃度為

(8)

其中ρT為T時刻xT處測得的污染元素的濃度.

對(5)和(8)分別作離散化，得到迭代格式

xm+1=xm+▽h(xm)Δt，

(9)

(10)

其中ρ0為T時刻污染元素的初始濃度.

若在某一給定點x0處的濃度為ρ0，則由(9)及(10)可以求出各時刻t(t≤T)沿著經過點xo處的特征線上的濃度值ρ(x,t;x0).這樣，當給定污染濃度的閾值U時，對給定時刻t1(t1≤T)，濃度超過閾值U的區(qū)域為

Bi={x(t1)|ρ(x(t1),t1;x0)≥U}，

此區(qū)域即為可能的污染源所在的區(qū)域.

根據特征線法我們畫出了污染元素濃度圖，如圖4；根據測得的污染物數據，我們畫出了污染元素濃度的測量值，如圖5.從圖中看，我們發(fā)現用特征線法求得的污染元素濃度與測量得到的濃度非常接近.同時，我們將坐標M(x,y)的值代入(7)式計算得到的濃度與測量得的濃度進行比較，發(fā)現兩值也非常接近.

圖4 用特征線法得到的污染元素的濃度 圖5 污染元素濃度的測量值

本模型適用于各種由于對流造成污染的重金屬，包括河流中的污染等，不僅可以根據測量的污染物濃度確定污染源，為治理污染省去大量的人力物力，加快治理污染的過程，而且可以根據污染源計算出距離污染源多少遠處污染物的濃度，為合理建造化工廠等提供相應的對策.

3 結論

我們利用流體力學中的質量守恒定律來建立城市大氣污染擴散模型和土壤污染擴散規(guī)律，得到污染物的濃度函數，由此得出污染源.這為城市氣體泄漏、瓦斯爆炸、土壤重金屬污染等突發(fā)事件提供了有效快速的解決途徑，有利于快速找到泄漏地點，解決突發(fā)事件，且省去了大量的人力、物力，節(jié)約了大量的社會成本.

致謝

本文是在陳志祥和胡金杰老師的指導下完成的，在此表示謝意.

參考文獻:

[1]姜啟源，謝金星，葉俊.數學模型[M].北京：高等出版社，2003.

[2]王冬琳.數學建模及實驗[M].北京：國防工業(yè)出版社，2004.

[3]劉振翼，周軼，黃平等.CO2管線泄漏擴散小尺度實驗研究[J].化工學報，2012年，63卷(5)：1655-1657.

[4]何龍慶，林繼成，石冰.菲克定律與擴散的熱力學理論[J].安慶師范學院學報，2006，12(4):38-39.

[5]石東偉，陳冬娜.高斯擴散模型在確定污染源位置中的應用[J].河南科技學院學報，2012年，40卷(2)：57-58.

[6]高存昌.沼氣在農業(yè)面源污染治理中的作用[J].河南農業(yè)，2012(5)：6-7.

[7]Ockendon J，Howison S，Lacey A，et al.應用偏微分方程[M].北京：科學出版社，2008.

[8]蔣玄葦，張野，王嘉舜.地下水在平面滲流狀態(tài)下的流速測定[J].科技資訊，2012(8)：100-101.

[9]馬莉.Matlab數學實驗與建模[M].北京：清華大學出版社，2010.

[10]Lin Changrong1,Wang Shangxu,Zhang Yong.Predicting the distribution of reservoirs by applying themethod of seismic data structure characteristics: Examplefrom the eighth zone in Tahe Oilfield[J].Applied Geophysics,2006(4)：234-240.