， ，

(1.中原工學(xué)院，鄭州 450007; 2.河南省計(jì)量科學(xué)研究院，鄭州 450008；3.中國(guó)礦業(yè)大學(xué)(北京)，北京 100083)

城市大氣可吸入顆粒物污染已經(jīng)成為大氣科學(xué)研究的熱點(diǎn)問題[1].通過對(duì)可吸入顆粒物濃度變化的分析，可以深入了解可吸入顆粒物污染的來源、變化規(guī)律及影響因素，從而采取有針對(duì)性的措施[2-3].受污染源、氣象條件等多種因素的影響，實(shí)際監(jiān)測(cè)出的顆粒物污染物濃度時(shí)間序列數(shù)據(jù)有很大的無規(guī)律變化，其中有許多是由于偶然因素造成的，這些偶然因素變化嚴(yán)重影響其變化規(guī)律及變化原因的分析[4].消除由于偶然因素變化造成的影響，分析污染物濃度的真實(shí)變化規(guī)律，從信號(hào)處理的角度來看是典型的去噪處理[5].小波分析技術(shù)出現(xiàn)后,已經(jīng)有一些研究用小波分析方法對(duì)大氣污染物的時(shí)間序列進(jìn)行去噪處理，取得了不錯(cuò)的效果[6-7].本文采用MATLAB小波分析對(duì)PM10質(zhì)量濃度時(shí)間序列數(shù)據(jù)進(jìn)行去噪處理,分析了去噪處理的具體過程，并對(duì)北京市2010-2011年度PM10質(zhì)量濃度隨時(shí)間變化的數(shù)據(jù)進(jìn)行了小波去噪分析.

1 時(shí)間序列數(shù)據(jù)小波去噪處理的原理

1.1 時(shí)間序列數(shù)據(jù)的去噪處理

信號(hào)處理時(shí)，一般認(rèn)為低頻信號(hào)是有用信號(hào)而把高頻信號(hào)看作噪聲，要求有用信號(hào)和噪聲的頻譜相互分開[5,8-9].

對(duì)時(shí)間序列數(shù)據(jù)去噪的傳統(tǒng)方法主要有移動(dòng)平均法、傳統(tǒng)濾波方法、卡爾曼濾波方法和維納濾波方法[10].傳統(tǒng)濾波方法的傅里葉(Fourier)變換是將時(shí)域信號(hào)變換到頻域，但由于傅里葉變換中采樣間隔都是常數(shù)，時(shí)間域與頻率域之間彼此是整體刻畫，不能用于局部分析，難以實(shí)現(xiàn)信噪的有效分離[11].

小波分析(或小波變換)具有自適應(yīng)性質(zhì)和數(shù)學(xué)顯微鏡性質(zhì),與傳統(tǒng)的傅里葉變換相比，能夠提供信號(hào)在時(shí)-頻域上的局部化特征，特別適合非平穩(wěn)、非線性信號(hào)的處理[10].將小波變換用于信號(hào)去噪，能在去噪的同時(shí)不損壞信號(hào)的突變部分，是一種對(duì)信號(hào)的局部頻譜進(jìn)行分析的比較理想的數(shù)學(xué)工具[12].

1.2 小波去噪處理的MATLAB實(shí)現(xiàn)

小波閾值消噪法是小波消噪方法中應(yīng)用較廣泛的一種方法,它具有方法簡(jiǎn)單、計(jì)算量小、去噪效果好的特點(diǎn).小波閾值消噪的主要理論依據(jù)為:屬于空間的信號(hào)在小波域內(nèi)的能量主要集中在有限的幾個(gè)系數(shù)中,而噪聲的能量卻分布于整個(gè)小波域內(nèi),因此經(jīng)過小波分解后得到小波變換系數(shù)m,信號(hào)的小波變換系數(shù)要大于噪聲的小波變換系數(shù),于是可以找到一個(gè)合適的數(shù)作為閾值(門限) ,當(dāng)m小于閾值時(shí),認(rèn)為這時(shí)的m主要是由噪聲引起的;當(dāng)m大于該閾值時(shí),認(rèn)為這時(shí)的m主要是由信號(hào)引起的,從而實(shí)現(xiàn)了信噪分離[12].

MATLAB將矩陣運(yùn)算、數(shù)值分析、圖形處理、編程技術(shù)等功能有機(jī)地結(jié)合在一起，具有強(qiáng)大的數(shù)學(xué)計(jì)算和分析功能、豐富的可視化圖形表現(xiàn)功能和方便的程序設(shè)計(jì)能力[13].

