基于模糊C均值的低空风切变预警算法

熊兴隆，杨立香，马愈昭，庄子波

(1.中国民航大学 天津市智能信号与图像处理重点实验室，天津 300300； 2.中国民航大学 民航气象研究所，天津 300300)

0 引言

低空风切变[1]是指600 m以下风速和风向在空间突然发生的变化，是一种严重影响飞机安全起降的天气现象。风切变[2]形成机制复杂和影响强度大，准确地判断风切变的位置和强度已成为航空领域的一个重要研究课题。与传统天气雷达相比，新一代多普勒天气雷达[3]不仅能监测大范围降水天气过程，还可以提供中小尺度的风场信息，已成为民航机场风切变探测的重要设备。在低空风切变检测方面，国内外学者进行了广泛的研究。1986年Uyeda等[4]利用径向和切向的组合梯度值来探测风切变，该方法适合有明显辐合线的强风切变识别。1989年Evans等[5]利用终端多普勒天气雷达(Terminal Doppler Weather Radar, TDWR)自动识别切变，该方法需采集大量不同地区的风切变数据集。1994年Wolfson等[6]运用智能图像处理以及数据融合技术识别风切变，该方法较适合微下击暴流引起的切变。2010年魏耀等[7]运用最小二乘法拟合径向和方位切变，改善了边缘识别效果，但该方法容易漏掉其他方向的切变。2013年Miller等[8]采用均值去噪，利用雷达的方位切变实现了气旋跟踪，但该方法对于其他对流天气引起的径向切变检测受限。2014年Zheng等[9]结合多普勒雷达资料提出了一种自动识别阵风锋的算法，该算法适合速度梯度值明显的切变识别。2015年Tagliaferri等[10]运用神经网络和支持向量机预报短时风向，该方法需要大量数据集。2016年赵丽艳[11]以多尺度进行最小二乘法径向数据拟合，该算法较适合单条径向的风切变识别。2017年Hwang等[12]使用模糊神经网络模型取得了较高的辐合线识别率，该方法需花费较长时间训练样本。目前大多数风切变识别算法在预处理时采用均值去噪，风切变搜索时是基于径向或切向的梯度值差异，或者是两者的结合，这适合强切变识别，但容易漏掉小切变。

针对小切变的检测问题，本文提出一种基于模糊C均值(Fuzzy C-Means, FCM)的低空风切变预警算法。速度基数据预处理时采用TV模型去噪，保持速度数据的细节信息，检测切变时采用8邻域系统获取4个方位梯度值，运用FCM思想将梯度值分为高低切变，采用武汉暴雨研究所提供的实测基数据进行验证。

1 雷达速度基数据的预处理

由于新一代多普勒天气雷达探测到的风场信息存在数据缺测、噪声问题，需要对雷达径向速度数据进行预处理。对于缺少的速度基数据，采用较为实用快速的线性内插法填补，然后对处理后的数据进行噪声的剔除。去噪方法一般为均值滤波方法(Mean Filtering Method, MFM)，其处理速度较快，但数据细节和边缘信息容易丢失，容易漏掉小切变。由于TV模型去噪[13]能保留边缘信息，本文将此方法首次应用在多普勒雷达速度基数据去噪中，同时保留速度数据的小切变细节。

TV模型由Rudin等[14]在1992年提出，用有界变差函数空间(Space of Functions for Bounded Variation, BV空间)对图像进行建模，并利用全变差方程实现。TV模型首先建立能量泛函及其约束条件；通过变分原理，求得欧拉-拉格朗日(Eulef-Lagrance, E-L)方程；在初值条件下求近似解，即求偏微分方程的数值解。令Ω为Rn中的有界开子集，u∈L1(Ω;R)为局部可积分函数，则TV模型可定义为：

TV[u(x,y)]=∬Ω|▽u(x,y)|dxdy

(1)

其中：u(x,y)为被处理的速度基数据，▽u(x,y)为速度基数据的梯度值。

建立能量泛函函数，其函数的约束优化表达式为：

(2)

通过E-L方程，式(2)可表示为：

(3)

由于本文处理的速度数据为离散的数值，将E-L方程离散化，其表达式为：

(4)

2 低空风切变预警算法

