城空气污染程度的分析报告和预测模型

1、实用标准文档城市空气污染程度的分析和预测摘 要 本文讨论了有关城市污染程度、污染因素及污染扩散的问题.对于问题一, 本文主要从大气污染、 噪声污染和水体污染这三个面选取主要 污染物,查阅北京、天津、上海、重庆和西安五座城市 2007-2021 年的年度平均 污染数据,采用降维的思想,运用主成分分析法减少变量个数,再借助 Matlab 软件计算各主成分的奉献率, 分析知可选取前三个主成分作为衡量污染程度的标 准,最后根据综合指标得到这五个城市的污染程度从高到低依次为： 重庆、上海、 北京、天津、西安.通过判断相关系数的大小,确定五个城市影响人们生活的主 要污染因素是水污染,其四项指标依次为化学需

2、氧量、总氮、总磷和氨氮.对于问题二,以北京市大气污染为例.首先,利用GPS记录北京市14个城区观测点的位置,并查阅2021年污染指标SO2、CO2、PM2.5与PM10的污染数 据,绘制出相应的空间浓度分布图, 估计这四种污染物的大致污染源位置依次为： (110,100)附近、 (130,83)附近, (125,85)附近和 (132,80)附近；其次,根据污染扩 散原理和方式,建立Cauchy污染传播模型,根据各地区空气污染物的浓度分布, 运用 Matlab 软件对数据非线性拟合,得出扩散模型各参数的值,计算得出各项 污染指标的污染源位置依次为： (115.5,97.3), (128.2,8

3、5.3), (129.1,80.8) 和 (125.5,87.6) ；最后,比拟污染物位置的计算值与实际值,发现误差相差较小, 故模型建立较为合理.对于问题三,分析西安市的主要污染大气污染.收集西安市 2021 年 4.1-7.31 日的空气污染数据,根据时间序列的平稳性特点及 AIC 定阶准那么选取 适宜的时间序列模型 ARMA (1,1) ,利用 Matlab 软件对序列模型的各项参数进行 估计并检验模型的合理性, 并将模型用于数据预报. 利用时间序列模型预测西安 市未来 10天的空气污染状况总体等级为良.对于问题四,基于问题一、二、三对污染因素的分析和污染扩散的特点,主 要从减少污染物的

4、产生和治理净化已产生的污染物两方面, 针对大气污染、 水体 污染和噪声污染为相关部门提供合理化防治建议.关键词 主成分分析；Cauchy污染传播模型；时间序列模型；Matlab软件文案大全一、问题重述环境质量的好坏对人们的生活有着最直接的影响, 包括身体健康、寿命长短、 生活质量以及经济开展水平、 生态环境等.因此保护环境的工作需要个人、 集体、 国家、乃至全球各国的共同努力.事实上,环境保护事业从 1962 年起步至今, 一直备受人们关注和重视, 环境污染的防治工作也一步步开展运营, 然而较经济 开展的速度和水平, 还有很多需要为之努力的工作要做. 环境保护工作一方面是 要从减少污染物的产生

5、方面, 预防和限制污染源头污染物的排放, 另一方面是治 理净化已产生的污染物.环境污染是指环境因受人类活动影响而改变了原有性质和状态的现象, 其实 质是人类活动中将大量的污染物排入环境, 影响其自净水平, 降低了生态系统的 功能.为考察某时期某区域的环境状况, 需统计该地区各项污染指标的数据, 根据 各污染指标的污染排放标准, 确定其污染程度, 并通过污染发生的原因、 途径和 后果制定合理的污染治理方案,加以实施运行,改善环境.统计北京、天津、上海、重庆和西安五座城市的污染数据,讨论如下问题： 问题一：比拟五座城市的污染程度及影响人们生活的主要污染因素； 问题二：建立城市污染扩散模型确定某城市

6、污染程度并检验； 问题三：给出西安市某方面污染状况的中长期预测数据； 问题四：给相关部门提供一篇治理污染的建议短文.二、问题分析 对于问题一,本文将查阅收集北京、上海、天津、重庆以及西安五座城市大 气污染、噪声污染和水体污染 2007-2021 年的数据,由于各项污染指标众多, 所 以选取具有代表性的污染物：可吸入颗粒物 PM10、SO2、CO2 、道路交通噪声、 区域环境噪声、 化学需氧量、 氨氮、总氮和总磷这 9 项指标来综合衡量各城市的 污染程度, 本文将运用降维的思想, 利用主成分分析法计算前几个主成分的奉献 率大小及通过线性组后的综合指标, 以此来衡量各城市的污染程度, 根据各主成

