2021年全国研究生数学建模竞赛华为杯B题空气质量预报二次建模求解全过程文档及程序

2021年全国研究生数学建模竞赛华为杯

B题空气质量预报二次建模

原题再现：

大气污染系指由于人类活动或自然过程引起某些物质进入大气中，呈现足够的浓度，达到了足够的时间，并因此危害了人体的舒适、健康和福利或危害了生态环境[1]。污染防治实践表明，建立空气质量预报模型，提前获知可能发生的大气污染过程并采取相应控制措施，是减少大气污染对人体健康和环境等造成的危害，提高环境空气质量的有效方法之一。
目前常用WRF-CMAQ模拟体系（以下简称WRF-CMAQ模型）对空气质量进行预报。WRF-CMAQ模型主要包括WRF和CMAQ两部分：WRF是一种中尺度数值天气预报系统，用于为CMAQ提供所需的气象场数据；CMAQ是一种三维欧拉大气化学与传输模拟系统，其根据来自WRF的气象信息及场域内的污染排放清单，基于物理和化学反应原理模拟污染物等的变化过程，继而得到具体时间点或时间段的预报结果。WRF和CMAQ的结构如图 1、图 2所示，详细介绍可以在附录提供的官网中进行查询。

请你们团队根据问题要求，基于一次预报数据及实测数据（见附件）进行空气质量预报二次数学建模，完成以下四个问题。请注意，实际工作中会遇到数据为空值或异常值的情况（见附录），故要求建立的模型具有一定的鲁棒性。

问题1. 使用附件1中的数据，按照附录中的方法计算监测点A从2020年8月25日到8月28日每天实测的AQI和首要污染物，将结果按照附录“AQI计算结果表”的格式放在正文中。
问题2. 在污染物排放情况不变的条件下，某一地区的气象条件有利于污染物扩散或沉降时，该地区的AQI会下降，反之会上升。使用附件1中的数据，根据对污染物浓度的影响程度，对气象条件进行合理分类，并阐述各类气象条件的特征。
问题3. 使用附件1、2中的数据，建立一个同时适用于A、B、C三个监测点（监测点两两间直线距离>100km，忽略相互影响）的二次预报数学模型，用来预测未来三天6种常规污染物单日浓度值，要求二次预报模型预测结果中AQI预报值的最大相对误差应尽量小，且首要污染物预测准确度尽量高。并使用该模型预测监测点A、B、C在2021年7月13日至7月15日6种常规污染物的单日浓度值，计算相应的AQI和首要污染物，将结果依照附录“污染物浓度及AQI预测结果表”的格式放在论文中。
问题4. 相邻区域的污染物浓度往往具有一定的相关性，区域协同预报可能会提升空气质量预报的准确度。如图 4，监测点A的临近区域内存在监测点A1、A2、A3，使用附件1、3中的数据，建立包含A、A1、A2、A3四个监测点的协同预报模型，要求二次模型预测结果中AQI预报值的最大相对误差应尽量小，且首要污染物预测准确度尽量高。使用该模型预测监测点A、A1、A2、A3在2021年7月13日至7月15日6种常规污染物的单日浓度值，计算相应的AQI和首要污染物，将结果依照附录“污染物浓度及AQI预测结果表”的格式放在论文中。并讨论：与问题3的模型相比，协同预报模型能否提升针对监测点A的污染物浓度预报准确度？说明原因。

整体求解过程概述(摘要)

大气污染对生态环境及人体健康有严重影响，因此准确预估未来大气污染物浓度具有十分重要的科学意义，预测未来污染物浓度情况有利于提前采取污染防治措施。目前常用的污染物预报模型 WRF-CMAQ 还存在一定缺陷，如受气象场模拟结果的制约等问题。为了提高现有预报模拟的准确性，本文首先根据六种常见大气污染物的实际观测浓度，计算空气质量指数及相应的首要污染物；然后针对不同的气象条件，分别讨论不同情况下对污染物浓度和空气质量指数的影响情况，进行气象条件的合理分类。此外，本文以 WRF-CMAQ模型模拟结果作为一次预报数据，结合实际观测数据，建立二次预报模型，预报未来三天三个监测站点 A、B、C 的常规污染物浓度、AQI 以及首要污染物；最后，考虑相邻区域污染物浓度的相关性，建立区域协同预报，结合反距离权重法建立 LSTM 模型，预测未来三天的污染物浓度、空气质量指数以及首要污染物。
针对问题一：空气质量指数常用于描述空气污染程度，通过六种常规大气污染物建立空气质量指数计算模型，求得给点监测点 A 未来四天的每日实测 AQI 和首要污染物。结果显示，这四天每日实测 AQI 先减小后增大，在 8 月 28 日达到 138，除 8 月 26 日外，其余三日的首要污染物均为 O3。
针对问题二：通过阅读相关参考文献和相关分析，发现空气质量指数及各污染物浓度主要受温度、湿度、压强和风速这些气象要素的影响。利用多元线性回归模型分别建立 AQI与各气象条件、各污染物与各气象条件之间的定量表达式。除 PM2.5 的拟合优度较小外，其他模型具有较好的解释能力。根据相关分析以及多元线性归回模型，将气象条件分为了高影响类气象条件及低影响类气象条件，湿度与风速属于高影响类气象条件，温度与压强属于低影响类气象条件。
针对问题三：WRF-CMAQ 模型输出的污染物浓度与真实值偏差较大，利用实际观测数据对一次预报模拟数据进行修正。本文将一次模拟的气象要素以及对应时间的实际观测污染物浓度数据进行数据预处理、归一化等操作，利用长短期记忆人工神经网络 LSTM 对时间序列数据集的高可预报性，建立污染物浓度的二次预报模型，结果表明该模型模拟效果良好。对三个监测站点的日 AQI 预测结果在 40-100 之间，首要污染物主要为 O3。
针对问题四：相邻区域污染物浓度具有一定的相关性，基于地理学第一定律，本文根据各监测点的地理位置进行反距离权重插值，结合问题三中的 LSTM 模型建立区域协同预报模型。结果表明模拟精度有所提升，并解决了部分在利用 lSTM 模型模拟问题三数据时出现的极大极小值问题。

模型假设：