新一代多普勒天气雷达速度基数据进行预处理之后，要实现低空风切变预警功能，必须获得速度的梯度值。首先采用每个速度基数据及其8邻域系统分别对应的速度值获取4个方位梯度值，提取满足阈值的风切变值;然后运用FCM思想将风切变值分为高切变和低切变两类，从而实现高低切变预警。

2.1 梯度值获取

新一代多普勒天气雷达通过在一系列固定仰角上扫描360°的方式获取基数据。在每个仰角上，以雷达为中心波束向外发射电磁波，当径向距离增加时，距离地面的高度也增加。新一代多普勒天气雷达径向速度是用脉冲对技术进行估计，其中脉冲与距离门之间复杂相关，然后以笛卡尔坐标[15]对数据进行处理。基于径向或者切向的风切变识别算法容易漏掉其他方向的切变信息，本文采用8邻域系统计算4个方位梯度值，8邻域系统如图1所示。

图1 中心点p及其8邻域系统 Fig. 1 Center p and its 8-neighborhood system

为去除地杂波[16]，一般当多普勒速度绝对值小于0.1 m/s时被视为无效数据。使每个速度基数据点p及其8邻域系统分别对应的速度值与4个方向模板依次进行卷积，获得4个方位差值，分别为V(N)-V(S)、V(NE)-V(SW)、V(NE)-V(SE)和V(E)-V(W)，其中方向模板如图2所示。

图2 4个方向模板 Fig. 2 Template of four directions

分别计算每个点对应的4个方位差值的绝对值，求取速度梯度值▽v，其中▽v定义如下:

▽v=mean[|V(N)-V(S)|/dN-S,|V(NE)-V(SW)|/dNE-SW,|V(NE)-V(SE)|/dNE-SE,|V(E)-V(W)|/dE-W]

(5)

其中：mean[a,b]=(a+b)/2，V(N)为当前点正北方向的近邻速度值，dN-S为正北方向与正南方向分别对应的近邻点之间的距离，其他三个方位的对应含义与V(N)及dN-S类似。

2.2 基于FCM的4个方位搜索风切变

由于不同对流引起的切变强度不同，本文运用FCM算法，根据风速梯度识别低切变和高切变，判别不同位置的切变相对强弱程度。FCM算法[17-18]是一种基于距离的聚类算法，通过模糊理论对数据进行分析和建模，建立数据类属的不确定性描述，具有较高的搜索速度。本文基于FCM的4个方位搜索风切变的流程如图3所示，描述如下。

图3 基于FCM的4个方位搜索风切变流程 Fig. 3 Flow chart of wind shear identification with FCM method

1)读取雷达速度基数据，采用较为实用的线性内插法添加缺测数据，利用TV模型去噪，同时保持细节信息。

2)运用4个方向模板与预处理数据中的每个点及其8邻域系统对应的速度值依次进行卷积，获得4个方位差值，分别计算绝对值，获得速度梯度值。

3)根据国际民航组织建议，水平风切变标准采用美国机场低空风切变报警系统[19]的报警标准。美国的报警标准规定水平风切变2.6 (m/s)km可作为对飞行构成危害的强度标准。由于新一代多普勒天气雷达距离分辨率为250 m，本文水平风切变的阈值为10.4 (m/s)km，根据阈值提取有效风切变值。

4)运用FCM思想对切变值分类，算法具体描述如下。令X={xi,i=1,2,…,n}为训练集，X⊆Rn，c为预定的类别数目，可取任意整数，vi(i=1,2,…,c)为第i个聚类的中心，uik为第k个样本对第i类的隶属度函数。本文设置c=2，最大迭代次数l=100，FCM算法的隶属矩阵可初始化为：

U={uik}； ∀k=1,2,…,n

(6)

通过式(5)，计算类别中心，计算式为：

(7)

其中:m为一个加权数，一般取2。

计算目标函数Jm(U,v)，计算式为：

(8)

(9)

其中∀i=1,2,…,c；∀k=1,2,…,n。

5)根据风速梯度，输出高低切变，进而实现不同强度的风切变预警功能。

3 算法对比及切变强度案例分析

