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

(安徽財經(jīng)大學 統(tǒng)計與應用數(shù)學學院，安徽 蚌埠 233030)

?

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

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

(安徽財經(jīng)大學 統(tǒng)計與應用數(shù)學學院，安徽 蚌埠 233030)

高斯模型；空氣污染；模糊聚類；綜合評價；F檢驗

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

0 引言

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

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

1.1研究思路

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

1.2模糊C均值聚類

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

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

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

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

1.2.1 構(gòu)建空氣質(zhì)量等級標準

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

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

表2　空氣質(zhì)量等級標準

再運用公式：

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

1.2.2 模型的檢驗

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

表3　國家標準和美國標準污染物濃度限值

2 不同污染源下空氣污染擴散研究

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

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

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

Q=CVΔt

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

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

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

煙羽抬升是指煙氣離開排放源后因初始動力和熱力作用而繼續(xù)上升并逐漸變平的過程。煙囪排出的煙氣既有一定的出口速度，煙氣的溫度也高于出口處的環(huán)境大氣。這樣煙流在動力和熱力的作用下，在一定條件下可抬至相當?shù)母叨?，后在空氣阻力和重力的作用下沿風的方向上變平。查閱有關(guān)文獻可知，其標準公式：

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

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

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

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

H=H0+Δh

2.2建立連續(xù)高架點源高斯擴散模型

2.2.1研究思路

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

2.2.2模型的建立

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

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

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

b.以有效污染源為坐標原點，風向所指方向為X軸正方向，建立三維空間坐標系。污染氣體沿著X軸橫向、垂直自由擴散，單位時間內(nèi)對下風向點濃度函數(shù)進行積分等于單位時間內(nèi)的污染氣體排放量：

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

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

(1)

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

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

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

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

(2)

2.3建立瞬時高架點源高斯擴散模型

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

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

(3)

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

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

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

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

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

(2)煙羽抬升高度估計

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

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

(3)各時點污染濃度函數(shù)

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

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

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

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

2.4.1 研究思路

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

2.4.2模型的構(gòu)建

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

利用點源的空氣擴散高斯模型，對Y積分得到無限長連續(xù)線源高斯模型模型數(shù)學表達式，如下：

再考慮風向，則可解得地面濃度為：

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

3 北京市空氣狀況分析

3.1問題的簡化

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

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

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

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

表4　北京天氣

3.2問題的的求解

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

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

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

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

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

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

令其中的常數(shù)項

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

[4]孫慶珍,董曉馬.大氣點源擴散模型的GIS可視化研究[J].測繪科學,2011,36: 24-26.

[5]劉克強.煙羽抬升估算和觀測方法[J].輻射防護通訊,1986,47-49.

[6]王文,于雷.基于高斯線源模式的主要尾氣擴散模型綜述[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