有研究表明, 受前体物排放^[1~3]和气象条件^[4~7]等因素共同驱动, 大气O₃已成为影响城市夏季环境空气质量的主要污染物^{[8, 9]}. 目前对于O₃污染驱动因素的研究, 主要应用物理化学机制驱动的演绎模型^{[10, 11]}和数据驱动的归纳模型^{[12, 13]}. 其中CMAQ和WRF-Chem等应用相对较多. 张迎春等^[14]采用CMAQ模型对成都市夏季O₃进行模拟评估, 并采用情景分析法得出夏季O₃对气象因子的灵敏度依次为: 气温、绝对湿度和大气压强. Hu等^[15]利用WRF-Chem模拟四川盆地O₃时空变化特征和传输路径. 但这些模型需要输入网格化排放清单和气象场等众多模型参数, 运算时效性较差^[16]. 机器学习ML模型善于处理非线性问题, 运算效率高^[17]. 蔡旺华^[18]和马润美等^[19]分别采用XGBoost和Random Forest模型预测O₃浓度, 精确度较高. 但这些模型的可解释性差, 都较难定量解析各因素对O₃的影响^[20]. 近年来, 基于可解释性数据驱动模型, 学者们将机器学习ML模块耦合可解释性SHAP模块, 为模型的合理性验证和模拟数据的定量解析提供接口^{[21, 22]}. ML-SHAP模型应用广泛, 如Xu等^[23]采用XGBoost-SHAP解释驾驶行为和机动车排放水平的关系, 廖彬等^[24]采用XGBoost-SHAP模型解释决定运动员身价的主要因素. XGBoost-SHAP模型在定量解析影响因素贡献方面已有一定应用, 但该类模型缺乏对各驱动因素和目标变量之间的关联性筛选和优化能力, 导致部分影响因素可解释性较差.

本研究建立可解释性数据驱动模型Correlation-ML-SHAP模型, Correlation模块挖掘O₃浓度的关联影响因素, 优化数据组合, 提高模型可解释性. 机器学习ML模块耦合可解释性SHAP模块计算各驱动因素对O₃浓度的贡献, 并以晋城市2021年夏季O₃污染为例开展应用研究, 旨在为夏季城市大气O₃污染防控提供基础依据.

Correlation-ML-SHAP研究方法耦合Correlation模块、ML模块和可解释性SHAP模块, 通过挖掘O₃浓度关联影响因素, 实现对模型结果可解释性的优化. 以上3个模块分别用于时间关联性挖掘、O₃浓度模拟和各驱动因素影响评估. 应用Correlation模块遍历不同时段驱动因素和模拟时段O₃浓度的Pearson相关系数, 搜索O₃浓度关联影响因素优化数据组合. 应用机器学习ML模块训练候选模型, 结合模型评价指标和可解释性SHAP特征图分别评估候选模型的精准度和可解释性, 确定最佳模型. 将优化数据组合和最佳模型输入可解释性SHAP模块计算任意时刻各驱动因素对O₃浓度的影响贡献. 应用于定量解析O₃浓度驱动因素和交互影响. 技术路线如图 1所示.

图 1

Fig. 1

图 1 技术路线 Fig. 1 Technology roadmap

对于数据驱动的机器学习ML模型, 设O₃污染数据集为{X}={x_i1, x_i2, x_i3, …, x_is, y_i}.对于i时刻的O₃浓度可由s个驱动因素输入候选机器学习ML模型计算得到^[25]. 如公式(1):

(1)

式中, O_3(i)为i时刻O₃浓度的模拟值(μg·m^-3); α为候选ML模型的优化参数配置; f为候选的机器学习ML模型.

各驱动因素对夏季城市O₃浓度的影响呈现非线性变化特征^{[26, 27]}. 处理此类数据特征的主流ML模型包括XGBoost、Random Forest、Adaboost、Decision Tree和KNN^[28~32]. 采用网格搜索验证法优化以上候选模型的参数配置, 如表 1所示. 以XGBoost模型为例, n_estimator为弱学习器数量, 通过生成新的学习器拟合上一步结果的残差以逐步提高模拟精度; max_depth、reg_alpha和reg_lambda分别为树的最大深度以及L1和L2正则化项, 这些参数主要用于避免模型发生过拟合.

表 1 (Table 1)

表 1 模型主要参数配置 Table 1 Configuration of the main parameters of the models 模型名称 参数配置 XGBoost n_estimators=210, max_depth=7, min_child_weight=6, subsample=1.0, colsample_bytree=1.0, gamma=0.6, reg_alpha=0.05, reg_lambda=3, learning_rate=0.15 Random Forest n_estimators=210, max_depth=22, max_features=5, max_features=5, min_samples_leaf=2 Adaboost n_estimators=177, learning_rate=1 Decision Tree criterion=‘mse’, min_samples_split=2 KNN n_neighbors=10, p=1, weights=distance

