引言

我国人工增雨业务缺乏有效的分析可播区域的方法，各种探测资料及研究成果在业务中应用能力差。目前，常用的方法是建立叠加分析平台，叠加各种资料粗略地分析可播区域，如河南“九五”期间开发的人影决策指挥地理信息平台^[1], 这种利用平台叠加分析可播区域的方法使用起来不太方便，在实际业务中发挥的作用有限，并且当资料种类比较多、资料分辨率差异较大时，该方法将无法准确判断可播区域的位置。许多省、市人工影响天气业务技术系统中都具有可播区域分析功能，山东采用冰水转化区诊断、卫星云图处理、探空资料处理、物理量诊断及降水预报、云水预报以及雷达资料处理等模块来分析可播区域^[2]，然而各模块之间相对独立; 黑龙江应用MICAPS平台分析、查询天气和雷达资料以及其它相关子系统分析可播区域^[3]，这些系统都无法进行自动分析可播区域。针对这一问题，本文提出一种方法，解决如何利用雷达、卫星、数值模式产品和探空等不同类型和不同分辨率的资料，应用计算机软件技术自动选择可播区域。

1 使用单一资料选择可播区域的方法

为描述问题的方便，假定使用单一的雷达回波资料来确定可播区域。假设雷达回波图，长为L，宽为D(L和D的单位均为像素)，该图每个像素点的回波强度可知，并且若强度大于20dBz为可播(当然可以附加其它各种条件来判断可播性，但这对本文方法研究不产生影响)，由于满足条件的可播区可能分散地有许多块，并且每块大小和形状各异，如何应用计算机软件技术自动识别?

1.1 方法的设计思路

由于L、D的值较大，为了简化问题，可采用分块法进行处理。把雷达回波图分解成边长为a的正方形小块(a的大小可以根据实际需要确定)，形成一个小块矩阵，用数组P[x][y]表示小块的分布矩阵(x和y为小块在矩阵中的行、列号，其中0≤x＜L/a, 0≤y＜D/a)。根据前面给定的可播条件，假设r为可播点占小块所有像素点的比率，若r大于60%，则该块可播，令P[x][y]等于r，否则P[x][y]=0。一般情况下, 在整个区域中可播区范围较小，因此可以把矩阵P[x][y]当作稀疏矩阵处理，可建立正交链表来存储小块矩阵^[4]。查找可播区方法如下：遍历正交链表的结点，先找出一个未被访问的结点，把其插入到一个单链表中，同时把正交链表中该结点标识为已访问，然后从单链表的此结点开始，在正交链表中查找所有与该结点相邻的未被访问的结点，找出后把其标识为已访问并插入到单链表的结尾，按上面步骤依次处理单链表中结点，直到处理完毕。这样最后得到的单链表就是一块可播区。重复上面操作直到正交链表中所有结点均已访问，这样就找出了所有可播区域。

1.2 方法的实现 1.2.1 数据结构设计

(1) 稀疏矩阵数据结构的设计

本文采用正交链表存储稀疏矩阵。稀疏矩阵的每一行设置一个带表头结点的循环链表，每一列也设置一个带表头结点的循环链表，行链表和列链表的表头可以共用一个结点。定义类MatrixNode，链表中的所有结点都属于其对象，此类包含head、down和right域。Head域用于区分该结点是表头结点还是可播小块结点，down和right域用于存放稀疏矩阵列链表和行链表的最前端第一个结点的地址。此外，类MatrixNode还包括一个联合体，若head=1(即为表头结点)时，联合体取值为一个next的指针域，该域用于链结表头结点；若head=0(即为非表头结点)时，其取值为一结构体SeedingNode。结构体SeedingNode包含可播小块的行号、列号、最小可播高度、最大可播高度和播撒剂量等信息，另外还包括一个next指针域，其用于链结同一可播区中其它小块结点。类MatrixNode完整的定义见附录。

(2) 线性链表数据结构的设计

线性链表用来存储查找得到的可播区域。该链表结点定义为SeedingFieldListNode，其包括两个指针域，分别为link和right，其中right指向由同一可播区域所有小块组成的单链表的表首结点，link用来链结不同的可播区域。SeedingFieldListNode还定义了一个count域，用来存放可播区域小块总个数。线性链表的类的完整定义见附录。

1.2.2 建立正交链表的算法

