近年来，大气污染模式的研究得到了很大发展。污染排放源信息对模式至关重要。然而，建立污染源清单是一个具有较大不确定性的复杂过程，涉及对经济发展状况、人口密度、能源结构、工厂分布、交通流量、车型比例等多方面的调查研究。因此，根据模式的模拟值和监测站的污染物浓度观测值来反演污染源信息对大气污染预报及控制均有重要意义。

大气污染预报是已知原因求结果，属于正问题，即给定了微分方程和初、边值条件，求满足给定条件的解；而大气污染控制是由结果求原因，属于反问题，即排放源等参数作为未知量出现在微分方程中，我们需利用观测值等附加信息反求未知参数(刘峰等，2003；杨一帆和张凯山，2013)。多种方法已被应用于反演污染源，比如逆向轨迹法(蔡旭晖等，2002)、遗传算法(陈军明等，2002)、伴随方程法(刘峰等，2003；郭少冬等，2010)、四维变分法(Mendoza-Dominguez and Russell，2001)、集合卡尔曼滤波法和集合卡尔曼平滑法(朱江和汪萍，2006；Tang et al.，2013)以及蒙特卡洛空间反演法(杨一帆和张凯山，2013)等。

目前，源反演普遍采用统计类方法，但其存在对观测误差、源清单先验估计误差敏感等弱点。而基于打靶法思想的算法恰好能弥补统计法的不足。打靶法思想的本质是将微分方程两时刻点间的演化问题转化为控制问题来求解，通过调整控制参量使观测时刻点的模拟值接近观测值。在源反演问题中，这个控制参量就是源排放速率。基于打靶法思想的各种算法以其精度高、程序简单、实用性强的特点被广泛应用于系统控制领域(向明华和刘云青，1990；屈香菊，1992；范慕辉和石铁君，1994；凌复华和李继跃，1988；Ning et al.，2004；吴文海等，2005；张洪波等，2008；齐笳羽等，2011)。

本文针对大气污染源反演问题提出一个基于打靶法思想的自适应算法。该算法在精度高、算法简明的基础上还弥补了许多统计方法的不足：能处理源清单中的大误差、初值大误差、观测值在个别时间点的大误差；无需先验分布假设及误差估计。本文还以简单模型的理想试验为例，给出了该自适应算法的计算效果。

从污染物浓度观测值估计污染物排放速率可作为系统控制问题进行考虑。其中，污染物排放速率为控制变量。关于打靶法思想介绍，见附录A。基于打靶法的思路，我们可以将观测值作为靶子，设法寻找使模拟值与观测值的距离小于容忍度的排放速率。

该算法设计的困难主要来自以下两方面：第一，直观上，当模拟值大于观测值时，我们需减小排放速率；当模拟值小于观测值时，我们需增大排放速率。那么如何增减能使射击次数尽量小且程序简明；第二，除了源清单的误差，实际中还不可避免地会遇到初值误差、观测误差。当这些误差较大时，算法能否及时纠错，使模拟值尽快主动回到正常位置(而不是被动地等待长时间模式预热，即Spin Up)。增强算法的健壮性，使算法具有纠错功能的必要性说明见附录B。

因此，采用何种打靶策略能够立刻对上述可能出现的大误差进行纠错并且较为快速、准确地进行源反演是本算法设计要解决的核心问题。这里，进行源反演可利用的唯一信息是误差(E_rr)＝模拟值－观测值。如何利用误差E_rr对排放速 率进行高效校正是十分重要的。判定何时需要对大误差及时纠错以及如何纠错也是要解决的关键问题。

若不考虑模式误差，大气污染物浓度在某一空间点上的抽象模型可以表示为

其中，C是污染物浓度的向量(其分量表示不同污染物种的浓度)；Q是污染源排放速率的向量；算子F表示污染物的输送、化学过程。

假定在t₁，...，t_k时刻有污染物的浓度观测向量 Y₁，...，Y_k，则在时间段［t_i，t_i+1］i=0，1，…，k−1的污染物浓度演化过程可视为打靶问题：在该时间段上控制Q，使得