7、分的相关系数大小即可确定影响人们生活的主要污染因素 (1) .对于问题二,本文将以北京市的大气污染为例,收集北京市 2021 年 14 个城 区的SO2、CO2、PM2.5和PM10的污染数据,绘制出这四种污染物浓度的空间 浓度分布图,估算出这四种污染物的大致污染源位置；建立Cauchy污染传播模型,根据各地区空气污染物的浓度分布,运用 Matlab 软件对数据 nlinfit 非线 性拟合,计算出扩散模型各参数值, 得出各空气污染物的污染源位置； 比拟污染 物位置的计算值与实际值,检验污染传播模型建立是否合理 (2) .对于问题三,本文将考虑西安市大气污染的各项污染指标,搜集 2021 年

8、4 月 1 日至 2021年 7 月 31 日的污染数据,根据时间序列的平稳性特点及 AIC 定阶 准那么选取适宜的时间序列模型, 利用 Matlab 软件(3) 估计序列模型的各项参数 并检验模型的合理性,考核通过时方可用于数据预报.对于问题四,基于问题一、二、三对污染因素的分析和污染扩散的特点,主 要从减少污染物的产生和治理净化已产生的污染物两方面, 针对大气污染、 水体 污染和噪声污染为相关部门提供合理化防治建议.三、模型假设1. 所有城市污染数据均真实有效,具有统计价值；2. 污染数据的取样为瞬时取样且在某个平面内, 不考虑取样前后污染物的浓度变 化；3. 绘制北京大气污染浓度分布图时

9、不考虑海拔因素的影响；4. 大气污染过程经历着扩散与衰减,且沿 x,y轴方向上的扩散系数分别为常数, 衰减与浓度成正比；5. 不考虑日照等自然因素对各项指标污染浓度的影响.四、符号说明4.1模型一XiX2X3X4X5X6X7X8X9Z可吸入颗粒物单位mg/m3; 二氧化硫单位mg/m3； 二氧化氮单位mg/m ;道路交通噪声单位分贝； 区域环境噪声单位分贝； 化学需氧量单位万吨； 氨氮单位万吨； 总氮单位万吨；总磷单位万吨；ei样本协方差矩阵； 协方差矩阵的特征值； 协方差矩阵的特征向量； 相关系数； 样本方差向量； 样本相关矩阵； 样本相关矩阵的特征值； 样本相关矩阵的特征向量；7:样本均值

10、向量；丫:主成分； 丫 :污染程度综合指标.其中 i =12,94.2模型二x :位置坐标/横轴； y :位置坐标/纵轴； z :位置坐标/竖轴； c :污染物浓度； m :污染物质量；M :污染物排放量;CT :扩散系数.4.3模型三XtSt序列；平稳白噪声；MA (q)滑动平均序列；AR(p)自回归序列；ARMA( p, q):自回归滑动平均序列c常数；自协方差函数；Pk自相关函数；f(x)概率密度；片未知参数向量?(m) k L最大似然估计；自回归参数.t =0,_1,_2,;i =1,2, ,k;l =1,2, ,P五、模型建立与求解5.1问题一的模型建立与求解环境因受人类活动影响而改

11、变了原有性质和状态的现象称为环境污染. 环境 污染的实质是人类活动中将大量的污染物排入环境, 影响其自净水平,降低了生 态系统的功能.然而生活中环境的污染也是多方面的,污染按环境要素分类,可 分为：大气污染、土壤污染和水体污染.为此,本文选取环境污染的主要污染影响指标, 建立污染评价模型,即主成 分分析模型,比拟北京、天津、上海、重庆和西安五座城市的污染程度,并得出 影响人们生活的主要污染因素.5.1.1建立主成分分析模型首先选取环境污染的局部污染指标：可吸入颗粒物x1、二氧化硫x2、二氧化氮X3、道路交通噪声X4、区域环境噪声X5、化学需氧量X6、氨氮X7、总氮X8 和总磷X9.北京、天津、

