本發(fā)明涉及氣象學領(lǐng)域，特別涉及一種基于光流法的單體跟蹤算法和一種基于移動立方體算法的單體三維重建算法。

背景技術(shù)：

冰雹、強降水、雷雨大風等強對流天氣，具有空間尺度小，生命周期短，突發(fā)性強，破壞力強等特點。同時，不同類型的強對流單體，其內(nèi)部結(jié)構(gòu)和粒子分布具備差異。例如冰雹云強回波的中心值、強回波核(45-65dbz)高度均高于強降水回波，易出現(xiàn)懸垂結(jié)構(gòu)^[1-4]，體現(xiàn)為冰雹云一般于降雹前在高、強、垂三個方面有別于單純性短時強降水；且降雹前vil值會不同程度持續(xù)躍增，降雹時vil值明顯大于短時強降水^[4]。

現(xiàn)有雷達產(chǎn)品存在以下缺點和不足：

(1)雷達圖像現(xiàn)采取的ppi顯示方式是在一系列固定仰角上掃描360°進行采樣的，即在某一仰角，雷達天線繞垂直軸z進行360°掃描^[5]，直觀性不高。

(2)組合反射率圖只顯示了最高反射率，無法觀察反射率在高度上的分布。

(3)剖面不易制作，人機交互體驗差。

(4)不提供強對流單體的結(jié)構(gòu)化特征，如描述強上升氣流造成的回波懸掛結(jié)構(gòu)的懸垂度，描述單體內(nèi)核體反射率強度的高度分布、梯度分布等；更談不上對它們演變規(guī)律的展示。

[參考文獻]

[1]李德俊，唐仁茂，熊守全等；強冰雹和短時強降水天氣雷達特征及臨近預(yù)警[j]，氣象，2011，37(4)：474-480。

[2]井喜，屠妮妮，井宇等；毛烏素沙漠東部邊緣中β尺度暴雨和冰雹綜合分析[j]，災(zāi)害學，2010，25(1)：73-77。

[3]李湘，張騰飛，胡娟等；云南冰雹災(zāi)害的多普勒雷達特征統(tǒng)計及預(yù)警指標[j]，災(zāi)害學，2015(2015年03)：88-93，98。

[4]張崇莉，向明堃，賴云華等；滇西北高原冰雹、短時強降水的多普勒雷達回波特征比較[j]，暴雨災(zāi)害，2011，30(1)：64-69。

[5]俞小鼎等；多普勒天氣雷達原理與業(yè)務(wù)應(yīng)用[m]，北京：氣象出版社，2009。

[6]李聰.強冰雹自動識別技術(shù)與防雹作業(yè)決策方法研究[d].天津：天津大學,2014：22-31

[7]王萍，潘躍.基于顯著性特征的大冰雹識別模型[j].物理學報,2013,62(6):515-524

[8]amburnsa,wolfpl.vildensityasahailindicator.weatherandforecasting,1997,12(3):473-478

[9]greenedr,clarkra.anindicatorofexplosivedevelopmentinseverestorms[c]7thconferenceofseverelocalstorms,missouri.1971.

[10]barronjl,fleetdj,beaucheminss.performanceofopticalflowtechniques.internationaljournalofcomputervision,1994,12(1):43-77

[11]lorensenwe,clinehe.marchingcubes:ahighresolution3dsurfaceconstructionalgorithm[c]acmsiggraphcomputergraphics.acm,1987,21(4):163-169.

[12]wangping,licong,zhangyuan.anadaptivesegmentationarithmeticadaptedtointertwinedirregularconvectivestormimages[c]，icml2013,vol.2,pp896-898

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

針對上述現(xiàn)有技術(shù)，本發(fā)明提供一種強對流單體結(jié)構(gòu)及結(jié)構(gòu)特征可視化時空剖析方法，可以用一組結(jié)構(gòu)化特征的變化曲線表現(xiàn)強對流單體的特征變化趨勢和變化規(guī)律；可以通過透明度的設(shè)置為用戶展現(xiàn)對流單體更好的三維透視性，通過透明度參數(shù)的設(shè)置為用戶展現(xiàn)特定反射率值域的三維構(gòu)成與分布；并改進了垂直剖面的人機交互體驗；通過交互式剖線角參量的設(shè)置，為用戶展現(xiàn)反射率強度在不同方位的垂直面上的分布；增加了水平剖面的人機交互體驗，通過交互式高度參量的設(shè)置，為用戶展現(xiàn)在不同高度的水平面上的分布。

