随着经济的发展，人类活动向大气环境排放的污染物越来越多，化学成分也越加复杂，使得大气污染成为世界上许多城市面临的最严重的环境问题之一.近年来，空气质量问题逐渐成为公众关注的焦点，从而对空气质量进行预报的需求也日益迫切(聂邦胜，2008).空气质量数值模式具备完善的理论基础、合理的模式设计和定时定量的预报结果，在计算机和数值计算技术飞速发展的今天，逐渐成为空气质量预报的主要发展趋势(房小怡等，2004).如今，空气质量数值模式已经发展到第三代，其在模拟大气污染物反应、输送和清除过程方面的能力有了长足发展，但化学过程的复杂性、多尺度耦合的影响及污染物排放源的不确定性等因素仍然给模拟准确度的进一步提高带来了巨大挑战(白晓平等，2007；Carmichael et al., 2008).

提高空气质量数值模式模拟能力的一种途径是利用数据同化方法，把观测资料和数值模式的信息结合起来获得对大气状态的最优估计，在数值模拟中为模式提供较优的初始场(王自发等，2006；Constantinescu et al., 2007)，从而最终改善模拟质量.随着越来越多可用的大气化学观测资料的出现，数据同化在空气质量数值模式中的应用前景也越来越广阔(Carmichael et al., 2008).

数据同化技术早在20 世纪60 年代就在气象领域得到应用(Charney et al., 1969)，随后发展到海洋领域，而针对大气化学资料的同化研究则是从20世纪90年代才开始的(Austin，1992).目前国际流行的同化方法主要是四维变分同化(4DVAR)和集合卡尔曼滤波(EnKF)(王自发等，2006)，但较难的编程实现和过大的计算量分别使这两种方法的应用受到限制，Evensen(2003)提出集合最优插值(EnOI)可作为EnKF的次优方案，能够在大大节省计算量的同时保持一定的有效性.万莉颖(2006)在海洋模式HYCOM的同化试验中对EnKF与EnOI的同化结果进行比较，得出了EnKF相对EnOI改变不是非常明显的结论.因而在计算能力受到制约的情况下，以EnOI进行同化试验是可以接受的替代方案.

与其他主流同化方法相比，国际上针对EnOI同化方法的应用研究还相对较少，且多集中在海洋模式中(Oke et al., 2005；2010；2002；Counillon et al., 2009)；国内的相关研究同样限于海洋领域，更是尚处于初步探讨阶段(万莉颖，2006；叶冬等，2011；闫长香等，2011；万莉颖等，2011；Wan et al., 2010)，研究还不够深入；而国内对EnOI在空气质量数值模式中的应用研究就更加缺乏.基于此研究现状，本文旨在探讨EnOI同化方法在珠三角地区空气质量数值模拟中应用的可行性，着重研究EnOI对MM5-STEM(Wang et al., 2005)模式污染物浓度预报场的修正能力，以期为集合同化方法在珠三角空气质量数值模拟中的应用奠定基础.

背景误差协方差矩阵(下称“B矩阵”)的确定是数据同化中的关键.同样是利用集合样本来构造B矩阵，EnOI与EnKF的最显著区别在于：前者使用由历史样本提供的静态样本构造B矩阵，B 矩阵不随模式积分更新，是静止的；而后者的B矩阵通常是通过将扰动叠加在某一给定模式状态的方法来构造(闫长香等，2011)，随着模式积分不断更新，是动态的.除此之外，两者的其余计算过程几乎一致.但与EnKF的实施需要并行多个算例不同，EnOI在同化过程中仅需运行一个算例，在计算资源需求上远小于EnKF(闫长香等，2011)，并且能较好地避免EnKF中潜在的滤波发散风险(叶冬等，2011).

EnOI同化更新公式(Evensen，2003)如下：

通常，B 矩阵以表示，而观测误差协方差矩阵(下称“R矩阵”)以表示.本研究中的同化试验均不通过扰动方式产生R矩阵，而是假定R矩阵为对角阵，观测误差为观测值的10%(崔应杰等，2006)，以观测误差平方作为R矩阵的对角线元素；同时假设不同站点的观测不相关，即R矩阵的非对角线元素设置为零.公式(1)中的γγ^T以N-1R代替(叶冬等，2011).

