甘建紅， 漆 慧,胡文東， 舒紅平， 羅 飛， 何童麗,尹麒名， 黎仁國(guó)

(1.成都信息工程大學(xué)， 成都 610200； 2.西華師范大學(xué)， 南充 637002)

1 引 言

大氣風(fēng)場(chǎng)急流現(xiàn)象最早在20世紀(jì)30年代被人們發(fā)現(xiàn)，急流與對(duì)強(qiáng)對(duì)流天氣形成發(fā)展、風(fēng)能轉(zhuǎn)化、空氣污染物的輸送和擴(kuò)散、森林火災(zāi)的蔓延等有著密切聯(lián)系[1-2].目前實(shí)際氣象業(yè)務(wù)中通過(guò)人工手動(dòng)方式在軟件上進(jìn)行識(shí)別和繪制急流，一定程度上限制了過(guò)程自動(dòng)化、預(yù)報(bào)智能化的快速發(fā)展.隨著智能天氣預(yù)報(bào)和格點(diǎn)預(yù)報(bào)興起，計(jì)算機(jī)自動(dòng)識(shí)別并標(biāo)識(shí)氣象風(fēng)場(chǎng)急流的研究逐漸受到人們重視.

Molnos等[3]以風(fēng)速、風(fēng)向和緯度導(dǎo)向函數(shù)構(gòu)造代價(jià)函數(shù)，通過(guò)Dijkstra方法尋找最短路徑方式計(jì)算副熱帶急流(Subtropical Jet Stream)和極鋒急流(Polar Jet Steam)，由于該急流線是兩條的閉合線，沒(méi)有考慮分叉和急流合流的現(xiàn)象，且算法需要已知起點(diǎn)，可見(jiàn)該算法不適合于非全球范圍風(fēng)場(chǎng)數(shù)據(jù)中計(jì)算具有多個(gè)急流的情況.

金宏憶[4]以VAD算法為基礎(chǔ)，用雷達(dá)探測(cè)資料以及多普勒天氣雷達(dá)二次產(chǎn)品，探討了機(jī)器識(shí)別低空急流的檢測(cè)方法.章晗等人[5]也利用雷達(dá)資料研究急流識(shí)別，該算法只針對(duì)西南風(fēng)低空急流自動(dòng)識(shí)別做出了研究，未對(duì)其它類型的急流提取回波特征.

路志英等[6]和王萍等[7]的方法類似，對(duì)站點(diǎn)風(fēng)數(shù)據(jù)計(jì)算急流線并用MICAPS進(jìn)行顯示，該方法利用地面站的觀測(cè)資料，經(jīng)過(guò)預(yù)處理、聚類及擬合實(shí)現(xiàn)了850 hPa 等壓面上中國(guó)東南部沿海地區(qū)低空急流的繪制.該方法根據(jù)該地區(qū)以東南風(fēng)和西南風(fēng)為主的特征進(jìn)行識(shí)別，對(duì)其它地方適用性需要進(jìn)一步優(yōu)化.

隨著智能網(wǎng)格預(yù)報(bào)的興起，不但預(yù)報(bào)數(shù)據(jù)采用格點(diǎn)形式，實(shí)況數(shù)據(jù)也逐漸采用格點(diǎn)數(shù)據(jù)形式進(jìn)行應(yīng)用[8-9].本文針對(duì)格點(diǎn)類型的風(fēng)場(chǎng)數(shù)據(jù)進(jìn)行急流軸計(jì)算，首先通過(guò)方向相干性計(jì)算急流區(qū)域，然后采用骨架提取算法獲取接近幾何中心的格點(diǎn)，基于獲取格點(diǎn)位置采用重心公式計(jì)算急流軸的點(diǎn)，最后連接急流軸的點(diǎn)得到急流軸.

2 急流計(jì)算方法

算法主要經(jīng)過(guò)三個(gè)階段實(shí)現(xiàn)：計(jì)算急流區(qū)域，計(jì)算中心軸(骨架)點(diǎn)及其微調(diào)，連接點(diǎn)得到急流線，具體算法流程如圖1所示.

3 計(jì)算急流區(qū)域