MATLAB中實(shí)現(xiàn)小波去噪功能的常用函數(shù)是wden[14],格式為:

[XD,CXD,LXD] = wden(X,TPTR,SORH,SCAL,N,′wname′)

式中：X是原始時(shí)間序列信號(hào)；TPTR表示閾值選擇的原則；SORH表示選擇軟閾值或硬閾值處理方法；SCAL表示小波分解時(shí)定義所用的閾值是否需要重新調(diào)整；N表示小波分解的層次；′wname′是所使用的正交小波名稱；XD、CXD、LXD分別是在上述參數(shù)作用下小波去噪后得到的信號(hào)及其附加小波分解結(jié)構(gòu).

2 MATLAB小波去噪函數(shù)參數(shù)及其選擇

在對(duì)PM10進(jìn)行小波分析時(shí),wden函數(shù)的參數(shù)選擇對(duì)去噪的結(jié)果有直接影響,現(xiàn)就不同參數(shù)對(duì)PM10時(shí)間序列去噪結(jié)果的影響分析如下：

2.1 閾值的幾種選擇原則

在MATLAB中有4 種閾值函數(shù)形式可以選用:rigrsure采用Stein的無偏似然估計(jì)原理進(jìn)行自適應(yīng)閾值選擇;sqtwolog表示采用固定的閾值形式; minimaxi表示采用極大極小原理選擇的閾值,它和sqtwolog一樣，也是一種固定的閾值, 它產(chǎn)生一個(gè)最小均方誤差的極值; heursure表示選擇啟發(fā)式閾值, 它是sqtwolog 和rigrsure 的綜合.具體的閾值公式可參考文獻(xiàn)[15].

rigrsure和minimaxi閾值選取原則較為保守(僅將部分系數(shù)置為零),可以將弱小的信號(hào)提取出來[15].通過對(duì)PM10的時(shí)間序列變化規(guī)律分析,并經(jīng)試驗(yàn)發(fā)現(xiàn),幾種選擇原則沒有太大區(qū)別.本文選取rigrsure作為閾值選取原則.

2.2 軟閾值和硬閾值

在對(duì)小波系數(shù)作門限閾值處理操作時(shí), 硬閾值處理只保留較大的小波系數(shù)并將較小的小波系數(shù)置零,軟閾值處理將較小的小波系數(shù)置零但對(duì)較大的小波系數(shù)向零作了收縮.經(jīng)試驗(yàn)發(fā)現(xiàn),本文分析的情況軟、硬閾值選擇沒有太大區(qū)別.本文選取軟閾值作為閾值處理方法.

2.3 閾值調(diào)整

SCAL參數(shù)定義所用的閾值是否需要重新調(diào)整,有3種情況：scal=′one′不進(jìn)行重新估計(jì); scal=′sln′只根據(jù)第一層小波分解系數(shù)估計(jì)噪聲水平;scal=′mln′在每個(gè)不同的小波分解層次估計(jì)噪聲水平. 圖1是對(duì)PM10進(jìn)行小波分解時(shí)用不同閾值重新調(diào)整方法得出的不同結(jié)果,使用函數(shù)為wden(X,′rigrsure′,′s′,SCAL,5,′sym8′).可以看出,對(duì)PM10進(jìn)行小波分解時(shí)應(yīng)該選擇′one′.

圖1 不同閾值重新調(diào)整方法對(duì)去噪結(jié)果的影響

2.4 分解層次的比較

圖2是采用不同分解層次的PM10小波去噪效果圖,使用函數(shù)為wden(X,′rigrsure′,′s′,′one′,level,′sym8′).從圖中可以看出，隨著小波分解層次的增加，去噪效果變好.但是分解層次增加到5層以上后去噪效果改善已經(jīng)不明顯.本研究分解層次取5層.

2.5 小波函數(shù)的選取

圖3是采用不同小波函數(shù)的PM10小波去噪效果圖,使用函數(shù)為wden(X,′rigrsure′,′s′,′one′,5,wname).小波函數(shù)(wname)可以從Daubechies (dbN)、symlets(symN)等常用的小波函數(shù)中選取一個(gè)正交小波,從圖中可以看出，小波函數(shù)Daubechies和symlets同層次的去噪效果基本相同,隨著小波函數(shù)序數(shù)的增加,去噪效果較好.本文選擇sym8作為去噪小波函數(shù).

3 實(shí)例分析

