引言

在20世纪中后期, 全球社会经济和工业化过程的飞速发展, 人口急剧的增长给大气环境带来了前所未有的巨大压力, 人们越来越认识到保护人类赖以生存的大气环境的重要性与紧迫性.从50年代后期起, 一些国家为了控制或减少危险大气污染事件的发生, 开始进行大气污染气象条件预报 (污染潜势预报), 例如WMO的121期技术报告中就介绍了70年代前后的许多污染潜势的预报方法, 其中有使用统计表对各种气象要素及其组合加权记分的预报法, 也有直接用箱模式导出的标准化浓度i/Q (i为混合层内污染物平均浓度, Q为面源源强) 作为污染潜势或空气停滞指标, 因为按标准化浓度计算公式i/Q=L/(2HV), 它仅与城市尺度L, 混合层厚H及混合层内的平均风速V有关.我国在1980年秋冬到1981年春也进行过一次北京的城市空气污染气象预报试验^[1], 这是一次以天气学方法为主, 兼采用经验统计表法的主观定性预报.后来的3个五年计划里, 我国污染气象条件方面的研究工作大多是在制定大气污染物排放标准、大气环境影响评价、酸雨成因研究、气候变化等重大项目中配合进行的.在这个阶段, 积累了许多国内研制和国外引进的大量各种尺度、各种类型的有关大气污染物传输、扩散和化学演变模式.例如用于区域酸沉降的RAMS系统, 用于区域和全球烟云扩散、轨迹计算的先进的Hysplit模式等等都已引入我国.同时也为这些模式积累了许多污染气象或气候参数, 加之近年来计算技术的高度发展, 为当前开展大气污染浓度预报和潜势预报都奠定了一定的基础.

大气污染物的浓度预报, 可分统计预报和数值预报两类, 统计预报是将历史上的污染物的浓度监测值与前期的气象条件 (气象要素及其时空分布、天气类型) 联系起来, 建立具有一定信度的统计关系并利用该关系对未来的污染物浓度进行预报.一般而言, 在有长期的监测数据的城市, 应用统计方法进行预报是不困难的.如果建立统计关系前后的城市大气污染物的排放源的源强及其时空分布没有太大的变化, 统计预报的效果也是相当好的.大气污染物浓度数值预报, 一般是在气象场预报模式的基础上数值求解污染物的有源汇传输扩散方程, 需要较详尽的源强及其时空分布资料和时空分辨率很高的气象预报模式, 在目前条件下, 用上述思路进行的浓度预报只能在有条件的个别城市进行.

在全国范围内则可进行污染气象条件或大气污染潜势指数预报, 它将在大尺度范围内给城市大气污染浓度预报提供有意义的气象背景.为了预报大气污染潜势, 本文对平流扩散方程积分后进行了物理分析, 建立了对实时监测浓度有记忆能力且能感受非临近地区浓度的箱格积分模型, 找出了能代表大气对污染物清除或稀释能力的量, 即大气污染潜势指数PPI, 并将它与可预报的气象要素联系起来.此外本文还在潜势预报的基础上发展了由前期污染物浓度及气象场监测值进行浓度预报的无源参数的箱格数值预报模式, 并给出了该模式在一些台站的业务预报试验的结果.

1 平流扩散方程的积分与箱格预报模型的建立 1.1 平流扩散方程的积分

如果大气污染物浓度c是由体积f内若干位于源强分别为q_i源产生, 风速为, 湍流交换系数为二阶张量, 干沉降速度为, 湿沉降速度为, 降水率为R, 降水清洗比为w_r, 那么大气中气载污染物不考虑化学反应的平流扩散方程, 按污染物守恒的原则, 可写为:

(1)

根据定义, 和只有铅直方向分量且指向地面. 的值等于降水率R和降水清洗比w_r的乘积

(2)

对式 (1) 在体积f内积分后再对体积平均可得:

(3)

其中

(4)

因为有, 在低速大气中, 空气十分近似于不可压流体, 故, 那么式 (3) 左方第二项可写为

(5)

式 (5) 的最右方为曲面积分, 为包围体积f的表面, 其法线方向指向体积外为正.式 (3) 右方第一项为体积f内单位时间大气污染物排放 (包括点源、面源与线源排放) 的总量Q, 即:

(6)

