本發(fā)明涉及氣象學(xué)領(lǐng)域，特別涉及一種基于多普勒天氣雷達資料的陣風(fēng)鋒自動識別方法。

背景技術(shù)：

天氣雷達是對強對流天氣進行監(jiān)測和預(yù)警的主要工具之一。天氣雷達發(fā)射脈沖形式的電磁波，當(dāng)電磁波遇到降水物質(zhì)(雪花、雨滴和冰雹等)，大部分能量繼續(xù)前進，而有一小部分能量被降水物質(zhì)向四面八方散射，向后散射的能量被雷達接收。

在對流單體的成熟階段，冷性下沉氣流作為一股冷空氣，在近地面的底層向外擴展，與單體運動前方的暖濕氣流交匯形成鋒面。這種現(xiàn)象反映在雷達反射率圖上，通常可以觀測到一條若隱若現(xiàn)、粗細不等、強度較弱且取值不定的弱脊帶狀區(qū)域，氣象上稱之為陣風(fēng)鋒或雷暴的出流邊界。陣風(fēng)鋒常引起氣壓突變、風(fēng)速風(fēng)向劇烈變化以及氣溫驟降等強烈天氣現(xiàn)象，會導(dǎo)農(nóng)作物倒伏、大樹或樹枝折斷、威脅飛機的起降，嚴(yán)重影響人民大眾的生命財產(chǎn)安全。

從圖像角度看，做陣風(fēng)鋒區(qū)域的橫剖面，其反射率值的分布呈脊?fàn)睿捎趨^(qū)域自身反射率值普遍較低，致使由區(qū)域兩側(cè)過渡到區(qū)域外圍的反射率值的梯度變化微弱，且時常出現(xiàn)部分區(qū)域段混于大面積的弱回波區(qū)域之中的現(xiàn)象，致使原本就不太強的梯度變化在這部分區(qū)域段基本消失。這就使常規(guī)的基于邊緣和區(qū)域圖像分割方法用于檢測這種陣風(fēng)鋒區(qū)域時奏效甚微,另外，弱回波區(qū)域中通常含有脊?fàn)罘顷囷L(fēng)鋒帶的干擾，這些干擾甚至與真正的陣風(fēng)鋒相鄰或相交。使識別陣風(fēng)鋒方法的空識率與擊中率之間的矛盾愈發(fā)突出。

在基于多普勒雷達數(shù)據(jù)的陣風(fēng)鋒識別算法中，Delanoy和Troxel^[1]較早提出了基于模糊理論的函數(shù)模板相關(guān)法,該方法將模板尺度固定為17×7，為使檢測算法與陣風(fēng)鋒的方向無關(guān)，該模板需圍繞其中心多次旋轉(zhuǎn)。Osama Alkhouli和Victor DeBrunner^[2]應(yīng)用熵函數(shù)模板法自動識別邊界層輻合線，且方法無需旋轉(zhuǎn)模板。鄭佳峰等^[3]結(jié)合速度場與強度場數(shù)據(jù)檢測陣風(fēng)鋒。在速度場數(shù)據(jù)中，檢測輻合帶；在反射率強度場中，該方法主要使用的是雙向梯度法尋找弱窄帶回波，兩者結(jié)合確定陣風(fēng)鋒。

發(fā)明人在實現(xiàn)本發(fā)明的過程中，發(fā)現(xiàn)現(xiàn)有技術(shù)中至少存在以下缺點和不足：

首先，人為觀測費時費力，會影響到預(yù)報的時效性。文獻[1]只適合檢測寬度不大于3與周邊區(qū)域存在較大梯度的陣風(fēng)鋒。文獻[2]應(yīng)用熵函數(shù)模板法提高了檢測速度，但僅對典型的陣風(fēng)鋒有效，對于部分陣風(fēng)鋒區(qū)域融入背景的情況，會使檢測結(jié)果出現(xiàn)斷裂。文獻^[3]的方法在速度場數(shù)據(jù)中出現(xiàn)無效速度的區(qū)域時無法發(fā)揮；在反射率強度場時，雙向梯度法限定窄帶回波與母體回波距離。至今為止，在對陣風(fēng)鋒檢測的文獻中尚未見到陣風(fēng)鋒與其他弱窄帶回波碰撞產(chǎn)生交疊現(xiàn)象的解決方案，也未見到在陣風(fēng)鋒的檢測中出現(xiàn)斷裂現(xiàn)象的處理方法。