為了解決上述技術(shù)問題，本發(fā)明提出的一種強對流單體結(jié)構(gòu)及結(jié)構(gòu)特征可視化時空剖析方法，包括獲得強對流單體特征的變化趨勢，獲得強對流單體三維圖像及垂直和水平的剖面；

一、獲得強對流單體特征的變化趨勢的步驟如下：

步驟1、選取強對流單體特征

以雷達基數(shù)據(jù)為原始數(shù)據(jù)，求取獲得有效厚度、高回波比、懸垂度、加權(quán)核高、峰度、組合反射率圖核區(qū)梯度高值、垂直累計液態(tài)水含量和累計液態(tài)水密度；

同時，通過對雷達基數(shù)據(jù)進行插值得到規(guī)則三維場數(shù)據(jù)，根據(jù)規(guī)則三維場數(shù)據(jù)得到各反射率的高度分布，并用一組柱形圖表示；

所述組合反射率圖核區(qū)梯度高值用于量化單體核區(qū)的變化程度，其求取方法如下：

在組合反射率圖中進行單體分割后，將單體中反射率r>40dbz的區(qū)域定義為單體核區(qū)，按照下述公式求取單體核區(qū)的中心p，

公式(1)中，pr>40dbz為所有反射率大于40dbz的點的坐標，rr>40dbz為對應(yīng)反射率；

以單體核區(qū)的中心p為圓心作圓覆蓋單體，并將該圓等分為八個扇形，選取三個相鄰扇形使之核區(qū)面積之和最小，并求取這三個扇形的平均梯度g：

公式(2)中，∑g為三個扇形區(qū)域所有像素點的梯度之和，像素點的梯度只考慮大小，不考慮方向，num為像素點總數(shù)，該平均梯度g即為組合反射率圖核區(qū)梯度高值；

步驟2、通過光流法實現(xiàn)前后體掃時刻的單體追蹤：

利用光流法求得光流數(shù)：將前后體掃時刻的兩幀組合反射率圖進行單體分割得到前體掃時刻單體和后體掃時刻單體；然后輸入該兩幀組合反射率圖，并通過光流法計算出光流場數(shù)據(jù)；遍歷光流場數(shù)據(jù)點，對于橫向速度或者縱向速度大于30的光流場數(shù)據(jù)點賦值為0；通過光流場數(shù)據(jù)計算得到前體掃時刻所有像素在后體掃時刻的流入位置；對比前后體掃時刻各單體的范圍，計算出光流數(shù)；

光流數(shù)記為sum(a→b)，用以表示由前體掃時刻單體a流入后體掃時刻單體b的像素數(shù)；用光流數(shù)量化前體掃時刻單體a和后體掃時刻單體b的關(guān)聯(lián)程度，光流數(shù)越高，前體掃時刻單體a和后體掃時刻單體b的相關(guān)聯(lián)程度就越高，反之越低；

結(jié)合前后體掃時刻單體的面積s，綜合考慮前后體掃時刻單體出現(xiàn)分裂或合并的情況，制定如下三個條件，若滿足其一則認為前體掃時刻單體a和后體掃時刻單體b存在關(guān)聯(lián)，

條件一、sum(a→b)>s(a)×50％；

條件二、sum(a→b)>s(b)×50％；

條件三、設(shè)sum(a→b)＝k，對前體掃時刻的l個單體流入后體掃時刻單體b的所有像素求和，記為sum，k/sum>50％，且k/s(b)>33％；

步驟3.根據(jù)步驟1選取的強對流單體特征和步驟2實現(xiàn)的前后體掃時刻的單體追蹤，用一組曲線表示強對流單體特征的變化趨勢。

進一步講，本發(fā)明中，通過結(jié)構(gòu)可視化剖析獲得強對流單體三維圖像及垂直和水平的剖面，步驟如下：

步驟1、構(gòu)建強對流單體三維圖像：首先，讀取雷達基數(shù)據(jù)獲得組合反射率圖，并進行單體分割，插值獲得強對流單體的規(guī)則三維場數(shù)據(jù)；然后，通過移動立方體算法提取等值面；最終，繪制顯示等值面；

步驟2、單體垂直及水平剖面顯示：根據(jù)步驟1中獲得的規(guī)則三維場數(shù)據(jù)，插值獲得垂直和水平剖面。

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

