GRAPES-CUACE伴随模式构建及其在北京地区黑碳和臭氧溯源分析中的应用

安兴琴，王 超，金 敏，李江涛

（1.中国气象科学研究院 大气成分研究所，北京 100081；2.复旦大学 大气与海洋科学系，上海 200438；3.武汉市气象局，武汉 430000）

在大气环境研究领域，敏感性分析受到越来越多的关注。常用的敏感性分析方法（如有限差分法）本质上都是基于正向模拟，需要多次运行大气化学模式，计算代价高，只适合输入变量较少的情况。伴随方法采用反向模拟的思想，通过建立大气化学模式的伴随模式，只需要一次伴随模拟，即可得到目标函数关于所有输入参数的敏感性，其计算效率是正向敏感性分析方法所不可比拟的[1]。

伴随方法提供的敏感性信息可以应用于多种非线性优化控制问题中，如最优污染源控制策略、模型参数估计、污染源反演和大气化学资料同化等[1-6]。这些非线性优化控制问题的求解途径主要为：构建合理的目标函数（如源反演问题的目标函数包含了观测浓度与模拟浓度的差异），利用伴随模式计算目标函数关于模式输入参数（如排放源）的梯度（敏感性），将此梯度信息作为优化算法中的下降步长，更新迭代获得新的输入参数值，如此反复迭代，直到目标函数达到极小值，此时得到的即为更加合理的输入参数值[3-5]。

伴随方法最早在20世纪70年代被引入大气环境领域[7-8]，经过数十年的发展，逐渐成为大气环境科学研究的重要工具[9-11]。目前，国内外不少学者基于国际上应用广泛的大气化学模式开发了相应的伴随模式，如美国学者Henze等[12]开发了全球化学传输模式GEOS-Chem的伴随模式，Sandu等[13]开发了空气质量模式STEM-Ⅲ的伴随模式，德国学者Elbern等[14]构建了欧洲空气污染扩散模式EURAD的伴随模式，国内学者刘峰等[3]建立了大气化学传输模式CAMx的伴随模式等。

GRAPES-CUACE是在中国研发的数值模式系统GRAPES（global-regional assimilation and prediction system）[15]和气溶胶模式CUACE（CMA unified atmospheric chemistry environmental forecasting system）[16]基础上开发的在线耦合的大气化学模式[17]，本文开发了GRAPES-CUACE气溶胶模块和气体模块的伴随模式[18-19]，为求解大气排放源与浓度之间的敏感性关系及下一步的污染源反演工作奠定基础。

1 GRAPES-CUACE大气化学模式介绍

1.1 GRAPES-CUACE气溶胶模块

GRAPES-CUACE气溶胶模块中涉及了六个物种的粒子，包括：黑碳、有机碳、海盐、硫酸盐、硝酸盐和沙尘。气溶胶模块将粒子粒径分成了12个粒径段[16-17]，粒子在各粒径段的质量守恒方程为：

其中，p为干颗粒要素，i为粒径范围。该方程描述了气溶胶在大气中的主要过程，如传输、排放、晴空过程、干沉降以及气溶胶-云相互作用等。此外，气溶胶模块也考虑了气溶胶的垂直扩散以及二次气溶胶的形成等机制。

1.2 GRAPES-CUACE气体模块

气体模块采用CBM-IV（carbon bond mechanism IV）化学机制，共包括36种化学物种和94个化学反应，根据有机化合物官能团及反应活性的不同，碳键共分为12类，即C—C饱和键（PAR）、乙烯（ETH）、C=C不饱和键（OLE）、甲苯（TOL）、二甲苯（XYL）、异戊二烯（ISOP）、甲醛（FORM）、乙醛（ALD2）、α-羰基醛（MGLY）、甲酚（CRES）、过氧乙酰硝酸酯（PAN）、有机硝酸酯（NTR）[20]。气体模块考虑了气体排放源及其垂直扩散，其核心是VOCs与NOx及O3之间一系列复杂的化学反应过程[19，21]。

2 GRAPES-CUACE伴随模式构建

2.1 伴随理论

假设W和G是两个不同的Hilbert空间，L是从G到W的连续线性算子，如果存在唯一一个从W到G的连续线性算子L*满足[1]：