[參考文獻]

[1]Delanoy R L,Troxel S W.A machine intelligent gust front algorithm for doppler weather radars[C].Contributions to the American Meteorological Society’s 26th International Conference on Radar Meteorology.1993:9.

[2]Alkhouli O,DeBrunner V.Gust front detection in weather radar images by entropy matched functional template[C].Image Processing,2005.ICIP 2005.IEEE International Conference on.IEEE,2005,1:I-645-8.

[3]鄭佳鋒,張杰,朱克云,等.陣風(fēng)鋒自動識別與預(yù)警[J].應(yīng)用氣象學(xué)報,2013,24(1):117.

[4]Ojala T,M,Harwood D.A comparative study of texture measures with classification based on featured distributions[J].Pattern recognition,1996,29(1):51-59.

[5]G.Two-frame motion estimation based on polynomial expansion[M].Image Analysis.Springer Berlin Heidelberg,2003:363-370.

技術(shù)實現(xiàn)要素：

本發(fā)明公開了一種基于多普勒天氣雷達資料的陣風(fēng)鋒自動識別方法，可以解決的技術(shù)問題包括：自動識別不同寬度的陣風(fēng)鋒；能識別與周邊區(qū)域的梯度較小或部分域融入背景的陣風(fēng)鋒；能分離與其他弱窄帶回波碰撞產(chǎn)生交疊的陣風(fēng)鋒；對陣風(fēng)鋒檢測中出現(xiàn)斷裂現(xiàn)象進行有效連接的處理；實現(xiàn)準(zhǔn)確、完整自動識別陣風(fēng)鋒的目的。

為了解決上述技術(shù)問題，本發(fā)明提出的一種基于多普勒天氣雷達資料的陣風(fēng)鋒自動識別方法，包括以下步驟：

步驟一、依據(jù)陣風(fēng)鋒的雷達表現(xiàn)特征，利用局部二值化算法提取出弱窄帶回波疑似區(qū)域；步驟如下：

1-1)設(shè)一幅低仰角雷達圖像的大小為N×N，對中間[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]區(qū)域內(nèi)的像素點p_ij進行其是否屬于弱脊?fàn)顓^(qū)域的判斷；

1-2)求出區(qū)域[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]內(nèi)大于等于35dBZ的連通區(qū)域ω_i，i＝1,2，…s，計算該連通區(qū)域ω_i的面積s_i及該連通區(qū)域ω_i外包矩形長度l_i，對于面積s_i小于S₁或外包矩形長度l_i小于L₁的連通區(qū)域進行標(biāo)記，得到ω_j′，j＝1,2，…m，m≤s；

1-3)若像素點p_ij的反射率值f(i,j)∈x₁或f(i,j)∈x₂且f(i,j)∈ω_j′，其中x₁＝[5,35)dBZ，x₂＝[35,40)dBZ，執(zhí)行1-4)，否則，執(zhí)行1-5)

1-4)在以該像素點p_ij為中心的區(qū)域根據(jù)式(1)做卷積運算，得到卷積運算結(jié)果g₁(i,j)和g₂(i,j)，當(dāng)f(i,j)≥g₁(i,j)且f(i,j)>g₂(i,j)，則認為像素點p_ij屬于弱窄帶回波疑似區(qū)域Ω，并將該像素點p_ij置為前景，否則，置為背景；

1-5)將像素點p_ij置為背景，至此將低仰角雷達圖像轉(zhuǎn)化為一張二值圖；

1-6)計算步驟1-5)所形成的二值圖中各個連通區(qū)域的面積和外接矩形長度，將連通區(qū)域面積小于S₂或外接矩形長度小于L₂的連通區(qū)域置為背景；從而提取出弱窄帶回波疑似區(qū)域；

步驟二、對弱窄帶回波疑似區(qū)域進行分割、連接和篩選處理，得到弱脊?fàn)顜?yīng)的骨架圖像；步驟如下：

2-1)對步驟一提取出的弱窄帶回波疑似區(qū)域的輪廓進行除毛刺處理，再進行細化得到區(qū)域的骨架圖像A；

2-2)在上述骨架圖像A中的骨架交叉點和折點處斷開骨架，得到骨架圖像B，包括：

依照折點的特點通過計算骨架上某點兩側(cè)切線夾角識別出折點，即沿骨架行進，若由某點到其前側(cè)第n個點的向量與該點到其后側(cè)第n個點的向量的夾角小于135度，則認為該點為折點；

