多普勒天氣雷達(dá)識(shí)別低空風(fēng)切變算法的改進(jìn)研究

王博

(中國(guó)民用航空西北地區(qū)空中交通管理局寧夏分局，寧夏 銀川 750009)

引言

多普勒天氣雷達(dá)是大氣風(fēng)場(chǎng)中尺度結(jié)構(gòu)的重要探測(cè)工具。在多普勒天氣雷達(dá)的徑向風(fēng)場(chǎng)信息中，大部分中尺度天氣系統(tǒng)都有其獨(dú)特的二維特征結(jié)構(gòu)[1]，可通過(guò)研究這些二維特征結(jié)構(gòu)識(shí)別和追蹤中尺度天氣系統(tǒng)[2-3]。我國(guó)內(nèi)陸目前使用較多的多普勒天氣雷達(dá)是CINRAD-CD，該款雷達(dá)有切變線識(shí)別產(chǎn)品，但它采用的算法并未把缺測(cè)區(qū)納入到分析范疇中，僅對(duì)探測(cè)數(shù)據(jù)進(jìn)行了簡(jiǎn)單的預(yù)處理，且識(shí)別產(chǎn)品僅限于平面位置顯示(plane position indicator，PPI)，沒有涉及到垂直切變。

自20世紀(jì)開始，國(guó)內(nèi)外的一些實(shí)驗(yàn)室、高校均已開展多普勒天氣雷達(dá)風(fēng)切變線識(shí)別領(lǐng)域的研究。CAMPBELL and OLSON[4]成功研發(fā)了可利用人工智能方法自動(dòng)識(shí)別風(fēng)切變的WX1系統(tǒng)。QIN et al.[5]基于多普勒天氣雷達(dá)回波數(shù)據(jù)，利用集合平方根濾波法對(duì)風(fēng)場(chǎng)變化進(jìn)行了預(yù)測(cè)，該方法得到的颮線傳播方向與觀測(cè)值一致，其傳播速度略高于觀測(cè)值，其有效預(yù)報(bào)時(shí)間約為5 h。ZHENG et al.[6]研究了陣風(fēng)鋒發(fā)生過(guò)程中回波的反射率、速度和譜寬等統(tǒng)計(jì)特征，提出了一種基于滑動(dòng)窗口型的陣風(fēng)鋒自動(dòng)識(shí)別算法。胡明寶等[7]利用單部多普勒天氣雷達(dá)探測(cè)低空風(fēng)切變，可計(jì)算出一維徑向風(fēng)切變值、一維方位風(fēng)切變值及二維復(fù)合風(fēng)切變值。魏耀和張興敢[8]基于多普勒天氣雷達(dá)速度場(chǎng)合成切變，首先取不同“擬合窗口”，然后使用最小二乘擬合法進(jìn)行運(yùn)算，該算法在定位精度、識(shí)別能力、邊緣切變識(shí)別等方面優(yōu)勢(shì)明顯。冷亮等[9]、吳舉秀等[10]提出了一種用于陣風(fēng)鋒自動(dòng)識(shí)別的數(shù)學(xué)形態(tài)學(xué)算法，通過(guò)計(jì)算陣風(fēng)鋒的反射率因子差、速度切變、速度梯度等三個(gè)參量場(chǎng)，得到一個(gè)二維特征量，采取數(shù)字形態(tài)學(xué)算法對(duì)該二維特征量進(jìn)行處理，從而獲取陣風(fēng)鋒框架結(jié)構(gòu)特性曲線。國(guó)內(nèi)外相關(guān)研究均集中在識(shí)別風(fēng)切變閾值及識(shí)別算法方面，并未關(guān)注探測(cè)資料的預(yù)處理，且未實(shí)際應(yīng)用這些研究成果。

目前，有些廠家的天氣雷達(dá)已可提供風(fēng)切變產(chǎn)品，這些風(fēng)切變產(chǎn)品包括徑向、仰角、方位角等多個(gè)方向的風(fēng)切變場(chǎng)[11-12]。他們采用的主流算法是：預(yù)處理天氣雷達(dá)的探測(cè)數(shù)據(jù)，并對(duì)缺測(cè)區(qū)域使用線性內(nèi)插技術(shù)填充，再分別計(jì)算風(fēng)場(chǎng)的切變信息[13]，但此前的風(fēng)場(chǎng)結(jié)構(gòu)可能會(huì)由于線性內(nèi)插而發(fā)生變化，從而影響風(fēng)切變的計(jì)算。

本文提出的處理方法可通過(guò)多普勒天氣雷達(dá)的徑向速度有效辨識(shí)風(fēng)切變線，顯示垂直切變圖，可為臨近預(yù)報(bào)提供更科學(xué)的參考信息。

1 計(jì)算方法