(1)本發(fā)明通過基于光流法的單體跟蹤算法實現(xiàn)單體跟蹤，用一組結(jié)構(gòu)化特征的變化曲線表現(xiàn)強對流單體的特征變化趨勢和變化規(guī)律。

(2)本發(fā)明通過移動立方塊算法實現(xiàn)單體的等值面提取后，采用面繪制的方法繪制單體的三維結(jié)構(gòu)，通過透明度的設(shè)置為用戶展現(xiàn)對流單體更好的三維透視性，通過透明度參數(shù)的設(shè)置為用戶展現(xiàn)特定反射率值域的三維構(gòu)成與分布。

(3)本發(fā)明通過連續(xù)調(diào)整角度，為用戶展現(xiàn)對流單體垂直剖面上的構(gòu)成及強度分布，改進了垂直剖面的人機交互體驗；通過高度參量的設(shè)置，為用戶展現(xiàn)任意高度上對流單體的水平構(gòu)成和強度分布，增加了水平剖面的人機交互體驗。

附圖說明

圖1是對前體掃時刻的組合反射率圖進行單體分割，單體分割結(jié)果如圖所示，矩形框為單體的最小外包矩形，矩形框左上角標號為對該單體的編號；

圖2是對后體掃時刻的組合反射率圖進行單體分割，單體分割結(jié)果如圖所示，矩形框為單體的最小外包矩形，矩形框左上角標號為對該單體的編號；

圖3(a)至圖3(e)是部分單體特征的變化趨勢；橫軸表示時間，縱軸表示特征的量值，虛線箭頭是對當前時刻的標識；其中：圖3(a)是懸垂度變化趨勢，圖3(b)是有效厚度，圖3(c)是加權(quán)核高，圖3(d)是高回波比，圖3(e)是各反射率在高度上的分布示意圖；

圖4表示單體中55dbz的等值面提取結(jié)果；

圖5表示各反射率等值面提取結(jié)果，包括：30dbz，35dbz，40dbz，45dbz，50dbz，55dbz的等值面提取結(jié)果；

圖6表示各反射率的等值面的繪制結(jié)果；

圖7是單體整體三維圖像加上水平剖面的結(jié)果。

具體實施方式

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

本發(fā)明提出的一種強對流單體結(jié)構(gòu)及結(jié)構(gòu)特征可視化時空剖析方法，主要包括獲得強對流單體三維圖像及垂直和水平的剖面、獲得強對流單體特征的變化趨勢。

一、獲得強對流單體特征的變化趨勢的步驟如下：

步驟1、選取強對流單體特征

以雷達基數(shù)據(jù)為原始數(shù)據(jù)，求取獲得有效厚度、高回波比、懸垂度、加權(quán)核高、峰度、組合反射率圖核區(qū)梯度高值、垂直累計液態(tài)水含量和累計液態(tài)水密度；其中，選取有效厚度^[7]、高回波比^[6]等特征表示單體的“高”、“強”特點，用懸垂度表示單體的強懸垂回波結(jié)構(gòu)^[6]，用加權(quán)核高表示水成物固態(tài)率^[6]，用峰度表示回波分布的緊密程度^[7]，用組合反射率圖核區(qū)梯度高值量化核區(qū)的變化程度，也選取常見的冰雹判別指標垂直累計液態(tài)水含量^[9]和累計液態(tài)水密度^[8]；

同時，通過對雷達基數(shù)據(jù)進行插值得到規(guī)則三維場數(shù)據(jù)，根據(jù)規(guī)則三維場數(shù)據(jù)得到各反射率的高度分布，并用一組柱形圖表示；

所選特征的具體定義如下：

(1)有效厚度

定義單體核反射率值大于45dbz的頂高和-20℃層高度之差h-20為有效厚度：

δh＝max{hr>45dbz}-h-20

(2)高回波比

強冰雹反射率圖通常出現(xiàn)大面積的高強反射率區(qū)域。求取反射率值在50dbz以上的區(qū)域面積sr>50占30dbz以上區(qū)域面積sr>30的比例ρ，并將其命名為高回波比；

(3)懸垂度

強對流云體中，強烈的上升氣流位于弱回波區(qū)，其上方出現(xiàn)強懸垂回波，懸掛區(qū)域與弱回波區(qū)域呈正相關(guān)關(guān)系，本文將0℃層之上的弱回波體積定義為回波的懸垂度；

