簡化VVP算法在一次雷暴風(fēng)場反演中的應(yīng)用

吳福浪

(中國民用航空寧波空中交通管理站，浙江 寧波 315000)

簡化VVP算法在一次雷暴風(fēng)場反演中的應(yīng)用

吳福浪

(中國民用航空寧波空中交通管理站，浙江 寧波 315000)

在研究Waldteufel等提出的VVP算法的基礎(chǔ)上，基于風(fēng)場在劃定的小積分體積內(nèi)是呈均勻分布的假設(shè)，提出簡化的VVP算法，并用簡化的VVP算法反演2015年04月02日發(fā)生在寧波機(jī)場的一次強(qiáng)雷暴的風(fēng)場；理論和反演結(jié)果表明，簡化的VVP算法能較好的反演出此次強(qiáng)雷暴的風(fēng)場分布。

簡化VVP算法；強(qiáng)雷暴；寧波機(jī)場；風(fēng)場反演

0 引 言

雷暴是一種災(zāi)害性天氣，強(qiáng)烈雷暴的發(fā)生，常伴隨大風(fēng)、大雨或者冰雹等強(qiáng)對流天氣出現(xiàn)，這種天氣狀況的出現(xiàn)對于飛機(jī)飛行來說，其影響是非常嚴(yán)重的。新一代多普勒天氣雷達(dá)具有高時空分辨率的特點(diǎn)，同時能夠獲取粒子的徑向速度和反射率因子等要素，成為研究災(zāi)害性天氣和超短期預(yù)報的重要手段。但是實際工作中，多普勒天氣雷達(dá)測得的三維風(fēng)場是徑向的，并不是實際的三維風(fēng)場矢量，這給準(zhǔn)確監(jiān)測以及預(yù)報雷暴天氣系統(tǒng)造成很大誤差。因此,多普勒雷達(dá)風(fēng)場反演受到普遍的關(guān)注。從20世紀(jì)60年代至今，直接從單部多普勒天氣雷達(dá)估計矢量風(fēng)場已經(jīng)取得相當(dāng)多的研究成果。

在國內(nèi)方面，忻翎艷等(1990)率先討論了利用VAD、VARD、VVP等技術(shù)獲取風(fēng)場信息的方法。陶祖鈺(1992)在假定相鄰方位角風(fēng)矢量相等的前提下，提出了一種反演風(fēng)場信息的速度方位處理(VAP)的新技術(shù)。姜海燕等(1997)從動力學(xué)角度出發(fā)，應(yīng)用中尺度渦度方程，經(jīng)過簡化導(dǎo)出風(fēng)速的徑向分量、切向分量和渦度、散度的關(guān)系。梁旭東等(2007)提出采用“按方位角均勻假定反演關(guān)系”進(jìn)行二維風(fēng)場反演的積分VAP技術(shù)(IVAP)，該技術(shù)與傳統(tǒng)VAP技術(shù)不同的是在整個選定的反演區(qū)間上進(jìn)行二維風(fēng)場反演，而不是僅考慮反演區(qū)間的兩個端點(diǎn)，因此IVAP技術(shù)具有濾波性能，而且可以根據(jù)需要調(diào)整反演區(qū)間以得到滿足一定平滑性并能反映中小尺度信息的風(fēng)場信息。周生輝等(2014)利用VVP算法，選取量級最大的3個主要參量進(jìn)行反演，引入隨機(jī)的觀測誤差，通過改變模擬風(fēng)速確定了反演算法的適用范圍。

國外方面，Lhermitte等(1961)提出了VAD方法，從圓周上反演風(fēng)場；VAD技術(shù)通常用于計算水平平均風(fēng)速的垂直廓線。在VAD假設(shè)的基礎(chǔ)上，Browning(1968)提出了線性風(fēng)場假設(shè)的反演方法，通過傅里葉變換將風(fēng)場的風(fēng)速、散度等參量分離了出來。Easterbrook等(1975)在局地均勻風(fēng)的假定條件下，提出了速度面積顯示(VARD)方法。Waldteufel等(1979)通過徑向風(fēng)場在空間呈線性分布的假設(shè)，提出VVP算法。VVP算法是將線性模式應(yīng)用到一個小體積內(nèi)的單多普勒徑向速度,估計出局地的散度和形變。Srivastava等(1986)根據(jù)距離圈的位置不同對傅里葉系數(shù)引入了權(quán)重系數(shù)來進(jìn)行訂正，并根據(jù)提供的垂直速度和散度的關(guān)系式，得到了垂直速度的廓線。Caya等(1992)在對線性風(fēng)場假設(shè)的討論中，指出通過傅里葉變換得到的各系數(shù)的物理含義并不準(zhǔn)確，這種數(shù)學(xué)上的處理得不到真實的風(fēng)場信息，特別是風(fēng)場特征為非線性時，為此提出了利用設(shè)定閾值和利用多項式擬合方法來得到近似解，并取得了較為可信的結(jié)果。

