1 引言

多尺度空气质量模型(CMAQ)作为MODES-3的核心组成部分，可用于局地、城市、区域和大陆等多尺度的空气质量预报研究。为了提高模式预报精度，使模型能够真实反映研究区域的大气污染过程，中外已经有许多学者开展了相关研究，如许建明等(2005)进行的CMAQ-MOS区域空气质量统计修正模型预报途径研究，通过观测气象要素值以及观测污染物浓度值对模型模拟结果进行统计订正，成功将概率预报方式与数值预报方式结合起来；Peng等(2014)基于CMAQ模型成功建立了东亚地区的碳同化反演系统；Chang等(2016)使用了同化后的MACC(Monitoring Atmospheric Composition and Climate)数据作为初始和边界条件，成功地对CMAQ模型进行优化；Gou等(2011)使用四维变分方法优化了CMAQ模型的化学传输过程；Zubrow等(2008)使用了集合调整卡尔曼滤波方法对CMAQ模型预报中CO的化学传输过程进行优化；针对CMAQ模型所使用的排放源清单，程兴宏等(2013, 2016)提出了使用自适应偏最小二乘回归法以及张弛逼近源同化反演对排放源清单进行订正的方法。

为了提高多尺度空气质量模型(CMAQ)的模拟准确度，本研究引入了结合WRF气象要素的CMAQ-MOS模型。在模拟试验过程中，发现传统的多元线性回归方法限制了MOS模型的优化效果，采用非线性的极端随机树方法和梯度提升回归树方法对模式输出统计模型进行优化，将两种方法应用到徐州空气质量预报业务中，并进行性能比较。

2 模式和方法 2.1 WRF-CMAQ-MOS模型

传统的多元线性回归方法限制了结合WRF气象要素的CMAQ-MOS模型的发挥，多元线性回归方法可较好地预测线性过程，但对于非线性过程，例如大气污染物预报过程，其效果显然会受到较大的限制。为了解决以上问题，提升结合WRF气象要素的CMAQ-MOS模型优化效果，文中选取较为新颖且易于实现的极端随机树(Extremely Randomized Trees)方法(Galelli, et al, 2013；Scalzo, et al, 2012)和梯度提升回归树(Gradient Boosted Regression Trees)方法(Carslaw, et al, 2009)，对结合WRF气象要素的CMAQ-MOS模型进行改进。传统的使用多元线性回归方法、基于WRF预报气象要素的MOS模型计算式为

(1)

(2)

式中，X_t|t-t0代表t|t-t₀时刻至t时刻某种污染物的浓度观测值矩阵, 代表t-t₀时刻至t时刻常数矩阵，代表t-t₀时刻至t时刻N列WRF大气要素预报值、L列空气质量模型污染物浓度模拟值，Ã_mt|t-t0^N+L为t-t₀时刻至t时刻相对应系数项矩阵。根据式(1)，求得和Ã_mt|t-t0^N+L后代入式(2)。在式(2)中，S_t+t1|t为t时刻至t+t₁时刻某种污染物的浓度预报值矩阵，本系统中，取t₀为1×24 h，取t₁为1×24 h，代表t时刻至t+t₁时刻WRF大气要素预报值、空气质量模型污染物浓度模拟值。

由式(1)、(2)可知，使用MOS模型综合考虑了污染物浓度观测数据、大气要素观测数据、WRF预报气象数据和空气质量模型预报数据进行污染物浓度预报订正，考虑了动力作用、化学作用和污染物排放的影响。

图 1为结合WRF预报气象要素的CMAQ-MOS模型运行流程，分为模拟阶段和预报阶段。在模拟阶段，模型读取t-1时段观测气象数据，t时段CMAQ模型预报污染物浓度数据，t时段WRF预报气象数据，t时段观测污染物浓度数据，建立基于多元线性回归的矩阵方程，其中，观测气象要素数据为4×25的矩阵，4列分别为观测温度、湿度、气压、风速，25行为从t-1天00时至t天00时的时间序列；WRF气象要素数据为4×25的矩阵，4列分别为WRF预报温度、湿度、气压、风速，25行为从t天00时至t+1天00时的时间序列；CMAQ模型预报污染物浓度数据为7×25的矩阵，7列分别为空气质量指数(AQI)、PM_2.5、PM₁₀、CO、NO₂、O₃、SO₂的模型模拟值，25行为从t天00时至t+1天00时的时间序列；观测污染物浓度数据为7×25的矩阵，7列分别为空气质量指数、PM_2.5、PM₁₀、CO、NO₂、O₃、SO₂的观测值，25行为从t天00时至t+1天00时的时间序列。通过以上数据组成多元线性回归模型，将所得模型代入预报阶段中，得到t+1时段MOS模型污染物预报浓度值，为7×25的矩阵，7列分别为空气质量指数、PM_2.5、PM₁₀、CO、NO₂、O₃、SO₂的污染物浓度预报值，25行为从t天00时至t+1天00时的时间序列。