利用局部二值模式算子檢測出端點和交叉點，即針對步驟2-1)得到的骨架圖像，骨架圖像中取值為1的點p為可能的端點或交叉點，以該點p為中心，考察其3×3區(qū)域及5×5區(qū)域邊界的取值分布；其中，5×5區(qū)域邊界點若取值為1，但在5×5范圍內(nèi)與區(qū)域中心不連通，則將該點置為0，分別從上述3×3區(qū)域及5×5區(qū)域的左上角開始沿逆時針方向形成8位01鏈碼和16位01鏈碼來描述所述兩個區(qū)域邊界的取值分布；然后分別沿兩個區(qū)域邊界循環(huán)一周，記錄依次取值變化的次數(shù)n₃(p)和n₅(p)，若n₃(p)＝2，則點p為端點；若n₃(p)≥6或n₅(p)≥6，則點p為交叉點；

在上述折點和交叉點處斷開骨架，從而得到骨架圖像B；

2-3)通過端點匹配的方法將骨架圖像B中任意兩段不連通的曲線連接；

根據(jù)陣風(fēng)鋒走向平緩的特點，設(shè)：曲線l_i的端點A和曲線l_j(j≠i)的端點B的匹配條件如下：

式(2)中，L＝30km，Φ₁＝Φ₂＝Φ₃＝-0.7，l_i和l_j表示曲線i和j的長度，指由端點A沿曲線i在端點A處切線的方向指向遠離該曲線的向量，指由端點B沿曲線j在端點B處切線的方向指向遠離該曲線的向量；

當(dāng)端點A僅與一個端點符合匹配條件時，則該端點為端點B，連接端點A和端點B；

當(dāng)端點A與多個端點符合匹配條件時，則將其中值最小的端點作為端點B，連接端點A和端點B，其中，端點C是與端點A符合匹配條件的任意一端點，指由端點C沿其所在曲線在端點C處切線的方向指向遠離該曲線的向量；

從而形成骨架圖像C；

2-4)判斷低仰角雷達圖像是否存在單側(cè)脊區(qū)域，若存在，則剔除單側(cè)脊區(qū)域的骨架；

對骨架圖像C中的每條曲線分別進行PCA處理，至少得到第二主成分方向的單位向量e₂；

當(dāng)距離k分別為2,3,4,5km時，設(shè)變量α_k的初值為0，對每條曲線上m個點中每一點p_i，i＝1,2...m,求得所述曲線一側(cè)的位置q_i，q_i＝p_i+k e₂對于每一點p_i計算：

式(3)中，f(x)表示點x的反射率值；

設(shè)，所述曲線一側(cè)的加權(quán)比例

當(dāng)距離k分別為-2,-3,-4,-5km時，設(shè)，所述曲線另一側(cè)的加權(quán)比例為β₀，根據(jù)上述過程求得加權(quán)比例β₀；

若α₀<Th且β₀<Th，則認為該曲線與弱窄帶回波疑似區(qū)域中對應(yīng)的區(qū)域是弱脊?fàn)顜?，否則，該曲線與弱窄帶回波疑似區(qū)域中對應(yīng)的區(qū)域是單側(cè)脊區(qū)域，濾除該曲線；其中，Th＝0.75；形成骨架圖像D；

2-5)判斷低仰角雷達圖像是否存在虛線回波，若存在，則剔除虛線回波的骨架；

用2-4)中對骨架圖像C中的每條曲線的PCA處理結(jié)果，能得到骨架圖像D中每條曲線的第一主成分特征值λ₁、第一主成分方向的單位向量e₁和第二主成分特征值λ₂，當(dāng)?shù)诙鞒煞痔卣髦郸?sub>2/第一主成分特征值λ₁<0.05，且雷達中心點到曲線所擬合的直線的距離小于d_l，則認為該曲線與低仰角雷達圖像中對應(yīng)的區(qū)域是虛線回波，濾除該曲線；其中d_l＝3km；最終得到的骨架圖像即為弱脊?fàn)顜?yīng)的骨架圖像；

