城市空氣污染程度的和預(yù)測模型

1、城市空氣污染程度的分析和預(yù)測摘 要 本文討論了有關(guān)城市污染程度、污染因素及污染擴(kuò)散的問題。對(duì)于問題一，本文主要從大氣污染、噪聲污染和水體污染這三個(gè)面選取主要污染物，查閱北京、天津、上海、重慶和西安五座城市2007-2012年的年度平均污染數(shù)據(jù)，采用降維的思想，運(yùn)用主成分分析法減少變量個(gè)數(shù)，再借助Matlab軟件計(jì)算各主成分的貢獻(xiàn)率，分析知可選取前三個(gè)主成分作為衡量污染程度的標(biāo)準(zhǔn)，最后根據(jù)綜合指標(biāo)得到這五個(gè)城市的污染程度從高到低依次為：重慶、上海、北京、天津、西安。通過判斷相關(guān)系數(shù)的大小，確定五個(gè)城市影響人們生活的主要污染因素是水污染，其四項(xiàng)指標(biāo)依次為化學(xué)需氧量、總氮、總磷和氨氮。對(duì)于問題二，以

2、北京市大氣污染為例。首先，利用GPS記錄北京市14個(gè)城區(qū)觀測點(diǎn)的位置，并查閱2013年污染指標(biāo)、與的污染數(shù)據(jù)，繪制出相應(yīng)的空間濃度分布圖，估計(jì)這四種污染物的大致污染源位置依次為：附近、附近，附近和附近；其次，根據(jù)污染擴(kuò)散原理和方式，建立Cauchy污染傳播模型，根據(jù)各地區(qū)空氣污染物的濃度分布，運(yùn)用Matlab軟件對(duì)數(shù)據(jù)非線性擬合，得出擴(kuò)散模型各參數(shù)的值，計(jì)算得出各項(xiàng)污染指標(biāo)的污染源位置依次為：，和；最后，比較污染物位置的計(jì)算值與實(shí)際值，發(fā)現(xiàn)誤差相差較小，故模型建立較為合理。對(duì)于問題三，分析西安市的主要污染大氣污染。收集西安市2014年4.1-7.31日的空氣污染數(shù)據(jù)，根據(jù)時(shí)間序列的平穩(wěn)性特點(diǎn)及

3、定階準(zhǔn)則選取合適的時(shí)間序列模型，利用Matlab軟件對(duì)序列模型的各項(xiàng)參數(shù)進(jìn)行估計(jì)并檢驗(yàn)?zāi)Ｐ偷暮侠硇?，并將模型用于?shù)據(jù)預(yù)報(bào)。利用時(shí)間序列模型預(yù)測西安市未來10天的空氣污染狀況總體等級(jí)為良。對(duì)于問題四，基于問題一、二、三對(duì)污染因素的分析和污染擴(kuò)散的特點(diǎn)，主要從減少污染物的產(chǎn)生和治理凈化已產(chǎn)生的污染物兩方面，針對(duì)大氣污染、水體污染和噪聲污染為相關(guān)部門提供合理化防治建議。關(guān)鍵詞 主成分分析；Cauchy污染傳播模型；時(shí)間序列模型；Matlab軟件30一、 問題重述環(huán)境質(zhì)量的好壞對(duì)人們的生活有著最直接的影響，包括身體健康、壽命長短、生活質(zhì)量以及經(jīng)濟(jì)發(fā)展水平、生態(tài)環(huán)境等。因此保護(hù)環(huán)境的工作需要個(gè)人、集體、

4、國家、乃至全球各國的共同努力。事實(shí)上，環(huán)境保護(hù)事業(yè)從1962年起步至今，一直備受人們關(guān)注和重視，環(huán)境污染的防治工作也一步步開展運(yùn)營，然而較經(jīng)濟(jì)發(fā)展的速度和水平，還有很多需要為之努力的工作要做。環(huán)境保護(hù)工作一方面是要從減少污染物的產(chǎn)生方面，預(yù)防和控制污染源頭污染物的排放，另一方面是治理凈化已產(chǎn)生的污染物。環(huán)境污染是指環(huán)境因受人類活動(dòng)影響而改變了原有性質(zhì)和狀態(tài)的現(xiàn)象，其實(shí)質(zhì)是人類活動(dòng)中將大量的污染物排入環(huán)境，影響其自凈能力，降低了生態(tài)系統(tǒng)的功能。為考察某時(shí)期某區(qū)域的環(huán)境狀況，需統(tǒng)計(jì)該地區(qū)各項(xiàng)污染指標(biāo)的數(shù)據(jù)，根據(jù)各污染指標(biāo)的污染排放標(biāo)準(zhǔn)，確定其污染程度，并通過污染發(fā)生的原因、途徑和后果制定合理的污染

5、治理方案，加以實(shí)施運(yùn)行，改善環(huán)境。統(tǒng)計(jì)北京、天津、上海、重慶和西安五座城市的污染數(shù)據(jù)，討論如下問題：問題一：比較五座城市的污染程度及影響人們生活的主要污染因素；問題二：建立城市污染擴(kuò)散模型確定某城市污染程度并檢驗(yàn)；問題三：給出西安市某方面污染狀況的中長期預(yù)測數(shù)據(jù)；問題四：給相關(guān)部門提供一篇治理污染的建議短文。二、 問題分析對(duì)于問題一，本文將查閱收集北京、上海、天津、重慶以及西安五座城市大氣污染、噪聲污染和水體污染2007-2012年的數(shù)據(jù)，由于各項(xiàng)污染指標(biāo)眾多，所以選取具有代表性的污染物：可吸入顆粒物、道路交通噪聲、區(qū)域環(huán)境噪聲、化學(xué)需氧量、氨氮、總氮和總磷這9項(xiàng)指標(biāo)來綜合衡量各城市的污染程度