式 (3) 右方第二项为体积f内大气污染物的干、湿清除项.

(7)

式 (3) 右方第三项湍流通量项, 在城市尺度、中等风速的条件下该项的大小约为平流项的百分之几, 可以忽略, 但在小风速或静风条件下该项必须保留.因为按原始定义, 湍流通量可表达为:

(8)

若将湍流通量用虚拟湍流输送速度表示, 即

(9)

那么式 (8) 可改写为:

(10)

这里要指出的是, 在式 (1) 的原来意义上, 的方向总是从浓度高值区指向低值区.在本文讨论的情况下, 垂直于体积f的表面, 其法线方向指向浓度低值区为正

(11)

于是式 (3) 可写为:

(12)

上式是平均浓度的预报方程, 它的最右方的积分项表示体积τ内大气对污染物的清除能力, 这些清除是由通风扩散稀释和干湿沉降过程进行的, 若设该积分值对平均浓度的比值为:

(13)

V_c代表着该时段的大气通风扩散稀释和干湿沉降的总能力的平均值.这时式 (12) 立即可解得:

(14)

其中积分常数表示初始平均浓度.在给定的时段W_T内, 若V_c和Q都与时间关联甚小, 上式可解出

(15)

由式 (15) 绘制的单位源强浓度曲线图 1 (这里已取/Q=1) 可见, 在大气对污染物的清除能力很大时, 如V_c=5.0, 那么该体积内的污染物将很快达到平衡浓度1/V_c=0.2;相反, 若大气对污染物的清除能力很小, 如V_c=0.25, 那么该体积内的污染物将会在初始浓度/Q的基础上缓慢地向平衡浓度1/V_c=4增加.这种特点很相似于实际大气中的过程:连续数天的静风条件, 将使大气污染状况日趋严重; 而一场持续数小时的强风, 就会使空气立即清新起来.当平均时段较长, 这时V_c将为时间的函数, 其平均浓度应按式 (14) 的积分计算.事实上, 上述积分是在一个单体积f中进行的, 只要V_c是常值, 时间间隔W_T足够长, 或者式 (12) 中的时间局地变化项为零, 它和经典的箱模型全同.在一般情况下多了非常重要的时间变化项, 从而能用来预报在小风条件下逐渐恶化的污染过程.此外平流扩散方程的差分模式只能描述短时间步长的过程且存在污染物不守衡问题, 而使用式 (15) 就可避免这些问题.在任何情形下表示积分的V_c总可由箱格模式求得.由式 (13) 见到, V_c不但依赖于箱体内的参数而且也依赖于箱体外参数.要准确计算某个箱体的V_c值就必须将所有有关箱体的浓度计算出来, 只有这样才能准确地计算出该箱体的浓度变化, 于是根据具体资料与计算条件设计不同计算模式就成了一个不可缺少的步骤.

1.2 箱格预报模式

在箱格预报模式中, 设中心位置在x=(i+1/2) WX, y=(j+1/2) WY, z=(k+1/2) WZ, i=1, 2, …, M; j=1, 2, …, N; K=1, 2, …, L (在黑圆点坐标网中x=iWX, y=jWY, z=kWZ) 的箱体abcd-efgh的体积f=WX×WY×WZ, 它沿X-方向水平体长为WX, Y-方向水平宽度为WY, 高为WZ, 见图 2a, 设箱格风速矢量坐标网为图 2b, 浓度标量坐标网为黑圆点网, 见图 2c.在t=nWT时体积内浓度平均值为c_i,j,k, 通过界面abcd、efgh、aecg、bfdh、cdgh、abef并与之垂直的平均风速分别为:

(16)

其中, u_i,j,k, v_i,j,k, w_i,j,k为c点上风速矢量的3个分量; u_i,j+1,k, v_i,j+1,k, w_i,j+1,k为d点上风速矢量的分量; u_i,j,k+1, v_i,j,k+1, w_i,j,k+1为a点风速矢量的分量; u_i,j+1,k+1,v_i,j+1,k+1, w_i,j+1,k+1为b点上风速矢量的分量; 其它类推.这些垂直通过界面的风速所携带的浓度分别为c_i、c_i+1、c_j、c_j+1、c_k、c_k+1, 可用黑圆点坐标网络上的浓度表示为:

(17)