由于集合样本数量的有限性，原始的B矩阵会包含许多虚假相关，容易造成同化结果的不合理.为避免这一现象，本研究采用了高斯分布函数(参看公式(6))对格点间的相关性作局地化约束，规定了其随距离增大而减小的特征，其中的r为格点或站点间的距离；L为局地化尺度或水平相关尺度，L越大则表示同化站点能够影响到的范围越大.

最终，公式(1)改写成以下形式：

本研究所采用的MM5-STEM空气质量模型系统包括气象条件模拟及化学传输、转化模拟等模块，其中，多尺度、多界面空气质量模型系统STEM(Sulfur Transport and dEposition Model)是模式系统的核心模块，它可以模拟气态、气溶胶等多种污染物的输送和转化过程(Tang et al., 2003a；2003b；Carmichael et al., 2003a；2003b).STEM已经成功应用于珠江三角洲及东亚地区污染物输送和转化特征的研究(Wang et al., 2005a；2005b；Wang et al., 2007；Liu et al., 2013)，并在多项国际大型野外观测计划(TRACE-P)中被用作预报模式，包括光化学反应、气溶胶、沉降等，所用化学反应机理为SAPRC99(Carter，2000)，通过在线TUV模块进行光解率的计算(Madronich et al., 1999)，能够较好地模拟区域空气质量的演变.

本研究所采用的污染源排放清单来自文献(Zheng et al., 2009；郑君瑜等，2009a；2009b；余宇帆等，2011；车汶蔚等，2009；尹沙沙，2011)，观测资料来源于广东省环境监测中心及广州市环境监测中心，包括珠三角地区合共57个污染物观测站点2010年11月的NO₂、SO₂及PM₁₀日均浓度值数据.57个站点中的13个为粤港珠三角区域空气监控网络(张远航等，2008)站点(见图 1同化站点)，分别为广州麓湖(LH)、广州万顷沙(WQS)、广州天湖(TH)、深圳荔园(LY)、惠州金果湾(JGW)、惠州下埔(XB)、东莞豪岗小学(HG)、佛山惠景城(HJC)、佛山金桔咀(JJJ)、肇庆城中(CZ)、中山紫马岭(ZML)、珠海唐家(TJ)、江门东湖(DH)；其余44个站点为珠三角城市省控站点(见图 1检验站点).

图 1(Fig. 1)

图 1 粤港珠三角区域空气监控网络站点(同化站点)及珠三角省控站点(检验站点)分布 Fig. 1 Distribution of PRD Regional Air Quality Monitoring Network(assimilation sites) and Guangdong provincial monitoring sites in the PRD(validation sites)

本研究旨在探讨集合最优插值法在珠三角空气质量数值模拟中修正模式预报场的效果，了解该方法对MM5-STEM模式中不同污染物预报场的修正能力.试验方式为对2010年11月30天的NO₂、SO₂和PM₁₀日均浓度值进行逐日同化，理论上可以用分析误差矩阵来检验同化效果，但直接求解分析误差矩阵非常困难，本研究将通过比较同化前后模拟值与观测值的均方差误差(RMSE)来检验同化效果(崔应杰等，2006).

对于进入同化系统观测站点的选择问题，考虑到粤港珠三角区域空气监控网络(张远航等，2008)在设计站点布局时已经得到一定的科学论证，该监控网络的站点也已在长期的业务运行中被用以代表整个珠三角的空气质量整体水平并一直受到严格的质控管理，综合数据代表性及数据质量两方面的因素，本研究直接以粤港珠三角区域空气监控网络位于广东省境内的13个站点作为进入同化系统的观测站点(下称“同化站点”，见图 1)，而把珠三角区域内的其他44个省控站点作为独立检验站点(下称“检验站点”，见图 1).大量检验站点的设置可以避免个别检验站点因代表性不足导致的对整体同化效果的错误分析，这也是本研究区别于前人研究(崔应杰等，2006)的一个特点.

