， ，

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

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

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

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

信號處理時，一般認為低頻信號是有用信號而把高頻信號看作噪聲，要求有用信號和噪聲的頻譜相互分開[5,8-9].

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

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

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

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

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

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

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

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

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

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

2.1 閾值的幾種選擇原則

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

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

2.2 軟閾值和硬閾值

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

2.3 閾值調整

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

圖1 不同閾值重新調整方法對去噪結果的影響

2.4 分解層次的比較

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

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

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

3 實例分析

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

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

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

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

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

表1 北京市標準氣候值(1971-2000年)

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

4 結 語

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

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

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

參考文獻：

[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時間序列及其成因分析[J].環(huán)境科學研究, 2007, 20 (6) : 83-86.

[3] 王宏，林長城，蔡義勇,等.福州市PM10突變特征與氣象條件的關系研究[J]. 熱帶氣象學報, 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] 張燕, 楊洋.基于小波分析的金融時間序列消噪方法及應用[J].寧波大學學報(理工版), 2010, 23 (3): 56-59.

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

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

[8] 錢舒. 小波在股市數(shù)據(jù)分析中的應用[J]. 經濟數(shù)學, 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] 張典,蔣勇軍,楊平恒,等.巖溶地下河水文時間序列小波消噪處理[J].人民長江, 2009, 40(21): 45-46.

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

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

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

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

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

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