王 萍，顧 愷，侯謹(jǐn)毅

(天津大學(xué)電氣自動(dòng)化與信息工程學(xué)院，天津 300072)

強(qiáng)對(duì)流天氣目前是指伴隨雷暴現(xiàn)象的對(duì)流性大風(fēng)(風(fēng)速≥17.2 m/s)、冰雹、短時(shí)強(qiáng)降水等．強(qiáng)對(duì)流天氣發(fā)生于中小尺度天氣系統(tǒng)，空間尺度小，其生命史短暫并帶有明顯的突發(fā)性，較短的僅有幾分鐘至一小時(shí)[1]. 雖然強(qiáng)對(duì)流天氣發(fā)生的時(shí)間短，但是破壞力卻非常強(qiáng)，它會(huì)引發(fā)洪澇、滑坡、泥石流等自然災(zāi)害，致使房倒屋毀、植被作物受災(zāi)、電信交通中斷甚至人員傷亡等，因此分析和預(yù)報(bào)強(qiáng)對(duì)流天氣對(duì)國(guó)民經(jīng)濟(jì)、民眾安全有著重要意義．強(qiáng)對(duì)流天氣是由空氣劇烈垂直運(yùn)動(dòng)導(dǎo)致的，以上升和下沉氣流為主．描述出包含上升和下沉氣流的三維流場(chǎng)有助于對(duì)強(qiáng)對(duì)流天氣的分析.

對(duì)強(qiáng)對(duì)流天氣的分析基本都依賴多普勒天氣雷達(dá)數(shù)據(jù)．多普勒天氣雷達(dá)作為一種有效、及時(shí)的觀測(cè)手段，時(shí)間、空間分辨率較高，推動(dòng)了短時(shí)臨近預(yù)報(bào)的發(fā)展，也是研究中小尺度天氣系統(tǒng)的重要資料來(lái)源．隨著全國(guó)天氣雷達(dá)布網(wǎng)的完善，雷達(dá)資料產(chǎn)品在天氣預(yù)報(bào)分析領(lǐng)域扮演著越來(lái)越重要的角色，新一代多普勒天氣雷達(dá)所采集到的是各仰角上的圓錐面上的信息[2]，傳統(tǒng)的強(qiáng)對(duì)流分析方法都是基于各仰角上的雷達(dá)數(shù)據(jù)進(jìn)行的，這對(duì)于由空氣強(qiáng)烈的垂直運(yùn)動(dòng)引發(fā)的強(qiáng)對(duì)流天氣而言，空間垂直維度連續(xù)性信息的欠缺使分析效果受到限制．

目前，國(guó)內(nèi)外一些科研人員主要通過(guò)三維風(fēng)場(chǎng)反演來(lái)達(dá)到描述強(qiáng)對(duì)流三維流場(chǎng)的目的．1975 年，Ray等[3]較早提出雙雷達(dá)風(fēng)場(chǎng)反演的直接合成法ODD(over determined dual-Doppler)技術(shù)，反演出了龍卷云的三維風(fēng)場(chǎng)結(jié)構(gòu)；1998 年，Bousquet 等[4]在ODD 技術(shù)的基礎(chǔ)上又提出了EODD(extended ODD)技術(shù)；2012 年，羅昌榮等[5]將雙雷達(dá)反演方法普遍采用的笛卡爾坐標(biāo)系改為地球坐標(biāo)系，并通過(guò)質(zhì)量連續(xù)方程，以及迭代的方法進(jìn)行估算從而使用雙雷達(dá)反演三維風(fēng)場(chǎng)．總之，通過(guò)三雷達(dá)或雙雷達(dá)的反演效果較好，但僅僅局限于多雷達(dá)探測(cè)范圍重疊的區(qū)域．因此也有些研究者嘗試通過(guò)單雷達(dá)反演三維風(fēng)場(chǎng)：1991年，Sun 等[6]提出一種全四維變分法，以Boussinesq模型為預(yù)報(bào)模型，具有反演三維風(fēng)場(chǎng)和熱力場(chǎng)的能力；莊子波等[7]提出了VAD 方法，但是VAD 方法對(duì)垂直速度的估計(jì)不準(zhǔn)確．總之，若僅使用單雷達(dá)進(jìn)行三維風(fēng)場(chǎng)反演，由于信息有限，反演的結(jié)果精確度有限，且計(jì)算復(fù)雜度較高．