雷達(dá)反演風(fēng)場在強(qiáng)對流應(yīng)用方面，梅玨等(2006)通過采用“給定方位角內(nèi)風(fēng)場均勻假定”提出積分速度—方位處理技術(shù)(IVAP)，并利用IVAP對2001年08月05日上海強(qiáng)雷雨時虹橋機(jī)場Doppler雷達(dá)速度資料進(jìn)行中尺度風(fēng)場反演試驗，反演風(fēng)的風(fēng)向風(fēng)速和地面實測風(fēng)的風(fēng)向風(fēng)速有很好的相關(guān)性。周鑫等(2009)提出利用交叉相關(guān)法反演高空風(fēng)場，并分析2006年07月06日發(fā)生在四川綿陽的一次強(qiáng)對流天氣，反演的風(fēng)場結(jié)構(gòu)和速度圖的關(guān)系是對應(yīng)的。夏網(wǎng)萍等(2010)針對單多普勒雷達(dá)風(fēng)場反演VVP方法中存在的線性方程組的病態(tài)矩陣問題，使用改進(jìn)的SVVP方法，通過多種模擬風(fēng)場信號對其進(jìn)行檢驗和分析，并給出兩個反演實例：2009年江蘇雨雪天氣和安徽冰雹個例。模擬和實例表明SVVP方法反演的風(fēng)場結(jié)構(gòu)合理，對于反演強(qiáng)對流天氣比較有效。周小剛等(2015)由單多普勒雷達(dá)反演水平風(fēng)場的基本原理人手，回顧了VAD產(chǎn)品在在暴雨與強(qiáng)對流臨近預(yù)報中的應(yīng)用，并重點(diǎn)對VAD產(chǎn)品的業(yè)務(wù)應(yīng)用誤區(qū)進(jìn)行了討論。

利用多普勒天氣雷達(dá)的徑向風(fēng)場信息反演的實際三維風(fēng)場有助于更加準(zhǔn)確研究雷暴天氣系統(tǒng)的風(fēng)場結(jié)構(gòu)與演變特征。三維風(fēng)場反演以Waldteufe等(1979)提出的VVP算法為代表。但是VVP算法在求解矩陣過程中造成較大擾動，引起大的舍入誤差，導(dǎo)致反演的風(fēng)速誤差很大；并且計算量過程復(fù)雜，耗時相對較多限制了它在實際工作中的應(yīng)用。本文意在通過對VVP算法進(jìn)行簡化，實現(xiàn)它在實際業(yè)務(wù)應(yīng)用中對雷暴天氣三維風(fēng)場的準(zhǔn)確反演，從而充分利用多普勒雷達(dá)的風(fēng)場信息，為促進(jìn)民航預(yù)報水平和多普勒雷達(dá)反演風(fēng)場在民航科研和業(yè)務(wù)中的應(yīng)用提供可靠的理論基礎(chǔ)。并且對于雷暴風(fēng)場的分析與探究有助于給飛機(jī)飛行帶來保護(hù)，其意義是重大的；充分的了解雷暴風(fēng)場，良好的掌握處理應(yīng)急狀況的方法，確保飛機(jī)在遇到狀況時能夠安全的飛行。

1 簡化的VVP算法