6、，本文將運(yùn)用降維的思想，利用主成分分析法計(jì)算前幾個(gè)主成分的貢獻(xiàn)率大小及通過線性組后的綜合指標(biāo)，以此來衡量各城市的污染程度，根據(jù)各主成分的相關(guān)系數(shù)大小即可確定影響人們生活的主要污染因素(1)。對(duì)于問題二，本文將以北京市的大氣污染為例，收集北京市2013年14個(gè)城區(qū)的、和的污染數(shù)據(jù)，繪制出這四種污染物濃度的空間濃度分布圖，估算出這四種污染物的大致污染源位置；建立Cauchy污染傳播模型，根據(jù)各地區(qū)空氣污染物的濃度分布，運(yùn)用Matlab軟件對(duì)數(shù)據(jù)nlinfit非線性擬合，計(jì)算出擴(kuò)散模型各參數(shù)值，得出各空氣污染物的污染源位置；比較污染物位置的計(jì)算值與實(shí)際值，檢驗(yàn)污染傳播模型建立是否合理(2)。對(duì)于問題

7、三，本文將考慮西安市大氣污染的各項(xiàng)污染指標(biāo)，搜集2014年4月1日至2014年7月31日的污染數(shù)據(jù)，根據(jù)時(shí)間序列的平穩(wěn)性特點(diǎn)及定階準(zhǔn)則選取合適的時(shí)間序列模型，利用Matlab軟件(3)估計(jì)序列模型的各項(xiàng)參數(shù)并檢驗(yàn)?zāi)Ｐ偷暮侠硇?，考核通過時(shí)方可用于數(shù)據(jù)預(yù)報(bào)。對(duì)于問題四，基于問題一、二、三對(duì)污染因素的分析和污染擴(kuò)散的特點(diǎn)，主要從減少污染物的產(chǎn)生和治理凈化已產(chǎn)生的污染物兩方面，針對(duì)大氣污染、水體污染和噪聲污染為相關(guān)部門提供合理化防治建議。三、 模型假設(shè)1.所有城市污染數(shù)據(jù)均真實(shí)有效，具有統(tǒng)計(jì)價(jià)值；2.污染數(shù)據(jù)的取樣為瞬時(shí)取樣且在某個(gè)平面內(nèi)，不考慮取樣前后污染物的濃度變化；3.繪制北京大氣污染濃度分布圖

8、時(shí)不考慮海拔因素的影響；4.大氣污染過程經(jīng)歷著擴(kuò)散與衰減，且沿軸方向上的擴(kuò)散系數(shù)分別為常數(shù)，衰減與濃度成正比；5.不考慮日照等自然因素對(duì)各項(xiàng)指標(biāo)污染濃度的影響。四、 符號(hào)說明4.1 模型一： 可吸入顆粒物（單位）；：二氧化硫（單位）；： 二氧化氮（單位）；： 道路交通噪聲（單位分貝）；： 區(qū)域環(huán)境噪聲（單位分貝）；： 化學(xué)需氧量（單位萬噸）；： 氨氮（單位萬噸）；： 總氮（單位萬噸）；： 總磷（單位萬噸） ；： 樣本協(xié)方差矩陣 ；： 協(xié)方差矩陣的特征值 ；： 協(xié)方差矩陣的特征向量； ： 相關(guān)系數(shù)；： 樣本方差向量 ；： 樣本相關(guān)矩陣；： 樣本相關(guān)矩陣的特征值；： 樣本相關(guān)矩陣的特征向量；： 樣