2.1.1 极端随机树方法

Geurts等(2006)提出了极端随机树方法，根据经典的自上而下的方法，极端随机树构建了一系列“自由生长”的回归树集合。与随机森林方法相似，极端随机树方法也是由多棵决策树构成，但与随机森林方法不同在于极端随机树方法是完全随机的得到分叉值，从而进行对回归树的分叉，不同于随机森林的在一个随机子集内得到最佳分叉属性。此外，极端随机树方法中的每一棵回归树用的都是全部训练样本。该方法的实现流程为：

(1) 根据观测污染物浓度数据建立随机拆分(S, a)

计算S的最大值和最小值，分别记为a_min^S和a_max^S；选择均匀切点为[a_min^S，a_max^S]中的a_c；返回拆分[a < a_c]。

(2) 结合WRF预报气象数据、模型预报数据建立一棵极端随机树t

在以下3种情况时，返回叶子结点：(ⅰ)|S| < n_min，即叶子节点的最小样例；(ⅱ)S中包含了所有的候选属性；(ⅲ)S中包含了输出变量。运行流程如下：在所有的属性中选择互不相同的随机K属性{a₁, …, a_K}; 利用随机拆分(S, a_i)，构建K拆分{S₁, …, S_K}；选择拆分S_*，其中，Score(S_*, S)=max_{i=1, …, K}Score(S_i, S)，Score为评分函数，文中，Score_R(S_i, S)=, var{y|S}为y在S中的方差；根据S_*，在分叉S_l和S_r中进行切片S；分叉S_l和S_r中分别形成下一步的回归树t_l和t_r；创建分叉点S_*将t_l和t_r作为左右子树并返回树t。

(3) 建立极端随机树集合

根据(2)构建极端回归树集合T={t₁, …, t_M}。

2.1.2 梯度提升回归树方法

梯度提升回归树是根据梯度提升算法构建的回归树模型，梯度提升是一种组合算法，它的基分类器是决策树，既可以用来回归，也可以用作分类。与随机森林相比，梯度提升算法在分类性能上能够与其媲美，甚至有可能表现得更好。梯度提升回归树具有以下三点优势：自然而然地处理混合类型的数据、预测能力强及在输出空间通过强大的损失函数对于异常值的鲁棒性强。

Friedman(2001)提出梯度提升算法，其核心就在于每棵树是从先前所有树的残差中来学习。利用的是当前模型中损失函数的负梯度值作为提升树算法中的残差的近似值，进而拟合一棵回归树。其运行流程为：

(1) 初始化，结合WRF预报气象数据、模型预报数据和观测污染物浓度数据通过函数

(3)

估计一个使损失函数极小化的常数值，式中，y_i为WRF预报气象数据和模型预报数据，c为观测污染物浓度数据，此时它是只有一个节点的树。

(2) 迭代法建立M棵提升树

在第1层循环中m=1，2，…，M；在第2层循环中i=1，2，…，N，计算损失函数的负梯度在当前模型的值，并将其作为残差的估计值。

(4)

对于r_{m, i}拟合一棵回归树，得到第m棵树的叶节点区域R_{m, j}，j=1，2，…，J；在第2层循环中j=1，2，…，J，计算

(5)

利用线性搜索估计叶节点区域的值，使损失函数极小化；之后，更新

(6)

式中，指示函数I在x∈R_{m, j}时，I=1，否则I=0。

(3) 最后得到的f_m(x)就是最终的梯度提升回归树模型

(7)

2.2 研究区域设置

以六大污染物(NO₂、PM_2.5、PM₁₀、SO₂、O₃、CO)及空气质量指数(AQI)为研究对象，利用空气质量模型CMAQ及Python、Fortran语言，shell脚本语言，进行CMAQ-MOS统计修正模型研究，用于徐州空气质量预报。

模型使用3层嵌套网格(图 2)，水平分辨率分别为27、9和3 km，第1层区域覆盖整个东亚、大部分东南亚和部分南亚地区，第2层区域覆盖大部分华东地区及部分华中地区，第3层区域覆盖徐州市全域，网格数分别为88×76，112×94和97×97；垂直方向采用η坐标系，共29层，模式层顶气压为50 hPa，模式层顶距地面20 km，边界层(约2 km)内分为10层。对于污染源排放数据，区域内层主要采用SMOKE模型输出，外层嵌套采用清华大学MEIC清单数据。

2.3 误差分析方法

采用的误差分析方法主要有4种，分别是均方根误差(RMSE)、相关系数(R)、平均相对误差(MFE)、平均相对偏差(MFB)(Boylan, et al, 2006)。当平均相对误差在±30%内，且平均相对偏差小于50%时，可认为预报是准确的。

(8)

(9)

(10)

(11)

