機載雷達風切變識別算法研究及在機場預報中的應用

魏鳴，張明旭，2，張培昌，郭巍，周生輝，2

(1．南京信息工程大學氣象災害預報預警與評估協同創新中心，江蘇南京210044;2．南京信息工程大學中國氣象局氣溶膠與云降水開放重點實驗室，江蘇南京210044;3．上海市衛星遙感與測量應用中心，上海201199)

0 引言

隨著我國航空事業的發展，減少危險天氣對飛行的威脅日益迫切。在愈加頻繁的氣象危害中，大氣的風切變特別是低空風切變對飛行的影響最大，造成的事故所占比例高(黃儀方，2002)。風切變形成機制復雜，持續時間短且強度較大(武虎子等，2013)。雷暴等強對流天氣在一定范圍內可產生以垂直風為主要特征的綜合風切變區域，尤其以雷暴云體的強下沉氣流區和冷性外流氣流前緣引起的陣風鋒更為嚴重，對飛行危害最大(張培昌等，2001;王叢梅等，2011)。鋒面天氣可產生以水平風的水平切變和垂直切變為主的風切變區域，對飛行的危害程度較之強對流天氣略弱(孫一昕和方娟，2012)。Williams(2005)提到輻射逆溫型的低空急流天氣可產生逆溫風切變，由于其強度通常較小，容易被忽視，但若處置不當也會發生飛行危險;機場周圍地物、建筑等也會引起風切變現象，且強度較大，這使機場的選址尤為重要(秦娟等，2011;遲繼峰等，2012)。因此及時準確地判斷風切變的位置、類型和強度，可有效減輕和避免危害，確保飛行及起降安全。

Chan and Shun(2010)研究認為風切變的強度標準涉及氣象條件、飛機性能及飛行員駕駛技術等多方面因素，因此以其對飛行的危害程度作為強度分類標準成為航空界共識。按照國際民航組織建議(章澄昌，2000)，水平風的垂直切變強度標準見表1，其中空氣層垂直厚度為30 m，時間尺度為2 min。當切變值達到0.1(1/s)以上時就威脅飛機安全。水平風的水平切變標準一般依據美國機場低空風切變報警系統所采用的報警標準(Britt et al．，1993)，水平風水平切變值在2.6(m/s)/km以上時極有可能對飛行構成危害。

表1 水平風垂直切變強度的標準Table 1 The standard ofverticalshear strength for horizontal wind

多普勒天氣雷達不僅具有多普勒徑向速度信息，還有反映風場變化的速度譜寬信息。在多普勒雷達資料的各種算法產品中，合成切變產品可反映大氣流場的脈動性，同時可以展示3～10 km尺度以內危險的陣風鋒、低空切變線等低空風切變現象以及與中尺度氣旋有關的風切變識別產品(張京英等，2009;梁建宇和孫建華，2012)。在風切變識別方面，國內外有許多研究進展。目前主要的識別算法有直接計算差值濾波法、最小二乘法等(魏耀和張興敢，2010)。Harris et al．(1985)首次提出了直接計算差值濾波合成切變算法。該算法在笛卡爾坐標系中分別計算一維徑向風切變和一維方位風切變，然后經過求平均及濾波等步驟，合成為二維風切變。合成切變產品用于監測陣風鋒、微下擊暴流以及冷鋒過境等風切變，并經逐步改進與機載雷達相結合，建立了早期的機載低空風切變預警系統。METSTAR公司生產的CINRAD-SA/SB等多普勒雷達均采用基于美國NEXTRAD的合成切變算法，但國內版權受限。在參考WSD-88D直接計算差值濾波合成切變算法的基礎上，王珊珊等(2008)進行了雷達回波自動識別鋒線的方法研究，產品能夠直觀識別風切變。胡明寶等(2000)用機場多普勒雷達資料根據最小二乘法計算合成切變值，定量描述機場附近的低空風切變區域，確定飛機起降的潛在危險區。2008年北京奧運會期間，氣象保障部門應用最小二乘法獲得的合成切變產品為人工影響天氣作業提供了重要參考。魏耀和張興敢(2010)對直接計算差值濾波和最小二乘法作對比分析，后者在定位精度、識別能力及邊緣識別能力等方面具有優勢。

