引言

污染物在大气中随风飘流，又因为湍流而扩散，这一过程，我们通常用平流-扩散方程来描述。在这里为了研究的方便，我们假定:方程是精确的；在给定了气象场、扩散系数、污染源及初、边值等条件后，方程的解是适定的。那么，我们可以把方程 (或者方程对应的数值模式) 看作是一个确定的系统，该系统在输入了上述这些条件 (原因) 后，就会得到一个确定的浓度场时空分布 (结果)。由原因求结果，称为正问题 (direct problem)。由结果求原因，称为反问题 (inverse problem)。求解反问题，称为反演^[1]。

大气边界层内，紊乱的气流造成气团的垂直混和，主宰了动量、热量、水汽和污染物等在垂直方向的输送，并进而影响着天气和气候系统。由于湍流的紊乱无序，我们通常不是研究其瞬间的性态，而是试图对它的总体效应 (扩散系数) 进行描述 (参数化)^[2，3]。关于大气边界层的参数化方案很多，但是要实际应用，必须先确定若干经验参数。这些参数通常是间接的、难以测量的，没有合理的参数数值，再好的方案也难以付诸实用。要从相对容易测量的物理量 (比如浓度) 反求扩散参数，就是一个典型的反演问题。

反演问题和正问题，初看起来都是已知一部分信息，通过信息间的联系，求未知信息，似乎正问题能解决了，反问题也不难。但是正问题是由原因求结果，顺着事物发展的趋势进行推演，是在模拟物理上可实现的过程。反问题是由结果回溯原因，这一过程，一般是物理上不可实现的 (不可逆)。正问题的解要求是存在、唯一且稳定，反问题的解则很难满足这些适定性要求。另外，即使正问题是线性的，其相应的反问题也完全可能是非线性的。因此反问题往往比正问题要困难得多，迄今为止反问题的研究还不是很成熟的。

本文用最优控制的观点来定义反问题的解，以“K理论”为例，用变分伴随方法求解最优控制问题，从而得到满足要求的最合理的解 (扩散系数)，并与扩散系数的“真实值”进行比较。

1 反问题的定义及求解方法

首先我们对反问题的解进行定义。

用最优化方法，得到一个扩散系数分布，使得预报值与观测值之差的平方和 (即目标函数D) 为最小。

(1)

其中Φ=Φ(x，y，z，t) 是浓度的预报值，Φ_O是观测值，p是权重因子。

且

(2)

边界条件:

其中f=f (x，y，z，t) 是污染源，K_H和K是水平和垂直扩散系数。

这里是通过一个最优控制问题来定义反问题的解。目标函数为D，控制变量 (即反问题的解) 为垂直扩散系数K，约束条件为方程 (2) 及其定解条件。不难看出，目标函数D是通过方程 (2) 而成为K的函数 (泛函)。我们要做的就是调整控制变量K，使得目标函数值最小，此时的K就是最合理的，因为由它预报的浓度最接近于实测值。

下面，我们用变分方法来求解上述优化问题。

假设对垂直扩散系数K作一个微小的调整K →K +δK，相应的，Φ→Φ+δΦ，目标函数D →D +δD，

则有：

略去高阶微量，可以导出

对方程 (2) 进行扰动，得到其切线方程:

这是当扩散系数发生微小的变动时，浓度的变化量所应满足的方程。

其中切线算子

相应的伴随算子

关于伴随算子的表达式，在参考文献[4]中第49~52页有详细推导。

定义切线方程的伴随方程及其定解条件^[4] :

下边界Σ₀处,

上边界Σ_H处,

侧边界Σ处, 当un < 0时, Φ^*=0

则由伴随算子的Lag range等式，有

如果把控制变量K看作是时间和空间里的广义向量，那么正是目标函数D关于控制变量 (向量) K的梯度，这个梯度是通过求解方程 (2) 和伴随方程得到的。

利用梯度信息，我们就可以应用各种最优化方法 (如最速下降法，共轭梯度法，拟牛顿方法等) 来使目标函数下降，从而使反问题得到解决^[5]。

优化步骤:

(1) 猜测一个扩散系数分布K (x，y，z);

(2) 把K代入正模式，计算Φ，

(3) 计算目标函数D，计算p (Φ-Φ_O)，代入伴随模式，计算Φ^*，

(4) 计算目标函数的梯度

(5) 用最优化算法更新K;

(6) 判断‖∇_KD ‖≤ε?如果是，则停机；否则返回 (2)。

2 反演结果及分析

为减少计算量，用二维模式 (水平坐标x，垂直坐标z) 进行了试验。水平方向风速均匀，垂直方向无平均风速。先根据方程 (2) 写出正模式^[6]，再用共轭码方法由正模式写出切线模式及其伴随模式^[7]。作为数值试验，为了验证公式推导及反演方法的正确性，先假定一个“真实”的扩散系数，该垂直扩散系数只是沿着垂直坐标变化，在水平方向是均匀的，即K=K (z)，如图 1所示。用这个“真实”的扩散系数模拟一个浓度分布作为“观测”值，图 2所示为测试源下风方向x=8处的浓度垂直分布，容易看到扩散系数数值较大的地方对应的浓度值比较低，浓度与扩散系数的相关是我们利用观测值来反演扩散系数的依据所在。反演时，先预估一个扩散系数分布，然后对估计的扩散系数进行优化，使得预报值充分逼近观测值，并将优化后的扩散系数值与“真实”值作比较，以检验方法的有效性。数值试验中对扩散项用了两种差分格式进行研究。