12、上海、重庆和西安五座城市的各项指标年平均数据(见 附录表1).问题中变量太多不但会增加计算的复杂性, 也给合理的分析问题和解决问题 带来很大困难；同时,这些变量之间存在的相关性,也使得这些变量所反映的信 息在一定程度上有所重叠.为了减少变量个数,本文采用了降维的思想,利用主 成分分析法得出污染空气的主要因素变量,提升问题研究的合理性,同时不会使 数据反映的信息量有大的损失.设V是X =(X1,X2,X3,Xp)T的协方差矩阵,V的特征值与正交化特征向量 分别为 - q - I -p - 0及e1,e2,e3/ ®,且x的第i个主成分为丫弋必 X2 $3X3 飞pXp,(i =1,2,

13、3, , p)( 1)实用标准文档1 n根据平均值计算公式Xj二丄Xj代入数据计算样本X = X!,X2,X3,Xp的均值向 n jj量 X =Xi, X2,X3,Xp为X= 0.0986 0.0354 0.0434 68.3480 54.4760 23.5925 3.2110 2.9390 0.3475根据协方差矩阵计算公式1 n* - XXj - XT 2n - 1 i 4利用Matlab软件代入数据可求得随机变量x = %,X2, X3,Xp p =9相应样本协方差矩阵为只写下三角局部0.00030.00010.0001000-0.0034-0.00640.00500.78200.002