则称L*是L伴随算子[1，3]。当L为矩阵型线性算子时，有(LX,Y)=YTLX=XTLTY=(X,LTY)，即矩阵型线性算子的伴随算子就是它的转置，即L*=LT；若矩阵的分量是复数，则伴随算子L*是其转置的共轭[1，3]。

2.2 伴随模式的构建

构造原模式的伴随模式时，可先将原模式分为多个小段程序，再将一小段程序抽象为一个向量函数F∶Rn→Rm[1，18]，即：

其中，X和Y分别为n维独立变量和m维依赖变量，对应原程序中的输入变量和输出变量[1，18]。假设F在给定点X0是连续可微的，则其切线性模式[1，18]为：

根据伴随理论式（2），（4）的伴随模式为：

其中dX*和dY*分别为n维和m维。对比式（4）和式（5）可知，两式交换了输入输出的维度，算子取了转置[1，18]。

这里以m=1，n=2为例，构建其对应的切线性语句和伴随语句。例如：

原语句：Y=X1·X2，式中X1、X2为输入变量，Y为输出变量，X1、X2和Y各自独立。

切线性语句：dY=X2-dX1+X1dX2，dX1、dX2、dY分别为X1、X2、Y的扰动态向量。

其中，X1*、X2*和Y*分别为X1、X2和Y的伴随变量。

由此可见，为计算输入变量的梯度，切线性程序需要运行两次，而伴随程序通过令Y*=1.0只需计算一次。切线性方法的计算代价随着输入变量的增多而增大，当m＜＜n时，伴随方法的计算效率是切线性方法所不可比拟的[1，18]。

2.3 伴随模式验证

构建完成后的伴随模式需要进行正确性测试以验证伴随代码的正确性。首先需要验证切线性代码的正确性，在此基础上才能对伴随代码进行正确性测试。

2.3.1 切线性验证

（1）验证方法

直接测试整体代码很难确定出错位置，为此可采取分段测试的方式，待局部代码通过正确性测试后再耦合起来测试整体程序[1，18]。假设每一小段代码可以视为Y=F(X),F(X+δX)在X点的泰勒展开式[1，18]为：

其中，分子为原模式输入参数加上扰动后的模式输出增量，分母为切线性模式的输出值。通过等比例减小δX的值，重复计算式（6），若式（6）右侧计算结果趋近于1.0，则说明构建的切线性代码正确[1，18]。

（2）验证结果

模型中各输入变量均能通过切线性模型验证。模型中输入变量较多，由于篇幅原因，这里只给出代表性变量（污染物浓度xrow以及粒子湿半径rhop）的验证结果，如表1所示。其中，扰动值δX=0.001，每次缩减比例a=0.1。可以看出，Index值随着a的减小而逐渐趋近于1.0，但由于机器舍入误差的存在，当a小到一定程度时，Index值又慢慢远离1.0，呈抛物线形[1，18]。

表1 切线性检验结果Table 1 Validation results of the tangent linear model

2.3.2 伴随验证

在验证切线性代码正确性的基础上，测试伴随代码的正确性。同样对伴随代码采取分段处理的方式。

（1）验证方法

伴随代码和切线性代码之间需要满足式（2）的关系：

(LX,Y)=(X,L*Y)

其中，L为切线性过程；L*为对应的伴随过程。令Y=L(X)，得到：

将dX作为切线性模式的输入项，得到切线性输出值∇F·dX，随后将切线性输出值∇F·dX作为伴随模式的输入项，得到伴随输出值∇TF(∇F·dX)，进而分别计算（7）式的左右端[1，18]。若左右两端在机器误差范围内近似相等，则认为构建的伴随代码是正确的。

（2）验证结果

表2以输入变量污染物浓度(xrow)为例，给出了伴随测试结果。表中，VALTGL为式（7）左端计算结果，VALADJ为式（7）右端计算结果，对比发现，VALTGL和VALADJ数值的一致性超过了14位有效数字。除变量xrow之外，模型中所有输入变量均进行了测试，测试结果验证了伴随代码的准确性及精确度，同时也验证了伴随代码与切线性代码之间的一致性。

表2 伴随检验结果Table 2 Validation results of the adjoint model

3 敏感性实验

目标函数可以定义为模式模拟值与观测值之间的差异，敏感性则体现了目标函数对输入参数（如排放源）扰动的响应大小。本研究分别利用已构建的CUACE气溶胶伴随模式和气体伴随模式，对北京地区黑碳气溶胶（BC）和臭氧（O3）进行数值模拟实验，分析污染物模拟浓度与观测浓度之差对于排放源的敏感性。