武汉暴雨研究所为中国气象局气象研究院所之一，主要任务为瞄准暴雨防灾减灾的国家目标和气象业务服务需求，开展中国暴雨的应用基础研究、应用研究和技术开发。本文采用武汉暴雨研究所提供的实测速度基数据进行实验，运用Microsoft Visual Studio 2012和Matlab 2013进行算法编程。风切变预警算法在预处理时一般采用MFM，在风切变搜索时一般采用径向或切向。本文算法在预处理时采用TV模型，在风切变搜索时采用4个方向模板。通过MFM算法与本文算法对比，分析两种算法的速度数据细节保留以及风切变识别效果。通过依次识别龙卷风和阵风锋引起的高切变和低切变，说明本文算法预警高低切变的可行性。

3.1 算法对比

在速度基数据预处理时，分别采用MFM与TV模型去噪，对比分析去噪效果和速度基数据细节信息的保留情况。在检测切变时，对于径向速度图，依次采用MFM预处理最小二乘法合成切变和本文算法检测，对于平面位置显示器(Plan Position Indicator, PPI)，分别将MFM预处理最小二乘法径向切变、MFM预处理最小二乘法切向切变以及MFM预处理最小二乘法合成切变与本文算法对比，分析基于径向或切向的切向识别算法和本文算法在定位精度和边缘检测两个方面的风切变识别效果。

速度基数据采用阜阳新一代多普勒天气雷达(CINRAD-SA)在2009年6月5日13:26时采集的仰角为0.57°的PPI扫描数据。受冷空气影响，阜阳西北干冷气流和西南暖湿气流交汇，从北向南出现大风、雷雨和冰雹等强对流天气，瞬时风力最大达到27 m/s。

图4为方位角270°的速度图和两种算法预处理图。图4(a)为径向速度图，从图中标出区域56～59 km处和65～68 km处可以看出速度有明显的波动。图4(b)为MFM预处理后径向速度图，图4(c)为TV模型预处理后径向速度图，结合图4(b)与图4(c)，可以看出MFM与TV模型的标出区域都滤去了原速度图的波动，但TV模型的标出区域具有较好的保真度。

图5为方位角270°的速度图和两种算法切变图。图5(a)为PPI径向速度图，可看出圆圈处风速有明显的突变。结合图5(b)的MFM预处理最小二乘法合成切变图和图5(c)的本文算法切变图，对比两幅图圆圈中超过阈值的点，MFM预处理最小二乘法合成切变识别个数为2，本文算法识别个数为5，可知本文算法识别的风切变点数较多。

图6雷达速度图和两种算法切变图。图6(a)为资料速度PPI图，显示多条径向速度数据，在雷达西侧圆圈处出现了清晰的速度辐合线。图6(b)为MFM预处理最小二乘法切向切变图，从A和B圆圈处可看出检测出了一些边缘切变。图6(c)为MFM预处理最小二乘法径向切变图，从圆圈处可看出检测出了一些径向切变。图6(d)的MFM预处理最小二乘法合成切变图和图6(e)的本文算法切变图中的圆圈处均识别出了切变，最小二乘法合成切变识别的风切变在C和D之间有明显间断，本文算法识别在边缘连续性优于前者。检测风切变时，TV模型去噪较好地保留了速度基数据特征，以及8邻域系统有效地检测出4个方位切变，识别的风切变在边缘连续性上优于基于径向或者切向的风切变识别算法。

图4 方位角270°的速度图和两种算法的预处理图 Fig. 4 Diagram of velocity and speed preprocessing with two algorithms at azimuth of 270°

图5 方位角270°速度图和两种算法切变图 Fig. 5 Diagram of velocity and wind shear with two algorithms at azimuth of 270°

表1为MFM预处理最小二乘法合成切变和FCM算法切变辐合辐散的定位及误差分析。速度图中辐合辐散位置、MFM最小二乘法合成切变对应位置以及FCM算法切变对应位置可分别参考图6(a)、图6(e)以及图6(f)。通过速度图中不同辐合辐散位置的对应误差分析，MFM最小二乘法合成切变对应误差最大为1.2°、0.75 km，FCM算法切变对应误差最大为0.5°、0.25 km，可得FCM算法切变的定位精度比MFM预处理最小二乘法合成切变定位精度高。

3.2 切变强度案例分析

目前大多数风切变算法只是预警，并没有进行切变强度分析，通过龙卷风和阵风锋的案例测试，分析不同天气现象引起高切变和低切变，说明本文算法预警高低切变的可行性。

个例分析1 速度基数据采用安徽泗县新一代多普勒天气雷达(CINRAD-SA)在2006年6月29日6:44时采集的仰角为0.53°的PPI扫描数据。