9、本均值向量；： 主成分；： 污染程度綜合指標(biāo)。(其中）4.2模型二： 位置坐標(biāo)/橫軸；： 位置坐標(biāo)/縱軸；： 位置坐標(biāo)/豎軸； ： 污染物濃度；： 污染物質(zhì)量； ： 污染物排放量；： 擴(kuò)散系數(shù)。4.3模型三： 序列；：平穩(wěn)白噪聲；： 滑動(dòng)平均序列；： 自回歸序列；: 自回歸滑動(dòng)平均序列 ；： 常數(shù)；： 自協(xié)方差函數(shù)；： 自相關(guān)函數(shù)；： 概率密度；： 未知參數(shù)向量: 最大似然估計(jì)；： 自回歸參數(shù)。五、 模型建立與求解5.1 問題一的模型建立與求解環(huán)境因受人類活動(dòng)影響而改變了原有性質(zhì)和狀態(tài)的現(xiàn)象稱為環(huán)境污染。環(huán)境污染的實(shí)質(zhì)是人類活動(dòng)中將大量的污染物排入環(huán)境，影響其自凈能力，降低了生態(tài)系統(tǒng)的功能。然

10、而生活中環(huán)境的污染也是多方面的，污染按環(huán)境要素分類，可分為：大氣污染、土壤污染和水體污染。為此，本文選取環(huán)境污染的主要污染影響指標(biāo)，建立污染評(píng)價(jià)模型，即主成分分析模型，比較北京、天津、上海、重慶和西安五座城市的污染程度，并得出影響人們生活的主要污染因素。5.1.1 建立主成分分析模型首先選取環(huán)境污染的部分污染指標(biāo)：可吸入顆粒物、二氧化硫、二氧化氮、道路交通噪聲、區(qū)域環(huán)境噪聲、化學(xué)需氧量、氨氮、總氮和總磷。北京、天津、上海、重慶和西安五座城市的各項(xiàng)指標(biāo)年平均數(shù)據(jù)（見附錄表1）。問題中變量太多不但會(huì)增加計(jì)算的復(fù)雜性，也給合理的分析問題和解決問題帶來很大困難；同時(shí)，這些變量之間存在的相關(guān)性，也使得這

11、些變量所反映的信息在一定程度上有所重疊。為了減少變量個(gè)數(shù)，本文采用了降維的思想，利用主成分分析法得出污染空氣的主要因素變量，提高問題研究的合理性，同時(shí)不會(huì)使數(shù)據(jù)反映的信息量有大的損失。設(shè)是的協(xié)方差矩陣，的特征值與正交化特征向量分別為及，且的第個(gè)主成分為 （1）根據(jù)平均值計(jì)算公式代入數(shù)據(jù)計(jì)算樣本的均值向量為根據(jù)協(xié)方差矩陣計(jì)算公式 （2）利用Matlab軟件代入數(shù)據(jù)可求得隨機(jī)變量相應(yīng)樣本協(xié)方差矩陣為（只寫下三角部分） 由于協(xié)方差矩陣對(duì)角線元素即為樣本的方差值，觀察矩陣易知其方差值相差較大，即各量綱間分散程度較大，不適用于環(huán)境污染模型中各成分指標(biāo)間的相關(guān)性分析，為消除由于量綱不同可能帶來的結(jié)果誤差，

12、故將樣本變量標(biāo)準(zhǔn)化，即令 （3）得到一組新的數(shù)據(jù)。其中為的平均值，為的方差。此時(shí)，的協(xié)方差矩陣即為樣本變量的相關(guān)矩陣,其中 （4） 為的協(xié)方差。則根據(jù)公式（4）運(yùn)用Matlab軟件代入數(shù)據(jù)即可得其樣本相關(guān)矩陣為由特征值計(jì)算公式 （5）第個(gè)主成分的貢獻(xiàn)率 （6）及前個(gè)主成分的累計(jì)貢獻(xiàn)率 （7）其中。計(jì)算的特征值及各標(biāo)準(zhǔn)化主成分的貢獻(xiàn)率和累計(jì)貢獻(xiàn)率如表1所示:表1 的特征值、貢獻(xiàn)率和累計(jì)貢獻(xiàn)率貢獻(xiàn)率(%)累計(jì)貢獻(xiàn)率(%)14.742052.6952.6922.588428.7681.4531.476816.4197.8640.19282.141005001006001007001008001009

13、00100由表1易知，第一、第二和第三主成分的累積貢獻(xiàn)率為97.86%，并且只有第一、第二和和第三主成分對(duì)應(yīng)的特征值大于1，所以選取第一主成分和第二主成分第三主成分作為標(biāo)準(zhǔn)化數(shù)據(jù)的主成分討論問題，剩余的主成分對(duì)本文問題的討論影響不大，因此可以忽略。由此可得的特征值所對(duì)應(yīng)的前三個(gè)正交單位化特征向量為：記第個(gè)主成分為,則標(biāo)準(zhǔn)化數(shù)據(jù)的第一主成分、第二主成分、第三主成分第九主成分。由與的相關(guān)系數(shù) （8）利用Matlab軟件代入數(shù)據(jù)計(jì)算出前三個(gè)主成分與各原始變量的相關(guān)系數(shù)如表2所示 ：表2 相關(guān)系數(shù)相關(guān)系數(shù)-0.3593-0.8135-0.40650.1677-0.93010.0655-0.65380.

14、4274-0.6164-0.67550.6799-0.2789-0.7407-0.19110.62780.96770.21370.13370.73410.56250.38020.9355-0.1263-0.32650.8955-0.0454-0.4317記樣本均值向量為，則再利用方差計(jì)算公式代入數(shù)據(jù)可以得到樣本向量的方差向量為：所以標(biāo)準(zhǔn)化數(shù)據(jù)的第一主成分、第二主成分和第三主成分可以用樣本數(shù)據(jù)表示為：由于、和特征值的大小決定它們貢獻(xiàn)率的大小，所以結(jié)合、和的特征值線性組合后可得衡量各地區(qū)消費(fèi)支出的綜合指標(biāo)： （9）把、和的值代入（9）式即可得到，按照綜合指標(biāo)的大小便可對(duì)5個(gè)城市進(jìn)行排序，從而得知各

15、個(gè)城市的污染程度：表3 三個(gè)主成分值、綜合指標(biāo)值及排名排名北京-0.75230.63240.6325-0.99683天津0.3239-1.3453-1.3453-3.93304上海-0.81752.40652.40655.90622重慶3.5358-0.1692-0.169216.03131西安-2.2798-1.5244-1.5244-17.00775綜合考慮大氣污染、水體污染和噪聲污染的各項(xiàng)指標(biāo)：可吸入顆粒物、二氧化硫、二氧化氮、道路交通噪聲、區(qū)域環(huán)境噪聲、化學(xué)需氧量、氨氮、總氮和總磷，五座城市中重慶市污染程度最為嚴(yán)重，上海次之，北京、天津污染程度較為接近，西安市污染程度較輕。5.1.2

16、結(jié)論分析通過對(duì)北京、天津、上海、重慶和西安五座城市各項(xiàng)指標(biāo)年平均數(shù)據(jù)統(tǒng)計(jì)和計(jì)算，得出結(jié)論：由于標(biāo)準(zhǔn)化后的第一主成分的貢獻(xiàn)率已達(dá)52.69%，大于其他主成分的貢獻(xiàn)率，可以認(rèn)為值越大，環(huán)境質(zhì)量越差。因此觀察與各項(xiàng)指標(biāo)間的相關(guān)系數(shù)可知，水體污染的四項(xiàng)指標(biāo)化學(xué)需氧量、總氮、總磷和氨氮對(duì)環(huán)境污染的影響較大，其中化學(xué)需氧量對(duì)人們生活環(huán)境造成的污染程度最大。化學(xué)需氧量是衡量水質(zhì)受有機(jī)物污染程度的綜合指標(biāo)，當(dāng)其值過高時(shí)，水體會(huì)發(fā)黑、發(fā)臭?？偟欠从乘w富營養(yǎng)化的主要指標(biāo)；過量磷可造成水體污穢異臭，湖泊發(fā)生富營養(yǎng)化和海灣赤潮現(xiàn)象；飲用水中氨氮過量在一定條件下可轉(zhuǎn)化為亞硝酸鹽，是強(qiáng)致癌物質(zhì)，對(duì)人體健康極為不利；水

17、體中氨氮過量則可導(dǎo)致水生物表現(xiàn)亢奮甚至死亡。綜合考慮大氣污染、水體污染和噪聲污染的各項(xiàng)指標(biāo)，得出結(jié)論：五座城市的總體污染程度依次為：重慶、上海、北京、天津、西安。重慶市由于三面環(huán)山的特殊地理位置和重工業(yè)發(fā)達(dá)、煤煙排放過量等原因，使得其總體污染程度在五座城市中居于首位；上海機(jī)動(dòng)車尾氣排放量常年居高不下，污染程度僅次于重慶；北京和天津地理位置接近，其工業(yè)污染物排放量較高，造成大量的霧霾揚(yáng)塵天氣；而西安汽車尾氣和工業(yè)廢氣排放量較大，因此空氣污染最為嚴(yán)重，但污染程度較其他城市略為輕微。5.2 問題二的模型建立與求解由于人類活動(dòng)的影響，大氣、水體及噪聲等污染物濃度及其分布發(fā)生了很大變化，了解各污染物的空

18、間分布及擴(kuò)散情況對(duì)評(píng)價(jià)其污染程度有很大幫助。一方面，確定污染擴(kuò)散的方式、速率和程度有助于控制其影響范圍，提早預(yù)防和改善，另一方面可確定污染源的位置，從根本上杜絕或減輕污染的發(fā)生和污染程度，為相關(guān)部門提供污染治理的意見和建議。以北京市為例，建立其污染擴(kuò)散模型并以空氣污染為檢驗(yàn)標(biāo)準(zhǔn)驗(yàn)證模型的合理性。5.2.1 估計(jì)北京市污染源位置查閱并搜集北京市污染空氣的、和（可吸入顆粒）的年平均濃度含量，取包含北京市在內(nèi)的長為183.4千米，寬為170.1千米的長方形區(qū)域（如圖1），建立二維坐標(biāo)系，得到北京市14個(gè)區(qū)的地理位置坐標(biāo)(單位：千米)如表5所示:圖1 北京市各城區(qū)圖表4 2013年北京市各項(xiàng)污染指標(biāo)值

19、及地理位置坐標(biāo)單位（/萬噸）橫坐標(biāo)縱坐標(biāo)朝陽區(qū)29.764112.491.3132.191.9海淀區(qū)26.963.611598.1126.576.8豐臺(tái)區(qū)28.157.5118.596.9135.373.7石景山區(qū)24.963.3116.492.8130.166.5門頭溝區(qū)24.651.8114.891.1128.361.2房山區(qū)31.261.9131.7106.8148.562.3通州區(qū)38.655.8123.5105.7130.7101.4順義區(qū)20.844.898.584.8104.7104.3大興區(qū)33.765.7130.3107.8149.778.8昌平區(qū)25.943.594.779

20、.294.969.3平谷區(qū)20.63598.784.8101.4145.4懷柔區(qū)22.337.995.376.183.1103.1密云縣21.343.685.971.676.9117.4延慶縣19.234.478.368.066.547.1根據(jù)上表數(shù)據(jù)利用Matlab軟件分別畫出、和的年平均濃度隨地理位置變化的空間濃度分布圖：圖2 的空間濃度分布圖 圖3 的空間濃度分布圖圖4 的空間濃度分布圖 圖5 的空間濃度分布圖觀察圖2、圖3、圖4和圖5估計(jì)污染指標(biāo)的污染源位置在附近，的污染源位置在附近，可吸入顆粒污染源位置在附近，細(xì)微顆粒的污染源位置在附近 。5.2.2 建立Cauchy擴(kuò)散模型設(shè)表示時(shí)

21、刻處某一種污染物的濃度。 在污染物傳播的區(qū)域內(nèi)，任意取一封閉曲面，它所圍區(qū)域?yàn)椋跁r(shí)間內(nèi)，由于擴(kuò)散，經(jīng)過這個(gè)區(qū)域內(nèi)的污染物的質(zhì)量為 （10）其中，分別代表沿軸方向的擴(kuò)散系數(shù)，是方向余弦。由Gauss-Green公式得擴(kuò)散過程中由于沉降等原因?qū)е聟^(qū)域內(nèi)質(zhì)量減少量為 （11）其中表示衰減系數(shù)。 時(shí)間區(qū)域內(nèi)污染物質(zhì)量的變化為時(shí)間區(qū)域內(nèi)污染物質(zhì)量的變化也等于這段時(shí)間內(nèi)增加的污染物質(zhì)量減去減少的污染物的質(zhì)量，即由于所選時(shí)間、區(qū)域的任意性，可以得到 (12)記污染源在點(diǎn)處，則在時(shí)刻有其中為污染物排放量。由以上分析可知此擴(kuò)散問題滿足Cauchy問題。 (13)對(duì)(13)進(jìn)行傅里葉(Fourier)變換，且令