其中，Tol为设定的允许误差界限。于是，源反演问题转化成了分时间段的打靶问题。Q 在每个时间段上为常值向量。

基于上面的叙述和分析，我们在本节给出各时间段［t_i，t_i+1］上的自适应算法流程。

假设共n个污染物种(即C、F、Q为n维向量)，其中s个物种有观测值(即Y₁，…，Y_k为s维向量)。s维向量Q_coeff 为经验参量，用于计算向量Q在每次射击后的校正幅度。不同的物种可使用不同的分量值。如果Q_coeff 的分量过大，则该物种浓度的模拟值在观测值上下频繁震荡；如果Q_coeff 的分量过小，则该物种浓度的模拟值变化慢，需要较大的射击次数。根据以上现象，我们可以通过数值试验找到恰当的经验参量Q_coeff 。

在每个时间段［t_i，t_i+1］上，算法分为三大步：初始化、若干次射击和保存本时间段的结果。若干次射击又分成4小步，其中步骤(c)(用Q的当前值在［t_i，t_i+1］上解初值问题)即射击。

在每次射击后更新s维误差向量E_rr＝C(t_i+1; Q)−Y_i+1，i=0，1，…，(y−1)(其中C只提取与Y相同的物种分量)，为下一次射击前校正向量Q作准备[见2.2.1节中步骤(d.1)]。

注意到只有通过有效的观测信息，我们才有可能得到有效的源反演值。当初值或观测值在个别时间点存在大误差时，反演信息必然无效。由于观测值在若干个时间点可能存在大误差，这使得系统在自然状态下会经历若干个Spin Up，导致反演信息在Spin Up期间均无效。因此，主动采取措施自动纠错、增强算法的健壮性、使污染物浓度和源反演结果立刻回到正常位置是十分必要的(纠错的必要性图例说明见附录B)。

当初值或观测值在个别时间点存在大误差时，我们的算法恰恰能自动纠错，使污染物浓度和源反演结果迅速回到正常位置。这就大大增强了算法的健壮性，从观测值中尽可能多地反演出有效源排放速率[见2.2.1节中步骤(d.2)]。实际上，为了得到最终的有效信息，我们将对结果进行后处理：去掉纠错阶段的源排放速率(见本文数值试验中图 1e-f’)。

图 1

Fig. 1

图 1 试验1 的结果：（a）浓度C1；（b）浓度C2；（b’）为（b）C2 的局部放大；（c）浓度C3；（d）浓度C4；（e）Q1 源排放速率；（e’）为（e）Q1 去掉纠错阶段；（f）Q 2 源排放速率；（f’）为（f）Q 2 去掉纠错阶段 Fig. 1 Results of Experiment 1: (a) Concentration C1, (b) concentration C2, (b’) partial enlargement of (b), (c) concentration of C3, (d) concentration of C4, (e) emissions rate of Q 1, (e’) panel (e) without error correction phase, (f) emissions rate of Q 2, (f’) panel (f) without error correction phase

(1)初始化Q和Q_coeff；用源清单中的值初始化Q(源清单中的量均为正)；Q的符号标记向量signQ各分量初始化为正；ShootingCount=0；(ShootingCount为射击次数的计数器)。

(2)当s维向量 | E_rr |存在分量大于允许误差界限Tol且ShootingCount小于MaxShooting(最大射击次数)时，步骤分为：(a)ShootingCount＝ ShootingCount＋1；(b)用设定的经验值初始化 Q_coeff；(c)用Q的当前值在［t_i，t_i+1］上解初值问题，得到各物种浓度C(t_i+1; Q)；(d)计算E_rr= C(t_i+1; Q)−Y_i+1，用E_rr更新Q，计算方法如下：

For j = 1到s

(d.1)计算E_rr(j)= C(t_i+1; Q; j)−Y_i+1(j)

(d.2)用E_rr(j)更新Q(j)：

[公式(3)为核心，其他均为了增强算法健壮性、为了纠错而设计]

Switch sign Q(j)

Case sign Q(j)为正：

如果Q(j)＜CheckCritical且E_rr(j)＞Tol，那么(CheckCritical为纠错阈值)