機載多普勒雷達的風切變識別算法需要具有快速、簡捷和高效的特點，本文根據機載氣象雷達的特點和應用需求，采用最小二乘法的識別原理，通過實例分析，將此算法與目前多普勒天氣雷達業務應用系統PUP上的風切變識別產品進行對比分析，討論它在機載雷達上應用的可行性，為改進機載雷達風切變預警系統提供參考。

1 最小二乘法的識別算法原理

1.1 一維徑向風切變

為了減少雷達測量噪聲和奇異點的影響，選取合適的擬合窗口，其中包含n個距離庫(同一徑向上)，即(v1，r1)…(vi，ri)…(vn，rn)(vi為窗口中第 i個點的徑向速度，ri為第i個點與雷達之間的距離)，然后利用最小二乘法計算雷達徑向速度沿雷達徑向的變化值CRs。設v*是擬合窗口的速度估計值，則線性回歸方程為v*=α+βr。其中α、β是待定系數，α為初始值，β為斜率。

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

根據最小二乘法原理，使上述誤差平方和最小，即令E最小，

則

解(3)式得

如圖1所示，A、C兩點之間的徑向切變為(vC－vA)與資料點總長度ΔR之比，而斜率β是(vC－vA)與A、C兩點之間距離Δr的比值。則計算一維徑向切變值:

其中:LR是資料點采樣體積在徑向上的單位長度。

圖1 徑向切變示意圖Fig．1 The illustration of radial shear

1.2 一維方位風切變

對于一維方位風切變(CAs)，在同一探測距離圓周上選取擬合窗口，計算過程中采取的擬合窗口內包含 n 個距離庫，即(v1，θ1)…(vi，θi)…(vn，θn)(vi為窗口內第i個點的徑向速度，θi為相應點的方位角)。由最小二乘法得擬合曲線的斜率β:

距離庫內的數據是以該點為中心的采樣體積內的平均值，采樣體積與波束寬度和距離庫長相關。如圖2所示，同一庫距內A、B兩點徑向速度的一維方位風切變值CAs(雷達的角度分辨率為1.5°)為:

圖2 方位切變示意圖Fig．2 The illustration of azimuth shear

由最小二乘法計算得 β=(vB－vA)/(rΔθ)，一維方位風切變值CAs為:

1.3 二維合成風切變

雷達徑向速度沿徑向、方位的綜合變化，即一維徑向風切變與一維方位風切變的合成，得到二維合成切變:

2 雷達資料處理

2.1資料質量控制

為保證識別效果，計算前需對雷達資料進行質量控制，包括雷達缺失數據填補及去除孤立噪聲點等。缺失數據主要是雷達掃描時數據重置和原始數據丟失所致(Lakshmanan et al．，2003;袁杲等，2007)。為填補缺失徑向數據，先判斷方位角索引，若與上次相同則另外保存。讀完后再判斷是否有突變，若沒有則不需進一步處理，否則就再對原始回波進行修正。去除孤立噪聲時，以當前回波點為中心，取5×5的模板，計算無效反射率因子值所占百分比P，若P大于75%，則判斷該點為孤立回波點，需去除(黃鈴光，2011)。

2.2 擬合窗口選取

選取擬合窗口(表2)時，要兼顧隨機載雷達探測距離遠近而變化的徑向切變和方位切變，以最大程度減少線性擬合的誤差，因此為使合成切變的擬合精度不隨距離而變化，規定組合擬合窗口面積保持不變。結合實驗分析，資料點的選取應滿足:

式中:Nr為徑向切變擬合窗口資料點;Nt為方位切變擬合窗口資料點;ε為常數。

表2 徑向切變、方位切變擬合窗口資料點數的選取Table 2 The data points of fitting window for radial shear and azimuth shear

