2021年华为杯数学建模竞赛B题论文和代码

空气质量预报二次建模

尽管目前已有WRF-CMAQ模拟体系对空气质量进行预报，但由于部分污染物生成机理不完全明晰以及排放清单不确定等因素，空气质量的预报结果并不理想。因此，在WRF-CMAQ 模型一次预报的基础上进行更加准确的二次预报对提前获知大气污染并采取相应控制措施具有深远的意义。

本文针对各监测点空气质量预报数据进行量化分析，对各监测点2021年7月13日至2021 年 7 月 15 日的污染物浓度值和温度等气象条件特征值进行预测，分析时间、气象条件特征值以及污染物浓度三者之间的关系，并构建相应的数学模型。文章综合采用了

K-means++聚类算法、BP 神经网络模型、随机森林算法、皮尔逊相关性分析等方法研究污染物浓度、气象条件特征值的分析预测问题。

针对问题一，首先提取监测点A从2020年8月25日至8月28日的逐日污染物浓度实测数据，根据附录中空气质量指数(AQI)的计算公式与首要污染物的选择标准，通过MATLAB 编程得到 8 月 25 日至 8 月 28 日的 AQI 分别是 60、46、109、138，除 8 月 26日无首要污染物外，其余三天的首要污染物均为𝑶𝟑。

针对问题二，首先定性分析气象条件与污染物浓度之间的关系，由于污染物浓度变化直接受限于当时的气象条件，而与监测时间并非直接相关，所以在数据预处理方面直接删除监测点A逐小时实测数据表中与异常数据同监测时间的所有数据，不考虑时间轴可能不完整的情况，仅分析气象条件与污染物浓度之间的相关性。此外，参考AQI计算公式计算各监测时间的实时 AQI，根据空气质量等级划分表得到各监测时间的实时空气质量等级，并赋予一定的分值得到实时空气质量等级分数，将实时AQI与实时空气质量等级分数纳入到气象条件与污染物浓度的相关性分析中，用以描述污染物扩散或沉降会导致AQI下降的现象，接着采用z-score标准化方法处理数据，并运用K-means++聚类算法进行聚类分析，最终得到5类气象条件。其中，对污染物浓度影响最大的是湿度与风速，湿度、风速越大，空气质量越好，各污染物浓度越小。

针对问题三，首先对数据表中的异常数据进行近邻均值填补，因为污染物浓度与监测时间并非直接相关，而是通过气象条件进行联系，所以可以将原先的二次预测问题变成一

个两阶段预测问题。第一阶段是根据预测气象条件的时序预测问题，采用各监测点逐小时气象实测数据与对应的时间数据构建数据集，建立 BP 神经网络模型进行训练，利用训练好的BP神经网络预测各监测点2021年7月13日至2021年7月15日的气象条件特征值。第二阶段先分析各监测点逐小时污染物浓度、气象条件实测与一次预报的差异性（一次预报值取当天预报值，数据从2020年7月23日开始）,发现误差虽大，但总体走向大致相同。接着，计算各监测点污染物浓度、气象条件的一次预报误差，构建随机森林预测模型，采用输入为气象条件一次预报误差、输出为污染物浓度一次预报误差的数据集进行训练，将第一阶段预测的气象条件视作实测值与相应时间的一次预报值作差后代入到训练好的随机森林预测模型中得到各监测点从7月13日到7月15日的污染物浓度修正预测值，再与相应时间的一次预报值累加得到相应的二次预报值。运用问题1中的程序计算得到监测点A三天的AQI分别是56、62、80，首要污染物均为𝑶𝟑；监测点B三天的AQI分别是20、22、23，均没有首要污染物；监测点C三天的AQI分别是77、74、91，首要污染物均为𝑶𝟑。