泗县遭遇龙卷风袭击，图7为雷达速度图和切变强度识别图。图7(a)为资料速度PPI图，从图中标出区域可以看出东南部冷暖气流进行交汇，出现风速较大的天气现象。图7(b)为部分速度放大图，其中A处仅出现简单的对流,其风向为东南风，B处有明显的气旋存在，其高空为东南风，低空为西北风，风向变化明显。图7(c)为总切变图，从图中标出区域可以看出切变发生位置与图7(b)的对流位置相互对应。图7(d)为低切变图，在方位角150°至180°处的高空处，由于仅是简单对流，引起了低切变，切变强度范围为10.406 9～24.316 9(m/s)km。图7(e)为高切变图，雷达西南部低空气旋，引起了高切变，切变强度范围为24.413 9～54.368 1(m/s)km。

个例分析2 速度基数据采用盐城新一代多普勒天气雷达(CINRAD-SA)在2007年7月25日16:36时采集的仰角为0.62°的PPI扫描数据。

盐城测站出现阵风锋，图8为雷达速度PPI图和切变强度识别图。图8(a)为盐城雷达资料图，从图中标出区域可以看出锋面过境时有一明显折角，测站风速变化明显。图8(b)为部分速度放大图，从图中标出区域可以看出冷暖气流交汇，此时风向为东北风。图8(c)为总切变图，从图中标出区域可以看出切变强度较大，是由于测站东南部大面积对流引起。图8(d)为低切变图，从图中标出区域可以看出对流引起的低切变，结合图8(c)可以看出大部分集中在低切变区，其切变强度范围为10.400 4～22.330 5(m/s)km。图8(e)为高切变图，除低切变之外，阵风锋引起少部分强切变，其切变强度范围为22.338 1～89.110 4(m/s)km。

通过切变强度案例分析，本文算法能够判断风切变的相对强弱程度，可以为多种天气现象引起的风切变提供指示，如龙卷风、阵风锋等，对于分析不同天气现象引起的风切变具有一定可行性。

图7 雷达速度图和切变强度识别图(泗县) Fig. 7 Diagram of radar speed and shear strength identification (Sixian)

图8 雷达速度图和切变强度识别图(盐城) Fig. 8 Diagram of radar speed and shear strength identification (Yancheng) 表1 两种算法辐合辐散的定位及误差分析 Tab. 1 Location and error analysis for convergence and divergence of two algorithms

速度图中辐合辐散位置MFM预处理最小二乘法合成切变对应位置合成切变对应误差FCM算法切变对应位置切变对应误差261.0°,59.00km260.0°,58.50km1.0°,0.50km261.1°,58.75km0.1°,0.25km267.2°,62.00km267.9°,62.25km0.7°,0.25km267.0°,62.00km0.2°,0.00km274.1°,65.75km275.3°,65.00km1.2°,0.75km273.8°,65.20km0.3°,0.50km282.0°,68.00km283.0°,69.00km1.0°,1.00km282.5°,68.25km0.5°,0.25km

4 结语

本文提出一种基于FCM的低空风切变预警算法，通过对龙卷风和阵风锋引起的风切变进行分析，实现了低空风切变预警功能。算法对比和案例分析表明：1)本文算法的切变识别效果在定位精度和边缘识别两个方面均优于基于径向或切向的风切变识别算法，对于准确地判断风切变位置有重要意义。2)本文算法能够判断风切变的相对强弱程度，对于分析不同天气现象引起的风切变具有一定可行性。采用实测速度基数据进行分析，对内部信号处理引起的杂波还有待进一步处理。若能对雷达基数据的质量控制进一步完善，必能提高风切变预警效率，下一步也将会研究这些内容。

参考文献(References)

[1] CHAN P W. Application of LIDAR-based F-factor in windshear alerting [J]. Meteorologische Zeitschrift, 2012, 21(2): 193-204.

[2] HUANG T J, TIAN Y H, LI J, et al. Salient region detection and segmentation for general object recognition and image understanding [J]. Science China Information Sciences, 2011, 54(12): 2461-2470.

[3] LI N, WEI M, NIU B, et al. Application of multiple wind retrieval algorithms in nowcasting [J]. Atmosphere, 2015, 6(6): 834-849.