22、，已知 ，故得常微分方程Cauchy問題唯一解 (14)對(duì)(14)求逆變換由于故得 為了簡化模型，考慮空氣污染的平面?zhèn)鞑ィ瑒t可得到不同空氣污染指標(biāo)的擴(kuò)散方程: (15)如果認(rèn)為經(jīng)過了相當(dāng)長的時(shí)間后，擴(kuò)散已經(jīng)終止，物質(zhì)分布處于平穩(wěn)狀態(tài)，則方程(12)中的，于是有線性橢圓方程的邊值問題也就是用傅里葉(Fourier)變換求解，其中是區(qū)域的邊界. 在Cauchy問題的解(14)中，有三個(gè)未知的參數(shù)，它們分別是擴(kuò)散與衰減過程中的擴(kuò)散系數(shù)與衰減系數(shù)的算術(shù)平方根，假設(shè)點(diǎn)源的質(zhì)量與位置已知，設(shè)觀測取樣為：取樣時(shí)刻為，或設(shè)，是取樣時(shí)間，則(12)變成：對(duì)而言，取樣時(shí)間為1，而方程形狀與(12)一致，把在點(diǎn)觀測

23、到物質(zhì)密度與公式(15)都取對(duì)數(shù)，令則有 (16)令，則(16)可以寫成。當(dāng)時(shí)，有。利用Matlab軟件對(duì)傳播模型的各項(xiàng)參數(shù)進(jìn)行估計(jì)并用于檢驗(yàn)?zāi)Ｐ偷暮侠硇裕】諝馕廴镜娜狱c(diǎn)數(shù)據(jù),運(yùn)用非線性擬合nlinfit函數(shù)對(duì)數(shù)據(jù)進(jìn)行擬合，得出函數(shù)參數(shù)即可得污染源位置。Matlab求得的空氣污染指標(biāo)污染源位置為，污染源位置為，可吸入顆粒的污染源位置為，的污染源位置為，與實(shí)際污染源位置很接近，故Cauchy污染傳播模型建立合理。5.3 問題三的模型建立與求解根據(jù)相關(guān)數(shù)據(jù)西安市在某個(gè)時(shí)間段內(nèi)PM2.5一直高居全國前列，因此空氣污染是西安市人民關(guān)注的主要方面。采用時(shí)間序列模型對(duì)先市空氣污染的各項(xiàng)指標(biāo)進(jìn)行預(yù)測。通