针对问题四，采用与问题三大致相同的预测流程，但在训练随机森林预测模型时进行改进，将监测点A、A1、A2、A3 的污染物浓度、气象条件的一次误差合并为一个数据集进行训练，构建随机森林协同预测模型，并按照问题三中的计算流程得到监测点A三天的AQI分别是48、53、52，7月13日没有首要污染物，14日、15日首要污染物均为𝑶𝟑；监测点A1三天的 AQI分别是 56、64、59，首要污染物均为𝑶𝟑；监测点 A2三天的AQI分别是 49、53、49，7 月 14 日首要污染物为𝑶𝟑，13 日、15 日均没有首要污染物；监测点A3三天的AQI分别是46、49、52，7月15日首要污染物为𝑶𝟑，13日、14日均没有首要污染物。最终比较问题三与问题四中监测点A的预测结果，比较了A、A1、A2、A3之间的距离大小与皮尔逊相关性大小，最终得出协同预报模型不能提升污染物浓度预测准确度的结论。

        关键词：空气质量预报；K-means++聚类算法；BP 神经网络；随机森林；修正预测值；协同预报

1.2 问题提出

          由于模糊的气象场以及排放清单的不确定性，另外对包括臭氧在内的污染物生成机理的不完全明晰，导致一次预报（WRF-CMAQ）预报模型的结果并不理想，因此有了二次建模，我们依据实测数据以及一次建模预测数据，进行二次建模，具体需依次解决以下四个问题：
          问题1：附件1给出了检测点A的6类污染物的日监测浓度，我们需要按照附录所给的计算方法，计算监测点A从2020年8月25日到8月28日每天实测的AQI和首要污染物。

问题2：不同的气象条件会对AQI值产生不同的影响，我们需要对所给的每天每小时的气象条件进行聚类处理，建立气象条件与 AQI 值变化的联系，并根据气象条件对 AQI值的影响程度，对气象条件进行合理的分类，阐述清楚各类气象条件的具体特征。

问题3：根据题中所给数据，我们需要建立一个同时适用于A、B、C三个监测点（监测点两两间直线距离>100km，忽略相互影响）的二次预报数学模型，用来预测未来三天6种常规污染物单日浓度值，题目要求二次预报模型预测结果中AQI预报值的最大相对误差应尽量小，且首要污染物预测准确度尽量高。之后我们需要利用所建立的二次预报模型预测监测点A、B、C在2021年7月13日至7月15日6种常规污染物的单日浓度值，计算相应的AQI和首要污染物。 问题4：考虑到相邻区域的污染物浓度往往具有一定的相关性，我们需要建立包含A、A1、A2、A3四个监测点的协同预报模型，题目要求二次模型预测结果中AQI预报值的最大相对误差应尽量小，且首要污染物预测准确度尽量高。之后我们需要利用所建立的二次协同预报模型预测监测点A、A1、A2、A3在2021年7月13日至7月15日6种常规污染物的单日浓度值，计算相应的AQI和首要污染物。并将此模型与前面所建立的二次预报模型作比较，判断其能否提高对监测点A的污染物浓度预报的准确度。

二、模型假设

          本文模型基于以下合理假设：
          假设一：本题中所提供的数据是基本真实可靠的，本题中坏数据与缺失的数据对预测结果的影响忽略不计；
          假设二：A、B、C三个监测点两两间直线距离>100km，忽略相互之间的影响。

三、符号说明

文中使用到的符号含义如下

符号	含义
IAQIP CP BPHi, BPLo IAQIHi, IAQILo 𝑥 𝜇 𝛿 𝑥∗ 𝑥𝑙𝑖𝑙 𝑥𝑙𝑎𝑥 𝑥𝑐𝑐𝑎𝑙𝑐𝑐	污染物P的空气质量分指数，结果进位取整数 污染物P的质量浓度值 与CP相近的污染物浓度限值的高位值与低位值 与BPHi,BPLo对应的空气质量分指数 待处理原始数据 原始数据均值 原始数据标准差 z-score标准化后的数据 原始数据最小值 原始数据最大值 最小-最大标准化后的数据

