基于多普勒天氣雷達(dá)的低空多普勒速度的切變識別算法研究

閆文輝，黃興友，李盈盈，楊敏，劉燕斐

(南京信息工程大學(xué)氣象災(zāi)害預(yù)報預(yù)警與評估協(xié)同創(chuàng)新中心，江蘇南京210044)

1 引 言

低空風(fēng)切變是指600 m高度以下空間兩點風(fēng)的矢量差，即風(fēng)向、風(fēng)速在空間兩點突然大幅度變化的一種現(xiàn)象[1]。低空風(fēng)切變是威脅航空安全的重要原因之一，特別是在飛機(jī)起飛和降落階段，由于飛機(jī)在低層大氣遭遇風(fēng)切變之后，空速、升力驟變，導(dǎo)致飛機(jī)失速和操縱困難，甚至?xí)l(fā)生墜機(jī)。國際民航組織2005年進(jìn)行的一項調(diào)查顯示，1964—1983年，低空風(fēng)切變至少造成了全球28次大型民航飛機(jī)事故，500多人死亡，200人受傷[2]。2015年2月臺灣復(fù)興航空公司的飛機(jī)在機(jī)場降落過程中，遭遇低空風(fēng)切變，導(dǎo)致機(jī)上乘客與機(jī)組人員48人死亡，10人受傷[3]。因此，研究出快速可靠的低空風(fēng)切變識別預(yù)警算法，是降低飛行風(fēng)險的一個有效措施。

在低空風(fēng)切變識別方面，國內(nèi)外開展了一系列研究。Wilson等[4]利用多普勒雷達(dá)數(shù)據(jù)分析了微下?lián)舯┝黠L(fēng)場的水平和垂直結(jié)構(gòu)，并提出利用多普勒雷達(dá)進(jìn)行微下?lián)舯┝髯R別和低空風(fēng)切變?yōu)暮︻A(yù)警的可能性。Uyeda等[5]第一次開發(fā)了一種基于徑向速度梯度的陣風(fēng)鋒識別算法，Campbell等[6]研制了一種基于人工智能的低空風(fēng)切變自動識別方法。Evans等[7]描述了集成的業(yè)務(wù)化風(fēng)切變預(yù)警系統(tǒng)，該系統(tǒng)一般由激光雷達(dá)風(fēng)切變預(yù)警系統(tǒng)(LIWAS)和終端多普勒天氣雷達(dá)系統(tǒng)(TDWRS)構(gòu)成，已在美國幾十個機(jī)場實施，效果明顯，香港機(jī)場也采用了類似的低空風(fēng)切變預(yù)警系統(tǒng)。此外，Boilley等[8]利用高分辨率的數(shù)值模式對Nice-Cote d'Azur機(jī)場的低空風(fēng)切變進(jìn)行預(yù)報。Augros等[9]將多普勒雷達(dá)徑向風(fēng)速投影到笛卡爾坐標(biāo)系中，對比鄰近格點風(fēng)速的差別來識別風(fēng)切變。