3 WRF-CMAQ-MOS统计修正模型的应用 3.1 极端随机树方法和梯度提升回归树方法比较

选取2016年9月13-24日徐州市黄河新村站点逐时污染物浓度观测数据、WRF气象预报数据、CMAQ模型污染物浓度预报数据，使用极端随机树方法和梯度提升回归树方法并结合WRF气象要素进行CMAQ-MOS模型试验。分别以WRFMOS1代表使用多元线性回归方法、结合WRF气象要素预报的CMAQ-MOS模型优化结果，EXTRA代表使用极端随机树方法、结合WRF气象要素预报的CMAQ-MOS模型优化结果，GBRT代表使用梯度提升回归树方法、结合WRF气象要素预报的CMAQ-MOS模型优化结果，试验结果如下：

由表 1可以看出，EXTRA方案优化效果最佳。结合观测污染物浓度值(OBV)、WRFMOS1方案结果、EXTRA方案结果、GBRT方案结果，绘制出2016年9月13-24日徐州市空气质量指数、CO、NO₂、O₃、PM₁₀、PM_2.5、SO₂浓度观测、WRFMOS1方案、EXTRA方案、GBRT方案时间变化和泰勒图(图 3、4)。

针对CMAQ模型，使用多元线性回归方法，结合WRF预报气象要素的CMAQ-MOS模型，其优化效果仍未达到预期。在选取极端随机树方法和梯度提升回归树方法，分别对结合WRF气象要素的CMAQ-MOS模型进行改进之后，经由试验结果发现，使用极端随机树方法比梯度提升回归树方法的改进效果更显著，其原因主要是因为其每一棵树都是用了全部的样本，并进行了完全随机的拆分。

3.2 极端随机树方法的应用

为了验证使用极端随机树方法并结合WRF气象要素的CMAQ-MOS模型优化效果并进行应用，进行了针对徐州地区2016年1-3月的逐时预报应用试验，选取黄河新村站点(1177A)观测数据进行评估，并针对7个站点的优化效果进行统计分析。模型运行效率较高，进行单站点一个月逐时预报试验耗时约15 min，业务化运行中每日对7个站点进行未来3 d逐时预报，耗时约3 min。

针对空气质量指数及六大污染物，使用极端随机树方法结合WRF气象要素预报的CMAQ-MOS模型进行试验，7个站点试验结果见表 2。可见，经过优化后的六大污染物及空气质量指数的准确度均得到了极大的提高。根据试验结果，针对黄河新村站点的空气质量指数与六种污染物分别绘制2016年1-3月徐州市空气质量指数、CO、NO₂、O₃、PM₁₀、PM_2.5、SO₂污染物浓度预报和订正误差时间变化、散点分布和泰勒图。从图 5-7可明显看出，针对空气质量指数、CO、NO₂、O₃、PM₁₀、PM_2.5、SO₂，使用极端随机树方法结合WRF气象要素预报的CMAQ-MOS模型具有良好的优化效果。

4 总结

为了提高CMAQ空气质量模型的预报精度，基于多元线性回归方法，引入了结合WRF气象要素的CMAQ-MOS模型，通过试验发现其优化效果一般，主要原因是多元线性回归方法限制了结合WRF气象要素的CMAQ-MOS模型的发挥。为了提升模型优化效果，选取极端随机树方法和梯度提升回归树方法，对结合WRF气象要素的CMAQ-MOS模型进行改进，通过试验发现，相对梯度提升回归树方法，极端随机树方法的改进效果更为显著。

采用极端随机树方法，并结合WRF预报气象要素的CMAQ-MOS模型在徐州地区进行2016年1-3月的模拟应用试验。结果表明，针对空气质量指数及六大污染物逐时预报，与未经优化的CMAQ模型相比，采用极端随机树方法并结合WRF预报气象要素的CMAQ-MOS模型在相关系数、均方根误差、平均相对误差、平均相对偏差方面都有较为明显的优化效果，NO₂及O₃两种污染物优化效果最为明显，其原因可能是O₃、NO₂与MOS模型引入的温度气象要素的相关性较高，这是因为温度的变化往往与太阳辐射的变化相关，而太阳辐射的变化会导致O₃光化学反应速率的变化。与此同时，NO₂和O₃的化学耦合较强，O₃浓度的变化会导致NO₂浓度也发生变化，因此MOS模型针对NO₂及O₃两种污染物优化效果最为明显。

综上所述，采用极端随机树方法并结合WRF预报气象要素的CMAQ-MOS模型预报准确度得到较大的提升，未来可以投入实际业务应用。

此外，采用极端随机树方法并结合WRF预报气象要素的CMAQ-MOS模型仍存在有一定的局限性，有待进一步改进，主要包括以下两个方面：(1)本方法优化的仅为点位数据，并不能完成对整个预报场的优化；(2)本方法针对空气质量模型预报较差时优化效果较好，但针对空气质量模型预报较好时优化效果并不明显。