大氣風(fēng)場(chǎng)急流具有風(fēng)速較大、局部區(qū)域內(nèi)風(fēng)向具有一致性特點(diǎn).對(duì)于站點(diǎn)數(shù)據(jù)可采用圖的方法計(jì)算點(diǎn)間的相似性，而本文輸入格點(diǎn)數(shù)據(jù)為棋盤(pán)格分布的數(shù)據(jù)，因此采用Hessian的方法計(jì)算相鄰點(diǎn)的相干系數(shù).本文分兩步計(jì)算激流區(qū)域：(1) 計(jì)算風(fēng)向一致性好的區(qū)域；(2) 從得到區(qū)域中去掉風(fēng)速小的地方，保留下來(lái)的區(qū)域即為急流區(qū)域.

為了便于使用算子進(jìn)行計(jì)算，在格點(diǎn)數(shù)據(jù)的上邊增加一行數(shù)據(jù)，其值等于最邊上一行數(shù)據(jù)，同理在格點(diǎn)的下邊增加一行數(shù)據(jù)，在左邊、右邊也各增加一列數(shù)據(jù)，將增加數(shù)據(jù)之前的格點(diǎn)稱為原格點(diǎn).

圖1 急流自動(dòng)識(shí)別流程圖Fig.1 Jet automatic identification

(1) 計(jì)算風(fēng)向一致性較好的區(qū)域.

首先，對(duì)全部原格點(diǎn)數(shù)據(jù)計(jì)算方向相干性系數(shù).在圖像處理等領(lǐng)域內(nèi)相干性系數(shù)計(jì)算主要面向標(biāo)量場(chǎng)，而風(fēng)場(chǎng)是一個(gè)矢量場(chǎng)，為了能讓其刻畫(huà)方向保持一致性特性，提出算子范圍內(nèi)水平、豎直方向的最大風(fēng)向夾角作為X、Y方向差分并構(gòu)建Hessian矩陣[10]的方法計(jì)算相干性系數(shù)，得到Hessian矩陣在此表達(dá)式為

(1)

其中，F(xiàn)x=max{θ|鄰域內(nèi)水平方向兩兩格點(diǎn)風(fēng)向夾角};Fy=max{θ|鄰域內(nèi)豎直方向兩兩格點(diǎn)風(fēng)向夾角}.

具體離散點(diǎn)中計(jì)算方法如下：假設(shè)當(dāng)前計(jì)算風(fēng)向相干系數(shù)格點(diǎn)位置為(x0,y0)，本文實(shí)驗(yàn)中鄰域尺寸為3×3，如圖2所示.計(jì)算(x0,y0)點(diǎn)風(fēng)向與(x0+1,y0)點(diǎn)風(fēng)向構(gòu)成的夾角，記為Fx1；計(jì)算(x0,y0)點(diǎn)風(fēng)向與(x0-1,y0)點(diǎn)風(fēng)向構(gòu)成的夾角，記為Fx2；計(jì)算(x0+1,y0)點(diǎn)風(fēng)向與(x0-1,y0)點(diǎn)風(fēng)向構(gòu)成的夾角，記為Fx3.于是,Fx1、Fx2、Fx3的最大值即為Fx.同理，得到y(tǒng)方向的差分值Fy.

然后，計(jì)算該Hessian矩陣的特征值λ1和λ2，并按照式(2)計(jì)算風(fēng)向相干性系數(shù)k.如果原格點(diǎn)的系數(shù)小于閾值T，則為急流區(qū)域候選點(diǎn).

k=|λ1-λ2|

(2)

(2) 提取急流區(qū)域.

按照公式(3)計(jì)算每個(gè)候選格點(diǎn)風(fēng)矢量的模值，即風(fēng)大小，如果模值大于閾值S，則為急流區(qū)域內(nèi)的格點(diǎn)，否則不是急流區(qū)域的格點(diǎn).至此可得到一個(gè)或多個(gè)急流區(qū)域.

(3)

其中,vx、vy分別為風(fēng)矢量在水平、豎直方向的分量.

(a) (b) (c)

4 計(jì)算急流區(qū)域幾何中心軸

《中尺度天氣圖技術(shù)規(guī)范(修改稿)》，低空急流的技術(shù)要求是：當(dāng)有兩個(gè)以上連續(xù)測(cè)站風(fēng)速超過(guò)12 m/s時(shí)， 沿12 m/s以上大風(fēng)區(qū)的幾何中心分析低空急流.基于該規(guī)范的思想，算法在計(jì)算出急流區(qū)域后利用骨架提取方法計(jì)算區(qū)域的骨架點(diǎn)，骨架點(diǎn)即也是格點(diǎn).