24、過平穩(wěn)性分析和模型定階的準(zhǔn)則確定預(yù)測模型，再利用最大似然法估計(jì)各級(jí)參數(shù)，最后進(jìn)行模型預(yù)報(bào)。5.3.1時(shí)間序列模型的建立1.平穩(wěn)性檢驗(yàn)（1）自相關(guān)函數(shù)序列（滑動(dòng)平均序列）的自相關(guān)函數(shù)：設(shè)是序列，。其中是零均值，即平穩(wěn)白噪聲，且，則是零均值平穩(wěn)序列，且其自協(xié)方差函數(shù)和自相關(guān)函數(shù)分別為 (17) (18)其自相關(guān)函數(shù)呈拖尾狀趨于零，偏自相關(guān)函數(shù)式非截尾但趨于零。序列（自回歸序列）的自相關(guān)函數(shù)：設(shè)是序列：，是平穩(wěn)白噪聲，且。則其自協(xié)方差相關(guān)函數(shù)和自相關(guān)函數(shù)為 (19) (20)其自相關(guān)函數(shù)呈拖尾狀趨于零，偏自相關(guān)函數(shù)非截尾但趨于零。而的自相關(guān)函數(shù)可以看做是序列的自相關(guān)函數(shù)和序列的自相關(guān)函數(shù)的混合物。計(jì)

25、算知本時(shí)間序列服從模型，下面考慮定階，即估計(jì)的值。（2）模型定階的準(zhǔn)則準(zhǔn)則又稱Akaile信息準(zhǔn)則，AIC準(zhǔn)則是信息論與統(tǒng)計(jì)學(xué)的重要研究成果，具有重要意義。準(zhǔn)則起源于Kullback-Lelibler信息量，設(shè)是隨機(jī)變量，它的概率密度是，其中含有個(gè)未知參數(shù)，設(shè)未知參數(shù)向量為。屬于分布族，其中。顯然。K-L信息量可以用來刻畫與的接近程度，其定義為則有，且有。K-L信息量是尋求最接近的參數(shù)概率密度，使得 （21）經(jīng)過理論分析，當(dāng)給定樣本觀測值（樣本容量為），設(shè)是模型參數(shù)（未知參數(shù)個(gè)數(shù)是）的最大似然估計(jì)。設(shè)是其對(duì)數(shù)似然函數(shù)，信息準(zhǔn)則是：使得（21）式中的滿足設(shè)是序列，其中未知參數(shù)的個(gè)數(shù)是個(gè)，包括自回

26、歸參數(shù)，滑動(dòng)平均參數(shù)及。對(duì)于對(duì)數(shù)似然函數(shù) （22）略去，得最小平方和估計(jì)對(duì)應(yīng)似然函數(shù) （23）由（23）式可得，的最大似然估計(jì) 將上式代入（23），得 因此，序列定階準(zhǔn)則為：選，使得 （24）其中是樣本容量，與有關(guān)。若當(dāng)，使得（24）達(dá)到最小值，則認(rèn)為序列為。當(dāng)序列含有未知均值時(shí)，模型為這時(shí)，未知參數(shù)個(gè)數(shù)為，準(zhǔn)則為：選取，使得 （25）其式（24）與式(25）有相同的最小值點(diǎn)。2.最大似然估計(jì)法（ML估計(jì)）估計(jì)參數(shù)假設(shè)是零均值正態(tài)序列，他的概率密度為其中，,表示參數(shù)向量是時(shí)取期望（因是序列，與有關(guān)）。令，則有 對(duì)數(shù)似然函數(shù)其中為樣本觀測值。下證僅與有關(guān)，而與無關(guān)，事實(shí)上的第行第列的元素是,設(shè)的

27、傳遞形式是，則。由此即知與無關(guān)，從而僅與有關(guān)，而與無關(guān)。因?yàn)閮H與有關(guān)，記 （26）則稱為平方和函數(shù)，進(jìn)一步記，即是關(guān)于的條件期望，也就是已知時(shí)關(guān)于的估計(jì)，則可證可見是一殘差平方和形式，若使（26）式達(dá)到最小，即，則稱為最小平方和估計(jì)?，F(xiàn)討論最大似然估計(jì)，對(duì)數(shù)似然函數(shù)是（27）因?yàn)?，均與無關(guān)，令解得 （28）將其帶入（27）中得 其中，略去常數(shù)，得 （29）的最大似然估計(jì)應(yīng)使（29）達(dá)到最大。若使得，則是的最大似然估計(jì)（ML估計(jì)），代入（28）即得的最大似然估計(jì)值。3.模型考核的檢驗(yàn) 若擬合模型的殘差記為，它是的估計(jì)值。例如對(duì)序列，設(shè)未知參數(shù)的估計(jì)是，則殘差記 Ljung-Box的檢驗(yàn)統(tǒng)計(jì)量是