如果有一種較為簡(jiǎn)捷的方法，能夠達(dá)到描述強(qiáng)對(duì)流三維流場(chǎng)的效果，那么不僅可以詳細(xì)、具體地表現(xiàn)出強(qiáng)對(duì)流天氣過(guò)程中各個(gè)時(shí)刻流場(chǎng)中氣流的情況，還有助于預(yù)報(bào)員分析強(qiáng)對(duì)流的發(fā)展趨勢(shì)，為強(qiáng)對(duì)流天氣的預(yù)報(bào)提供依據(jù)．

本文基于多普勒天氣雷達(dá)基數(shù)據(jù)，對(duì)雷達(dá)基數(shù)據(jù)進(jìn)行雙線性插值的預(yù)處理．設(shè)置反射率因子閾值，獲得反射率因子大于等于閾值的三維格點(diǎn)數(shù)據(jù)，對(duì)三維格點(diǎn)數(shù)據(jù)通過(guò) Marching Cubes 方法進(jìn)行曲面重建．基于相鄰體掃時(shí)刻同一數(shù)據(jù)點(diǎn)的反射率因子基本不變的假設(shè)，設(shè)計(jì)了一種非剛性的、基于特征融合的曲面匹配方法，對(duì)相鄰體掃時(shí)刻的曲面間進(jìn)行曲面匹配．通過(guò)曲面匹配可獲得三維流場(chǎng)中各數(shù)據(jù)點(diǎn)的位移向量，從而近似描述出強(qiáng)對(duì)流天氣的三維流場(chǎng)．最后，通過(guò)模擬數(shù)據(jù)進(jìn)行曲面匹配方法的效果測(cè)試．

1 數(shù)據(jù)準(zhǔn)備

新一代多普勒天氣雷達(dá)掃描所采集到的是圓錐面上的信息[2]，掃描方式如圖1 所示．對(duì)于多普勒雷達(dá)的VCP21 體掃模式[8]，雷達(dá)在6 min 間隔下掃描9個(gè)仰角，獲取9 個(gè)仰角上的圓錐面數(shù)據(jù)．其中9 個(gè)仰角具體為 0.5°、1.5°、2.4°、3.4°、4.3°、6.0°、9.9°、14.6°、19.5°．實(shí)際上，預(yù)報(bào)員分析所依賴的雷達(dá)圖(如圖2 所示)就是這些仰角上圓錐面的俯視平面圖．

圖1 多普勒雷達(dá)掃描示意Fig.1 Schematic of Doppler radar scanning

圖2 雷達(dá)反射率因子回波Fig.2 Echo of radar reflectivity factor

對(duì)于同一組體掃數(shù)據(jù)的9 張錐面圖，對(duì)應(yīng)二維水平面上的點(diǎn)，選取不同高度9 個(gè)數(shù)據(jù)中的最大反射率因子合成一張組合反射率圖．本文采用的單體分割算法[9]首先根據(jù)對(duì)流單體從內(nèi)到外反射率因子強(qiáng)度逐級(jí)減弱的分布規(guī)律，通過(guò)區(qū)域生長(zhǎng)法得到一系列獨(dú)立的對(duì)流單體核區(qū)，再使用膨脹避讓算法[9]，自適應(yīng)地在雷達(dá)組合反射率圖上分割出對(duì)流單體．通過(guò)這種方法分割的單體，不受單體尺度、雷達(dá)回波強(qiáng)度的限制，也較好地解決了外圍相互交織的多單體對(duì)象的自適應(yīng)分割問(wèn)題．

為了獲得對(duì)流單體的三維描述，本文利用9 個(gè)仰角的反射率因子數(shù)據(jù)經(jīng)過(guò)雙線性插值[10]得到高度分辨率為0.25 km、水平分辨率為1 km×1 km 的70 層平行于水平面的反射率因子數(shù)據(jù)，雷達(dá)射線與對(duì)流單體三維格點(diǎn)生成數(shù)據(jù)的關(guān)系如圖3 所示．最終得到70 層三維格點(diǎn)數(shù)據(jù)如圖4 所示，它們可以表現(xiàn)出對(duì)流單體三維區(qū)域幾乎每一處的信息，作為本文研究的輸入數(shù)據(jù)．

圖3 70層水平格點(diǎn)數(shù)據(jù)與雷達(dá)射線的關(guān)系示意Fig.3 Schematic of the relationship between 70 layers of horizontal grid data and radar rays

圖4 三維格點(diǎn)數(shù)據(jù)示意Fig.4 Schematic of three-dimensional grid point data