3.1 黑碳气溶胶伴随敏感性分析

3.1.1 数据资料及实验设计

南郊观象台（简称南郊站，NJ）和上甸子站（SDZ）分别是北京市的城区和郊区代表站（图1），本文使用的BC观测数据选取了两个站2010年1、4、7和10月的BC浓度小时均值。

BC的模拟实验在箱模式中进行，不考虑气溶胶的传输过程，且固定气象条件不变，以2007年估算出的BC排放源[22]作为先验源驱动大气化学模式CUACE，得到南郊站和上甸子站2010年1、4、7和10月的BC模拟浓度；定义目标函数为BC观测浓度与模拟浓度之差，随后逆时序运行CUACE气溶胶伴随模式，得到目标函数关于排放源的敏感性。

图1 站点位置Fig.1 The locations of observation sites

3.1.2 正向模式模拟结果

图2和图3分别给出了上甸子站和南郊站2010年各季节BC模拟浓度与观测浓度的日变化曲线。可以看出，CUACE模式模拟的北京城区和郊区站点四个季节的BC浓度变化趋势与观测数据较为一致，能够反映BC的日变化过程，但是BC模拟浓度与观测值仍存在一定差异。

图2 2010年上甸子站黑碳气溶胶模拟浓度与观测浓度对比Fig.2 Comparison between simulated and observed concerations of biack carbon aerosol at Shangdianzi station in 2010

图3 2010年南郊站黑碳气溶胶模拟浓度与观测浓度对比Fig.3 Comparison between simulated and observed concentrations of black aerosol at Nanjiao station in 2010

3.1.3 敏感性分析

图4与图5分别给出了上甸子站和南郊站目标函数（观测浓度与模拟浓度差值）关于BC排放源的敏感性日变化曲线。可以看出，各图中敏感性曲线与浓度差值曲线的变化趋势较为一致，两者之间成正比关系。这是由于BC不参与大气中的其他化学反应，且在气象条件一定以及箱模式不与外界发生物质交换的情况下，模拟的BC浓度与输入的BC源排放正相关。

敏感性数值的符号反映了模拟浓度相对观测浓度的偏离情况，敏感性数值为正（负），说明模拟浓度较观测浓度高（低），此时的BC排放源高（低）于实际排放源强。从图4、图5可以看出，上甸子站和南郊站各季节的敏感性值均为负值，说明两个站点各季节模拟的BC浓度较观测浓度低，这与图2、图3中的模拟浓度与观测浓度的对比情况是一致的。

敏感性绝对值的大小反映了模式模拟的准确度，敏感性绝对值越小（大），说明模拟浓度与观测浓度差异越小（大），此时的BC排放源与实际排放情况差异越小（大）。可以看出，上甸子站（图4）和南郊站（图5）各季节敏感性较大绝对值出现的时刻主要为夜间和傍晚，说明夜间和傍晚的BC模拟浓度与观测浓度的差异较其他时刻大，这与图2、图3中的模拟浓度与观测浓度对比情况是一致的。

敏感性的符号和绝对值大小，可以为BC排放源的反演提供初步方向。例如，南郊站（图5）1月00:30的敏感性数值为-8×10-11，说明此时需要对BC排放源强进行较大幅度的增加；而上甸子站（图4）7月12:30的敏感性数值为-2×10-12，说明此时只需要对BC排放源强进行小幅度增加。更进一步的源反演工作需要结合最优化算法，以敏感性值为依据确定下降方向，对排放源强不断更新迭代，直到目标函数取得极小值，此时得到的即为排放源强的最优解。

图4 2010年上甸子站黑碳气溶胶目标函数（观测浓度与模拟浓度差值）关于排放量的敏感性日变化Fig.4 Diurnal variation of the sensitivity for the objective function(difference between observed and simulated concentrations)of black carbon aerosol with respect to emissions at Shangdianzi station in 2010

图5 2010年南郊站黑碳气溶胶目标函数（观测浓度与模拟浓度差值）关于排放量的敏感性日变化Fig.5 Diurnal variation of the sensitivity for the objective function(difference between observed and simulated concentrations)of black carbon aerosol with respect to emissions at Nanjiao station in 2010

