引言

Fujita(1979)根据国际上第一次针对下击暴流开展的外场观测计划(the Northern Illinois Meteorological Research on Downburst)观测到的50个下击暴流的统计特征，正式定义了下击暴流：雷暴天气中在地面或近地面附近产生风速＞17.9 m·s^-1的灾害性辐散风的强下沉气流。Fujita(1985)按下击暴流地面出流的灾害范围大小，把下击暴流分为宏下击暴流(＞4 km)和微下击暴流(≤4 km)。Wilson et al(1984)基于多普勒雷达的观测对微下击暴流进行了补充定义：穿过辐散中心的多普勒速度差≥10 m·s^-1，初始的最大入流和出流中心(正负速度对)之间的距离＜4 km。Hjelmfelt(1988)把下击暴流按出流形态分为孤立的微下击暴流和下击暴流线。

多普勒天气雷达是监测预报下击暴流的重要手段。下击暴流在多普勒天气雷达观测的低仰角径向速度中表现为小尺度径向辐散特征(雷达径向上的正负速度对)或大风核区(因外部环境因子影响和离散的小尺度下沉气流叠加引起的非对称性出流)(Fujita，1979；Wilson and Schreiber, 1986；Potts，1989)。弓状回波也是一种在较大范围内产生较强下击暴流的回波形态(Fujita，1979)。Roberts and Wilson(1989)使用单个和多个多普勒天气雷达研究了一些云内下击暴流发生的先兆，他们确定了在下击暴流(初始地面出流)发生前2~6 min普遍出现的4个主要特征：①下降的反射率因子核；②云底附近或之上高度增加的径向辐合；③沿风暴边缘腐蚀的反射率因子(反射率因子槽口)；④绕垂直轴的旋转。

基于对下击暴流雷达特征的一些认知，下击暴流识别和临近预报方法先后被研发。早在20世纪80年代末和90年代初，Merritt(1987)基于下击暴流在低仰角径向速度中的小尺度径向辐散特征，研发了一些下击暴流的识别算法。由于这种特征只出现在距地面1 km以内的高度，雷达用最低扫描仰角(一般为0.5°)探测到这种特征的有效半径也只有60 km左右，而且当雷达探测到这种特征时，地面下击暴流大风即将或正在发生，几乎没有预报时效。因此在认知下击暴流雷达先兆特征的基础上，一些下击暴流临近预报算法被提出。最初的算法仅基于下击暴流在雷达反射率因子及其反演产品中的先兆特征(例如风暴质心高度和垂直液态水含量VIL的快速下降)而设计(Roberts and Wilson, 1989；Wolfson et al，1994)，预报成功率非常有限。Smith et al(2004)综合使用雷达反射率因子和径向速度等数据得到风暴单体的26个雷达特征量来预报下击暴流，该算法(Damaging Downburst Prediction and Detection Algorithm)在美国的WSR-88D系统中运行。

我国新一代多普勒天气雷达网的建设始于1998年。俞小鼎等(2006)首次利用我国新一代多普勒天气雷达资料对2003年6月6日发生在安徽的一次系列下击暴流过程进行了详细分析。在这之后，有许多气象学者基于多普勒天气雷达观测到的下击暴流资料进行了个例分析，例如：毕旭等，2007；吴芳芳等，2009；程月星等，2018；郭英莲和孙继松，2019；沈杭锋等，2019；龙柯吉等，2020；俞小鼎和郑永光，2020。杨璐等(2018)分析了北京地区雷暴大风不同生命期内的雷达统计特征及预警提前量。

我国的一些气象工作者也在下击暴流的识别和临近预报技术方面开展了一些研发工作。基于低仰角径向速度中的小尺度径向辐散特征的下击暴流识别算法被研发(张钢等，2011；陶岚和戴建华，2011；杜牧云等，2015; 2017)。李国翠等(2014)利用雷达反射率因子三维拼图数据开展了雷暴大风识别技术研究。罗辉等(2015)依据反射率因子核的快速下降和中层径向辐合特征，设计了下击暴流出流强度公式。王兴等(2019)基于风暴核顶高的快速下降和中层径向辐合特征来实现下击暴流的智能预报。王萍和牛智勇(2014)和肖艳姣(2018)使用不同的方法研发了中层径向辐合(mid-altitude radial covergence, MARC)特征自动识别算法。周康辉等(2017)在地面气象观测站大风记录的基础上，结合多源数据，利用模糊逻辑算法，实现雷暴大风与非雷暴大风的有效识别，可对雷暴大风进行实时监测。孙京等(2019)通过云模型和支持向量机等方法从诸多风暴单体的雷达特征量中筛选出9个下击暴流预报因子，利用Bayes和BP神经网络两种方法建立了下击暴流临近预报模型。

