基于高斯模型對(duì)京津冀地區(qū)空氣污染的研究

張德鑫，朱家明，朱晗，余湉

(安徽財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030)

?

基于高斯模型對(duì)京津冀地區(qū)空氣污染的研究

張德鑫，朱家明，朱晗，余湉

(安徽財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030)

高斯模型；空氣污染；模糊聚類；綜合評(píng)價(jià)；F檢驗(yàn)

針對(duì)京津冀地區(qū)的空氣污染，使用模糊C均值聚類、模糊綜合評(píng)價(jià)、F檢驗(yàn)，分類討論等方法，分別構(gòu)建空氣質(zhì)量優(yōu)劣衡量體系、連續(xù)高架點(diǎn)源高斯擴(kuò)散、瞬時(shí)高架點(diǎn)源高斯擴(kuò)散、連續(xù)線源高斯擴(kuò)散等模型，使用MATLAB、SPSS軟件，得到影響空氣質(zhì)量的關(guān)鍵參數(shù)。研究結(jié)果表明，京津冀地區(qū)污染的主要原因分別是機(jī)動(dòng)車的大量使用、工業(yè)廢氣排放以及煤炭與燃油等化石能源燃燒等。

0 引言

人類生產(chǎn)生活排放大量大氣污染物，導(dǎo)致城市空氣污染日益加重，而京津冀地區(qū)的空氣污染問題尤為嚴(yán)重，據(jù)我國環(huán)保部發(fā)布的報(bào)告2013年6月京津冀地區(qū)空氣質(zhì)量重度污染以上天數(shù)占21.2%。空氣質(zhì)量級(jí)別的明確有利于人們的出行安排，所以需要通過建立一個(gè)空氣質(zhì)量評(píng)價(jià)數(shù)學(xué)模型來評(píng)價(jià)各地區(qū)、各時(shí)間段的空氣質(zhì)量狀況，分析其主要污染源的性質(zhì)和種類，以便做出正確的出行選擇，并可在一定程度上改善和治理空氣污染[1]。

1 空氣質(zhì)量優(yōu)劣程度等級(jí)的判別

1.1研究思路

以京津冀地區(qū)為研究對(duì)象，參考國家標(biāo)準(zhǔn)和美國標(biāo)準(zhǔn)建立衡量空氣質(zhì)量優(yōu)劣程度等級(jí)的數(shù)學(xué)模型。以河北石家莊作為代表，從氣象部門網(wǎng)站調(diào)取2014年7月、10月和2015年1月、4月共4個(gè)月份的6號(hào)、16號(hào)、26號(hào)北京、天津、石家莊地區(qū)的空氣質(zhì)量情況相關(guān)指數(shù)共36組數(shù)據(jù)作為樣本并依次排序，對(duì)這些樣本進(jìn)行聚類分析得到衡量空氣質(zhì)量優(yōu)劣程度等級(jí)的數(shù)學(xué)模型。

1.2模糊C均值聚類

首先運(yùn)用MATLAB軟件對(duì)樣本中的數(shù)據(jù)進(jìn)行模糊C均值聚類，結(jié)果如表1中所示。

表1　MATLAB模糊C均值聚類結(jié)果

其中，t1～t6表示將樣本分為6類，其后的數(shù)值代表樣本的序號(hào)。

然后運(yùn)用SPSS軟件對(duì)分類的結(jié)果進(jìn)行Friedman檢驗(yàn)，對(duì)于檢驗(yàn)水平為α=0.04，分別采取各類別之間檢驗(yàn)與其中2類合并以后檢驗(yàn)，顯然各類別之間的檢驗(yàn)接受原假設(shè)比其中2類合并后檢驗(yàn)接受原假設(shè)更多，這說明各類之間差距較大，同一類中個(gè)體之間差異較小，結(jié)果符合聚類的基本原則。

1.2.1 構(gòu)建空氣質(zhì)量等級(jí)標(biāo)準(zhǔn)