（1）假设各监测站点位于同一海拔；
（2）假设各监测站点位于平原；
（3）假设各监测站点污染物浓度不受风向影响。

问题分析

问题一分析

本题需要我们计算监测点 A 从 2020 年 8 月 25 日到 8 月 28 日每日实测的 AQI 和首要污染物。监测点 A 逐日污染物浓度实测数据包括从 2019 年 4 月 16 日到 2021 年 7 月 12 日各污染物浓度的逐日监测浓度，共 819 条数据。由于本题只需求解四日的实测 AQI 和首要污染物，在实测数据中这四天不存在数据异常情况，因此不需要进行数据预处理操作。
针对附件 1 给出的主要污染物：二氧化硫（SO2）、二氧化氮（NO2）、粒径小于 10μm的颗粒物（PM10）、粒径小于 2.5μm的颗粒物（PM2.5）、臭氧（O3）、一氧化碳（CO），有实测的逐小时浓度数据和实测的逐日浓度数据。根据 AQI 的计算方法，我们可以直接采用已经把逐小时浓度数据平均了的 SO2 的 24 小时平均浓度、NO2 的 24 小时平均浓度、CO 的 24 小时平均浓度、PM2.5 的 24 小时平均浓度、PM10 的 24 小时平均浓度和 O3 的最大 8 小时滑动平均浓度。利用以上数据计算各污染物的空气质量分指数（IAQI），再将IAQI 的最大值作为当日的 AQI 值，且 IAQI 最大的污染物即为当日的首要污染物，若 AQI小于 50，则空气质量为优，当日无首要污染物。下面根据以上思路建立模型进行求解。

问题二分析

除受到污染物排放变化的影响，当地的气象条件也会导致 AQI 的变化。气象条件通过使污染物发生扩散或者沉降从而导致该地区 AQI 的下降，反之则会引起 AQI 的上升。
附件 1（监测点 A 空气质量预报基础数据.xlsx）提供了逐日污染物浓度实测数据、逐小时的五项气象指标的实测数据，结合问题一计算得到的逐日 AQI 实测数据，将时期统一为 2019 年 4 月 16 日到 2021 年 7 月 12 日，共 819 条数据。经过异常数据处理等工作后得到本题所需数据。
为确定不同气象条件对 AQI 的影响，首先建立相关分析模型检验各气象条件与实测AQI 之间的关系，结合显著性分析情况，建立各气象条件与 AQI 之间的多元回归模型。同时，对六种主要污染物也进行上述操作，分析各气象条件与各污染物浓度的相关性，通过多元线性回归模型建立实测 AQI 与各气象条件之间的回归方程。结合以上两部分的模型结果，对气象条件进行分类，从而分析各类气象条件的特征。

问题三分析

由于一次预报模型存在不确定性等问题，加入实测数据可以对一次预报模型进行修正从而建立二次预报模型，以提高预报精度。本题给出 A、B、C 三个监测点，在此不考虑三个监测点之间的相互影响，利用三个监测点的一次预报模型数据和实际观测数据，建立一个能同时满足三个站点的二次预报模型。通过该二次预报数学模型来预报未来三天（2021年 7 月 13 日至 2021 年 7 月 15 日）各站点 6 种污染物的单日浓度值。根据各站点的一次预报的污染物浓度数据和气象数据和各站点的实测数据，我们建立二次预报数学模型，通过模型模拟未来三天各污染物的逐小时浓度，经过平均得到单日浓度值，再根据 AQI 计算公式计算得到监测点 A、B、C 未来三日的 AQI 及首要污染物。
需要注意的是，根据题目要求，模型结果应使单日 AQI 预报值尽可能小而首要污染物预测浓度尽可能准确。为此我们利用附件材料中提供的实测污染物浓度数据与一次预报的气象要素数据进行模拟，建立一种基于长短期记忆人工神经网络的预测模型（LSTM），预测未来逐小时的各站点 6 种污染物浓度值，再计算得到每日 AQI 及首要污染物。

问题四分析

在问题三中，我们没有考虑到 A、B、C 三个监测站点之间的相互影响。而在实际情况下，地理事物之间是存在空间关联性的，根据地理学第一定律可知，任何事物之间都是相互关联的，距离较近的事物之间的联系更紧密。因此本题中，在对空气质量预报进行二次建模时，还需要考虑各站点之间的相互影响，建立区域协同预报模型，可能会提高预报结果的准确性。
本题给出了监测点 A 附近的三个监测点 A1、A2、A3 的一次预报气象数据与污染物浓度数据，实测气象数据与污染物浓度数据。在此需考虑 A、A1、A2、A3 之间的相互影响，建立包含四个监测站点的协同预报模型，因此我们基于反距离权重法对各站点的预报数据重新补充，再利用预报数据和实测数据建立适用于该站点二次预报模型。
通过建立的模型，我们需要对四个监测点未来三天的常规污染物单日浓度进行预测，同时计算单日 AQI 以及首要污染物。将该模型得到的预报结果与问题 3 中的没有考虑相互影响的二次预报模型的预报结果进行对比，根据对比结果，研究分析协同预报模型是否能起到提高未来污染物浓度预报精度的作用。

模型的建立与求解

在计算空气质量分指数（IAQI）时，要对使用的数据进行异常值的剔除。例如当 PM2.5浓度普遍小于 200μg∕m3 时，出现一个大于 400 的值，则该值为异常值，需要将它去除在外不进行 IAQI 的计算。
使用空气质量分指数（IAQI）计算得到空气质量指数（AQI）后，再根据计算得到的空气质量指数（AQI）划分空气质量等级，具体的等级及对应空气质量指数（AQI）范围如下表 2 所示。

由表 2 可知，当 AQI 小于或等于 50 时，空气质量等级为优，当天无首要污染物；而当 AQI 大于 50 时，IAQI 最大的污染物项目为首要污染物。IAQI 最大的污染物为首要污染物。若 IAQI 最大的污染物为两项或两项以上时，并列为首要污染物。IAQI 大于 100 的污染物为超标污染物。

根据由上述 AQI 及首要污染物的定义，在 MATLAB 软件中建立计算模型。计算得到2019 年 4 月 16 日到 2021 年 7 月 12 日监测点 A 每天实测的 AQI 和首要污染物。根据问题一要求，提取出 2020 年 8 月 25 日到 8 月 28 日的数据，结果如下表 3 所示：

