數(shù)值分析實(shí)驗(yàn)報(bào)告-12

成都空氣質(zhì)量指數(shù)擬合與預(yù)測(cè)——基于正交最小二乘法擬合Xxxxx一、問題敘述空氣質(zhì)量的好壞直接反映了空氣的污染程度，它是依據(jù)空氣中污染物濃度的高低來判斷的，空氣中的污染物主要有大氣中直徑小于或等于2.5微米的顆粒物（PM2.5）、可吸入懸浮顆粒物\浮塵（PM10）、二氧化硫（SO2）、二氧化氮（NO2）、一氧化碳（CO）等。目前，城市空氣質(zhì)量污染指數(shù)的分級(jí)標(biāo)準(zhǔn)是根據(jù)空氣質(zhì)量指數(shù)（AQI）的取值界定的，并代表著空氣污染程度和空氣質(zhì)量狀況。成都是中國(guó)重要的西南部城市，其工業(yè)化與城市化的快速發(fā)展對(duì)城市環(huán)境產(chǎn)生了重要影響。對(duì)空氣質(zhì)量進(jìn)行預(yù)測(cè)并采取合適的控制措施就有著至關(guān)重要的作用，因此，本文根據(jù)搜集的成都2013年11月-2014年11月的空氣質(zhì)量的數(shù)據(jù)資料，基于正交最小二乘法建立適當(dāng)?shù)臄?shù)學(xué)模型，選出影響空氣質(zhì)量的主要污染物并對(duì)成都市空氣質(zhì)量變化以及主要污染物的來源進(jìn)行分析，并根據(jù)擬合模型預(yù)測(cè)成都短期的空氣質(zhì)量指數(shù)。為了便于后續(xù)分析對(duì)AQI的解釋如下：AQI：空氣質(zhì)量指數(shù)（AirQualityIndex）是定量描述空氣質(zhì)量狀況的指數(shù)，其數(shù)值越大說明空氣污染狀況越嚴(yán)重，對(duì)人體健康的危害也就越大。二、問題分析由于matlab沒有提供現(xiàn)成的函數(shù)用于實(shí)現(xiàn)正交最小二乘擬合，所以需要用matlab自定義編寫正交最小二乘擬合程序，然后再用matlab畫出污染物AQI、PM2.5、PM10、CO、NO2、SO2的離散數(shù)據(jù)圖以及AQI、PM2.5、PM10、CO、NO2、SO2的擬合曲線圖，然后通過這些圖分別對(duì)各污染物隨時(shí)間的變化特征進(jìn)行分析總結(jié)，最后根據(jù)擬合的模型帶入數(shù)據(jù)預(yù)測(cè)未來的AQI、PM2.5、PM10、CO、NO2、SO2指數(shù)值。由于影響空氣質(zhì)量的因素較多以及有很多的不確定性，因此我們對(duì)該擬合模型做以下假設(shè)：假設(shè)搜集數(shù)據(jù)均有效，忽略空氣中在指標(biāo)控制范圍內(nèi)的其他污染物，空氣質(zhì)量相同等級(jí)的污染程度相同，各種影響空氣質(zhì)量的因素保持平穩(wěn)變化并且各種因素對(duì)環(huán)境的影響最終主要表現(xiàn)在PM2.5、PM10、CO、NO2、SO2上，不考慮其他隨機(jī)因素的影響。運(yùn)用matlab對(duì)成都過去11個(gè)月以及月內(nèi)的主要污染物擬合分析發(fā)現(xiàn)在多項(xiàng)式擬合中四和五次曲線模型的較接近實(shí)際情況，其曲線模型擬合度較高，因此我們分別運(yùn)用四次和五次曲線模型進(jìn)行擬合及預(yù)測(cè)分析。三、實(shí)驗(yàn)程序程序一：基于正交最小二乘擬合的matlab程序%******************************************%%基于正交最小二乘擬合的matlab程序%x輸入數(shù)據(jù)向量，y實(shí)際輸出數(shù)據(jù)向量，m擬合階次%******************************************%functionf=zhengjiaomin(x,y,m)%基于正交最小二乘擬合的matlab程序if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等');return;end%檢查維數(shù)symsv;d=zeros(1,m+1);q=zeros(1,m+1);alpha=zeros(1,m+1);fork=0:mpx(k+1)=power(v,k);end%x的冪多項(xiàng)式H2=[1];d(1)=n;forh=1:nq(1)=q(1)+y(h);alpha(1)=alpha(1)+x(h);endq(1)=q(1)/d(1);alpha(1)=alpha(1)/d(1);f(1)=q(1);H1=[-alpha(1)1];forh=1:nd(2)=d(2)+(x(h)-alpha(1))^2;q(2)=q(2)+y(h)*(x(h)-alpha(1));alpha(2)=alpha(2)+x(h)*(x(h)-alpha(1))^2;endq(2)=q(2)/d(2);alpha(2)=alpha(2)/d(2);f(1)=f(1)+q(2)*(-alpha(1));f(2)=q(2);beta=d(2)/d(1);fori=3:(m+1)H=zeros(1,i);H(i)=H1(i-1);H(i-1)=-alpha(i-1)*H1(i-1)+H1(i-2);forj=2:i-2;H(j)=-alpha(i-1)*H1(j)+H1(j-1)-beta*H2(j);endH(1)=-alpha(i-1)*H1(1)-beta*H2(1);HF=H*transpose(px(1:i));forh=1:nW=subs(HF,'v',x(h));d(i)=d(i)+(W)^2;q(i)=q(i)+y(h)*W;alpha(i)=alpha(i)+x(h)*(W)^2;endalpha(i)=alpha(i)/d(i);q(i)=q(i)/d(i);beta=d(i)/d(i-1);fork=1:i-1f(k)=f(k)+q(i)*H(k);endf(i)=q(i)*H(i);H2=H1;H1=H;end程序二：2013.11月和2014.11月內(nèi)各日AQI、PM2.5、PM10數(shù)據(jù)擬合%2013.11和2014.11月AQI、PM2.5、PM10數(shù)據(jù)擬合functionkongqidate%定義m函數(shù)xone=[1:1:29];%某月的日期%2014.11月與2013年11月各日的AQI、PM2.5、PM10數(shù)據(jù)從excel成都空氣%質(zhì)量數(shù)據(jù)中的14.11和13.11單元表讀取，數(shù)據(jù)見附件“成都空氣質(zhì)量數(shù)據(jù)”aqione=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','14.11','B2:B30');%2014.11AQI數(shù)據(jù)m25one=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','14.11','E2:E30');m10one=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','14.11','F2:F30');aqitwo=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','13.11','B2:B30');%2013.11AQI數(shù)據(jù)m25two=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','13.11','E2:E30');m10two=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','13.11','F2:F30');none=5;%某月各天數(shù)據(jù)的擬合階次%2014.11的正交最小二乘擬合多項(xiàng)式系數(shù)計(jì)算，調(diào)用自編matlabm函數(shù)faqione=zhengjiaomin(xone,aqione,none);%2014.11的AQI數(shù)據(jù)擬合fm25one=zhengjiaomin(xone,m25one,none);fm10one=zhengjiaomin(xone,m10one,none);faqitwo=zhengjiaomin(xone,aqitwo,none);fm25two=zhengjiaomin(xone,m25two,none);fm10two=zhengjiaomin(xone,m10two,none);%系數(shù)按降冪排列并輸出faqione=fliplr(faqione)fm25one=fliplr(fm25one)fm10one=fliplr(fm10one)faqitwo=fliplr(faqitwo)fm25two=fliplr(fm25two)fm10two=fliplr(fm10two)%計(jì)算各日擬合指標(biāo)值yaqione=polyval(faqione,xone);ym25one=polyval(fm25one,xone);ym10one=polyval(fm10one,xone);yaqitwo=polyval(faqitwo,xone);ym25two=polyval(fm25two,xone);ym10two=polyval(fm10two,xone);%根據(jù)數(shù)據(jù)擬合的模型預(yù)測(cè)2014年和2013年11月30日的各項(xiàng)指標(biāo)，并輸出yaqi1one=polyval(faqione,30)ym251one=polyval(fm25one,30)ym101one=polyval(fm10one,30)yaqi1two=polyval(faqitwo,30)ym251two=polyval(fm25two,30)ym101two=polyval(fm10two,30)holdon%用不同符號(hào)畫出各指標(biāo)的離散數(shù)據(jù)圖plot(xone,aqione,'*',xone,m25one,'o',xone,m10one,'v',xone,aqitwo,'+',xone,m25two,'h',xone,m10two,'s')holdon%不同顏色畫出各指標(biāo)的擬合曲線%實(shí)線代表2014.11AQI藍(lán)色，PM2.5綠色，PM10紅色%虛線代表2013.11AQI藍(lán)綠色，PM2.5紫紅色，PM10黃色plot(xone,yaqione,'b',xone,ym25one,'g',xone,ym10one,'r',xone,yaqitwo,'c--',xone,ym25two,'m--',xone,ym10two,'y--')holdonlegend('2014AQI指數(shù)','2014PM2.5','2014PM10','2013AQI指數(shù)','2013PM2.5','2013PM10');%圖例標(biāo)識(shí)xlabel('xnoe時(shí)間/天');ylabel('y各日指數(shù)值');%標(biāo)定坐標(biāo)軸及名稱程序三：2013.11-2014.10月的AQI、PM2.5、PM10、CO、NO2、SO2數(shù)據(jù)擬合%2013.11-2014.10AQIPM2.5PM10CONO2SO2數(shù)據(jù)擬合functionkongqimonth%定義m函數(shù)x=[1:1:12];%2013.11-2014.10各月，按序分別代表2013.11，2013.12，2014.01......2014.10%2013.11-2014.10AQIPM2.5PM10CONO2SO2數(shù)據(jù)從成都空氣質(zhì)量數(shù)據(jù)excel%中month單元表讀取，數(shù)據(jù)見附件“成都空氣質(zhì)量數(shù)據(jù)”aqi=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','month','B2:B13');%2013.11-2014.10AQI數(shù)據(jù)從excelmonth單元表讀取m25=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','month','D2:D13');m10=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','month','E2:E13');co=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','month','F2:F13');no2=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','month','G2:G13');so2=xlsread('成都空氣質(zhì)量數(shù)據(jù).xlsx','month','H2:H13');n=4;%各月數(shù)據(jù)的擬合階次%2013.11-2014.10的正交最小二乘擬合多項(xiàng)式系數(shù)計(jì)算，調(diào)用自編matlab函數(shù)faqi=zhengjiaomin(x,aqi,n);%2013.11-2014.10的AQI數(shù)據(jù)擬合fm25=zhengjiaomin(x,m25,n);fm10=zhengjiaomin(x,m10,n);fco=zhengjiaomin(x,co,n);fno2=zhengjiaomin(x,no2,n);fso2=zhengjiaomin(x,so2,n);%系數(shù)按降冪排列并輸出faqi=fliplr(faqi)fm25=fliplr(fm25)fm10=fliplr(fm10)fco=fliplr(fco)fno2=fliplr(fno2)fso2=fliplr(fso2)%計(jì)算各月擬合指標(biāo)值yaqi=polyval(faqi,x);ym25=polyval(fm25,x);ym10=polyval(fm10,x);yco=polyval(fco,x);yno2=polyval(fno2,x);yso2=polyval(fso2,x);%根據(jù)數(shù)據(jù)擬合的模型預(yù)測(cè)2014年11月的各項(xiàng)指標(biāo)并輸出yaqi1=polyval(faqi,13)%x=13代表2014年11月AQI預(yù)測(cè)值ym251=polyval(fm25,13)ym101=polyval(fm10,13)yco1=polyval(fco,13)yno21=polyval(fno2,13)yso21=polyval(fso2,13)%用不同符號(hào)畫出各指標(biāo)的離散數(shù)據(jù)圖plot(x,aqi,'*',x,m25,'o',x,m10,'v',x,co,'+',x,no2,'h',x,so2,'s')holdon%不同顏色畫出各指標(biāo)的擬合曲線%AQI藍(lán)色，PM2.5綠色，PM0紅色，CO藍(lán)綠色，NO2紫紅色，SO2黃色plot(x,yaqi,'b',x,ym25,'g',x,ym10,'r',x,yco,'c',x,yno2,'m',x,yso2,'y')xlabel('x時(shí)間/月');ylabel('y各月指數(shù)值');%標(biāo)定坐標(biāo)軸及名稱legend('AQI指數(shù)','PM2.5','PM10','Co','No2','So2');%圖例標(biāo)識(shí)holdon程序四：在matlabCommandwindow分別調(diào)用kongqidate（程序二）、kongqimonth（程序三）就可得到擬合結(jié)果數(shù)據(jù)及圖形。四、實(shí)驗(yàn)數(shù)據(jù)及結(jié)果分析1、對(duì)2013年11月至2014年10月的各項(xiàng)指標(biāo)采用四次曲線的擬合結(jié)果為：yaqi=-0.1732x4+5.3046x3-54.2104x2+196.1882x-44.4040ypm25=-0.1489x4+4.5790x3-46.8634x2+169.1943x-51.9116ypm10=-0.1611x4+5.0701x3-52.7769x2+190.1411x-9.8662yco=-0.0011x4+0.0327x3-0.3250x2+1.1131x+0.3962yno2=-0.0193x4+0.6098x3-6.1883x2+19.9566x+48.5707yso2=-0.0114x4+0.3579x3-3.5240x2+9.6499x+23.7071以上各式分別代表2013年11月到2014年10月的AQI、PM2.5、PM10、CO、NO2、SO2的擬合模型，擬合曲線如圖1所示。圖1-2013.11-2014.10的各指標(biāo)擬合曲線圖通過擬合的模型預(yù)測(cè)2014年11月的各項(xiàng)指標(biāo)數(shù)據(jù)如下：表1-2014年11月各項(xiàng)指標(biāo)預(yù)測(cè)值2014.11AQIPM2.5PM10CONO2SO2預(yù)測(cè)數(shù)據(jù)51.313134.138980.63890.957950.520215.2929實(shí)際數(shù)據(jù)5741891.064815分析表1可以發(fā)現(xiàn)2014年11月的AQI、PM2.5、PM10、CO、NO2、SO2預(yù)測(cè)數(shù)據(jù)與實(shí)際數(shù)據(jù)保持較小的誤差，處于可接受范圍。當(dāng)繼續(xù)通過上述模型預(yù)測(cè)發(fā)現(xiàn)，該模型預(yù)測(cè)值隨著時(shí)間的逐漸增加將降低至負(fù)值，而根據(jù)實(shí)際情況，各項(xiàng)指標(biāo)值不可能低至負(fù)值，因此這個(gè)模型只適合短時(shí)間預(yù)測(cè)，對(duì)于長(zhǎng)期的預(yù)測(cè)會(huì)出現(xiàn)較大的誤差。為了更好的評(píng)價(jià)成都市空氣質(zhì)量，查得空氣污染指數(shù)范圍及相應(yīng)的空氣質(zhì)量類別如表2。表2-空氣污染指數(shù)范圍及相應(yīng)的空氣質(zhì)量類別空氣污質(zhì)量指數(shù)AQI空氣質(zhì)量狀況對(duì)健康的影響建議采取的措施0～50優(yōu)空氣質(zhì)量令人滿意，基本無空氣污染各類人群可正?；顒?dòng)51～100良101～150輕度污染易感人群癥狀有輕加劇，健康人群出現(xiàn)刺激癥狀心臟病和呼吸系統(tǒng)疾病患者應(yīng)減少體力消耗和戶外活動(dòng)151～200中度污染201～300重度污染心臟病和肺病患者癥狀顯著加劇，運(yùn)動(dòng)耐受力降低，健康人群中普遍出現(xiàn)癥狀老年人和心臟病、肺病患者應(yīng)在停留在室內(nèi)，并減少體力活動(dòng)>300嚴(yán)重污染健康人運(yùn)動(dòng)耐受力降低，有明顯強(qiáng)烈癥狀，提前出現(xiàn)某些疾病老年人和病人應(yīng)當(dāng)留在室內(nèi)，避免體力消耗，一般人群應(yīng)避免戶外活動(dòng)從圖1中可分析得出，從2013年11月到2014年10月，成都的AQI指數(shù)徘徊于50～200之間，前6個(gè)月（2013年11月-2014年4月）AQI指數(shù)徘徊在100～200之間，其中2013年12月到2014年2月污染較重，屬于中度污染，而2014年1月污染最重，并且其主要污染物為PM2.5和PM10，后6個(gè)月（2014年5月-2014年10月）AQI指數(shù)徘徊在50～100之間，結(jié)合表2發(fā)現(xiàn)，成都空氣質(zhì)量由輕度和中度污染狀態(tài)逐漸轉(zhuǎn)變?yōu)榱?，說明成都的環(huán)境治理工作得到了一定成果，空氣質(zhì)量正在朝好的方向不停發(fā)展。PM2.5：主要來源是汽車尾氣的排放、集中供暖產(chǎn)生的廢氣。PM10主要來源是汽車尾氣的排放、不合格煙塵排放。由圖1不難看出PM2.5和PM10從2013年11月到2014年1月上升的，2014年1月達(dá)到這幾個(gè)月的最大值，然后一直到2014年7月都在下降，7月的時(shí)候PM2.5和PM10均達(dá)到最小值，但是到8月份的時(shí)候又開始有上升的趨勢(shì)，但是整體逐漸成下降趨勢(shì)。NO2：主要來源是汽車尾氣的排放。分析NO2的趨勢(shì)線可知，2013年11月份至2014年10月份，相對(duì)于趨于平緩趨勢(shì)。但在2013年11月份至2014年1月份略有增加，可能春節(jié)用車量的增加有關(guān)。SO2：主要來源是集中供暖產(chǎn)生的廢氣。分析SO2的趨勢(shì)線可知，2013年11月份至2014年1月份，SO2相對(duì)處于高位，這與冬季供暖期污染源增加相對(duì)應(yīng)。預(yù)計(jì)，每年的SO2污染濃度最大值與當(dāng)年的最冷月相對(duì)應(yīng)。2、由上述擬合發(fā)現(xiàn)PM2.5、PM10對(duì)空氣質(zhì)量影響最大，因此下面對(duì)于月份內(nèi)各天的數(shù)據(jù)擬合僅對(duì)AQI、PM2.5、PM10這三項(xiàng)指標(biāo)進(jìn)行分析，對(duì)2014年11月和2013年11月的數(shù)據(jù)采用五次曲線模型擬合的AQI、PM2.5、PM10指標(biāo)結(jié)果為：yaqi14=0.0003x5-0.0234x4+0.7178x3-9.6652x2+54.0198x-14.3235ypm2514=0.0002x5-0.0203x4+0.6254x3-8.4612x2+47.6679x-29.0132ypm1014=0.0004x5-0.0338x4+1.0335x3-13.8206x2+76.9661x-37.2057yaqi13=0.0007x5-0.0500x4+1.3074x3-14.1673x2+59.8361x-3.6069ypm2513=0.0007x5-0.0481x4+1.2484x3-13.4747x2+56.8556x-19.9148ypm1013=0.0007x5-0.0539x4+1.3682x3-14.3064x2+58.4476x+6.1893以上前三式分別代表2014年11月的AQI、PM2.5、PM10的擬合模型，后三式分別代表2013年11月的AQI、PM2.5、PM10的擬合模型，擬合曲線如圖2所示。圖2-2013.11和2014.11月的各指標(biāo)擬合曲線圖通過擬合的模型預(yù)測(cè)2014年和2013年的11月30日的AQI、PM2.5、PM10指標(biāo)數(shù)據(jù)如下：表3-2014年11月30日與2013年11月30日的主要指標(biāo)預(yù)測(cè)預(yù)測(cè)數(shù)據(jù)AQIPM2.5PM102014.11.3081.187754.4425104.32212013.11.30167.2111151.1519278.2431表4--2014年11月30日與2013年11月30日的主要指標(biāo)實(shí)際值實(shí)際數(shù)據(jù)AQIPM2.5PM102014.11.307349912013.11.30154145278分析表3與表4可以發(fā)現(xiàn)2014年和2013年11月30日的AQI、PM2.5、PM10預(yù)測(cè)數(shù)據(jù)與實(shí)際數(shù)據(jù)保持較小的誤差，處在可接受范圍。分析圖2可以發(fā)現(xiàn)2014年11月空氣質(zhì)量變化平穩(wěn)且基本處于良，而2013年同期相比之下污染稍重，這說明2014年空氣質(zhì)量得到了較好的提升。 通過查閱2013年11月和2014年11月降雨數(shù)據(jù)（見附錄）發(fā)現(xiàn)，2013年11月后半月降雨天數(shù)相對(duì)前半月較少因而AQI指數(shù)較前半月增加，污染加重，而在后半月21日至26日又有降雨使得空氣質(zhì)量指數(shù)下降，空氣質(zhì)量逐漸趨于良好。2014年11月整體降雨天數(shù)較2013年同期多，空氣質(zhì)量指數(shù)也較2013年低，空氣質(zhì)量整體趨于良好?？梢娍諝赓|(zhì)量指數(shù)與降雨有很大關(guān)系，降雨對(duì)空氣污染有很大緩解。五、實(shí)驗(yàn)結(jié)論在數(shù)據(jù)不太少的情況下，該模型有一定的誤差，但短期預(yù)測(cè)可以得到誤差不太大的預(yù)測(cè)值，而長(zhǎng)期預(yù)測(cè)則會(huì)出現(xiàn)負(fù)值等不合理數(shù)據(jù)，因此，不適合做長(zhǎng)期的預(yù)測(cè)或分析。另外，用正交最小二乘法擬合來求解問題時(shí)要選取合適的階次，當(dāng)階次過大或過小時(shí)均會(huì)出現(xiàn)不符合實(shí)際情況的數(shù)據(jù)。六、注記本模型還有比較大的改進(jìn)的空間，比如改成雙指數(shù)模型，誤差應(yīng)該會(huì)減少。還可以考慮利用更多的數(shù)據(jù)來進(jìn)行預(yù)測(cè)以及對(duì)數(shù)據(jù)進(jìn)行篩選，去除奇異數(shù)據(jù)，再進(jìn)行擬合，這樣得到的擬合及預(yù)測(cè)結(jié)果會(huì)更準(zhǔn)確，還可對(duì)各污染物的變化規(guī)律對(duì)AQI影響權(quán)值做分析。七、附錄1、