鑒于機載雷達掃描方式、掃描半徑與地基雷達的差別，在處理機載雷達資料時，擬合窗口的選取需要根據實際參數調整。

3 實例分析

3.1 石家莊機場風切變

1)個例1:2010年6月16—17日

2010年6月16日傍晚至17日下半夜，石家莊機場本場及周邊地區出現大范圍雷雨、短時暴雨、大風和冰雹等強對流天氣，石家莊、滄州以北地區28個縣(市)出現大風，萬全地區瞬時風力最大達24 m/s;河北省除邢臺東部、邯鄲東部外均出現雷雨天氣，中北部地區有33個縣市降水量達25～106 mm;有8個測站出現冰雹，寬城測站冰雹最大直徑為14 mm。這也是當年出現的范圍最大、強度最強的對流天氣過程。

圖3是2010年6月16日20時(北京時，下同)的850 hPa天氣圖?？梢姡仪f地區處于北部低壓中心的東南側，位于槽前脊后，降水量在20時前后并不多，但風速較大，東北西南走向的切變較明顯。由于石家莊沒有探空站，只能參考周邊山西太原(53772)和河北邢臺(53798)的探空資料。圖4是2010年6月16日20時太原和邢臺的溫度對數壓力(T-lgp)圖，可見兩站的對流有效位能CAPE值并不大，且水汽不充足，但此時850、700、500 hPa的風速較大，太原依次為 10、8、13 m/s，邢臺依次為11、6、8 m/s，且邢臺的 925 hPa 風速達 19 m/s，由此推斷石家莊地區風速較大，風切變明顯。

圖3 6月16日20時850 hPa天氣圖Fig．3 850 hPa synoptic chart at 20:00 BST 16 June

圖5為石家莊2010年6月16日23:54多普勒雷達徑向速度資料圖及PUP合成切變識別圖和最小二乘法合成切變識別圖?？梢姡瑘D5b圓圈區域識別得模糊，而圖5c相應區域則識別得清晰連續。

2)個例2:2010年7月10日

圖46月16日20時太原(a)和邢臺(b)的T-lgp圖Fig．4 T-lgp plots of(a)Taiyuan and(b)Xingtai at 20:00 BST 16 June

圖5 6月16日23:54雷達徑向速度及合成切變圖(半徑:115 km;仰角:1.5°) a．PPI的雷達徑向速度;b．PUP合成切變圖;c．最小二乘法合成切變圖Fig．5 Radar radial velocity and composite shear at 23:54 BST 16 June(radius:115 km;elevation:1.5°) a．radial velocity on PPI;b．composite shear of PUP product;c．composite shear of least square method

2010年7月10日石家莊及周邊地區出現大范圍雷雨、短時暴雨、大風和冰雹等強對流天氣，石家莊機場本場20時左右遭遇大風天氣，持續近半小時，最強陣風達20 m/s，降雨強、風速大給當地經濟和人員造成重大損失。

圖6是2010年7月10日08時850 hPa天氣圖?？梢?，河北省處于低壓中心東南側，西北干冷空氣與東南暖濕空氣在石家莊附近相遇并出現降水。20時低槽略東移，石家莊處于槽前脊后，降水增加，強對流明顯(圖略)。圖7是20時太原和邢臺的T-lgp圖。可見，太原的對流有效位能 CAPE為658.3 J/kg，對流抑制能CIN為149.1 J/kg，邢臺的對流有效位能 CAPE為256.9 J/kg，對流抑制能CIN 為39.4 J/kg。850、700、500 hPa的風速依次為:太原 4、7、16 m/s，邢臺 1、10、15 m/s。分析表明石家莊將會繼續存在風切變、降水和強對流天氣。

圖8為石家莊2010年7月10日21:24多普勒雷達徑向速度資料圖及PUP合成切變識別圖和最小二乘法合成切變圖?？梢?，圖8b圓圈所示的弱切變區識別得較零散，而圖8c則識別得較連續。