根据计算结果可知，除 2020 年 8 月 26 日 AQI 小于 50，无首要污染物外，其余三天AQI 均大于 50，且首要污染物均为 O3。其中 2020 年 8 月 25 日空气质量等级为良，2020年 8 月 27 日空气质量等级为轻度污染，2020 年 8 月 28 日空气质量也为轻度污染。

对实测逐日气象数据、AQI、各气象要素建立的多元回归方程模型结果如下表 5 所示。根据模型参数，如 P 值、标准估算误差等，各模型均通过了显著性检验。

从上表中可以看出，AQI 对应的 R^2 值大约为 0.445，拟合效果较好，可以解释一部分气象要素变化与 AQI 的变化的关系，并且 P 值约等于 0，说明接受原假设的概率很低，推翻了原假设，有总体显著性差异，可以接受各气象要素与对应参数组成的回归曲线与 AQI存在线性关系。

计算结果表明，压强对应的回归系数的 P 值大于 0.05，说明该值没有通过显著性检验，即气压对于空气质量指数的改变可能没有明显的线性关系，所以气压对于空气质量指数的变化没有明确的作用。
而剩下的温度、湿度、风速的都通过显著性检验，表明温度湿度和风速对 AQI 有明显的线性关系，对 AQI 的影响大小排序为湿度＞风速＞温度，其中，温度与AQI 呈正相关，但影响较小；风速和湿度与 AQI 呈负相关关系，对 AQI 的影响也较大。
对污染物 SO2 的浓度而言，对应的 R^2 为 0.434，这些气象条件可以解释近一半的 SO2浓度。各自变量项的 P 值都小于 0.01，结果显著，具有明显的线性关系。对比湿度与风速的回归系数，温度和压强的回归系数都很小，表明于 SO2 而言，湿度和风速的影响较大，且都呈负相关关系，湿度对 SO2 浓度的影响最大；而温度和压强对 SO2 的影响较小。
对污染物 NO2 来说，对应 R^2 为 0.468，气象条件对它的解释程度略高于对 AQI 值的解释程度。所有气象条件都通过了显著性检验，且所有气象条件都与 NO2 呈负相关关系，其中，各参数影响大小排序为风速＞湿度＞温度＞气压。
对污染物 PM10 而言，它的对应 R^2 为 0.53，有超过一半的污染物浓度可以通过这四个气象参数来解释，是这些组里解释程度最高的，整体显著性水平较高。所有气象参数中，温度没有通过显著性检验，在此不做讨论。气压、湿度与风速对 PM10 浓度的影响大小为风速＞湿度＞气压，风速与湿度的回归系数大致相同，属于对 PM10 影响较大的条件，且都为负影响；气压对 PM10 浓度影响较小，为正影响。
PM2.5 的 R^2 值只有 0.245，接近四分之一的该值可以被这些气象条件解释到，说明还有很多其他的因素会影响到 PM2.5 的变化。整体显著性水平高，但温度和气压的没有通过显著性检验，从回归系数上看也对 PM2.5 影响较小。湿度和气压对 PM2.5 的影响较大，且都是负影响，风速对 PM2.5 浓度的影响最大。
O3 的回归方程和其他物种较为不同。它的 R^2 为 0.478，解释程度大致处于平均水平上。所有气象参数均通过了显著性检验。对 O3 的影响大小排列为湿度＞温度＞风速＞气压。湿度和温度的影响程度基本一致，但温度对 O3 的影响为明显的正影响，而湿度和风速是负影响。气压也呈正影响，但影响较小。
对 CO 而言，R^2 为 0.328，约三分之一的该浓度值可以被气象条件解释。气压项没有通过显著性检验。温度与风速对 CO 浓度的变化有较大影响，湿度对 CO 浓度也略有影响，都呈负相关关系。
综上所述，在上面各气象要素对于六种污染物的影响的统计下，发现：
无论是六种污染物中的哪一种，湿度与风速都对于污染物浓度变化有明显的负相关关系，并且在一些污染物中，只有这两个气象要素通过了显著性检验。说明湿度变化与风速变化对于污染物浓度变化影响极大，属于高影响类气象条件，即湿度与风速增大，污染物浓度极有可能下降，而对应两因素减少时，污染物浓度极有可能上升。
而对于气压和温度来说，在对于六种污染物的影响的统计下，观察到，两个气象要素各有三次在总体显著性较高的情况下，而该要素的显著性水平没有通过显著性检验，并且通过显著性检验的情况下，对于一些污染物的浓度的变化存在负影响，而另一些则为正影响，说明气压变化与温度变化对于污染物浓度变化影响较小，属于低影响类气象条件。
并且发现高影响类气象要素具有对所有污染物变化有一致的影响作用，而低影响类污染物对特定的污染物的浓度变化可能有不同的影响作用，可能为正也可能为负，并且有些情况可能都无法通过显著性检验，而无法判定对于污染物的浓度变化有影响。

问题 3 代码