在体积界面abcd、efgh、aecg、bfdh、cdgh、abef上虚拟湍流输送速度分别为u_t,i, u_t,i+1, v_t,j, v_t,j+1, w_t,k, w_t,k+1输送的浓度分别为:

(18)

此外在k≤0的平面上, w_k=w_t,k=0, 干、湿沉降速度v_d,i,j,k、v_w,i,j,k都只有k方向分量, 其方向向下 (其中对非重力沉降物v_d只有在与地面相贴的箱表面上有有限值, 在其他箱面上均为零值), 并认为干沉降在时间间隔WT内, 在形式上仅将污染物降到临近的下一个箱体.那么按式 (13), V_c的平均值就可以写为:

(19)

式 (19) 仅考虑两个相邻体积间的浓度的扩散过程, 如果平均的时距较长, 湍涡的尺度谱分布很广且与单个箱尺度相比, 尺度很大的湍涡的能量占总湍能的比例不能忽略, 这时应该按文献[2][3]的概念将u_t, v_t, w_t按拉氏自相关尺度分解为:

(20)

显然有:

(21)

考虑到大湍涡的穿越 (transilient) 扩散效应^[4], 式 (19) 可写为:

(22)

其中

(23)

由式 (15) 还可定义出每个箱体内污染物变化的特征 (惯性) 时间尺度T_i,j,k, 它的倒数应满足下式:

(24)

按式 (19) 很容易求得:

(25)

(26)

式 (25) 说明, 箱体内惯性时间尺度T_i,j,k的倒数为箱体三个方向的惯性时间分尺度的倒数和.在其他条件相同的情况下, 同体积长条形的箱体比正方形的箱体的惯性时间短.这也说明, 圆形或正方形的城市与长条形城市相比, 大气污染物在市内滞留的时间要长.以长边垂直于盛行风且上风方污染浓度低的城市的惯性时间为最短.

箱格预报模式的计算流程见图 3, 由气象预报模式给出t=nWT时刻的风场、降水量分布 (计算湿沉降速度)、各种尺度的虚拟扩散速度场等初始场, 还要给出各箱体内的源强及浓度值 (要注意的是, 本模式在短时间内对初始浓度依赖性较大, 在长时间后则对源强依赖性大), 通过式 (19) 或式 (22) 求得各箱体的通风扩散稀释和干、湿沉降的总能力的平均值V_c (同时也给出污染物的干、湿沉降量) 后, 再按式 (15) 计算出在一个时间步长δT即t=(n+1) δT时刻的各箱体的浓度值, 该值也就是下一个时间步长开始时的浓度初值.这里的时间步长可以与气象模式中的步长一致, 也可远大于气象模式中的时间步长, 因为它们之间并无任何直接联系.

一般地说, δT选取得较小时模式对前期浓度记忆强: δT选取得大时, 模式对源强依赖大, 在使用中可选δT使其满足.在求取单个时段内的V_c时, 本箱格模式的作用与传统箱格模式的差异在于各箱体对局地浓度有记忆并且对远方箱体的浓度有感受能力, 能同化实时监测浓度.箱体可堆满对流层 (这时浓度c应以浓度与大气密度的积代替), 也可仅堆置于大气边界层 (这时箱体堆积总高度为边界层高度, 因此箱体侧边界面积的高度依赖于当地边界层高H, 例如取δZ=H/L, L为垂直方向箱体个数), 甚至可以将箱体离散地放在任意指定的地点.

2 空气污染潜势指数与污染指数的确定 2.1 空气污染潜势指数的确定

在大气边界层顶和地面之间建立一个箱体, 设u=u_i=u_i+1 > 0, u_t,i=u_t,i+1, v_t,j=v_t,j+1, v_i=v_i+1=w_i=w_i+1=0, WZ=H为混合层深度, 上风箱体的平均浓度c_i+1不为零, 其它相邻箱体内的平均浓度都保持为零, 再设:箱底面积S=WX×WY, 那么按式 (19) 可得:

(27)

在小风情形下, 湍流输送虚拟速度u_t, v_t, w_t, 以及干、湿沉降将成为V_c的主要项. H为混合层深度时, 垂直湍流输送虚拟速度w_t=0, 再设上风箱体的平均浓度c_i-1 < < c_i, 平流风速远大于湍流输送虚拟速度, 式 (27) 即为