國內(nèi)學(xué)者對風(fēng)切變的識別和預(yù)警也進(jìn)行了不少研究，段鶴等[10]詳細(xì)介紹了典型颮線發(fā)生時多普勒速度回波特征及風(fēng)切變的位置。王珊珊[11]基于直接計算差值濾波算法識別了低空風(fēng)切變，并以動態(tài)鏈接庫的方式實現(xiàn)了工程設(shè)計。王楠等[12]采用最小二乘法識別了低空風(fēng)切變和輻合線，并討論了擬合“窗口”大小(即參與計算的資料點數(shù)n)對資料預(yù)處理效果和梯度計算的影響。魏耀等[13]提出了擬合“窗口”隨探測距離遠(yuǎn)近而變的最小二乘法低空風(fēng)切變識別算法，發(fā)現(xiàn)其優(yōu)于直接計算差值濾波算法。胡琦[14]利用香港國際機(jī)場多普勒激光雷達(dá)數(shù)據(jù)，使用灰色預(yù)測、BP神經(jīng)網(wǎng)絡(luò)等4種算法，得到了預(yù)測風(fēng)速圖，并用于低空風(fēng)切變的預(yù)警。Wang等[15]基于區(qū)域增長算法識別了不同高度水平風(fēng)的水平和垂直切變。周生輝等[16]利用VVP風(fēng)場反演的方法進(jìn)行了低空風(fēng)切變識別分析。蔣立輝等[17]針對新一代多普勒天氣雷達(dá)單尺度的風(fēng)切變識別方法容易產(chǎn)生漏報的問題，提出了基于自適應(yīng)多尺度的梯度搜索風(fēng)切變預(yù)警算法。杜愛軍[18]采用最小二乘法和退模糊技術(shù)，計算了基于南京信息工程大學(xué)C波段雙偏振雷達(dá)資料的二維合成風(fēng)切變。國內(nèi)還有較多相應(yīng)的研究[19-23]。

低空風(fēng)切變是在不同的天氣背景和環(huán)境條件下形成的，因而具有不同的空間(水平)尺度。但國內(nèi)的低空風(fēng)切變識別算法都是通過單一尺度的一維徑向風(fēng)切變、一維切向風(fēng)切變和二維合成風(fēng)切變來判別低空風(fēng)切變的位置、強(qiáng)度。這些算法不能很好地識別各種尺度的低空風(fēng)切變，容易產(chǎn)生漏報現(xiàn)象。而且擬合“窗口”大小的選取帶有人為主觀性，會影響低空風(fēng)切變識別結(jié)果。此外，一維徑向切變值隨著擬合“窗口”的增大而減小，因為窗口越大，平均作用越強(qiáng)。擬合“窗口”過小，計算結(jié)果會受到預(yù)處理后殘留噪聲的影響；擬合“窗口”過大，可能會將某些大切變區(qū)漏掉。本文的自適應(yīng)多尺度低空風(fēng)切變識別算法，可有效解決“窗口”大小的問題。

本文對已有風(fēng)切變識別算法中存在的單一尺度問題進(jìn)行了分析，根據(jù)風(fēng)切變強(qiáng)度因子，自適應(yīng)選擇合適大小的擬合“窗口”，采用最小二乘法，自動處理二維合成風(fēng)切變。并與我國新一代S波段多普勒天氣雷達(dá)(CINRAD/SA)的二維合成風(fēng)切變產(chǎn)品進(jìn)行了對比分析。此外，本文還計算了垂直風(fēng)切變和時間風(fēng)切變，實現(xiàn)了低空風(fēng)切變的四維識別預(yù)警，為飛行安全提供保障。

2 低空風(fēng)切變識別算法

2.1 自適應(yīng)擬合“窗口”選取

風(fēng)切變強(qiáng)度因子是Woodfield等[24]提出的，它是衡量低空風(fēng)切變強(qiáng)度的有效標(biāo)準(zhǔn)之一，其計算公式如下：

式中dv/dt是風(fēng)速變化率(沿著擬合路徑)，Δv是風(fēng)速變化量(沿著擬合路徑)，Vapp是飛機(jī)的正常進(jìn)近速度，H表示擬合路徑尺度。最大風(fēng)切變強(qiáng)度因子所在的局部區(qū)域一般最容易發(fā)生低空風(fēng)切變。

根據(jù)飛機(jī)的響應(yīng)時間，風(fēng)切變強(qiáng)度因子的計算“窗口”分別選為0.5 km、1 km、2 km和4 km四個尺度。在同一距離庫處，沿徑向分別計算四個尺度下的風(fēng)切變強(qiáng)度因子值。比較四個風(fēng)切變強(qiáng)度因子值，取最大因子值所對應(yīng)的尺度作為該距離庫一維徑向風(fēng)切變計算的自適應(yīng)擬合“窗口”。以此得到每個距離庫的一維徑向低空風(fēng)切變計算的自適應(yīng)擬合尺度。同理，在同一探測距離上、得到距離庫的一維切向低空風(fēng)切變計算的自適應(yīng)擬合窗口尺度。