28、假設(shè)檢驗(yàn)。在成立時(shí)，若充分大，近似于分布，其中是估計(jì)的模型參數(shù)個(gè)數(shù)。檢驗(yàn)法：給定顯著水平，設(shè)實(shí)際算得的值是，值為，則當(dāng)，拒絕，即認(rèn)為非白噪聲，模型考核不通過，需重新估計(jì)參數(shù)；而當(dāng)，接受，認(rèn)為是白噪聲，模型考核通過，可以進(jìn)行預(yù)報(bào)。4.的預(yù)報(bào)時(shí)間序列的步預(yù)報(bào)，是根據(jù)是取值對(duì)未來時(shí)刻的隨機(jī)變量做出估計(jì)，估計(jì)量記作，它是的線性組合。是零均值正態(tài)序列，則下列預(yù)報(bào)差分方程成立 （30）只需知道，就可以遞推算得，因此定義預(yù)報(bào)變量，令可證下列遞推預(yù)報(bào)公式：式中第三項(xiàng)當(dāng)時(shí)為零。且當(dāng)較小時(shí)，令初值。5.3.2模型求解及結(jié)論考慮數(shù)據(jù)的連貫性，選取2014.4.1日至2014.7.31日的空氣質(zhì)量指標(biāo)數(shù)據(jù)(見附錄表2

29、)進(jìn)行時(shí)間序列分析，運(yùn)用Matlab軟件作出相應(yīng)時(shí)序圖，進(jìn)行時(shí)間序列的驗(yàn)證：圖6 空氣質(zhì)量時(shí)序圖以及相應(yīng)的自相關(guān)系數(shù)圖與偏自相關(guān)系數(shù)圖： 圖7 自相關(guān)系數(shù)圖 表8 偏自相關(guān)系數(shù)圖觀察圖像明顯可知該序列不平穩(wěn)，故作差分處理。圖9為作一階差分后所得序列圖。圖9 一階差分后的空氣質(zhì)量時(shí)序圖為進(jìn)一步驗(yàn)證平穩(wěn)性，考察差分后序列自相關(guān)圖和偏自相關(guān)圖。圖10 一階差分后自相關(guān)圖 圖11 一階差分后偏自相關(guān)圖觀察圖像可知一階差分后的序列較平穩(wěn)，故可確定模型。利用準(zhǔn)則定階：當(dāng)時(shí)，計(jì)算得到值如下表所示：表5 值值1146.961128.281123.941120.141139.551115.621117.4811

30、18.161129.591117.561118.421120.421129.081117.291119.281119.83由表中數(shù)據(jù)知當(dāng)，時(shí)，其值最小，故取模型。又由檢驗(yàn)法，取顯著性水平,運(yùn)用Matlab軟件計(jì)算得到,故該模型通過考核檢驗(yàn)，可以用來預(yù)報(bào)。進(jìn)一步計(jì)算相關(guān)參數(shù)得到，故得到模型為其十步預(yù)報(bào)值見下表：表6 8月1-10日空氣質(zhì)量預(yù)測值與實(shí)際值比較8月1日2日3日4日5日6日7日8日9日10日預(yù)報(bào)值85.9274.4074.2768.2570.8666.2369.5265.2568.7064.54真實(shí)值12311311293797447396383進(jìn)一步對(duì)8月1-10日各項(xiàng)污染指標(biāo)、和的

31、濃度進(jìn)行預(yù)測，預(yù)測值如下表：表7 8月1-10日各項(xiàng)污染指標(biāo)的預(yù)測值（單位：）8月1日59.02117.978.4829.811.178月2日44.7999.238.0929.611.138月3日39.2388.467.7929.411.108月4日35.8882.127.5529.201.078月5日34.1078.267.3729.001.058月6日33.0875.797.15281.038月7日32.4274.086.9728.791.018月8日31.9472.806.8128.5918月9日31.5571.766.6528.380.988月10日31.1970.856.4927.9

32、70.97查閱8月1-10日西安市各項(xiàng)污染指標(biāo)、和的實(shí)際觀測數(shù)據(jù)，發(fā)現(xiàn)與本模型預(yù)測值之間存在誤差，當(dāng)誤差較小，故可認(rèn)為本模型采用模型進(jìn)行預(yù)測具有實(shí)際參考價(jià)值。5.4 污染治理建議地球是人類長期賴以生存的家園，然而目前嚴(yán)峻的環(huán)境污染態(tài)勢使得污染防治工作迫在眉睫，環(huán)境的長期污染不僅危害人類的健康，也對(duì)各動(dòng)植物物的生存繁衍造成威脅。而環(huán)境的保護(hù)工作不僅要對(duì)已產(chǎn)生的污染進(jìn)行治理，還應(yīng)從根本上預(yù)防污染的發(fā)生。因此，本文就空氣污染、水體污染和噪聲污染的預(yù)防和治理手段給出如下建議：環(huán)境污染的防治工作應(yīng)主要從減少污染物的產(chǎn)生和治理凈化已產(chǎn)生的污染物兩方面著手。空氣污染：一方面政府通過出臺(tái)相關(guān)的空氣污染指標(biāo)排放

33、標(biāo)準(zhǔn)和檢查制度，嚴(yán)格檢測控制各工業(yè)污染區(qū)及汽車尾氣、城鎮(zhèn)住宅區(qū)等各項(xiàng)指標(biāo)的排放量和區(qū)域密度，布局合理的工業(yè)區(qū)，調(diào)整產(chǎn)業(yè)結(jié)構(gòu)，整治燃煤企業(yè)，緊抓“禁燃區(qū)”工作，控制污染源；此外還應(yīng)調(diào)整能源結(jié)構(gòu)，改革生產(chǎn)工藝，采用無污染或低污染能源和工藝，如太陽能、風(fēng)能，并合理利用能源，對(duì)燃料進(jìn)行預(yù)處理，改進(jìn)燃燒設(shè)備，提高燃燒效率。另一方面加強(qiáng)廢氣的循環(huán)及回收利用，使用凈化裝置去除煙塵和各種工業(yè)粉塵，利用物理、化學(xué)方法凈化尾氣中的有害成分，引入潔凈煤技術(shù)。此外還應(yīng)加大環(huán)境保護(hù)的宣傳力度，提高公民素質(zhì)和環(huán)保意識(shí)。水體污染：其首要任務(wù)是保護(hù)水資源，調(diào)整缺水地區(qū)的產(chǎn)業(yè)結(jié)構(gòu)，控制高耗水、高耗能和重污染的建設(shè)項(xiàng)目，其次，資

