引言

人工增雨效果检验是人工影响天气研究的重要课题，人工增雨效果统计检验是人工增雨效果评估的主要方法之一。对于“人工增雨效果”如何评估这一问题，朴素的想法是选择作业影响区以外的区域或者天气相似的非作业日的降水情况与影响区的降水情况进行比较。从科学的角度实现这种方法，一般是根据一定假设条件的历史样本和作业观测样本，建立对比区与影响区的历史回归方程，进而结合统计学方法给出评估的结果。早在1982年叶家东等(1982)在其著作中就系统地阐述了人工影响天气效果检验的各种统计数学方法。胡志晋(1979)、程克明等(1980)介绍了协变量(气象参量)在人工增雨效果统计评估分析中应用的方法。近年来，随着新一代多普勒天气雷达、机载云物理探测仪器等设备的发展，还开展了利用新一代雷达、机载仪器等多种探测资料对人工增雨效果物理检验的研究(Woodley et al, 2004; 唐仁茂等，2009; 2010；张瑞波等, 2010；于丽娟等，2009；刘晴等，2013；陈小敏等，2014；李红斌等，2014)。

在建立人工增雨效果统计检验历史回归方程时，显然涉及以下三个基本问题：影响区的确定、对比区(或对比站点)的选择、回归方程及显著性检验。对于影响范围较小的高炮、火箭作业，可以指定其风向下游以外的区域作为对比区，建立站点对站点或者区域对区域的回归方程(叶家东等，1982；李书严等，2006)；对于大范围的飞机作业，也可以根据FCM方法(浮动对比区历史回归统计检验方法)(Abbas et al, 1999；Koloskov et al, 1999)或者CA-FCM方法(基于聚类的浮动对比区历史回归统计检验方法)(房彬等，2006；翟雨等，2008)基于一定的统计学依据浮动地选择对比站点，建立站点数据对区域雨量的回归方程。

值得注意的是，在人工增雨特别是飞机增雨的效果统计检验工作中，人们往往主要关注“对比区选择”与“回归方程及显著性检验”的各种方法，而认为“影响区的确定”是理所当然的。但事实上，该项工作仍然存在一些亟待解决的技术问题：第一，在统计检验工作中只注重统计学方法，而没有很好地考虑播云催化剂的扩散，没有很好地结合已有的一些云中催化剂扩散的物理学与数值模拟的理论与结果；第二，在日常效果统计检验方法研究和业务应用中，多数仍然采用主观圈定或者仅根据飞机路径大致划定影响区的方式来确定影响区。这些问题的原因在于，大范围的飞机作业所构成的影响区有着很强的不规则性，而这些不规则区域本身以及区域里的降雨量均需要准确地确定。站点雨量观测值实际上是站点一定面积上的降水量，要获得一个不规则区域的降水量首先需要有该区域准确的面积资料，为此，以往在实际工作当中人们通常用相对规则的形状(譬如规则的多边形、圆形、椭圆形等)(吴兑, 1980；叶家东等, 1981)来表示影响区。但是，如果得到的影响区范围不够客观，势必会使统计检验结果的科学性与客观性大打折扣。

我们注意到，虽然现在的GIS二次开发技术已经具备了较强大的功能，但是，目前许多人工增雨统计检验产品中只应用了其显示功能(黄彦彬等, 2011)，而没有使其与统计检验方法相结合并参与到相关的运算中。如果能充分发挥GIS二次开发技术的功能，那么，前述的不规则影响区准确确定问题应该会得到较好的解决。为此，本文介绍一种基于Visual Basic与MapObjects (简称VB+MO)技术, 结合二维线源扩散方程，并且可在人工增雨效果统计评估业务化系统中实现自动地确定不规则影响区范围、准确地计算其面积和体积降水量的技术方法。

1 原理与方法 1.1 飞机播云催化剂的扩散方程

为了便于讨论并且不失一般性，这里假定飞机飞行作业线源是瞬时形成的，催化源强定常，不考虑作业区域的云层条件(如冰晶浓度、过冷液态水含量等)对线源扩散的影响。根据国内外对播云催化剂扩散的研究(申亿铭, 1982；申亿铭等，1986；1996；Gagin et al, 1985)，一种常用的模型是将飞机播云视为线源扩散，并可由三维点源扩散方程得到如下二维的线源扩散方程：