CINRAD-CD多普勒天氣雷達(dá)切變線識(shí)別產(chǎn)品的運(yùn)算方法[13]為：先順著徑向選取反射率達(dá)到閾值的像素點(diǎn)，并滑動(dòng)平均原始徑向數(shù)據(jù)，隨后運(yùn)算這部分?jǐn)?shù)據(jù)的徑向切變和方位切變，最后在笛卡爾坐標(biāo)系上對(duì)其進(jìn)行梯度插值且對(duì)其濾波[14]。這種算法的局限之處在于：未將缺測(cè)區(qū)域和回波邊緣的數(shù)據(jù)有效處理，且僅有平面位置顯示(PPI)的產(chǎn)品，沒有垂直切變的產(chǎn)品。

本文使用多普勒天氣雷達(dá)徑向速度的體積掃描數(shù)據(jù)，在極坐標(biāo)上進(jìn)行運(yùn)算。

1.1 探測(cè)數(shù)據(jù)預(yù)處理

為降低徑向速度噪聲，且有效保留徑向速度的中尺度結(jié)構(gòu)，對(duì)探測(cè)資料實(shí)施預(yù)處理。預(yù)處理步驟是：

1)中值濾波

分別沿著方位角和徑向選取n1×n2個(gè)點(diǎn)，對(duì)其排序，但不需選取所有目標(biāo)點(diǎn)的數(shù)值，僅使用中間值代替即可。

2)滑動(dòng)平均

對(duì)探測(cè)數(shù)據(jù)實(shí)施中值濾波后，選取沿著方位角方向和徑向的n3和n4個(gè)連續(xù)像素點(diǎn)，分別計(jì)算其均值作為中間點(diǎn)數(shù)值。針對(duì)缺少探測(cè)數(shù)據(jù)的區(qū)域及附近區(qū)域，對(duì)其采用的處理原則是：唯有滿足一定數(shù)量的有效數(shù)據(jù)才進(jìn)行處理，但是對(duì)其不采取外插值處理；如若有效數(shù)據(jù)不能達(dá)到一定數(shù)量，則把無(wú)數(shù)據(jù)區(qū)周邊的點(diǎn)以無(wú)探測(cè)數(shù)據(jù)處理。如此則可有效規(guī)避邊緣區(qū)域探測(cè)資料不足而導(dǎo)致的誤差，切變線通常是連續(xù)的，但缺測(cè)區(qū)是孤立的，因此去除這部分區(qū)域?qū)τ谧罱K計(jì)算效果影響有限。

參與計(jì)算的探測(cè)數(shù)據(jù)點(diǎn)數(shù)n1、n2、n3、n4不同(即計(jì)算“窗口”大小)，會(huì)極大影響預(yù)處理效果。為降低差別，沿方位角方向的n1、n3通常取值為3。本文重點(diǎn)分析沿著徑向的n2、n4變化時(shí)對(duì)預(yù)處理結(jié)果造成的影響。

1.2 計(jì)算探測(cè)數(shù)據(jù)徑向速度的徑向切變和方位切變

對(duì)任意一組隨機(jī)數(shù)據(jù)(y1,x1)，(y2,x2)，…，(yn,xn)，通過(guò)最小二乘擬合法，可得到利用這些點(diǎn)的最佳直線，且可得出這條直線的斜率。

(1)

式(1)中vi表示第i個(gè)點(diǎn)的徑向速度，ri表示該點(diǎn)距天氣雷達(dá)的距離，n表示參與運(yùn)算的探測(cè)數(shù)據(jù)點(diǎn)數(shù)，θ表示方位角，此處為常數(shù)。

(2)

式(2)中vj表示第j個(gè)點(diǎn)的徑向速度，r在此處為常數(shù)，表示該點(diǎn)距天氣雷達(dá)的距離，n表示參與運(yùn)算的探測(cè)數(shù)據(jù)點(diǎn)數(shù)，θj表示相應(yīng)的方位角。

1.3 計(jì)算垂直切變

(3)

式(3)中v1、v2為上下兩層的徑向速度，并對(duì)其進(jìn)行平滑擬合(對(duì)應(yīng)的區(qū)間范圍均為1)；r1、r2為上下兩層的徑向間距：α1、α2為上下兩層PPI的仰角。僅當(dāng)上下兩層數(shù)據(jù)同時(shí)有效時(shí)才進(jìn)行運(yùn)算，否則作為無(wú)探測(cè)數(shù)據(jù)區(qū)域處理，隨后對(duì)其采取中值濾波。

1.4 計(jì)算組合切變

2 天氣雷達(dá)探測(cè)數(shù)據(jù)