步驟三、由當(dāng)前時刻和前一時刻的兩幅低仰角雷達圖像得到光流場，將步驟二得到的弱脊?fàn)顜?yīng)的骨架圖像中前后時刻匹配的骨架擬定為疑似陣風(fēng)鋒，根據(jù)該疑似陣風(fēng)鋒的位置和速度與風(fēng)暴單體的位置和速度的關(guān)系及該疑似陣風(fēng)鋒的走向與速度的關(guān)系識別出陣風(fēng)鋒；步驟如下：

3-1)利用光流法由當(dāng)前時刻和前一時刻的兩幅低仰角雷達圖像得到光流場；

3-2)利用光流場信息將前一時刻的弱脊?fàn)顜?yīng)的骨架圖像的曲線移動到由步驟二得到的當(dāng)前時刻骨架圖像中的對應(yīng)位置，同時滿足下述條件一和條件二的曲線為當(dāng)前時刻和前一時刻的同一條弱窄帶回波區(qū)域?qū)?yīng)，即為疑似陣風(fēng)鋒；

條件一：前后時刻弱窄帶回波區(qū)域的交疊長度大于30％；

條件二：利用PCA處理后得到的前后時刻弱窄帶回波區(qū)域的第一主成分方向夾角小于30度；

3-3)若前一時刻弱窄帶回波某端比當(dāng)前時刻弱窄帶回波長10km以上時，則利用多出的此段曲線對步驟3-2)中認定的疑似陣風(fēng)鋒進行延長補全；

3-4)根據(jù)該疑似陣風(fēng)鋒的位置和速度與風(fēng)暴單體的位置和速度的關(guān)系及該疑似陣風(fēng)鋒的走向與速度的關(guān)系識別是否具有陣風(fēng)鋒，包括

判斷當(dāng)前時刻疑似陣風(fēng)鋒走向與疑似陣風(fēng)鋒速度方向銳角大于45度；

同時，當(dāng)前時刻疑似陣風(fēng)鋒與風(fēng)暴單體的速度方向與位置關(guān)系滿足下述條件(1)至條件(5)中的一條；則認定當(dāng)前時刻疑似陣風(fēng)鋒為陣風(fēng)鋒；并針對下一時刻低仰角雷達圖像，直接將步驟3-3)補全后的該時刻的疑似陣風(fēng)鋒直接認定為陣風(fēng)鋒；

(1)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒僅位于風(fēng)暴單體運動的前側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向與風(fēng)暴單體移動方向夾角小于30度；

(2)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒僅位于風(fēng)暴單體運動的右側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向與風(fēng)暴單體移動方向右轉(zhuǎn)90度后的方向夾角小于30度；

(3)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒僅位于風(fēng)暴單體運動的左側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向與風(fēng)暴單體移動方向左轉(zhuǎn)90度后的方向夾角小于30度；

(4)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒既位于風(fēng)暴單體前側(cè)40km內(nèi)且位于風(fēng)暴單體右側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向位于風(fēng)暴單體移動方向與風(fēng)暴單體移動方向右轉(zhuǎn)90度后的方向之間；

(5)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒既位于風(fēng)暴單體前側(cè)40km內(nèi)且位于風(fēng)暴單體左側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向位于風(fēng)暴單體移動方向與風(fēng)暴單體移動方向左轉(zhuǎn)90度后的方向之間；

上述步驟3-2)中，利用光流場信息將前一時刻的弱脊?fàn)顜?yīng)的骨架圖像的曲線移動到由步驟二得到的當(dāng)前時刻骨架圖像中的對應(yīng)位置后，將沒有同時滿足條件一和條件二的前一時刻的弱脊?fàn)顜?yīng)的骨架圖像的曲線保留兩個提掃疊加到下一循環(huán)的疑似陣風(fēng)鋒匹配過程中的前一時刻的骨架圖像中；

上述步驟3-3)補全后的疑似陣風(fēng)鋒用于替換下一循環(huán)中前一時刻對應(yīng)位置的曲線。

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是：

根據(jù)陣風(fēng)鋒在雷達強度場中的回波特性和幾何特征，首先通過局部區(qū)域二值化等方法提取弱窄帶回波，再根據(jù)陣風(fēng)鋒與風(fēng)暴單體之間的關(guān)系從弱窄帶回波中篩選出陣風(fēng)鋒，實現(xiàn)了準(zhǔn)確、完整自動檢測陣風(fēng)鋒，對災(zāi)害進行及時的預(yù)警，減少了經(jīng)濟損失和人員傷亡；并通過實驗驗證了本方法的有效性。

附圖說明