$ q = \frac{Q}{{4{\rm{\pi }}Kt}}{{\rm{e}}^{\frac{{ - {w^2}t}}{{4K}} - \frac{{w\left({z - H} \right)}}{{2K}} - \frac{{{{\left({x - ut} \right)}^2}}}{{4Kt}}}}\left[ {{{\rm{e}}^{ - \frac{{{{\left({z - H} \right)}^2}}}{{4Kt}}}} - {{\rm{e}}^{ - \frac{{{{\left({z + H} \right)}^2}}}{{4Kt}}}}} \right] $

(1)

式中，q为(x, z)处的浓度，H为播撒高度，Q为源强，K为湍流扩散系数。

本文研究的播云影响区范围实际上是催化剂有效扩散的水平范围，这既包括由风速u导致的平流输送作用，也包括由湍流扩散系数K导致的湍流扩散作用。由式(1) 不难发现平流输送作用对于水平扩散距离具有相对独立的贡献。假定K定常，设在水平环境风u_a作用下扩散的浓度为q_a(x_a, z, t)，而无水平风时扩散的浓度为q_b(x_b, z, t)，在t时刻，取它们的顺风向水平扩散前沿处具有相同浓度(即q_a＝q_b)的点分别为x_a和x_b，然后，根据式(1) 很容易求得x_a与x_b的关系式：

$ {x_a} = {u_a}t + {x_b} $

(2)

由式(2) 可见，x_b是在无风情况下只靠湍流作用扩散出的水平距离，而有风时的扩散距离x_a等于平流输送距离u_at与无风速时的湍流扩散距离x_b之和，即在该二维线源扩散模型下平流输送与湍流扩散作用是两个相对独立的过程，可以分开讨论。

对于湍流扩散系数K，除了利用动力学和热力学参量代入式(1) 计算外，通常在实际应用中人们也倾向于通过总结归纳不同地区不同天气条件下K的取值，并依据经验使用(申亿铭等，1996；王俊等，1998)。根据模拟研究(申亿铭等，1986；1993)，以扩散方程[式(1)]为基础选取合理范围的不同K、Q、u、w等参数值进行数值实验，结果表明，有效扩散时间[催化剂扩散稀释到临界值(取为10l^-1)或是催化粒子下落到0℃层的时间]的量级在10³~10⁴ s，在这种时间尺度下，由式(1) 计算得到不计风速的湍流有效水平扩散范围的量级在10⁰ km左右，而高空风速的量级一般最大为10¹ m·s^-1，平流有效扩散范围的量级在10¹~10² km。因此，在水平方向上可近似认为催化剂扩散范围等于其随高空风的平流范围。对于直线飞行的播云作业来说，只要给出有效扩散时间，结合作业层高空风速风向，就可以计算出催化剂的平流扩散范围。鉴于湍流扩散作用导致的水平扩散距离比平流距离小1到2个量级以及湍流作用本身的不确定性，本文在平流范围的基础上将其边缘向外扩展一定的小范围作为修正，就可以近似视之为该线源的播云影响区范围。

1.2 CA-FCM方法

在国家“十五”科技攻关项目支持下，中国科学院大气物理研究所建立了基于聚类的浮动对比区历史回归人工增雨效果统计检验方法(简称CA-FCM方法)(房彬等，2006；翟雨等，2008)，根据历史雨量将站点进行聚类分析并在显著性检验指标内分为m类，飞机播云的影响区包含的站点分别属于其中的k个分区(1≤k≤m)，由此将总影响区分为k个分区，分别对每个分区挑选影响区以外历史相关性好的站点作为对比站点建立多元预报回归方程：

$ {\widehat V^{\left(k \right)}} = b_o^k + \sum\limits_{i = 1}^n {{a_i}^{\left(k \right)}{r_i}^{\left(k \right)} + \sum\limits_{i = 1}^m {{b_i}^{\left(k \right)}{w_i}^{\left(k \right)}} } $

(3)