表 1 模型主要参数配置 Table 1 Configuration of the main parameters of the models

为定量评价模型模拟结果的准确度, 分别采用平均值绝对误差(MAE)、均方根误差(RMSE)和决定系数(R²)对候选ML模型进行定量评价. 如公式(2)~(4):

(2)

(3)

(4)

式中, O_{3(i true)}为i时刻O₃浓度的真实值(μg·m^-3); O_3(i)为i时刻O₃浓度的模拟值(μg·m^-3). MAE和RMSE用于评价模拟值和真实值的差异, 两种指标越低越准确; R²用于评价回归性能, 取值越大效果越好. 表 2为候选的不同机器学习ML模型模拟结果的评估指标.

表 2 (Table 2)

表 2 模型评价指标结果 Table 2 Model evaluation index results 模型名称 MAE RMSE R2 XGBoost 15.4 21.7 0.82 Random Forest 16.0 21.9 0.81 Adaboost 22.2 28.2 0.69 Decision Tree 21.8 30.4 0.64 KNN 22.2 30.1 0.65

表 2 模型评价指标结果 Table 2 Model evaluation index results

KNN、Decision Tree和Adaboost模型模拟结果的准确度较低, 主要原因在于KNN对多维数据的回归和可解释性较差; Decision Tree和Adaboost模型易受异常数据干扰, 稳定性较弱. XGBoost和Random Forest模型模拟结果的准确度相对较高, 但由于Random Forest模型max_depth参数值较大, 导致模型复杂度高, 容易发生过拟合. 采用XGBoost模型时MAE、RMSE和R²这3项指标分别为15.4、21.7和0.82, 相比其它ML模型模拟结果的准确度最高, 原因在于该模型通过引入损失函数的二阶泰勒展开和正则项, 优化模型复杂度, 避免过拟合^{[33, 34]}.

基于合作博弈中SHAP理论的归因分析方法, 计算机器学习ML模型中各驱动因素对O_3(i)的边际贡献值, 结合SHAP可加性解释模型量化不同时刻各驱动因素的贡献大小^[35], 计算方法如公式(5)~(7):

(5)

(6)

(7)

式中, O_3(base)为O_3(i)的期望值(μg·m^-3); shap(x_{i, j})为i时刻因素j对O_3(i)的贡献值(μg·m^-3), 当shap(x_{i, j})>0时, 对O_3(i)呈正贡献; shap(x_{i, j}) < 0时, 对O_3(i)呈负贡献; O_3(M)为ML模型在子集M条件下的模拟值(μg·m^-3); M为不包括因素j的特征子集; β_(M)为特征子集M的权重.

环境中不同驱动因素对O_3(i)的贡献权重不同, 定量解析因素j对O_3(i)的贡献权重, 可以通过计算因素j的SHAP特征重要性与全部因素重要性的比值得到^[36]. 计算方法如公式(8)和公式(9):

(8)

(9)

式中, I(x_j)为因素j对O_3(i)的贡献权重(%); α(x_j)为因素j的SHAP特征重要性, 以因素j的|shap(x_{i, j})|均值计(μg·m^-3), α(x_j)取值越大, 因素j对O_3(i)的贡献权重越大.

环境中驱动因素通常在其他因素交互影响下共同驱动O_3(i)变化, 分析O_3(i)驱动因素的交互影响^[37]. 计算方法如公式(10)~(12):

(10)

(11)

(12)

式中, shap(x_i, x_j)为因素i和j对O_3(i)的交互影响值(μg·m^-3), O_3(M^*)为ML模型在子集M^*条件下的模拟值(μg·m^-3); M^*为不包括因素i和j的特征子集; β^*为特征子集M^*的权重.

(a) 2021年各季节臭氧日最大8 h第90百分位数; (b)2021年各月臭氧日最大8 h第90百分位数; (c)2021年5~9月每小时O₃浓度均值

晋城市作为“2+26”城市之一, 位于山西省东南部, 全境被山脊环绕形成盆地地形. 如图 2所示, 2021年臭氧日最大8 h第90百分位数[ρ(O₃-8h)]为180 μg·m^-3, 位居“2+26”城市倒数. ρ(O₃-8h)呈“春夏高、秋冬低”分布, 即: 夏季(212 μg·m^-3)>春季(172 μg·m^-3)>秋季(153 μg·m^-3)>冬季(96 μg·m^-3). 其中5~9月是污染高发期, ρ(O₃-8h)分别为185、225、192、179和177 μg·m^-3. 2021年5~9月每小时O₃浓度均值变化特征呈单峰分布, 07:00开始浓度持续升高, 12:00~18:00处于高值区. 主要原因可能是夏季昼间有利的气象条件促进光化学反应向O₃生成方向进行^[38~42].