通過對(duì)聚類結(jié)果的檢驗(yàn)和分析，將分成的6類結(jié)果人為定義為優(yōu)、良、輕度污染、中度污染、重度污染和嚴(yán)重污染6個(gè)空氣質(zhì)量優(yōu)劣程度等級(jí)，再對(duì)其每個(gè)等級(jí)的各參數(shù)指標(biāo)數(shù)據(jù)進(jìn)行如下處理：

在a、b、c中選取合適的數(shù)值代表每個(gè)空氣質(zhì)量等級(jí)的各指標(biāo)所在范圍，注意同一指標(biāo)的上下范圍要連貫，空氣質(zhì)量優(yōu)劣程度等級(jí)劃分如表2中所示。

表2　空氣質(zhì)量等級(jí)標(biāo)準(zhǔn)

再運(yùn)用公式：

其中，I為空氣質(zhì)量指數(shù)，即FQI，是輸出值；C為污染物濃度，輸入值；Clow表示小于或等于C的濃度限值；Chigh為大于或等于C的濃度限值；Ilow是對(duì)應(yīng)于Clow的指數(shù)限值；Ihigh是對(duì)應(yīng)于Chigh的指數(shù)限值，Clow、Chigh、Ilow、Ihigh皆為常量。利用該公式，根據(jù)污染物濃度C，可以方便地計(jì)算出空氣質(zhì)量指數(shù)I，再與表1相對(duì)照，判斷其空氣質(zhì)量屬于何種等級(jí)。

1.2.2 模型的檢驗(yàn)

假若給定某城市某日24時(shí)各污染物平均濃度，即可按照已建立的模型算出其FAQI，再根據(jù)空氣質(zhì)量等級(jí)標(biāo)準(zhǔn)表對(duì)照判斷出其空氣質(zhì)量狀況如何，可通過按照國家標(biāo)準(zhǔn)和美國標(biāo)準(zhǔn)計(jì)算其空氣質(zhì)量指數(shù)AQI與表3中的國標(biāo)和美標(biāo)濃度限值相對(duì)照，檢驗(yàn)其結(jié)果的正確性。

表3　國家標(biāo)準(zhǔn)和美國標(biāo)準(zhǔn)污染物濃度限值

2 不同污染源下空氣污染擴(kuò)散研究

2.1相關(guān)指標(biāo)公式

(1)公式I 瞬時(shí)污染氣體排放量

瞬時(shí)污染氣體排放量主要針對(duì)的是排放一段時(shí)間污染氣體后短時(shí)間內(nèi)不再排放的情況，假設(shè)排放的氣體不發(fā)生任何化學(xué)反應(yīng)且沒有其它氣體的混入時(shí)，將之前排放量看作是接下來一段時(shí)間內(nèi)的瞬時(shí)排放量，即：

Q=CVΔt

(2)公式Ⅱ單位時(shí)間污染氣體排放量

點(diǎn)源污染氣體的排放連續(xù)，用單位時(shí)間污染氣體的排放量來研究對(duì)濃度的影響，即：

(3)公式Ⅲ 煙羽抬升高度估計(jì)

煙羽抬升是指煙氣離開排放源后因初始動(dòng)力和熱力作用而繼續(xù)上升并逐漸變平的過程。煙囪排出的煙氣既有一定的出口速度，煙氣的溫度也高于出口處的環(huán)境大氣。這樣煙流在動(dòng)力和熱力的作用下，在一定條件下可抬至相當(dāng)?shù)母叨龋笤诳諝庾枇椭亓Φ淖饔孟卵仫L(fēng)的方向上變平。查閱有關(guān)文獻(xiàn)可知，其標(biāo)準(zhǔn)公式：

Δh=2(1.5VD+0.04QH)u-1

式中，V表示煙氣出口速度(立方米每時(shí)),D表示煙囪口內(nèi)徑，QH為排熱量(kcal /s)，u表示煙囪頂處的平均風(fēng)速(m·s-1)。

(4)公式Ⅳ污染源有效高度