2.2 基于最小二乘法的二維合成風(fēng)切變算法

2.2.1 一維徑向風(fēng)切變

對于一維徑向風(fēng)切變(RS)，自適應(yīng)擬合“窗口”包含n個距離庫，即：

式中vi是“窗口”中第i個距離庫的徑向速度，ri是相應(yīng)點與雷達(dá)之間的徑向距離，利用最小二乘法計算一維徑向風(fēng)切變的方法如下。

設(shè)線性回歸方程為：

式中v*是擬合“窗口”速度估計值，a為初始值，b為斜率，r是徑向距離。

則估計值與實測值之間的誤差為：

根據(jù)最小二乘法原理，使上述誤差平方的和最小，即使下式取最小值。

分別對a、b求偏導(dǎo)并令其等于0，即：

解式(6)得：

如圖1所示，A點和D點之間的一維徑向風(fēng)切變是(vD-vA)與資料點總長度ΔR之比。而按照式(5)，斜率b是(vD-vA)與A、D兩點之間的距離Δr的比值。則一維徑向風(fēng)切變的大小為：

式中LR是資料點距離庫在徑向上的單位長度，即雷達(dá)徑向分辨率。

圖1 一維徑向風(fēng)切變原理示意圖

2.2.2 一維切向風(fēng)切變

對于一維切向風(fēng)切變(AS)，自適應(yīng)擬合“窗口”包含n個距離庫，即：

式中vi是“窗口”內(nèi)第i個距離庫的徑向速度，θi是相應(yīng)點的方位角。設(shè)v*為擬合“窗口”的速度估計值，則線性回歸方程為v*=a+b(rθ)，r為該距離庫與雷達(dá)之間的距離，仍利用最小二乘法得擬合曲線的斜率b為：

距離庫內(nèi)的數(shù)據(jù)是以該點為中心的采樣體積內(nèi)的平均值，采樣體積與雷達(dá)波束寬度和脈沖寬度相關(guān)。如圖2所示，相同探測距離上的A、D兩點的徑向速度的一維切向風(fēng)切變?yōu)?

式中θR為雷達(dá)的角度分辨率(CINRAD/SA雷達(dá)的角度分辨率約為1°)，Δθ為弧度。

由于最小二乘法計算得到的斜率b為：

則一維切向低空風(fēng)切變值A(chǔ)S為：

2.2.3 二維合成風(fēng)切變

二維合成風(fēng)切變定義為徑向速度在徑向、切向的綜合變化，即一維徑向風(fēng)切變和一維切向風(fēng)切變的合成。

2.3 垂直風(fēng)切變算法

垂直風(fēng)切變定義為多普勒雷達(dá)徑向速度在距雷達(dá)相同距離，相同方位角處(即r、θ都是常數(shù))上下兩層PPI上隨垂直高度的變化[12]，即：

式中，vp、vq為相鄰上下兩層PPI的徑向速度，并進(jìn)行了“窗口”為n的平滑擬合，Rp、Rq分別為上下兩層的徑向距離；αp、αq分別為上下兩層PPI的仰角。當(dāng)上下兩層仰角的徑向速度數(shù)據(jù)有任一缺失時不進(jìn)行垂直切變計算，該點作為無資料處理。本文中“窗口”n取為3點，上下兩層PPI徑向速度分別取最低兩層仰角(1.45 °和 0.53 °)。

2.4 時間風(fēng)切變算法

基于多普勒雷達(dá)徑向速度的時間風(fēng)切變是指同一空間位置 (即徑向距離R、仰角α和方位角θ都為常數(shù))處時間間隔為Δt的前后兩次低仰角PPI掃描徑向速度隨時間的變化，單位為m/(s·min)，即：