34、源的開發(fā)利用堅(jiān)持開源節(jié)流并舉的方針，完善工農(nóng)業(yè)用水管理措施，建立用水考核制度，結(jié)合低資源消耗、高效益的高新技術(shù)，推行清潔生產(chǎn)。此外還應(yīng)建立循環(huán)用水制度和污水處理系統(tǒng)，尤其是有毒污染物的排放，建立資源更新補(bǔ)償機(jī)制，提高排污費(fèi)標(biāo)準(zhǔn)，使之高于污染治理成本等。噪聲污染：嚴(yán)格合理的控制機(jī)動(dòng)車數(shù)量、流量及汽車?yán)鹊镍Q笛，加強(qiáng)城市道路建設(shè)（多孔隙瀝青路面可降噪），施工區(qū)也應(yīng)注意施工時(shí)間和噪聲標(biāo)準(zhǔn)。設(shè)置部分聲屏障設(shè)施、綠化均能降噪。六、 模型評(píng)價(jià)與推廣6.1 模型評(píng)價(jià)：6.1.1模型一（主成分分析模型）優(yōu)點(diǎn):(1)本文采用降維的思想，用較少的變量反映原始變量，簡化模型的同時(shí)結(jié)合了其各項(xiàng)指標(biāo)的綜合權(quán)重，對(duì)考察城

35、市污染程度具有實(shí)際意義；(2)本文選擇先主成分分析后評(píng)價(jià)污染程度的思想，通過忽略一些次要因素的影響，考察某一綜合指標(biāo)的實(shí)際意義以便衡量五座城市的污染程度，具有現(xiàn)實(shí)意義。缺點(diǎn)：(1)由于所需土壤污染各指標(biāo)數(shù)據(jù)的缺少，比較五座城市的污染程度時(shí)未考慮土壤污染，可能使得主要污染因素與實(shí)際略有偏差或不足；(2) 降維后的新變量只是反映原始變量所提供的絕大部分信息，在實(shí)際問題的研究中，不能以此作為衡量城市污染程度的唯一量度。6.1.2模型二（污染傳播模型）優(yōu)點(diǎn)：(1)采用定性和定量結(jié)合的方法確定并檢驗(yàn)污染源位置，增強(qiáng)了模型的合理性；(2)各污染指標(biāo)的空間濃度分布圖較直觀的反映了其指標(biāo)污染程度集中區(qū)，使估計(jì)

36、的污染源位置更接近實(shí)際位置。缺點(diǎn)：(1)由于忽略空氣的縱向擴(kuò)散和風(fēng)速等對(duì)污染傳播的影響，使得計(jì)算結(jié)果與實(shí)際測量值有些出入；(2)由于所需水體污染和其他污染的數(shù)據(jù)缺少，只考慮了空氣污染的擴(kuò)散模型，對(duì)考察整體環(huán)境質(zhì)量不能提供完善的數(shù)據(jù)支持；(3)由于人的主觀影響，北京市地理位置坐標(biāo)及各項(xiàng)污染指標(biāo)的污染源估計(jì)值與實(shí)際存在誤差。6.1.3模型三（時(shí)間序列模型）優(yōu)點(diǎn)：(1) 對(duì)具有時(shí)序性、隨機(jī)性、前后時(shí)刻具有相依性、呈現(xiàn)某種趨勢或周期性的數(shù)據(jù)序列做了很好的處理，并做出了較準(zhǔn)確的預(yù)測；(2)生活中按時(shí)間順序排列、隨時(shí)間變化且相互關(guān)聯(lián)的例子較多，如工程，經(jīng)濟(jì)等，因此，本模型為其他部分領(lǐng)域的預(yù)測提供了參考。缺

37、點(diǎn)：(1)未考慮某些污染指標(biāo)可能存在的季節(jié)性和周期性變化，只運(yùn)用平穩(wěn)時(shí)間序列對(duì)各項(xiàng)污染指標(biāo)進(jìn)行大致預(yù)測，預(yù)測值與實(shí)際值之間可能存在誤差；(2)時(shí)間序列模型運(yùn)用不成熟，未能對(duì)模型殘差進(jìn)行精確分析和異常數(shù)據(jù)的剔除。6.2 模型推廣考慮到西安市每年11月15日至次年3月15日的冬日大面積供暖設(shè)備可能存在某些排放指標(biāo)過高現(xiàn)象，故考慮空氣重可吸入顆粒物的季節(jié)性變化，建立季節(jié)性序列模型：考慮時(shí)間序列，則階差分為，其中 程為階差分算子。則稱為序列。當(dāng)通過平穩(wěn)性檢驗(yàn)時(shí)便可建立乘積型季節(jié)性模型：進(jìn)一步利用最大似然估計(jì)（ML估計(jì)法）估計(jì)其相應(yīng)參數(shù)，通過檢驗(yàn)的考核通過時(shí)便可進(jìn)行預(yù)報(bào)。另外，本模型對(duì)同類型：具有時(shí)序性