建立正交链表过程中，设计一个长度为max{L/a, D/a}的辅助数组H，H[i]是一个指向第i个列链表和行链表的指针。在该算法中，首先建立所有的表头结点(开始时每个表头结点的next域指向它本身)，然后，按列读取矩阵的每个可播小块P[i][j]，定义一个指针last，每开始读取新的一列时，last指向H[j]，为每个读入的可播小块建立相应的结点，并把其链结到last->down后面，重新调整last，让其指向新链入的结点，以便链结后续读入的该列的其它结点。由于表头结点H[i]的next域中暂时存放着该行链表中最后一个结点的地址，我们可以方便地利用它在行链表中链结新结点, 但链结完新结点后还必须重新调整H[i]的next域，使其存放新链入的结点的地址。每读完一列，把last->down指向该列表头，以便封闭列链表。完成建立所有列链表后，还必须封闭所有的行链表。最后，使用next域把所有表头结点相互链结起来，即完成了正交链表的建立。建立正交链表的算法见附录。

1.2.3 可播区的求解算法

利用前面定义的数据结构和存储结构，可以非常方便的求解和存储可播区域。在求解算法中，采用一个双重嵌套遍历正交链表, 外层依次访问表头结点，内层依次访问列链表，每找到一个未访问的结点(该结点的Visited域的值为False)，则开始处理一块新的可播区域。在求解新的可播区域时，首先在线性表中建立一个新结点，把前面在正交链表中找到的未被访问的结点链入到线性表中新结点的right域后面，依次处理right域链结的单链表(开始时right后面只链结了一个结点)，找出与单链表结点所表示的小块行、列号差值的绝对值均小于等于1的未被访问的可播小块结点(即相邻的未被访问的小块)，然后把其链入单链表表尾。处理完成单链表中的所有结点，所求的单链表即为一块可播区。遍历完正交链表所有结点，得到的线性表就是所求的所有可播区域。可播区的求解算法见附录(略)。

2 综合利用各种资料选择可播区域的方法

在使用多种资料选择可播区域时，把前面数组P[x][y]的元素定义为结构体，这种对数据结构的修改，可给处理小块矩阵带来方便。

2.1 利用卫星资料选择可播区的方法

以GOES卫星为例，其分辨率大概9km左右，若雷达资料分辨率为1km，则划分小块时其边长a可取10，这样能保证每个小块至少有一个卫星资料的值点。在定义数组P[x][y]元素的结构体时可增设Satellite_Seeding域，用其存放该小块中根据卫星资料的值点(或反演产品)判断的可播性。如果某个小块包含多个值点，可取其平均值或根据可播值点占总的值点的比率来判断小块的可播性。最终确定小块可播性时，可根据各种资料判断的结果按一定加权系数来确定。确定小块可播性后，就可以利用上面的算法求解所有可播区域。

2.2 利用数值模拟产品和探空资料选择可播区域的方法

在利用数值产品和探空资料选择可播区域时，首先找出数值产品格点或探空站所处的小块，根据格点资料或探空资料判断该小块的可播属性，再由此确定其相邻小块的可播属性。具体方法为：在设计数组P[x][y]元素的结构体时可增设七个域，分别用来存放根据数值产品判断的可播性、可播最小高度、可播最大高度和播撒剂量以及根据探空资料判断的可播性、可播最小高度和可播最大高度。在使用数值产品选择可播区域时，假设资料的水平格距为30km，若格点所处的小块可播，则与该小块行、列号差值的绝对值小于等于1的8个小块均可播并且可播高度和播撒剂也相同。对于探空资料，可以插值到每个小块，然后判断小块可播性。最后，综合各种单一资料判断的结果，得出小块可播属性，这样就可以根据前面描述的算法求解所有可播区域。

至于使用更多的其它资料，其方法与上面的分析类似。

3 该方法在飞机和高炮(火箭)人工增雨中的应用

用前面算法求解得到的所有可播区域以线性表形式存储。线性表结点的count域记录了每块可播区域的小块数量，根据count值的大小可直接计算该可播区面积(count×a²)。根据飞机性能指标、可催化时间以及机场与可播区之间的距离等因素，可自动选择一个最佳大小和位置的可播区域。把存储该最佳可播区域的单链表按第一、二关键字分别为row、col进行排序，经处理可得到该可播区边界，可以直接显示可播区域，并能详细标明各小块可播高度以及播撒剂量等信息。根据最佳可播高度的风向风速情况，能够使用计算机自动预设航线，并能详细标明每段轨迹的播撒高度以及轨迹长度等信息，这样就可以直接用于指挥飞机进行增雨作业。

对于使用高炮(火箭)人工增雨，根据可播区信息和炮点的位置，可直接确定进行增雨作业的炮点，并能根据小块可播高度确定作业的仰角。

4 个例应用

以2004年5月2日的郑州714CD雷达资料和探空资料为例，讨论该方法初步应用的一个个例。雷达取14：56的资料，可播性判据为强度大于20dBz。探空取早晨8点的资料，可播性判据为：T在-5℃至-15℃之间; T-Td≤2℃; e-Ei≥-0.1。

