基于模糊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.