(28)

按式 (15) 可得平衡浓度:

(29)

WX按文献[5]取为, 那么, 上式可写为:

(30)

式 (30) 即为中国国家标准GB/T3840-91制定SO₂等气态大气污染物总允许排放量的基本公式原型, 在那里已令干、湿沉降速度为零, 且按地区给出了控制参数A_i=u_i×H×c/2, 其中u_i×H为各地多年平均通风量.由此, 对考虑干、湿沉降的大气污染物的允许排放总量由下式计算:

(31)

式中、为多年平均的干、湿沉降速度.对不考虑干、湿沉降的大气污染物总允许排放量的基本公式为:

(32)

由式 (31) 或 (32) 按面积S、标准浓度及所在地区的参数A_i、、取得Q_{总允许排放量}后, 代入式 (15) 得到

(33)

该式给出箱内初始浓度为且存在标准源强时的浓度在一个时间步长上的浓度变化.如果直接将式 (31) 或式 (32) 代入式 (33) 便得到标准源强在实际气象条件下生成的浓度与在标准气象条件下产生的浓度的比值, 这个比值本文称为污染潜势指数PPI.该指数事实上表示的是实际气象条件下的通风扩散稀释和干、湿沉降清除大气污染物的总能力与标准气象条件下这种能力的比值.考虑沉降时有:

(34)

其中就是标准气象条件下的大气稀释和沉降清除污染物能力V_c标准或者写为:

(35)

不考虑沉降时, 只要将上式中的干、湿沉降速度置零即可.

从式 (35) 可见, PPI的预报和源强与浓度初值场的关联不大, 只是在模式第一次运行时需提供, 这时可假定=1, 经过数小时或数十小时运行后, 其结果将与初值无关.这两个式子的其它自变量都是由气象模式提供如V_c (见式 (28))、H, 或由地理参数决定如面积S、A_i.由于模式的运行完全脱离污染源和污染浓度监测场, 它预测的只是大气通风扩散稀释和干、湿沉降清除大气污染物总能力的大小, 并不能给出大气的实际污染状态.如果要预测有限空间f的平均浓度, 那么在各个f体积的污染源源强显然是不可缺少的.为了弥补源强资料中通常存在的不确定性, 污染浓度监测场的使用将会有效地提高预报效果.

2.2 空气污染指数的确定

污染指数PSI, 对单种污染物而言, 它定义为真实源强在真实气象条件下产生的平均浓度c对标准浓度c_标准的比值即: PSI=c/c_标准.例如要计算某个箱体的污染指数PSI, 则仅需将式 (33) 中的Q_{总允许排放量}更换为Q_{真实排放总量}并除以c_标准, 按定义得:

(36)

在仅有监测浓度的季平均值而无真实排放总量时, 按式 (29) 可得预报PSI的近似式 (37):

(37)

如果存在着较好的大气污染物监测网, 能按一定的时间间隔提供各种原生的大气污染物的PSI值, 在不考虑化学过程情况下, 其污染指数可近似用下式计算.

(38)

式 (38) 中的PSI_t=nWT可随时用相应时刻的实测值取代.无实测值时则用前一步长的计算值迭代计算.

3 空气污染潜势指数与污染指数的试预报

采用京津冀中尺度试验基地的MOMS业务模式作为气象预报模式^[6], 格点数为61 (纬向)×46 (经向)×10 (高度), 水平格距60 km, 顶层气压为150 hPa.在使用总体边界层时, 利用e=1和e=0.95两层的温度差及文献[7]给出的判据, 判别大气稳定度, 再由稳定度类别按文献[8]的公式与参数确定大气边界层高度, 再以边界层高度以内的平均风速、当地降水量计算V_c值, 其中干沉降速度是作为参数输入的.当使用高分辨率大气边界层时则直接利用模式输出的边界层高度, 进行上述计算, 不过对边界层高度作了不得低于200 m的限制.

