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的污染物浓度预报准确度？说明原因。

  题中所给附件和附录我就不放这上面了

解题

  问题一主要是通过所提供的计算公式进行简单的数据验算。首先根据各项污染物的空气质量分指数的计算公式求出每一类污染物的 IAQI ，公式如下：

  其中臭氧（O3）浓度需要用其最大 8 小时滑动平均值来计算：臭氧滑动平均指一个自然日内 8 时至 24 时的所有 8 小时滑动平均浓度中的最大值，其中 8 小时滑动平均值指连续 8 小时平均浓度的算术平均值。其计算公式如下：

  其中c+为臭氧在某日𝑡 − 1时至𝑡时的平均污染物浓度。

  通过对比空气质量分指数（IAQI）与对应的污染物项目浓度限值表 4-1 来对计算出来的空气质量分指数，从而选择某污染物浓度相近的污染物浓度限值的高位值与低位值，接着通过空气质量指数公式取各分指数中的最大值，即


  问题二要求处理两个方面：（1）对所提供的六种污染物浓度、每一种气象条件与 AQI 进行灰色关联度分析，筛选出与影响污染物浓度相关的关键因素，忽略关联性不高的次要因素，从而达到对数据进行降维的目的；（2）要求根据的降维的特征对气象条件进行合理分类，此处用聚类实现分类，并阐述各类气象条件的特征。
  附件 1 中实测和预报的数据，都因为各种原因，导致存在缺失值，而时间序列的数据， 不能简单的删除去处理缺失值，因此要进行时间序列缺失值的填充插值处理。一般而言， 有以下几种方法：对时间序列的缺失值插值法有零填充、前向填充、后向填充、knn 最近邻均值法填充以及单线性插值填充等。本论文中我们选择采用 knn 最近邻均值法填充和单线性插值填充法进行缺失值处理。


  六种污染物的拟合曲线图

  Calinski-Harabasz 分数值 s 的数学计算公式是：

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

程序代码

# k-means 聚类
from numpy import unique
from numpy import where
from sklearn.datasets import make_classification
from sklearn.cluster import KMeans,MiniBatchKMeans,DBSCAN,AgglomerativeClustering,Birch
from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn import metrics
import pandas as pd
import numpy as np 
import scipy.io as io
# 引入必要的第三方库及包
# 数据处理后训练数据集
# 采用kmeans聚类
# from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import Axes3D
def minibatchKmeans(X):
    for index, k in enumerate((3,4,5,6)):
        plt.subplot(2,2,index+1)
        y_pred = MiniBatchKMeans(n_clusters=k, batch_size = 200, random_state=9).fit_predict(X)
        score= metrics.calinski_harabasz_score(X, y_pred)  
        plt.scatter(X[:, 0], X[:, 1], c=y_pred)
        plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                    transform=plt.gca().transAxes, size=10,
                    horizontalalignment='right')
    plt.show()

def Kmeans3D(data):

    estimator = KMeans(n_clusters=3)  
    # 构造聚类器
    y = estimator.fit_predict(data)  
    # 聚类
    label_pred = estimator.labels_  
    # 获取聚类标签
    centroids = estimator.cluster_centers_  
    # 获取聚类中心
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(data[:, 0], data[:, 1], data[:, 2], c=y, marker='*')
    ax.scatter(centroids[:, 0], centroids[:, 1], centroids[:, 2], marker='>')
    plt.axis([0, 1, 0, 1])
    plt.show()
def minibatchKmeans(X,k,p=False):
    if p==True:
        for index, kk in enumerate((3,4,5,6)):
            plt.subplot(2,2,index+1)
            y_pred = MiniBatchKMeans(n_clusters=kk).fit_predict(X)
            score= metrics.calinski_harabasz_score(X, y_pred)  
            plt.scatter(X[:, 0], X[:, 1], c=y_pred)
            plt.text(.99, .01, ('k=%d, score: %.2f' % (kk,score)),
                        transform=plt.gca().transAxes, size=10,
                        horizontalalignment='right')
        plt.show()
    else:
        y_pred = KMeans(n_clusters=k).fit_predict(X)
        clusters = unique(y_pred)
        indexdict ={}
        for i in clusters:
            indexarr = np.where(y_pred[:]==i)
            indexdict[i] = indexarr[0]
        return indexdict