图 2

Fig. 2

(a)2021年各季节臭氧日最大8 h第90百分位数; (b)2021年各月臭氧日最大8 h第90百分位数; (c)2021年5~9月每小时O3浓度均值 图 2 晋城市2021年O3浓度时间变化特征 Fig. 2 Temporal variation characteristics of O3 concentration in Jincheng City in 2021

以晋城市2021年夏季O₃污染过程为例, 应用建立的Correlation-ML-SHAP模型定量解析夏季O₃污染驱动因素影响贡献.

选取晋城市2021年夏季5~9月昼间07:00~18:00作为模拟时段. 收集O₃、NO₂、CO、PM₁₀和PM_2.5浓度的小时数据以及气温、湿度、日照强度、风向、风速、大气压和能见度小时数据作为O₃污染数据集. 小时浓度数据来源于大气环境科学综合数据采集与共享平台(http://106.38.83.12:9002/ahpcp/login/index.html); 小时气象数据来源于国家气象科学数据中心(http://data.cma.cn/). 数据集中存在的部分异常、缺失值, 采用基于时间序列的拉格朗日插值法填充处理.

如表 3所示, 应用Correlation-XGBoost-SHAP模型的Correlation模块, 遍历不同时段各驱动因素与模拟时段O₃浓度之间的相关性系数, 得到关联因素为: 提前2 h的气温、日照强度、湿度以及同时刻的NO₂、CO、大气压、风速和能见度, Pearson相关系数分别为0.777、0.770、-0.611、-0.548、-0.387、-0.346、0.326和0.331, 在0.01水平上显著相关. 数据特征方面, 气温、日照强度、湿度的相关性在提前2 h达到最大, 表明O₃浓度变化滞后于上述特征. 姜峰等^[43]在研究日照强度对O₃浓度影响时也发现了这种关系. 可知应用数据驱动的方法对O₃浓度进行模拟时, 模型可能更关注提前2 h的气温、日照强度和湿度变化.

表 3 (Table 3)

表 3 驱动因素与O3浓度的相关性系数1) Table 3 Correlation coefficient between driving factors and O3 concentration1) 同时段 提前1 h 提前2 h 提前3 h 提前4 h 提前5 h 气温 0.694** 0.762** 0.777** 0.751** 0.694** 0.614** 日照强度 0.539** 0.711** 0.770** 0.757** 0.707** 0.642** 湿度 -0.514** -0.583** -0.611** -0.600** -0.558** -0.494** 风速 0.326** 0.319** 0.287** 0.231** 0.162** 0.109** 大气压 -0.346** -0.323** -0.286** -0.245** -0.211** -0.188** 能见度 0.331** 0.344** 0.339** 0.324** 0.291** 0.235** NO2 -0.548** -0.475** -0.372** -0.273** -0.177** -0.058 CO -0.387** -0.278** -0.141** -0.012 0.096** -0.090** PM10 -0.020 -0.012 0.015 0.050 0.094** -0.024 PM2.5 -0.017 0.028 0.083** 0.141** 0.195** 0.023 1)**表示在0.01水平上显著相关

表 3 驱动因素与O3浓度的相关性系数1) Table 3 Correlation coefficient between driving factors and O3 concentration1)

如图 3所示, 由于XGBoost-SHAP模型缺乏对各驱动因素和O₃浓度之间的时间关联性筛选和优化能力, 导致部分驱动因素的可解释性较差. 建立的Correlation-XGBoost-SHAP模型关注数据之间的时间关联性, 通过Correlation模块挖掘O₃浓度的关联驱动因素, 优化数据组合, 使模型模拟结果的合理性得以验证和改善. 从模拟结果的可解释性来看, 该模型优化了影响O₃浓度变化的驱动因素贡献特征, 如该模型的日照强度贡献排名提升至第2, 对O₃污染呈正贡献. 从模拟结果的准确度来看, 该模型的MAE、RMSE和R²分别为15.4、19.8和0.85, 在XGBoost-SHAP模型的基础上准确度有所提高. 由此可知Correlation-ML-SHAP模型在解析影响O₃污染的各驱动因素方面具备可行性, 其中ML模块采用XGBoost模型模拟准确度最佳.

图 3

Fig. 3

(a)模型可解释性; (b)模型评价指标 图 3 优化数据组合的模型性能评估 Fig. 3 Model performance evaluation of optimized data composition