本小节将简要展示MM5-STEM在本研究中的模拟情况，关于MM5-STEM的模拟性能及其在珠三角地区应用的更多细节，可参考前人的研究(Wang et al., 2005a； 2005b；Wang et al., 2007; Liu et al., 2013).

表 1 列出了2010年11月30天的NO₂、SO₂和PM₁₀平均浓度在全部57个站点的模拟效果的统计分析.从观测均值、模拟均值及MNB均值上看，NO₂、SO₂和PM₁₀的模拟值较观测值均有不同程度的偏低，其中，NO₂和SO₂的偏离程度相对较小，其MNB均值分别为-18%和-22%，PM₁₀则偏低较多，MNB均值达-70%.MNB的最小最大值体现了各污染物在不同站点模拟值偏离观测值的极端情况，可见NO₂和SO₂的模拟值同时存在偏高和偏低的现象，而PM₁₀则表现出整体偏低，表明PM₁₀的源清单可能被较大程度地低估.从FAC2值来看，NO₂和SO₂有超过或接近一半的模拟值与观测值的比例在0.5~2之间，而PM₁₀则只有15%满足此条件.

表1(Table.1)

表1 NO2、SO2和PM10日均浓度模拟效果统计分析 Table.1 Statistical analysis of NO2，SO2 and PM10 daily mean concentration simulation performance

表1 NO2、SO2和PM10日均浓度模拟效果统计分析 Table.1 Statistical analysis of NO2，SO2 and PM10 daily mean concentration simulation performance 污染物 平均值/(μg·m-3) MNB FAC2 观测 模拟 AVE MIN MAX 注：观测平均值，模拟平均值，平均标准偏差; 均方根误差，FAC2为满足的数据比例，其中，Oi为观测值，Mi为模拟值. NO2 55 39 -18% -93% 286% 59% SO2 30 20 -22% -89% 270% 47% PM10 87 24 -70% -96% -29% 15%

考虑到本研究所采用的污染源排放清单在广州地区的精细化程度优于广州以外的珠三角地区，参照文献(Wu et al., 2006)的做法，将广州地区及广州以外珠三角地区的站点分别进行评估，其中，广州市10个国控点为一组(以GZ表示)，粤港珠三角区域空气监控网络中广州以外珠三角地区的9个站点(不含区域站)为另一组(以PRD表示).图 2以泰勒图(Taylor，2001)的形式给出了综合了相关系数、标准差、中心均方根误差3个统计参数的评估结果，其中，OBS为参考点，越接近参考点(即相关系数与标准差接近1，中心均方根误差RMSD接近0)时模拟效果越好.从图可知，GZ站点中NO₂及SO₂的大部分站点均位于相关系数0.3~0.8、标准差0.5~1及中心均方根误差0.5~1的范围内，模拟相对较好；而GZ站点中PM₁₀则集中在相关系数 -0.2~0.3，标准差0~0.5及中心均方根误差为1左右的范围内，模拟相对较差.PRD站点在相关系数及中心均方根误差方面均与GZ站点差异不大，但标准差普遍小于GZ站点，表明总体上GZ站点模拟效果优于PRD站点.

图 2(Fig. 2)

图 2 广州及广州以外珠三角地区的NO2、SO2和PM10日均浓度模拟泰勒图(黑线代表归一化标准差，蓝线代表相关系数，绿线代表归一化的中心均方根误差) Fig. 2 Taylor plot of NO2，SO2 and PM10 daily mean concentration at Guangzhou stations and other PRD cities stations(dark line denotes normalized st and ard deviation，blue line denotes correlation coefficient， and green line denotes normalized centered root mean square error)

综上所述，本研究中SO₂和NO₂的模拟效果相对较好，而PM₁₀偏低较多；广州地区站点的模拟效果普遍优于广州以外珠三角地区站点.

局地化尺度L决定着同化站点的影响范围，在大气污染物观测站点分布密度较低情况下，若局地化尺度L设置得太小，大部分区域将得不到更新，也达不到对数据的有效利用；若L设置得太大，每个观测点都在较大的区域内产生影响，则可能会造成不合理的分析增量，反而降低同化效果.因而可能存在一个合适的L，在这个尺度下同化效果将达到最优.另外，从公式(7)来看，经验系数α的不同取值也将对同化效果造成一定影响，本小节将探讨不同的L和α组合对同化效果的影响.