圖6 7月10日08時850 hPa天氣圖Fig．6 850 hPa synoptic chart at 08:00 BST 10 July

圖7 2010年7月10日20時太原(a)和邢臺(b)的T-lgp圖Fig．7 T-lgp plots of(a)Taiyuan and(b)Xingtai at 20:00 BST 10 July 2010

3)個例3:2009年6月5日

2009年6月5日午后安徽省大部出現雷雨大風及冰雹等強對流天氣，至晚間上半夜全省40余縣市出現8級以上大風，其中淮南20時53分最大瞬時風速達36 m/s，阜陽地區監測到強烈的陣風鋒，懷遠、淮南分別出現8、10 mm直徑的冰雹。強對流引起大風天氣在當地造成生命財產的嚴重損失。

圖8 7月10日21:24雷達徑向速度及合成切變圖(半徑:115 km，仰角:1.5°) a．PPI的雷達徑向速度;b．PUP合成切變圖;c．最小二乘法合成切變圖Fig．8 Radar radial velocity and composite shear at 21:24 BST 10 July(radius:115 km;elevation:1.5°) a．radial velocity on PPI;b．composite shear of PUP product;c．composite shear of least square method

圖9是2009年6月5日08時和20時的天氣圖。可見，西北干冷氣流和西南暖濕氣流在安徽交匯，該地區出現強對流并伴有雷暴天氣。阜陽處于安徽北部，有西北—東南方向的強切變區，兩時次天氣圖顯示切變區穩定并維持，雷達速度圖上切變帶清晰可見。圖10是20時阜陽T-lgp圖?？梢姡瑘D中紅色區域所占比例較大，對流有效位能CAPE達1 782.7 J/kg;對流抑制能CIN為171.7 J/kg;溫度露點僅為 3 ℃;925、850、700、500 hPa的風速依次為13、13、12、21 m/s，從地面到高空風速都很大，使得阜陽地區維持強對流天氣和風切變。

圖11為阜陽2009年6月5日21:44多普勒雷達徑向速度資料圖及PUP合成切變圖和最小二乘法合成切變圖。由圖11b(圓圈處)和圖11c均可見雷達西南側存在陣風鋒，而圖11c識別的陣風鋒邊界比圖11b更加清晰，說明最小二乘法識別的梯度強，邊緣信息保留完整。

圖10 6月5日20時阜陽T-lgp圖Fig．10 T-lgp plot in Fuyang at 20:00 BST 5 June

圖96月5日08時(a)和20時(b)850 hPa天氣圖Fig．9 850 hPa synoptic charts at(a)08:00 and(b)20:00 BST 5 June

3.2 識別結果分析

為進一步定量分析兩種方法的識別效果，選取輻合輻散性氣流(合成切變值點大于8.0×10－4s－1)在兩種方法中對應的典型位置，進行定位誤差分析。選取2010年6月16日23:54和2010年7月10日21:24的PUP合成切變圖與最小二乘法合成切變圖，分別進行對比統計，結果見表3和表4。選取2009年6月5日21:44的PUP合成切變圖與最小二乘法合成切變圖，對比統計陣風鋒典型位置，結果見表5。

通過以上比較，得出如下結論:

圖11 6月5日21:44雷達徑向速度及合成切變圖(半徑:115 km，仰角:1.5°) a．PPI的雷達徑向速度;b．PUP合成切變圖;c．最小二乘法合成切變圖Fig．11 Radar radial velocity and composite shear at 21:44 BST 5 June(radius:115 km;elevation:1.5°) a．radial velocity on PPI;b．composite shear of PUP product;c．composite shear of least square method

1)識別個數。對比分析PUP合成切變圖和最小二乘法的合成切變圖可知，在識別合成切變點個數時，最小二乘法識別得較多，且邊緣切變的識別較突出，合成切變數值在8.0×10－4s－1以上的顯示數量較多，參考價值高，在方位和距離上與速度圖的輻散輻合位置較吻合;PUP產品顯示的合成切變值個數較少且都在大值區，忽略了邊緣切變的識別，參考價值不高。對比表明最小二乘法的識別能力優于PUP產品。