本文提出一种综合使用雷达和探空观测资料的下击暴流临近预报算法。在对雷达基数据进行地物杂波抑制和径向速度退模糊以及对探空资料进行处理得到0℃、-20℃和最小相当位温高度的基础上，该算法首先进行风暴单体识别追踪和冰雹指数的计算；然后进行MARC特征和中气旋识别，并使之与识别的风暴单体相关联；最后提取诸多风暴单体的雷达特征量，经批量下击暴流和非下击暴流个例统计分析后挑选出9个下击暴流的雷达先兆因子作为模糊逻辑法的输入，建立下击暴流临近预报方程。

1 资料和预处理

文中用于风暴单体雷达特征量统计分析的资料包含2010-2016年的CINRAD-SA/SB(China New Generation Weather Radar-S Band A/B Type) 型雷达观测的基数据(武汉、襄阳、荆州、宜昌、恩施、十堰和随州)、SWAN系统(Severe Weather Automatic Nowcast System)雷达拼图数据(27°~35°N、105°~117°E)、湖北省区域自动气象站观测的地面风数据和探空观测数据(安康、南阳、恩施、宜昌、武汉、长沙)；用于2015年6月1日个例分析的雷达基数据来自岳阳雷达，探空资料来自当日20时的长沙站；用于评估检验的是2019年6-8月湖北省的雷暴大风过程个例资料(包括雷达基数据、区域自动气象站地面大风数据和探空数据)。

雷达基数据采用VCP21体扫模式(6 min扫描9个固定仰角，分别为0.5°、1.5°、2.4°、3.4°、4.3°、6.0°、9.9°、14.6°、19.5°)采样，雷达波束宽度约为1°，反射率因子数据库长为1 km，径向速度数据库长为250 m。对雷达基数据的预处理包含地物杂波抑制(吴涛等，2013)和径向速度退模糊(肖艳姣等，2012；2016)处理。

下击暴流样本的选取结合地面测风数据和雷达观测数据，范围限定在距离雷达10~150 km内，主要考虑两个因素: 一是近距离处雷达探测不到风暴的中上部信息; 二是新一代天气雷达的径向速度最大探测距离约为150 km; 虽然通过距离退折叠后，其探测范围可扩展到230 km，但150 km以外的有效速度信息很少。针对地面大风数据，首先剔除了 < 17.9 m·s^-1的数据；然后用概率统计分析方法剔除因地形和设备异常等原因引起的大风站点，即剔除大风出现频次远大于站均频次的站点；最后结合SWAN系统的雷达拼图数据剔除冷空气大风和脱离风暴母体一定距离的阵风锋大风，具体方法为：对于出现了大风的站点，如果其周围40 km范围内有面积大于100 km²的40 dBz以上的回波和面积大于40 km²的50 dBz以上的回波，则认为该站是因雷暴引起的大风，否则被剔除。将保留的地面大风与大风发生前6 min最低仰角质心距离站点10 km范围内的风暴进行匹配，能与风暴单体匹配上的地面大风被确定为下击暴流大风。

对于探空观测数据，首先计算各高度的环境相当位温(刘健文等，2005)，找到最小相当位温所在高度，其次使用线性插值处理获取0℃、-20℃的高度。

2 下击暴流临近预报方法

图 1给出了下击暴流临近预报算法流程图，其中下击暴流临近预报方程是在使用批量下击暴流和非下击暴流个例对提取的雷达特征量进行统计分析后使用模糊逻辑法建立的，下面对几个主要模块分别进行介绍。

2.1 风暴单体识别、追踪和冰雹指数计算