4.1 計(jì)算急流軸附近的格點(diǎn)

骨架提取算法[11-12]是提取幾何中心軸常用的方法，而查表法是骨架提取的一個(gè)重要方法.查表法將一系列約束條件用周圍像素權(quán)重乘積之和通過(guò)查表方式進(jìn)行表達(dá)，比如內(nèi)部點(diǎn)不能刪除、孤立點(diǎn)不能刪除、直線端點(diǎn)不能刪除等約束條件.在此，每個(gè)格點(diǎn)的值只能為1或0，如果為急流區(qū)域的點(diǎn)，則為1，否則為0.權(quán)重分布如圖3所示，對(duì)于任意一格點(diǎn)，取其周圍緊鄰8個(gè)格點(diǎn)為鄰域，并按照一定順序給與權(quán)重，權(quán)重分別為1、2、4、8、16、32、64、128，中間方格表示當(dāng)前格點(diǎn).

圖3 鄰域權(quán)重分布圖

計(jì)算方法為將當(dāng)前格點(diǎn)周圍8個(gè)格點(diǎn)的值分別乘以對(duì)應(yīng)位置的權(quán)重，并對(duì)乘積之和除以255得到最終結(jié)果，在表1中查詢?cè)撟罱K結(jié)果對(duì)應(yīng)的值，對(duì)應(yīng)值如果是1，則表示這個(gè)點(diǎn)應(yīng)該刪掉，反之若值為0，表示這個(gè)點(diǎn)應(yīng)該保留.

4.2 微調(diào)急流軸點(diǎn)

通過(guò)上述算法獲取的急流軸位于格點(diǎn)之上，但在實(shí)際情況中，風(fēng)場(chǎng)的急流軸剛好經(jīng)過(guò)格點(diǎn)的概率非常小.采用重心公式對(duì)骨架點(diǎn)進(jìn)行微調(diào)，使得骨架點(diǎn)接近風(fēng)速最大位置，與實(shí)際情況一致.具體方法為，設(shè)第i個(gè)骨架點(diǎn)為當(dāng)前待微調(diào)的骨架點(diǎn)，將當(dāng)前骨架點(diǎn)及其左、右、上、下、左上、右上、左下、右下共9個(gè)格點(diǎn)的風(fēng)矢量的模帶入式(4)和(5)，9個(gè)格點(diǎn)中如果某個(gè)格點(diǎn)不在急流區(qū)域?qū)⑵滹L(fēng)大小F設(shè)置為0進(jìn)行計(jì)算，得到微調(diào)后的骨架點(diǎn)位置(Xi，Yi).

(4)

(5)

其中,xij，yij分別是當(dāng)前骨架點(diǎn)的橫坐標(biāo)與縱坐標(biāo);Fij是格點(diǎn)的風(fēng)矢量的模.圖4(a)中圓點(diǎn)標(biāo)注的是計(jì)算得出的急流區(qū)域，五角星標(biāo)注的是微調(diào)后的骨架點(diǎn).

最后,將計(jì)算的急流軸點(diǎn)(微調(diào)后的骨架點(diǎn))進(jìn)行連接得到急流線，如圖4所示，通常需要對(duì)其進(jìn)行平滑處理，平滑處理不是本文研究重點(diǎn)且方法較為成熟，不在贅述.

(a)

(b)

5 算法測(cè)試結(jié)果及分析

算法實(shí)驗(yàn)采用MeteoInfo開(kāi)源組件讀取數(shù)據(jù)，在VisualStudio2015開(kāi)發(fā)環(huán)境中使用C#語(yǔ)言實(shí)現(xiàn)，輸入數(shù)據(jù)為MICAPS第11類數(shù)據(jù)，程序結(jié)果為急流軸的經(jīng)緯度坐標(biāo).

5.1 統(tǒng)計(jì)分析