sign Q(j)设为负；[使得下一次能进入Case sign Q(j)为负分支纠错]

Q(j)=－10；[用负值重新初始化Q(j)]

否则

Case sign Q(j)为负：[当初值或观测值在个别时间点远大于真值时，使用负值Q(j)能使污染物浓度和源反演结果立刻纠错，回到正常位置，而不是被动地等待长时间Spin Up]

当E_rr(j)·Q_coeff(j)≤－1时，作：

减小Q_coeff(j)，使得公式(4)中的1+E_rr(j)·Q_coeff(j)＞0。

防止本需要为负的Q(j)由于Q_coeff(j)设定得大而在正、负间来回震荡

Q_coeff(j)＝Q_coeff(j)×0.5;

End of 当 [E_rr(j)·Q_coeff(j)≤－1]

End of Switch sign Q(j)

End of For j =从1到s

(3)保存本时间段［t_i，t_i+1］的计算结果，为下一时间段作准备；保存经过若干次射击最终接受 的Q和C([t_i，t_i+1]; Q)；用C(t_i+1; Q)给下一时间段赋初值。

(1)源清单中各物种的排放速率为非负值。在运行此算法之前，我们需要对源清单进行预处理：如果某物种有观测值且该物种在源清单中某时刻的排放速率为零，则用一个小的正数代替零。这样才能使步骤(d.2)的Q(j)更新公式起作用。

(2)在试验中，我们对各物种使用了相同的允许误差界限Tol。实际上，为了进一步减小射击次数，不同的物种可以使用不同的Tol：值域范围大的物种取值相对大，值域范围小的物种取值相对小。Tol 不宜过小。

(3)2.2.1节中步骤(d.2)的CheckCritical 为设定的较小正值。当源排放速率已经非常小时，若仍有Err(j)＞Tol，那么我们决定使用负值纠错。实际上，这负值在后处理中作为纠错阶段的源排放速率将被去掉。(见图 1e-f’)

在本节中，我们采用下面描述空气质量的简单无量纲模型进行试验：

其中，C_i是第i个物种的浓度，Q₁、Q₂是随时间变化的污染源排放速率。试验在在t∈[0,12]上进行，共3个试验，具体设置见表 1。在试验中，我们假设真实初始浓度C_1true(0) =C_2true(0) =C_3true(0) =C_4true(0) =1，在试验1、试验2中，假设真实排放速率Q_1true(t)=10+ sin(2t)，Q_2true(t)=0.3+0.2cos(t)；在试验3中，假设真实排放速率Q_1true(t)=[1/1.2+sin(t−7.2)](t−7.2)²+3，Q_2true(t)=[2/1.3+cos(t−1)](t−7)²+6。在我们的试验中不考虑模式误差。试验各时间段的射击次数见表 2。

表 1 (Table 1)

表 1 试验列表 Table 1 The detailed setup of the experiments 试验 观测时间间隔 其他设置 试验目的 试验1 1 Tol=1 MaxShooting=10 CheckCritical=0.01 检验本算法对污染物浓度的反 演效果和污染源反演能力 试验2 0.5 Qcoeff=[0.1，0.1] C2(0)=C1(0)=00 C3(0)=C4(0)=2 考察观测时间加密对数值结果的影响 试验3 1 Tol=1 MaxShooting=20 CheckCritical=0.01 Qcoeff=[0.02，0.02] C1(0)=C2(0)=C3(0) =C4(0)=2 给出初值误差、观测误差较小情况下的计算效果

表 1 试验列表 Table 1 The detailed setup of the experiments

表 2 (Table 2)

