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

城市空气污染程度的分析和预测摘 要 本文讨论了有关城市污染程度、污染因素及污染扩散的问题。对于问题一，本文主要从大气污染、噪声污染和水体污染这三个面选取主要污染物，查阅北京、天津、上海、重庆和西安五座城市2007-2012年的年度平均污染数据，采用降维的思想，运用主成分分析法减少变量个数，再借助Matlab软件计算各主成分的贡献率，分析知可选取前三个主成分作为衡量污染程度的标准，最后根据综合指标得到这五个城市的污染程度从高到低依次为：重庆、上海、北京、天津、西安。通过判断相关系数的大小，确定五个城市影响人们生活的主要污染因素是水污染，其四项指标依次为化学需氧量、总氮、总磷和氨氮。对于问题二，以北京市大气污染为例。首先，利用GPS记录北京市14个城区观测点的位置，并查阅2013年污染指标、与的污染数据，绘制出相应的空间浓度分布图，估计这四种污染物的大致污染源位置依次为：附近、附近，附近和附近；其次，根据污染扩散原理和方式，建立Cauchy污染传播模型，根据各地区空气污染物的浓度分布，运用Matlab软件对数据非线性拟合，得出扩散模型各参数的值，计算得出各项污染指标的污染源位置依次为：，和；最后，比较污染物位置的计算值与实际值，发现误差相差较小，故模型建立较为合理。对于问题三，分析西安市的主要污染大气污染。收集西安市2014年4.1-7.31日的空气污染数据，根据时间序列的平稳性特点及定阶准则选取合适的时间序列模型，利用Matlab软件对序列模型的各项参数进行估计并检验模型的合理性，并将模型用于数据预报。利用时间序列模型预测西安市未来10天的空气污染状况总体等级为良。对于问题四，基于问题一、二、三对污染因素的分析和污染扩散的特点，主要从减少污染物的产生和治理净化已产生的污染物两方面，针对大气污染、水体污染和噪声污染为相关部门提供合理化防治建议。关键词 主成分分析；Cauchy污染传播模型；时间序列模型；Matlab软件30一、 问题重述环境质量的好坏对人们的生活有着最直接的影响，包括身体健康、寿命长短、生活质量以及经济发展水平、生态环境等。因此保护环境的工作需要个人、集体、国家、乃至全球各国的共同努力。事实上，环境保护事业从1962年起步至今，一直备受人们关注和重视，环境污染的防治工作也一步步开展运营，然而较经济发展的速度和水平，还有很多需要为之努力的工作要做。环境保护工作一方面是要从减少污染物的产生方面，预防和控制污染源头污染物的排放，另一方面是治理净化已产生的污染物。环境污染是指环境因受人类活动影响而改变了原有性质和状态的现象，其实质是人类活动中将大量的污染物排入环境，影响其自净能力，降低了生态系统的功能。为考察某时期某区域的环境状况，需统计该地区各项污染指标的数据，根据各污染指标的污染排放标准，确定其污染程度，并通过污染发生的原因、途径和后果制定合理的污染治理方案，加以实施运行，改善环境。统计北京、天津、上海、重庆和西安五座城市的污染数据，讨论如下问题：问题一：比较五座城市的污染程度及影响人们生活的主要污染因素；问题二：建立城市污染扩散模型确定某城市污染程度并检验；问题三：给出西安市某方面污染状况的中长期预测数据；问题四：给相关部门提供一篇治理污染的建议短文。二、 问题分析对于问题一，本文将查阅收集北京、上海、天津、重庆以及西安五座城市大气污染、噪声污染和水体污染2007-2012年的数据，由于各项污染指标众多，所以选取具有代表性的污染物：可吸入颗粒物、道路交通噪声、区域环境噪声、化学需氧量、氨氮、总氮和总磷这9项指标来综合衡量各城市的污染程度，本文将运用降维的思想，利用主成分分析法计算前几个主成分的贡献率大小及通过线性组后的综合指标，以此来衡量各城市的污染程度，根据各主成分的相关系数大小即可确定影响人们生活的主要污染因素(1)。对于问题二，本文将以北京市的大气污染为例，收集北京市2013年14个城区的、和的污染数据，绘制出这四种污染物浓度的空间浓度分布图，估算出这四种污染物的大致污染源位置；建立Cauchy污染传播模型，根据各地区空气污染物的浓度分布，运用Matlab软件对数据nlinfit非线性拟合，计算出扩散模型各参数值，得出各空气污染物的污染源位置；比较污染物位置的计算值与实际值，检验污染传播模型建立是否合理(2)。对于问题三，本文将考虑西安市大气污染的各项污染指标，搜集2014年4月1日至2014年7月31日的污染数据，根据时间序列的平稳性特点及定阶准则选取合适的时间序列模型，利用Matlab软件(3)估计序列模型的各项参数并检验模型的合理性，考核通过时方可用于数据预报。对于问题四，基于问题一、二、三对污染因素的分析和污染扩散的特点，主要从减少污染物的产生和治理净化已产生的污染物两方面，针对大气污染、水体污染和噪声污染为相关部门提供合理化防治建议。三、 模型假设1.所有城市污染数据均真实有效，具有统计价值；2.污染数据的取样为瞬时取样且在某个平面内，不考虑取样前后污染物的浓度变化；3.绘制北京大气污染浓度分布图时不考虑海拔因素的影响；4.大气污染过程经历着扩散与衰减，且沿轴方向上的扩散系数分别为常数，衰减与浓度成正比；5.不考虑日照等自然因素对各项指标污染浓度的影响。四、 符号说明4.1 模型一： 可吸入颗粒物（单位）；：二氧化硫（单位）；： 二氧化氮（单位）；： 道路交通噪声（单位分贝）；： 区域环境噪声（单位分贝）；： 化学需氧量（单位万吨）；： 氨氮（单位万吨）；： 总氮（单位万吨）；： 总磷（单位万吨） ；： 样本协方差矩阵 ；： 协方差矩阵的特征值 ；： 协方差矩阵的特征向量； ： 相关系数；： 样本方差向量 ；： 样本相关矩阵；： 样本相关矩阵的特征值；： 样本相关矩阵的特征向量；： 样本均值向量；： 主成分；： 污染程度综合指标。(其中）4.2模型二： 位置坐标/横轴；： 位置坐标/纵轴；： 位置坐标/竖轴； ： 污染物浓度；： 污染物质量； ： 污染物排放量；： 扩散系数。4.3模型三： 序列；：平稳白噪声；： 滑动平均序列；： 自回归序列；: 自回归滑动平均序列 ；： 常数；： 自协方差函数；： 自相关函数；： 概率密度；： 未知参数向量: 最大似然估计；： 自回归参数。五、 模型建立与求解5.1 问题一的模型建立与求解环境因受人类活动影响而改变了原有性质和状态的现象称为环境污染。环境污染的实质是人类活动中将大量的污染物排入环境，影响其自净能力，降低了生态系统的功能。然而生活中环境的污染也是多方面的，污染按环境要素分类，可分为：大气污染、土壤污染和水体污染。为此，本文选取环境污染的主要污染影响指标，建立污染评价模型，即主成分分析模型，比较北京、天津、上海、重庆和西安五座城市的污染程度，并得出影响人们生活的主要污染因素。5.1.1 建立主成分分析模型首先选取环境污染的部分污染指标：可吸入颗粒物、二氧化硫、二氧化氮、道路交通噪声、区域环境噪声、化学需氧量、氨氮、总氮和总磷。北京、天津、上海、重庆和西安五座城市的各项指标年平均数据（见附录表1）。问题中变量太多不但会增加计算的复杂性，也给合理的分析问题和解决问题带来很大困难；同时，这些变量之间存在的相关性，也使得这些变量所反映的信息在一定程度上有所重叠。为了减少变量个数，本文采用了降维的思想，利用主成分分析法得出污染空气的主要因素变量，提高问题研究的合理性，同时不会使数据反映的信息量有大的损失。设是的协方差矩阵，的特征值与正交化特征向量分别为及，且的第个主成分为 （1）根据平均值计算公式代入数据计算样本的均值向量为根据协方差矩阵计算公式 （2）利用Matlab软件代入数据可求得随机变量相应样本协方差矩阵为（只写下三角部分） 由于协方差矩阵对角线元素即为样本的方差值，观察矩阵易知其方差值相差较大，即各量纲间分散程度较大，不适用于环境污染模型中各成分指标间的相关性分析，为消除由于量纲不同可能带来的结果误差，故将样本变量标准化，即令 （3）得到一组新的数据。其中为的平均值，为的方差。此时，的协方差矩阵即为样本变量的相关矩阵,其中 （4） 为的协方差。则根据公式（4）运用Matlab软件代入数据即可得其样本相关矩阵为由特征值计算公式 （5）第个主成分的贡献率 （6）及前个主成分的累计贡献率 （7）其中。计算的特征值及各标准化主成分的贡献率和累计贡献率如表1所示:表1 的特征值、贡献率和累计贡献率贡献率(%)累计贡献率(%)14.742052.6952.6922.588428.7681.4531.476816.4197.8640.19282.14100500100600100700100800100900100由表1易知，第一、第二和第三主成分的累积贡献率为97.86%，并且只有第一、第二和和第三主成分对应的特征值大于1，所以选取第一主成分和第二主成分第三主成分作为标准化数据的主成分讨论问题，剩余的主成分对本文问题的讨论影响不大，因此可以忽略。由此可得的特征值所对应的前三个正交单位化特征向量为：记第个主成分为,则标准化数据的第一主成分、第二主成分、第三主成分第九主成分。由与的相关系数 （8）利用Matlab软件代入数据计算出前三个主成分与各原始变量的相关系数如表2所示 ：表2 相关系数相关系数-0.3593-0.8135-0.40650.1677-0.93010.0655-0.65380.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记样本均值向量为，则再利用方差计算公式代入数据可以得到样本向量的方差向量为：所以标准化数据的第一主成分、第二主成分和第三主成分可以用样本数据表示为：由于、和特征值的大小决定它们贡献率的大小，所以结合、和的特征值线性组合后可得衡量各地区消费支出的综合指标： （9）把、和的值代入（9）式即可得到，按照综合指标的大小便可对5个城市进行排序，从而得知各个城市的污染程度：表3 三个主成分值、综合指标值及排名排名北京-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%，大于其他主成分的贡献率，可以认为值越大，环境质量越差。因此观察与各项指标间的相关系数可知，水体污染的四项指标化学需氧量、总氮、总磷和氨氮对环境污染的影响较大，其中化学需氧量对人们生活环境造成的污染程度最大。化学需氧量是衡量水质受有机物污染程度的综合指标，当其值过高时，水体会发黑、发臭。总氮是反映水体富营养化的主要指标；过量磷可造成水体污秽异臭，湖泊发生富营养化和海湾赤潮现象；饮用水中氨氮过量在一定条件下可转化为亚硝酸盐，是强致癌物质，对人体健康极为不利；水体中氨氮过量则可导致水生物表现亢奋甚至死亡。综合考虑大气污染、水体污染和噪声污染的各项指标，得出结论：五座城市的总体污染程度依次为：重庆、上海、北京、天津、西安。重庆市由于三面环山的特殊地理位置和重工业发达、煤烟排放过量等原因，使得其总体污染程度在五座城市中居于首位；上海机动车尾气排放量常年居高不下，污染程度仅次于重庆；北京和天津地理位置接近，其工业污染物排放量较高，造成大量的雾霾扬尘天气；而西安汽车尾气和工业废气排放量较大，因此空气污染最为严重，但污染程度较其他城市略为轻微。5.2 问题二的模型建立与求解由于人类活动的影响，大气、水体及噪声等污染物浓度及其分布发生了很大变化，了解各污染物的空间分布及扩散情况对评价其污染程度有很大帮助。一方面，确定污染扩散的方式、速率和程度有助于控制其影响范围，提早预防和改善，另一方面可确定污染源的位置，从根本上杜绝或减轻污染的发生和污染程度，为相关部门提供污染治理的意见和建议。以北京市为例，建立其污染扩散模型并以空气污染为检验标准验证模型的合理性。5.2.1 估计北京市污染源位置查阅并搜集北京市污染空气的、和（可吸入颗粒）的年平均浓度含量，取包含北京市在内的长为183.4千米，宽为170.1千米的长方形区域（如图1），建立二维坐标系，得到北京市14个区的地理位置坐标(单位：千米)如表5所示:图1 北京市各城区图表4 2013年北京市各项污染指标值及地理位置坐标单位（/万吨）横坐标纵坐标朝阳区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 的空间浓度分布图 图3 的空间浓度分布图图4 的空间浓度分布图 图5 的空间浓度分布图观察图2、图3、图4和图5估计污染指标的污染源位置在附近，的污染源位置在附近，可吸入颗粒污染源位置在附近，细微颗粒的污染源位置在附近 。5.2.2 建立Cauchy扩散模型设表示时刻处某一种污染物的浓度。 在污染物传播的区域内，任意取一封闭曲面，它所围区域为，在时间内，由于扩散，经过这个区域内的污染物的质量为 （10）其中，分别代表沿轴方向的扩散系数，是方向余弦。由Gauss-Green公式得扩散过程中由于沉降等原因导致区域内质量减少量为 （11）其中表示衰减系数。 时间区域内污染物质量的变化为时间区域内污染物质量的变化也等于这段时间内增加的污染物质量减去减少的污染物的质量，即由于所选时间、区域的任意性，可以得到 (12)记污染源在点处，则在时刻有其中为污染物排放量。由以上分析可知此扩散问题满足Cauchy问题。 (13)对(13)进行傅里叶(Fourier)变换，且令，已知 ，故得常微分方程Cauchy问题唯一解 (14)对(14)求逆变换由于故得 为了简化模型，考虑空气污染的平面传播，则可得到不同空气污染指标的扩散方程: (15)如果认为经过了相当长的时间后，扩散已经终止，物质分布处于平稳状态，则方程(12)中的，于是有线性椭圆方程的边值问题也就是用傅里叶(Fourier)变换求解，其中是区域的边界. 在Cauchy问题的解(14)中，有三个未知的参数，它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根，假设点源的质量与位置已知，设观测取样为：取样时刻为，或设，是取样时间，则(12)变成：对而言，取样时间为1，而方程形状与(12)一致，把在点观测到物质密度与公式(15)都取对数，令则有 (16)令，则(16)可以写成。当时，有。利用Matlab软件对传播模型的各项参数进行估计并用于检验模型的合理性，取空气污染的取样点数据,运用非线性拟合nlinfit函数对数据进行拟合，得出函数参数即可得污染源位置。Matlab求得的空气污染指标污染源位置为，污染源位置为，可吸入颗粒的污染源位置为，的污染源位置为，与实际污染源位置很接近，故Cauchy污染传播模型建立合理。5.3 问题三的模型建立与求解根据相关数据西安市在某个时间段内PM2.5一直高居全国前列，因此空气污染是西安市人民关注的主要方面。采用时间序列模型对先市空气污染的各项指标进行预测。通过平稳性分析和模型定阶的准则确定预测模型，再利用最大似然法估计各级参数，最后进行模型预报。5.3.1时间序列模型的建立1.平稳性检验（1）自相关函数序列（滑动平均序列）的自相关函数：设是序列，。其中是零均值，即平稳白噪声，且，则是零均值平稳序列，且其自协方差函数和自相关函数分别为 (17) (18)其自相关函数呈拖尾状趋于零，偏自相关函数式非截尾但趋于零。序列（自回归序列）的自相关函数：设是序列：，是平稳白噪声，且。则其自协方差相关函数和自相关函数为 (19) (20)其自相关函数呈拖尾状趋于零，偏自相关函数非截尾但趋于零。而的自相关函数可以看做是序列的自相关函数和序列的自相关函数的混合物。计算知本时间序列服从模型，下面考虑定阶，即估计的值。（2）模型定阶的准则准则又称Akaile信息准则，AIC准则是信息论与统计学的重要研究成果，具有重要意义。准则起源于Kullback-Lelibler信息量，设是随机变量，它的概率密度是，其中含有个未知参数，设未知参数向量为。属于分布族，其中。显然。K-L信息量可以用来刻画与的接近程度，其定义为则有，且有。K-L信息量是寻求最接近的参数概率密度，使得 （21）经过理论分析，当给定样本观测值（样本容量为），设是模型参数（未知参数个数是）的最大似然估计。设是其对数似然函数，信息准则是：使得（21）式中的满足设是序列，其中未知参数的个数是个，包括自回归参数，滑动平均参数及。对于对数似然函数 （22）略去，得最小平方和估计对应似然函数 （23）由（23）式可得，的最大似然估计 将上式代入（23），得 因此，序列定阶准则为：选，使得 （24）其中是样本容量，与有关。若当，使得（24）达到最小值，则认为序列为。当序列含有未知均值时，模型为这时，未知参数个数为，准则为：选取，使得 （25）其式（24）与式(25）有相同的最小值点。2.最大似然估计法（ML估计）估计参数假设是零均值正态序列，他的概率密度为其中，,表示参数向量是时取期望（因是序列，与有关）。令，则有 对数似然函数其中为样本观测值。下证仅与有关，而与无关，事实上的第行第列的元素是,设的传递形式是，则。由此即知与无关，从而仅与有关，而与无关。因为仅与有关，记 （26）则称为平方和函数，进一步记，即是关于的条件期望，也就是已知时关于的估计，则可证可见是一残差平方和形式，若使（26）式达到最小，即，则称为最小平方和估计。现讨论最大似然估计，对数似然函数是（27）因为，均与无关，令解得 （28）将其带入（27）中得 其中，略去常数，得 （29）的最大似然估计应使（29）达到最大。若使得，则是的最大似然估计（ML估计），代入（28）即得的最大似然估计值。3.模型考核的检验 若拟合模型的残差记为，它是的估计值。例如对序列，设未知参数的估计是，则残差记 Ljung-Box的检验统计量是 假设检验。在成立时，若充分大，近似于分布，其中是估计的模型参数个数。检验法：给定显著水平，设实际算得的值是，值为，则当，拒绝，即认为非白噪声，模型考核不通过，需重新估计参数；而当，接受，认为是白噪声，模型考核通过，可以进行预报。4.的预报时间序列的步预报，是根据是取值对未来时刻的随机变量做出估计，估计量记作，它是的线性组合。是零均值正态序列，则下列预报差分方程成立 （30）只需知道，就可以递推算得，因此定义预报变量，令可证下列递推预报公式：式中第三项当时为零。且当较小时，令初值。5.3.2模型求解及结论考虑数据的连贯性，选取2014.4.1日至2014.7.31日的空气质量指标数据(见附录表2)进行时间序列分析，运用Matlab软件作出相应时序图，进行时间序列的验证：图6 空气质量时序图以及相应的自相关系数图与偏自相关系数图： 图7 自相关系数图 表8 偏自相关系数图观察图像明显可知该序列不平稳，故作差分处理。图9为作一阶差分后所得序列图。图9 一阶差分后的空气质量时序图为进一步验证平稳性，考察差分后序列自相关图和偏自相关图。图10 一阶差分后自相关图 图11 一阶差分后偏自相关图观察图像可知一阶差分后的序列较平稳，故可确定模型。利用准则定阶：当时，计算得到值如下表所示：表5 值值1146.961128.281123.941120.141139.551115.621117.481118.161129.591117.561118.421120.421129.081117.291119.281119.83由表中数据知当，时，其值最小，故取模型。又由检验法，取显著性水平,运用Matlab软件计算得到,故该模型通过考核检验，可以用来预报。进一步计算相关参数得到，故得到模型为其十步预报值见下表：表6 8月1-10日空气质量预测值与实际值比较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日各项污染指标、和的浓度进行预测，预测值如下表：表7 8月1-10日各项污染指标的预测值（单位：）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.970.97查阅8月1-10日西安市各项污染指标、和的实际观测数据，发现与本模型预测值之间存在误差，当误差较小，故可认为本模型采用模型进行预测具有实际参考价值。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日的冬日大面积供暖设备可能存在某些排放指标过高现象，故考虑空气重可吸入颗粒物的季节性变化，建立季节性序列模型：考虑时间序列，则阶差分为，其中 程为阶差分算子。则称为序列。当通过平稳性检验时便可建立乘积型季节性模型：进一步利用最大似然估计（ML估计法）估计其相应参数，通过检验的考核通过时便可进行预报。另外，本模型对同类型：具有时序性、随机性、前后时刻具有相依性、呈现某种趋势或周期性的数据序列，如天文、电力、经济等其他领域的行业预测也具有一定参考价值。参考文献1杨启帆等.数学建模M.杭州：浙江大学出版社.20062石辛民等.计算方法及Matlab实现M.北京：清华大学出版社.20133范金城等.数据分析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<fx1=d(:,m);a=mean(x1);b=std(x1);n=length(x1)+1;c=1;while c<nx1(c)=(x1(c)-a)/b;c=c+1;endd(:,m)=x1;m=m+1;endx=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.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)axis(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.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)%函数拟合确定二氧化硫污染源位置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.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,190,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 95.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 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,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 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.时间序列模型a=data(152:273,1);t=1:1:122;r11=autocorr(a) %画出数据的自相关图r12=parcorr(a) %画出数据的偏自相关图da=diff(a);r21=autocorr(da) %画出一阶差分数据的自相关图r22=parcorr(da) %画出一阶差分数据的偏自相关图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准则确定模型的阶数spec2=garchset('R',1,'M',1,'Display','off');coeffX,errorsX,LLFX=garchfit(spec2,da)sigmaForecast,w_Forecast=garchpred(coeffX,da,10)x_pred=a(end)+cumsum(w_Forecast)%确定模型参数并预测二、 数据表1 各城市各项污染指标年平均数据北京天津上海重庆西安（）0.1090.1050.0710.0900.118（）0.0290.0480.0230.0370.040（）0.0520.0420.0460.0350.042（分贝）69.267.969.367.268.2（分贝）54.054.354.754.055.3（万吨）18.6522.9424.2640.2811.83（万吨）2.052.544.745.341.38（万吨）3.263.291.555.411.18（万吨）0.440.360.180.640.12表2 西安市2014.4.1日至2014.7.31日各项指标污染数据日期AQI指数质量等级当天AQI排名PM2.5PM10CONO2SO22014/4/1108轻度污染143761491.8757262014/4/2110轻度污染155771671.8562322014/4/3103轻度污染156761371.8351272014/4/4102轻度污染140761361.8952302014/4/592良119581301.948232014/4/6108轻度污染166521641.7951272014/4/791良111521311.8848192014/4/8107轻度污染97621632.0862232014/4/9117轻度污染112791802.1959242014/4/10122轻度污染141891852.2561212014/4/1182良9059912.23792014/4/1276良9755902.1645112014/4/1390良127651222.2344122014/4/1485良80581131.6535152014/4/1574良5949901.293482014/4/1671良7947901.433482014/4/1783良97571121.6447132014/4/1888良120641181.784292014/4/1955良6934591.733282014/4/2081良125561061.5239112014/4/2162良7337741.223682014/4/2266良6838791.1437112014/4/2399良124681451.3247152014/4/24112轻度污染160831631.3740162014/4/25145轻度污染172701870.922282014/4/26112轻度污染172531681.0531122014/4/27183中度污染185793121.2642172014/4/28177中度污染188803021.5847182014/4/29142轻度污染185712331.748172014/4/30131轻度污染169662111.7650202014/5/1145轻度污染164602181.3937162014/5/2103轻度污染130471551.2433152014/5/3114轻度污染172531781.4342192014/5/4103轻度污染187491561.4532162014/5/599良184501471.3841242014/5/6107轻度污染165601621.541232014/5/7103轻度污染157531541.4441182014/5/8110轻度污染142551691.4333202014/5/9128轻度污染171701871.317132014/5/1060良9740661.341882014/5/1162良13125731.3323102014/5/1279良113401061.4739162014/5/13