式中，历史样本中副影响区的降水体积V_e^(k)为回归方程的因变量，对比站点的降水量r_i^(k)和其他物理协变量(如整层大气可降水量)w_i^(k)为回归方程的自变量。关于该方法的详细介绍可参见房彬等(2006)和翟雨等(2008)。

1.3 VB＋MO技术简介

MO全称为MapObjects，是ESRI公司的产品(薛伟，2004)，包括了一个ActiveX控件(OCX)也就是Map控件以及45个ActiveX自动对象，适用于符合IT行业标准的Windows编程环境。

VB+MO(Visual Basic与MapObjects)在地理信息系统二次开发当中虽然不是最先进的技术，但是，其开发构架简单，语言通俗易懂，非常适合快速开发，并适应业务人员的日常应用。通过VB可以非常方便地调用MO的各种功能，如图形的构建叠加与缓冲、投影变换以及面积计算等，并且可以动过调用动态链接(dll)结合Fortran语言进行复杂的数学运算。

2 飞机播云影响区计算 2.1 影响区范围

在一次实际的飞机作业过程中，飞机作业高度是变化的，因此，首先根据飞机飞行航迹的特征拐点将播撒路径分解成几条特征线段，并将每条特征线段近似为直线段，分别计算每条线段平流而成的四边形范围，再将它们叠加，并在其边缘添加一定的范围作为湍流扩散而成的范围，构成总影响区。

MO的基本显示功能不在此赘述，这里仅介绍下列关键步骤：

(1) 四边形构建

假设，在水平方向上，特征线源AB的催化剂在有效扩散时间内由高空风V平流到达A′B′位置。A′与B′的坐标分别依据A与B的坐标和风速风向计算得到，那么，用这4个点按照一定的顺序(顺时针或逆时针)可构造出一个四边形(图 1)。假定A、B、B′、A′的经纬度依次分别是lon(i)与lat(i)(i=1, 2, 3, 4)，相应的VB代码如下：

Dim poly As MapObjects2.Polygon

  %定义多边形对象

Dim pp As New MapObjects2.Point

  %定义点对象

Dim pts As New MapObjects2.Points

  %定义点集对象

For i＝1 to 4

  pp.x＝lon(i) : pp.y＝lat(i)

  pts.add pp  %将4个点依次加入点集

Next i

Poly..Parts.Add pts

  %由上述点集组成多边形

(2) 多边形叠加

假定按照上一步产生了n个四边形ploym(i) (i＝1, 2…, n)，将它们叠加至总影响区poly上，相应的VB代码如下：

For i＝1 to n

  Set poly＝poly.Union(polym(k))

    %叠加多边形

Next i

(3) 缓冲区操作

如果不计风速影响，催化剂靠湍流作用也能扩散，因此，需要在平流作用构成的影响区基础上增加一定范围的缓冲区。鉴于湍流扩散参数存在不确定性，为了简单起见，这里通过使其边缘向外扩展一定距离来表征湍流造成的扩散范围，具体代码如下：

Set poly＝poly.Buffer(ks)

  % ks代表缓冲半径

以吉林省2010年6月27日的飞机作业为例，计算得到的影响区如图 2所示。飞机作业高度3100 m，风速6.2 m·s^-1，风向105°。播撒路线为白城－四平－梅河口－九台－白城，有效扩散时间假定为3 h。

由图 2可见，影响区是由飞机路线分解成的4条特征线段经平流作用叠加构成的，与传统的主观划定影响区不同，图 2中的影响区是空心的，更加客观地表征了影响区的实际范围。

在实际作业中播撒高度可能起伏较大，有时还在不同高度进行播撒，且不同高度上风速风向变化较大，在这种情况下，可以根据上述方法将特征线源线段作细分，并分别计算每条播撒段的扩散范围后再叠加，以体现播撒高度变化对扩散范围的影响。

2.2 影响区面积与降水量

要建立对比站点对影响区体积降水量的历史回归方程，关键是要求得影响区更加准确的体积降水量，为此，这里先求出影响区的面积，再对影响区进行雨量插值，最后获得总的体积降水量。

(1) 多边形投影变换