38、、隨機(jī)性、前后時(shí)刻具有相依性、呈現(xiàn)某種趨勢或周期性的數(shù)據(jù)序列，如天文、電力、經(jīng)濟(jì)等其他領(lǐng)域的行業(yè)預(yù)測也具有一定參考價(jià)值。參考文獻(xiàn)1楊啟帆等.數(shù)學(xué)建模M.杭州：浙江大學(xué)出版社.20062石辛民等.計(jì)算方法及Matlab實(shí)現(xiàn)M.北京：清華大學(xué)出版社.20133范金城等.數(shù)據(jù)分析M.北京：科學(xué)出版社.2000附錄一、 主要程序1.主成分分析模型x=data;n=corrcoef(x);j,k=eig(n);p1=sqrt(k(1,1)*j(:,1)p2=sqrt(k(2,2)*j(:,2)p3=sqrt(k(3,3)*j(:,3) %計(jì)算相關(guān)矩陣和相關(guān)系數(shù)d=data;f=size(d,2)+1;m

39、=1;while mfx1=d(:,m);a=mean(x1);b=std(x1);n=length(x1)+1;c=1;while cnx1(c)=(x1(c)-a)/b;c=c+1;endd(:,m)=x1;m=m+1;endx=d; %標(biāo)準(zhǔn)化變量y1=x*j(:,1)y2=x*j(:,2)y3=x*j(:,2)y=k(1,1)*y1+k(2,2)*y2+k(3,3)*y3 %計(jì)算綜合指標(biāo)2.Cauchy污染傳播模型x=132.1 126.5 135.3 130.1 128.3 148.5 130.7 104.7 149.7 94.9 101.4 83.1 76.9 66.5;y=91.9

40、 76.8 73.7 66.5 61.2 62.3 101.4 104.3 78.8 69.3 145.4 103.1 117.4 47.1;z=29.7 26.9 28.1 24.9 24.6 31.2 38.6 20.8 33.7 25.9 20.6 22.3 21.3 19.2;xmin=min(x)xmax=max(x)ymin=min(y)ymax=max(y)x0=linspace(xmin,xmax,20)y0=linspace(ymin,ymax,20)x0,y0=meshgrid(x0,y0)z0=griddata(x,y,z,x0,y0,v4)surf(x0,y0,z0)a

41、xis(0,190,0,180,15,45)colormap(gray)xlabel(x)ylabel(y)zlabel(z) %畫出二氧化硫的平面分布圖x1=132.1 126.5 135.3 130.1 128.3 148.5 130.7 104.7 149.7 94.9 101.4 83.1 76.9 66.5;x2=91.9 76.8 73.7 66.5 61.2 62.3 101.4 104.3 78.8 69.3 145.4 103.1 117.4 47.1;x=x1 x2y=64 63.6 57.5 63.3 51.8 61.9 55.8 44.8 65.7 43.5 35 37

42、.9 43.6 34.4;myfunc=(beta,x)beta(1)-log(beta(2)*beta(3)-(x(:,1)-beta(4)/(4*(beta(2)2)+(x(:,2)-beta(5)/(4*(beta(3)2)+beta(6)beta=nlinfit(x,y,myfunc,1000,0.103,0.103,149700,78800,0.01)%函數(shù)擬合確定二氧化硫污染源位置x=132.1 126.5 135.3 130.1 128.3 148.5 130.7 104.7 149.7 94.9 101.4 83.1 76.9 66.5;y=91.9 76.8 73.7 66.

43、5 61.2 62.3 101.4 104.3 78.8 69.3 145.4 103.1 117.4 47.1;z=112.4 115 118.5 116.4 114.8 131.7 123.5 98.5 130.3 94.7 98.7 95.3 85.9 78.3;xmin=min(x)xmax=max(x)ymin=min(y)ymax=max(y)x0=linspace(xmin,xmax,20)y0=linspace(ymin,ymax,20)x0,y0=meshgrid(x0,y0)z0=griddata(x,y,z,x0,y0,v4)surf(x0,y0,z0)axis(0,19

44、0,0,180,70,135)colormap(gray)xlabel(x)ylabel(y)zlabel(z) %PM10的空間分布圖x1=132.1 126.5 135.3 130.1 128.3 148.5 130.7 104.7 149.7 94.9 101.4 83.1 76.9 66.5;x2=91.9 76.8 73.7 66.5 61.2 62.3 101.4 104.3 78.8 69.3 145.4 103.1 117.4 47.1;x=x1 x2y=112.4 115 118.5 116.4 114.8 131.7 123.5 98.5 130.3 94.7 98.7 9

45、5.3 85.9 78.3;myfunc=(beta,x)beta(1)-log(beta(2)*beta(3)-(x(:,1)-beta(4)/(4*(beta(2)2)+(x(:,2)-beta(5)/(4*(beta(3)2)+beta(6)beta=nlinfit(x,y,myfunc,1000,0.103,0.103,148000,62000,0.01) %確定PM10的污染源位置x=132.1 126.5 135.3 130.1 128.3 148.5 130.7 104.7 149.7 94.9 101.4 83.1 76.9 66.5;y=91.9 76.8 73.7 66.5

46、 61.2 62.3 101.4 104.3 78.8 69.3 145.4 103.1 117.4 47.1;z=91.3 98.1 96.9 92.8 91.1 106.8 105.7 84.8 107.8 79.2 84.8 76.1 71.6 68.0;xmin=min(x)xmax=max(x)ymin=min(y)ymax=max(y)x0=linspace(xmin,xmax,20)y0=linspace(ymin,ymax,20)x0,y0=meshgrid(x0,y0)z0=griddata(x,y,z,x0,y0,v4)surf(x0,y0,z0)axis(0,190,0,

47、180,65,115)colormap(gray)xlabel(x)ylabel(y)zlabel(z) %畫出PM2.5的平面分布圖x1=132.1 126.5 135.3 130.1 128.3 148.5 130.7 104.7 149.7 94.9 101.4 83.1 76.9 66.5;x2=91.9 76.8 73.7 66.5 61.2 62.3 101.4 104.3 78.8 69.3 145.4 103.1 117.4 47.1;x=x1 x2y=91.3 98.1 96.9 92.8 91.1 106.8 105.7 84.8 107.8 79.2 84.8 76.1

48、71.6 68.0;myfunc=(beta,x)beta(1)-log(beta(2)*beta(3)-(x(:,1)-beta(4)/(4*(beta(2)2)+(x(:,2)-beta(5)/(4*(beta(3)2)+beta(6)beta=nlinfit(x,y,myfunc,1000,0.103,0.103,143000,79000,0.01) %確定PM2.5污染源位置3.時(shí)間序列模型a=data(152:273,1);t=1:1:122;r11=autocorr(a) %畫出數(shù)據(jù)的自相關(guān)圖r12=parcorr(a) %畫出數(shù)據(jù)的偏自相關(guān)圖da=diff(a);r21=auto

49、corr(da) %畫出一階差分?jǐn)?shù)據(jù)的自相關(guān)圖r22=parcorr(da) %畫出一階差分?jǐn)?shù)據(jù)的偏自相關(guān)圖n=length(da); for i=0:3 for j=0:3 spec=garchset(R,i,M,j,Display,off);coeffX,errorsX,LLFX=garchfit(spec,da);num=garchcount(coeffX);aic,bic=aicbic(LLFX,num,n);fprintf(R=%d,M=%d,AIC=%f,i,j,aic,bic)end end %利用AIC準(zhǔn)則確定模型的階數(shù)spec2=garchset(R,1,M,1,Display,off);coeffX,errorsX,LLFX=garchfit(spec2,da)sigmaFor