使用SCIT(storm cell identification and tracking)算法识别风暴单体(Johnson et al，1998)，但对其中的反射率因子阈值进行了调整，只使用了6个反射率因子阈值，分别为60、55、50、45、40、35 dBz，其目的是为了提高强风暴单体的识别率。

使用冰雹探测算法(hail detection algorithm，HDA)(Witt et al，1998)计算强冰雹指数(SHI)和强冰雹概率(POSH)，公式如下：

$ S H I=0.1 \int_{H_{0}}^{H_{\mathrm{T}}} W_{\mathrm{T}}(H) E \mathrm{~d} H $

(1)

$ P O S H=29 \ln \left(\frac{S H I}{W T}\right)+50 $

(2)

其中，

$ E=5 \times 10^{-6} \times 10^{0.084 Z} W(Z) $

(3)

$ W_{\mathrm{T}}(H)= \begin{cases}0 & H \leqslant H_{0} \\ \frac{H-H_{0}}{H_{\mathrm{m} 20}-H_{0}} & H_{0}<H<H_{\mathrm{m} 20} \\ 1 & H \geqslant H_{\mathrm{m} 20}\end{cases} $

(4)

$ W T=57.5 H_{0}-121 $

(5)

$ W(Z)=\left\{\begin{array}{lc} 0 & Z \leqslant Z_{\mathrm{L}} \\ \frac{Z-Z_{\mathrm{L}}}{Z_{\mathrm{U}}-Z_{\mathrm{L}}} & Z_{\mathrm{L}}<Z<Z_{\mathrm{U}} \\ 1 & Z \geqslant Z_{\mathrm{U}} \end{array}\right. $

(6)

式中：Z_L=40 dBz，Z_U=50 dBz，H₀和H_m20分别为0℃和-20℃高度，Z为雷达反射率因子，单位为dBz。

E为冰雹动能通量，单位为J·m^-2·s^-1。W_T(H)为基于温度的权重函数，W(Z)为用于定义降雨和冰雹反射率因子转换区的权重函数，WT为与零度层高度有关的冰雹警报阈值，单位为J·m^-1·s^-1。

对于识别出来的风暴单体，如果其VIL < 5 kg·m^-2则被舍弃，余下的风暴单体根据VIL值从大到小排序，用字母A~Z加数字0~19(比如A0，B19)标记识别出来的风暴单体。

基于移动路径最短原则对相邻体扫时次的风暴单体进行匹配，其中使用最大反射率因子、风暴单体底高、移动速度和方向作为约束，即相邻体扫时次的同一风暴单体的最大反射率因子的变化不能超过12 dBz，移动速度不能超过120 km·h^-1，移动方向与前一体扫时次所有风暴单体的平均移动方向的偏差不能大于60°，底高变化不能超过2 km。Han et al(2009)讨论了移动速度动态阈值S_max选择方法，根据风暴回波面积大小设置了三组阈值，当面积 < 300 km²时，S_max=100 km·h^-1；当300≤面积 < 500 km²时，S_max=150 km·h^-1；当面积≥500 km²时，S_max=200 km·h^-1。由于本文是识别风暴核，面积不是特别大，故设置S_max=120 km·h^-1。在匹配过程中，对风暴合并和分裂进行了处理。如果当前时次的一个风暴单体同时匹配上前一时次的两个风暴单体，即出现风暴合并，那么把最大反射率因子变化小的前一时次的风暴单体的标号赋给当前时次的合并单体，前一时次的另一个单体就被认为消散了；如果当前时次的两个风暴单体同时匹配上前一时次的同一个单体，即出现了风暴分裂，则把前一时次的风暴单体序号赋给当前时次最大反射率因子变化小的风暴单体，另一个单体当成新生单体，赋予新的标号。

2.2 径向散度和方位涡度切变估计

使用二维局部线性最小二乘导数(LLSD)方法估计径向散度和方位涡度切变(Elmore et al，1994)。在极坐标下，取一个(2N+1)×(2M+1)大小的二维局地窗口，设窗口内各方位-斜距库的径向速度为u_i，j(i=-N，…，0，1, …, N；j=-M，…，0，1, …, M)。以窗口中心为基准(i=0，j=0)，沿径向中心以远，i为正，沿切向顺时针方向，j为正，则窗口中心的径向散度切变u_r和方位涡度切变u_a分别为：

$ \begin{gathered} u_{\mathrm{r}}=\frac{\sum\limits_{j} \sum\limits_{i} i u_{i, j} w_{i, j}}{\Delta r \sum\limits_{j} \sum\limits_{i} i^{2} w_{i, j}} \\ (i=-N, \cdots, 0, 1, \cdots, N ; j=-M, \cdots, 0, 1, \cdots, M) \end{gathered} $

(7)

$ \begin{gathered} u_{\mathrm{a}}=\frac{\sum\limits_{j} \sum\limits_{i} S_{i, j} u_{i, j} w_{i, j}}{\sum\limits_{j} \sum\limits_{i} S_{i, j}{ }^{2} w_{i, j}} \\ (i=-N, \cdots, 0, 1, \cdots, N ; j=-M, \cdots, 0, 1, \cdots, M) \end{gathered} $

(8)

式中：S_i，j=(r+jΔr)·i·Δφ，r为局地二维窗口中心离雷达的斜距(单位：km)，Δr为径向速度数据库库长(单位：km)，Δφ为波束宽度(单位：弧度)。w_i，j为权重，对于所有的i和j，都应有w_ij>0，w_i，j=w_-i，j，w_i，j=w_i，-j，本文取w_i，j=1。如果速度场中有径向辐合，则u_r < 0，u_r越小，径向辐合切变越大；如果速度场中有气旋性旋转，则u_a>0，u_a越大，正方位涡度切变越大。

由于使用LLSD方法计算径向散度和方位涡度切变时使用的是极坐标径向速度数据，其方位分辨率随着离开雷达的距离r而变化，故所取局地窗口大小也随着r变化。设径向辐合核和中气旋核的大小是某个常数D(缺省1.25 km)，那么N，M的取值分别为：

$ N=(\text { int })\left[\left(\frac{D}{\Delta r}-1\right) / 2\right] $

(9)

$ M=(\text { int })\left[\left(\frac{D}{r \Delta \varphi}-1\right) / 2\right] $

(10)

式中：(int)表示四舍五入取整。如果M < 1，则令M=1；如果M>13，则令M=13。

使用二维局部LLSD方法估计径向散度和方位涡度切变后再基于这两个参量去识别MARC特征和中气旋，相比直接使用相邻距离或方位库的速度差的识别方法，前者可以减轻径向速度的质量和小尺度自然脉动对识别的影响。

2.3 MARC特征和中气旋识别

肖艳姣(2018)给出了基于极坐标三维径向散度切变数据的MARC特征自动识别算法(storm cell radial convergence zone identification，SCRCZI)。本文使用的中气旋识别算法与SCRCZI算法类似，区别在于中气旋识别是基于极坐标三维方位涡度切变数据，识别也分为三个部分，即一维正方位涡度切变段、二维正方位涡度切变特征和三维正方位涡度切变特征(以下简称为一维切变段、二维切变特征和三维切变特征)的识别，步骤如下：

(1) 一维切变段识别

一维切变段的识别是沿径向方向搜寻方位涡度切变大于预设阈值(缺省分7级，分别是70、60、50、40、30、20、10，单位为10^-4s^-1)的连续库。当沿径向从雷达处开始往外第一次搜索到大于阈值的方位涡度切变库时，那么其后面大于该阈值的连续方位涡度切变库被聚在一起，直到小于该阈值的方位涡度切变库被搜索到；如果随后库的方位涡度切变值与该阈值的差≤方位涡度切变差阈值(缺省值为4，单位为10^-4s^-1)，且满足条件的连续库的数量小于中断计数阈值(缺省值为2个)，那么该切变段的搜索继续，否则切变段截止在第一个小于阈值的库之前。使用这样的搜索方法可以减轻速度的小尺度自然脉动对一维切变段搜索的影响。如果切变段的长度大于段长度阈值(缺省值为0.95 km)，那么该段就被保存，否则被剔除。每个切变段的开始和结束的距离库、方位、最大和最小速度等参量被记录。

(2) 二维切变特征识别

当一个仰角的最后一根径向的一维切变段被搜索完后，基于空间邻近原则，一维切变段被合并为二维切变特征。相邻段的合并需满足两个标准：一是2个段的方位差小于方位间隔阈值(缺省值为1.5°)，二是2个段在径向上的重叠距离要大于段重叠阈值(缺省值为0.45 km)。一个二维切变特征包含的一维切变段数量必须大于段数阈值(缺省值为2)，其面积要大于二维切变特征面积阈值(缺省值为3 km²)。在完成所有等级阈值的二维切变特征识别后，为了提取最强方位涡度切变区信息，需要舍弃包裹着高阈值二维切变特征的低阈值二维切变特征。如果高阈值的二维切变特征的中心落在低阈值二维切变特征面积内，那么低阈值的二维切变特征被舍弃，但是在被保留的最高等级阈值的二维切变特征中增加记录包裹它的最低等级阈值的二维切变特征的参量信息。对二维切变特征做尺寸大小、尺寸对称性和切变对称性检查，即使用椭圆拟合二维切变特征，计算其长、短轴半径a和b以及最大、最小速度库的连线与连线中点所在雷达径向的小夹角α。当a≤15 km、a/b≤4、45°≤α≤135°时，该二维切变特征被保留，否则被舍弃。计算并保存二维切变特征的一些参量，其中包含中心位置(斜距、方位和仰角)、开始和结束的方位和斜距、最大方位涡度切变、最大最小速度及其所在的方位和斜距及对应的差分方位切变等。

(3) 三维切变特征识别

把保留的二维切变特征按最大方位涡度切变大小排序，之后进行垂直关联。每个被识别的三维切变特征至少包含连续仰角的两个二维切变特征。垂直关联是从最低仰角开始的一个反复过程。首先垂直关联相邻仰角的中心距离小于阈值(缺省值2.5 km)的二维切变特征，如果有多个二维切变特征可以关联，那么只选择垂直积分方位涡度切变最大的那个二维切变特征进行关联。如果第一次垂直关联结束后还有未被垂直关联的二维切变特征存在，那么把搜索半径增加到5 km，对未被垂直关联的所有二维切变特征重复上面的步骤。如果第二次垂直关联结束后仍然有未被垂直关联的二维切变特征存在，那么把搜索半径增加到7.5 km后再次重复上面步骤进行第三次垂直关联。

如果以三维切变特征为中心的20 km半径内没有风暴单体存在，那么该三维切变特征被删除。余下的三维切变特征被认定为中气旋，之后计算中气旋的一些特征参量，包含垂直积分最大方位涡度切变、6 km以下的最大方位涡度切变和旋转速度(最大正负速度绝对值之和的一半)等，把中气旋按垂直积分最大方位涡度大小排序。

2.4 风暴单体雷达特征参量提取

用距离最近原则，把识别出来的MARC特征和中气旋与风暴单体相关联，然后提取风暴单体的一些雷达特征量。从风暴识别算法中可提取风暴单体质心位置、顶高、底高、VIL及其密度、最大反射率因子及其所在高度、最小相当位温高度附近及其以上高度的最大反射率因子、1 km高度以下的最大绝对径向速度和最大径向辐散等；从冰雹探测算法中可提取强冰雹指数SHI和概率POSH等；从风暴追踪算法中可提取风暴移动速度和方向、VIL的变化、风暴顶高和质心高度变化；从MARC特征识别算法中可提取MARC特征的质心位置、顶高和底高、垂直积分径向辐合值、最强径向辐合分量的强度等级、1~6 km高度的最大径向辐合值和径向方向上的最大最小速度差、最小相当位温高度附近的最大径向辐合值和径向方向上的最大最小速度差等；从中气旋识别算法中可提取中气旋的质心位置、顶高和底高、直径、最大旋转速度及其所在高度、最大方位切变及其所在高度、最低4个方位涡度切变分量的旋转速度和方位切变等。其中有一些特征参量用到了基于探空资料计算的0℃、-20℃和最小位温高度，探空资料来自最接近雷达体扫时间且距离风暴质心最近的探空观测。需要说明的是：与MARC特征和中气旋有关的一些雷达特征量本可以通过识别风暴单体后直接从对应的极坐标径向速度、三维径向散度切变和方位涡度切变数据中计算得到，但是考虑到风暴单体的识别是尽可能识别强风暴核，其对应的径向速度场中可能不包含完整的MARC特征或中气旋结构(MARC特征可能部分位于风暴后侧中层的入流槽口中，中气旋核区可能部分位于弱回波或有界弱回波区内)，所以本文采取先识别MARC特征和中气旋后再提取相关雷达特征参量的方式。

使用2010-2016年湖北省所有新一代天气雷达观测的降水天气过程提取了下击暴流和非下击暴流风暴单体的18个雷达特征量(表略)，并进行统计分析。对于下击暴流风暴单体，选取地面大风发生之前最多1小时的风暴单体雷达特征量；对于非下击暴流风暴单体，选取其整个生命史的雷达特征量。选取较长时间的雷达特征量的目的是为了延长下击暴流预报时效。结果表明，下击暴流和非下击暴流风暴的这18个特征量都有重叠隶属部分，筛选出重叠隶属部分相对较少的9个雷达特征量作为下击暴流预报因子(见表 1)。此外，分析了产生下击暴流的风暴单体核心高度、顶高和VIL的时间变率，发现这三个量的快速下降通常发生在地面大风发生之前的一个体扫，使用这几个量预报下击暴流并没有多少预报时效，所以没有被作为候选预报因子。

2.5 下击暴流临近预报方程建立

本文使用模糊逻辑法建立下击暴流临近预报方程，公式(11)给出了产生下击暴流的风暴单体的雷达特征量x的隶属函数F(x)，其中T_U和T_L分别为高、低临界值。

$ F(x)=\left\{\begin{array}{lc} 0 & x \leqslant T_{\mathrm{L}} \\ \frac{x-T_{L}}{T_{\mathrm{U}}-T_{\mathrm{L}}} & T_{\mathrm{L}}<x<T_{\mathrm{U}} \\ 1 & x \geqslant T_{\mathrm{U}} \end{array}\right. $

(11)

把雷达特征量作为输入变量，使用模糊逻辑法建立下击暴流临近预报概率方程：

$ P=\frac{\sum\limits_{i=1}^{n} F\left(x_{i}\right) W\left(x_{i}\right)}{\sum\limits_{i=1}^{n} W\left(x_{i}\right)} $

(12)

式中：P为下击暴流发生的概率，F(x_i)和W(x_i)分别为雷达特征量x_i的隶属函数和权重，n为输入变量个数。当某风暴单体的P值大于预设阈值(缺省值为50)时就认为该风暴单体将在其余下的生命期内产生下击暴流。

表 1给出了用于下击暴流临近预报的风暴单体的各雷达特征量隶属函数的高、低临界值和权重，其中隶属函数的高、低临界值是使用2010-2016年发生在湖北省的下击暴流个例和非下击暴流个例统计分析各雷达特征量的概率分布情况而得到的，权重是根据使用单个雷达特征量预报下击暴流的击中率得到的。

3 个例分析

2015年6月1日21:28(北京时，下同)左右，“东方之星”客轮在长江航道湖北监利段开始侧翻，21:31倾覆，导致442人死亡。两个气象专家组进行了现场调查和多种探测资料综合分析，得出相似结论。孟智勇等(2016)研究表明，客轮在倾覆时处于一条飑线弓形回波顶点处，根据树木折断情况估计遭遇了至少31 m·s^-1的强风，该强风可能是由微下击暴流里的直线大风或其中的小尺度涡旋(非龙卷)造成。郑永光等(2016)研究表明，客轮倾覆江段和周边区域发生了显著微下击暴流导致的强烈大风灾害，估计最强风力超过12级，并具有空间分布不连续、多尺度和强灾害时空尺度小等特征。事发周边区域北部的四台村养猪场附近树林中同时发生了多条相邻的微下击暴流条迹，呈现出辐散和辐合交替分布的特征。由于微下击暴流尺度非常小，现有地面自动气象观测站网没有捕捉到这次微下击暴流大风过程。

本文使用岳阳雷达观测的体扫资料分析了该下击暴流临近预报算法的预报效果，客轮开始侧翻的地点在岳阳雷达极坐标中的位置为(339.3°、49.1 km)。使用该算法从北京时19:31(体扫开始时间，下同)的体扫资料开始处理，从长沙20时的探空资料获取0℃、-20℃和最小相当位温的所在高度分别为5.2、8.7和6.4 km(长沙探空站与客轮侧翻点之间的距离为180 km)。算法从20:11-21:09连续11个体扫识别到标号为A17的风暴单体，其中20:52预报该单体将产生下击暴流；20:35-21:27连续10个体扫识别到标号为B5的风暴单体，其中20:40-21:21预报该单体将产生下击暴流；21:21-21:38连续4个体扫识别到标号为C1的风暴单体，其中21:21和21:27预报该单体将产生下击暴流。图 2给出了20:35-21:38时间段的下击暴流临近预报结果。在21:21，算法首次识别到单体C1，此时预报B5和C1都将产生下击暴流，到21:27算法没有再预报B5将产生下击暴流(下击暴流概率为41%)，只预报C1将产生下击暴流，之后单体B5消散，2个体扫后单体C1消散。在21:21，单体B5和C1距离很近，它们的极坐标质心位置分别为(337.3°、56.6 km)和(346.5°、59.9 km)，单体B5正好位于弓形回波的顶点处，离客轮倾覆地点更近，且移动方向正朝向客轮倾覆地点。Wilson et al(1984)认为典型的导致微下击暴流的下沉气流强度达到峰值后的0~10 min内会影响到地面，考虑到一个体扫周期为6 min，我们认为是单体B5产生的微下击暴流导致了“东方之星”客轮的侧翻。图 3对单体B5的雷达特征量随时间的变化进行了分析。

从图 3a和3b可看出，该单体的最大反射率因子(ZMAX)、最小相当位温高度附近的最大反射率因子(ZTH)、顶高(HT)、底高(HB)、质心高度(HC)、垂直积分液态水含量(VIL)和强冰雹指数(SHI)都有一个先增加达到最大值后逐渐减小的趋势，且在最后一个体扫，HT和VIL快速下降。从图 3c可看出，MARC特征的垂直积分径向辐合量IMARC和最小相当位温高度附近的径向辐合量CMTH都有一个先增加后减小的趋势，但是它们达到峰值的时间要落后于基于反射率因子得到的一些雷达特征量。从图 3d可看出，在单体B5生命周期的大部分时间里，中气旋的旋转速度都只维持弱中气旋，只有最后3个体扫才达到中等强度中气旋；1~6 km高度径向上的最大、最小速度差几乎都大于20 m·s^-1，最大径向辐合等级RCL都≥40×10^-4 s^-1，表明存在显著的中层径向辐合特征。从图 3e可看出，只有在单体B5生命周期后半段的5个体扫时间，风暴底位于1 km以下；在21:21，风暴底的最大绝对径向速度MRV01达到最大(22.5 m·s^-1)，径向辐散值MDIV也达到最大(104×10^-4 s^-1)，且出现一个陡升现象。从图 3f可看出，从20:41-21:21，单体B5产生强冰雹的概率都不大于16%，产生下击暴流的概率都大于60%，最早预报单体B5将产生下击暴流的时间比“东方之星”客轮倾覆时间提前了47 min。

图 4给出了21:21-21:32时段3个体扫时次的0.5°仰角的径向速度，从图中可看出，在21:21，弓形回波顶点后侧有一块负速度大风核区，但是不存在辐散式的正负速度对；到21:27，弓形回波顶点后侧出现3块相邻的正负速度对，不过正速度大小和面积都不大(见图中虚线椭圆区)；到了21:32，弓形回波顶点后侧的3个辐散式的正负速度对非常清晰，表明发生了相邻的3个下击暴流。

4 预报效果评估

使用湖北省2019年6-8月的17个雷暴大风过程(表 2)对下击暴流预报击中率进行了评估，评估范围为距离雷达10~150 km。由于本算法是在识别风暴单体的基础上进行下击暴流临近预报，在挑选个例时，未考虑冷空气叠加雷暴大风的个例；评估时，尽量不考虑阵风锋(对应弱回波或无回波)引起的大风。评估方法为：在区域自动气象站出现大风之前1小时内和站点周围40 km范围内(以统计风暴平均移速40 km·h^-1为依据)，如果有风暴单体被识别且其最大反射率因子≥35 dBz，那么当下击暴流概率≥50%时，则判断该站点预报下击暴流正确，否则为漏报；如果没有风暴单体被识别或者有风暴单体被识别但是其最大反射率因子 < 35 dBz时，则认为该站点大风为非下击暴流大风，不参加评估。这种评估方法不能完全排除脱离风暴母体40 km内的出流大风的影响，但是考虑到这种出流也是从母体脱离的，如果预报了下击暴流也算正确。

以各雷达为中心，统计10~150 km范围内的区域地自动气象站观测到的雷暴大风共有353站次(当相邻雷达间距小于300 km，且选取的个例落在其观测重叠区域时，会重复统计)，其中识别出风暴单体并正确预报有下击暴流大风的有242站次，平均预报时效为39 min；识别出风暴单体但没有预报有下击暴流大风的有38站次(下击暴流发生概率位于35%~49%)，其中回波面积很小的有6站次，在下击暴流发生后预报有下击暴流的有9站次；未识别出风暴单体由阵风锋(有窄带弱回波)产生的大风有40站次；未识别出风暴单体且几乎无回波的大风有33站次(距离雷达较远的阵风锋雷达探测不到)。去掉弱回波、无回波的情况，共有280站次的雷暴大风可认为是下击暴流大风，正确预报242站次，击中率为86.4%。按照回波形态分类评估，飑线类、线状对流类(飑线除外)和非线状对流类风暴的下击暴流预报击中率分别为93.2%、90.5%和75.6%。对于非线状对流类风暴，MARC或后侧入流急流特征没有飑线类和线状对流类的那么明显，这可能是其预报正确率较低的原因之一。

5 结论与讨论

(1) 综合使用雷达反射率因子、径向速度和基于探空观测的0℃、-20℃和最小相当位温高度等资料，基于风暴单体识别追踪、冰雹指数计算、MARC特征识别和中气旋识别等算法，提取风暴单体的诸多雷达特征量，经统计分析挑选出下击暴流的雷达先兆因子9个作为模糊逻辑法的输入，建立了下击暴流临近预报方程。

(2) 使用2015年6月1日发生在监利导致“东方之星”客轮侧翻的下击暴流个例对本算法进行了测试，结果表明本算法从20:41-21:21共有8个体扫时次预报了引起沉船事件的那个风暴单体将会产生下击暴流，首次预报时间比客轮侧翻时间早47 min。21:27-21:32，0.5°仰角的径向速度中出现了辐散式的正负速度对，表明下击暴流大风的发生。

(3) 使用2019年6-8月发生在湖北省的雷暴大风个例对本文提出的下击暴流临近预报算法进行了效果评估，结果表明本算法对下击暴流临近预报的击中率为86.4%，平均预报时效为39 min；按回波形态分类评估，飑线类、线状对流类(飑线除外)和非线状对流类风暴的下击暴流临近预报击中率分别为93.2%、90.5%和75.6%。

在本算法中没有使用识别的风暴单体区域内的径向速度直接提取与MARC特征和中气旋有关的雷达特征量，而是单独识别了MARC特征和中气旋后，再把它们与风暴单体关联。这样做的原因是考虑到风暴单体的识别是尽可能识别强风暴核，其对应的径向速度场中可能不包含完整的MARC特征或中气旋结构，MARC特征可能部分位于风暴后侧中层的入流槽口中，中气旋核区可能部分位于弱回波或有界弱回波区内。本文没有评估虚警率，今后将继续研发1 km高度以下的辐散式速度对和大风核的自动识别算法，在此基础上，结合地面自动站观测进行下击暴流临近预报虚警率的批量评估。该算法模块已集成到中国气象局武汉暴雨研究所研发的分类强对流天气自动识别预警系统中(崔春光等，2021)，并于2019年开始投入业务运行，今后将在业务应用中不断地优化算法，特别是需要使用更多的产生下击暴流的孤立单体风暴来统计各雷达特征量，修改下击暴流临近预报方程，因为对于线状对流和团状对流来说，很难把地面发生的大风与风暴单体进行准确匹配。对于已出现的下击暴流大风，将加入1 km以下的大风核区的识别和0.5°仰角的径向速度出现的辐散式的正负速度对或强辐散区的识别。对于提取的诸多雷达特征量，今后也将使用机器学习方法进行下击暴流的临近预报技术研究。