2 曲面重建

2.1 三維點(diǎn)云的獲取與去噪

從內(nèi)到外反射率因子逐層降低的分布特點(diǎn)是對(duì)流單體的共性．因此，對(duì)于前一節(jié)獲取的70 層水平格點(diǎn)數(shù)據(jù)，設(shè)定一個(gè)反射率因子閾值 Rth后，所有反射率因子大于等于該閾值的三維格點(diǎn)組成了一個(gè)實(shí)心的點(diǎn)云核．

對(duì)于該三維點(diǎn)云核，本文采取半徑濾波方法[11]濾除點(diǎn)云中的噪聲．具體而言，對(duì)于三維點(diǎn)云數(shù)據(jù)，如果一個(gè)點(diǎn)在設(shè)定搜索半徑r范圍內(nèi)的鄰近點(diǎn)數(shù)量小于給定的閾值t，則判斷該點(diǎn)為遠(yuǎn)離主體點(diǎn)云的稀疏離群點(diǎn)，做刪除處理．經(jīng)過(guò)測(cè)試與統(tǒng)計(jì)分析，設(shè)置r為3 個(gè)單位長(zhǎng)度，閾值t為3．

2.2 通過(guò)Marching Cubes方法重建曲面

拓?fù)浣Y(jié)構(gòu)對(duì)于點(diǎn)云的匹配、形狀分析等應(yīng)用具有重要作用．而上文提取的空間點(diǎn)云僅僅表示物體的幾何形狀，無(wú)法表達(dá)內(nèi)部拓?fù)浣Y(jié)構(gòu)．本文采用一種常見(jiàn)的曲面重建方法——Marching Cubes 方法[12]．這種方法本質(zhì)上是從一個(gè)三維空間數(shù)據(jù)場(chǎng)中抽取出一個(gè)等值面．

Marching Cubes 方法采用分治的策略，基本思想是把三維點(diǎn)云中相鄰的8 個(gè)體元看作一個(gè)六面立方體的8 個(gè)頂點(diǎn)．對(duì)于圖4 中H×W×70 的三維格點(diǎn)數(shù)據(jù)而言，六面立方體個(gè)數(shù)為(H－1)×(W－1)×69.把每個(gè)六面立方體的8 個(gè)頂點(diǎn)的值與設(shè)定的等值面的閾值進(jìn)行比較，而每個(gè)頂點(diǎn)的值與閾值相比有大于閾值或不大于閾值兩種情況，因此一共有28＝256 種情況，可通過(guò)建立一個(gè)長(zhǎng)度為256 的表，供算法查詢調(diào)用．

圖5 對(duì)其中15 種典型的情況給出了說(shuō)明. 對(duì)于每個(gè)立方體，標(biāo)記有實(shí)心點(diǎn)的頂點(diǎn)表示該點(diǎn)的值大于設(shè)定的閾值，其余頂點(diǎn)數(shù)值未達(dá)到閾值. 以圖5 中的2 號(hào)子圖為例，頂點(diǎn)A是該立方體8 個(gè)頂點(diǎn)中唯一大于所設(shè)閾值的點(diǎn)，與其相鄰的頂點(diǎn)有B、C、D．在AB、AC、AD這3 條線段上，分別根據(jù)相應(yīng)線段的兩個(gè)頂點(diǎn)的數(shù)值做線性插值，找到每條線段上數(shù)值等于閾值的點(diǎn)的位置．對(duì)于這3 條線段，設(shè)線段的端點(diǎn)坐標(biāo)為1P、2P ，它們的值分別為R1和R2．等值面的閾值為R，P為待求的坐標(biāo)，其計(jì)算式為

本文曲面重建的具體步驟如下．

步驟1建立長(zhǎng)度為256 的索引表．

步驟2對(duì)三維格點(diǎn)數(shù)據(jù)，取出與其相鄰的7 個(gè)數(shù)據(jù)格點(diǎn)，構(gòu)成一個(gè)正六面體，然后根據(jù)這8 個(gè)數(shù)據(jù)點(diǎn)與設(shè)定的等值面的閾值的大小關(guān)系，得出對(duì)應(yīng)的索引值．

步驟3根據(jù)索引值搜索索引表，得到需要進(jìn)行插值的邊，進(jìn)而得到三角面．

步驟4遍歷所有正六面體，重復(fù)步驟2 和步驟3，進(jìn)而得到所有三角面的輸出．