预报污染潜势指数PPI时是使用式 (34) 进行的.计算V_c值时是在气象模式的e=1平面的每一个格点进行的, 在这里放置的箱体的底面积相当于全国50个大中城市的平均面积81 km², 并按国家标准GB/T 3840-91给出A_i, v_d只有在与地面相贴的箱表面上是按文献[9]给出, 其他箱面上均取零值.平均降水量按气候图内插值给出.虚拟湍流输送速度各分量的绝对值的大小则按稳定度的级别 (由箱格的总体理查逊数决定), 取为平均风速的一定百分数, 静风 ( < 0.5 m/s) 时, 则取为各稳定度条件下的e_u、e_v、e_w的典型值.其它地区未置箱体, 其实质相当于在其它地区设置了浓度永为零的箱体.WT取10分钟.

1997年3月1日以来, 在京津冀中尺度试验基地, 使用实时气象资料对全国大气污染潜势指数进行了试预报.图 4就是1998年10月18、20、21和22日我国华北的污染潜势指数 (图右色标数即为PPI指数) 分布与500 hPa等压面形势图, 黑色为等高线 (gpm), 红色为等温线 (K).图中PPI=1, 表示当地大气的通风扩散稀释和干、湿沉降清除大气污染物的总能力与标准 (多年平均) 气象条件的相同. PPI=0.5表示当地大气清除污染物的总能力比标准条件的强1倍. PPI=2.0, 表示当地大气清除污染物的总能力弱1倍, 此地即使按标准源强排放, 由于当天的气象条件不利于大气污染物的稀释扩散与清除, 惰性污染物的浓度将为标准浓度的两倍.以北京为例, 18日 (图 4a) 由于有新鲜的冷空气南下, PPI在0.25与0.5之间, 大气的污染潜势很低.到20日, 冷平流已经消失, 北京的PPI在1.5与2.0之间, 空气污染潜势已较高. 22日北京的污染潜势指数接近2.5, 显著高于邻近的河北与天津.实际上这一周北京监测站发布的污染指数是该月最高的周均值. 1998年10月模拟的结果与北京监测站发布的污染指数的对应关系很好, 已另外撰文叙述.

为了对上述计算的效果作出估计, 于1998年8月至12月, 在上海市使用本文模式和式 (38), 对SO₂、TSP和NO_x的日平均PSI值, 逐日进行了预报, 同时有预报值与监测值的日数为118天, 并将PSI预报值与相应的监测值做了相关对比.在这里PSI指数定义为:

(39)

SO₂、TSP和NO_x的日平均PSI预报值与监测值的相关图见图 5~7.其中SO₂的PSI日平均预报值与监测值的相关系数为63.5%, TSP为61.8%, NO_x为45.0% (样本数都为118).

NO_x的预报水平明显低于SO₂及TSP的预报水平.这是因为NO_x主要是城市车辆等流动源排放的, 其源强的时间分布很不均匀, 在预报中未考虑排放的时间系数, 估计在考虑时间系数以后, NO_x的预报效果会有一定改善.

鉴于污染潜势预报的基本原理与PSI指数的预报原理相似, 且不依赖于局地源强的大小, 其预报准确性不会低于PSI的预报准确率.

4 结论

(1) 由平流扩散方程积分后再离散化而建立的箱格预报模型能克服经典箱模型对时间步长δT远大于f/V_c的要求和不能预测小风条件下逐渐恶化的污染状态的缺点; 同平流扩散方程的差分模式相比较, 箱格预报模型对时间步长和空间步长的要求特别灵活且不存在不守衡问题; 此外箱格预报模型对初始浓度有很好的记忆力, 并对非相邻地区的浓度有感受力, 因而特别适用于各种时间、空间尺度或多尺度的空气污染预报.

(2) 本文空气污染潜势预报所使用的PPI指数与源强和浓度初值场的关联不大, 但能很好地反映实际气象条件下的通风扩散稀释和干、湿沉降清除大气污染物的总能力, 与经典的空气污染潜势预报方法相比, 结构严谨、物理意义更清晰明确.

(3) 本文给出的空气污染指数PSI的预报公式和方法对源强和浓度监测的要求灵活, 既可以使用源强进行预报也可以使用浓度监测数据进行预报, 或同时使用二者进行预报.从初步试预报的检验结果来看, 其效果还是可接受的.

(4) 要进一步提高预报准确率, 首先要提高气象场, 特别是风场的预报效果, 还应该在模式中放入一些重要的化学过程.当然, 在进行O₃预报时必须考虑化学反应过程和太阳辐射强度的预报模式的引入.

致谢 本文的上海市污染指数实测与预报值由邵德民高级工程师提供, 在此表示感谢.