四、问题一的分析与求解

          从监测点A中提取2020年8月25日到8月28日每天实测的6种污染物的监测浓度数据。

计算空气质量分指数（IAQI）

IAQIP =	IAQIHi−IAQILo	∙ (CP − BPLo) + IAQILo                                     （4-1）
	BPHi−BPLo

计算空气质量指数（AQI）

 AQI = max{IAQI1, IAQI2, IAQI3, . . . , IAQIn}                                       （4-2）

根据公式（4-1）、（4-2），利用 MATLAB 编写相应的程序，再导入所提取的日监测浓度数据，得到监测点 A 从2020 年 8 月25 日到 8月 28 日每天实测的 AQI和首要污染物，程序见附录，计算结果表4-1所示：

表4-1：实测的AQI和首要污染物

监测日期	地点	AQI计算
		AQI	首要污染物
2020/8/25	监测点A	60	O3
2020/8/26	监测点A	46	无
2020/8/27	监测点A	109	O3
2020/8/28	监测点A	138	O3

五、问题二的分析与求解

气象条件往往指不同天气现象下的水热情况，衡量特征主要包括：温度、相对湿度、

气压、风速、风向、降水量、蒸发量、冻土深度、积雪深度、日照时长等。不同的气象条件会对污染物的排放情况有着不同的影响，当气象条件利于污染物扩散或沉降时，AQI会下降，空气质量明显改善。基于这种特征，本章对附件1中监测点A逐小时污染物浓度与气象实测数据进行特性分析，并在此基础上对气象条件进行分类。

5.1问题二的分析与思路

本题要求根据气象条件对污染物浓度的影响程度，合理划分气象条件，同时还需要阐

述分类后各类气象条件的特征。换句话说，主要就是要探究不同的污染物浓度下气象条件各因素的分布情况。

          本文所提及的污染物包括SO2、NO2 、PM10、PM2.5、O3、CO。在不同的气象条件下，上述污染物会有着不同的扩散、聚集、沉降、上升等情况，所以它们的浓度实测值的变化

其实是受制于气象条件的变化，也就是说气象条件可以看作成自变量，污染物浓度可以看

作成因变量。整个问题就抽象成根据因变量划分自变量，将自变量定义域划分成若干个区间，同一区间的自变量、因变量关系相对一致，不同区间的自变量、因变量关系差异较大。

由于数据中所给的气象条件实测特征有5类，包括温度、相对湿度、气压、风速、风向，因此自变量可以看成 1×5 的向量；污染物有 6 类，所以因变量可以看成 1×6 的向量。但在实际生活中，我们的评判标准并不是污染物的浓度值，往往是根据空气质量或者天气

情况等综合评价指标来反映气象条件的分布。因此，我们希望借助不同污染物浓度所对应

的AQI值来对气象条件进行划分，得到与AQI相关的气象条件分类情况。

5.2数据处理

5.2.1预处理

在监测点A逐小时污染物浓度与气象实测数据表中，存在样本数据缺失、数据不合理、时间轴不完整的情况，而污染物的浓度变化与气象条件是直接相关的，但与时间轴的关系

是借助气象条件这一中间因素进行挂钩，其与时间轴是一个间接关系。因此，在本题探究

气象条件与污染物浓度的关系问题中，时间的相关性较低，可以删去数据表中所有存在异

常数据的样本，仅保留数据合理且完整的样本数据用于分类。

5.2.2数据标准化

          由于5类气象条件特征值的量级不同，采用z-score标准化方法将5类气象条件特征值化为同一量级下,标准化公式如下：

表5-1 五类特征平均水平

特征	平均水平
温度	25.130 ℃
湿度	68.193 %
气压	1010.607 MBar
风速	1.408 m/s
风向	155.645°

表5-2 数据处理结果