本文利用的雷達(dá)探測(cè)數(shù)據(jù)來(lái)自銀川新一代天氣雷達(dá)(CINRAD-CD)對(duì)2019年6月26日暴雨過(guò)程、2019年8月2日帶狀強(qiáng)降水過(guò)程的探測(cè)。CINRAD-CD雷達(dá)是一種C波段全相參體制的多普勒天氣雷達(dá)，工作頻率為5 300～5 500 MHz，天線直徑4.3 m，增益43 dB，最大可測(cè)半徑為200 km，天線的方位、俯仰定位精度均為0.2°[15]，銀川雷達(dá)饋源海拔高度為1 070.439 m。

3 不同預(yù)處理方法試驗(yàn)對(duì)比

對(duì)2019年6月26日多普勒天氣雷達(dá)的探測(cè)數(shù)據(jù)進(jìn)行預(yù)處理，探測(cè)時(shí)選擇0.6°仰角，45°方位角的徑向，以此次天氣過(guò)程為例說(shuō)明探測(cè)數(shù)據(jù)預(yù)處理的效果。

3.1 中值濾波

首先分別令n2=5，n2=10，n2=20，n2=40，然后對(duì)該探測(cè)數(shù)據(jù)實(shí)施中值濾波，最后把結(jié)果和原始資料比較，將其繪成圖1。

圖1 中值濾波效果對(duì)比

圖1中的原始探測(cè)數(shù)據(jù)存在的庫(kù)間脈動(dòng)比較明顯，經(jīng)過(guò)中值濾波處理之后，庫(kù)間脈動(dòng)基本消失，但是曲線平滑度依然不足。圖1中A—B段、B—C段梯度較大。A—B段，當(dāng)n2=5和n2=10時(shí)，濾波結(jié)果和原始資料較為相近，保持了大梯度值；當(dāng)n2=20和n2=40時(shí)，梯度明顯降低，甚至幾乎消失。B—C段，當(dāng)n2=5，n2=10，n2=20，n2=40的濾波結(jié)果和大梯度區(qū)重合度較高，與原始資料非常接近，可保持大梯度特征。由以上分析可知，對(duì)大尺度的梯度區(qū)域而言，挑選合理的濾波參數(shù)，可在消除庫(kù)間脈動(dòng)的情況下仍然保持探測(cè)數(shù)據(jù)的大梯度特征。

3.2 滑動(dòng)平均

首先分別令n4=5，n4=10，n4=20，n4=40，然后對(duì)該探測(cè)數(shù)據(jù)實(shí)施徑向滑動(dòng)平均，最后把結(jié)果與原始數(shù)據(jù)對(duì)比，將其繪成圖2。

圖2 滑動(dòng)平均效果對(duì)比

由圖2分析可得，當(dāng)n4取值較小時(shí)，可以保持大的梯度區(qū)間，但不能有效去除庫(kù)間脈動(dòng)；當(dāng)n4取值增大時(shí)，庫(kù)間脈動(dòng)的去除效果顯著，但是大梯度區(qū)間逐漸降低。由此可見，若只對(duì)探測(cè)數(shù)據(jù)進(jìn)行滑動(dòng)平均，將無(wú)法滿足去除庫(kù)間脈動(dòng)且保留大梯度區(qū)間的要求。

3.3 中值濾波后滑動(dòng)平均

首先對(duì)原始探測(cè)數(shù)據(jù)實(shí)施中值濾波(n1×n2=3×10)，然后分別令n4=5，n4=10，n4=20，n4=40，對(duì)其作徑向滑動(dòng)平均，最后把結(jié)果與原始數(shù)據(jù)對(duì)比，將其繪成圖3。

圖3 先實(shí)施中值濾波再實(shí)施滑動(dòng)平均的效果對(duì)比

由圖3分析可得，當(dāng)n4=5，n4=10，n4=20(即n4較小時(shí))，3個(gè)處理結(jié)果較為接近，不僅過(guò)濾了庫(kù)間脈動(dòng)，而且保留了原有的中尺度結(jié)構(gòu)。但是當(dāng)n4=40時(shí)，該處理方法會(huì)影響中尺度結(jié)構(gòu)。

3.4 切變計(jì)算

對(duì)探測(cè)數(shù)據(jù)實(shí)施切變計(jì)算的過(guò)程中，探測(cè)數(shù)據(jù)點(diǎn)數(shù)的不同表示辨識(shí)風(fēng)切變的尺度不同。通過(guò)對(duì)比預(yù)處理結(jié)果發(fā)現(xiàn)：在預(yù)處理尺度和切變計(jì)算尺度相差不大(即都較大或較小)時(shí)，風(fēng)切變的辨識(shí)效果不良；預(yù)處理尺度和切變計(jì)算尺度相差較大時(shí)可得到良好的辨識(shí)效果。

4 實(shí)例識(shí)別結(jié)果