圖5 Marching Cubes 15種典型情況示意Fig.5 Schematic of 15 typical Marching Cubes cases

通過(guò)以上4 步，把這些生成的面相連就得到了最終的等值面．這種方法就好像移動(dòng)立方塊一樣，Marching Cubes 算法因此得名．總的來(lái)說(shuō)，Marching Cubes 算法原理簡(jiǎn)單、效果較好．圖 6 是通過(guò)Marching Cubes 算法對(duì)對(duì)流單體點(diǎn)云數(shù)據(jù)重建表面的一個(gè)示例．

圖6 Marching Cubes重建效果Fig.6 Reconstructed rendering of Marching Cubes

3 非剛性曲面匹配方法設(shè)計(jì)

曲面匹配是逆向工程的關(guān)鍵技術(shù)之一．一般曲面匹配是通過(guò)求解一個(gè)剛性坐標(biāo)變換矩陣(對(duì)于三維空間，空間變換矩陣是4×4 的)，將不同三維數(shù)據(jù)點(diǎn)集統(tǒng)一到一個(gè)坐標(biāo)系下．最為常見(jiàn)的曲面匹配方法為基于最小二乘原理進(jìn)行多次迭代的ICP(iterative closest point)方法[13]，該類方法的最終目標(biāo)是求解出兩個(gè)曲面間最佳的剛性坐標(biāo)變換矩陣．因此主要適用于兩個(gè)外形完全一致的曲面，如不同攝像頭視角下的同一物體．但是本文涉及的前后時(shí)刻反射率因子大于設(shè)定閾值的三維格點(diǎn)重建的一組曲面形狀雖然相似但并非完全一致，為此本文設(shè)計(jì)了一種精確至點(diǎn)到點(diǎn)、非剛性的曲面匹配方法．

3.1 曲面點(diǎn)特征向量設(shè)計(jì)

三維曲面上各點(diǎn)的法向量指的是垂直于該點(diǎn)切平面的向量，通過(guò)三維曲面上某點(diǎn)的法向量，可以獲得該點(diǎn)處局部曲面的朝向信息．另外，曲率是曲面彎曲程度的一種度量，是一種重要的曲面局部特征，通過(guò)三維曲面上某點(diǎn)的曲率，可以得到該點(diǎn)處局部曲面的凹凸信息．綜合曲面上一點(diǎn)的法向量與曲率兩種信息，就可以相當(dāng)程度上描述出該點(diǎn)處的曲面特征．本文具體的法向量與曲率的求法如下所述．

1) 曲面點(diǎn)的法向量

面模型是由一定數(shù)目的面片來(lái)逼近的，面片越多則面模型越精細(xì)．對(duì)于各面片的法向量，通過(guò)組成各面片的任意兩條邊的叉乘向量并單位化來(lái)表示．這里以三角形面片為例，給出具體計(jì)算方法．設(shè)圖7 中的F是一個(gè)三角形面片，P0( x0, y0, z0)、P1( x1, y1, z1)、P2( x2, y2, z2)是這三角形面片的3 個(gè)頂點(diǎn)．G是F所在平面．n為位于P0處的該三角形面片的法向量，則有

圖7 三角形面片法向量求解示意Fig.7 Schematic of normal vector solution for triangular patch

對(duì)于曲面任意點(diǎn)P的法向量，設(shè)為該點(diǎn)周圍所有面片的法向量的平均．假設(shè)NP為該點(diǎn)周圍面片數(shù)，in 為該點(diǎn)周圍第i個(gè)面片的法向量，nP為P點(diǎn)的法向量，則有

向量示意如圖8 所示．

2) 曲面點(diǎn)的曲率

曲面上某一點(diǎn)的曲率，是基于該點(diǎn)法向量進(jìn)行求解的．經(jīng)過(guò)曲面上點(diǎn)P的法線(法向量)的某一平面與曲面相交可得到一條曲線，這條曲線在點(diǎn)P處的曲率即為曲面上的點(diǎn)P的曲率，如圖9 所示．經(jīng)過(guò)頂點(diǎn)P的法線的平面可以有無(wú)數(shù)個(gè)(可任意旋轉(zhuǎn))，這也意味著相交的曲線也有無(wú)數(shù)條，而每條曲線在頂點(diǎn)P都有一個(gè)曲率．設(shè)所有這些曲率中最大的曲率為κH，最小的曲率為κL．其他常見(jiàn)曲率如高斯曲率、平均曲率基于最大曲率和最小曲率求得．最大曲率與最小曲率求法分別為