污染源的有效高度包括煙囪的實(shí)際高度和煙羽抬升高度：

H=H0+Δh

2.2建立連續(xù)高架點(diǎn)源高斯擴(kuò)散模型

2.2.1研究思路

假設(shè)污染氣體排放的過程中不被吸收且未發(fā)生化學(xué)反應(yīng)。以有效污染源為坐標(biāo)原點(diǎn)，風(fēng)向所指方向?yàn)閄軸正方向，構(gòu)建空間坐標(biāo)系。污染氣體沿X軸向橫側(cè)和垂直方向自由擴(kuò)散，假設(shè)Y軸、Z軸上濃度分布服從二維正態(tài)分布。先從不考慮地面障礙物著手，對(duì)下風(fēng)向點(diǎn)濃度建立高架點(diǎn)源高斯擴(kuò)散模型，再結(jié)合形象圖利用鏡像對(duì)模型進(jìn)行完善，構(gòu)建下風(fēng)向點(diǎn)濃度函數(shù)。而在點(diǎn)源的實(shí)際擴(kuò)散過程中，污染氣體可能會(huì)受到地面障礙物的影響。在凌晨5點(diǎn)至上午8點(diǎn)和下午4點(diǎn)至晚上9點(diǎn)這2個(gè)時(shí)間段上不再排放污染氣體，可將其看成瞬時(shí)排放，構(gòu)建瞬時(shí)高架點(diǎn)源高斯擴(kuò)散模型。而早上9點(diǎn)至下午3點(diǎn)、晚上10點(diǎn)至凌晨4點(diǎn)2個(gè)時(shí)間段上污染源連續(xù)排放污染氣體，對(duì)其構(gòu)建連續(xù)高架點(diǎn)源擴(kuò)散模型。最終將數(shù)據(jù)代入模型進(jìn)行求解，根據(jù)各個(gè)下風(fēng)向點(diǎn)濃度可知空間污染濃度分布，再對(duì)照空氣質(zhì)量指標(biāo)可判斷出空氣質(zhì)量等級(jí)。

2.2.2模型的建立

(1)不考慮地面障礙物的影響

a.下風(fēng)向點(diǎn)濃度分布和空間坐標(biāo)以及單位時(shí)間排放量有關(guān)，構(gòu)建濃度函數(shù)：

C=C(X,Y,Z,t)

b.以有效污染源為坐標(biāo)原點(diǎn)，風(fēng)向所指方向?yàn)閄軸正方向，建立三維空間坐標(biāo)系。污染氣體沿著X軸橫向、垂直自由擴(kuò)散，單位時(shí)間內(nèi)對(duì)下風(fēng)向點(diǎn)濃度函數(shù)進(jìn)行積分等于單位時(shí)間內(nèi)的污染氣體排放量：

?C(X,Y,Z,t)dYdZ=q(t)

c.假設(shè)污染氣體沿X軸，隨風(fēng)向Y軸、Z軸方向上自由擴(kuò)散相互獨(dú)立且服從二維正態(tài)分布。在不考慮地面障礙物的影響下μy=0,μz=0,可構(gòu)建二維正態(tài)分布密度函數(shù)，并化簡可得風(fēng)向點(diǎn)濃度函數(shù)為：

(1)

(2)考慮地面障礙物的影響

a.假設(shè)有效污染源位于與地面相距H處，以有效污染源為坐標(biāo)原點(diǎn)，風(fēng)向所指方向?yàn)閄軸正方向，建立三維空間坐標(biāo)系。假設(shè)下風(fēng)向點(diǎn)K位于離地面z處，根據(jù)不考慮地面障礙物影響下風(fēng)點(diǎn)濃度密度函數(shù)可知：