表 2 射击次数列表 Table 2 Shooting times over each time span 试验 各时间段的射击次数 试验1 [0,1]:4,[1,2]:5,[2,3]:5,[3,4]:5,[4,5]:5,[5,6]:9,[6,7]:6,[7,8]:4,[8,9]:5,[9,10]:5,[10,11]:5,[11,12]:5 试验2 [0,0.5]:10,[0.5,1]:10,[1,1.5]:10,[1.5,2]:10,[2,2.5]:7,[2.5,3]:10,[3,3.5]:10,[3.5,4]:10,[4,4.5]:10,[4.5,5]:10,[5,5.5]:10,[5.5,6]:9,[6,6.5]:8,[6.5,7]:8,[7,7.5]:9,[7.5,8]:10,[8,8.5]:10,[8.5,9]: 10,[9,9.5]:10,[9.5,10]:10,[10,10.5]:10,[10.5,11]:10,[11,11.5]: 10,[11.5,12]: 10 试验3 [0,1]:5,[1,2]:14,[2,3]:18,[3,4]:20,[4,5]:18,[5,6]:17,[6,7]:20,[7,8]:20,[8,9]:20,[9,10]:20,[10,11]:17,[11,12]:8 注：表 2给出了各时间段的实际射击次数。在试验1和2中,最大射击次数MaxShooting=10,试验1中的射击次数大多约为5次,试验2(加密试验1的观测时间,总观测次数为试验1的两倍)中的射击次数大多为最大射击次数10。我们看到,加密观测时间不仅使时间段数增多,而且每个时间段需要的射击次数也明显增大。因此,加密观测时间付出的时间代价是显然的(由图 2我们将看到,加密观测时间得到的反演结果也不一定更精确)。

表 2 射击次数列表 Table 2 Shooting times over each time span

在本文的数值试验中，经验参量Q_coeff第一个分量对应Q₁，第二个分量对应Q₂。源清单的排放速率均为Q₁(t)≡1,Q₂(t)≡50,t∈[0,12]。#

试验1的目的是检验本算法对污染物浓度的反演效果和污染源反演能力。从试验1的结果我们看到：本算法可使用较少射击次数反演源(表 2和表 3)，并且能处理源清单中的大误差、初值大误差、观测值在个别时间点的大误差(图 1；图 1b’在时刻6的观测值存在大误差)。反演后，没有用观测值反演的物种浓度也有可能得到改善(图 1d)。在源排放速率变化很小的情况下，我们从反演结果仍能够看到源变化趋势(图 1f’)。另外，试验1与附录B展示了负值纠错功能的必要性。

表 3 (Table 3)

表 3 试验1 的源反演过程 Table 3 The shooting process of emission in Experiment 1 时间段 各时间段的源反演射击过程 Q1 Q2 [0,1] 1、0.0001、－10、－90.5470 50、0.0000、－10、－96.3230 [1,2] 1、2.4883、5.6614、9.3787、10.7136 50、0.3490、0.3291、0.3109、0.2943 [2,3] 1、2.1170、4.3012、7.0244、8.7374 50、0.3483、0.3240、0.3022、0.2824 [3,4] 1、2.5105、5.8458、9.7524、11.0087 50、0.3457、0.3176、0.2927、0.2703 [4,5] 1、2.3320、5.1421、8.5613、10.1266 50、0.3470、0.3196、0.2951、0.2732 [5,6] 1、2.0201、3.9626、6.3092、7.7880、8.3781、8.4318、8.4863、 8.4605 50、2.0652、9.5764、21.2562、14.9676、19.5630、16.2727、18.7078、16.9272 [6,7] 1,2.2451、4.6349、7.5352、9.1663、9.5831 50、0.0707、0.0142、0.0028、－10.0000、－16.1352 [7,8] 1、2.9597、7.5530、12.1761 50、0.3641、0.3618、0.3595 [8,9] 1、2.1208、4.2072、6.7743、8.4383 50、0.3568、0.3485、0.3407、0.3333 [9,10] 1、2.4418、5.3846、8.8472、10.2822 50、0.3469、0.3308、0.3159、0.3021 [10,11] 1、2.4833、5.5287、9.0774、10.4518 50、0.3417、0.3220、0.3041、0.2877 [11,12] 1、2.2929、4.7884、7.7914、9.3893 50、0.3436、0.3268、0.3113、0.2969 注：表 3 给出了试验1的源反演过程。在试验1中，第一个物种和第二个物种源清单的排放速率分别为Q1(t) ≡ 1,Q2(t) ≡ 50, t ∈[0,12] 。故在每个时间段，Q1 和Q2 分别从1 和50 起始。初值大误差使得Q1和Q2在[0,1]用负值纠错；观测值大误差使得Q2在[6,7] 用负值纠错。Q1源清单的值小于实际值，从表3 看出，Q1的调整是稳步增长的过程；Q2源清单的值远远大于实际值，从表 3看出，Q2的调整是先迅速降到实际值附近，然后微调的过程。