圖8 通過(guò)P點(diǎn)周圍面片法向量求解該點(diǎn)法向量示意Fig.8 Schematic for calculating the normal vector of point P by the normal vectors of the surrounding sur-faces

圖9 曲面曲率示意Fig.9 Schematic of surface curvature

綜上，曲面上某一點(diǎn)的法向量nP=( nx, ny, nz)代表該處的朝向，曲率又代表了該處的凹凸程度．本文把法向量與最大、最小曲率進(jìn)行融合[14-15]，形成一個(gè)5 維的特征向量，以表示曲面上某一點(diǎn)的形狀信息，即

考慮到這個(gè)特征向量I 的5 元素量綱并不一致，為了消除各參數(shù)量綱的影響，本文對(duì)法向量進(jìn)行單位化處理，對(duì)最大、最小曲率進(jìn)行歸一化處理．具體公式分別為

式中：為所有特征向量中最大曲率κH的平均值；為所有特征向量中最大曲率κH的最大值；為所有特征向量中最大曲率κH的最小值；為所有特征向量的最小曲率κL的平均值；為所有特征向量中最小曲率κL的最大值；為所有特征向量中最小曲率κL的最小值．通過(guò)式(7)～(11)，原特征向量I變?yōu)樾碌奶卣飨蛄科涓髟氐哪Ｖ稻怀^(guò)1，消除了量綱的影響．

3.2 匹配區(qū)域的限定與“相對(duì)坐標(biāo)向量”

在使用特征向量進(jìn)行曲面匹配時(shí)，還需要限定匹配區(qū)域．例如，對(duì)于單體曲面東南角的某點(diǎn)，它不應(yīng)該與下一時(shí)刻單體曲面西北角的某點(diǎn)匹配，即便它們的5 維特征向量非常相似．所以對(duì)于曲面上的每一點(diǎn)，除了該點(diǎn)處曲面的形狀和朝向信息外，還應(yīng)考慮該點(diǎn)在曲面的位置信息來(lái)限定匹配區(qū)域．對(duì)于每個(gè)點(diǎn)，它們都有一個(gè)3 維空間的坐標(biāo)，可以通過(guò)這個(gè)3維坐標(biāo)來(lái)直接表示出一個(gè)點(diǎn)的位置信息．然而，對(duì)流單體內(nèi)部的氣流不僅有上升與下沉的運(yùn)動(dòng)，還有幅度更大的東南西北方向的水平移動(dòng)，對(duì)于同一降水目標(biāo)物粒子，它在相鄰體掃時(shí)刻下的兩個(gè)位置的3 維坐標(biāo)可能相差很大，因此直接用3 維空間的坐標(biāo)來(lái)表示位置信息并不能很好地達(dá)到匹配區(qū)域限定的目的．

為了解決這個(gè)問(wèn)題，本文提出了一種新的位置信息表達(dá)方法，以達(dá)到表示某點(diǎn)在3 維曲面上的“相對(duì)位置”的效果．假設(shè)反射率因子大于等于設(shè)定閾值的3 維格點(diǎn)的外邊界曲面包含的外邊界點(diǎn)有N 個(gè)．在如圖10 所示的以單體框?yàn)榈酌?、垂直地面向上的射線作為z 軸建立的坐標(biāo)系下，設(shè)曲面上的點(diǎn)P的坐標(biāo)為( xP, yP, zP)．對(duì)于這N 個(gè)外邊界點(diǎn)，坐標(biāo)x 小于 xP的點(diǎn)的個(gè)數(shù)設(shè)為Nx，坐標(biāo)y 小于 yP的點(diǎn)的個(gè)數(shù)設(shè)為Ny，坐標(biāo)z 小于 zP的點(diǎn)的個(gè)數(shù)設(shè)為Nz．進(jìn)而用3 維向量作為一種新的坐標(biāo)信息表達(dá)方式，并稱之為頂點(diǎn)P在重建曲面上的“相對(duì)坐標(biāo)向量”．以曲面最高點(diǎn)為例(假設(shè)最高點(diǎn)唯一)，那么它的“相對(duì)坐標(biāo)向量”中的第 3 維的值為(N－1)/N．“相對(duì)坐標(biāo)向量”很好地表現(xiàn)了一個(gè)點(diǎn)在曲面上的相對(duì)位置．對(duì)于空氣中同一個(gè)降水物粒子而言，雖然它在相鄰體掃時(shí)刻間隔時(shí)間內(nèi)位移的絕對(duì)距離可能很大，但它在重建表面上的相對(duì)位置變化應(yīng)該較小，同時(shí)，“相對(duì)坐標(biāo)向量”變化也會(huì)較?。ㄟ^(guò)多次