b.由于考慮地面障礙物的影響，某部分氣體從始源出發(fā)擴(kuò)散碰到地面發(fā)生全反射，又折回地面上空，剛好到達(dá)下風(fēng)向點(diǎn)K。2次濃度發(fā)生疊加，也就是區(qū)別于不考慮地面障礙物的影響。而疊加上去的濃度就相當(dāng)于將有效污染源關(guān)于地面的鏡像點(diǎn)為污染源，不考慮地面障礙物影響自由擴(kuò)散到下風(fēng)向點(diǎn)K的濃度，K點(diǎn)在Z軸方向上坐標(biāo)為(z+H)，可得疊加上去的濃度為：

由此可知，考慮地面障礙物影響下風(fēng)向點(diǎn)K的濃度函數(shù)為：

(2)

2.3建立瞬時(shí)高架點(diǎn)源高斯擴(kuò)散模型

2.3.1模型的分析與構(gòu)建

與上述模型的構(gòu)建思路相同，唯一的區(qū)別就在于對(duì)某一時(shí)間點(diǎn)上的下風(fēng)向點(diǎn)濃度進(jìn)行積分等于瞬時(shí)污染氣體排放量且煙羽抬升高度為零，因此按照上述模型的思路可構(gòu)建瞬時(shí)高架點(diǎn)源高斯擴(kuò)散模型，求得考慮地面障礙物影響下風(fēng)向點(diǎn)K的濃度函數(shù)為：

(3)

2.3.2模型的參數(shù)估計(jì)與求解

(1)擴(kuò)散參數(shù)σy,σz的估計(jì)

查閱有關(guān)文獻(xiàn)，我國GB3804-91《制定地方大氣污染物排放標(biāo)準(zhǔn)的技術(shù)方法》采用如下經(jīng)驗(yàn)公式確定擴(kuò)散參數(shù)σy,σz：

σy=γ1Xα1, σz=γ2Xα2

式中，γ1,α1,γ2,α2稱為擴(kuò)散系數(shù)。這些系數(shù)由實(shí)驗(yàn)確定，在一個(gè)相當(dāng)長的X距離內(nèi)為常數(shù)，可從GB3804-91文獻(xiàn)的表中查取。

(2)煙羽抬升高度估計(jì)

假定該工廠的煙囪內(nèi)徑為d(m),排熱量為qh(kcal/s)，以煙囪出口處的環(huán)境風(fēng)速us( m·s-1)估計(jì)煙囪頂處的平均風(fēng)速。問題三中要求解的各個(gè)時(shí)點(diǎn)上的煙氣出口速度各不相同,由于早上8時(shí)、中午12時(shí)、晚上9時(shí)的煙氣出口速度分別為V1=0m3,V2=1 200m3/h，V3=0m3/h，代入公式可得早上8時(shí)、中午12時(shí)、晚上9時(shí)的煙羽抬升高度分別為：

Δh1=0，Δh2=2(1 800d+0.04qh)us-1，Δh3=0

(3)各時(shí)點(diǎn)污染濃度函數(shù)

早上8點(diǎn)屬于凌晨5時(shí)至上午8時(shí)，符合瞬時(shí)高架點(diǎn)源高斯擴(kuò)散模型,排放濃度C1=1 160 mg/m3，排放速度為V1=5 700 m3/h，間隔時(shí)間Δt=9 h，在下風(fēng)距離X=51 000 m，煙囪高H0=50 m，將數(shù)據(jù)代入公式(3)即可求解出下風(fēng)向點(diǎn)濃度函數(shù)，如下：

中午12點(diǎn)是屬于上午9時(shí)至下午3時(shí)，符合連續(xù)高架點(diǎn)源高斯擴(kuò)散模型，這時(shí)的排放濃度為C2=406.92 mg/m3,排放速度為V2=1 200 m3/h，在下風(fēng)距離X=51 000 m，污染源的有效高度為H=H0+Δh2，將數(shù)據(jù)代入公式(2)可求解出中午12點(diǎn)下風(fēng)向點(diǎn)濃度函數(shù)，形式見上。晚上9點(diǎn)屬于下午4點(diǎn)至晚上9點(diǎn)這段期間，符合瞬時(shí)高架點(diǎn)源高斯擴(kuò)散模型，求解過程同早上8點(diǎn)。