14、00.000300.08910.2917-0.1114-0.0031-0.0409-5.0768-3.8197110.5272-0.0281-0.0062-0.0051-0.3351-0.388116.06303.00581-0.00280.0039-0.0049-0.9387-0.789514.80031 .4350 2.8425-0.00030.0003-0.0005-0.0963-0.10291.75740.1690 0.34850.0437由于协方差矩阵对角线元素即为样本X = (%,X2, X3,Xp的方差值,观察矩11 =阵易知其方差值相差较大,即各量纲间分散程度较大,不适用于环境

15、污染模型中 各成分指标间的相关性分析,为消除由于量纲不同可能带来的结果误差,故将样本变量X=X1,X2,X3, ,Xp p =9标准化,即令X 4.X- i i ,i =1,2, ,93g得到一组新的数据X =X,X2,X9.其中Ji为Xi的平均值,二ii为Xi的方 差.此时,X = X1 , X2,X9 的协方差矩阵即为样本变量X=X1,X2X9的 相关矩阵'=ijp p,其中?j询Xi ,Yj )=Cov(Xi,Xj)(4)文案大全运用Matlab软件代入数据即可得其样本相关矩阵为CovXj,Xj为 Xi,Xj 的协方差. 那么根据公式4-10.602610.1171-0.5158

16、1-0.2097-0.74460.910210.19650.04850.00130.186611-0.5744-0.0229-0.6245-0.5461-0.67271-0.8787-0.3709-0.4726-0.2186-0.41440.881311-0.09080.2378- 0.4691-0.6296-0.86700.83500.49091.0.08870.1325-0.3486-0.5213-0.91140.79980.46640.9891P =由特征值计算公式止冃=0(5)第i个主成分的奉献率丸i / P(6)及前k个主成分的累计奉献率送入/Ps4/(7)其中p = 9.计算的特征

17、值及各标准化主成分的奉献率和累计奉献率如表1所示：表1 r的特征值、奉献率和累计奉献率i奉献率(%)累计奉献率(%)14.742052.6952.6922.588428.7681.4531.476816.4197.8640.19282.14100500100600100700100800100900100由表1易知,第一、第二和第三主成分的累积奉献率为97.86%,并且只有第一、第二和和第三主成分对应的特征值大于 1,所以选取第一主成分Y；和第二 主成分丫2第三主成分丫3作为标准化数据X* =(X；,X；,X8)的主成分讨论问 题,剩余的主成分对本文问题的讨论影响不大,因此可以忽略.由此可得T

18、的特征值,(1, ；,Q所对应的前三个正交单位化特征向量为:8=(-0.1650, 0.0770, -0.3002, -0.3102, -0.3402, 0.4444, 0.3371, 0.4296,0.4112) e； =(-0.505Q -0.5781 0.2657, 0.4226 -0.1188 0.1328 0.3496 -0.0785 -0.0282)仓=(-0.3345,0.0539, -0.5073, -0.2295,0.5166,0.1100, 0.3129, -0.2686, - 0.3552)记第i个主成分为Yi (i =1,2,3),那么标准化数据X*的第一主成分=X*

19、$、 第二主成分丫2二X* e2、第三主成分Y3 = X* ©3第九主成分= X* 9 o 由Y 与 Xi的相关系数(8)利用Matlab软件代入数据计算出前三个主成分与各原始变量的相关系数如表2相关系数卩丫恐严Y/y2*丫3*X/-0.3593-0.8135-0.4065X2*0.1677-0.93010.0655X3*-0.65380.4274-0.6164x4*-0.67550.6799-0.2789X5*-0.7407-0.19110.6278X60.96770.21370.1337X70.73410.56250.3802X80.9355-0.1263-0.3265X90.8

20、955-0.0454-0.4317所示：表2相关系数记样本均值向量X =Xi,X2,X3,乂9为亠,18,%,那么再利用方差计算公式1 n 261 a M -叫2代入数据可以得到样本向量n -1二=(0.0986 0.0354 0.0434 68.3480 54.4760 23.5925 3.2110 2.9390 0.3475)X =X1, X2, X3,X9的方差向量,- 88 ,99 )为:11?二 22,二=0.0003 0.0001 0 0.7820 0.2917 110.5272 3.0058 2.8425 0.0437所以标准化数据X*的第一主成分Y、第二主成分丫2和第三主成分弋

21、可以用样本数据X =X1,X2, X3, ,X9表示为：Y =xe2_( X1 X2- J 01V °2X9- 9 , )e2.°99e3X1 X2- 珂J'. 02,孕)e3. 099由于丫1、丫2和Y/特征值的大小决定它们奉献率的大小,所以结合Y1、Y2和丫3 的特征值线性组合后可得衡量各地区消费支出的综合指标丫 ：丫 =咯 比3=7 4.7420< 2.5884Y21.4768丫39把Y1、丫2“和Y/的值代入9式即可得到Y,根据综合指标Y的大小便可对5个城市进行排序,从而得知各个城市的污染程度：表3三个主成分值、综合指标值及排名Yi丫2丫3丫排名北京-

22、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综合考虑大气污染、水体污染和噪声污染的各项指标： 可吸入颗粒物、二氧 化硫、二氧化氮、道路交通噪声、区域环境噪声、化学需氧量、氨氮、总氮和总 磷,五座城市中重庆市污染程度最为严重,上海次之,北京、天津污染程度较为 接近,西安市污染程度较轻.5.1.2结论分析通过对北京、天津、上海、重庆和西安五座城市各项指标

23、年平均数据统计和 计算,得出结论：由于标准化后的第一主成分 丫的奉献率已达52.69%,大于其 他主成分的奉献率,可以认为丫1值越大,环境质量越差.因此观察 丫与各项指 标间的相关系数可知,水体污染的四项指标化学需氧量x6、总氮x8、总磷x9和氨氮X7对环境污染的影响较大,其中化学需氧量对人们生活环境造成的污染程 度最大.化学需氧量是衡量水质受有机物污染程度的综合指标,当其值过高时,水体会发黑、发臭.总氮是反映水体富营养化的主要指标； 过量磷可造成水体污秽异 臭,湖泊发生富营养化和海湾赤潮现象；饮用水中氨氮过量在一定条件下可转化 为亚硝酸盐,是强致癌物质,对人体健康极为不利；水体中氨氮过量那么

24、可导致水 生物表现亢奋甚至死亡.综合考虑大气污染、水体污染和噪声污染的各项指标, 得出结论：五座城市 的总体污染程度依次为：重庆、上海、北京、天津、西安.重庆市由于三面环山 的特殊地理位置和重工业兴旺、煤烟排放过量等原因,使得其总体污染程度在五 座城市中居于首位；上海机动车尾气排放量常年居高不下,污染程度仅次于重庆； 北京和天津地理位置接近,其工业污染物排放量较高,造成大量的雾霾扬尘天气； 而西安汽车尾气和工业废气排放量较大,因此空气污染最为严重,但污染程度较其他城市略为稍微.5.2问题二的模型建立与求解由于人类活动的影响,大气、水体及噪声等污染物浓度及其分布发生了很大 变化,了解各污染物的空

25、间分布及扩散情况对评价其污染程度有很大帮助.一方面,确定污染扩散的方式、速率和程度有助于限制其影响范围,提早预防和改善, 另一方面可确定污染源的位置,从根本上杜绝或减轻污染的发生和污染程度,为相关部门提供污染治理的意见和建议. 以北京市为例,建立其污染扩散模型并以 空气污染为检验标准验证模型的合理性.5.2.1估计北京市污染源位置查阅并搜集北京市污染空气的SO2、NO 2和PM10 可吸入颗粒的年平均浓度含量,取包含北京市在内的长为183.4千米,宽为170.1千米的长方形区域 如图1,建立二维坐标系,得到北京市 14个区的地理位置坐标x, y单位:千米如表5所示：表4 2021年北京市各项污

26、染指标值及地理位置坐标单位mg /m3/万吨SO2NO 2PM10PM2.5横坐标x纵坐标y朝阳区29.764112.491.3132.191.9海淀区26.963.611598.1126.576.8丰台区28.157.5118.596.9135.373.7石景山区24.963.3116.492.8130.166.5门头沟区24.651.8114.891.1128.361.2房山区31.261.9131.7106.8148.562.3通州区38.655.8123.5105.7130.7101.4顺义区20.844.898.584.8104.7104.3大兴区33.765.7130.3107.8

27、149.778.8昌平区25.943.594.779.294.969.3平谷区20.63598.784.8101.4145.4怀柔区22.337.995.376.183.1103.1密云县21.343.685.971.676.9117.4延庆县19.234.478.368.066.547.1实用标准文档根据上表数据利用Matlab软件分别画出SO2、NO2和PM10的年平均浓度 随地理位置变化的空间浓度分布图：(11)文案大全图2 SO2的空间浓度分布图z图3 NO2的空间浓度分布图150o o o o o o O3 2 1 0 9 8 7图5 PM2.5的空间浓度分布图图4 PM10的空间浓

28、度分布图观察图2、图3、图4和图5估计污染指标SO2的污染源位置在(110,100)附 近,NO 2的污染源位置在(130,83)附近,可吸入颗粒PM10污染源位置在(132,80) 附近,细微颗粒PM2.5的污染源位置在(125,85)附近.5.2.2 建立Cauchy扩散模型设c(x, y, z,t)表示t时刻(x, y,z)处某一种污染物的浓度. 在污染物传播的区 域内,任意取一封闭曲面 S,它所围区域为门,在氏时间内,由于扩散,经过 这个区域内的污染物的质量'=m1为2ccos-：i、jy2 ccos" 1 - z2 c cos )dSdt ( 10):x: y: z

29、x, y,z轴方向的扩散系数,cos,cos ： ,cos是方向Hii 二.t . .(:-S 其中,m分别代表沿 余弦.由Gauss-Green公式得t过：Q-2 -22 ；一 C 2 ：一 Ccr +cr +crx - 2 y -:x: y-22 C2 )dxdydt :z扩散过程中由于沉降等原因导致 区域内质量减少量m2为=m2k c(x, y, z,t)dxdydzdtt Q其中k表示衰减系数.t时间i "区域内污染物质量的变化m为实用标准文档文案大全m 11 ic(x, y, z,t t) 一c(x, y,乙t)dxdydz Q=“眾 fdxdydzdtQ ctt时间O区域

30、内污染物质量的变化 Am也等于这段时间内增加的污染物质 量减去减少的污染物的质量,二m - . ：mi| -. ：m2t沁2 22 ：'C2 ；'C2；y r+c；x：y2 ? C 22 -k c(x, y, z,t)dxdydzdt :z.(Q_c兰 dxdydzdtt、区域门的任意性,可以得到2 2 2:C2 ；- C 2 C 2 ；- C 2x2 Yy22-k c A:x；y: z记污染源在点(x0,y0,zo)处,那么在t = 0时刻有c(x,y, z,0) =M、(x -x.)、(y - y.)、(z 二)其中M为污染物排放量.由以上分析可知此扩散问题满足 Cauch

31、y问题由于所选时间t、(12)O-2 -22 C 2 C 2 y 7z 2 _ k C5-.2CC2 C C77 =7TT + 口 y - 2 . - Z . 2ytexcyczC(x, y,z,0) =M 6(x x0)6(y y0)6(zZ0)对(13)进行傅里叶(Fourier)变换,且令 =(1, 2, 3),2u(13)t, )= ：u(x, y, z,t),_;c珂故得常微分方程唯一解对(14)求逆变换由于_ x ； C(x,y, z,0) 1二 M ； L (x - X.)！ ； L (y - y°)! ； L (y - yo) 1二 Me 4( 1 x0 2yo 3z

32、0)Cauchy问题眩+&皆+叩計或鶯+k2)C = 0* dt&0j) = Me4"% 旳申C(t巧=Me&邛弋少£92*2)1(叫勺0吧G-aC,(14)12le益CXx2二e故得c(x, y,z,t) mjT(C)M(2石T)'店；晞;Mexp_ k2t丄4 二 yt4；zt4 ft2 2 2: exp(8 二 t；x；y；z * 二t)4；xt4；yt4；zt为了简化模型,考虑空气污染的平面传播,那么可得到不同空气污染指标的扩 散方程：c(x,y,t)=-exp -3-M-k2t(15)(8肮6巧、4't加 yt如果认为经过了