1. # -*- coding: utf-8 -*-
2. import torch
3. import numpy as np
4. import pandas as pd
5.
6. from torch.utils.data import Dataset
7. from sktime.transformations.series.impute import Imputer
8. from sklearn.model_selection import train_test_split
9.
10.# 数据预处理
11.path = './resources/附件 1 监测点 A 空气质量预报基础数据.xlsx'
12.na_values = ['NAN', 'NA', '—']
13.df_sheet = pd.read_excel(path, sheet_name='监测点 A 逐日污染物浓度实测数据',
na_values=na_values)
14.
15.
16.def preprocess_data(dataframe, interp='drift', interval=8):
17. y = dataframe.copy()
18.
19. date = [y['监测日期'].dt.month,
20. y['监测日期'].dt.isocalendar().week,
21. y['监测日期'].dt.day,
22. y['监测日期'].dt.quarter]
23.
24. del y['监测日期']
25. del y['地点']
26.
27. if interp:
28. transformer = Imputer(method=interp)
29. y = transformer.fit_transform(y, y)
30.
31. # 切分时间
32. y['月'] = date[0]
33. y['周'] = date[1]
34. y['日'] = date[2]
35. y['季节'] = date[3]
36.
37. mean, std = y.mean(axis=0), y.std(axis=0)
38. y = y.apply(lambda x: (x - x.mean()) / x.std(), axis=0)
39.
40. y, y_ = y[:-3], y[-3:] # predict last three day, skip it
41.
42. return [y[i:i + interval] for i in range(0, len(y) - interval + 1)],
y_, mean.to_numpy(), std.to_numpy()
43.
44.
45.dataset, time_encodding, mean, std = preprocess_data(df_sheet)
46.train_data, test_data = train_test_split(dataset, test_size=0.2,
random_state=42)
47.
48.
49.# 组件
50.class MovingAverageLogger:
51.
52. def __init__(self, beta=0.98):
53. self.b = beta
54.
55. self.v = None
56.
57. def update(self, v):
58.
59. if self.v is None:
60. self.v = v
61. else:
62. self.v = self.v * self.b + v * (1. - self.b)
63.
64. def value(self):
65. return self.v
66.
67.
68.class TimeSeriesDataset(Dataset):
69.
70. def __init__(self, data):
71. super(TimeSeriesDataset, self).__init__()
72. self.data = data
73.
74. def __len__(self):
75. return len(self.data)
76.
77. def __getitem__(self, x):
78. X = self.data[x].to_numpy().astype(np.float32)[:-1]
79. Y = self.data[x].to_numpy().astype(np.float32)[-1]
80.
81. return torch.from_numpy(X), torch.Tensor(Y[:6])
82.
83.
84.class Baseline(torch.nn.Module):
85.
86. def __init__(self, in_dim, out_dim, head_dim=64, dropout=0.2):
87. super(Baseline, self).__init__()
88.
89. self.bilstm = torch.nn.LSTM(in_dim,
90. head_dim,
91. 2,
92. batch_first=True,
93. dropout=dropout,
94. bidirectional=True)
95.
96. self.mlp = torch.nn.Sequential(
97. torch.nn.Dropout(dropout),
98. torch.nn.Linear(head_dim * 2 * 7, head_dim),
99. torch.nn.ReLU(),
100. torch.nn.Dropout(dropout),
101. torch.nn.Linear(head_dim, out_dim)
102. )
103.
104. def forward(self, x):
105. return self.mlp(self.bilstm(x)[0].reshape(len(x), -1))
106.
107. def forecasting(self, x, z, mean, std):
108. x = x.to_numpy().astype(np.float32)
109. z = z.to_numpy().astype(np.float32)[:, 6:]
110.
111. ys = []
112. for i in range(len(z)):
113. X = torch.Tensor([x])
114. Y = self(X)
115. ys.append(Y)
116.
117. y = Y.detach().numpy()[0]
118. x = np.concatenate([x[1:], np.concatenate([y,
z[i]])[np.newaxis, :]], axis=0)
119.
120. return torch.cat(ys, dim=0).detach().numpy() * std[:6] + mean[:6]
121.
122.
123.from torch.utils.data import DataLoader
124.
125.train_dataset = DataLoader(TimeSeriesDataset(train_data),
126. batch_size=16,
127. shuffle=True, drop_last=True)
128.
129.test_dataset = DataLoader(TimeSeriesDataset(test_data),
130. batch_size=32,
131. shuffle=False)
132.
133.lr = 3e-4
134.model = Baseline(10, 6)
135.optimizer = torch.optim.AdamW(model.parameters(), lr=lr)
136.criterion = torch.nn.MSELoss()
137.
138.if __name__ == '__main__':
139.
140. from tqdm import tqdm
141.
142. for i in range(200):
143.
144. model.train()
145. logger = MovingAverageLogger()
146. with tqdm(total=len(train_dataset), desc=f'Epoch {i}:
Training ...') as t:
147. for X, y in train_dataset:
148. optimizer.zero_grad()
149. Y = model(X)
150. loss = criterion(Y, y)
151. loss.backward()
152. optimizer.step()
153.
154. logger.update(loss.item())
155. t.update()
156. t.set_postfix(loss=logger.value())
157.
158. model.eval()
159. logger = MovingAverageLogger()
160. mse = 0
161. count = 0
162. with tqdm(total=len(test_dataset), desc=f'Epoch{i}: Testing ...')
as t:
163. for X, y in test_dataset:
164. Y = model(X)
165. loss = criterion(Y, y)
166.
167. count += len(Y)
168. mse += loss.item() * count
169.
170. logger.update(loss.item())
171. t.update()
172. t.set_postfix(loss=logger.value())
173. print(mse / count)
174.
175. print(model.forecasting(dataset[-1][1:], time_encodding, mean,
std))

论文缩略图：

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码：

问题 1 代码

1. import pandas as pd
2. import matplotlib.pyplot as plt
3. import numpy as np
4. import math
5. %matplotlib inline
6.
7. import warnings
8. warnings.filterwarnings("ignore")
9.
10.plt.rcParams['font.sans-serif'] = ['SimHei'] # 步骤一（替换 sans-serif 字体）
11.plt.rcParams['axes.unicode_minus'] = False # 步骤二（解决坐标轴负数的负号显
示问题）
12.
13.O3_df=pd.read_excel("../附件 1 监测点 A 空气质量预报基础数
据.xlsx",sheet_name=1,index_col=0,usecols=[0,6])
14.O3_df_range=O3_df.loc["2020-8-25":"2020-8-28"]
15.
16.O3_rs_dict={}
17.for i in range(4):
18. print(f"第{25+i}天的数据情况如下：")
19. a=O3_eight_hours.iloc[i*24:(i+1)*24,:]
20. print(a)
21. max_value=a.max()
22. print(max_value)
23. day=25+i
24. key="2020-8-"+str(day)
25. O3_rs_dict[key]=max_value
26.
27.
28.day_ground=pd.read_excel("../附件 1 监测点 A 空气质量预报基础数
据.xlsx",sheet_name=2,index_col=0)
29.day_ground.drop(["地点"],axis=1,inplace=True)
30.day_ground.columns=[i.split("(")[0] for i in list(day_ground.columns)]
31.
32.q1_data=day_ground["2020-8-25":"2020-8-28"]
33.q1_data.iloc[:,4]=[112.25,92.125,169.125,200.75]
34.IAQI_table={
35. "IAQI":[0,50,100,150,200,300,400,500],
36. "CO":[0,2,4,14,24,36,48,60],
37. "SO2":[0,50,150,475,800,1600,2100,2620],
38. "NO2":[0,40,80,180,280,565,750,940],
39. "O3":[0,100,160,215,265,800,np.nan,np.nan],
40. "PM10":[0,50,150,250,350,420,500,600],
41. "PM2.5":[0,35,75,115,150,250,350,500]
42. }
43.IAQI_DF=pd.DataFrame(data=IAQI_table,index=None)
44.IAQI=IAQI_DF.iloc[:,0].tolist()
45.
46.def calculate_AQI(C_p):
47. """
48. 每次计算一种指标的四天的值，返回一个字典，最后放一起比一下。每次传入 q1_data
的列名，算一个的
49. """
50. C_p_value=q1_data.loc[:,C_p].tolist()
51. C_p_range=IAQI_DF.loc[:,C_p].tolist()
52.
53. C_p_dict={}
54. C_p_list=[]
55.
56. for num in C_p_value:
57. for i in range(7):
58. if (num > C_p_range[i]) and (num < C_p_range[i+1]):
59. Lo=i # 低位值和高位值序号只差 1，只记录一个就可以
60. break
61. Lo_value=C_p_range[Lo]
62. Hi_value=C_p_range[(Lo+1)]
63. IAQI_Lo_value=IAQI[Lo]
64. IAQI_Hi_value=IAQI[(Lo+1)]
65.
66. IAQI_p=math.ceil((IAQI_Hi_value-IAQI_Lo_value)/(Hi_valueLo_value)*(num-Lo_value) + IAQI_Lo_value)
67.
68. C_p_list.append({"IAQI":IAQI_p,"AQI_lo":IAQI_Lo_value})
69. return C_p_list
70.
71.rs_dict={}
72.for i in ["SO2","NO2","PM10","PM2.5","O3","CO"]:
73. a=calculate_AQI(i)
74. print(f"{i}：{a}")
75. a_value=list(i['IAQI'] for i in a)
76. rs_dict[i]=a_value
77. print(f"\n 对于{i}来说：单纯 IAQI 值：{a_value}\n")
78.
79.
80.rs_df=pd.DataFrame(data=rs_dict)
81.rs_df.max(axis=1)

问题 2 代码

1. import pandas as pd
2. import matplotlib.pyplot as plt
3. import numpy as np
4. import seaborn as sn
5. %matplotlib inline
6.
7. import warnings
8. warnings.filterwarnings("ignore")
9.
10.plt.rcParams['font.sans-serif'] = ['SimHei'] # 步骤一（替换 sans-serif 字体）
11.plt.rcParams['axes.unicode_minus'] = False # 步骤二（解决坐标轴负数的负号显
示问题）
12.
13.day_ground=pd.read_excel("../附件 1 监测点 A 空气质量预报基础数
据.xlsx",sheet_name=2,index_col=0)
14.day_ground.drop(["地点"],axis=1,inplace=True)
15.day_ground.columns=[i.split("(")[0] for i in list(day_ground.columns)]
16.
17.plt.figure(figsize=(15,15))
18.ax1=plt.subplot(3,1,1)
19.ax1.set_title("原数据")
20.plt.plot(day_ground['PM10'][:"2020-04-30"], color="c")
21.
22.ax2=plt.subplot(3,1,2)
23.ax2.set_title("填充插值")
24.plt.plot(day_ground["PM10"][:"2020-04-
30"].interpolate(method='pad', limit=2),color="r")
25.
26.ax3=plt.subplot(3,1,3)
27.ax3.set_title("线性插值")
28.plt.plot(day_ground["PM10"][:"2020-04-
30"].interpolate(method='linear', limit=2),color="b")
29.
30.plt.savefig("插值 1.png", dpi=150)
31.
32.plt.figure(figsize=(20,10))
33.plt.xticks(fontsize=20)
34.plt.yticks(fontsize=20)
35.plt.plot(day_ground['PM10']["2019-05-31":"2019-12-
15"], color="c",linewidth=4)
36.plt.plot(day_ground["PM10"]["2019-05-31":"2019-12-
15"].interpolate(method='linear', limit=2),'o',color='lightcoral',linewidt
h=3)
37.plt.plot(day_ground["PM10"]["2019-05-31":"2019-12-
15"].interpolate(method='pad', limit=2),'--',color='steelblue')
38.plt.legend(['row data', 'linear', 'pad'], loc='best')
39.
40.plt.savefig("插值 3.png", dpi=150)
41.
42.# 未插值的图像
43.plt.figure(figsize=(15,15))
44.ax1=plt.subplot(3,1,1)
45.ax1.set_title("原数据")
46.plt.plot(day_ground['PM10'][:"2020-04-30"], color="c")
47.
48.ax2=plt.subplot(3,1,2)
49.ax2.set_title("时间插值")
50.plt.plot(day_ground["PM10"][:"2020-04-
30"].interpolate(method='time', limit=2),color="r")
51.
52.ax3=plt.subplot(3,1,3)
53.ax3.set_title("最近邻插值")
54.plt.plot(day_ground["PM10"][:"2020-04-
30"].interpolate(method='nearest', limit=2),color="b")
55.
56.plt.savefig("插值 2.png", dpi=150)
57.
58.corr_day = day_ground.corr()
59.corr_day
60.
61.mask = np.array(corr_day)
62.mask[np.tril_indices_from(mask)] = False
63.# 绘制图像
64.
65.fig,ax = plt.subplots()
66.fig.set_size_inches(20,10)
67.# Spectral bwr twilight Pastel1
68.plt.xticks(fontsize=14)
69.plt.yticks(fontsize=14)
70.
71.sn.heatmap(corr_day,mask=mask,vmax=0.9,square=True,annot=True)
72.
73.plt.yticks(rotation=360)
74.plt.savefig("相关 1.png", dpi=150)
75.
76.weather_air=pd.read_excel("../附件 1 监测点 A 空气质量预报基础数
据.xlsx",sheet_name=1,index_col=0,na_values=['NAN', 'NA', '—'])
77.
78.weather_air.drop("地点",axis=1,inplace=True)
79.weather_air.columns=[i.split("
（")[0].split("(")[0] for i in list(weather_air.columns)]
80.weather_air.fillna(value=0,inplace=True)
81.
82._, ax = plt.subplots(figsize=(12, 10))
83.corr = weather_air.corr(method='pearson')
84.# corr = df.corr(method='kendall')
85.# corr = df.corr(method='spearman')
86.cmap = sns.diverging_palette(220, 10, as_cmap=True) # 在两种 HUSL 颜色之间制
作不同的调色板。图的正负色彩范围为 220、10，结果为真则返回 matplotlib 的 colormap
对象
87._ = sns.heatmap(
88. corr, # 使用 Pandas DataFrame 数据，索引/列信息用于标记列和行
89. cmap=cmap, # 数据值到颜色空间的映射
90. square=True, # 每个单元格都是正方形
91. cbar_kws={'shrink': .9}, # `fig.colorbar`的关键字参数
92. ax=ax, # 绘制图的轴
93. annot=True, # 在单元格中标注数据值
94. annot_kws={'fontsize': 12}) # 热图，将矩形数据绘制为颜色编码矩阵
95.
96.plt.xticks(rotation=30)
97.plt.savefig("相关 2.png", dpi=150)
98.plt.show()
99.
100.once_predict=pd.read_excel("../附件 1 监测点 A 空气质量预报基础数
据.xlsx",sheet_name=0,index_col=0,na_values=['NAN', 'NA', '—'])
101.once_predict_air=once_predict.iloc[:,2:]
102.once_predict_air.columns=[i.split("
（")[0].split("(")[0] for i in list(once_predict_air.columns)]
103.once_predict_air.head()
104.once_predict_air.fillna(value=0,inplace=True)
105.
106.plt.figure(figsize=(14,12))
107.corr=once_predict_air.corr()
108.mask = np.array(corr)
109.mask[np.tril_indices_from(mask)] = False
110.
111.sns.heatmap(corr,mask=mask,linewidths=0.1,vmax=1.0,square=True,linecolor=
'white',annot=True,cmap="Spectral")
112.plt.xticks(rotation=30)
113.
114.plt.savefig("相关 3.png", dpi=150)
115.
116.once_predict_air[["O3 小时平均浓度","地表温度"]]["2020-09-01":"2020-12-31"]
117.
118.Week_data=once_predict_air.resample("W").mean()
119.
120.fig, ax1 = plt.subplots(figsize=(16,8))
121.plt.xticks(rotation=45,fontsize=14)
122.plt.yticks(fontsize=14)
123.
124.ax1.bar(Week_data.index, Week_data["O3 小时平均浓度
"], color="blue",label="O3 小时平均浓度",alpha=0.5, width=3)
125.ax1.set_xlabel("月份")
126.ax1.set_ylabel("O3 小时平均浓度")
127.
128.ax2 = ax1.twinx()
129.ax2.plot(Week_data.index, Week_data["地表温度"], color="red", label="地表温
度")
130.ax2.set_ylabel("地表温度")
131.
132.fig.legend(loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax1.t
ransAxes)
133.plt.savefig("温度和 O3.png", dpi=150)
134.plt.show()
135.
136.Week_data_once=weather_air.resample("d").mean()["2020-1-1":"2020-3-1"]
137.
138.fig, ax1 = plt.subplots(figsize=(16,8))
139.plt.xticks(rotation=45,fontsize=14)
140.plt.yticks(fontsize=14)
141.
142.
143.ax1.plot(Week_data_once.index.date,Week_data_once["SO2 监测浓度"],"+--
",label="SO2 监测浓度",color="cornflowerblue",alpha=0.5,linewidth=3)
144.
145.ax1.set_xlabel("月份")
146.ax1.set_ylabel("O3 小时平均浓度")
147.
148.
149.ax2 = ax1.twinx()
150.ax2.plot(Week_data_once.index.date, Week_data_once["温度"], "o--",label="
温度",color="c")
151.ax2.plot(Week_data_once.index.date, Week_data_once["湿度"], ">--",label="
湿度",color="b")
152.ax2.plot(Week_data_once.index.date, Week_data_once["风向"], "*--
",color="darksalmon",label="风向")
153.
154.
155.ax2.set_ylabel("气象条件")
156.
157.fig.legend(loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax1.t
ransAxes)
158.plt.savefig("SO2 和其他.png", dpi=150)
159.plt.show()
160.
161.Week_data_once=weather_air.resample("d").mean()["2020-1-1":"2020-3-1"]
162.
163.fig, ax1 = plt.subplots(figsize=(16,8))
164.plt.xticks(rotation=45,fontsize=14)
165.plt.yticks(fontsize=14)
166.
167.ax1.plot(Week_data_once.index.date,Week_data_once["NO2 监测浓度
"],"+:",label="NO2 监测浓度",color="cornflowerblue",alpha=0.5,linewidth=3)
168.ax1.plot(Week_data_once.index.date,Week_data_once["PM10 监测浓度
"],"+-.",label="PM10 监测浓度",color="#fd4659",alpha=0.5,linewidth=3)
169.ax1.plot(Week_data_once.index.date,Week_data_once["PM2.5 监测浓度"],"+-
",label="PM2.5 监测浓度",color="#fdaa48",alpha=0.5,linewidth=3)
170.ax1.set_xlabel("月份")
171.ax1.set_ylabel("O3 小时平均浓度")
172.
173.
174.ax2 = ax1.twinx()
175.ax2.plot(Week_data_once.index.date, Week_data_once["温度"], "o--",label="
温度",color="c")
176.ax2.plot(Week_data_once.index.date, Week_data_once["湿度"], ">--",label="
湿度",color="b")
177.ax2.plot(Week_data_once.index.date, Week_data_once["风向"], "*--
",color="darksalmon",label="风向")
178.
179.ax2.set_ylabel("气象条件")
180.
181.fig.legend(loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax1.t
ransAxes)
182.plt.savefig("NO2PM10 和其他.png", dpi=150)
183.plt.show()

问题 4 代码

1. import tensorflow as tf
2. def placeholder(P, Q, N):
3. X = tf.compat.v1.placeholder(shape=(None, P, N), dtype=tf.float32)
4. TE = tf.compat.v1.placeholder(shape=(None, P + Q, 2), dtype=tf.int32)
5. label = tf.compat.v1.placeholder(shape=(None, Q, N), dtype=tf.float32)
6. is_training = tf.compat.v1.placeholder(shape=(), dtype=tf.bool)
7. return X, TE, label, is_training
8.
9. def FC(x, units, activations, bn, bn_decay, is_training, use_bias=True):
10. if isinstance(units, int):
11. units = [units]
12. activations = [activations]
13. elif isinstance(units, tuple):
14. units = list(units)
15. activations = list(activations)
16. assert type(units) == list
17. for num_unit, activation in zip(units, activations):
18. x = tf_utils.conv2d(
19. x, output_dims=num_unit, kernel_size=[1, 1], stride=[1, 1],
20. padding='VALID', use_bias=use_bias, activation=activation,
21. bn=bn, bn_decay=bn_decay, is_training=is_training)
22. return x
23.
24.def STEmbedding(SE, TE, T, D, bn, bn_decay, is_training):
25. '''
26. spatio-temporal embedding
27. SE: [N, D]
28. TE: [batch_size, P + Q, 2] (dayofweek, timeofday)
29. T: num of time steps in one day
30. D: output dims
31. retrun: [batch_size, P + Q, N, D]
32. '''
33. # spatial embedding
34. SE = tf.expand_dims(tf.expand_dims(SE, axis=0), axis=0)
35. SE = FC(
36. SE, units=[D, D], activations=[tf.nn.relu, None],
37. bn=bn, bn_decay=bn_decay, is_training=is_training)
38. # temporal embedding
39. dayofweek = tf.one_hot(TE[..., 0], depth=7)
40. timeofday = tf.one_hot(TE[..., 1], depth=T)
41. TE = tf.concat((dayofweek, timeofday), axis=-1)
42. TE = tf.expand_dims(TE, axis=2)
43. TE = FC(
44. TE, units=[D, D], activations=[tf.nn.relu, None],
45. bn=bn, bn_decay=bn_decay, is_training=is_training)
46. return tf.add(SE, TE)
47.
48.def spatialAttention(X, STE, K, d, bn, bn_decay, is_training):
49. '''
50. spatial attention mechanism
51. X: [batch_size, num_step, N, D]
52. STE: [batch_size, num_step, N, D]
53. K: number of attention heads
54. d: dimension of each attention outputs
55. return: [batch_size, num_step, N, D]
56. '''
57. D = K * d
58. X = tf.concat((X, STE), axis=-1)
59. # [batch_size, num_step, N, K * d]
60. query = FC(
61. X, units=D, activations=tf.nn.relu,
62. bn=bn, bn_decay=bn_decay, is_training=is_training)
63. key = FC(
64. X, units=D, activations=tf.nn.relu,
65. bn=bn, bn_decay=bn_decay, is_training=is_training)
66. value = FC(
67. X, units=D, activations=tf.nn.relu,
68. bn=bn, bn_decay=bn_decay, is_training=is_training)
69. # [K * batch_size, num_step, N, d]
70. query = tf.concat(tf.split(query, K, axis=-1), axis=0)
71. key = tf.concat(tf.split(key, K, axis=-1), axis=0)
72. value = tf.concat(tf.split(value, K, axis=-1), axis=0)
73. # [K * batch_size, num_step, N, N]
74. attention = tf.matmul(query, key, transpose_b=True)
75. attention /= (d ** 0.5)
76. attention = tf.nn.softmax(attention, axis=-1)
77. # [batch_size, num_step, N, D]
78. X = tf.matmul(attention, value)
79. X = tf.concat(tf.split(X, K, axis=0), axis=-1)
80. X = FC(
81. X, units=[D, D], activations=[tf.nn.relu, None],
82. bn=bn, bn_decay=bn_decay, is_training=is_training)
83. return X
84.
85.def temporalAttention(X, STE, K, d, bn, bn_decay, is_training, mask=True):
86. '''
87. temporal attention mechanism
88. X: [batch_size, num_step, N, D]
89. STE: [batch_size, num_step, N, D]
90. K: number of attention heads
91. d: dimension of each attention outputs
92. return: [batch_size, num_step, N, D]
93. '''
94. D = K * d
95. X = tf.concat((X, STE), axis=-1)
96. # [batch_size, num_step, N, K * d]
97. query = FC(
98. X, units=D, activations=tf.nn.relu,
99. bn=bn, bn_decay=bn_decay, is_training=is_training)
100. key = FC(
101. X, units=D, activations=tf.nn.relu,
102. bn=bn, bn_decay=bn_decay, is_training=is_training)
103. value = FC(
104. X, units=D, activations=tf.nn.relu,
105. bn=bn, bn_decay=bn_decay, is_training=is_training)
106. # [K * batch_size, num_step, N, d]
107. query = tf.concat(tf.split(query, K, axis=-1), axis=0)
108. key = tf.concat(tf.split(key, K, axis=-1), axis=0)
109. value = tf.concat(tf.split(value, K, axis=-1), axis=0)
110. # query: [K * batch_size, N, num_step, d]
111. # key: [K * batch_size, N, d, num_step]
112. # value: [K * batch_size, N, num_step, d]
113. query = tf.transpose(query, perm=(0, 2, 1, 3))
114. key = tf.transpose(key, perm=(0, 2, 3, 1))
115. value = tf.transpose(value, perm=(0, 2, 1, 3))
116. # [K * batch_size, N, num_step, num_step]
117. attention = tf.matmul(query, key)
118. attention /= (d ** 0.5)
119. # mask attention score
120. if mask:
121. batch_size = tf.shape(X)[0]
122. num_step = X.get_shape()[1].value
123. N = X.get_shape()[2].value
124. mask = tf.ones(shape=(num_step, num_step))
125. mask = tf.linalg.LinearOperatorLowerTriangular(mask).to_dense()
126. mask = tf.expand_dims(tf.expand_dims(mask, axis=0), axis=0)
127. mask = tf.tile(mask, multiples=(K * batch_size, N, 1, 1))
128. mask = tf.cast(mask, dtype=tf.bool)
129. attention = tf.compat.v2.where(
130. condition=mask, x=attention, y=-2 ** 15 + 1)
131. # softmax
132. attention = tf.nn.softmax(attention, axis=-1)
133. # [batch_size, num_step, N, D]
134. X = tf.matmul(attention, value)
135. X = tf.transpose(X, perm=(0, 2, 1, 3))
136. X = tf.concat(tf.split(X, K, axis=0), axis=-1)
137. X = FC(
138. X, units=[D, D], activations=[tf.nn.relu, None],
139. bn=bn, bn_decay=bn_decay, is_training=is_training)
140. return X
141.
142.def gatedFusion(HS, HT, D, bn, bn_decay, is_training):
143. '''
144. gated fusion
145. HS: [batch_size, num_step, N, D]
146. HT: [batch_size, num_step, N, D]
147. D: output dims
148. return: [batch_size, num_step, N, D]
149. '''
150. XS = FC(
151. HS, units=D, activations=None,
152. bn=bn, bn_decay=bn_decay,
153. is_training=is_training, use_bias=False)
154. XT = FC(
155. HT, units=D, activations=None,
156. bn=bn, bn_decay=bn_decay,
157. is_training=is_training, use_bias=True)
158. z = tf.nn.sigmoid(tf.add(XS, XT))
159. H = tf.add(tf.multiply(z, HS), tf.multiply(1 - z, HT))
160. H = FC(
161. H, units=[D, D], activations=[tf.nn.relu, None],
162. bn=bn, bn_decay=bn_decay, is_training=is_training)
163. return H
164.
165.def STAttBlock(X, STE, K, d, bn, bn_decay, is_training, mask=False):
166. HS = spatialAttention(X, STE, K, d, bn, bn_decay, is_training)
167. HT = temporalAttention(X, STE, K, d, bn, bn_decay, is_training,
mask=mask)
168. H = gatedFusion(HS, HT, K * d, bn, bn_decay, is_training)
169. return tf.add(X, H)
170.
171.def transformAttention(X, STE_P, STE_Q, K, d, bn, bn_decay, is_training):
172. '''
173. transform attention mechanism
174. X: [batch_size, P, N, D]
175. STE_P: [batch_size, P, N, D]
176. STE_Q: [batch_size, Q, N, D]
177. K: number of attention heads
178. d: dimension of each attention outputs
179. return: [batch_size, Q, N, D]
180. '''
181. D = K * d
182. # query: [batch_size, Q, N, K * d]
183. # key: [batch_size, P, N, K * d]
184. # value: [batch_size, P, N, K * d]
185. query = FC(
186. STE_Q, units=D, activations=tf.nn.relu,
187. bn=bn, bn_decay=bn_decay, is_training=is_training)
188. key = FC(
189. STE_P, units=D, activations=tf.nn.relu,
190. bn=bn, bn_decay=bn_decay, is_training=is_training)
191. value = FC(
192. X, units=D, activations=tf.nn.relu,
193. bn=bn, bn_decay=bn_decay, is_training=is_training)
194. # query: [K * batch_size, Q, N, d]
195. # key: [K * batch_size, P, N, d]
196. # value: [K * batch_size, P, N, d]
197. query = tf.concat(tf.split(query, K, axis=-1), axis=0)
198. key = tf.concat(tf.split(key, K, axis=-1), axis=0)
199. value = tf.concat(tf.split(value, K, axis=-1), axis=0)
200. # query: [K * batch_size, N, Q, d]
201. # key: [K * batch_size, N, d, P]
202. # value: [K * batch_size, N, P, d]
203. query = tf.transpose(query, perm=(0, 2, 1, 3))
204. key = tf.transpose(key, perm=(0, 2, 3, 1))
205. value = tf.transpose(value, perm=(0, 2, 1, 3))
206. # [K * batch_size, N, Q, P]
207. attention = tf.matmul(query, key)
208. attention /= (d ** 0.5)
209. attention = tf.nn.softmax(attention, axis=-1)
210. # [batch_size, Q, N, D]
211. X = tf.matmul(attention, value)
212. X = tf.transpose(X, perm=(0, 2, 1, 3))
213. X = tf.concat(tf.split(X, K, axis=0), axis=-1)
214. X = FC(
215. X, units=[D, D], activations=[tf.nn.relu, None],
216. bn=bn, bn_decay=bn_decay, is_training=is_training)
217. return X
218.
219.def GMAN(X, TE, SE, P, Q, T, L, K, d, bn, bn_decay, is_training):
220. '''
221. GMAN
222. X： [batch_size, P, N]
223. TE： [batch_size, P + Q, 2] (time-of-day, day-of-week)
224. SE： [N, K * d]
225. P： number of history steps
226. Q： number of prediction steps
227. T： one day is divided into T steps
228. L： number of STAtt blocks in the encoder/decoder
229. K： number of attention heads
230. d： dimension of each attention head outputs
231. return： [batch_size, Q, N]
232. '''
233. D = K * d
234. # input
235. X = tf.expand_dims(X, axis=-1)
236. X = FC(
237. X, units=[D, D], activations=[tf.nn.relu, None],
238. bn=bn, bn_decay=bn_decay, is_training=is_training)
239. # STE
240. STE = STEmbedding(SE, TE, T, D, bn, bn_decay, is_training)
241. STE_P = STE[:, : P]
242. STE_Q = STE[:, P:]
243. # encoder
244. for _ in range(L):
245. X = STAttBlock(X, STE_P, K, d, bn, bn_decay, is_training)
246. # transAtt
247. X = transformAttention(
248. X, STE_P, STE_Q, K, d, bn, bn_decay, is_training)
249. # decoder
250. for _ in range(L):
251. X = STAttBlock(X, STE_Q, K, d, bn, bn_decay, is_training)
252. # output
253. X = FC(
254. X, units=[D, 1], activations=[tf.nn.relu, None],
255. bn=bn, bn_decay=bn_decay, is_training=is_training)
256. return tf.squeeze(X, axis=3)
257.
258.def mae_loss(pred, label):
259. mask = tf.not_equal(label, 0)
260. mask = tf.cast(mask, tf.float32)
261. mask /= tf.reduce_mean(mask)
262. mask = tf.compat.v2.where(
263. condition=tf.math.is_nan(mask), x=0., y=mask)
264. loss = tf.abs(tf.subtract(pred, label))
265. loss *= mask
266. loss = tf.compat.v2.where(
267. condition=tf.math.is_nan(loss), x=0., y=loss)
268. loss = tf.reduce_mean(loss)
269. return loss

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可