地点	温度(pu)	湿度(pu)	气压(pu)	风速(pu)	风向(pu)	AQI
监测点A	-0.299	0.504	1.505	-0.544	1.458	94
监测点A	-0.624	0.813	1.076	-0.874	-0.871	80
监测点A	-0.809	0.915	0.290	-0.874	-0.531	75
监测点A	-0.902	0.915	-0.354	0.445	-0.443	49
监测点A	-0.856	0.813	-0.640	0.445	-0.478	55
监测点A	-0.809	0.813	-0.640	1.765	-0.456	55
监测点A	-0.717	0.504	-0.426	-0.214	-0.108	57
监测点A	-0.717	0.607	0.433	-0.544	-0.929	67
监测点A	-0.809	0.813	1.291	1.765	1.530	72
监测点A	-0.948	0.813	1.577	1.435	-0.853	73
监测点A	-0.763	0.710	1.291	0.115	-0.440	63
监测点A	-0.392	0.401	1.076	-0.874	0.693	67
监测点A	-0.067	0.195	0.075	-0.544	0.764	73
监测点A	0.489	-0.525	-0.068	0.115	-0.880	64
监测点A	1.185	-1.142	-0.783	-0.214	-0.963	50
监测点A	1.788	-1.656	-1.212	0.445	-0.957	74
监测点A	1.788	-1.862	-1.284	0.445	-0.947	75
监测点A	1.556	-1.553	-1.141	-0.214	1.469	95
监测点A	1.139	-1.245	-1.141	0.445	1.518	70
监测点A	0.768	-0.833	-0.926	-2.524	1.423	78

注：表中的正、负值分别代表高于、低于平均水平

5.3模型建立与求解

5.3.1相关性分析

          针对5类标准化后的气象条件特征，进行相关性分析，计算其相关性矩阵，结果如表5-3。

表5-3 五类特征相关性分析

特征	温度	湿度	气压	风速	风向
温度	1.000	0.141	-0.828	0.076	0.146
湿度	0.141	1.000	-0.404	-0.240	0.044
气压	-0.828	-0.404	1.000	-0.030	-0.195
风速	0.076	-0.240	-0.030	1.000	-0.050
风向	0.146	0.044	-0.195	-0.050	1.000

8

从表中可以发现温度与气压有着极强的负相关，温度越高，气压越低；湿度与气压也

有着较强的负相关，但相关性程度比不上温度与气压的关系程度。这两个结论对后续的分类结果有一定的参考意义。

5.3.2空气质量指数(AQI)的处理

          空气质量等级范围根据AQI进行划分，同时对进行打分，并将相关的空气质量分数纳入到气象条件的分类依据当中，相关数据见表5-4。

表5-4 空气质量等级及对应的分数值与AQI范围

空气质量等级	优	良	轻度污染	中度污染	重度污染	严重污染
分数	100	80	60	40	20	0
AQI	[0,50]	[51,100]	[101,150]	[151,200]	[201,300]	[301,+∞)

5.3.3求解

样本的聚类特征值包括：空气质量分数、AQI、温度、湿度、气压、风速、风向，借助 matlab 采用 kmeans++算法进行聚类分析（聚类结果中的负数代表低于平均水平，正值表示高于平均水平），处理聚类结果并将其可视化，如图5-1： 类别(从左到右空气质量依次变差)

图5-1 气象条件分类结果可视化

分析上图，得到以下气象条件分类特征：

类别1：高温；湿度大；气压低；风速大；西北风、西南风、西风

类别2：低温；湿度小；气压高；风速小；东北风、东南风、东风

类别3：高温；湿度小；气压高；风速小；东北风、东南风、东风

类别4：高温；湿度小；气压低；风速小；东北风、东南风、东风

类别5：高温；湿度小；气压低；风速小；西北风、西南风、西风

在该分类结果中，温度、气压、风向与空气质量指数的相关性不强；而湿度、风速与

空气质量指数呈现了较强的相关性。湿度、风速越大，空气质量越好，这可能是因为湿度

大、风速大的气象条件有利于污染物扩散、沉淀，从而降低了大气中污染物浓度，改善了

空气质量。但总体来看，该分类的结果不是很理想，未能找到温度、气压、风向对污染物浓度影响的明确关系。