33、相当长的时间后,扩散已经终止,物质分布处于平稳状态,c那么方程(12)中的-=0,于是有线性椭圆方程的边值问题-2-2 -2 J C 2 J C . 2Qrx7TT+Qry7-k C = 0«cXcy&(x, y)倉®(x, y)也就是用傅里叶(Fourier)变换求解,其中 戸丨是区域i丨的边界.在Cauchy问题的解(14)中,有三个未知的参数二x&y,k,它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根,假设点源的质量M与位置(xo,y.),设观测取样为：(捲,口1),区,y?, m?),川(x., yng),取样时刻为t =1,或设t二t&

34、#176;. , t.是取样时间,那么(12)变成:Ct = t.xCxx t0 yCyy -tok C对.而言,取样时间为1,而方程形状与(12) 一致,把在(Xj,yj点观测到物质密度mi与公式(15)都取对数,令t =1那么有In c(x, y,t) =ln3 In ；x；y2(x-X.)十(y-y.)2-k2(16)令u =So£ ,V 5oX2 ?-x12cy2ln(2厂厂1,那么(16)可以写成 In c(x, y,1)=:U :V ;.当 x = X., y = y.时,有 In c(x., y.,1)=；.利用Matlab软件对传播模型的各项参数进行估计并用于检验模型