Waldteufel和Corbin等于1979年提出了VVP算法(Volume Velocity Processing)，通過徑向風(fēng)場在空間呈線性分布的假設(shè)，理論上可以獲得三維風(fēng)場的12個變量，求解出它的各項動力因子，比如風(fēng)向，風(fēng)速，輻散場等。VVP算法是基于風(fēng)場線性分布假設(shè)的，在仰角層數(shù)足夠多的觀測資料中取一個多仰角多方位角多距離庫長的三維小體積做積分，并在分析體積中假設(shè)：1)風(fēng)場在三維小體積內(nèi)呈線性變化；2)掃描期間，風(fēng)場不隨時間變化。

由于VVP算法分析體積較大，各個元素在量級上存在較大差異，在編程求解矩陣過程中造成較大擾動，引起大的舍入誤差，導(dǎo)致矩陣的嚴(yán)重病態(tài)，徑向速度經(jīng)過病態(tài)系數(shù)矩陣的放大，導(dǎo)致反演的風(fēng)速誤差很大(馮果忱等,1991)；并且實際工作中電腦性能一般，因此計算速度較慢，無法及時的把三維風(fēng)場反演用于實際預(yù)報工作。

針對VVP算法求解中遇到的困難，一般的做法是減少反演參數(shù)的數(shù)量來降低計算困難(Xin等，1998)。魏鳴等(1998)通過對系數(shù)矩陣的處理降低了條件數(shù)，并用共軛梯度法減少求解難度。Li等(2007)利用分布計算的方法，先求出均勻風(fēng)場假設(shè)下的反演結(jié)果，再通過風(fēng)場參量之間的聯(lián)系獲得其他的量。

根據(jù)上述風(fēng)場反演中的問題和研究進(jìn)展，并結(jié)合多年預(yù)報工作總結(jié)發(fā)現(xiàn)，三維風(fēng)場在一個較大的局部體積內(nèi)線性分布的可能性較??；如果把局部積分體積劃分盡可能的大，那么計算量過程變得復(fù)雜，耗時相對較多限制了使用。實際上在一個很小的積分體積中，三維風(fēng)場分布更接近于均勻分布(余艷梅等，2008)。Holleman(2005)通過和探空資料的對比，發(fā)現(xiàn)VVP算法在省略部分待求參量的情況下，仍可以得到較好的反演結(jié)果，基于此，我們?yōu)榱藴p少計算量，提高三維風(fēng)場反演精度，提出簡化的VVP算法。

簡化的VVP算法基于風(fēng)場在劃定的小積分體積內(nèi)是呈均勻分布的假設(shè)，即在劃分的三維小積分體積內(nèi)，各個點(diǎn)的三維風(fēng)場矢量是相同的，并且在多普勒天氣雷達(dá)一個體掃時間內(nèi)，風(fēng)場不隨時間變化。

根據(jù)以上假設(shè)，假如積分體積內(nèi)有m個數(shù)據(jù)點(diǎn)，那么有：

V1=V2=V3=…=Vm

(1)

每一個點(diǎn)的三維風(fēng)場矢量Vp在徑向上的速度分量Vrp為：

Vrp=usinθpcosφp+νcosθpcosφp+wsinφp

(2)

其中φp為雷達(dá)仰角，θp為方位角，rp為格點(diǎn)距離雷達(dá)中心的位置。從方程(2)可以看出，理論上只要有3個點(diǎn)，就能直接解出風(fēng)場矢量Vp。但是實際中并不是所有的點(diǎn)都滿足均勻假設(shè)，因此造成方程(2)無定解。為獲得與雷達(dá)資料同分辨率的三維風(fēng)場，本文采用最小二乘法方法來估計三維風(fēng)場矢量。

根據(jù)最小二乘法原理，有：

(3)

其中N是小積分體積的數(shù)據(jù)點(diǎn)個數(shù)，為了使方程(3)方差和最小，分別對u、ν、w求偏導(dǎo)數(shù)，并使得偏導(dǎo)數(shù)等于零，則有：

(4)

解方程(4)可以得到3個速度分量u、ν、w。若要反演出三維風(fēng)場風(fēng)矢量，要求解u、ν、w3個變量，則至少需要3個格點(diǎn)。這要求小積分體積必須大于等于2x2x2。經(jīng)過實際預(yù)報工作多次測試，本文把分析體積定為4x4x4。

2 雷暴三維風(fēng)場反演個例