由于垂直扩散系数不是常数，为了能得到它在垂直方向的分布，数值试验中，沿一条垂直线等距离地分布多个点源，使得点源下风方向各点的浓度都不太小，而彼此之间又有差异。实践表明，这样的浓度分布，有利于反演出各点的扩散系数。

在数值计算中，为了方便起见，对各个物理量都作了无量纲化处理，这并不改变问题的实质和反演方法的有效性。为了精细地反映垂直扩散系数的效果，计算网格在z方向是75个格点，而在x方向是25个格点。测试源分布在x=5的垂线上。在下游一定距离后，物质浓度在垂直方向趋于均匀化而对扩散系数不敏感，在数值误差的限制下，这样的浓度分布不利于精确反演扩散系数。因此，数值试验中采用的是x=5~10之间的浓度分布作为“观测值”来进行反演。

2.1 格式一的结果

格式一:

在一条垂直线上，每5个格点设置一个点源。反演得到的扩散系数垂直分布如图 3所示。它是一个锯齿状分布，看起来很不理想，但是当把“真实”的扩散系数叠加到图线上后发现，在其它的点上吻合得不错，“真实”曲线正是反演曲线的包络线。

为了探究“锯齿”现象的原因，画出目标函数梯度值的垂直分布，如图 4所示。可以发现，图线上周期地出现数值很小几乎为0的点，对应的正是那些反演不好的点。由于算法中调整扩散系数都是在原估计值上加上一个梯度的倍数，若该点的梯度值一直为0，则该点一直得不到调整，当然也就反演不出来了。之所以在这些点上梯度很小，可以大致地解释如下:为了能反演出扩散系数的分布曲线，测试源间隔不能太大，格距也不能太大，使得浓度的变化率远大于扩散系数的变化率，这样才能捕捉到扩散系数的变化趋势。相对于浓度的变化来说，扩散系数的变化是缓慢的，在一段小距离内，可近似地看作是常数。又因为点源之间的距离是奇数倍格距，因此，两个点源之间的浓度极小值正好位于两个相邻格点之间，这两个相邻格点上的浓度值也就几乎相等，差分后就几乎为零。这就解释了该点处梯度值一直很小的原因。

改变点源的分布，每四个格点设置一个点源，可以避免出现梯度值始终很小的点。果然，反演的扩散系数曲线与“真实”曲线吻合得很好，如图 5所示。

进一步的试验表明，当点源的间隔为奇数，反演后图线就有“锯齿”，偶数时则不出现。点源间隔太大，则一些点处浓度太小，反演效果不好；反之，点源太密，浓度太均匀，扩散效应弱，反演效果也不好。

由于这一格式只用到了相邻两个点的扩散系数，只具有一阶精度，下面尝试一种精度更高的格式。

2.2 格式二的结果

格式二:

该格式是二阶精度。

先取点源间隔为5(奇数)，反演的结果出现意想不到的振荡现象，如图 6所示。进一步分析发现，振荡解的平均值与”真实”值曲线相吻合，把这一振荡的扩散系数曲线代入到模式中，计算的浓度与“观测”值是一致的，如图 7所示。这说明了用这一格式反演的结果不是唯一的。实际上，该格式是用相邻格点扩散系数的平均值来计算的，因此，只要在“真实”值上叠加一个等幅波 (波长为二倍格距，在格点处取极大极小值，振幅可以任取)，就可以得到一致的浓度分布。

但是取点源间隔为4后，却不出现上述的振荡现象，最后收敛到“真实”值。而且，离源较远，靠近边界的点也反演得不错，如图 8所示。

多次试验结果表明，虽然反演问题的可能解不是唯一的，但是，最终得到哪一个解，却具有必然性。当点源间隔为偶数，最终收敛到真实值。当点源间隔为奇数，肯定得到振荡解，而且振荡的幅度，正比于初估值与“真实”值的差值。

3 结论

综上所述，我们用最优化原理定义了垂直扩散系数反演问题，并在二维情形下，就两种差分格式进行了多次数值试验，分析了结果并探讨了问题出现的原因。我们可以看到，反问题和正问题是有着本质区别的，无论是在问题的表述，解的定义，影响因素和求解方法上都有着很大的不同。但是，反问题的解决，又是以正问题的完整解决为基础的，正问题模式的不断完善是顺利求解反演问题的重要前提。

我们从大气边界层垂直扩散系数的反演这一具体问题的数值试验中，可以看出，差分格式对于反演的结果起着非常重要的作用。一些精度高、在正问题中用得很好的格式，到了反问题中就可能不适用，或者适用条件受到限制。另外，反演效果的好坏，还与试验测量的设置有着直接的关系，数值试验有助于找到合理的试验设置方案，为真实测量提供理论依据。