本研究循环测试了在不同α(2^-3~2³，间隔为2倍)和L(10~100 km，间隔10 km)设置下的NO₂、SO₂和PM₁₀的30天逐日同化实验，分别考查13个同化站点及44个检验站点的平均RMSE值以检验同化效果.不同污染物的集合样本分别由相应污染物2010年1、4、7、11月4个月的120个日均浓度模拟值组成，以代表全年各个季节的状态，EnOI的同化计算过程参照公式(1)~(7).

经过共70个算例的试验，求出使各污染物检验站点平均RMSE最小的L和α组合，定义为“最优L和α组合”，结果列于表 2.可见不同污染物由于其本身浓度分布特征及模拟效果的差异，其最优L和α组合也不一样，结合其他时段(1、4、7月)的试验结果，总体而言有如下规律：①相对于检验站点，同化站点本身受同化影响更大，且对α的敏感性比L大，一般情况下，无论L取何值，α取值越大，则同化站点的模拟值越向观测值靠近，RMSE值越小；②对于检验站点，当α取值较大时往往存在一个最优的L值，L取此值时RMSE值最小；③对于模拟值整体偏离观测值较多的污染物(如本研究中的PM₁₀)，在一定的L范围内，α及L的增大均有利于检验站点RMSE值的减小；④无论是同化站点还是检验站点，当L足够大时，不同污染物的RMSE值均随L的增大而逐渐稳定.

表2(Table.2)

表2 各污染物的最优L和α组合 Table.2 Optimal settings of L and α of different pollutants

表2 各污染物的最优L和α组合 Table.2 Optimal settings of L and α of different pollutants 污染物 L/km α NO2 100 0.5 SO2 20 4 PM10 40 8

应该指出的是，本小节中所称的“最优L和α组合”只是在特定试验设置下的计算结果，在不同的样本来源及试验时段中可能获得不同结果，但由此获得的规律是普遍适用的.

为探讨EnOI对不同污染物的同化能力，本节以4.2节中NO₂、SO₂和PM₁₀各自的“最优L和α组合”方案的计算结果为例，从同化站点及检验站点同化前后RMSE值对比、同化前后月均浓度分布及同化分析增量与模式误差对比等角度进行讨论分析.

由表 3可知，各污染物在所有同化站点经过同 化后RMSE值均有不同程度的下降，RMSE下降站点比例均达100%;RMSE平均降幅的排序为PM₁₀>SO₂>NO₂，平均下降比例的排序为SO₂>PM₁₀>NO₂.由于PM₁₀同化前的RMSE均值远大于NO₂及SO₂，故其同化后的RMSE平均降幅也最大.而虽然在两次排序中SO₂均比NO₂靠前，但若把L和α均设置在同一水平(L＝50 km，α＝1)，SO₂和NO₂的位置将对调(图表略)，表明对于同化前RMSE均值相当的不同污染物，同化站点的RMSE下降情况受α值的影响较大，这与4.2节中第1个结论是一致的，即对于同化站点本身，其同化效果对α的取值较为敏感.

表3(Table.3)

表3 同化站点及检验站点同化前后RMSE值对比 Table.3 Comparison of RMSE before and after assimilation at assimilation sites and validation sites

表3 同化站点及检验站点同化前后RMSE值对比 Table.3 Comparison of RMSE before and after assimilation at assimilation sites and validation sites 数据类型 单位 NO2 SO2 PM10 同化站 检验站 同化站 检验站 同化站 检验站 同化前RMSE均值 μg·m-3 32 30 25 23 80 66 同化后RMSE均值 μg·m-3 16 20 5 16 26 38 RMSE平均降幅 μg·m-3 17 10 20 7 54 28 RMSE平均下降比例   52% 33% 80% 32% 68% 42% RMSE最小降幅 μg·m-3 2 -25 5 -13 38 -22 RMSE最大降幅 μg·m-3 32 32 35 27 92 59 RMSE最小下降比例   12% -215% 70% -107% 43% -48% RMSE最大下降比例   98% 77% 98% 83% 94% 83% RMSE下降站点比例   100% 86% 100% 84% 100% 91%