(4)垂直累計液態(tài)水含量，累計液態(tài)水密度與加權(quán)核高

垂直累計液態(tài)水含量vil通過將雷達反射率因子z轉(zhuǎn)換成等價的液態(tài)水，積分得到垂直柱體內(nèi)液態(tài)水總量：

vil＝∫mdhm＝3.44x10^-3z^4/7

設(shè)∫mdh的柱底面積為s，則進一步可以計算出累計液態(tài)水密度vild：

考慮到云中水成物粒子在零度層以上會越來越多的轉(zhuǎn)變成固態(tài)，加入衰減系數(shù)k，形成加權(quán)垂直累積液態(tài)水含量密度：

式中，h0為回波底，h1是回波頂高；k取為高度的函數(shù)，且：

若h0<h<h0℃，則k＝1；

若h0℃<h<h-20℃，則k由1線性衰減到0.1；

若h1>h>h-20℃，則k由0.1線性衰減0。

設(shè)的柱底面積為s，得體積∑(h1-h0)*s，則加權(quán)垂直累積液態(tài)水含量密度為：

只考慮反射率大于45dbz的單體部分，且不區(qū)分核區(qū)內(nèi)的反射率值的差異，則加權(quán)垂直累積液態(tài)水含量密度簡化為下式，并將該特征命名為加權(quán)核高；

(5)峰度

峰度是一種突出直方圖厚重尾部的四階統(tǒng)計量，在一定程度上反映出單體低層反射率圖所呈現(xiàn)的冰雹單體核偏向高梯度一側(cè)、反射率分布緊致的特點：

其中，是均值，s是標準差；

(6)組合反射率圖核區(qū)梯度高值

在組合反射率圖中進行單體分割后，將單體中反射率大于40dbz的區(qū)域定義為單體核區(qū)，按照下述公式求取單體核區(qū)的中心p，

公式(1)中，pr>40dbz為所有反射率大于40dbz的點的坐標，rr>40dbz為對應(yīng)反射率；

以單體核區(qū)的中心p作為圓心作圓覆蓋單體，并將該圓等分為八個扇形，選取三個相鄰扇形使之核區(qū)面積之和最小，并求取這三個扇形的平均梯度g：

公式(2)中，∑g為三個扇形區(qū)域所有像素點的梯度之和，像素點的梯度只考慮大小，不考慮方向，num為像素點總數(shù)，該平均梯度g即為組合反射率圖核區(qū)梯度高值；

步驟2、通過光流法實現(xiàn)前后體掃時刻的單體追蹤：

利用光流法^[10]求得光流數(shù)：將前后體掃時刻的兩幀組合反射率圖進行單體分割^[12]得到前體掃時刻單體和后體掃時刻單體；然后輸入該兩幀組合反射率圖，并通過光流法計算出光流場數(shù)據(jù)；遍歷光流場數(shù)據(jù)點，對于橫向速度或者縱向速度大于30的光流場數(shù)據(jù)點賦值為0；通過光流場數(shù)據(jù)計算得到前體掃時刻所有像素在后體掃時刻的流入位置；對比前后體掃時刻各單體的范圍，計算出光流數(shù)；

光流數(shù)記為sum(a→b)，用以表示由前體掃時刻單體a流入后體掃時刻單體b的像素數(shù)；用光流數(shù)量化前體掃時刻單體a和后體掃時刻單體b的關(guān)聯(lián)程度，光流數(shù)越高，前體掃時刻單體a和后體掃時刻單體b的相關(guān)聯(lián)程度就越高，反之越低；

結(jié)合前后體掃時刻單體的面積s，綜合考慮前后體掃時刻單體出現(xiàn)分裂或合并的情況，制定如下三個條件，若滿足其一則認為前體掃時刻單體a和后體掃時刻單體b存在關(guān)聯(lián)，

條件一、sum(a→b)>s(a)×50％；

條件二、sum(a→b)>s(b)×50％；

條件三、設(shè)sum(a→b)＝k，對前體掃時刻的l個單體流入后體掃時刻單體b的所有像素求和，記為sum，k/sum>50％，且k/s(b)>33％；

消亡單體和新生單體均不會滿足上述任一條件；

若前體掃時刻單體a既不分裂也不合并，且大部分被后體掃時刻單體b繼承，此時對應(yīng)條件一；

若后體掃時刻單體b大部分繼承于前體掃時刻單體a，此時對應(yīng)條件二；

