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

城市空气污染程度的分析和预测 摘 要 本文讨论了有关城市污染程度 污染因素及污染扩散的问题 对于问题一 本文主要从大气污染 噪声污染和水体污染这三个面选取主 要污染物 查阅北京 天津 上海 重庆和西安五座城市 2007 2012 年的年度 平均污染数据 采用降维的思想 运用主成分分析法减少变量个数 再借助 Matlab 软件计算各主成分的贡献率 分析知可选取前三个主成分作为衡量污染 程度的标准 最后根据综合指标得到这五个城市的污染程度从高到低依次为 重庆 上海 北京 天津 西安 通过判断相关系数的大小 确定五个城市影 响人们生活的主要污染因素是水污染 其四项指标依次为化学需氧量 总氮 总磷和氨氮 对于问题二 以北京市大气污染为例 首先 利用 GPS 记录北京市 14 个城 区观测点的位置 并查阅 2013 年污染指标 与的污染 2 SO 2 CO5 2PM10PM 数据 绘制出相应的空间浓度分布图 估计这四种污染物的大致污染源位置依 次为 附近 附近 附近和附近 其次 根 100 110 83 130 85 125 80 132 据污染扩散原理和方式 建立 Cauchy 污染传播模型 根据各地区空气污染物的 浓度分布 运用 Matlab 软件对数据非线性拟合 得出扩散模型各参数的值 计 算得出各项污染指标的污染源位置依次为 3 97 5 115 3 85 2 128 和 最后 比较污染物位置的计算值与实际值 发现误 8 80 1 129 6 87 5 125 差相差较小 故模型建立较为合理 对于问题三 分析西安市的主要污染 大气污染 收集西安市 2014 年 4 1 7 31 日的空气污染数据 根据时间序列的平稳性特点及定阶准则选取AIC 合适的时间序列模型 利用 Matlab 软件对序列模型的各项参数进行 11 ARMA 估计并检验模型的合理性 并将模型用于数据预报 利用时间序列模型预测西 安市未来 10 天的空气污染状况总体等级为良 对于问题四 基于问题一 二 三对污染因素的分析和污染扩散的特点 主要从减少污染物的产生和治理净化已产生的污染物两方面 针对大气污染 水体污染和噪声污染为相关部门提供合理化防治建议 关键词 主成分分析 Cauchy 污染传播模型 时间序列模型 Matlab 软件 一 问题重述 环境质量的好坏对人们的生活有着最直接的影响 包括身体健康 寿命长 短 生活质量以及经济发展水平 生态环境等 因此保护环境的工作需要个人 集体 国家 乃至全球各国的共同努力 事实上 环境保护事业从 1962 年起步 至今 一直备受人们关注和重视 环境污染的防治工作也一步步开展运营 然 而较经济发展的速度和水平 还有很多需要为之努力的工作要做 环境保护工 作一方面是要从减少污染物的产生方面 预防和控制污染源头污染物的排放 另一方面是治理净化已产生的污染物 环境污染是指环境因受人类活动影响而改变了原有性质和状态的现象 其 实质是人类活动中将大量的污染物排入环境 影响其自净能力 降低了生态系 统的功能 为考察某时期某区域的环境状况 需统计该地区各项污染指标的数据 根 据各污染指标的污染排放标准 确定其污染程度 并通过污染发生的原因 途 径和后果制定合理的污染治理方案 加以实施运行 改善环境 统计北京 天津 上海 重庆和西安五座城市的污染数据 讨论如下问题 问题一 比较五座城市的污染程度及影响人们生活的主要污染因素 问题二 建立城市污染扩散模型确定某城市污染程度并检验 问题三 给出西安市某方面污染状况的中长期预测数据 问题四 给相关部门提供一篇治理污染的建议短文 二 问题分析 对于问题一 本文将查阅收集北京 上海 天津 重庆以及西安五座城市 大气污染 噪声污染和水体污染 2007 2012 年的数据 由于各项污染指标众多 所以选取具有代表性的污染物 可吸入颗粒物 道路交通噪10PM 2 SO 2 CO 声 区域环境噪声 化学需氧量 氨氮 总氮和总磷这 9 项指标来综合衡量各 城市的污染程度 本文将运用降维的思想 利用主成分分析法计算前几个主成 分的贡献率大小及通过线性组后的综合指标 以此来衡量各城市的污染程度 根据各主成分的相关系数大小即可确定影响人们生活的主要污染因素 1 对于问题二 本文将以北京市的大气污染为例 收集北京市 2013 年 14 个 城区的 和的污染数据 绘制出这四种污染物浓度的空 2 SO 2 CO5 2PM10PM 间浓度分布图 估算出这四种污染物的大致污染源位置 建立 Cauchy 污染传播 模型 根据各地区空气污染物的浓度分布 运用 Matlab 软件对数据 nlinfit 非 线性拟合 计算出扩散模型各参数值 得出各空气污染物的污染源位置 比较 污染物位置的计算值与实际值 检验污染传播模型建立是否合理 2 对于问题三 本文将考虑西安市大气污染的各项污染指标 搜集 2014 年 4 月 1 日至 2014 年 7 月 31 日的污染数据 根据时间序列的平稳性特点及定AIC 阶准则选取合适的时间序列模型 利用 Matlab 软件 3 估计序列模型的各项 参数并检验模型的合理性 考核通过时方可用于数据预报 对于问题四 基于问题一 二 三对污染因素的分析和污染扩散的特点 主要从减少污染物的产生和治理净化已产生的污染物两方面 针对大气污染 水体污染和噪声污染为相关部门提供合理化防治建议 三 模型假设 1 所有城市污染数据均真实有效 具有统计价值 2 污染数据的取样为瞬时取样且在某个平面内 不考虑取样前后污染物的浓度 变化 3 绘制北京大气污染浓度分布图时不考虑海拔因素的影响 4 大气污染过程经历着扩散与衰减 且沿轴方向上的扩散系数分别为常数 yx 衰减与浓度成正比 5 不考虑日照等自然因素对各项指标污染浓度的影响 四 符号说明 4 1 模型一 可吸入颗粒物 单位 1 x 3 mmg 二氧化硫 单位 2 x 3 mmg 二氧化氮 单位 3 x 3 mmg 道路交通噪声 单位分贝 4 x 区域环境噪声 单位分贝 5 x 化学需氧量 单位万吨 6 x 氨氮 单位万吨 7 x 总氮 单位万吨 8 x 总磷 单位万吨 9 x 样本协方差矩阵 协方差矩阵的特征值 i 协方差矩阵的特征向量 i e 相关系数 ii XY 样本方差向量 i 样本相关矩阵 样本相关矩阵的特征值 i 样本相关矩阵的特征向量 i e 样本均值向量 i 主成分 i Y 污染程度综合指标 Y 其中 9 2 1 i 4 2 模型二 位置坐标 横轴 x 位置坐标 纵轴 y 位置坐标 竖轴 z 污染物浓度 c 污染物质量 m 污染物排放量 M 扩散系数 4 3 模型三 序列 t X 平稳白噪声 t 滑动平均序列 MA q 自回归序列 AR p 自回归滑动平均序列 ARMAqp 常数 c 自协方差函数 k 自相关函数 k 概率密度 xf 未知参数向量 i 最大似然估计 m k 自回归参数 l plkit 2 1 2 1 2 1 0 五 模型建立与求解 5 1 问题一的模型建立与求解 环境因受人类活动影响而改变了原有性质和状态的现象称为环境污染 环 境污染的实质是人类活动中将大量的污染物排入环境 影响其自净能力 降低 了生态系统的功能 然而生活中环境的污染也是多方面的 污染按环境要素分 类 可分为 大气污染 土壤污染和水体污染 为此 本文选取环境污染的主要污染影响指标 建立污染评价模型 即主 成分分析模型 比较北京 天津 上海 重庆和西安五座城市的污染程度 并 得出影响人们生活的主要污染因素 5 1 1 建立主成分分析模型 首先选取环境污染的部分污染指标 可吸入颗粒物 二氧化硫 二氧 1 x 2 x 化氮 道路交通噪声 区域环境噪声 化学需氧量 氨氮 总氮 3 x 4 x 5 x 6 x 7 x 和总磷 北京 天津 上海 重庆和西安五座城市的各项指标年平均数据 8 x 9 x 见附录表 1 问题中变量太多不但会增加计算的复杂性 也给合理的分析问题和解决问 题带来很大困难 同时 这些变量之间存在的相关性 也使得这些变量所反映 的信息在一定程度上有所重叠 为了减少变量个数 本文采用了降维的思想 利用主成分分析法得出污染空气的主要因素变量 提高问题研究的合理性 同 时不会使数据反映的信息量有大的损失 设是的协方差矩阵 的特征值与正交化特征向量 T 321 p xxxxx 分别为及 且的第 个主成分为 0 321 p p eeee 321 xi 3 2 1 332211 pixexexexeY pipiiii 1 根据平均值计算公式代入数据计算样本的均值向 n 1j i 1 ij x n x 321p xxxxx 量为 321p xxxxx 3475 0 9390 2 2110 3 5925 234760 543480 680434 0 0354 0 0986 0 x 根据协方差矩阵计算公式 T 1 1n 1 xxxx ii n i 2 利用 Matlab 软件代入数据可求得随机变量 相应 321p xxxxx 9 p 样本协方差矩阵为 只写下三角部分 0437 0 3485 0 1690 07574 1 1029 00963 0 0005 00003 0 0003 0 8425 2 4350 18003 147895 0 9387 00049 00039 0 0028 0 0058 30630 163881 0 3351 00051 00062 00281 0 5272 1108197 30768 5 0409 00031 01114 0 2917 0 0891 0 00003 0 0020 0 7820 0 0050 0 0064 00034 0 000 0001 0 0001 0 0003 0 由于协方差矩阵对角线元素即为样本的方差值 观察矩 321p xxxxx 阵易知其方差值相差较大 即各量纲间分散程度较大 不适用于环境污染模型 中各成分指标间的相关性分析 为消除由于量纲不同可能带来的结果误差 故 将样本变量标准化 即令 321p xxxxx 9 p 3 9 2 1 i x X ii ii i 得到一组新的数据 其中为的平均值 为的方 921 XXXX i i x ii i x 差 此时 的协方差矩阵即为样本变量 921 XXXX 的相关矩阵 其中 921 xxxx ppij 4 jjii ji jiij XX YX Cov E 为的协方差 ov j XXC ij XXi 则根据公式 4 运用 Matlab 软件代入数据即可得其样本相关矩阵为 19891 04664 0 7998 09114 0 5213 0 3486 0 1325 00887 0 14909 0 8350 0 8670 0 6296 0 4691 0 2378 00908 0 18813 0 4144 0 2186 0 4726 0 3709 0 8787 0 16727 0 5461 06245 0 0229 0 5744 0 11866 00013 00485 01965 0 19102 0 7446 0 2097 0 15158 0 1171 0 16026 0 1 由特征值计算公式 5 0 E 第 个主成分的贡献率 i 6 p i 及前个主成分的累计贡献率 k 7 p k s s 1 其中 计算的特征值及各标准化主成分的贡献率和累计贡献率如表 1 所9 p 示 表 1 的特征值 贡献率和累计贡献率 i i 贡献率 累计贡献率 14 742052 6952 69 22 588428 7681 45 31 476816 4197 86 40 19282 14100 500100 600100 700100 800100 900100 由表 1 易知 第一 第二和第三主成分的累积贡献率为 97 86 并且只有 第一 第二和和第三主成分对应的特征值大于 1 所以选取第一主成分和第 1 Y 二主成分第三主成分作为标准化数据的主成分讨论 2 Y 3 Y 8 2 1 XXXX 问题 剩余的主成分对本文问题的讨论影响不大 因此可以忽略 由此可得的特征值所对应的前三个正交单位化特征向 921 量为 4112 0 4296 0 3371 0 4444 0 3402 0 3102 0 3002 0 0770 0 1650 0 1 e 0282 0 0785 0 3496 0 1328 0 1188 0 4226 0 2657 0 5781 0 5056 0 2 e 3552 0 2686 0 3129 0 1100 0 5166 0 2295 0 5073 0 0539 0 3345 0 3 e 记第 个主成分为 则标准化数据的第一主成分 i 3 2 1 iYi X 1 1 eXY 第二主成分 第三主成分第九主成分 2 2 eXY 3 3 eXY 9 9 eXY 由与的相关系数 i Y i X 8 iji XY e ii 利用 Matlab 软件代入数据计算出前三个主成分与各原始变量的相关系数如表 2 所示 表 2 相关系数 相关系数 iI XY 1 Y 2 Y 3 Y 1 X 0 3593 0 8135 0 4065 2 X 0 1677 0 93010 0655 3 X 0 65380 4274 0 6164 4 X 0 67550 6799 0 2789 5 X 0 7407 0 19110 6278 6 X 0 96770 21370 1337 7 X 0 73410 56250 3802 8 X 0 9355 0 1263 0 3265 9 X 0 8955 0 0454 0 4317 记样本均值向量为 则 9321 xxxxx 9821 3475 0 9390 22110 3 5925 234760 543480 680434 00354 00986 0 再利用方差计算公式代入数据可以得到样本向量 n iiii x n 1j 2 j 1 1 的方差向量为 9321 xxxxx 99882211 0437 0 8425 2 0058 3 5272 1102917 0 7820 0 00001 0 0003 0 所以标准化数据的第一主成分 第二主成分和第三主成分可以用 X 1 Y 2 Y 3 Y 样本数据表示为 9321 xxxxx 1 99 99 22 22 11 11 11 e x x x eXY 2 99 99 22 22 11 11 22 e x x x eXY 3 99 99 22 22 11 11 33 e x x x eXY 由于 和特征值的大小决定它们贡献率的大小 所以结合 1 Y 2 Y 3 Y 1 Y 和的特征值线性组合后可得衡量各地区消费支出的综合指标 2 Y 3 YY 9 321332211 4768 1 5884 2 7420 4YYYYYYY 把 和的值代入 9 式即可得到 按照综合指标的大小便可 1 Y 2 Y 3 YYY 对 5 个城市进行排序 从而得知各个城市的污染程度 表 3 三个主成分值 综合指标值及排名 1 Y 2 Y 3 Y Y 排名 北京 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 结论分析 通过对北京 天津 上海 重庆和西安五座城市各项指标年平均数据统计 和计算 得出结论 由于标准化后的第一主成分的贡献率已达 52 69 大于 1 Y 其他主成分的贡献率 可以认为值越大 环境质量越差 因此观察与各 1 Y 1 Y 项指标间的相关系数可知 水体污染的四项指标化学需氧量 总氮 总磷 6 x 8 x 和氨氮对环境污染的影响较大 其中化学需氧量对人们生活环境造成的污 9 x 7 x 染程度最大 化学需氧量是衡量水质受有机物污染程度的综合指标 当其值过高时 水 体会发黑 发臭 总氮是反映水体富营养化的主要指标 过量磷可造成水体污 秽异臭 湖泊发生富营养化和海湾赤潮现象 饮用水中氨氮过量在一定条件下 可转化为亚硝酸盐 是强致癌物质 对人体健康极为不利 水体中氨氮过量则 可导致水生物表现亢奋甚至死亡 综合考虑大气污染 水体污染和噪声污染的各项指标 得出结论 五座城 市的总体污染程度依次为 重庆 上海 北京 天津 西安 重庆市由于三面 环山的特殊地理位置和重工业发达 煤烟排放过量等原因 使得其总体污染程 度在五座城市中居于首位 上海机动车尾气排放量常年居高不下 污染程度仅 次于重庆 北京和天津地理位置接近 其工业污染物排放量较高 造成大量的 雾霾扬尘天气 而西安汽车尾气和工业废气排放量较大 因此空气污染最为严 重 但污染程度较其他城市略为轻微 5 2 问题二的模型建立与求解 由于人类活动的影响 大气 水体及噪声等污染物浓度及其分布发生了很 大变化 了解各污染物的空间分布及扩散情况对评价其污染程度有很大帮助 一方面 确定污染扩散的方式 速率和程度有助于控制其影响范围 提早预防 和改善 另一方面可确定污染源的位置 从根本上杜绝或减轻污染的发生和污 染程度 为相关部门提供污染治理的意见和建议 以北京市为例 建立其污染 扩散模型并以空气污染为检验标准验证模型的合理性 5 2 1 估计北京市污染源位置 查阅并搜集北京市污染空气的 和 可吸入颗粒 的年平 2 SO 2 NOPM10 均浓度含量 取包含北京市在内的长为 183 4 千米 宽为 170 1 千米的长方形 区域 如图 1 建立二维坐标系 得到北京市 14 个区的地理位置坐标 单位 千米 如表 5 所示 yx 图 1 北京市各城区图 表 4 2013 年北京市各项污染指标值及地理位置坐标 单位 3 mmg 万吨 2 SO 2 NO PM10PM2 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 8149 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 软件分别画出 和的年平均浓度 2 SO 2 NOPM10 随地理位置变化的空间浓度分布图 0 50 100 150 0 50 100 150 15 20 25 30 35 40 45 x y z 0 50 100 150 0 50 100 150 30 40 50 60 70 x y z 图 2 的空间浓度分布图 图 3 的空间浓度分布图 2 SO 2 NO 0 50 100 150 0 50 100 150 70 80 90 100 110 120 130 x y z 0 50 100 150 0 50 100 150 70 80 90 100 110 x y z 图 4 的空间浓度分布图 图 5 的空间浓度分布图PM10PM2 5 观察图 2 图 3 图 4 和图 5 估计污染指标的污染源位置在附 2 SO 100 110 近 的污染源位置在附近 可吸入颗粒污染源位置在 2 NO 83 130 PM10 附近 细微颗粒的污染源位置在附近 80 132 PM2 5 85 125 5 2 2 建立 Cauchy 扩散模型 设表示 时刻处某一种污染物的浓度 在污染物传播的 c x y z tt x y z 区域内 任意取一封闭曲面 它所围区域为 在时间内 由于扩散 经S t 过这个区域内的污染物的质量为 1 m 10 tt t S zyx tS z c y c x c mdd coscoscos 222 1 其中 分别代表沿轴方向的扩散系数 是方向 222 xyz x y zcos cos cos 余弦 由 Gauss Green 公式得 tyx z c y c x c m tt t zyx ddd 2 2 2 2 2 2 2 2 2 1 扩散过程中由于沉降等原因导致区域内质量减少量为 2 m 11 2 2 d d d d tt t mk c x y z tx y z t 其中表示衰减系数 2 k 时间区域内污染物质量的变化为t m d d dmc x y z ttc x y z tx y z d d d d tt t c x y z t t 时间区域内污染物质量的变化也等于这段时间内增加的污染物质t m 量减去减少的污染物的质量 即 12 mmm 222 2222 z 222 d d d d tt xy t ccc k c x y z tx y z t xyz d d d d tt t c x y z t t 由于所选时间 区域的任意性 可以得到tt 12 222 2222 222 xyz cccc k c txyz 记污染源在点处 则在时刻有 000 xyz0t 000 0 c x y zMxxyyzz 其中为污染物排放量 M 由以上分析可知此扩散问题满足 Cauchy 问题 13 222 2222 222 000 0 xyz cccc k c txyz c x y zMxxyyzz 对 13 进行傅里叶 Fourier 变换 且令 123 c tu x y z t 已知 2 2 1 2 u c x 2 2 2 2 u c y 2 2 3 2 u c z 1 023 0 000 0 o ixyz c x y z Mxxyyyy Me 故得常微分方程 Cauchy 问题 1 023 0 2222222 123 d 0 d 0 o xyy ixyz c kc t cMe 唯一解 14 2222222 1231 0203 0 xyz kt ixyz c tMe 对 14 求逆变换由于 1 2 1 222 41 xxx x ee 2 1 222 1 00 4 1 xx ixx x x eee 故得 1 222 2 000 222 3222 222 2 000 222 exp 444 2 exp 444 8 xyz xyz xyz xyz c x y z tc xxyyzzM k t ttt t xxyyzzM k t ttttt 为了简化模型 考虑空气污染的平面传播 则可得到不同空气污染指标的 扩散方程 15 22 2 00 22 exp 44 8 xy xy xxyyM c x y tk t tttt 如果认为经过了相当长的时间后 扩散已经终止 物质分布处于平稳状态 则方程 12 中的 于是有线性椭圆方程的边值问题0 c t 22 222 22 0 xy cc k c xy c x yx y 也就是用傅里叶 Fourier 变换求解 其中是区域的边界 在 Cauchy 问题的解 14 中 有三个未知的参数 它们分别是扩散 xy k 与衰减过程中的扩散系数与衰减系数的算术平方根 假设点源的质量与位置M 已知 设观测取样为 取样时刻为 00 xy 111 x y m 222 nnn xy mxy m 或设 是取样时间 则 12 变成 1t 0 tt 0 t 222 000txxxyyy CtCtCt k C 对而言 取样时间为 1 而方程形状与 12 一致 把在点观测到物 ii x y 质密度与公式 15 都取对数 令则有 i m1t 16 22 2 00 22 3 ln lnln 44 2 xy xy xxyyM c x y tk 令 2 0 4 xx U 2 0 4 yy V 2 1 x 2 1 y 则 16 可以写成 2 3 lnln 2 xy M k ln 1 c x yUV 当时 有 00 xxyy 00 ln 1 c xy 利用 Matlab 软件对传播模型的各项参数进行估计并用于检验模型的合理性 取空气污染的取样点数据 运用非线性拟合 nlinfit 函数对数据进行拟合 得出 函数参数即可得污染源位置 Matlab 求得的空气污染指标污染源位置为 污染源位 2 SO 3 97 5 115 2 NO 置为 可吸入颗粒的污染源位置为 的污 3 85 2 128 PM10 8 80 1 129 PM2 5 染源位置为 与实际污染源位置很接近 故 Cauchy 污染传播模型 6 87 5 125 建立合理 5 3 问题三的模型建立与求解 根据相关数据西安市在某个时间段内 PM2 5 一直高居全国前列 因此空气 污染是西安市人民关注的主要方面 采用时间序列模型对先市空气污染的各项指标进行预测 通过平稳性分析 和模型定阶的准则确定预测模型 再利用最大似然法估计各级参数 最后AIC 进行模型预报 5 3 1 时间序列模型的建立 1 平稳性检验 1 自相关函数 序列 滑动平均序列 的自相关函数 MA q 设是序列 t XXX 21 MA q qtqttttt BX 2211 其中是零均值 即平稳白噪声 且 则是零均值平稳序列 t 2 tt Var t X 且其自协方差函数和自相关函数分别为 17 0 1 0 1 11 2 22 2 2 1 2 qk qk k qkqkk q k 18 0 1 1 0 1 22 1 11 qk qk k q qkqkk k 其自相关函数呈拖尾状趋于零 偏自相关函数式非截尾但趋于零 序列 自回归序列 的自相关函数 AR p 设是序列 是平稳白t X AR p tptptt XXXX 22111t 噪声 且 则其自协方差相关函数和自相关函数为 2 tt Var 19 22111tptpttktk XXXXE 20 pkpkkk 2211 其自相关函数呈拖尾状趋于零 偏自相关函数非截尾但趋于零 而的自相关函数可以看做是序列的自相关函数和 ARMAqp MA q 序列的自相关函数的混合物 AR p 计算知本时间序列服从模型 下面考虑定阶 即估计的 ARMAqpqp 值 2 模型定阶的准则AIC 准则又称 Akaile 信息准则 AIC 准则是信息论与统计学的重要研究成AIC 果 具有重要意义 准则起源于 Kullback Lelibler 信息量 设是随机变量 它的概率AIC t X 密度是 其中含有个未知参数 设未知参数向量为 xfk T 21 k 属于分布族 其中 显然 xf xg T 21 k 0 xgxf K L 信息量可以用来刻画与的接近程度 其定义为 xg xf x xg xf xfgfId ln 则有 且有 0 gfI0 0 gfI K L 信息量是寻求最接近的参数概率密度 使得 xf xg 21 min gfI 经过理论分析 当给定样本观测值 样本容量为 设 21n xxxx n 是模型参数 未知参数个数是 的最大似然估计 设 m k T 21 k k 是其对数似然函数 信息准则是 使得 21 式中的满足 ln LAICk min2 ln 2 AIC kLk m k 设是序列 其中未知参数的个数是个 包括自 t X ARMA p q1 qpk 回归参数 滑动平均参数及 对于对数似 T 21 p T 21 q 2 然函数 22 2 22 2 ln 2 2ln 2 ln 2 1 Snn MxL nn 略去 得最小平方和估计对应似然函数 2ln 2 ln 2 1 n M 23 2 22 2 ln 2 Sn xL 由 23 式可得 的最大似然估计 2 1 2 S n 将上式代入 23 得 2 ln 2 22 nn xL 因此 序列定阶准则为 选 使得 ARMA p qAICqp 24 min 1 2lnAIC 2 qpn 其中是样本容量 与有关 若当 使得 24 达到最小n 2 qp pp qq 值 则认为序列为 ARMA q p 当序列含有未知均值时 模型为 ARMA q p tt BXB 这时 未知参数个数为 准则为 选取 使得2 qpkAICqp 25 min 2 2lnAIC 2 qpn 其式 24 与式 25 有相同的最小值点 q p 2 最大似然估计法 ML 估计 估计参数 假设是零均值正态序列 他的概率密度为 t XXX 21 ARMAqp 2 1 exp 2 1 2 1 2 2 nn T nn n tn xxxf 其中 表 T nn T qp xxxx 21121 T nnn XXE E 示参数向量是时取期望 因是序列 与有关 t X ARMAqp n X 令 则有 12 M 2 exp 2 2 T 2 22 Mxx Mxf n 对数似然函数 2 ln 2 1 2ln 2 ln 2 1 ln 2 T 2 22 Mxxn M xfxL 其中为样本观测值 21n xxxx 下证仅与有关 而与无关 事实上的第 行第列的元素是MxxT 2 ij 设的传递形式是 则 jijiij rXXEr t X ji j jt GX 0 0 2 k jikkji GGr 由此即知与无关 从而仅与有关 而与无关 12 M 2 MxxT 2 因为仅与有关 记 26 MxxT MxxS T 则称为平方和函数 进一步记 即是关于 S 21ntt xxxE t t 的条件期望 也就是已知时关于的估计 则可证 n xxx 21 n xxx 21 t n t t S 2 可见是一残差平方和形式 若使 26 式达到最小 即 S S min S 则称为最小平方和估计 S 现讨论最大似然估计 对数似然函数是 27 2 ln 2 1 2ln 2 ln 2 1 ln 2 222 Sn MxfxL 因为 均与无关 令M S 2 0 2 2 2222 SnL 解得 28 n Mxx n S T 2 将其带入 27 中得 cMxx n Mxf ln 2 ln 2 1 ln T2 其中 略去常数 得 ln1 2 ln 2 n n c c 29 ln 2 ln 2 1 TMx x n ML 的最大似然估计应使 29 达到最大 若使得 则是的最大似然估计 ML 估计 代入 m max m L m 28 即得的最大似然估计值 2 3 模型考核的检验 2 若拟合模型的残差记为 它是的估计值 例如对序列 设未知 t t AR p 参数的估计是 则残差 p 21 ntXXX ptpttt 2 1 11 记 mk n t t kn t kttk 2 1 1 2 1 Ljung Box 的检验统计量是 2 m k k kn nn 1 2 2 2 假设检验 在成立时 0 0 10 kk mkHmkH 对某些当 0 H 若充分大 近似于分布 其中 是估计的模型参数个数 n 2 2 rm r 检验法 给定显著水平 设实际算得的值是 值为 2 2 2 0 p 则当 拒绝 即认为非白噪声 模型考核不通过 2 0 2 Pp p 0 H t 需重新估计参数 而当 接受 认为是白噪声 模型考核通过 可 p 0 H t 以进行预报 4 的预报 ARMAqp 时间序列的 步预报 是根据是取值对未来时刻的随机变l 1 kk XXlk 量做出估计 估计量记作 它是的线性组合 0 lX lk lXk 1 kk XX 是零均值正态序列 则下列预报差分方程 t XXXX 210 ARMAqp 成立 2 1 21 plqlXlXlXlX kpkkk 30 只需知道 就可以递推算得 因此定义预 2 1 qXXX kkk pllXk 报变量 令 T kkk q k qXXXX 2 1 0 2 1 pj pj j j 可证下列递推预报公式 1221 1 2 1 1 1000 0010 0001 q k qqqp q q k X G G G G X p qj jpkj k q q X X G G G G 1 1 1 1 2 1 0 0 0 式中第三项当时为零 且当较小时 令初值 qp 0 k0 1 q k X 5 3 2 模型求解及结论 考虑数据的连贯性 选取 2014 4 1 日至 2014 7 31 日的空气质量指标 数据 见附录表 2 进行时间序列分析 运用 Matlab 软件作出相应时序图 AQI 进行时间序列的验证 020406080100120140 0 50 100 150 200 250 t AQI 图 6 空气质量时序图 以及相应的自相关系数图与偏自相关系数图 02468101214161820 0 2 0 0 2 0 4 0 6 0 8 Lag Sample Autocorrelation Sample Autocorrelation Function 02468101214161820 0 2 0 0 2 0 4 0 6 0 8 Lag Sample Partial Autocorrelations Sample Partial Autocorrelation Function 图 7 自相关系数图 表 8 偏自相关系数图 观察图像明显可知该序列不平稳 故作差分处理 图 9 为作一阶差分后所 得序列图 020406080100120140 100 50 0 50 100 150 200 t AQI 图 9 一阶差分后的空气质量时序图 为进一步验证平稳性 考察差分后序列自相关图和偏自相关图 02468101214161820 0 4 0 2 0 0 2 0 4 0 6 0 8 1 Lag Sample Autocorrelation Sample Autocorrelation Function 02468101214161820 0 4 0 2 0 0 2 0 4 0 6 0 8 1 Lag Sample Partial Autocorrelations Sample Partial Autocorrelation Function 图 10 一阶差分后自相关图 图 11 一阶差分后偏自相关图 观察图像可知一阶差分后的序列较平稳 故可确定模型 ARMAqp 利用准则定阶 当时 计算得到值如下表所示 AIC1 dAIC 表 5 值AIC 值AIC0 q1 q2 q3 q 0 p 1146 961128 281123 941120 14 1 p 1139 551115 621117 481118 16 2 p 1129 591117 561118 421120 42 3 p 1129 081117 291119 281119 83 由表中数据知当 时 其值最小 故取模型 又1 p1 qAIC 1 1 ARMA 由检验法 取显著性水平 运用 Matlab 软件计算得到 2 05 0 故该模型通过考核检验 可以用来预报 05 0 72 0 p 进一步计算相关参数得到 故得到模型为4891 0 1 1 1 1 1 ARMA tt BXBB 1 1 4891 0 1 其十步预报值见下表 表 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 5PM10 2 SO 2 NOCO 浓度进行预测 预测值如下表 表 7 8 月 1 10 日各项污染指标的预测值 单位 3 mmg PM2 5PM10 2 SO 2 NO CO 8 月 1 日 59 02117 978 4829 811 17 8 月 2 日 44 7999 238 0929 611 13 8 月 3 日 39 2388 467 7929 411 10 8 月 4 日 35 8882 127 5529 201 07 8 月 5 日 34 1078 267 3729 001 05 8 月 6 日 33 0875 797 15281 03 8 月 7 日 32 4274 086 9728 791 01 8 月 8 日 31 9472 806 8128 591 8 月 9 日 31 5571 766 6528 380 98 8 月 10 日 31 1970 856 4927 970 97 查阅 8 月 1 10 日西安市各项污染指标 和PM2 5PM10 2 SO 2 NO 的实际观测数据 发现与本模型预测值之间存在误差 当误差较小 故可CO 认为本模型采用模型进行预测具有实际参考价值 1 1 ARMA 5 4 污染治理建议 地球是人类长期赖以生存的家园 然而目前严峻的环境污染态势使得污染 防治工作迫在眉睫 环境的长期污染不仅危害人类的健康 也对各动植物物的 生存繁衍造成威胁 而环境的保护工作不仅要对已产生的污染进行治理 还应 从根本上预防污染的发生 因此 本文就空气污染 水体污染和噪声污染的预 防和治理手段给出如下建议 环境污染的防治工作应主要从减少污染物的产生和治理净化已产生的污染 物两方面着手 空气污染 一方面政府通过出台相关的空气污染指标排放标准 和检查制度 严格检测控制各工业污染区及汽车尾气 城镇住宅区等各项指标 的排放量和区域密度 布局合理的工业区 调整产业结构 整治燃煤企业 紧 抓 禁燃区 工作 控制污染源 此外还应调整能源结构 改革生产工艺 采 用无污染或低污染能源和工艺 如太阳能 风能 并合理利用能源 对燃料进 行预处理 改进燃烧设备 提高燃烧效率 另一方面加强废气的循环及回收利 用 使用净化装置去除烟尘和各种工业粉尘 利用物理 化学方法净化尾气中 的有害成分 引入洁净煤技术 此外还应加大环境保护的宣传力度 提高公民 素质和环保意识 水体污染 其首要任务是保护水资源 调整缺水地区的产业结构 控制高 耗水 高耗能和重污染的建设项目 其次 资源的开发利用坚持开源节流并举 的方针 完善工农业用水管理措施 建立用水考核制度 结合低资源消耗 高 效益的高新技术 推行清洁生产 此外还应建立循环用水制度和污水处理系统 尤其是有毒污染物的排放 建立资源更新补偿机制 提高排污费标准 使之高 于污染治理成本等 噪声污染 严格合理的控制机动车数量 流量及汽车喇叭的鸣笛 加强城 市道路建设 多孔隙沥青路面可降噪 施工区也应注意施工时间和噪声标准 设置部分声屏障设施 绿化均能降噪 六 模型评价与推广 6 1 模型评价 6 1 1 模型一 主成分分析模型 优点 1 本文采用降维的思想 用较少的变量反映原始变量 简化模型的同时结合了 其各项指标的综合权重 对考察城市污染程度具有实际意义 2 本文选择先主成分分析后评价污染程度的思想 通过忽略一些次要因素的影 响 考察某一综合指标的实际意义以便衡量五座城市的污染程度 具有现实 意义 缺点 1 由于所需土壤污染各指标数据的缺少 比较五座城市的污染程度时未考虑土 壤污染 可能使得主要污染因素与实际略有偏差或不足 2 降维后的新变量只是反映原始变量所提供的绝大部分信息 在实际问题的 研究中 不能以此作为衡量城市污染程度的唯一量度 6 1 2 模型二 污染传播模型 优点 1 采用定性和定量结合的方法确定并检验污染源位置 增强了模型的合理性 2 各污染指标的空间浓度分布图较直观的反映了其指标污染程度集中区 使估 计的污染源位置更接近实际位置 缺点 1 由于忽略空气的纵向扩散和风速等对污染传播的影响 使得计算结果与实际 测量值有些出入 2 由于所需水体污染和其他污染的数据缺少 只考虑了空气污染的扩散模型 对考察整体环境质量不能提供完善的数据支持 3 由于人的主观影响 北京市地理位置坐标及各项污染指标的污染源估计值与 实际存在误差 6 1 3 模型三 时间序列模型 优点 1 对具有时序性 随机性 前后时刻具有相依性 呈现某种趋势或周期性的 数据序列做了很好的处理 并做出了较准确的预测 2 生活中按时间顺序排列 随时间变化且相互关联的例子较多 如工程 经 济等 因此 本模型为其他部分领域的预测提供了参考 缺点 1 未考虑某些污染指标可能存在的季节性和周期性变化 只运用平稳时间序列 对各项污染指标进行大致预测 预测值与实际值之间可能存在误差 2 时间序列模型运用不成熟 未能对模型残差进行精确分析和异常数据的剔除 6 2 模型推广 考虑到西安市每年 11 月 15 日至次年 3 月 15 日的冬日大面积供暖设备可能 存在某些排放指标过高现象 故考虑空气重可吸入颗粒物的季节性变化 PM10 建立季节性序列模型 ARIMA 考虑时间序列 则阶差分为 其中 2 1 0 tXtd t d t d XBX 1 dd B 1 dddd B d B d B d B d 1 1 1 1 21 1 112 程为阶差分算子 则称为序列 d d t Xq d ARIMA p 当通过平稳性检验时便可建立乘积型季节性模型 t S t s BBWBB 进一步利用最大似然估计 ML 估计法 估计其相应参数 通过检验的 2 考核通过时便可进行预报 ARIMA 另外 本模型对同类型 具有时序性 随机性 前后时刻具有相依性 呈 现某种趋势或周期性的数据序列 如天文 电力 经济等其他领域的行业预测 也具有一定参考价值 参考文献 1 杨启帆等 数学建模 M 杭州 浙江大学出版社 2006 2 石辛民等 计算方法及 Matlab 实现 M 北京 清华大学出版社 2013 3 范金城等 数据分析 M 北京 科学出版社 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 计算相关矩阵和相关系数 d data f size d 2 1 m 1 while m f x1 d m a mean x1 b std x1 n length x1 1 c 1 while c n x1 c x1 c a b c c 1 end d m x1 m m 1 end x d 标准化变量 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 计算综合指标 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 76 8 73 7 66 5 61 2 62 3 101