根據(jù)國(guó)家環(huán)境保護(hù)部網(wǎng)站公布的北京市2010-2011年度大氣空氣質(zhì)量監(jiān)測(cè)結(jié)果，將用污染指數(shù)表示的當(dāng)天的可吸入顆粒物的污染情況轉(zhuǎn)換為質(zhì)量濃度值,具體轉(zhuǎn)換方法見文獻(xiàn)[16]，其中部分缺省數(shù)據(jù)根據(jù)前后兩天線性插值得到.

將北京市2010-2011年度PM10質(zhì)量濃度數(shù)據(jù)導(dǎo)入MATLAB，對(duì)數(shù)據(jù)進(jìn)行小波去噪分析，北京市2010-2011年度PM10質(zhì)量濃度隨時(shí)間的變化曲線和小波分析結(jié)果如圖4所示.圖5是按月份平均的柱狀圖.可見，經(jīng)過小波分析去噪后，污染物的變化規(guī)律比較明顯地顯現(xiàn)出來了，小波去噪處理在分析時(shí)間序列數(shù)據(jù)中確實(shí)有明顯的效果.

從圖4中PM10濃度變化可以看出，北京市PM10濃度變化在一年中呈現(xiàn)一定的規(guī)律性，每年年初有一下降的過程，隨后又有上升，4、5月份開始下降，在7月份有一小峰值，隨后逐漸下降，一般到8、9月份達(dá)到一年中的低點(diǎn)，隨后逐漸上升到11月份達(dá)到較高值，以后逐漸下降到次年.分析可知，造成這種變化的原因,除了造成顆粒物污染的工業(yè)生產(chǎn)等活動(dòng)的季節(jié)性波動(dòng)外,也與影響顆粒物污染的氣候氣象條件及人類活動(dòng)有關(guān).每年的3-4月份，PM10濃度較高，這可能和春季的沙塵入侵有關(guān)，此時(shí)北京市氣溫逐漸升高,土壤解凍,風(fēng)力較大,植被尚未完全長(zhǎng)出,容易形成沙塵天氣,造成顆粒物濃度較高; 隨后風(fēng)力變小,植被完全恢復(fù),對(duì)抑制沙塵天氣形成有利; 7月份是北京市的雨季，天氣溫度較高,風(fēng)力較弱(表1所示的北京市1971-2000年的氣候值)，空氣的擴(kuò)散能力減小，造成污染物積累，同時(shí)大量排放的機(jī)動(dòng)車尾氣在高溫和強(qiáng)烈紫外線作用下發(fā)生二次反應(yīng)，會(huì)生成較多的細(xì)顆粒物，這些因素都有可能使顆粒物濃度增加[17]；9月份秋高氣爽,風(fēng)力增加，污染物得以擴(kuò)散,是一年中空氣質(zhì)量最好的月份；進(jìn)入10月份后，溫度降低，風(fēng)力較小,霧天和逆溫天氣出現(xiàn)較多,同時(shí)降水減少，由降水造成的顆粒物清除減少，造成污染物積累，大氣顆粒物濃度開始增加；到年底PM10濃度開始降低，主要原因可能是進(jìn)入冬季后，冷空氣頻繁入侵，使顆粒物得到擴(kuò)散、凈化，同時(shí)北京冬季天寒地凍,局地土壤不容易被風(fēng)吹起,使顆粒物濃度降低.

圖4 北京市2010-2011年度PM10濃度隨時(shí)間變化曲線和小波去噪結(jié)果

圖5 北京市2010-2011年度PM10月平均變化圖

表1 北京市標(biāo)準(zhǔn)氣候值(1971-2000年)

分析北京市大氣污染物質(zhì)量濃度的變化情況,不難看出，大氣污染物濃度的影響因素非常復(fù)雜，不同年份相同月份的變化有所差異，如2011年1月份PM10污染特別小，而2月份PM10較高,這主要是由于1月冷空氣頻繁,有利于大氣污染物的清除和擴(kuò)散，且恰遇農(nóng)歷春節(jié)放假,工業(yè)排放和機(jī)動(dòng)車尾氣排放減少,而2月16日至24日，出現(xiàn)了典型的靜風(fēng)穩(wěn)定天氣，逆溫較強(qiáng)，濕度上升，有時(shí)有霧，污染物易積累、難擴(kuò)散，空氣質(zhì)量轉(zhuǎn)差[17]，這與北京所處的地理位置和自然環(huán)境有關(guān)，也和氣候、氣象條件有關(guān).總的說來，北京市大氣質(zhì)量較好的季節(jié)在5-10月份.