图 3展示了NO₂在各个同化站点同化前后的RMSE对比情况.以RMSE下降比例最高的广州天湖站(98%，图 3中的TH)和最低的中山紫马岭站(12%，图 3中的ZML)为例进行分析.因广州天湖站位于氮氧化物源排放相对较少的清洁区，且试验时段(11月份)外来源影响微弱，故NO₂观测值较低(平均仅10 μg · m^-3)，按照取观测误差为观测值10%的假设，该站点被认为具有极小的观测误差，即同化过程对观测值的信任程度很大；另一方面，集合样本因包含了不同季节的波动使得计算出的背景误差不可忽略，两种因素的叠加使同化分析值极大程度地向自身观测值靠近.而对于中山紫马岭站，因该站点位于4个同化站点的包围之中，而L又被设置为较大的100 km，其同化结果体现了相邻同化站点的共同作用，导致同化分析值与自身观测值仍有较大偏差.

图 3(Fig. 3)

图 3 同化站点NO2同化前后RMSE值对比 Fig. 3 Comparison between RMSE before and after assimilation of NO2 at assimilation sites

对于检验站点，不难发现表 3中表征RMSE下降情况的各项指标均小于同化站点，且各污染物均出现了少量同化后RMSE反而上升的站点，占全部检验站点的9%~16%.原因可能是个别受局地污染源影响的检验站点因缺乏区域代表性而表现出与邻近同化站点差异显著的浓度特征，故同化效果不理想.另外，即使考虑了把L和α均设置在同一水平的情况，从RMSE平均降幅、RMSE平均下降比例及RMSE下降站点比例3项指标来看，均保持了PM₁₀>NO₂>SO₂的排序，可能反映了不同污染物扩散能力的差异.

图 4展示了NO₂在各个检验站点同化前后的RMSE对比情况.在44个检验站点中，除广州、佛山及肇庆的6个站点表现出同化后RMSE值反而上升 的趋势外，其余站点的RMSE值均有不同程度的下降，占检验站点数的86%，所有检验站点平均下降比例为32%，整体同化效果较好.其中，最大RMSE下降幅度(32 μg · m^-3)及最大RMSE下降比例(77%)均出现在肇庆.图 5给出了检验站点同化前后NO₂模拟均值与观测均值的时间序列，可见经过同化修正的模拟值已与观测值相当接近.

图 4(Fig. 4)

图 4 检验站点NO2同化前后RMSE值对比 Fig. 4 Comparison between RMSE before and after assimilation of NO2 at validation sites

图 5(Fig. 5)

图 5 检验站点同化前后NO2模拟均值与观测均值时间序列(2010年) Fig. 5 Time series of average observation and simulation before and after assimilation at validation sites(2010)

从图 6可以直观地对比同化前后的NO₂浓度分布差异，同化前的NO₂高值区呈零散的斑点状分布，主要对应于排放源强度较大的点源；同化后的NO₂高值区把广州城区及佛山中部连在一起，呈带状分布，更真实地反映了汽车保有量更大的城市中心区的NO₂浓度分布特征，表明EnOI能够有效地结合观测值和模拟值，提供更符合实际情况的污染物浓度预报场.

图 6(Fig. 6)

图 6 NO2同化前后月均浓度分布图(a：同化前；b：同化后) Fig. 6 Monthly average concentration field of NO2 before and after assimilation(a: before assimilation; b: after assimilation)

为进一步了解同化前后整个区域的浓度变化细节及同化效果，用图 6中同化后模拟值减去同化前模拟值，得出模拟区域内平均分析增量，附以各检验站点的平均模拟误差进行对比，具体如图 7所示.其中，平均分析增量以背景颜色表示，平均模拟误差(观测值-同化前模拟值)以方格颜色表示，当代表某站点的方格颜色与该站点附近的背景颜色较接近时，可认为同化产生的分析增量较好地填补了模拟误差，即同化效果较好.

图 7(Fig. 7)

图 7 检验站点NO2同化分析增量与模式误差对比 Fig. 7 Comparison between assimilation analysis increment and model error of NO2 at validation sites