式中，vt1、vt2分別為同一空間位置處時間間隔Δt的前后兩次PPI掃描的徑向速度。當(dāng)時間間隔為Δt的兩次PPI數(shù)據(jù)中徑向速度有任一缺失時不進(jìn)行時間切變計算，該點作為無資料處理。由于新一代多普勒天氣雷達(dá)VCP21模式下一次完整體掃需要6 min，為了計算第n-2次和第n次體掃的時間風(fēng)切變，時間間隔Δt取12 min。

3 實例分析

本文利用兩次天氣過程的新一代多普勒天氣雷達(dá)(CINRAD/SA)資料及MICAPS高空資料，進(jìn)行風(fēng)切變識別算法研究和驗證分析。CINRAD/SA雷達(dá)的波束寬度為1.0°、探測距離為230 km、反射率因子數(shù)據(jù)的庫長為1 km、徑向速度數(shù)據(jù)的庫長為250 m。在計算風(fēng)切變之前，需對雷達(dá)基數(shù)據(jù)進(jìn)行質(zhì)量控制，以消除并糾正基數(shù)據(jù)中的錯誤。采用的方法是孤立噪聲點[25]濾除、缺測數(shù)據(jù)填補(bǔ)和濾波平滑。本文使用鄭佳鋒[26]提出的“鄰域頻數(shù)法”的噪聲點濾除與缺測值填補(bǔ)算法對反射率因子、平均多普勒速度進(jìn)行處理?？紤]到氣象回波的連續(xù)性，為了進(jìn)一步減少徑向速度噪聲，同時有效地保留徑向速度中的小尺度結(jié)構(gòu)，先對徑向速度數(shù)據(jù)進(jìn)行了3×3中值濾波，再沿徑向和切向分別作5點、3點滑動平均。

3.1 實例1

2009年6月3日傍晚—4日凌晨，徐州及周邊地區(qū)颮線過境，遭遇大范圍雷雨、大風(fēng)和冰雹等強(qiáng)對流天氣，19個縣市出現(xiàn)了17 m/s以上的短時大風(fēng)。沙氏指數(shù)SI為-10.02℃，氣團(tuán)指標(biāo)K為38℃，大氣低層溫度露點差較大，十分有利于強(qiáng)對流的發(fā)生，地面可能出現(xiàn)強(qiáng)低空風(fēng)切變。

圖3為2009年6月3日23:30(北京時間，下同)徐州站多普勒天氣雷達(dá)回波圖、PUP合成風(fēng)切變識別產(chǎn)品及本文方法得到的時空四維風(fēng)切變(二維合成風(fēng)切變、垂直風(fēng)切變及時間風(fēng)切變)識別圖。從0.53°仰角雷達(dá)回波圖(圖3a、3c)中可知，“人”字形颮線系統(tǒng)正在過境，颮線附近表現(xiàn)為活躍的雷暴帶，這一區(qū)域反射率因子最大(45～60 dBz)，雷達(dá)站附近徑向速度大于20 m/s。圖中右邊無層狀云區(qū)的暖式颮線從雷達(dá)站周圍向東延伸，左下方有層狀云區(qū)的冷式颮線從雷達(dá)站向南延伸，相對于右邊的暖式颮線，左下方的冷式颮線具有更加明顯的颮線特征，反射率因子回波圖中有明顯的陣風(fēng)鋒細(xì)線回波，徑向速度圖中的陣風(fēng)鋒附近有明顯的強(qiáng)切變線存在。