4 結(jié) 語

(1)小波去噪分析方法能較好地用來分析PM10的時(shí)間序列變化規(guī)律，可以消除偶然因素對(duì)大氣顆粒物濃度變化的影響，更好地顯示PM10的時(shí)間變化規(guī)律.

(2)在運(yùn)用MATLAB提供的小波去噪分析函數(shù)時(shí)，不同參數(shù)對(duì)小波去噪的效果有不同的影響，應(yīng)該根據(jù)時(shí)間序列數(shù)據(jù)的特點(diǎn)進(jìn)行選擇，以達(dá)到較好的去噪效果.

(3)對(duì)北京市PM10質(zhì)量濃度數(shù)據(jù)進(jìn)行小波去噪分析表明，北京市PM10濃度在一年中變化呈現(xiàn)一定的規(guī)律性，造成這種變化的原因除了造成顆粒物污染的工業(yè)生產(chǎn)等活動(dòng)的季節(jié)性波動(dòng)外,也與影響顆粒物污染的氣候、氣象條件及人類活動(dòng)有關(guān).

參考文獻(xiàn)：

[1] Buseck P R，Posfai M. Airborne Minerals and Related Aerosol Particles: Effects on Climate and the Environment [J].P. Natl. Acad. Sci., 1999, 96:3372-3379.

[2] 孫杰，高慶先，周鎖銓.2002年北京PM10時(shí)間序列及其成因分析[J].環(huán)境科學(xué)研究, 2007, 20 (6) : 83-86.

[3] 王宏，林長(zhǎng)城，蔡義勇,等.福州市PM10突變特征與氣象條件的關(guān)系研究[J]. 熱帶氣象學(xué)報(bào), 2008, 24(5): 564-568.

[4] Li S T, Shu L Y. Data Mining to Aid Policy Making in Air Pollution Management[J]. Expert Systems with Applications, 2004, 27(3): 331-340.

[5] 張燕, 楊洋.基于小波分析的金融時(shí)間序列消噪方法及應(yīng)用[J].寧波大學(xué)學(xué)報(bào)(理工版), 2010, 23 (3): 56-59.

[6] 王海鵬，張斌，劉祖涵，等． 基于小波變換的蘭州市近十年空氣污染指數(shù)變化[J]． 環(huán)境科學(xué)學(xué)報(bào)，2011, 31(5): 1070-1076.

[7] 陳柳,馬廣大.小波變換在大氣污染物時(shí)間序列分析中的應(yīng)用[J].西安科技大學(xué)學(xué)報(bào), 2006, 26(1): 58-61.

[8] 錢舒. 小波在股市數(shù)據(jù)分析中的應(yīng)用[J]. 經(jīng)濟(jì)數(shù)學(xué), 2002, 19(4): 80-84.

[9] Sang Y, Wang D, Wu J, et al. Entropy-Based Wavelet De-noising Method for Time Series Analysis[J].Entropy, 2009, 11:1123-1147.

[10] De Gooijer J G，Hyndman R J. 25 Years of Time Series Forecasting[J]. International Journal of Forecasting, 2006, 22: 443-473.

[11] 張典,蔣勇軍,楊平恒,等.巖溶地下河水文時(shí)間序列小波消噪處理[J].人民長(zhǎng)江, 2009, 40(21): 45-46.

[12] 范曉志. 小波變換的信號(hào)去噪應(yīng)用[J]. 武漢科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2004, 27(3): 286-288.

[13] 李海濤，鄧櫻.MATLAB程序設(shè)計(jì)教程[M].北京：高等教育出版社，2002.

[14] 杜浩藩,叢爽.基于MATLAB 小波去噪方法的研究[J].計(jì)算機(jī)仿真, 2003, 20(7): 119-123.

[15] 飛思科技產(chǎn)品研發(fā)中心.MATLAB 6.5輔助小波分析與應(yīng)用[M].北京:電子工業(yè)出版社, 2003.

[16] 楊書申,邵龍義,楊園園.北京、上海兩地2004 和2005 年大氣污染特征對(duì)比分析[J].長(zhǎng)江流域資源與環(huán)境, 2008, 17(2): 0323-0327.

[17] 北京市環(huán)境保護(hù)局.北京市空氣質(zhì)量月報(bào)[EB/OL].[2013-03-10].http://www.bjepb.gov.cn/portal0/tab223 /info8206.htm.