把整个区域划分为34×34的矩阵，每个小块边长为14，像素点分辨率为1km。按可播点占小块所有像素点的比率大于60%来判断小块可播性，根据雷达回波资料，共获得可播小块68块。对于探空资料，由于整个河南省只有两个探空站[郑州站(位置在东北方)和南阳站(位置在西南方)]，资料分辨率较粗，在使用其确定小块可播性时，我们以经过这两个站的中心点，斜率为-1的直线为分界线，小块左上角在直线上方(包括在直线上)时，使用郑州站探空资料来计算，在直线下方时，使用南阳站探空资料来计算。这样结合前面的探空资料可播性判据，所求得的68个可播小块播撒高度均在5.23km与6.45km之间。建立正交链表，并使用上述求解算法共求得大小不等的7块可播区域，所求线性表结构图如图 1，图中线性链表结点count域为right域所指可播区小块个数，link域指向下一可播区表头结点，right域指向该可播区的第一个小块。描述小块的结构体中col、row为小块所在矩阵的行、列号，amount域用来存放播撒剂量(本例没有计算播撒剂量), Height_Botton、Height_Top域用来存放该小块可播的最小高度和最大高度。从图中可以看出，7块可播区包含的小块数量分别为：6、30、4、9、15、3和1(按求解所得的顺序排列)，输出所求可播区域边界如图 2。

若使用飞机增雨，飞机飞行时速平均约360km·h^-1，一架次携带的催化剂量可催化1~2小时，则7块可播区中只有第二块可播区满足条件(主要从可播区大小考虑)，其可播面积约为5880km²。若播撒轨迹与风向垂直，由探空资料可知，高度在5.23km~6.45km时，风向为220°左右，则播撒轨迹斜率应为tan(2π-220π/180)。在预设可播轨迹时，通过与可播区上边界和左边界相邻的小块的左上角(若轨迹斜率为正值，则需通过与可播区上边界和右边界相邻的小块的右上角)，且以tan(2π-220π/180)为斜率作直线，则直线在该可播区内的部分即为所求的可播轨迹，相邻两条轨迹间距约为10.7km[acos(220π/180-π)]，输出所有的可播轨迹可得到飞机人工增雨的预设航线(如图 3)。按此航线进行播撒，其总的可播长度约为545km(总可播长度为前面所求可播线段的和，而可播线段长度d由表达式①、②求得), 飞机播撒约耗时1.5小时左右。

$\begin{array}{l} {\rm{Temp}} = \sin \left( {{\rm{lat}}1{\rm{ \mathit{ π} }}/180} \right)\sin \left( {{\rm{lat}}2{\rm{ \mathit{ π} }}/180} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\cos \left( {{\rm{lat}}1{\rm{ \mathit{ π} }}/180} \right)\cos \left( {{\rm{lat}}2{\rm{ \mathit{ π} }}/180} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\cos \left( {{\rm{lon}}2{\rm{ \mathit{ π} }}/180 - {\rm{lon}}1{\rm{ \mathit{ π} }}/180} \right) \end{array} $

(1)

$d = {\rm{acos}}\left( {{\rm{temp}}} \right){\rm{R}} $

(2)

其中(lat1, lon1)和(lat2, lon2)分别是每段可播轨迹的两个端点，R是地球半径，取值为6370km。

在该个例中，使用本文方法选取的飞机增雨可播区回波较强，并且连续成片，说明该区域有较大范围的稳定的层状云或混合云降水，其常常具有较大的增雨潜力。同时，从当日14~20点的雨量资料看，该个例预设航线区域(豫北)普遍降水在4~8mm，其降水量明显高于其它地区，且降水较均匀稳定，这说明该区域较有利于增雨作业。

5 结语

本文先通过讨论利用单一雷达回波资料，使用分块法，把整个区域划分为小块矩阵，确定各小块可播性，然后把小块矩阵作为稀疏矩阵来处理，采用正交链表存储方式，使用相邻搜索算法，自动找出所有可播区域，并详细描述了该方法的实现算法以及采用的相应数据结构和存储结构。然后解决了综合使用雷达、卫星、数值模式产品和探空等各种不同类型、不同分辨率的资料如何自动选择可播区的问题，提出了对于分辨率小于小块尺度的资料直接用与小块相应的值点的平均值判断小块的可播性，对于分辨率大于小块尺度的资料采用由资料值点确定与之相邻的小块可播性。该方法实现了根据各种资料利用计算机软件技术自动选择可播区，并直接用于指挥飞机和高炮(火箭)人工增雨，其对人工影响天气业务现代化建设具有重要的意义。文章最后介绍了使用该方法的一个初步应用的例子，其结果表明该方法选择可播区域较为合理。