圖10 三維笛卡爾坐標(biāo)系與曲面上的點(diǎn)Fig.10 Three-dimensional Cartesian coordinate system and points on surface

3.3 相似性度量的選擇

對(duì)于不同點(diǎn)的特征向量，評(píng)價(jià)它們的相似性程度時(shí)，需要用到相似性度量．常見(jiàn)的相似性度量[16]有歐氏距離、曼哈頓距離、切比雪夫距離、馬氏距離、余弦相似度等等．本文選擇歐式距離．

考慮到對(duì)于前一時(shí)刻曲面上的點(diǎn)P，后一時(shí)刻曲面上可能有多個(gè)點(diǎn)的5 維特征向量與P的5 維特征向量間的歐式距離都很小甚至相同的情況，構(gòu)建附加匹配條件如下：設(shè)前一體掃時(shí)刻曲面上的點(diǎn)為P，設(shè)后一時(shí)刻曲面上與點(diǎn)P的5 維特征向量之間的歐氏距離最小的點(diǎn)為Q1，歐式距離次近的點(diǎn)為Q2，即當(dāng)最近歐式距離與次近歐式距離的比值小于設(shè)定閾值時(shí)，才認(rèn)為這兩點(diǎn)組成點(diǎn)對(duì)是匹配的．這里的閾值設(shè)為0.5(即

3.4 方法流程

綜合第3.1、3.2、3.3 節(jié)，本文提出的曲面匹配方法的具體步驟可以總結(jié)如下．

步驟1前一時(shí)刻的曲面選取一點(diǎn)P并作標(biāo)記，求它的5 維特征向量以及“相對(duì)坐標(biāo)向量”．

步驟2對(duì)后一時(shí)刻的曲面，限定匹配區(qū)域?yàn)椤跋鄬?duì)坐標(biāo)向量”與P相差10%以內(nèi)的區(qū)域．

步驟3在后一時(shí)刻的曲面的限定匹配區(qū)域下，通過(guò)遍歷找出與P點(diǎn)的5 維特征向量間歐式距離最近的點(diǎn)1Q以及次近點(diǎn)2Q ．

步驟4判斷歐式最近距離與次近距離的比值是否小于設(shè)定閾值(本文取為0.5)，如果小于0.5 則把P與Q設(shè)為“匹配點(diǎn)對(duì)”，否則不建立“匹配點(diǎn)對(duì)”. 然后返回步驟1，步驟1 選取點(diǎn)時(shí)跳過(guò)所有已標(biāo)記的點(diǎn).

不停地迭代以上步驟，直至遍歷完前一時(shí)刻曲面上的所有點(diǎn)．

通過(guò)上述步驟就獲得了前一時(shí)刻點(diǎn)云重建曲面上的點(diǎn)在下一時(shí)刻對(duì)應(yīng)的匹配點(diǎn)，前后時(shí)刻點(diǎn)之間的坐標(biāo)差就被認(rèn)定是大氣中這些降水物粒子的位移向量，用以反映氣流場(chǎng)中的氣流流向．

3.5 三維流場(chǎng)的描述

通過(guò)本文提出的曲面匹配方法，可以得到前一時(shí)刻曲面的所有外邊界點(diǎn)在下一時(shí)刻的匹配點(diǎn)，前后時(shí)刻點(diǎn)之間的坐標(biāo)差就是大氣中粒子的位移向量．將這些位移向量通過(guò)可視化函數(shù)庫(kù)VTK 進(jìn)行3 維可視化，如圖11 所示，其中每個(gè)箭頭的起點(diǎn)代表了前一體掃時(shí)刻的曲面上的某一點(diǎn)P的位置，箭頭的末端(尖頭處)代表后一體掃時(shí)刻的曲面上與點(diǎn)P組成“匹配點(diǎn)對(duì)”的點(diǎn)Q．箭頭代表了在一個(gè)體掃時(shí)刻間隔(6 min)內(nèi)點(diǎn)P到點(diǎn)Q的相對(duì)位移．由此就獲得了這個(gè)曲面上所有點(diǎn)在時(shí)間間隔6 min 內(nèi)的位移，并獲得強(qiáng)對(duì)流三維流場(chǎng)的近似描述．