在均勻假設(shè)風(fēng)場模型中，舍棄較小量級的參量之后，簡化的VVP算法反演結(jié)果會好于選擇基于線性風(fēng)場假設(shè)的復(fù)雜風(fēng)場模型時的結(jié)果(周生輝等，2014)，因此，我們根據(jù)簡化的VVP算法來反演2015年4月2日(世界時，下同)發(fā)生在寧波機(jī)場一次強(qiáng)雷暴過程的三維風(fēng)場。

強(qiáng)雷雨發(fā)生時段為2015年04月02日17:30

UTC至19:00 UTC。在強(qiáng)雷雨發(fā)生前，我國北方有弱冷空氣滲透南下，15:00 UTC地面圖上(圖1a)氣旋主體位于朝鮮半島，冷鋒由東北沿華東中部向西南位置延伸至我國華南地區(qū)。寧波市剛好處在地面冷鋒附近；12:00 UTC的各層高空圖上(圖略)，寧波市處于高空槽前西南氣流之中，中低層明顯的西南急流帶來充足的暖平流和水汽輸送。隨著高空槽系統(tǒng)東移南壓，高空槽逼近寧波市，高低層配置觸發(fā)強(qiáng)雷暴的發(fā)生發(fā)展。從雷達(dá)資料來看，從杭州灣到浙西南有東北—西南向?qū)拸V的雨帶分布，雨帶走向與地面鋒面走向一致，在寬廣的降水帶中嵌有較強(qiáng)的對流性降水回波，其中南北向的細(xì)長帶狀回波為影響本場的對流云系(圖略)。前期，該對流云帶自西向東快速向?qū)幉C(jī)場靠近，移速約50 km/h。在對流云帶距寧波機(jī)場10 km左右時，寧波機(jī)場出現(xiàn)雷雨天氣。對流云帶影響寧波機(jī)場時出現(xiàn)強(qiáng)降水及大風(fēng)天氣，最大風(fēng)速5～6級，陣風(fēng)8～9級，1 h累計降水量超過30 mm；強(qiáng)降水及大風(fēng)出現(xiàn)在18:06 UTC左右，對比此時雷達(dá)反射率(圖1b)，對流云帶此時正好位于寧波機(jī)場上空。

本文在Visual Studio 2008的C++環(huán)境下編程實現(xiàn)對簡化VVP算法的計算，并對單多普勒雷達(dá)基數(shù)據(jù)進(jìn)行了反演。在對此次雷暴三維風(fēng)場進(jìn)行反演前，先利用線性插值算法和補(bǔ)零算法填補(bǔ)風(fēng)場缺測數(shù)據(jù)，再用9點(diǎn)中值濾波濾去較大噪音點(diǎn)。速度退模糊采用“二維多路多普勒雷達(dá)風(fēng)場自動退模糊算法”(蔡親波等，2000)進(jìn)行退模糊，最后再根據(jù)我們提出的簡化VVP算法反演三維風(fēng)場矢量并繪圖。

圖1 雷暴發(fā)生前天氣形勢和發(fā)生時雷達(dá)反射率回波

為了說明反演風(fēng)場資料風(fēng)場的可靠性，本文選擇雷雨影響寧波機(jī)場前17:32 UTC和影響最強(qiáng)時刻18:06 UTC 1.5 km～3 km高度的反演風(fēng)場進(jìn)行檢驗：在1.5 km高度上，距離雷達(dá)中心50 km范圍圈內(nèi)，簡化VVP算法反演的風(fēng)場與雷達(dá)探測徑向風(fēng)速分布有較好的對應(yīng)關(guān)系(圖略)；在3.0 km高度上，距離雷達(dá)中心100 km范圍圈內(nèi)，簡化VVP算法反演的風(fēng)場與雷達(dá)探測徑向風(fēng)速分布有較好的對應(yīng)關(guān)系(圖略)；綜上可知，簡化VVP算法反演的風(fēng)場資料是可靠的。