应用Correlation-XGBoost-SHAP模型解析晋城市夏季O₃污染因素. 如图 4所示, 通过可解释性SHAP模块公式(8)和公式(9)计算各驱动因素的贡献权重, 用NO₂和CO对O₃的贡献权重和代表前体物排放水平的贡献权重. 影响O₃浓度变化的各驱动因素贡献权重依次为: 32.1%(气温)、21.3%(日照强度)、16.5%(湿度)、15.6%(前体物的排放水平)、8.5%(风向风速)和6.0%(其它). 气温、日照强度、湿度以及前体物对O₃污染的贡献权重均高于传输扩散相关的风向风速, 表明晋城市夏季O₃浓度变化主要受本地光化学反应生成影响. 其中前体物的排放贡献权重在污染天排名提升至第3, 其它驱动因素的贡献权重排名不变, 但贡献程度不同. 气温、日照强度和前体物排放水平的贡献权重在污染天分别为34.0%、21.9%和16.3%, 高于优良天3.4%、1.2%和1.2%; 湿度的贡献权重在污染天为14.3%, 低于优良天4.0%. 可知夏季城市O₃浓度受气温和日照强度影响较大, 污染天前体物排放对O₃浓度的影响也不容忽视.

图 4

Fig. 4

图 4 晋城市夏季O3浓度的驱动因素贡献权重 Fig. 4 Contribution weights of driving factors for O3 concentration in Jincheng City in summer

通过Correlation-XGBoost-SHAP模型可解释性SHAP模块的交互值计算, 进一步探究各驱动因素对O₃浓度变化的交互影响. 如图 5所示, 气象条件方面, 气温、日照强度和O₃浓度呈非线性正相关, 随着日照强度的增强, 气温对O₃污染的正贡献提升, 当气温超过24℃时在较强日照的交互影响下有94.9%的概率对O₃污染呈正贡献; 湿度低于70%时有94.1%的概率对O₃污染呈正贡献, 但超过70%后在低温的交互影响下对O₃污染呈负贡献. 传输扩散方面, 东南风风速低于5.8m·s^-1或南风风速低于5.3m·s^-1时, 风速增加对O₃污染的正贡献提升. 由此可知, 晋城市夏季O₃污染受多方面的气象因素交互影响, 强日照条件下高温和低湿是引发O₃污染的原因, 同时还受到东南部和南部区域的传输影响. 与周边城市临汾^[44]、郑州^[45]的O₃污染条件基本一致. 前体物方面, NO₂、CO和O₃浓度的变化相近. 高温(≥24℃)且低湿(≤70%)的气象条件下, ρ(NO₂)≥9 μg·m^-3或ρ(CO)≥0.7 mg·m^-3, 分别有94.9%和99.3%的概率对O₃污染呈正贡献.

图 5

Fig. 5

散点为横坐标上因素j在其他因素共同驱动下对O3的贡献值; 色柱为产生交互影响的因素取值; (d)和(e)用虚线将反应条件加以划分: Ⅰ阶段为午后的光化学反应, Ⅱ阶段为日出后的光化学反应 图 5 晋城市夏季O3浓度的驱动因素交互影响 Fig. 5 Interaction of driving factors on O3 concentration in Jincheng City in summer

(1) 基于可解释性数据驱动思想, 以城市O₃污染过程为研究对象, 建立Correlation-ML-SHAP模型, 模型能够挖掘并利用强驱动因素模拟O₃浓度和量化影响贡献, 其中机器学习ML模块采用XGBoost模型模拟准确度最佳, 模型可为夏季城市大气O₃污染防控提供基础依据.

(2) 晋城市2021年夏季5~9月ρ(O₃-8h)分别为185、225、192、179和177 μg·m^-3. 影响夏季O₃浓度的强驱动因素依次为气温、日照强度、湿度和前体物排放水平, 贡献权重分别为32.1%、21.3%、16.5%和15.6%.

(3) 气温、日照强度、前体物排放水平在污染天对O₃浓度的贡献权重分别为34.0%、21.9%和16.3%, 相比优良天分别提升3.4%、1.2%和1.2%, 前体物排放贡献权重排名提升至第3.

(4) 各驱动因素对O₃浓度呈非线性交互影响, 气温超24℃或湿度低于70%时分别有94.9%和94.1%的概率对O₃污染呈正贡献, 这种气象条件下ρ(NO₂)≥9 μg·m^-3或ρ(CO)≥0.7 mg·m^-3, 分别有94.9%和99.3%的概率对O₃污染呈正贡献. 东南风风速低于5.8m·s^-1或南风风速低于5.3m·s^-1, 均对O₃污染呈正贡献.