圖11 位移向量三維可視化示意Fig.11 Schematic of three-dimensional visualization of displacement vectors

具體地，可令上升箭頭為白色，下沉箭頭為黑色，這樣可以很清晰地、定性地觀察出流場(chǎng)各處上升與下沉氣流的情況，如圖12 所示．通過(guò)對(duì)這些箭頭(向量)進(jìn)行求均值、計(jì)數(shù)等數(shù)據(jù)處理，也可以定量地表示出當(dāng)前時(shí)刻流場(chǎng)整體的上升與下沉氣流的對(duì)比．

另外，在雷達(dá)徑向速度圖中，正速度指的是在雷達(dá)波束的徑向上遠(yuǎn)離雷達(dá)(即徑向速度為正)，一般用紅色等暖色表示；負(fù)速度指的是在雷達(dá)波束的徑向上靠近雷達(dá)(即徑向速度為負(fù))，一般用綠色等冷色表示．相應(yīng)地，也可以通過(guò)紅色箭頭和綠色箭頭(對(duì)應(yīng)雷達(dá)徑向速度圖的正速度與負(fù)速度)來(lái)分別代表徑向上遠(yuǎn)離雷達(dá)和靠近雷達(dá)的位移向量，如圖 13 所示．而位移向量是遠(yuǎn)離或靠近雷達(dá)是通過(guò)位移向量與該點(diǎn)到雷達(dá)的連接向量(即雷達(dá)掃描徑線)的點(diǎn)乘的正負(fù)判斷．

圖12 上升、下沉位移向量示意Fig.12 Schematic of rising and sinking displacement vectors

圖13 徑向正、負(fù)位移向量示意Fig.13 Schematic of radial positive and negative displacement vectors

4 曲面匹配效果測(cè)試

在實(shí)際應(yīng)用中，對(duì)于兩個(gè)相似但不完全一致的曲面，無(wú)法得到客觀上“正確”的匹配結(jié)果．而如果是對(duì)于兩個(gè)形狀一致但角度不同的兩個(gè)曲面(比如第2 個(gè)曲面是由第1 個(gè)曲面進(jìn)行旋轉(zhuǎn)后得到的)，可以知道這兩個(gè)曲面上點(diǎn)的“正確”的匹配結(jié)果，很容易對(duì)算法的結(jié)果進(jìn)行評(píng)價(jià)．比如把通過(guò)本文算法得到的正確匹配數(shù)除以全部匹配數(shù)，得到的百分比值就可以作為一個(gè)很好的評(píng)價(jià)曲面匹配效果標(biāo)準(zhǔn)，稱之為“匹配正確率”．因此本文通過(guò)模擬數(shù)據(jù)進(jìn)行測(cè)試，對(duì)于一個(gè)已知的曲面S，對(duì)其施加平移、旋轉(zhuǎn)、縮放、仿射變換后，獲得新的曲面T．運(yùn)用本文方法建立S 與T 的匹配，并統(tǒng)計(jì)“匹配正確率”，給出對(duì)本文匹配方法的效果評(píng)測(cè)．

4.1 三維平移變換

將曲面上一點(diǎn)P( x0, y0, z0)在x、y、z 3 個(gè)坐標(biāo)軸方向上分別平移 a、b 和 c，得到一個(gè)新的點(diǎn)P′ ( x ′, y ′, z′)．設(shè)R為點(diǎn)P到點(diǎn) P ′的轉(zhuǎn)移矩陣，并設(shè)P′=RP，那么4×4 的轉(zhuǎn)移矩陣R為

4.2 三維旋轉(zhuǎn)變換

三維空間下的旋轉(zhuǎn)變換為繞某一軸進(jìn)行旋轉(zhuǎn)，最基本的旋轉(zhuǎn)變換為繞坐標(biāo)軸x、y、z 的旋轉(zhuǎn)．本文僅通過(guò)基本的繞坐標(biāo)軸的旋轉(zhuǎn)變換來(lái)測(cè)試曲面匹配效果．設(shè)曲面上原本的點(diǎn)為P( x0, y0, z0)，旋轉(zhuǎn)變換后得到的點(diǎn)為P′ ( x′ , y ′, z′)，設(shè)P′=RP，繞x、y、z 3 個(gè)坐標(biāo)軸旋轉(zhuǎn)α角的4×4 的變換矩陣分別為Rx、Ry、Rz，如表1 所示．

表1 三維基本旋轉(zhuǎn)變換Tab.1 Three-dimensional basic rotation transform