35、的合理 性,取空气污染的取样点数据,运用非线性拟合nlinfit函数对数据进行拟合,得出函数参数即可得污染源位置.Matlab求得的空气污染指标SO2污染源位置为(115.5,97.3),NO2污染源位 置为(128.2,85.3),可吸入颗粒PM10的污染源位置为(129.1,80.8),PM2.5的污 染源位置为(125.5,87.6),与实际污染源位置很接近,故Cauchy污染传播模型建 立合理.5.3问题三的模型建立与求解根据相关数据西安市在某个时间段内 PM2.5 一直高居全国前列,因此空气污染是西安市人民关注的主要方面.采用时间序列模型对先市空气污染的各项指标进行预测. 通过平稳性

36、分析和 模型定阶的AIC准那么确定预测模型,再利用最大似然法估计各级参数,最后进行 模型预报.531时间序列模型的建立1. 平稳性检验(1) 自相关函数MA(q)序列(滑动平均序列)的自相关函数：设1 , X2,X t f 是 MA (q)序列,Xi = (B)-. = ；t jR j Vt _2 -31 q 討 _q.其中；t是零均值,即平稳白噪声,且Var(；t)=：；t2,那么Xt是零均值平稳序列,且 其自协方差函数和自相关函数分别为石；(1中日12 +日22十+材),k二二h 可入1皿),0/k = 0,1 _ k _ q, k q；(17)k =0,(18)(19)(20)1,-日k

37、 +日1日心+3q"q2 21 +01 + +0qP,其自相关函数呈拖尾状趋于零,偏自相关函数式非截尾但趋于零.AR(p)序列(自回归序列)的自相关函数：t是平稳白噪设 Xt是AR(p)序列：X1 lXt'Xt,：pXt斗；t,声,且Var( ；t)=门.那么其自协方差相关函数和自相关函数为% =E(X* 化餐1X+%Xt“ + 哥)其自相关函数呈拖尾状趋于零,偏自相关函数非截尾但趋于零.而ARMA (p, q)的自相关函数可以看做是MA (q)序列的自相关函数和AR(p)序列的自相关函数的混合物.计算知本时间序列服从ARMA (p,q)模型,下面考虑定阶,即估计p、q的值(