表3 2010年6月16日23:54兩種方法識別結果的輻合輻散定位及誤差分析(P為速度圖中輻合輻散點的位置，P1為PUP產品中對應的位置，P2為最小二乘法對應的位置，D1為PUP產品定位誤差，D2為最小二乘法對應誤差)Table 3 Error analysis and location for convergence and divergence of two kinds results at 23:54 BST 16 June 2010(P is the location of convergence and divergence in radial velocity;P1 is the location of PUP product;P2 is the location of least square method;D1 is the location error of the PUP product;D2 is the error of least square method)

表4 2010年7月10日21:24兩種方法識別結果的輻合輻散定位及誤差分析(P為速度圖中輻合輻散點的位置，P1為PUP產品中對應的位置，P2為最小二乘法對應的位置，D1為PUP產品定位誤差，D2為最小二乘法對應誤差)Table 4 Error analysis and location for convergence and divergence of two kinds results at 21:24 BST 10 July 2010(P is the location of convergence and divergence in radial velocity;P1 is the location of PUP product;P2 is the location of least square method;D1 is the location error of the PUP product;D2 is the error of least square method)

表5 2009年6月5日21:44兩種方法識別結果的輻合輻散定位及誤差分析(P為速度圖中輻合輻散點的位置;P1為PUP產品中對應的位置;P2為最小二乘法對應的位置;D1為PUP產品定位誤差;D2為最小二乘法對應誤差)Table 5 Error analysis and location for convergence and divergence of two kinds results at 21:44 BST 5 June 2009(P is the location of convergence and divergence in radial velocity;P1 is the location of PUP product;P2 is the location of least square method;D1 is the location error of the PUP product;D2 is the error of least square method)

2)定位精度。表3—5表明，最小二乘法的定位與輻合輻散等合成切變大值點的實際位置差別較小，PUP產品定位誤差較大。表3中D1值最大為0.8°、0.9 km，D2 值相應為 0.3°、0.1 km。表4 中，D1 值最大為 1.7°、0.5 km，D2 值相應為 0.1°、0.1 km;D1值最大為2.0°、0.8 km，D2值相應為0.4°、0.1 km。在定位精度上，最小二乘法較準確。

3)邊緣信息。由對比分析看出，PUP合成切變產品的邊緣信息不完整，邊緣點的切變識別較少或直接過濾掉。最小二乘法的合成切變擬合效果較好，邊緣信息提取完整，切變識別個數較多。分析2009年6月5日21:44的PUP合成切變圖和最小二乘法切變圖，可直觀地看出PUP合成切變圖中陣風鋒較零散，出現較明顯的鋸齒狀，而最小二乘法合成切變圖上識別的陣風鋒邊緣清晰。

4)連續性。分析2009年6月5日阜陽資料，最小二乘法擬合特性保證了信息的完整，陣風鋒的位置、平滑度、連續性等方面明顯好于PUP合成切變產品，在分析對比連續時次的徑向速度時，能夠找到速度資料變化的依據，這為預報未來時次的危險切變區域提供了參考，有助于改進機載雷達風切變預警。

4 結論

為探索機載氣象雷達的風切變識別算法，本文用地基多普勒雷達速度資料研究了風切變識別的最小二乘法，通過分析石家莊兩次大風過程及阜陽一次陣風鋒的實例，對比了PUP合成切變產品與最小二乘法合成切變的識別效果，結論如下:

1)在對多普勒雷達速度資料進行識別之前，結合相應時次的天氣圖和T-lgp圖對天氣形勢做出預測非常重要。天氣形勢是產生局地強對流的背景條件，可大致確定對流性天氣現象的位置，再結合雷達速度資料進行識別分析，有助于對風切變危險區域做出預警。

2)最小二乘法擬合效果較好，邊緣信息較完整，徑向速度資料的一維方位切變和一維徑向切變識別與徑向速度圖中的切變、輻合、輻散信息相吻合。其二維合成切變在輻合或輻散等現象出現時，都會產生大值區域，這是對一維切變風場的必要補充。

