尽管目前已有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，首要污染物均为𝑶𝟑。

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中提取2020年8月25日到8月28日每天实测的6种污染物的监测浓度数据。

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

计算空气质量指数（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和首要污染物

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

表5-2 数据处理结果

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

5.3模型建立与求解

5.3.1相关性分析

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

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

8

六、问题三分析与求解

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神经网络原理介绍

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

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

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

4）得出实验结果

6.3.3数据预处理

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

6.3.4神经网络训练

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

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

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

(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一次预报图