六、问题三分析与求解

6.1问题分析

本题根据附件 1 和附件 2 中的数据，建立一个同时适用于 A、B、C 三个监测点的二次预报模型，用来预测未来3天内6种常规污染物的单日浓度值。然而附件中的实测数据存在时间轴缺失以及数据错误的情况，为减少其影响，需对原始数据进行预处理，从而保证经过处理修正后污染物浓度数据能够具有良好的连续性，使得实测数据与一次预测数据时间上能够对应，方便后续的进一步处理。

6.2研究方法

本节根据已有的 A、B、C 三个监测点污染物浓度与气象一次预报数据和实测数据，首先分析了气象条件随时间的变化规律，通过分析气象条件预报值与实测值的数据，我们发现气象条件随时间有一定的变化规律，而污染物浓度的数据随时间无明显变化规律，但却随气象条件的变化有一定的规律，相关的对比图见附录。因此，将原先需要解决的二次预报问题转化为两阶段预测问题，第一阶段是气象条件的时序预测问题，我们采用神经网络模型预测未来三天的气象条件数据；第二阶段是根据气象条件一次预报误差数据与污染物浓度一次预报误差进行数据进行训练，并根据第一阶段预测出的气象条件变化值采用随机森林预测模型预测出其对应的污染物浓度变化值，然后修正污染物二次预报数据得到二次预报结果具体的预测流程如图6-1。

图6-1 空气质量二次预报流程

6.3利用BP神经网络预测气象条件

6.3.1 BP神经网络原理介绍

人工神经网络有很多强大的功能，它可以使用提供给它的相应数据来挖掘两者之间的未知规律，最后在新输入数据已知时使用挖掘出的规律作为工具来推断相应的结果。BP（Back Propagation）神经网络是按照误差逆向传播算法训练的多层前馈神经网络，是一种典型的深度学习模型，是近些年来应用最广泛的深度学习模型。其通过反向传播算法和最

速下降法不断地调整网络中的参数，以使网络的实际输出值和期望输出值的误差达到最小，其突出优点是具有很强的非线性映射能力和柔性的网络结构。

BP神经网络由一个输入层、一个或多个隐含层和一个输出层组成，每层含有若干个节点，各层节点通过加权路径来与相邻层的节点进行链接。当预测目标变量为分类变量时，输出层包含多个输出结点，但是当目标变量为区间型变量时，则只包含一个输出节点。

BP神经网络的工作方法主要分为两个过程。首先是信号的前向传播，信号从输入层输入到隐含层，通过隐含层的计算输出新的权重，最后再到输出层; 其次是误差的反向传播，BP神经网络模型根据输出层的输出结果计算预测数据与实际数据的误差，该误差用来依次调节从隐含层到输出层的权重和偏差以及输入层到隐含层的权重和偏差。如此反复进行，直到输出结果与期望的输出结果均方误差达到一个合理的范围内则结束训练。

          总而言之，BP神经网络的核心思想就是根据得到的结果计算误差，并通过反馈误差，不断修改权重和阈值，从而使输出结果的误差最小。

图6-2 BP神经网络模型

          BP神经网络的步骤主要为
          1）通过分配各层中每一个神经元的随机权重和偏差值，对前向传播过程初始化。      2）在得到每个神经元和偏差的和之后，在每个神经元中应用激活函数，将输出结果传递到下一层的神经元。

          3）若该层是输出层，则在此处应用激活函数，即可得到该模型的输出数据。           4）基于训练集，通过采用实际输出结果减去激活函数的输出来识别每个输出神经元的误差，并计算总误差。

5）对其计算偏导，将总的误差反向传递。

6）根据偏导数对权重进行更新，对更新权重后的模型重复步骤1）到步骤5）。

6.3.2建模步骤

          1）模型设计
          本文采用三层结构建立模型对天气状况进行预测，采用附件中的气象实测值作为输出、对应的时间信息作为输入进行训练，最终对7月13日至15日的气象条件进行预测。