若后體掃時刻單體b大部分來源于單體合并，合并部分的主要貢獻者是前體掃時刻單體a，此時對應(yīng)條件三；

步驟3.根據(jù)步驟1選取的強對流單體特征和步驟2實現(xiàn)的前后體掃時刻的單體追蹤，用一組曲線表示強對流單體特征的變化趨勢。

二、通過結(jié)構(gòu)可視化剖析獲得強對流單體三維圖像及垂直和水平的剖面，步驟如下：

步驟1、構(gòu)建強對流單體三維圖像：首先，讀取雷達基數(shù)據(jù)獲得組合反射率圖，并進行單體分割^[12]，插值獲得強對流單體的規(guī)則三維場數(shù)據(jù)；然后，通過移動立方體算法^[11]提取等值面；最終，繪制顯示等值面；

步驟2、單體縱向及橫向顯示：根據(jù)步驟1中獲得的規(guī)則三維場數(shù)據(jù)，插值獲得垂直和水平剖面。

綜上，本發(fā)明提出的一種強對流單體結(jié)構(gòu)及結(jié)構(gòu)特征可視化時空剖析方法，其中，通過光流法實現(xiàn)前后體掃時刻的單體追蹤，是在實現(xiàn)單體跟蹤的基礎(chǔ)上，選取一組強對流單體的結(jié)構(gòu)化特征，用折線圖的方式體現(xiàn)它們各自的變化趨勢及變化規(guī)律，可為氣象工作者的預(yù)報與科研提供方便。

研究材料：數(shù)據(jù)來源是天津市塘沽雷達站2015年8月18日的雷達數(shù)據(jù)。

圖1和圖2分別是前體掃時刻和后體掃時刻的單體分割結(jié)果，單體標號如圖所示，根據(jù)光流法得到圖1，圖2的光流數(shù)如表1所示，前后體掃時刻單體面積如表2所示：

表1光流數(shù)表

表2單體面積表

根據(jù)上表得到sum(0→0)＝290>s(后體掃時刻0號單體)×50％＝339×50％＝169.5，故后體掃時刻0號單體與前體掃時刻0號單體存在聯(lián)系，同理可判斷后體掃時刻1，2號單體與前體掃時刻滿0號單體存在聯(lián)系，后體掃時刻3號單體來與前體掃時刻1號單體存在聯(lián)系。以上結(jié)果表明前體掃時刻0號單體分裂為后體掃時刻的0，1，2號單體，后體掃時刻3號單體主要繼承了前體掃時刻的1號單體。

在實現(xiàn)單體跟蹤的基礎(chǔ)上，計算每個體掃時刻所有單體的特征，用圖表的形式展現(xiàn)單體的特征變化趨勢，部分特征的變化趨勢如圖3(a)至圖3(e)所示，橫軸表示時間，縱軸表示特征的量值，虛線箭頭是對當前時刻的標識；其中：圖3(a)是懸垂度變化趨勢，圖3(b)是有效厚度，圖3(c)是加權(quán)核高，圖3(d)是高回波比，圖3(e)是各反射率在高度上的分布示意圖；

單體三維重建的實現(xiàn)主要分為以下三個步驟：

(1)在讀取雷達基數(shù)據(jù)→單體分割→插值成規(guī)則三維場數(shù)據(jù)。

(2)通過移動立方體算法提取特定等值面，圖4為55dbz等值面提取結(jié)果，可見其提取結(jié)果是由多個三角面構(gòu)成的網(wǎng)格面。

(3)圖5為30dbz，35dbz，40dbz，45dbz，50dbz，55dbz等值面提取結(jié)果，在提取等值面后，將相應(yīng)的網(wǎng)格數(shù)據(jù)繪制為面，并設(shè)置合適的透明度。由于反射率高的等值面一般被反射率低的等值面包含，故反射率高的等值面的透明度低，反射率低的等值面反之，以此帶來較好的透視性，圖6為最后繪制結(jié)果，圖7是單體三維圖像加上水平剖面的效果圖。

盡管上面結(jié)合附圖對本發(fā)明進行了描述，但是本發(fā)明并不局限于上述的具體實施方式，上述的具體實施方式僅僅是示意性的，而不是限制性的，本領(lǐng)域的普通技術(shù)人員在本發(fā)明的啟示下，在不脫離本發(fā)明宗旨的情況下，還可以做出很多變形，這些均屬于本發(fā)明的保護之內(nèi)。