在2015年04月02日17:32 UTC的1.5 km反射率回波上(圖略)，強(qiáng)對流區(qū)分為兩個區(qū)域：位于A位置西側(cè)約40 km(記為TS1)和位于B位置西側(cè)約45 km(記為TS2)；對應(yīng)的徑向速度圖和反演風(fēng)場(圖略)，(1.5 km高度)以及圖(圖略)，(3 km高度)。由圖(圖略)可知，1.5 km反演風(fēng)場有明顯的風(fēng)切變輻合區(qū)，分別位于A位置西北15 km和B位置西南10 km。對比與TS1和TS2區(qū)域，反演風(fēng)場得出的風(fēng)切變輻合區(qū)超前于強(qiáng)對流區(qū)30 km左右。此時寧波櫟社機(jī)場實況為東北風(fēng)，風(fēng)速2～3 m/s。對比同時刻3 km高度反演風(fēng)場(圖略)，反演的風(fēng)切變輻合區(qū)位置與1.5 km高度反演位置相差不大，但是長度比1.5 km高度的長。

18:06 UTC的1.5 km反射率回波上(圖1b)，強(qiáng)對流區(qū)TS1移動至寧波機(jī)場上空，位于A位置西側(cè)2 km左右位置，強(qiáng)度加強(qiáng)；TS2移動至B區(qū)域，強(qiáng)度略有減弱。同時刻1.5 km徑向速度圖上，A位置有一個明顯的中氣旋，與圖1b的強(qiáng)雷暴位置相對應(yīng)；B位置附近，50～100 km為負(fù)速度，雷達(dá)中心至50 km范圍為正速度；這是明顯的速度輻合速度圖，對應(yīng)圖1b雷達(dá)中心北側(cè)反射率回波30～35dBz區(qū)域。同時刻1.5 km反演的風(fēng)場如圖2d，此時反演的風(fēng)切變輻合區(qū)超前強(qiáng)對流區(qū)域10～20 km。此時寧波機(jī)場實況出現(xiàn)大雷雨天氣，西北風(fēng)，平均風(fēng)5～6級，陣風(fēng)8～9級，給飛機(jī)安全飛行帶來極大影響。對比同時刻3 km高度反演風(fēng)場(圖略)，反演出的風(fēng)切變輻合區(qū)位置更加遠(yuǎn)離強(qiáng)雷暴區(qū)域，并且長度比1.5 km高度的長。

根據(jù)1.5 km高度反演風(fēng)場不同時刻風(fēng)切變輻合區(qū)演變特征可知，17:32 UTC至18:06 UTC之間，B位置附近風(fēng)切變輻合區(qū)移動緩慢，速度大約為20 km/h；而A位置由于雷暴較強(qiáng)，移動速度也快，17:32 UTC至18:06 UTC之間A位置附近風(fēng)切變輻合區(qū)移動速度大約為40 km/h。在3 km高度反演風(fēng)場不同時刻風(fēng)切變輻合區(qū)演變特征上，由于環(huán)境風(fēng)的引導(dǎo)，17:32 UTC至18:06 UTC之間，風(fēng)切變輻合區(qū)移動了25 km，速度約50 km/h。

3 結(jié) 語

本文對VVP算法以及簡化VVP算法進(jìn)行了介紹分析，并利用簡化VVP算法對2015年04月02日發(fā)生在寧波機(jī)場的一次強(qiáng)雷暴過程進(jìn)行風(fēng)場反演，主要結(jié)論有：

1)簡化VVP算法能較好的反演此次強(qiáng)雷暴過程風(fēng)場。在低層1.5 km高度，距離雷達(dá)中心50 km范圍圈內(nèi)反演效果較好；在3 km高度，距離雷達(dá)中心100 km范圍圈內(nèi)反演效果較好。

2)1.5 km高度上，反演風(fēng)場的風(fēng)切變輻合區(qū)超前于強(qiáng)對流系統(tǒng)10～30 km。風(fēng)切變輻合區(qū)未移動到機(jī)場時，寧波機(jī)場風(fēng)速較小；移動到機(jī)場時，機(jī)場出現(xiàn)大風(fēng)天氣，對飛機(jī)安全起飛降落有很大影響。

3)1.5 km高度，如果對應(yīng)強(qiáng)對流強(qiáng)度較強(qiáng)，風(fēng)切變輻合區(qū)移動速度也相應(yīng)較快，大約為40 km/h。在3 km高度，由于環(huán)境風(fēng)引導(dǎo)氣流的疊加，速度較快約50 km/h。