3.2 O3伴随敏感性分析

3.2.1 数据资料及实验设计

图6 2015年7月7日08:00—8日07:00顺义新城站O3模拟浓度与观测浓度对比Fig.6 Comparison of simulated and observed O3concentration at Shunyi station during 08:00 7th July 2015 to 07:00 8th July 2015

观测数据选取北京市顺义新城监测站（图1）2015年7月7日08:00至8日07:00的O3浓度小时值。通过运行GRAPES-CUACE正向模式和伴随模式，分别模拟得到顺义新城站2015年7月7日08:00至8日07:00的O3浓度以及目标函数（观测浓度与模拟浓度差值）关于NOx和VOCs排放源的敏感性。

3.2.2 正向模式模拟结果

2015年7月7日08:00至8日07:00顺义新城站O3模拟浓度与观测浓度的时间序列如图6所示。可以看出O3模拟浓度与观测浓度的变化趋势较为一致，能够反映出O3的日变化过程，但O3的模拟值较观测值偏低。

3.2.3 敏感性分析

图7给出了目标函数（观测浓度与模拟浓度差值）关于NOx和VOCs排放源的敏感性分布。可以看出，NOx和VOCs敏感性高值（绝对值）主要位于顺义、北京西南部的市区、廊坊以及天津，这与当天北京地区受东南风影响有关。在源反演的优化迭代中，通常以负梯度（敏感性）方向作为下降方向，对比NOx和VOCs的敏感性分布可以看出，NOx的敏感性在区域内为正值，而VOCs的敏感性在区域内主要为负值，这表明，若要减小2015年7月7日08:00至8日07:00北京顺义站O3模拟浓度与观测浓度的差异，需要对NOx和VOCs排放源分布进行调整，即在当前状态减小区域内的NOx排放源以及增大对应网格点的VOCs排放源。

图7 2015年7月7日08:00至8日07:00顺义新城站O3目标函数（观测浓度与模拟浓度差值）关于NOx和VOCs排放源的敏感性分布Fig.7 Distribution of the sensitivities of the O3 objective function(difference between observed and simulated concentrations)with respect to NOx and VOCs emissions at Shunyi station during 08:00 7th July 2015 to 07:00 8th July 2015

值得注意的是，由于O3浓度与前体物NOx和VOCs之间复杂的非线性关系，同时敏感性是随着排放源变化而变化的[3]，此处讨论的NOx和VOCs敏感性仅针对当前状态，不能仅仅从这一个状态的敏感性来确定最终的排放源分布，而是需要结合优化算法，利用当前状态的敏感性计算得到一个新的排放源分布，在新的状态下，重新计算敏感性，再进行下一步优化，如此反复迭代，直到满足收敛条件，才能给出更为合理的排放源分布。

4 结论

（1）依据伴随理论，基于GRAPES-CUACE大气化学模式，建立了该模式气溶胶模块和气体模块的伴随模式，并对其进行了正确性测试。结合BC和O3浓度观测数据，分别利用气溶胶和气体伴随模块进行了数值模拟及源-浓度的敏感性实验。

（2）气溶胶伴随模型的模拟结果表明：目标函数关于BC排放源的敏感性日变化趋势与观测浓度和模拟浓度差值的日变化趋势较为一致，两者之间成正比关系。敏感性绝对值的大小反映了模式模拟的准确度，敏感性数值的正（负）则反映了模拟值较观测值偏高（低），可依据敏感性确定下降方向，结合最优化算法，反演排放源强。

（3）气体伴随模型的模拟结果表明：若要减小2015年7月7日08:00至8日07:00北京顺义站O3模拟浓度与观测浓度的差异，需要对NOx和VOCs排放源分布进行调整，即在当前状态减小区域内的NOx排放源以及增大对应网格点的VOCs排放源，在得到新的排放源分布情况下，重新计算敏感性，再进行下一步优化，如此反复迭代，当满足收敛条件时即可给出合理的排放源分布。

（4）本文开发的GRAPES-CUACE气溶胶伴随和气体伴随能够有效地针对气溶胶和气体进行敏感性分析，为下一步的污染源反演工作奠定了基础。未来将继续开发优化算法模块，构建完整的GRAPES-CUACE大气化学四维变分同化系统，为大气污染优化控制提供科学依据。