PUP合成風(fēng)切變產(chǎn)品是利用1.45°仰角徑向速度數(shù)據(jù)，基于直接計算差值濾波算法得到的。該算法先用差值法計算徑向風(fēng)切變和切向風(fēng)切變，再經(jīng)過濾波得到合成風(fēng)切變[13]。為了作對比，本文的自適應(yīng)多尺度合成風(fēng)切變也采用1.45°仰角徑向速度數(shù)據(jù)。圖3e和圖3f分別為同一時刻PUP合成風(fēng)切變產(chǎn)品和本文的自適應(yīng)多尺度合成風(fēng)切變識別結(jié)果，顯示閾值都為2.0×10-3s-1?？梢姡铣娠L(fēng)切變大值區(qū)都出現(xiàn)在冷式颮線的前沿，在暖式颮線附近都沒有明顯的合成風(fēng)切變，這與1.45°仰角徑向速度圖(圖3d)中由陣風(fēng)鋒引起的風(fēng)切變位置相吻合。相對于PUP合成風(fēng)切變產(chǎn)品，本文的自適應(yīng)多尺度合成風(fēng)切變識別結(jié)果能有效地減少線性擬合的誤差，濾除由于噪聲導(dǎo)致的虛假合成風(fēng)切變，回波邊緣風(fēng)切變識別更加合理，識別的風(fēng)切變也更加連續(xù)。

圖3g是0.53°仰角徑向速度與1.45°仰角徑向速度間的垂直風(fēng)切變?？梢?，在左下方的冷式颮線后部出現(xiàn)負(fù)垂直切變最大值，冷式颮線前部出現(xiàn)正垂直切變最大值。說明在低層冷式颮線后部存在和颮線移動方向相同的徑向速度，冷式颮線前部存在和颮線移動方向相反的徑向速度。并且在陣風(fēng)鋒附近垂直風(fēng)切變的水平梯度很大，在一定程度上也能反映颮線的位置。圖3h是2009年6月3日23:18和2009年6月3日23:30之間0.53°仰角徑向速度的時間風(fēng)切變，顯示閾值為0.5 m/(s·min)。可見，時間風(fēng)切變大值區(qū)出現(xiàn)在陣風(fēng)鋒在12 min內(nèi)移動經(jīng)過的區(qū)域，所以時間切變可反映颮線在12 min內(nèi)的移動情況。

綜合分析圖3可看到，合成風(fēng)切變、垂直風(fēng)切變和時間風(fēng)切變構(gòu)成的時空四維風(fēng)切變識別圖中，切變發(fā)生的位置與雷達(dá)徑向速度圖中颮線的位置有較好的對應(yīng)關(guān)系，風(fēng)切變識別結(jié)果較理想，可為民航機(jī)場低空風(fēng)切變預(yù)警提供一定的幫助。

圖3 2009年6月3日23:30徐州颮線過程多普勒天氣雷達(dá)回波及低空風(fēng)切變識別結(jié)果

3.2 實例2

2017年6月5日西南低空急流造成南京及周邊地區(qū)出現(xiàn)大范圍大風(fēng)、降水天氣。圖4a是6月5日08時850 hPa天氣圖?？梢姡暇┑貐^(qū)位于西南急流中心的左前方，該區(qū)域由于渦度平流引起的垂直運動造成上升運動最強(qiáng)，并且該區(qū)域與溫度場中的暖舌相結(jié)合。圖4b是6月5日08時南京站單站高空風(fēng)圖，從圓心到該點的距離表示風(fēng)速大小，從圓心指向該點的矢量表示風(fēng)的來向，可見在850 hPa附近存在風(fēng)速的極大值，大于20 m/s，達(dá)到低空急流的風(fēng)速標(biāo)準(zhǔn)，風(fēng)向隨高度增加順轉(zhuǎn)，出現(xiàn)暖平流。

圖4 2017年6月5日08:00 850 hPa天氣圖(a)和南京站單站高空風(fēng)圖(b)

圖5 (見下頁)為2017年6月5日12:01南京站多普勒天氣雷達(dá)回波圖、PUP合成風(fēng)切變識別產(chǎn)品和本文方法得到的時空四維風(fēng)切變 (二維合成風(fēng)切變、垂直風(fēng)切變及時間風(fēng)切變)識別圖。強(qiáng)降水區(qū)域的回波強(qiáng)度為20～45 dBz，沒有明顯的強(qiáng)對流。徑向速度圖(圖5b)中，零徑向速度線呈正“S”形，風(fēng)向隨高度順轉(zhuǎn)，風(fēng)速隨高度增加先增加后減小，在雷達(dá)站西北方向46～92 km附近存在正徑向速度中心，在雷達(dá)站東南方向46～92 km附近存在負(fù)徑向速度中心，在距雷達(dá)站80 km(高度約2 km，即850 hPa等壓面)附近存在明顯的風(fēng)切變。這與上述的天氣學(xué)分析結(jié)果基本吻合。