4)飛機(jī)在飛行過程中，不僅應(yīng)該避開強(qiáng)雷暴區(qū)域，也應(yīng)避免在風(fēng)切變輻合區(qū)中飛行。簡化VVP算法為促進(jìn)民航預(yù)報水平和多普勒雷達(dá)反演風(fēng)場在民航科研和業(yè)務(wù)中的應(yīng)用提供可靠的理論基礎(chǔ)。

[1] 蔡親波,顏瓊丹.二維多路多普勒雷達(dá)風(fēng)場自動退模糊算法應(yīng)用研究[J].氣象科學(xué),2000,29(5):625-632.

[2] 馮果忱,于庚蒲.矩陣分析導(dǎo)論[M].長春:吉林大學(xué)出版社,1991.

[3] 姜海燕,葛潤生.一種新的單部多普勒天氣雷達(dá)反演技術(shù).應(yīng)用氣象學(xué)報,1997,8(2):219-223.

[4] 梁旭東,王斌.多普勒雷達(dá)積分VAP反演技術(shù).氣象學(xué)報,2007,65(2):261-271.

[5] 梅玨,梁旭東,呂環(huán)宇.單Doppler雷達(dá)IVAP風(fēng)場反演試驗.北京大學(xué)學(xué)報(自然科學(xué)版),2006,3(1).

[6] 陶祖鈺.從單Doppler速度場反演風(fēng)矢量的VAP方法.氣象學(xué)報,1992,50(1):81-90.

[7] 魏鳴.1998.單多普勒天氣雷達(dá)資料的變分同化三維風(fēng)場反演和VVP三維風(fēng)場反演.南京:南京大學(xué):97-105.

[8] 忻翎艷,湯達(dá)章.由VVP探測技術(shù)提取風(fēng)場信息的方法.氣象科學(xué),1990,10(4):401-408.

[9] 夏網(wǎng)萍,魏鳴,李南.單多普勒雷達(dá)SVVP風(fēng)場反演及應(yīng)用.氣象科學(xué),2010,30(4):503-509.

[10] 余艷梅,馮啟航,謝明,等.VVP三維反演方法研究.四川大學(xué)學(xué)報,2008,45(3):321-326.

[11] 周鑫,何建新,劉艷,等.多普勒雷達(dá)高空風(fēng)場反演方法研究.現(xiàn)代雷達(dá),2009,5(31):18-20.

[12] 周生輝,魏鳴,張培昌,等.簡化VVP反演算法在臺風(fēng)風(fēng)場反演中的應(yīng)用.遙感學(xué)報,2014,18(5):1128-1137.

[13] 周小剛,費(fèi)海燕,王秀明.基于多普勒雷達(dá)算法的業(yè)務(wù)應(yīng)用討論氣象.2015,41(1):113-120.

[14] Browning K A, Wexler R. The determination of kinematic properties of a wind field using Doppler radar[J]. Journal of Applied Meteorology,1968,7(1):105-113.

[15] Caya A, Laroches S, Zawadzki I, et al. Using single-Doppler data to obatain a mesoscale environmental filed[J]. Journal of Atmospheric & Oceanic Technology,2002,19(1):21-36.

[16] Esterbrook C C. Estimating Horizontal Wind Fields by Two-Dimensional Curve Fitting of Single-Doppler Radar Measurements. Preprints 16th Radar Meteorology Conf., Houston: Amer. Meteor. Sci.,1975.214-219.

[17] Holleman I. 2005.Quality control and verification of weather radar wind profiles.Journal of Atmospheric and Oceanic Technology,22(10):1541-1550.

[18] Lhermitte R M, Atlas D. Precipitation motion by pulse doppler[C]. Boston: Ninth Weather Radar Conference Amer Meteor Soc,1961.

[19] Li N, Wei M, Tang X W and Pan Y J. 2007.An improved velocity volume processing method.Advances in Atmospheric Sciences,24(5):893-906.

[20] Srivastava R C, Matejka T J, Lorello T J. Doppler radar study of the trailing anvile region associated with a squall line[J]. Journal of the atmospheric sciences,1986,43(4):356-377.

[21] Waldteufel Pand Corbin H. On the analysis of single-Doppler radar data. J. Appl. Meteor., 1979,18:532-542.

[22] Xin L G, Reuter W.1998.VVP technique applied to an Alberta storm.Journal of Atmospheric and Oceanic Technology,15(2):587-592.

2016-12-07