[4] UYEDA H, ZRNIC D S. Automatic detection of gust fronts [J]. Journal of Atmospheric and Oceanic Technology, 1986, 3(1): 36-50.

[5] EVANS J, TURNBULL D. Development of an automated windshear detection system using Doppler weather radar [J]. Proceedings of the IEEE, 1989, 77(11): 1661-1673.

[6] WOLFSON M M, DELANOY R L, FORMAN B E, et al. Automated microburst wind-shear prediction [J]. Lincoln Laboratory Journal, 1994, 7(2): 399-426.

[7] 魏耀,张兴敢.多普勒天气雷达合成切变算法及改进方法的研究[J].电子与信息学报,2010,32(1):41-47.(WEI Y, ZHANG X G. Research and improve on a method of radial velocity azimuth composite shear using Doppler weather radar [J]. Journal of Electronics & Information Technology, 2010, 32(1): 41-47.)

[8] MILLER M L, LAKSHMANAN V, SMITH T M. An automated method for depicting mesocyclone paths and intensities [J]. Weather and Forecasting, 2013, 28(3): 570-585.

[9] ZHENG J F, ZHANG J, ZHU K Y, et al. Gust front statistica characteristics and automatic identification algorithm for CINRAD [J]. Acta Meteorologica Sinica, 2014, 28(4): 607-623.

[10] TAGLIAFERRI F, VIOLA I M, FLAY R G J. Wind direction forecasting with artificial neural networks and support vector machines [J]. Ocean Engineering, 2015, 97(15): 65-73.

[11] 赵丽艳.基于多普勒天气雷达的风切变预警算法研究[D].天津:中国民航大学,2016:37-38.(ZHAO L Y. Wind shear forecasting algorithm based on Doppler weather radar [D]. Tianjin: Civil Aviation University of China, 2016:37-38.)

[12] HWANG Y, YU T Y, LAKSHMANAN V, et al. Neuro-fuzzy gust front detection algorithm with S-band polarimetric radar [J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(3): 1618-1628.

[13] DENG H Y, ZHU Q X, SONG X L. Decision-based marginal total variation diffusion for impulsive noise removal in color images [J]. Journal of Sensors, 2017: Article ID 7635641.

[14] RUDIN L I, OSHER S. Total variation based image restoration with free local constraints [C]// Proceedings of the 1994 IEEE International Conference Image Processing. Piscataway, NJ: IEEE, 1994, 1: 31-35.

[15] TABARY P, GUIBERT F, PERIER L, et al. An operational triple-PRT Doppler scheme for the French radar network [J]. Journal of Atmospheric and Oceanic Technology, 2006, 23(12): 1645-1656.

[16] LIM E, SUN J. A velocity dealiasing technique using rapidly updated analysis from a four-dimensional variational Doppler radar data assimilation system [J]. Journal of Atmospheric and Oceanic Technology, 2010, 27(7): 1140-1152.

[17] JAIN A K. Data clustering: 50 years beyond K-means [J]. Pattern Recognition Letters, 2010, 31(8): 651-666.

[18] 周双,冯勇,吴文渊,等.一种基于模糊C均值聚类小数量计算最大Lyapunov指数的新方法[J].物理学报,2016,65(2):020502.(ZHOU S, FENG Y, WU W Y, et al. A novel method based on the fuzzy C-means clustering to calculate the maximal Lyapunov exponent from small data [J]. Acta Physica Sinica, 2016, 65(2): 42-48.)

[19] 黄仪方,朱志愚.航空气象[M].成都:西南交通大学出版社,2002:131-134.(HUANG Y F, ZHU Z Y. Aviation Weather[M]. Chengdu: Southwest Jiaotong University Press, 2002:131-134.)

This work is partially supported by the National Natural Science Foundation of China (U1533113, U1433202), the Fundamental Research Funds for the Central Universities (3122016B001).

XIONGXinglong, born in 1962, M. S., professor. His research interests include signal and information processing, LIDAR (Laser Intensity Direction And Ranging) weather detection.

YANGLixiang, born in 1990, M. S. candidate. Her research interests include weather radar weather detection, image processing.

MAYuzhao, born in 1978, Ph. D., associate professor. Her research interests include fiber optics, aviation weather detection, electromagnetic computation.

ZHUANGZibo, born in 1980, M. S., lecturer. His research interests include aviation weather detection, data processing.