前述步骤都是在以经纬度为单位的地理坐标下进行的，要求其面积，需要将地理坐标转换成投影坐标(假设前述步骤已经计算出影响区范围poly)，相应的代码如下：

Dim pcs As New MapObjects2.ProjCoordSys

  %定义投影坐标系

Dim pcs2 As New MapObjects2.GeoCoordSys

  %定义地理坐标系

Dim newpoly As MapObjects2.Polygon

  %定义一个新的多边形

pcs2.Type＝moGeoCS_WGS1984

  %原图地理坐标系

pcs.Type＝moProjCS_AsiaLambertConformalConic

  %设为MO预置的亚洲%兰伯特投影

Set newpoly＝pcs.Transform(pcs2, poly)

  %由poly经过投影得到%投影坐标下的新多边%形newpoly。

(2) 多边形面积计算

投影坐标系下多边形newpoly的面积S为：

式中，S的单位为m²。

(3) 影响区内格点划分

在规则的区域(譬如矩形)内，可以很容易布置均匀的格点以供雨量插值，但是，现在影响区是不规则的区域，因此，需要在更大的矩形范围内布置均匀格点，再通过MO规则判断出哪些格点在影响区范围内，进而得到影响区内均匀分布的格点数。对于某一格点(lon1, lat1)，相应的代码为：

Dim pt As MapObjects2.Point

  %定义MO点对象

Dim bl As Boolean

  %定义布尔型变量

pt.x＝lon1 : pt.y＝lat1

bl＝poly.IsPointIn(pt)

  %判断该点是否在影响区内

(4) 格点雨量插值

除了较容易实现的距离平方反比方法以外，VB还可以调用Fortran、Surfer等编程语言实现更为复杂的插值计算，以便得到影响区内分布密集的格点雨量。

(5) 体积降水量计算

降水量的观测值实际上是单位面积降水量，需要乘上对应的面积才能得到体积降水量。在前述步骤中，已经求得了影响区总面积S以及当中n个格点的雨量R_i (i＝1, 2, …, n)。由于格点是均匀分布的，当其足够密集时，可将S等分成n份：

$ {S_0} = S/n $

(4)

并近似认为R_i对应的面积是S₀，那么，影响区的体积降水量为：

$ V = \sum\limits_{i = 1}^n {{R_i} \cdot {S_0}} $

(5)

(6) 影响区分区的体积降水量计算

在CA-FCM方法中，要根据式(3) 对影响区中由聚类分析确定的每个分区分别建立回归方程，这就需要计算每个分区的体积降水量。假设影响区中包含m个分区，需要将格点也做相应划分，则第k个分区的体积降水量为：

$ {V^{\left(k \right)}} = \sum\limits_{i = 1}^{N\left(k \right)} {{R_i}^{\left(k \right)} \cdot {S_0}} $

(6)

式中N(k)为第k个分区的格点数(k＝1, 2, …, m)。

3 应用举例

将上述技术应用于中国科学院大气物理研究所云降水物理与强风暴重点实验室(LACS)开发的人工增雨效果统计评估技术系统(V3.2) 中，并计算吉林省2010年4月18日飞机播云增雨的效果。

4月18日受蒙古气旋影响，吉林省出现大范围降水云系，适合于开展飞机播云作业。吉林省人工影响天气办公室进行了飞机播云增雨作业，飞行播云航线为长春－长岭－通榆－乾安－农安－长春。该次播云催化剂为AgI焰条4根，液态二氧化碳46 kg。飞机播云时间是09:12－11:15，播云高度主要在3 600 m高度，根据当天的探空资料，相应的温度为-9℃，风速12.7 m·s^-1，风向208°。这次播撒作业高度基本能达到碘化银的播云温度窗要求，相应的播撒AgI活化成核率约为10¹³ g^-1。利用上述飞机播云不规则影响区计算方法得到影响区范围如图 3a所示。