表 3 试验1 的源反演过程 Table 3 The shooting process of emission in Experiment 1

试验2除了观测时间间隔不同，其他均与试验1相同(表 1)。对比试验1与试验2，其他设置不变，只在时间上加密观测数据不一定得到更准确的反演源(图 2)，而且射击次数可能会明显增加(表 2)。

图 2

Fig. 2

图 2 试验2 的结果：（a）浓度C1；（b）浓度C2；（b’）为（b）C1 的局部放大；（c）Q1源排放速率；（c’）为（c）Q1去掉纠错阶段；（d）Q2源排 放速率；（d’）为（d）Q 2 去掉纠错阶段 Fig. 2 Results of Experiment 2: (a) Concentration C1, (b) concentration of C2, (b’) partial enlargement of (b), (c) emissions rate of Q1, (c’) panel (c) without error correction phase, (d) emissions rate of Q2, (d’) panel (d) without error correction phase

试验3给出了在初值、观测值误差不太大情况下的算法效果。反演后的源较好地与真值吻合(图 3)。

图 3

Fig. 3

图 3 试验3 的结果：（a）浓度C1；(b) 浓度C2；（c）Q1源排放速率；（d）Q2源排放速率 Fig. 3 Results of Experiment 2: (a) Concentration C1, (b) concentration C2, (c) emissions rate of Q1, (d) emissions rate of Q2

污染源反演问题常作为统计估计问题处理。比如，四维变分法、集合卡尔曼滤波是人们常用于反演的统计方法。四维变分法本质上是使目标函数极小化的泛函问题。集合卡尔曼滤波是一种用蒙特卡罗短期集合预报来估计预报误差协方差的四维变分方法。集合卡尔曼滤波无需模式反向积分，无需预报模式的切线性模式和伴随模式，因此与四维变分法相比，更容易实现、可移植性强。然而，四维变分、集合卡尔曼滤波等统计方法存在下面几个弱点：先验分布的高斯假设；线性的分析方程在模式为强非线性时，对结果造成不利影响；对观测误差较敏感；污染源排放的先验误差估计对反演结果有较大影响(朱江等，2006)。

本文提出的基于打靶法思想的大气污染源反演自适应算法本质上是系统控制问题，与统计量无关，且计算过程中用到的主要信息为观测值。只要观测值本身合理、有效，就可以反演出可信的源排放速率。如果我们具备可靠的观测误差信息，则可利用该信息先对观测数据作前处理校正，然后再应用本文设计的打靶算法进行源反演。

该基于打靶法思想的大气污染源反演算法弥补了四维变分、集合卡尔曼滤波等统计方法的不足，具有如下特点：(1)能处理源清单中的大误差、初值大误差、观测值在个别时间点的大误差；(2)无需先验分布假设及误差估计。在实际中，对污染源进行反演是一个很复杂的问题。污染源与观测站点的空间分布不均匀性、可观测的物种少于需要反演的污染源物种等局限性会使反演结果不唯一。在今后的工作中，我们会考虑通过有效地增加约束条件或信息来处理这一问题。

气象场误差、模式分辨率误差、化学机制误差等因素均可造成模拟与观测不一致。如果忽略这些模式误差、仅通过调整排放源来减小模拟和观测之间的不一致，那么源反演结果将会出现偏差。集合卡尔曼滤波/平滑考虑了各种误差对源反演的 影响(Tang et al.，2013)。如何将集合卡尔曼滤波/平滑与本文提出的基于打靶法思想的自适应算 法结合起来解决源反演问题也是我们今后将要探讨的。

致谢 感谢中国科学院大气物理研究所朱江研究员及审稿专家们对本文提出的宝贵意见。