圖1(a)至圖1(e)為天氣雷達反射率圖圖例與低仰角反射率圖上的陣風(fēng)鋒，其中，圖1(a)至圖1(e)的左側(cè)顯示了圖例(下同)，右側(cè)顯示了低仰角反射率圖上的陣風(fēng)鋒；

圖2為卷積模板，其中每個方格代表一個像素點，也就是雷達圖中1km*1km的區(qū)域(下同)，圖中的黑色點代表要輸出的像素方格，大的加粗框表示一個卷積操作數(shù)矩陣，小的加粗框代表另一個卷積操作數(shù)矩陣。

圖3(a)至圖3(f)為局部二值化算法提取的前景區(qū)域示意圖，其中，圖3(a)顯示了單側(cè)脊，圖3(b)顯示了與其他弱脊?fàn)顓^(qū)域交疊的陣風(fēng)鋒，圖3(c)顯示了虛線回波，圖3d)顯示了非陣風(fēng)鋒的邊界層輻合線，圖3(e)顯示了非邊界層輻合線，圖3(f)顯示了真正的陣風(fēng)鋒；

圖4(a)至圖4(c)為骨架的局部區(qū)域，其中，圖4(a)中3×3區(qū)域的中心點為骨架的端點，圖4(b)中5×5區(qū)域的中心點為骨架的交叉點，圖4(c)中5×5區(qū)域的中心點既不是骨架的端點也不是骨架的交叉點；

圖5為本發(fā)明提供的一種用于氣象學(xué)中尋找陣風(fēng)鋒方法的流程圖；

圖6(a)和圖(b)為本發(fā)明提供的測試效果圖。

具體實施方式

下面結(jié)合附圖和具體實施例對本發(fā)明技術(shù)方案作進一步詳細描述，所描述的具體實施例僅對本發(fā)明進行解釋說明，并不用以限制本發(fā)明。

本發(fā)明提供了一種基于多普勒天氣雷達資料的陣風(fēng)鋒自動識別方法，本方法能自動檢測出陣風(fēng)鋒，對災(zāi)害進行及時的預(yù)警，減少經(jīng)濟損失和人員傷亡。

實施例：在多普勒雷達低仰角反射率圖中觀測到的陣風(fēng)鋒如圖1(a)至圖1(e)所示，其圖像特征包括但不限于：

1)陣風(fēng)鋒通常表現(xiàn)為在強回波外圍的一條弱窄帶回波，一般位于雷暴母體運動方向的前端或側(cè)面，它與雷暴母體脫離或部分粘連，見圖1(a)；

2)寬度變化范圍：2km-10km；

3)長度變化范圍：40km-300km；

4)反射率取值不定：在同一時刻，一條陣風(fēng)鋒區(qū)域上的回波強度通常不同，見圖1(b)，不同時刻，同一條陣風(fēng)鋒上的回波強度取值和分布也會變化(如圖1(b)和圖1(c)，日期為2006.06.14其中圖1(b)時間為09:49和圖1(c)時間為09:31)，但多集中在10dBz—30dBz這個范圍；

5)脊特征：陣風(fēng)鋒上的反射率值整體上高于其兩側(cè)的反射率值，但有時差值不大，呈現(xiàn)很弱的脊特征，且不排除局部區(qū)域或局部點這種特征消失；

6)不連續(xù)性：當(dāng)陣風(fēng)鋒與非降水回波混合在一起時，可能會被較高反射率值區(qū)域或者較低的反射率值區(qū)域“污染”而發(fā)生斷裂，見圖1(d)；

7)交疊性：一條陣風(fēng)鋒可能與其他弱窄帶回波相交甚至發(fā)生碰撞，使兩者相連或相交，見圖1(e)。

本發(fā)明基于多普勒天氣雷達資料的陣風(fēng)鋒自動識別方法，如圖5所示，包括以下步驟：

101、依據(jù)陣風(fēng)鋒的雷達表現(xiàn)特征，利用局部二值化算法提取出弱窄帶回波疑似區(qū)域；具體內(nèi)容如下：

1-1)設(shè)一幅低仰角雷達圖像的大小為N×N，對中間[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]區(qū)域內(nèi)的像素點p_ij進行其是否屬于弱脊?fàn)顓^(qū)域的判斷；

