2011, Vol. 37 
当前国内人工影响天气作业参数的制定，主要是依据天气雷达产品进行指标识别或寻找依据的，通过分析识别，来给出人工影响天气的作业时机、作业部位 (地面作业的方位和仰角)、用弹量等。从20世纪90年代开始，随着数字化天气雷达的普遍应用，通过雷达建立人工影响天气作业参数的业务开展，以此为核心，建立了人工影响天气的作业指挥系统，使人工影响天气作业进入了科学化、定量化、客观化的阶段[1-9]。
进入21世纪，全国S波段多普勒雷达布局逐步完成，全国大多数省份都均匀分布多个S波段多普勒雷达。目前的实时人工增雨作业指挥系统均采用的是单雷达参数指挥。由于单雷达在离开雷达较远的方向空间分辨率较低，垂直部位缺测点较多。因此通过雷达组网技术，建立各种插值方法[10-13]，获得各水平层的CAPPI资料，是区域型雷达组网技术应用发展的趋势。本文利用安徽省及其周边省的S波段多普勒雷达强度场资料，利用雷达拼图技术，对雷达回波重叠区进行插值，可以获得雷达在垂直方向上的缺测资料，结合火箭弹道曲线和用弹量公式，再加上实时作业参数的算法方法，可建立人影作业指挥系统，通过业务软件，在人影作业时可精细化进行作业指挥。
1 作业参数的系统特点和总体框架作业参数依托于安徽省火箭人影作业指挥系统进行分析。系统总体框架分为3个部分：雷达组网拼图、火箭弹道曲线和作业参数精细化。雷达组网拼图：安徽省及周边省8部S波段多普勒雷达进行组网，并用插值方法进行强度场资料的拼图，生成的图叠加在使用GIS数据处理过的地图上，可以进行缩放、移动等矢量处理；火箭弹道曲线：提供火箭的AgI含量、火箭运动轨迹，与雷达组网拼图综合分析，可生成实时作业参数；作业参数精细化：通过作业参数算法方法，考虑各种条件及限制，根据平均风场，可以精细化地对人影作业进行实时指挥。
2 GIS数据处理和雷达组网拼图方法 2.1 GIS数据及处理方法 2.1.1 数据来源GIS数据来源于国家基础地理信息中心的《全国1:25万地形数据库、数字高程模型》。坐标系统采用1954北京坐标系，高程基准采用1956黄海高程系。网格尺寸：3″×3″。等高线的等高距，东部平原丘陵区为50 m。
2.1.2 数据处理使用1:25万安徽省GIS边界数据，建立以110°E为参考轴的6°带的高斯正形投影。可以建立以任何作业点为中心的任何半径的矢量地图。误差小于1 m。
高斯-克吕格投影直角坐标的计算公式[14]如下：a为椭球体长半轴；b为椭球体短半轴; f为扁率 (a-b)/a; e为第一偏心率,
| $\begin{array}{l} {X_N} = {k_0}\left\{ {M + NtgB\left[ {\frac{{{A^2}}}{2} + (5 - T + 9C + 4{C^2})\frac{{{A^4}}}{{24}}} \right] + } \right.\\ \left. {(61 - 58T + {T^2} + 270C - 330TC)\frac{{{A^6}}}{{720}}} \right\}\\ {Y_F} = FE + {k_0}N\left[ {A + (1 - T + C)\frac{{{A^3}}}{6} + } \right.\\ \left. {(5 - 18T + {T^2} + 14C - 58TC)\frac{{{A^5}}}{{120}}} \right]\\ T = t{g^2}B\\ C = {e′^{2}}{\cos ^2}B\\ A = (L - {L_0})\cos B\\ M = a\left[ {(1 - \frac{{{e^2}}}{4} - \frac{{3{e^4}}}{{64}} - \frac{{5{e^6}}}{{256}})B - (\frac{{3{e^2}}}{8} - } \right.\\ \frac{{3{e^4}}}{{32}} - \frac{{45{e^6}}}{{1024}})\sin 2B + (\frac{{15{e^4}}}{{256}} + \frac{{45{e^6}}}{{1024}}) \times \\ \left. {\sin 4B - \frac{{35{e^6}}}{{3072}}\sin 6B} \right]\\ N = \frac{a}{{\sqrt {1 - {e^2}{{\sin }^2}B} }} = \frac{{({a^2} + b)}}{{\sqrt {1 - {e′^{2}}{{\sin }^2}B} }}\\ \end{array}$ |
上面公式中东纬偏移FE=500000 m+带号×1000000；高斯-克吕格投影比例因子k0=1；高斯-克吕格投影反解公式：(X, Y) →(B, L)，原点纬度0，中央经度L0。
| $\begin{array}{l} B = {B_f} - \frac{{{N_f}tg{B_f}}}{{{R_f}}}\left[ {\frac{{{D^2}}}{2} - (5 + 3{T_f} + {C_f} - } \right.\\ \left. {9{T_f}{C_f})\frac{{{D^4}}}{{24}}} \right]\left. { + (61 + 90{T_f} + 45T_f^2)\frac{{{D^6}}}{{720}}} \right]\\ L = {L_0} + \frac{1}{{\cos {B_f}}}\left[ {D - 1 + 2{T_f} + {C_f})\frac{{{D^3}}}{6}} \right] + \\ \left. {(5 + 28{T_f} + 6{C_f} + 8{T_f}{C_f} + 24T_f^2)\frac{{{D^5}}}{{120}}} \right]\\ {N_f} = \frac{{({a^2}/b)}}{{\sqrt {1 + {e′^{2}}{{\cos }^2}{B_f}} }} = \frac{a}{{\sqrt {1 + {e^2}{{\sin }^2}{B_f}} }}\\ {R_f} = \frac{{a(1 - {e^2})}}{{{{(1 - {e′^{2}}{{\sin }^2}{B_f})}^{3/2}}}}\\ {B_f} = \varphi + (3{e_1}/2 - 27e_1^3/32)\sin 2\varphi + (21e_1^2/16 - \\ 55e_1^4/32)\sin 4\varphi + (151e_1^3/96)\sin 6\varphi \\ {e_1} = \frac{{1 - b/a}}{{1 + b/a}}\\ \varphi = \frac{{{M_f}}}{{a(1 - {e^2}/4 - 3{e^4}/64 - 5{e^6}/256)}}\\ {M_f} = ({X_N} - FN)/{k_0}\\ {T_f} = {\rm tg}^2{B_f}\\ {C_f} = {e′^{2}}{\cos ^2}{B_f}\\ D = \frac{{{Y_F} - FE}}{{{k_0}{N_f}}} \end{array}\]$ |
上面公式中东纬偏移FE=500000 m；北纬偏移FN北半球=0，FN南半球=10000000 m。
然后将平面矢量地图进行平移和缩放，在屏幕指定大小的范围内显示 (1个像素点代表 1 km)。
| $\begin{array}{l} X - {X_0} = k(X′ - X′{_0})\\ Y - {Y_0} = k(Y′ - Y′{_0}) \end{array}$ |
(X，Y)、(X′, Y′) 为平移前的大地坐标以及平移后的屏幕坐标，k为放缩倍数。
2.2 组网雷达拼图方法 2.2.1 数据来源
|
图 1 合肥、阜阳、蚌埠、黄山、南京、徐州、驻马店、九江8部雷达覆盖范围示意图 Fig. 1 Areas covered by the 8 radars at Hefei, Fuyang, Bengbu, Huangshan, Nanjing, Xuzhou, Zhumadian and Jiujiang |
|
表 1 8部雷达的位置 Table 1 The location of the 8 radars |
把来自安徽省附近合肥、阜阳、蚌埠、黄山、南京、徐州、驻马店、九江8部雷达的反射率场插值到统一的笛卡尔坐标网格，8个雷达的格点反射率场拼接起来形成3D拼图网格。由于国家气象业务网中S波段雷达只有0.5°、1.45°、2.4°三个仰角的反射率场资料, 因此本文利用多普勒雷达的0.5°、1.45°、2.4°三个仰角的反射率场资料进行插值，这些反射率场资料已经进行了质量控制，分辨率都为1 km。根据人工增雨冷云催化剂的高度要求，形成2 km、3 km、4 km、5 km、6 km、7 km、8 km高度的反射率场格点数据并进行绘图。
|
图 2 VCP21扫描模式下合肥和阜阳雷达各仰角及观测高度 Fig. 2 Radar elevations and observation heights at Stations Hefei and Fuyang in VCP21 scan mode |
(1) 单个雷达在2~8 km的不同高度上选取上下最近的2个点进行插值。
| $dBz = dB{\mathcal{z}_1} \times {d_1}/({d_1} + {d_2}) + dB{\mathcal{z}_2} \times {d_2}/({d_1} + {d_2})$ |
式中
(2) 对于多雷达重叠区的拼图算法，国内普遍采取三种方式[15]：第一种是选取多部雷达中回波最大值作为重叠区域雷达回波值；第二种是根据不同雷达距离重叠区的距离权重来计算重叠区的雷达回波值；第三种是直接取多部雷达回波平均值作为重叠区的雷达回波值。由于本文采取的方案中各个雷达都已插值到相同水平面的相同格点，因此选取多部雷达中回波最大值作为重叠区域雷达回波最大反射率值的方法。
3 弹道曲线和用弹量 3.1 火箭作业的弹道轨迹由于碘化银的成冰阈温低于-5 ℃，所以作业的部位应在云中0 ℃层以上负温区，而且还要把碘化银送到含水量丰富、上升气流较大的地方。所以火箭发射仰角是由目标云作业部位高度计算出的作业点与目标云水平距离及火箭弹道飞行参数综合计算而确定。火箭弹道轨迹 (图 3) 结合雷达的强回波区与作业点的距离和高度可计算得到火箭发射仰角[16]。
|
图 3 WR-1B型火箭弹弹道曲线 Fig. 3 Ballistic curves of the WR-1B rocket |
对流云人工催化要达到动力催化效果，人工冰核浓度应达到300~500个/L[4]。而WR-1B型火箭每弹种装有10克AgI催化剂，在-10 ℃时冰核数可达1.8×1016个[17]。若作业对象为体积100 km3的对流云体，需3~5枚火箭可使云中冰核达300~500个/L。结合安徽多年实践经验和每次作业云的尺度及强度等特征：对流云用弹量选择3~5枚/次；层状云催化人工冰核浓度应达到20个/L，对射程内层状云作业需2~3枚/次。
4 精细化的作业参数单多普勒雷达，由于离地面远处探测资料稀少的原因，等高面CAPPI的高度一般取2 km左右，最高不超过4 km。但由于火箭人工增雨作业催化剂的播撒区域都在4 km以上，因此了解4 km以上高分辨率的拼图CAPPI对火箭人工增雨就很有必要。
作业参数算法：
(1) 由于安徽地处南方，作业时零度线大都在5 km左右，同时考虑弹道曲线，取作业点周围8 km范围内回波值，4 km和5 km两层的资料。
(2) 考虑作业安全，避开人口稠密区，对比作业点周围村庄分布，过滤掉有村庄的点。
(3) 根据安徽经验，大于30 dBz的回波有降水潜力。挑选出回波周围4个点都大于30 dBz的点，计算各点到作业点的方位角。
(4) 使用MM5数值模式的UV风场资料，对作业点周围8 km范围内的水平风场进行平均，得到平均风向。
(5) 对比各个点方位角和平均风向的角度差，根据迎风作业方法，按照角度差从小到大排序。
(6) 按照序列，选择某点计算和作业点的仰角，与火箭射击抛物线比较，如果符合，此点就是作业目标点。如果不符合，就依次查找下一个点，直到找到作业目标点。
图 4、图 5是2008年11月7日12时5分安徽组网雷达人工增雨防雹作业指挥系统软件界面截图。作业点是在安徽省宿州市符离镇作业点。可以看到界面右下方是目标作业云的火箭作业参数，并且以作业点为起点，绘制出了火箭作业轨迹，细线表示火箭上升阶段，粗线表示火箭播撒阶段。通过这种图形的方式，直观的表示出火箭发射的方位和作业不同阶段的距离。
|
图 4 2008年11月7日12时05分安徽组网雷达人工增雨防雹作业指挥系统软件界面截图 (4~5 km) Fig. 4 Diagram for technique parameters of weather modification given by the radar network of Anhui Province at 12:05 BT 7 November 2008 (4~5 km) |
|
图 5 2008年11月7日12时05分安徽组网雷达人工增雨防雹作业指挥系统软件界面截图 (6~8 km) Fig. 5 As in Fig. 4, but for heights of 6-8 km |
在4 km高度图上，绘出了火箭作业轨迹示意图，并标有作业指导参数：火箭仰角45°~70°，对该点云体火箭用弹量2~3枚。在5 km图上 (图 4)，绘有火箭作业轨迹示意图，并标有作业指导参数：火箭仰角50°~74°，对该点云体火箭用弹量2~3枚。
而6 km (图 5)、7 km、8 km的图，显示云体面积太小，不适合作业。
不同的高度，用不同的仰角；回波强度不同，用弹量也不一样。雷达拼图CAPPI分层显示的数据，结合火箭参数，可以精细化地获得作业指导参数。
5 总结(1) 多部雷达组网拼图可以较细致地反映离雷达较远处回波的垂直变化情况。
(2) 能够建立笛卡尔坐标平面中3D网格的反射率场，可以以水平和垂直方面的分辨率为1 km进行精细化的图像显示。
(3) 在可作业层能够精细化地进行计算并获得火箭作业参数。
(4) 雷达参数只考虑了反射率场的组网拼图，没有考虑VIL、回波顶高、风场等参量。
(5) 组网拼图的插值在远离雷达而没有雷达重叠区的地方，有较大的误差。