根据吉林省50个雨量站点1997—2007年共11年的4月雨量数据，采用聚类分析方法对这些站点进行聚类分析，进而得到图 3中的5个聚类分区，各分区之间的相关系数见表 1。影响区中包含了第一、第三、第五这3个聚类分区的站点，据此将总影响区划分为3个分区，并分别根据相关性的大小依次挑选缓冲区以外的站点作为对比站点(如图 3b~3d)，譬如，第三分区的对比站点选择与第三分区相关性最好的前三类分区(3、4、5) 的站点，在历史资料中筛选出没有受作业影响并且有一定降水的数据作为历史样本建立多元逐步回归方程。然后，根据2010年4月18日飞机作业信息和降雨量等数据，可以得到各分区的实测降雨量、估计自然降雨量、绝对增雨量和相对增雨率，结果如表 2。对3个分区的结果进行累加，得到飞机作业增加的总降水量为30.92×10⁶ m³，相对增雨效果为11.24%，影响区面积为33425.15 km²。

从表 2注意到，总影响区中三个分区的增雨率差异较大，偏东部的第3分区的增雨率为44.85%，明显高于其他2个分区。下面从播云物理角度论证这个结果的合理性。

由作业时段的卫星红外云图(图 4)看到，在吉林大部地区上空有大片云层覆盖，并且在长春地区的云层比白城附近的云层要密实，说明长春地区(即第三分区所在地)播云时的云层发展较强。而理论研究(胡志晋，2001；Dinevich et al, 2005)指出了当云降水处于发展或持续阶段，云中有比较深厚的上升气流，云厚较大，过冷云层较厚，云底较低等人工增雨条件时有利于催化作业增雨。因此，第三分区很可能比第一分区拥有更有利的催化增雨条件，相应的增雨效果也更大。

由白城市雷达看到，4月18日上午，有一大范围的降水层状云系在白城南部自西向东移动，其云层底高大约几百米，顶高约5 km，云层的最大回波强度达到30~35 dBz(图略)。根据探空资料，云顶对应的环境温度约－20℃，0℃层高度约在1 km。飞机播云作业时段播云影响区正好位于该降水层状云系范围内，作业高度主要在3600 m，根据飞行观测记录，在作业时作业高度上出现少量积冰现象。这些观测说明该次飞机作业是在降水云层的过冷区中进行的，表明本文评估给出的播云正效果是有物理基础的，结果是合理的。

从上述白城雷达观测看到，飞机播云区位于持续发展的层状雨区中, 云顶、云底的高度、温度和过冷层厚度均符合播云增雨的条件。从卫星云图(图 4)也看到，长春地区附近的第三分区云层要比第一分区密实，说明那里催化增雨条件更有利，因此，第三分区的增雨效果应当比第一分区的高。

由于当日上午作业时段长春站雷达因故没有开机，因此，以上分析主要是从播云条件和云降水宏观特征来定性判断作业的有效性，如果同时有云微物理观测资料来进一步证明本文通过统计评估效果的合理将更有说服力。

4 结论与讨论

本文提出了更为客观地计算飞机作业影响区的适用方法，该方法将VB＋MO地理信息系统二次开发技术与二维线源扩散方程相结合，更加真实地划定出飞机播云的不规则影响区，更为精确地实现了不规则区域内体积降水量的计算，并且结合统计学方法，成功地应用于飞机人工增雨效果统计评估系统中实现了自动化。通过吉林省飞机增雨作业实例介绍了该方法的具体应用，便于进一步的学术探讨与业务应用。

最后，值得说明的是，本文所用的二维线源扩散模型在估算飞机播云催化剂扩散范围时具有一定的合理性，但是，这里仍然包含了一些假设，譬如，为简单起见，文中假定飞机飞行作业线源是瞬时形成的，并且不考虑作业区域的云层条件(如冰晶浓度、过冷液态水含量等)对线源扩散的影响，而这些假定与实际情况并不完全一致。在实际中飞机播云是按照连续源方式进行的，并且播云粒子还受自然云粒子捕获清除的影响。为此，下一步将根据更加接近实际的连续点源扩撒方程对飞机播云催化剂扩散范围进行计算，并引进云粒子的捕获清除作用，以获得更符合实际的播云扩散结果。此外，更为重要的是，今后还应致力于发展把统计评估与物理响应分析和人工增雨数值模拟结合的人工增雨效果综合检验评价方法，使统计评估结果更加令人信服。