观察图 7可得，由于L被设置为较大的100 km，分析增量超过±10 μg · m^-3的覆盖范围也较大，包括肇庆、佛山、广州、江门、东莞、惠州和深圳的大部分区域，其中，肇庆南部、佛山北部及广州中部片区的分析增量超过了20 μg · m^-3，江门部分地区的分析增量则超过了-20 μg · m^-3.总体而言，大多数站点的分析增量均与其模拟误差有较好的吻合，尤其是肇庆地区的模拟低估与江门地区的模拟高估都得到了较好的修正.

SO₂的同化效果总体情况与NO₂类似，部分特征略有区别.84%的SO₂检验站点表现出同化后RMSE下降，平均下降比例为32%，但不同站点间的同化效果差异小于NO₂.从图 8可以看出，分析增量绝对值的较大值仅分布在肇庆、江门及广州与东莞交界处，量值超过10 μg · m^-3的区域也仅在部分同化站点附近出现.这是由于L被设置为较小的20 km，以致同化站点的影响范围较小；同时因为SO₂在许多同化站点中均模拟较好，故产生的分析增量不大.除了佛山的几个站点，其余大多数站点的分析增量较好地修正了各自的模拟误差.

图 8(Fig. 8)

图 8 NO2检验站点SO2同化分析增量与模式误差对比 Fig. 8 Comparison between assimilation analysis increment and model error of SO2 at validation sites

由模式性能评估一节可知，相对于SO₂和NO₂而言，PM₁₀的模拟值较观测值偏低较多，理论上当模拟值误差较大时应先通过其他手段(如调整源清单)使其达到较合理范围才适宜进行同化处理，本研究在未对PM₁₀源清单进行调整的前提下进行同化试验，可看作对EnOI处理极端情况能力的考察.结果显示，有91%的检验站点PM₁₀表现出同化后RMSE下降，平均下降比例为42%.虽然对PM₁₀的同化设置了最大的α，但PM₁₀同化站点的RMSE值下降比例整体上并不及SO₂，这是由于PM₁₀的观测值较大，从而人为假设的观测误差较大，导致同化结果较不信任观测值的缘故.图 9中PM₁₀的分析增量基本为正，增量较大值基本覆盖了除江门以外的其它珠三角区域，表明在L=40 km的设置下珠三角区域的大部分地区都能够得到观测值的同化，但也同时也反映了模式PM₁₀源清单存在整体被低估的可能.虽然从数据上看，经EnOI同化后的PM₁₀预报场较同化前更加接近观测值，但由于PM₁₀原始模拟值过低，可能会导致EnOI对PM₁₀预报场的过度修正，从而破坏了模式的平衡，对最终的模拟结果产生负面影响.

图 9(Fig. 9)

图 9 检验站点PM10同化分析增量与模式误差对比 Fig. 9 Comparison between assimilation analysis increment and model error of PM10 at validation sites

由于数据的限制，本研究仅进行了二维(地面)的同化试验，也未对近年来公众比较关注O₃和PM_2.5进行EnOI的同化应用研究.进一步的工作将考虑把O₃和PM_2.5列入同化研究对象，并尝试在模式动态运行中探讨EnOI的同化效果.另外，因EnOI采用静态样本而存在天然缺陷，只适合于在计算量受到限制的情况下作为EnKF的次优方法实施，EnKF与4DVAR等方法仍然是数据同化的主流方向，故本研究可视作EnKF在珠三角空气质量模拟中开展的基础工作.

本研究就集合最优插值同化方法(EnOI)在珠三角空气质量模拟中修正不同污染物预报场的应用进行了探讨，结果发现：EnOI在NO₂、SO₂及PM₁₀的同化试验中均取得了较好的效果，检验站点均方根误差(RMSE)的平均下降比例分别可达33%、32%和42%，RMSE值下降的站点占检验站点总数的比例分别为86%、84%和91%.表明该方法能够有效地应用在珠三角空气质量模拟中，为模式提供与实际更为接近的预报场.

致谢： 华南理工大学郑君瑜教授为本文提供了珠三角网格化排放源清单，广东省环境监测中心为本文提供了粤港珠三角区域空气监控网络数据，在此一并表示感谢.