最后，將求得數(shù)據(jù)轉(zhuǎn)化為AQI指數(shù)，對(duì)照空氣質(zhì)量指標(biāo)判斷空氣質(zhì)量等級(jí)。

2.4連續(xù)高架線源高斯擴(kuò)散模型

2.4.1 研究思路

多污染源是指n個(gè)點(diǎn)污染源在n→∞時(shí)的情況，此時(shí)可將無窮多個(gè)密集的點(diǎn)看成連續(xù)狀態(tài)下的污染發(fā)生源，即從點(diǎn)到線，構(gòu)造線源排放模型。利用此模型對(duì)北京環(huán)城高速的汽車污染排放情況進(jìn)行評(píng)判，作出時(shí)序圖分析空氣污染濃度梯度變化，并利用PAQI評(píng)判空氣質(zhì)量等級(jí)，最后基于謹(jǐn)慎性原則再利用AQI二次評(píng)判。

2.4.2模型的構(gòu)建

假設(shè)氣象條件恒定不變且風(fēng)向水平。污染氣體沿X軸向橫側(cè)和垂直方向自由擴(kuò)散，假設(shè)Y軸、Z軸上濃度分布服從二維正態(tài)分布。假設(shè)相對(duì)于觀測(cè)點(diǎn)，高速公路平直且無限長，則利用積分可得高斯線源排放模型。并利用模型，綜合統(tǒng)計(jì)數(shù)據(jù)，求解分析2015年1月16日北京二環(huán)、四環(huán)(5 km)六環(huán)路(20 km)早上8點(diǎn)、中午12點(diǎn)、晚上9點(diǎn)時(shí)空氣污染濃度梯度變化及空氣質(zhì)量等級(jí)。

利用點(diǎn)源的空氣擴(kuò)散高斯模型，對(duì)Y積分得到無限長連續(xù)線源高斯模型模型數(shù)學(xué)表達(dá)式，如下：

再考慮風(fēng)向，則可解得地面濃度為：

其中ρ為任意一點(diǎn)污染物的濃度，q(kg/h)為排放量，H為有效源高H=HS+ΔH,HS為幾何高度，ΔH為煙羽抬升高度，β為風(fēng)向與線源的夾角。

3 北京市空氣狀況分析

3.1問題的簡化

由于高速公路相對(duì)于觀測(cè)點(diǎn)為無限長，不存在線源末端，故選擇采用無限長連續(xù)線源高斯模型模型求解汽車尾氣排放污染濃度。

首先，基于求解需要，對(duì)北京二環(huán)、四環(huán)、六環(huán)圖形進(jìn)行簡化如圖1所示。

圖1　北京二、四、六環(huán)城高速公路(從下向上)簡化圖

收集北京市2015年1月13日至1月16日的天氣情況整理如表4所示。

表4　北京天氣

3.2問題的的求解

3.2.1不考慮公路之間的相互影響

當(dāng)在環(huán)城高速上測(cè)量的時(shí)候，可將與污染物的距離視為0，構(gòu)造極限如下所示：

對(duì)原初式運(yùn)用洛必達(dá)法則求極限值，由于分子分母均為初等函數(shù)，排放量q與x無關(guān)，所以，其極限存在。令

3.2.2考慮公路之間的相互影響

先假定下風(fēng)向?qū)ι巷L(fēng)向的影響忽略不計(jì)?？紤]各公路之間的影響，需要引入各環(huán)城高速之間的距離。通過查閱文獻(xiàn)資料得知，北京市第二環(huán)城高速與第四、第六環(huán)城高速的平均水平距離為6 km與21 km。依照距離以及公式組可以算出公路之間污染物擴(kuò)散的影響濃度。

二環(huán)對(duì)四環(huán)的影響(6km)，用公式(b)和公式(e)，即

令其中的常數(shù)項(xiàng)