在對(duì)銀川市多普勒天氣雷達(dá)的探測(cè)數(shù)據(jù)處理過(guò)程中，不僅要有效過(guò)濾庫(kù)間脈動(dòng)，盡可能保留天氣過(guò)程原有的中尺度結(jié)構(gòu)，而且要盡可能多地使用探測(cè)數(shù)據(jù)，因此選用n1=3，n2=10，n3=3，n4=20。

4.1 暴雨天氣過(guò)程識(shí)別結(jié)果

2019年6月26日銀川市發(fā)生一次暴雨天氣過(guò)程，與中尺度切變線有關(guān)(圖4)。分別對(duì)比同一仰角上的反射率因子和徑向速度，繼續(xù)對(duì)比同時(shí)次、不同仰角間的垂直切變。分析可得：此天氣過(guò)程強(qiáng)降水區(qū)域的反射率因子為25～40 dBZ。負(fù)垂直切變和正垂直切變分別出現(xiàn)在切變線的北側(cè)和南側(cè)，該現(xiàn)象說(shuō)明垂直切變的計(jì)算結(jié)果與實(shí)際徑向風(fēng)場(chǎng)的變化是一致的。

圖4 銀川市2019年6月26日天氣過(guò)程識(shí)別結(jié)果(a. 0.6°仰角反射率因子；b. 0.6°仰角徑向速度；c. 0.6°仰角、1.4°仰角間垂直切變；d. 0.6°仰角徑向切變，顯示閾值：0.5 m·s-1·km-1；e. 0.6°仰角方位切變，顯示閾值：0.8 m·s-1·km-1；f. 0.6°仰角組合切變，顯示閾值：0.8 m·s-1·km-1)

4.2 帶狀強(qiáng)降水天氣過(guò)程識(shí)別結(jié)果

2019年8月2日銀川市發(fā)生一次帶狀強(qiáng)降水天氣過(guò)程(圖5)。分析可見，圖5a中強(qiáng)降水區(qū)的反射率因子為40～50 dBZ，正垂直切變最大值出現(xiàn)在降水區(qū)前部，負(fù)垂直切變出現(xiàn)在降水區(qū)后部。帶狀回波強(qiáng)度略有減弱，切變線識(shí)別結(jié)果強(qiáng)度減弱且不連續(xù)。

圖5 銀川市2019年8月2日天氣過(guò)程識(shí)別結(jié)果(a. 0.6°仰角反射率因子；b. 0.6°仰角徑向速度；c. 0.6°仰角、1.4°仰角間垂直切變；d. 0.6°仰角徑向切變，顯示閾值：1.2 m·s-1·km-1；e. 0.6°仰角方位切變，顯示閾值：1.7 m·s-1·km-1；f. 0.6°仰角組合切變，顯示閾值：1.7 m·s-1·km-1)

5 試驗(yàn)結(jié)果分析

本文通過(guò)對(duì)多普勒天氣雷達(dá)徑向速度的探測(cè)資料進(jìn)行預(yù)處理，得出垂直切變圖，探究了識(shí)別風(fēng)切變線的新方法，并將該方法分別應(yīng)用于暴雨和帶狀強(qiáng)降水兩個(gè)案例，證明此方法可有效識(shí)別風(fēng)切變。識(shí)別結(jié)果表明：

1)首先對(duì)天氣雷達(dá)的探測(cè)數(shù)據(jù)實(shí)施預(yù)處理，然后挑選合理的運(yùn)算“窗口”，通過(guò)最小二乘擬合法得到徑向速度沿著兩個(gè)方向的切變現(xiàn)象，辨識(shí)結(jié)果與徑向速度在圖中的切變信息保持一致，證明了該方法的有效性。

2)由于各降水過(guò)程中顯示閾值存在著區(qū)別，垂直切變能夠提供徑向風(fēng)場(chǎng)的高低層配置信息，徑向速度水平切變的強(qiáng)弱和降水強(qiáng)度具有一定聯(lián)系。

3)通過(guò)合理運(yùn)用多普勒天氣雷達(dá)探測(cè)到的徑向速度信息，提供垂直切變圖，能夠有效識(shí)別風(fēng)切變，為災(zāi)害氣象預(yù)報(bào)提供更科學(xué)的參考信息。

然而，需要特別指出的是，本文提供的方法可有效辨識(shí)徑向速度的切變，但與真實(shí)風(fēng)場(chǎng)的切變尚有區(qū)別。對(duì)于算法中顯示閾值的設(shè)置，需根據(jù)各個(gè)地區(qū)的現(xiàn)實(shí)狀況，對(duì)一系列天氣案例采取試驗(yàn)、分類和歸納，方可獲得最佳結(jié)果，不能將本文的閾值直接使用。