1-2)求出區(qū)域[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]內(nèi)大于等于35dBZ的連通區(qū)域ω_i，i＝1,2，…s，計算該連通區(qū)域ω_i的面積s_i及該連通區(qū)域ω_i外包矩形長度l_i，對于面積s_i小于S₁或外包矩形長度l_i小于L₁的連通區(qū)域進行標(biāo)記，得到ω_j′，j＝1,2，…m，m≤s；

1-3)若像素點p_ij的反射率值f(i,j)∈x₁或f(i,j)∈x₂且f(i,j)∈ω_j′，其中x₁＝[5,35)dBZ，x₂＝[35,40)dBZ，執(zhí)行1-4)，否則，執(zhí)行1-5)

1-4)則用圖2所示的模板在在以該像素點p_ij為中心的區(qū)域根據(jù)式(1)做卷積運算，得到卷積運算結(jié)果g₁(i,j)和g₂(i,j)，當(dāng)f(i,j)≥g₁(i,j)且f(i,j)>g₂(i,j)，則認為像素點p_ij屬于弱窄帶回波區(qū)域Ω，并將該像素點p_ij置為前景，否則，置為背景；

1-5)將像素點p_ij置為背景，至此將低仰角雷達圖像轉(zhuǎn)化為一張二值圖；

1-6)計算步驟1-5)所形成的二值圖中各個連通區(qū)域的面積和外接矩形長度，將連通區(qū)域面積小于S₂或外接矩形長度小于L₂的連通區(qū)域置為背景；其中參數(shù)k＝3，S1＝15，L1＝5，S2＝45，L2＝15；圖3(a)、圖3(b)、圖3(c)、圖3(d)、圖3(e)和圖3(f)顯示了局部二值化算法提取的前景區(qū)域示意圖，圖3(a)顯示了單側(cè)脊，圖3(b)顯示了與其他弱脊?fàn)顓^(qū)域交疊的陣風(fēng)鋒，圖3(c)顯示了虛線回波，圖3(d)顯示了非陣風(fēng)鋒的邊界層輻合線，圖3(e)顯示了非邊界層輻合線，圖3(f)顯示了真正的陣風(fēng)鋒；從而提取出弱窄帶回波疑似區(qū)域。

102、對弱窄帶回波疑似區(qū)域進行分割、連接和篩選處理，得到弱脊?fàn)顜?yīng)的骨架圖像；具體內(nèi)容如下：

2-1)對步驟一提取出的弱窄帶回波疑似區(qū)域的輪廓進行除毛刺處理，再進行細化得到區(qū)域的骨架圖像A；

2-2)在上述骨架圖像A中的骨架交叉點和折點處斷開骨架，得到骨架圖像B，包括：

依照折點的特點通過計算骨架上某點兩側(cè)切線夾角識別出折點，即沿骨架行進，若由某點到其前側(cè)第n個點的向量與該點到其后側(cè)第n個點的向量的夾角小于135度，則認為該點為折點，本例中n＝3；

利用局部二值模式(Local Binary Pattern，LBP)^[4]算子檢測出端點和交叉點，即針對步驟2-1)得到的骨架圖像，骨架圖像中取值為1的點p為可能的端點或交叉點，以該點p為中心，考察其3×3區(qū)域及5×5區(qū)域邊界的取值分布；其中，5×5區(qū)域邊界點若取值為1，但在5×5范圍內(nèi)與區(qū)域中心不連通，如圖4(c)中覆蓋較深(第4行最右邊)的方格，則將該點置為0，分別從上述3×3區(qū)域及5×5區(qū)域的左上角開始沿逆時針方向形成8位01鏈碼和16位01鏈碼來描述所述兩個區(qū)域邊界的取值分布；然后分別沿兩個區(qū)域邊界循環(huán)一周，記錄依次取值變化的次數(shù)n₃(p)和n₅(p)，若n₃(p)＝2，如圖4(a)，則點p為端點；若n₃(p)≥6或n₅(p)≥6，如圖4(b)，則點p為交叉點；

在上述折點和交叉點處斷開骨架，從而得到骨架圖像B；

2-3)通過端點匹配的方法將骨架圖像B中任意兩段不連通的曲線連接；

根據(jù)陣風(fēng)鋒走向平緩的特點，設(shè)：曲線l_i的端點A和曲線l_j(j≠i)的端點B的匹配條件如下：

式(2)中，L＝30km，Φ₁＝Φ₂＝Φ₃＝-0.7，l_i和l_j表示曲線i和j的長度，指由端點A沿曲線i在端點A處切線的方向指向遠離該曲線的向量，指由端點B沿曲線j在端點B處切線的方向指向遠離該曲線的向量；

