2013, 24 (1): 117-125
2. 成都军区空军气象中心，成都 610041
2. Air Force Meteorological Center of Chengdu Military Region, Chengdu 610041
一般强风暴或飑线中处于成熟阶段的单体下沉气流在近地面处沿水平方向推进，常常形成辐散性的阵风，当阵风前沿的辐合线达到一定强度时，就可称为阵风锋；由于在天气雷达的强度回波图上阵风锋表现为窄带回波，因此在雷达气象学中，也将阵风锋称为窄带回波。风暴中的强烈出流不仅温度低而且具有较高的水平动量，故阵风锋常引起气压突变、风向突变、风速突变、温度急降和相对湿度上升等强烈天气现象。一般情况下阵风锋不伴随有降水发生，但常常造成严重的地面风灾和强烈的风切变，因此，对阵风锋的研究已成为短时临近预报关注的重点之一。
目前，国内外对阵风锋的研究较多侧重于阵风锋的形成机制探讨和天气过程分析，Simpson[1]指出阵风锋的动力机制和密度流类似，造成阵风锋的那股冷空气是由降水拖曳和雨滴冷却蒸发综合作用生成的。Wakimoto[2]，Wilson等[3]、葛润生[4]、陈明轩等[5-6]和王彦等[7]均利用实测资料分析了两种类型 (运动型和静止型) 的阵风锋，结果表明：运动型的阵风锋常与猛烈发展的强风暴相伴随，它的出现预示着风暴将持续、猛烈地发展；静止型的阵风锋对应风暴发展后期，它的出现加速了风暴的消亡；两类不同的阵风锋反映了下沉气流与环境的不同作用，沿阵风锋抬升的上升气流可以使母体风暴得以维持，有可能形成新的对流风暴。葛润生[4]、黄璇璇等[8]先后提出雷达强度回波图上细长的窄带回波是由于冷暖气团交界面两侧的折射指数的突变, 而导致对电磁波的散射或反射形成的阵风锋回波。毕旭等[9]、刘勇等[10]研究表明:窄带回波的强度代表其中降水粒子的贡献，阵风锋的真实强度和持续时间取决于母体风暴的强度和后部冷空气堆中积蓄的能量。何彩芬等[11]和姚建群等[12]指出阵风锋在低空引起强烈的风切变，会造成飞行器在起飞或降落时发生事故。
国外较早开展阵风锋识别算法研究，美国联邦航空局 (FAA) 机场多普勒天气雷达系统 (TDWR) 通过对径向速度辐合线识别，并结合反射率因子和方位角切变形成了阵风锋探测算法[13]，随后Delanoy等[14]和Troxel等[15]提出了函数模版相关法识别窄带回波，该方法应用于ASR-9和TDWR雷达中，随后该方法又应用于NEXRAD中，宗蓉[16]分析了阵风锋的速度场特征，通过搜索速度辐合区来识别阵风锋，取得了一定的效果。李劲[17]借鉴了文献[14-15]提出的方法设计了窄带模版，并结合反射率利用函数评分实现对窄带回波的识别。陈刚[18]提出了在频域找出窄带回波区域区别于其他区域的形状特征的算法, 观察窄带回波与其他回波的功率谱，发现窄带回波区域的功率密度相对其他闭合区域更集中于峰值附近。王楠等[19]利用径向切变和切向切变计算出组合切变以识别切变线和辐合线。
目前，阵风锋识别算法尚未成熟，识别效果不甚理想，国内布网的多普勒天气雷达用户终端尚未开发出阵风锋的识别产品，本文基于对雷达原始基数据进行一定质量控制，充分利用强度场窄带回波的结构特征、几何特征，设计了双向梯度算法来识别窄带回波；在速度场中通过径向辐合来识别风的突变；最后综合双向梯度和切变辐合得到结果，确定出阵风锋的准确位置。文中还利用河南商丘、郑州及安徽阜阳3个雷达站探测的阵风锋98个体扫资料，结合地面自动气象站要素资料对本文的识别算法进行效果检验，并利用临界成功指数进行了效果评估。
1 数据的预处理多普勒天气雷达能够获得较高时间和空间分辨率的体扫数据。但在阵风锋识别前，必须先对原始数据进行一定的质量控制，避免后续工作引起较大的估算误差和错误，也有利于提高算法的准确性。
本文对强度数据的质量控制包括强度数据的填补、噪声滤除和数据平滑。数据填补是针对单个缺值点或者单个缺测径向；噪声滤除主要针对孤立点和奇异点的滤除；数据平滑采用3×3网格的中值滤波。对速度数据的质量控制包括速度数据的填补 (处理同强度数据)、数据平滑和速度退模糊，数据平滑为沿径向进行九点滑动平均；速度退模糊的方法类似于WSR_88D雷达的处理算法，即根据连续性原则, 将每个速度初猜值与其周围的相邻速度值相比较，先计算平均值，用相同径向上 (靠近雷达) 的4个邻近距离库的速度值和紧邻的已用算法处理过的前一个径向上的5个距离库的速度值计算平均速度，然后将正在进行退模糊处理的距离库的第1猜速度与这个平均值作比较。如果其差值大于最大可测速度，那么该点的可替代值为第1猜测速度减去 (或加上)2倍的最大可测速度，如果该可替代值大于最大可测速度，去除该点，否则，保留该点。
经过质量控制后，强度场图像中的窄带特征更加显著，轮廓更加分明 (噪声被滤除、缺测点和缺测径向数据被补齐、窄带回波周围的毛糙点被滤除)；速度场图像中阵风锋对应的辐合 (切变) 区域更加清晰规整 (噪声被滤除、缺测点和缺测径向被补齐、径向数据被平滑、速度模糊区域被纠正)。因此，可以说明一定的质量控制对实现阵风锋识别是有利的，可为后续算法准确识别打下基础。
2 辐合线 (切变线) 识别阵风锋在速度图上表现为径向速度的辐合 (切变), 可以分为风向切变和风速辐合两种情况，这两种情况的径向速度都在锋线附近有明显变化。第1种是径向速度的风向表现为由正速度向负速度的转变，第2种是径向速度的速度值表现为由较高的值向较低的值转变。在锋线两侧，两种情况对应的速度值均呈一致减小，因此，根据该特征，首先沿着径向距离递增的方向寻找径向速度值连续减小的一组相邻距离库，可得到一维的辐合段。但考虑到雷达接收到的速度数据有很大的脉动性，在经过径向上九点滑动平均之后，相邻两个距离点上的速度大小仍很难代表某一段距离上速度的增加和减少的趋势；因此，对径向上的每一个点 (指有速度资料的点) 都计算一个辐合 (切变) 参量S[13]，
|
(1) |
式 (1) 中，Vi表示在与计算点同一径向上，计算点前后各3个点的任一速度值，Δr为径向速度的距离分辨率。通过计算每个计算点上S的大小、正负来判断在此点附近速度是增大还是减小的趋势。在具体算法中为S设定一个阈值，当某一计算点上的S值小于此阈值，就将其保留，归于一维辐合段中的一点。图 1给出了2009年6月3日14：15(世界时，下同) 商丘雷达的一条阵风锋锋线径向数据的处理结果。由图 1可以看出，原始径向速度经过质量控制后大部分脉动被消除，曲线被平滑；再次经过S＜-2选定后的曲线的l1段和l3段因下降程度不满足阈值而被排除，而l2段的下降程度满足阈值被保留，即为选定的辐合段。
|
|
| 图 1. 2009年6月3日14:15商丘雷达速度场锋线上的辐合段选取(仰角为0.5°, 方位角为160°) Fig 1. Convergence segment search across the front line of Shangqiu radar at 1415 UTC 3 June 2009(elevation:0.5°; azimuth:160°) | |
为了描述每个辐合段的强度，计算了如下物理量：
|
(2) |
|
(3) |
|
(4) |
|
(5) |
式 (2)～(5) 中，vb表示辐合段起点处的径向速度；ve表示辐合段终点处的径向速度；rb表示辐合段起点离雷达站的距离；re则表示辐合段终点离雷达站的距离。因阵风锋必须在速度梯度和动量上满足一定条件，故引入了4个阈值：速度梯度下限 (Lg)，动量下限 (Lf)，速度梯度上限 (Hg)，动量上限 (Hf)。如果某一辐合段的速度梯度g或者动量f低于其下限值, 那么，该辐合段将被删除，否则，再对其进行高值检验；如果该辐合段的g或f中任一个值大于上限, 则此辐合段将被保存。在所有保存的辐合段中，判断当前辐合段与下一辐合段之间的径向中心和所在方位之差，如果差值大于径向距离阈值或方位阈值，即|si-si+1|＞ΔS，|azi-azi+1|＞Δaz (Δs为径向距离阈值，Δaz为方位阈值) 时，该辐合段被删除，否则保留。将所有满足条件的辐合段合并成二维特征，并滤除辐合段数量过少的特征，形成最终识别的辐合线 (切变线)。
3 窄带回波识别阵风锋在强度场上表现为细长的带状，通常称之为窄带回波。天气雷达探测到的阵风锋常常出现在低仰角层，回波高度一般不超过3 km，强度一般为10～35 dBZ，更高的强度通常是降水粒子的贡献，不代表阵风锋回波的真实强度。
本文先对所有回波强度值进行分级处理，即将强度值按递增划分成若干个区间，每个区间用该区间最高值代替。如表 1所示，其中-999代表无效值。因雷达图像在绘制时，是将不同区间的反射率因子值用特定的色标表示，分类后的强度场的回波值虽然被改变，但所属区间并未改变，即颜色不发生变化，所以分级处理前后的图像完全相同，图像保持了原有特征，因此可用分级后的强度场代替原始强度场作为窄带回波的识别算法的处理对象。
|
表 1 反射率因子分级表 Table 1 The reflectivity classification table |
窄带回波一般出现在强风暴回波移动方向的前沿，与强风暴回波之间存在一定间距，不同阶段间距不等。在初始阶段，窄带离母体风暴较近，随着母体风暴的发展，窄带与母体逐步分离，它们之间的距离越来越大，且窄带回波的强度较其前后方一般要高5 dBZ或以上。利用这个特点本文通过计算每个回波点的双向梯度来保留线状的回波，如图 2所示，点O代表当前所要计算的回波点，点A，B为与O同根径向相距5个距离库的点 (考虑一般的窄带回波宽度＜10 km)，点C，D为与O同距离圈相距5个方位的点，R(O)，R(A)，R(B)，R(C)，R(D) 为对应点处的反射率因子值。计算了如下两个方向的梯度，横向梯度：
|
(6) |
|
|
| 图 2. 双向梯度示意图 Fig 2. Bilateral grads schematic diagram | |
纵向梯度：
|
(7) |
对两个梯度进行判断：当GOC和GOD均为0，且GOA和GOB均大于等于5 dBZ时保留该点；如果不满足，则删除该点。
实际工作中无法得知窄带回波与径向线的夹角，因此在具体算法中设计了12个不同夹角的模型，即每隔15°建立1个模型。经过双向梯度处理后的强度场大部分降水回波被滤除，窄带回波被完整保留，但图像中仍然存在若干短线。为了去除这些短线，先将图像细化，即只保留每个强度段的中心点 (强度段是指从雷达站沿径向穿越某一回波内部的段称为强度段)，并记录下段宽，随后利用递归算法计算每根短线的长度，将不满足长度阈值的短线滤除。在此过程中保存每个有效点的方位角和径向库数。最后只保留一定长度的短线 (阈值可调)，并利用先前记录的段宽还原成带状，形成最终的窄带回波识别图像。
4 锋线定位与预警理论上，同一条阵风锋的窄带回波和辐合线 (切变线) 在空间上是完全对应的，但从实际得到的阵风锋回波图中会发现，窄带回波和辐合线 (切变线) 往往只有大部分重合，甚至会出现基本不重合的情况。因此，为了较为完整地定位出锋线，本文将识别的窄带回波和辐合线 (切变线) 综合，即某点对应位置处在窄带回波图像或辐合线图像中有一方存在有效值，则保留该点，作为锋线中的一点。再次，将综合得到的图像细化，并去除不满足长度阈值的短线 (方法如强度场)，保留下来的线为最终识别的锋线。在实际业务中，预报员更关心阵风锋在速度场引起的风切变，因此为了定量描述阵风锋的强度，将如下物理量作为输出：锋区最大正径向速度VP(单位：m·s-1)；锋区最大负径向速度VN(单位：m·s-1)；锋区两侧平均切变强度 (单位：m·(s·km)-1)：
|
(8) |
锋线上最大径向速度 (单位：m·(s·km)-1)：
|
(9) |
锋线长度：
|
(10) |
其中，Vn为辐合段中任意点的径向速度，Rn为该点离雷达站的距离；m为辐合段的数目。计算得到的5个物理量均可在相应的软件终端显示，以定量化所识别的阵风锋。为了更容易引起预报人员的注意，程序为识别的图像设定了2 s的定时器，即每隔2 s闪烁1次；预报员可以根据锋线的闪烁报警和界面提示的物理量，结合多个时次识别的锋线位置对下游区域发布大风预警。值得注意的是, 不同地域的阵风锋，算法采用的相关阈值不一样，因此在业务应用中还需对阈值做本地化处理，调整阈值以达到最优的识别效果。
5 识别效果利用2009年6月3日河南商丘、郑州和2009年6月5日安徽阜阳3个雷达站探测到的阵风锋雷达资料对本文的识别算法进行效果检验，并利用临界成功指数对98个体扫样本的识别效果进行评估。
5.1 地面要素变化实况2009年6月3日阵风锋于10：00左右开始影响郑州站，10：00—13：00测站记录表明，气温逐渐下降7.5℃，气压上升4 hPa，2 min平均风速从4.6 m·s-1增强至7.2 m·s-1，随后又逐渐减小。14：00左右阵风锋到达商丘站，14：00—15：00商丘站气温骤降9.4℃，气压上升3.6 hPa，2 min平均风速从1.7 m·s-1增强为15.6 m·s-1。2009年6月5日，阵风锋于12：00左右开始影响阜阳站，12：00—13：00阜阳站气温骤降9.7℃，气压上升3.7 hPa，2 min平均风速从3.1 m·s-1增强到6.5 m·s-1。以上地面要素变化说明，阵风锋过境时，3个站的气象要素均发生了变化，即地面气温骤降、气压上升、风速突增。
5.2 识别效果分析2009年6月3日从郑州雷达上可初次判断出阵风锋的时刻为11:07，11:13在雷达站东南测开始可以看到弱小的窄带回波 (如图 3所示)，对应速度场中有辐合线，但位置较窄带回波偏西，算法成功识别出了阵风锋；11:19—11:44窄带回波较窄，数据密度不够，辐合线消失，算法未能识别出阵风锋；11:50—13:21窄带回波加强，辐合线重新出现并进一步加强，算法均成功识别出阵风锋；13:27—13:33阵风锋处于消亡阶段，窄带回波和辐合线出现断裂并逐渐消失，算法未能识别出消亡时刻的阵风锋。
|
|
| 图 3. 2009年6月3日郑州雷达阵风锋识别图 Fig 3. Identification image of gust front of Zhengzhou radar on 3 June 2009 | |
商丘雷达在12:32—15:46也观测到了阵风锋 (如图 4所示)，12:32窄带回波较弱且与前端回波混杂，表现得并不明显，速度场中出现了辐合线，但由于部分辐合线与径向线平行，因此所识别的锋线较短。13:08窄带回波的强度值不统一，大部分为20 dBZ (归类后)，但内部为25 dBZ，影响了窄带回波识别。此外，速度场中辐合线基本与径向平行，辐合线也未能成功识别。13:14—15:22窄带回波和辐合线加强，均表现得较为明显，算法均成功识别出阵风锋。算法识别的过程中出现了偏短的情况，主要有两个原因：首先由于双向梯度法与窄带回波到母体回波之间的间距大小密切相关，当窄带回波有一段与母体回波相连时，窄带回波识别将受到一定影响，但一般该相连段也较短或仅限于阵风锋初生阶段；其次由于本文辐合线的识别子算法基于径向，当辐合线与径向线近似平行时，辐合线将无法识别，但该情况较少出现。
|
|
| 图 4. 2009年6月3日商丘雷达阵风锋识别图 Fig 4. Identification image of gust front of Shangqiu radar on 3 June 2009 | |
图 5给出了2009年6月5日安徽阜阳雷达探测到的阵风锋3个时次识别图，该过程阵风锋初始时刻为11:05，11:05—12:25窄带回波与降水回波严重混杂，辐合线偏弱或与径向近似平行，算法未能识别出阵风锋；12:31—14:39窄带回波和辐合线加强，算法均成功识别出阵风锋；14:45—14:51阵风锋处于消亡阶段，窄带回波和辐合线断裂并逐渐消亡，算法失效。
|
|
| 图 5. 2009年6月5日阜阳雷达阵风锋识别图 Fig 5. Identification image of gust front of Fuyang radar on 3 June 2009 | |
5.3 识别效果评估
本文选取这3个雷达站探测到的98个阵风锋体扫描样本 (从锋线初生到消亡)，其中，商丘33个，郑州25个，阜阳40个，算法识别情况如表 2所示。商丘阵风锋过程成功识别的样本数为26个，未能识别的样本数为7个；郑州阵风锋过程成功识别的样本数为17个，未能识别的样本数为8个；阜阳阵风锋过程成功识别的样本数为24个，未能识别的样本数为16个；3个站探测的阵风锋资料误识别样本数均为0。
|
|
表 2 阵风锋过程样本数和识别情况 Table 2 The number of gust front processes, samples and identification |
利用临界成功指数ICS、命中率RH、漏报率RM和虚警率RFA对总样本的识别效果进行评估，
|
(11) |
其中，x为有阵风锋且成功识别的样本数，y为有阵风锋但未能识别的样本数，z为无阵风锋但误识别为阵风锋的样本数。算法对98个样本的临界成功指数达到了0.684，成功识别率为68.4%，表明本方法对阵风锋识别效果较好 (表 3)。
|
|
表 3 总样本的临界成功指数ICS、命中率RH、漏报率RM和虚警率RFA Table 3 ICS, RH, RM and RFA of total samples |
5 结论
本文在对天气雷达体扫基数据进行一定质量控制的基础上，根据阵风锋在强度场和速度场所表现的回波特征和几何特征，通过识别窄带回波和辐合线 (切变线) 实现阵风锋的识别，最后利用3个雷达站探测的阵风锋资料，对识别算法进行效果检验，并利用临界成功指数进行评估。得到如下结论：
1) 根据窄带回波与母体回波的间隔，且窄带回波强度较前后高5 dBZ或以上，提出了双向梯度法，该方法可以有效滤除大范围降水回波而保留窄带回波。
2) 结合组合反射率因子产品，对比各个仰角层的回波发现，高仰角层对应位置处均不存在窄带回波或不存在强于窄带回波的反射率因子。因此, 算法只需考虑较低仰角层，大大提高识别效率。
3) 针对速度场采用的方法能有效识别出径向辐合线，同时也适用于低空径向风切变和辐合线的识别，但当辐合线与径向线平行或近似平行时，识别将被干扰；强度场中将图像沿径向细化可以消除类似带状回波的干扰。
4) 利用临界成功指数对98个体扫样本进行识别效果评估可知，本识别算法识别率达到68.4%，表明利用本算法对阵风锋自动识别有一定效果。
| [1] | Simpson J E A. Comparison between laboratory and atmospheric density currents. Quart J Roy Meteor Soc, 1969, 95: 578–765. |
| [2] | Wakimoto R M. The life cycle of thunderstorm gust front as viewed with Doppler Radar and Rawinsonde Data. Mon Wea Rev, 1982, 110: 1060–1082. DOI:10.1175/1520-0493(1982)110<1060:TLCOTG>2.0.CO;2 |
| [3] | Wilson J W, Schreiber W E. Initiation of convective storms at radar observer boundary layer convergence lines. Mon Wea Rev, 1986, 114: 2516–2536. DOI:10.1175/1520-0493(1986)114<2516:IOCSAR>2.0.CO;2 |
| [4] | 葛润生. 阵风锋的雷达探测和研究. 气象科学研究院院刊, 1986, 1, (2): 113–121. |
| [5] | 陈明轩, 愈小鼎, 谭晓光, 等. 对流天气临近预报系统技术的发展与研究进展. 应用气象学报, 2004, 15, (6): 754–766. |
| [6] | 陈明轩, 高峰, 孔荣, 等. 自动临近预报系统及其在北京奥运期间的应用. 应用气象学报, 2010, 21, (4): 394–404. |
| [7] | 王彦, 于莉莉, 李艳伟, 等. 边界层辐合线对强对流系统形成和发展的作用. 应用气象学报, 2011, 22, (6): 724–731. DOI:10.11898/1001-7313.20110610 |
| [8] | 黄旋旋, 何彩芬, 徐迪峰, 等. 阵风锋过程形成机制探讨. 气象, 2008, 34, (7): 380–387. |
| [9] | 毕旭, 刘慧敏, 赵榆飞. 陕北系列阵风锋天气过程分析. 陕西气象, 2008, (2): 23–26. |
| [10] | 刘勇, 王楠, 刘黎平. 陕西两次阵风锋的多普勒雷达和自动气象站资料分析. 高原气象, 2007, 26, (2): 380–387. |
| [11] | 何彩芬, 姚秀萍, 胡春蕾, 等. 一次台风前部龙卷的多普勒天气雷达分析. 应用气象学报, 2006, 17, (3): 370–375. |
| [12] | 姚建群, 戴建华, 姚祖庆. 一次强飑线的成因及维持和加强机制分析. 应用气象学报, 2005, 16, (6): 746–753. DOI:10.11898/1001-7313.20050615 |
| [13] | Uyeda H, D Zrnic S. Automatic detection of gust front. J Atmos Oceanic Technol, 1985, 3: 36–50. |
| [14] | Delanoy R L, Troxel S W. The Machine Intelligent Gust Front Algorithm. MIT Lincoln Laboratory Project Report ATC-196, 1993. |
| [15] | Troxel S W, DelanoyR L, Pmorgan J P, et al. Machine Intelligent Gust Front Algorithm for the Terminal Doppler Weather Radar (TDWR) and Integrated Terminal Weather System (ITWS). AMS Workshop on Wind Shear and Wind Shear Alert Systems, 1996: 70–79. |
| [16] | 宗蓉. 多普勒天气雷达的阵风锋识别方法探索. 南京: 南京信息工程大学, 2009. |
| [17] | 李劲. 利用多普勒天气雷达自动识别阵风锋方法研究. 南京: 南京信息工程大学, 2010. |
| [18] | 陈刚. 阵风锋的检测与识别. 西安: 西安电子科技大学, 2009. |
| [19] | 王楠, 刘黎平, 徐宝祥, 等. 利用多普勒雷达资料识别低空风切变和辐合线方法研究. 应用气象学报, 2007, 18, (3): 314–320. |