38、2) 模型定阶的AIC准那么AIC准那么又称Akaile信息准那么,AIC准那么是信息论与统计学的重要研究成 果,具有重要意义.AIC准那么起源于Kullback-Lelibler 信息量,设Xt是随机变量,它的概率 密度是f(x),其中含有k个未知参数,设未知参数向量为一：=('2,-k)T. f(x)属于分布族 g(x| '),其中 2匕,-k)T.显然 f(x)二 g(x|").K-L信息量可以用来刻画g(x| ：)与f (x)的接近程度,其定义为I(f(),g(|» f(x)lnfhdx那么有 i(f(),g()一0,且有 i(f(),g( I J)

39、lj = o.K-L信息量是寻求最接近f(x)的参数概率密度g(x| -),使得I (f (), g( | 0) =min(21)经过理论分析,当给定样本观测值X =(为,乂2,xn)(样本容量为n ),设k '-?(m) 是模型参数=(2,GT(未知参数个数是k )的最大似然估计.设In (L(J) 是其对数似然函数,AIC信息准那么是：使得(21)式中的k满足AIC(k) »2Inga) 2k=min设Xt是ARMA( p,q)序列,其中未知参数的个数是 p q 1个,包括自回(22)归参数,=(1,2,;：p)T,滑动平均参数二=厲门2厂Jq)T及匚2.对于对数似然 函

40、数1Xn) = TnMn2L(点2-(2二)- 二2 S2 21n略去2ln |M丨一,n(2二),得最小平方和估计对应似然函数L(点2|x)n 二2S(：)(23)由(23)式可得,：?2的最大似然估计；?2S(?)n将上式代入(23),得,',：2, 、 n,2 nL(十；|x) 才n 一?因此,ARMA( p,q)序列AIC定阶准那么为：选p,q,使得AIC = nln 貯十 2( p+q+1) = min( 24)其中n是样本容量,er2与p,q有关.假设当p=p,q=q,使得(24)到达最小值, 那么认为序列为ARMA( ?,(?).当ARMA( ?,?)序列含有未知均值时,

41、模型为这时,未知参数个数为p q 2,AIC准那么为：选取p,q,使得(25)2AIC 二 nln 二；2(p q 2) =min其式(24)与式(25)有相同的最小值点p,q.2. 最大似然估计法(ML估计)估计参数假设乂必2,Xt ?是零均值正态ARMA (p,q)序列,他的概率密度为fXn Pn62 =2兀尸厶力冷X" nxn)其中=(1 ,2,p,R,如丁风二(Xi,X2 /人)丁,、' n 二E-(XnXn)E-表示参数向量是1时取期望(因Xt是ARMA (p,q)序列,Xn与有关)令M二匚2广,那么有对数似然函数仙点2-)*®誓)L点2；|x= In f

42、x| 二二2；xtMx1 n1ln | M |ln2二ln 二 22 22"其中X =(为,X2, Xn)为样本观测值.下证xtMx仅与1有关,而与2无关,事实上的第i行第j列的元素是0oO5 =E(XiXJ=r4|,设Xt的传递形式是XZ Gjgg,那么帝|=(送GkG讯2 j =0k=0由此即知M =O"2区与石2无关,从而xtMx仅与卩有关,而与石2无关.由于xtMx仅与：有关,记S(J二xtMx(26)那么称S( J为平方和函数,进一步记？=E(JX1,X2,Xn),即?是,关于X1,X2/ ,Xn的条件期望,也就是X1,X2/ ,Xn时关于；t的估计,那么可证nS

43、(J 八?t=jQO可见S( J是一残差平方和形式,假设?S使(26)式到达最小,即S( J = min,那么?S称为最小平方和估计.现讨论最大似然估计,对数似然函数是1 n1S( P)L(:,二2 | x) =l n f (x| :卢2) =T n |M |- l n(2二)- l n 2 - 一 ？)( 27)2 222.f实用标准文档由于M , S(P)均与er：无关,令解得2 S( ) xT Mx将其带入(27)中得 In f(x|宾升夕nlMIln(xTMx) c 其中-nln(2二)1lnn,略去常数c,得(29)1nL( -) In | M | 一才n(xTMx)-的最大似然估计

44、应使(29)到达最大.假设?m使得L(?m) =max,那么?m是一：的最大似然估计(ML估计),代入(28) 即得:?2的最大似然估计值.3. 模型考核的2检验假设拟合模型的残差记为？,它是引的估计值.例如对AR(p)序列,设未知参数的估计是鶴角,陽,那么残差？t =Xt -昭-购pX j,t = 1,2,nk =1,2,m500文案大全Ljung-Box的2检验统计量是假设检验H o：k = 0,当k乞m i比：对某些k乞m,k = 0.在H 0成立时,假设 n充分大,2近似于2(m-r)分布,其中r是估计的模型参数个数.2检验法：给定显著水平：,设实际算得的2值是02 ,P值为p=p&#

45、39; 2一 02 I 那么当P ： :,拒绝H0,即认为；t非白噪声,模型考核不通过,需重新估计参数； 而当p -,接受H.,认为“是白噪声,模型考核通过,可以进行预报.4. ARMA (p,q)的预报时间序列的l步预报,是根据X k, X k 4厂 是取值对未来k l时刻的随机变 量Xk .(l 0)做出估计,估计量记作X?k(l),它是txk,xk？的线性组合.*0,X1,X2, ,x是零均值正态ARMA (p,q)序列,那么以下预报差分方程成 立江1忍(-1) ；*k(l -2) pXUl -q),(l p) (30)只需知道Xki,Xk2, >?kq,就可以递推算得z?kl,l

46、P,因此定义预报变量江(q) =(X,X(2),X(q)T,令可证以下递推预报公式:仰.=<0,= 1,2/ , p；P._G1100G20a1 亠a0a-Gq000-Gp +%(p * 十q A.(p *屮q_2cp 十2G -00G20 101 /;X?k(q)GqGq0Xk“1式中第三项当p乞q时为零.且当k.较小时,令初值X.T =0.5.3.2模型求解及结论考虑数据的连贯性,选取2021.4.1日至2021.7.31日的空气质量指标AQI数 据见附录表2进行时间序列分析,运用 Matlab软件作出相应时序图,进行时 间序列的验证：Q丰250200150100M VIi| i11

47、:1'I1 % JfiV1 A u. Ai/ I匚u '20406080100120140t图6空气质量时序图实用标准文档以及相应的自相关系数图与偏自相关系数图:-0.20 2468101214161820LagSample Autocorrelation Function8iT1I0 ;246581012416 -8 2-0.2LagSample Partial Autocorrelation Function8G 4-2o a a a a、/fa OCQ - an Arp mas图7自相关系数图表8偏自相关系数图观察图像明显可知该序列不平稳,序列图.Q肌200 -故作差分处

48、理.图9为作一阶差分后所得nnrLeFOCOIUA toHpmao-100020406080 100 1201曲图9 一阶差分后的空气质量时序图为进一步验证平稳性,考察差分后序列自相关图和偏自相关图.Sample Autocorrelation FunctionSample Partial Autocorrelation Function1 . 1 -8-0.402468101214161820Lag图10 一阶差分后自相关图8aeFOC nrurthpLrar Bwpmao2-0-0.402468101214161820Lag图11 一阶差分后偏自相关图文案大全观察图像可知一阶差分后的序列较

49、平稳,故可确定ARMA (p,q)模型利用AIC准那么定阶：当d =1时,计算得到AIC值如下表所示：表5 AIC值AIC值q =0q = 1q=2q = 3p =01146.961128.281123.941120.14p =11139.551115.621117.481118.16p=211129.591117.5611118.421120.42p = 31129.081117.291119.281119.83由表中数据知当P =1 , q =1时,其AIC值最小,故取ARMA (1,1)模型.又 由2检验法,取显著性水平 0.05 ,运用Matlab软件计算得到 p =0.72=0.05

50、,故该模型通过考核检验,可以用来预报.进一步计算相关参数得到？1二0.4891,:?1,故得到ARMA (1,1)模型为(1 _0.4891B)(1 _ B)Xt =(1 B) ；t其十步预报值见下表：表6 8月1-10日空气质量AQI预测值与实际值比拟8月1日2日3日4日5日6日7日8日9日10日预报值85.9274.4074.2768.2570.8666.2369.5265.2568.7064.54真实值12311311293797447396383进一步对8月1-10日各项污染指标PM2.5、PM10、SO2、NO2和CO的浓 度进行预测,预测值如下表：表7 8月1-10日各项污染指标的预测值(单位：mg/m3)PM2.5PM10SO2NO 2CO8月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日