4.3 三維縮放變換

對(duì)于曲面而言，進(jìn)行整體的縮小或放大，設(shè)縮小或放大的因數(shù)為k．設(shè)曲面上原本的點(diǎn)為( x0, y0, z0)，縮放變換后得到的點(diǎn)為P′( x′ , y ′, z′)，設(shè)R為點(diǎn)P到點(diǎn)P′的轉(zhuǎn)移矩陣，并設(shè)P ′=RP ．那么轉(zhuǎn)移矩陣R的形式為

4.4 三維仿射變換

仿射變換指的是一次非均勻的(即不是縮放變換)線性變換加上一個(gè)平移變換．通過(guò)仿射變換，球可以轉(zhuǎn)換為橢球．設(shè)P′=RP，三維空間下的仿射變換轉(zhuǎn)移矩陣的形式可以表示為

式中 tx、ty、tz分別為平移變換在x、y、z 3 個(gè)坐標(biāo)軸方向上分別平移的距離．

最終，對(duì)于各種3 維變換后分別得到的各種曲面Ti與原曲面S 進(jìn)行匹配，各選取10 組樣例，平均匹配正確率如表2 所示．從表2 中可知對(duì)于平移變換與縮放變換(表2 的前3 行)，匹配平均正確率分別達(dá)到了99.6%、92.0%、89.7%，這是由于這些變換下，曲面的改變相對(duì)較小，并且均為剛性變化．而對(duì)于旋轉(zhuǎn)變換(表2 的第4 行)，匹配平均正確率為77.6%，匹配正確率較低的原因主要是對(duì)于原曲面S 中的一點(diǎn)，經(jīng)過(guò)旋轉(zhuǎn)變換后它的“相對(duì)坐標(biāo)向量”改變較大，導(dǎo)致與之對(duì)應(yīng)的變換后曲面T 上的點(diǎn)可能并不在本文曲面匹配方法限定的匹配區(qū)域中，這是本文曲面匹配方法的一個(gè)不足之處．對(duì)于仿射變換(表2 的第5 行)，匹配平均正確率為76.3%，仿射變換與其他變換相比，形變更大，因此匹配平均正確率最低．最后，對(duì)這5 種變換的正確率求平均，綜合平均匹配正確率為87.0%．

表2 曲面匹配平均正確率Tab.2 Average accuracies of surface matching cases

5 結(jié) 語(yǔ)

本文為了描述強(qiáng)對(duì)流天氣的三維流場(chǎng)，在相鄰體掃時(shí)刻同一降水物粒子的反射率因子基本不變的假設(shè)下，設(shè)定反射率因子閾值，對(duì)相鄰體掃時(shí)刻反射率因子大于等于閾值的點(diǎn)云構(gòu)建出一組相似但不完全一致的曲面．本文提出了一種新的基于法向量與曲率特征融合、點(diǎn)到點(diǎn)的非剛性曲面匹配的方法，并且設(shè)計(jì)了“相對(duì)坐標(biāo)向量”進(jìn)行匹配區(qū)域的限制，設(shè)置附加條件解決一點(diǎn)對(duì)多點(diǎn)匹配的不確定性問(wèn)題．最后，通過(guò)模擬數(shù)據(jù)進(jìn)行測(cè)試，經(jīng)測(cè)試，本文的曲面匹配方法綜合平均正確率達(dá)到了87.0%．借助雷達(dá)反射率因子數(shù)據(jù)的插值三維格點(diǎn)數(shù)據(jù)，通過(guò)本文的曲面匹配方法可獲得三維流場(chǎng)中各數(shù)據(jù)點(diǎn)的位移向量，從而近似描述出強(qiáng)對(duì)流天氣的三維流場(chǎng)．

當(dāng)然，本文方法尚存在一些不足之處．本文使用的雙線性插值方法，隨著數(shù)據(jù)點(diǎn)到雷達(dá)間距離的增加，雷達(dá)射線變得更稀疏，插值效果也會(huì)受到影響，因此本文方法更適用于距離雷達(dá)較近的區(qū)域，比如水平距離100 km 內(nèi)．另外，本文方法的計(jì)算量較大，可以通過(guò)增大反射率因子閾值 Rth或提升硬件配置實(shí)現(xiàn)縮減運(yùn)算時(shí)間的效果，如何通過(guò)改進(jìn)曲面匹配方法以提升運(yùn)算效率值得進(jìn)一步的研究．