def P_Kmeans(X,k,p=False):
    if p==True:
        for index, kk in enumerate((3,4,5,6)):
            plt.subplot(2,2,index+1)
            y_pred = KMeans(n_clusters=kk).fit_predict(X)
            score= metrics.calinski_harabasz_score(X, y_pred)  
            plt.scatter(X[:, 0], X[:, 1], c=y_pred)
            plt.text(.99, .01, ('k=%d, score: %.2f' % (kk,score)),
                        transform=plt.gca().transAxes, size=10,
                        horizontalalignment='right')
        plt.show()
    else:
        y_pred = KMeans(n_clusters=k).fit_predict(X)
        clusters = unique(y_pred)
        indexdict ={}
        for i in clusters:
            indexarr = np.where(y_pred[:]==i)
            indexdict[i] = indexarr[0]
        return indexdict
def Gau(X):
    for index, k in enumerate((3,4,5,6)):
        plt.subplot(2,2,index+1)
        y_pred = GaussianMixture(n_components=k).fit_predict(X)
        score= metrics.calinski_harabasz_score(X, y_pred)  
        plt.scatter(X[:, 0], X[:, 1], c=y_pred)
        plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                    transform=plt.gca().transAxes, size=10,
                    horizontalalignment='right')
    plt.show()
    y_pred = GaussianMixture(n_components=k).fit_predict(X)
   

def P_DBSCAN(X):
    for index, k in enumerate((3,4,5,6)):
        plt.subplot(2,2,index+1)
        y_pred = DBSCAN(eps=0.30, min_samples=k).fit_predict(X)
        score= metrics.calinski_harabasz_score(X, y_pred)  
        plt.scatter(X[:, 0], X[:, 1], c=y_pred)
        plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                    transform=plt.gca().transAxes, size=10,
                    horizontalalignment='right')
    plt.show()
def Agg(X):
    # 聚合聚类
    for index, k in enumerate((3,4,5,6)):
        plt.subplot(2,2,index+1)
        y_pred = AgglomerativeClustering(n_clusters=k).fit_predict(X)
        score= metrics.calinski_harabasz_score(X, y_pred)  
        plt.scatter(X[:, 0], X[:, 1], c=y_pred)
        plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                    transform=plt.gca().transAxes, size=10,
                    horizontalalignment='right')
    plt.show()
def P_Birch(X):
    # BIRCH 聚类（ BIRCH 是平衡迭代减少的缩写，聚类使用层次结构)包括构造一个树状结构，从中提取聚类质心。
    for index, k in enumerate((3,4,5,6)):
        plt.subplot(2,2,index+1)
        y_pred = Birch(threshold=0.01, n_clusters=k).fit_predict(X)
        score= metrics.calinski_harabasz_score(X, y_pred)  
        plt.scatter(X[:, 0], X[:, 1], c=y_pred)
        plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                    transform=plt.gca().transAxes, size=10,
                    horizontalalignment='right')
    plt.show()
if __name__ =="__main__":
    # 定义数据集
    data_1_predict_raw = pd.read_excel('./data/data_1_knn_IAQI.xlsx',sheet_name='监测点A逐小时污染物浓度与气象实测数据')
    # X_out = data_1_predict_raw[['so2','no2','pm10','pm2.5','o3']]#,'co']]#,'direction']]#,'pressure'
    # X_out = data_1_predict_raw[['so2','no2','pm10','pm2.5','o3',]]#,'IAQI_SO_2','IAQI_NO_2','IAQI_PM_10','IAQI_PM_25','IAQI_O_3']]#,'co']]#,'direction']]#,'pressure'
    # X_out = data_1_predict_raw[['so2','no2','pm10','pm2.5','o3','temperature','humidity','pressure','wind','direction']]
    # X_out = data_1_predict_raw[['temperature','humidity','wind','direction','IAQI_SO_2','IAQI_NO_2','IAQI_PM_10','IAQI_PM_25','IAQI_O_3']]
    # X_out = data_1_predict_raw[['IAQI_SO_2','IAQI_NO_2','IAQI_PM_10','IAQI_PM_25','IAQI_O_3']]
    
    X_out = X_out.values

    scaler = StandardScaler()      
    X = scaler.fit_transform(X_out)  
    # Kmeans3D()
    # minibatchKmeans(X)
    # Gau()
    # Agg(X)
    # P_Birch(X)
    k = 3
    indexdict = P_Kmeans(X,k,True)
    # indexdict = minibatchKmeans(X,k,True)
    X_all = data_1_predict_raw[['so2','no2','pm10','pm2.5','o3','co','temperature','humidity','pressure','wind','direction','aqi']].values
    # 生成mat文件到MATLAB中绘图
    io.savemat('./IAQI_data_3f/X_3all.mat', {'X_all': X_all})
    io.savemat('./IAQI_data_3f/index1.mat', {'index1': indexdict[0]})
    io.savemat('./IAQI_data_3f/index2.mat', {'index2': indexdict[1]})
    io.savemat('./IAQI_data_3f/index3.mat', {'index3': indexdict[2]})

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