3)三次個例分析表明，最小二乘法在識別效果、定位精度、邊緣識別等方面均優于PUP合成切變產品。

4)從2009年6月5日阜陽多普勒雷達速度資料合成切變識別圖上可以明顯看出，在陣風鋒的識別方面，最小二乘法具有良好的擬合特性，識別出的陣風鋒在邊界清晰度、平滑性及連續性等方面均具優勢。

5)針對機載雷達資料的寬波束、風速等級少及向下掃描時遠處距庫高度低及易平滑的特點，探索機載雷達風切變算法時擬合窗口的選取尤為關鍵，需綜合考慮機載雷達參數特點，以達到最佳識別和預警效果。

遲繼峰，鐘天宇，劉慶超，等．2012．復雜地形多測風塔綜合地貌及風切變擬合修正的風資源評估方法研究［J］．華電技術，34(11):75-78．

胡明寶，談曙青，湯達章，等．2000．單部多普勒天氣雷達探測低空風切變方法［J］．南京氣象學院學報，23(1):113-119．

黃鈴光．2011．基于雷達數據的強對流天氣的識別算法及實現［D］．南京:南京信息工程大學．

黃儀方．2002．航空氣象［M］．成都:西南交通大學出版社:131-138．

梁建宇，孫建華．2012.2009年6月一次颮線過程災害性大風的形成機制［J］．大氣科學，36(2):316-336．

秦娟，吳仁彪，蘇志剛，等．2011．機載前視氣象雷達地雜波仿真方法［J］．現代雷達，33(8):72-75．

孫一昕，方娟．2012.2010年5月6日重慶強對流過程的天氣學分析［J］．氣象科學，32(6):609-621．

王叢梅，景華，王福俠，等．2011．一次強烈雹暴的多普勒天氣雷達資料分析［J］．氣象科學，31(5):659-665．

王珊珊，黃興友，蘇磊，等．2008．利用雷達回波自動識別鋒線的方法［J］．南京氣象學院學報，31(4):563-573．

魏耀，張興敢．2010．多普勒天氣雷達合成切變算法及改進方法的研究［J］．電子與信息學報，32(1):43-48．

武虎子，耿建中，羅振誼，等．2013．大型飛機在風切變環境中飛行安全尺度分析［J］．飛行力學，31(4):305-308．

袁杲，張志強，謝明元．2007．一種CINRAD雷達基數據讀取丟失的改進方法［J］．成都信息工程學院學報，22(2):178-181．

張京英，孫成武，王慶華，等．2009．一次颮線大風的多種資料分析和臨近預報［J］．氣象科學，29(1):126-132．

張培昌，杜秉玉，戴鐵丕．2001．雷達氣象學［M］．北京:氣象出版社:122-128．

章澄昌．2000．飛行氣象學［M］．北京:氣象出版社:179-187．

Britt C L，Harrah S D，Crittenden L H．1993．Microburst hazard detection performance of the NASA experimental windshear radar system［C］//AIAA93-3943:AIAA Aircraft Design，Systems and Operations Meeting．Monterey，CA:AIAA:1-13．

Chan P W，Shun C M．2010．Latest developments of wind-shear alerting services at the Hong Kong International Airport［C］//14th Conference on Aviation，Range，and Aerospace Meteorology．Atlanta，Georgia:American Meteorological Society:17-21．

Harris F I，Glover K M，Smythe G R．1985．Gust front detection and prediction［C］//Preprints，14th Conference on Severe Local Storms．Boston:Bulletin of the American Meteorological Society:342-345．

Lakshmanan V，Hondl K，Tumpf G，et al．2003．Quality control of WSR-88D data［C］//Preprints，31st International Conference on Radar Meterology．Seattle，WA:American Meteorological Society:522-525．

Williams P．2005．Aircraft trajectory planning for terrain following incorporating actuator constraints［J］．Journal Aircraft，42(5):1358-1362．