當(dāng)端點A僅與一個端點符合匹配條件時，則該端點為端點B，連接端點A和端點B；

當(dāng)端點A與多個端點符合匹配條件時，則將其中值最小的端點作為端點B，連接端點A和端點B，其中，端點C是與端點A符合匹配條件的任意一端點，指由端點C沿其所在曲線在端點C處切線的方向指向遠離該曲線的向量；

從而形成骨架圖像C；

2-4)判斷低仰角雷達圖像是否存在單側(cè)脊區(qū)域，若存在，則剔除單側(cè)脊區(qū)域的骨架；

對骨架圖像C中的每條曲線分別進行PCA處理，至少得到第二主成分方向的單位向量e₂；

當(dāng)距離k分別為2,3,4,5km時，設(shè)變量α_k的初值為0，對每條曲線上m個點中每一點p_i，i＝1,2...m,求得所述曲線一側(cè)的位置q_i，q_i＝p_i+k e₂對于每一點p_i計算：

式(3)中，f(x)表示點x的反射率值；

設(shè)，所述曲線一側(cè)的加權(quán)比例

當(dāng)距離k分別為-2,-3,-4,-5km時，設(shè)，所述曲線另一側(cè)的加權(quán)比例為β₀，根據(jù)上述過程求得加權(quán)比例β₀；

若α₀<Th且β₀<Th，則認為該曲線與弱窄帶回波疑似區(qū)域中對應(yīng)的區(qū)域是弱脊?fàn)顜?，否則，該曲線與弱窄帶回波疑似區(qū)域中對應(yīng)的區(qū)域是單側(cè)脊區(qū)域，濾除該曲線；其中，Th＝0.75；形成骨架圖像D；

2-5)判斷低仰角雷達圖像是否存在虛線回波，若存在，則剔除虛線回波的骨架；

用2-4)中對骨架圖像C中的每條曲線的PCA處理結(jié)果，能得到骨架圖像D中每條曲線的第一主成分特征值λ₁、第一主成分方向的單位向量e₁和第二主成分特征值λ₂，當(dāng)?shù)诙鞒煞痔卣髦郸?sub>2/第一主成分特征值λ₁<0.05，且雷達中心點到曲線所擬合的直線的距離小于d_l，則認為該曲線與低仰角雷達圖像中對應(yīng)的區(qū)域是虛線回波，濾除該曲線；其中d_l＝3km；最終得到的骨架圖像即為弱脊?fàn)顜?yīng)的骨架圖像。

103、由當(dāng)前時刻和前一時刻的兩幅低仰角雷達圖像得到光流場，將弱脊?fàn)顜?yīng)的骨架圖像中前后時刻匹配的骨架擬定為疑似陣風(fēng)鋒，根據(jù)該疑似陣風(fēng)鋒的位置和速度與風(fēng)暴單體的位置和速度的關(guān)系及該疑似陣風(fēng)鋒的走向與速度的關(guān)系識別出陣風(fēng)鋒，具體內(nèi)容如下：

3-1)利用光流法^[5]由當(dāng)前時刻和前一時刻的兩幅低仰角雷達圖像得到光流場，

3-2)利用光流場信息將前一時刻的弱脊?fàn)顜?yīng)的骨架圖像的曲線移動到由步驟二得到的當(dāng)前時刻骨架圖像中的對應(yīng)位置，同時滿足下述條件一和條件二的曲線為當(dāng)前時刻和前一時刻的同一條弱窄帶回波區(qū)域?qū)?yīng)，即為疑似陣風(fēng)鋒；

條件一：前后時刻弱窄帶回波區(qū)域的交疊長度大于30％；

條件二：利用PCA處理后得到的前后時刻弱窄帶回波區(qū)域的第一主成分方向夾角小于30度；

3-3)若前一時刻弱窄帶回波某端比當(dāng)前時刻弱窄帶回波長10km以上時，則利用多出的此段曲線對步驟3-2)中認定的疑似陣風(fēng)鋒進行延長補全；

3-4)根據(jù)該疑似陣風(fēng)鋒的位置和速度與風(fēng)暴單體的位置和速度的關(guān)系及該疑似陣風(fēng)鋒的走向與速度的關(guān)系識別是否具有陣風(fēng)鋒，包括