綜合自身污染狀況與其它路徑干擾，得到北二環(huán)、四環(huán)、六環(huán)的污染濃度表達(dá)式如下：

由于汽車是北京市氮氧化物的主要排放者，所以直接用氮氧化物的指標(biāo)來計(jì)量空氣質(zhì)量等級(jí)，也就是這里的排污量q。并用氮氧化物排放的指標(biāo)變化來反映北京市空氣污染濃度梯度。通過查找相關(guān)文獻(xiàn)，了解到每千輛汽車每天排出一氧化碳約3 000 kg,碳?xì)浠衔?00～400 kg,氮氧化合物50～150 kg。調(diào)查得北京市二、四、六環(huán)限行前后車流量情況如表5所示。

表5　北京市二、四、六環(huán)車流量/(千輛/h)

表6　北京市二四六環(huán)q值

表7　未限行北京市二四六環(huán)ρ值

表8　限行后北京市二四六環(huán)ρ值

將此值與國家標(biāo)準(zhǔn)比較即可得到空氣質(zhì)量指數(shù)，繼而通過分級(jí)表查詢得到空氣質(zhì)量等級(jí)。

將參數(shù)視為常參數(shù)，考察北京市二、四、六環(huán)限行前后值表與部分?jǐn)?shù)據(jù)，可得到：從時(shí)序上看，濃度從8點(diǎn)到21點(diǎn)順次遞減；從空間上看，四環(huán)的污染物濃度最高，其次是二環(huán)，最末是六環(huán)。

限行前后比對(duì)，所有路線污染物濃度均顯著下降，該政策效果顯著。

4 結(jié)論

(1)通過對(duì)模型的求解結(jié)果分析得知影響空氣質(zhì)量的關(guān)鍵參數(shù)指標(biāo)有：風(fēng)速、大氣情況、污染源模式、污染源排污狀態(tài)(排放濃度與排放速度)、地面是否有障礙物等。

(2)空氣質(zhì)量的優(yōu)劣與人們的生活息息相關(guān)，治理京津冀地區(qū)的空氣污染迫在眉睫。

[1]http://www.saikr.com/wuyi/2/1427.html.

[2]王寒梅,賈秩然. 大氣污染物擴(kuò)散模型應(yīng)用研究[J].遼寧師專學(xué)報(bào),2006,08(02): 78-79.

[3]鄒本東,徐子優(yōu).因子分析法解析北京市大氣顆粒物PM10的來源[J].中國環(huán)境監(jiān)測(cè)，2007, 23(02)：79-85.

[4]孫慶珍,董曉馬.大氣點(diǎn)源擴(kuò)散模型的GIS可視化研究[J].測(cè)繪科學(xué),2011,36: 24-26.

[5]劉克強(qiáng).煙羽抬升估算和觀測(cè)方法[J].輻射防護(hù)通訊,1986,47-49.

[6]王文,于雷.基于高斯線源模式的主要尾氣擴(kuò)散模型綜述[J].交通環(huán)保, 2004,25(05):7-10.

Research on Air Pollution in Beijing-Tianjin-Hebei Region Based on Gaussian Model

ZHANG De-xin, ZHU Jia-ming, ZHU Han, YU Tian

(College of Statistics and Applied mathematics, Anhui University of Finance and Economics, Bengbu Anhui 233030, China)

Gaussian model; air pollution; fuzzy clustering; comprehensive evaluation; and F-test

This article is about the air pollution, using fuzzy c-means clustering, fuzzy comprehensive evaluation, F-test, the methods of classification discussion, each building air quality measurement system, continuous Gaussian diffusion overhead Gaussian diffusion, instantaneous point source elevated point sources, such as continuous line source Gaussian diffusion model, the use of MATLAB, the SPSS software, get the key parameters affecting the quality of the air, and it is concluded that the main causes of pollution in Beijing-Tianjin-Hebei Region are the use of motor vehicles, industrial emissions and conclusions such as the burning of fossil fuels such as coal and oil.

2095-2716(2016)03-0124-09

X823

A