將算法應(yīng)用于2010年1月1日至7月1日及2017年7月1日至12月31日08～20時(shí)500 hPa風(fēng)場(chǎng)數(shù)據(jù)中的352組數(shù)據(jù)進(jìn)行測(cè)試.將測(cè)試結(jié)果交若干位專家進(jìn)行核實(shí)，如果專家認(rèn)為與手動(dòng)繪制較吻合則定義為準(zhǔn)確擊中；如果與人工繪制偏差較大但是流線方向與風(fēng)向趨勢(shì)一致則定義為偏差擊中；若本算法識(shí)別存在低空急流但人工分析不存在低空急流則認(rèn)為空?qǐng)?bào)；反之，若人工分析存在低空急流但本算法并未識(shí)別出，則認(rèn)為漏報(bào).實(shí)驗(yàn)結(jié)果如表2所示.

圖5是2010年1月1日至4日08時(shí)高度為500 hPa的UV數(shù)據(jù)實(shí)驗(yàn)結(jié)果，紅色線條為計(jì)算的急流線，藍(lán)色為表示風(fēng)向的風(fēng)廓線.算法獲得的急流走向和風(fēng)向走勢(shì)具有良好的一致性.

表2 算法應(yīng)用于2012年和2017年部分?jǐn)?shù)據(jù)結(jié)果

(a)

(b)

(c)

(d)

5.2 個(gè)例分析

將本算法應(yīng)用于2017年7月14日8時(shí)150 hPa的UV數(shù)據(jù),結(jié)果如圖6所示.圖b中圓形點(diǎn)為急流區(qū)域內(nèi)的格點(diǎn)，圓形格點(diǎn)組成的區(qū)域即為急流區(qū)域.五角星為微調(diào)后得到的骨架點(diǎn)，粗線即為急流軸.通過(guò)自動(dòng)檢測(cè)能快速檢測(cè)出急流軸，對(duì)于小短的急流軸也能找到，而這些通常是人工識(shí)別容易忽略的急流軸.急流的合并與分流是描述風(fēng)場(chǎng)變化的重要特性，算法計(jì)算的急流軸也保持了合并與分流的情況.

文獻(xiàn)[7]對(duì)2009年3月21日8時(shí)、2009年4月19日8時(shí)地面站數(shù)據(jù)進(jìn)行了急流線分析，本文針對(duì)2009年3月21日8時(shí)、2009年4月19日8時(shí)的實(shí)況格點(diǎn)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，本文算法的急流線在北海附近更靠近最大風(fēng)速、更靠近風(fēng)力大的站點(diǎn)，與急流線的定義更一致.由于文獻(xiàn)[7]基于站點(diǎn)數(shù)據(jù)計(jì)算急流，因此本文算法急流經(jīng)過(guò)海洋較多，文獻(xiàn)[7]的急流主要在陸地，由于數(shù)據(jù)源的不同，結(jié)果存在一定差異.

圖6 高空風(fēng)場(chǎng)計(jì)算結(jié)果Fig.6 Result of high level UV

一定程度上，風(fēng)場(chǎng)本身也是流體，根據(jù)流線的性質(zhì)，二維平面內(nèi)同一時(shí)刻的不同流線不能相交,因此，本實(shí)驗(yàn)結(jié)果沒(méi)有出現(xiàn)急流相交的情況，更具有可解釋性，符合流體力學(xué)性質(zhì).

6 結(jié) 論

本文根據(jù)大氣風(fēng)場(chǎng)格點(diǎn)矢量數(shù)據(jù)，利用相干性系數(shù)獲取特征值，將格點(diǎn)矢量數(shù)據(jù)轉(zhuǎn)為標(biāo)量.再利用骨架提取的方法可以自動(dòng)識(shí)別急流區(qū)域并繪制急流線.相比人工處理效率更高、準(zhǔn)確度更高，為天氣預(yù)報(bào)自動(dòng)化、智能化提供方法參考.為了得到更好的結(jié)果，可對(duì)原始格點(diǎn)數(shù)據(jù)進(jìn)行插值預(yù)處理，比如基于徑向基的插值方法.

根據(jù)《中尺度天氣圖技術(shù)規(guī)范(修改稿)》，低空急流的技術(shù)要求是：當(dāng)有兩個(gè)以上連續(xù)測(cè)站風(fēng)速超過(guò)12 m/s時(shí)， 沿12 m/s以上大風(fēng)區(qū)的幾何中心分析低空急流.可見(jiàn)算法首先計(jì)算急流區(qū)域然后通過(guò)骨架提取幾何中心軸與該規(guī)范思路一致.