判斷當(dāng)前時刻疑似陣風(fēng)鋒走向與疑似陣風(fēng)鋒速度方向所夾銳角大于45度；

同時，當(dāng)前時刻疑似陣風(fēng)鋒與風(fēng)暴單體的速度方向與位置關(guān)系滿足下述條件(1)至條件(5)中的一條；則認定當(dāng)前時刻疑似陣風(fēng)鋒為陣風(fēng)鋒；并針對下一時刻低仰角雷達圖像，直接將步驟3-3)補全后的該時刻的疑似陣風(fēng)鋒直接認定為陣風(fēng)鋒；

(1)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒僅位于風(fēng)暴單體運動的前側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向與風(fēng)暴單體移動方向夾角小于30度；

(2)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒僅位于風(fēng)暴單體運動的右側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向與風(fēng)暴單體移動方向右轉(zhuǎn)90度后的方向夾角小于30度；

(3)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒僅位于風(fēng)暴單體運動的左側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向與風(fēng)暴單體移動方向左轉(zhuǎn)90度后的方向夾角小于30度；

(4)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒既位于風(fēng)暴單體前側(cè)40km內(nèi)且位于風(fēng)暴單體右側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向位于風(fēng)暴單體移動方向與風(fēng)暴單體移動方向右轉(zhuǎn)90度后的方向之間；

(5)當(dāng)當(dāng)前時刻疑似陣風(fēng)鋒既位于風(fēng)暴單體前側(cè)40km內(nèi)且位于風(fēng)暴單體左側(cè)40km內(nèi)時，疑似陣風(fēng)鋒速度方向位于風(fēng)暴單體移動方向與風(fēng)暴單體移動方向左轉(zhuǎn)90度后的方向之間；

上述步驟3-2)中，利用光流場信息將前一時刻的弱脊?fàn)顜?yīng)的骨架圖像的曲線移動到由步驟二得到的當(dāng)前時刻骨架圖像中的對應(yīng)位置后，將沒有同時滿足條件一和條件二的前一時刻的弱脊?fàn)顜?yīng)的骨架圖像的曲線保留兩個提掃疊加到下一循環(huán)的疑似陣風(fēng)鋒匹配過程中的前一時刻的骨架圖像中；

上述步驟3-3)補全后的疑似陣風(fēng)鋒用于替換下一循環(huán)中前一時刻對應(yīng)位置的曲線。

下面以具體的實驗來驗證本發(fā)明實施例提供的一種基于多普勒天氣雷達資料的陣風(fēng)鋒自動識別方法的可行性，測試樣本由中國天津氣象臺提供，分兩部分驗證，詳見下文描述：

圖6(a)和圖6(b)為本發(fā)明提供的測試效果圖，其中，圖6(b)標(biāo)注的白色區(qū)域即最后自動檢測的真實陣風(fēng)鋒，第一部分是天津地區(qū)含有陣風(fēng)鋒的6個過程115組雷達基數(shù)據(jù)。第二部分是天津2012年6月含有強對流天氣的7天1680個基數(shù)據(jù)，其中包含5個陣風(fēng)鋒過程的88個基數(shù)據(jù)。利用擊中率POD、虛警率FAR和臨界成功指數(shù)CSI對檢驗結(jié)果進行評價(見表1和表2)。

表1 第一部分樣本數(shù)和識別情況

表1分別描述了6個陣風(fēng)鋒過程的識別情況，對其進行統(tǒng)計后得到樣本總數(shù)為115個，成功識別樣本數(shù)93個，未識別樣本數(shù)22，誤識別樣本數(shù)3，得出擊中率、虛警率和臨界成功指數(shù)分別為80.87％,3.13％,78.81％。

表2 第二部分2012年6月天津樣本數(shù)和識別情況

表2分別描述了7天強對流天氣的識別情況，擊中了88個陣風(fēng)鋒中的51個，在1592個數(shù)據(jù)中空報出了15個陣風(fēng)鋒，得出擊中率、虛警率和臨界成功指數(shù)分別為57.95％,22.73％,49.51％。

本領(lǐng)域技術(shù)人員可以理解附圖只是一個優(yōu)選實施例的示意圖，上述本發(fā)明實施例序號僅僅為了描述，不代表實施例的優(yōu)劣。以上所述僅為本發(fā)明的較佳實施例，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進等，均應(yīng)包含在本發(fā)明的保護范圍之內(nèi)。