圖5c和圖5d分別是1.45°徑向速度對應(yīng)的PUP合成風(fēng)切變產(chǎn)品和本文的自適應(yīng)多尺度合成風(fēng)切變識別結(jié)果，顯示閾值都為2.0×10-3s-1。切變大值區(qū)與徑向速度圖中切變線位置吻合良好，本文方法得到的合成風(fēng)切變相對于PUP的合成風(fēng)切變產(chǎn)品，有效地濾除徑向速度自身波動引起的虛假風(fēng)切變，回波邊緣識別結(jié)果更可靠。在切變線西側(cè)存在負(fù)垂直切變大值區(qū)(圖5e)，說明上層負(fù)徑向速度增強(qiáng)，正徑向速度減弱。圖5f是2017年6月5日11:49和2017年6月5日12:01之間1.45°仰角徑向速度的時間風(fēng)切變，顯示閾值為0.5 m/(s·min)?？梢?，時間風(fēng)切變大值區(qū)同樣和徑向速度圖中切變線的位置對應(yīng)較好，相對于颮線過程，低空急流引起的時間風(fēng)切變值較小，說明該天氣過程相對于颮線系統(tǒng)移動較慢。

4 結(jié)論與討論

本文研究了利用地基多普勒天氣雷達(dá)基數(shù)據(jù)識別低空風(fēng)切變的方法，并通過徐州一次颮線過程及南京一次混合云降水過程的實例從時空四維空間上驗證了低空風(fēng)切變識別算法的有效性。

(1)對雷達(dá)徑向速度資料進(jìn)行一定尺度的濾波和平滑，可有效地消除徑向速度自身波動產(chǎn)生的噪聲，并保留產(chǎn)生低空風(fēng)切變的中小尺度風(fēng)信息。

圖5 2017年6月5日12:01南京混合云降水過程多普勒天氣雷達(dá)回波及低空風(fēng)切變識別結(jié)果

(2)基于風(fēng)切變強(qiáng)度因子的最小二乘法自適應(yīng)擬合“窗口”的選取可有效地減少人為主觀性，降低單一窗口尺度造成的低空風(fēng)切變漏報現(xiàn)象。

(3)相對于PUP二維合成風(fēng)切變產(chǎn)品，本文方法得到的合成風(fēng)切變在切變連續(xù)性、切變識別效果和邊緣虛假切變處理等方面更優(yōu)，與雷達(dá)徑向速度圖中風(fēng)切變的位置相吻合。

(4)垂直風(fēng)切變可反映徑向速度的高低空配置情況，時間風(fēng)切變可提供徑向速度隨時間的變化情況。

利用多普勒雷達(dá)徑向速度數(shù)據(jù)從合成風(fēng)切變、垂直風(fēng)切變和時間風(fēng)切變四維空間實現(xiàn)對低空風(fēng)切變的自動識別，可為民航機(jī)場低空風(fēng)切變預(yù)警提供有效的輔助工具。需要指出的是，由于多普勒雷達(dá)探測的是取樣體積內(nèi)的平均徑向風(fēng)，并且空間分辨率有限，很難準(zhǔn)確反映真實的三維風(fēng)場特性，因此，這里識別的低仰角徑向速度的低空風(fēng)切變是真實風(fēng)場的低空風(fēng)切變的近似，兩者之間可能稍有差異。另外，本文未考慮速度模糊現(xiàn)象，算法的預(yù)處理過程還值得進(jìn)一步改進(jìn)和優(yōu)化。