2）数据归一化处理

首先对数据采取归一化处理，以消除各因素因量级不同而对预测精度产生影响。此处采用最小-最大标准化方法，其公式为：

𝑥𝑐𝑐𝑎𝑙𝑐𝑐 =

𝑥−𝑥𝑚𝑖𝑚 𝑥𝑚𝑎𝑥−𝑥𝑚𝑖𝑚                                                           （6-1）

处理后所有指标的数值都在0-1之间。

          3）模型训练
          在对数据归一化处理后，在同一训练样本下逐个进行模型训练，得到模型在取不同隐含层节点数时的误差。

4）得出实验结果

6.3.3数据预处理

分析附件中的数据，发现个别数据存在异常或丢失的情况，因此为了保证分析的结果合理性，需要对异常数据进行识别删除，而对缺失数据进行插值处理。实测数据中个别污染物浓度的实测值小于0，考虑实际情况可以判断其为异常数据，因此需要对其进行剔除，对于缺失数据的处理，常见的处理方法有移动平均法、直接删除法和指数平滑法等。由于实测数据按时间排布，如果不考虑缺失数据，则会对后续的处理带来影响，本文的处理方法为取相邻时间的几组数据取平均值近似代替所缺数据。

6.3.4神经网络训练

          采用MATLAB神经网络工具箱进行气象条件的时序预测，分别对监测点A、B、C的温度、湿度、气压、风速以及风向的实测值进行了训练拟合，得到如下结果。

（1）以温度实测为输出的神经网络训练结果

（a）监测点A温度训练情况

（b）监测点B温度训练情况

（c）监测点C温度训练情况

图6-3 以温度为特征的神经网络训练情况

如图所示为监测点 A、B、C 的温度实测值随时间的变化情况，其中，外围波动较大的曲线为温度的实测值，而中央的红色曲线为根据实测值拟合出的等效值，其可作为今后预报时的参考，整体的训练效果比较良好。从图中可以发现，实测数据曲线波动较大，呈不规则的毛刺状。但通过分析曲线走势可以发现，温度的实测数据大体呈周期性变化，且波动周期在一年左右，数据起始时间为四月份，根据其走向可得，在八月份左右温度达到最大，而在十二月左右温度达到最低，这与实际情况相吻合。这表明，同一地点每年同一时期的温度大体相当，即使略有误差也不会有太大的变动。温度随时间规律性的变化可为之后污染物浓度二次预报提供一定的参考。

（2）以湿度实测为输出的神经网络训练结果

代码

(1)AQI计算程序

function [AqI,k] = AQI_calculate(fname)

%% 计量标准设定

IAQI = [0,50,100,150,200,300,400,500];

SO2 = [0,50,150,475,800,1600,2100,2620];

NO2 = [0,40,80,180,280,565,750,940];

PM10 = [0,50,150,250,350,420,500,600];

PM25 = [0,35,75,115,150,250,350,500];

O3 = [0,100,160,215,265,800,800,800];

CO = [0,2,4,14,24,36,48,60];

chart = [SO2;NO2;PM10;PM25;O3;CO];

%% 读取数据

filename = fname;

sheet = 1;

xlRange = 'C2:H4';

A = xlsread(filename,sheet,xlRange);

[m,n] = size(A);

for j = 1:n

  for i = 1:m

  Chart = chart(j,:);

  B = [Chart,A(i,j)];

  B = sort(B);

  k = find(B == A(i,j));

  BP1 = Chart(k-1);

  BP2 = Chart(k);

  IA1 = IAQI(k-1);

  IA2 = IAQI(k);

  C(i,j) = ceil(((IA2-IA1)/(BP2-BP1))*(A(i,j)-BP1)+IA1);

  end         

end

for i = 1:m

  AqI(i) = max(C(i,:));

  k(i) = find(C(i,:)==AqI(i));

end

end

(2)监测点A一次预报图