王 興， 錢代麗， 朱 彬，， 安俊琳， 呂晶晶

（南京信息工程大學(xué) a.大氣科學(xué)與環(huán)境氣象國家級(jí)實(shí)驗(yàn)教學(xué)示范中心，b.大氣科學(xué)與氣象信息國家級(jí)虛擬仿真實(shí)驗(yàn)教學(xué)中心，南京 210044）

0 引 言

曲面重建是對離散的空間數(shù)據(jù)進(jìn)行三維渲染的主要技術(shù)手段，其基本思想是按照一組給定的閾值從三維點(diǎn)云數(shù)據(jù)中提取相應(yīng)的等值面，利用曲面渲染算法生成一組由不同顏色相互疊加的等值面［5-6］。曲面重建的算法很多，主要可分為兩類：一是基于斷層輪廓線的表面重建，即將三維數(shù)據(jù)分割成若干個(gè)二維切面，分別計(jì)算各切面的等值線，再將這些等值線連接成等值曲面；二是基于體素的等值面重建，即對三維點(diǎn)云構(gòu)成的體數(shù)據(jù)場中的各個(gè)體素加以分析，提取等值面信息，再將這些等值面連接成等值曲面［7-8］。其中，移動(dòng)立方體法（Marching Cubes，MC）［9］是沿用至今且頗具影響力的體素等值面重建算法之一。

曲面重建技術(shù)在醫(yī)療、建筑等工程領(lǐng)域的研究最為深入，近年來，也有一些研究將該技術(shù)應(yīng)用于天氣雷達(dá)、衛(wèi)星和數(shù)值模式等資料，實(shí)現(xiàn)氣象數(shù)據(jù)的三維可視化交互［10-11］。由于雷達(dá)探測數(shù)據(jù)的空間分布極不規(guī)則，距離雷達(dá)近端與遠(yuǎn)端的相鄰離散點(diǎn)之間，間隔距離相差超過100倍，在重構(gòu)均勻三維格點(diǎn)場雷達(dá)回波時(shí)，傳統(tǒng)空間插值算法并不適用。雷達(dá)遠(yuǎn)端的插值點(diǎn)不僅難以還原真實(shí)的探測值，還極大地增加了點(diǎn)云的密度，造成面片的計(jì)算規(guī)模陡增，網(wǎng)絡(luò)結(jié)構(gòu)復(fù)雜化，大幅增加了三維渲染的時(shí)長，影響三維可視化的實(shí)時(shí)性和可交互性。

為解決上述技術(shù)難題，針對多普勒天氣雷達(dá)探測方式及其數(shù)據(jù)特征，在移動(dòng)立方體算法基礎(chǔ)上，提出一系列改進(jìn)，并對改進(jìn)前、后的渲染效果做出比對分析。

1 基本原理與技術(shù)路線

1.1 移動(dòng)立方體（MC）算法

MC算法最早是Lorensen等［9］提出，用于解決三維物體的面繪制問題。其本質(zhì)是將三維坐標(biāo)系下的物體分解為體素集合，每個(gè)體素包含8個(gè)空間上緊相鄰的頂點(diǎn)，構(gòu)成一個(gè)正立方體，將該正立方體視為基本單位。曲面重建時(shí)，先設(shè)定一組大小不同的閾值，通過計(jì)算每個(gè)體素內(nèi)的梯度變化，并與各個(gè)閾值分別進(jìn)行比較，找出所有被閾值穿過的正立方體，再利用插值算法計(jì)算出這些表面（即等值面），再通過繪制一系列閾值的等值面，實(shí)現(xiàn)三維矢量圖像的繪制。

受制于20世紀(jì)80年代有限的計(jì)算機(jī)圖形處理能力，MC算法巧妙地利用三角面片近似表示待求取的等值面，極大地提升了計(jì)算效率。由于等值面與體素邊界面的交線是雙曲線，在特定情況下對相同的等值點(diǎn)可以采取不同的連接方式，這種近似表示的算法無法保證三角片所構(gòu)成的等值面的拓?fù)湟恢滦?，造成很多三角面片銜接處出現(xiàn)孔隙和不連續(xù)現(xiàn)象。為克服這一問題，很多學(xué)者對MC算法提出了改進(jìn)，有的是利用插值和雙曲漸近線將體素進(jìn)一步拆分［12-13］，有的是將MC的正立方體進(jìn)一步簡化，采用四面體剖分［14］等。

1.2 MC應(yīng)用于雷達(dá)回波渲染的技術(shù)缺陷

氣象雷達(dá)探測得到的數(shù)據(jù)，其空間分布是一個(gè)不規(guī)則的三維點(diǎn)云結(jié)構(gòu)，隨著探測距離的增加，點(diǎn)云愈發(fā)稀疏，距離雷達(dá)近端與遠(yuǎn)端的相鄰離散點(diǎn)之間，間隔距離相差可達(dá)100倍以上。在重構(gòu)均勻三維各點(diǎn)場雷達(dá)回波時(shí)，諸如雙線性插值、最近鄰插值、三次內(nèi)插等傳統(tǒng)空間插值算法并不適用。

宴會(huì)中，馬老當(dāng)場宣布，這次活動(dòng)的成功舉辦，在座的各位功不可沒。公司將論功行賞，給予在座的各位同仁以物質(zhì)獎(jiǎng)勵(lì)，對公司的高管層將另行獎(jiǎng)賞。獎(jiǎng)勵(lì)田卓小姐、高潮先生赴巴厘島六日游，馮可兒小姐本來也在去巴厘島之列，但考慮到公司不可一日無主，就留下來代理田卓行使管理權(quán)，但公司會(huì)額外補(bǔ)償馮可兒小姐不能成行的損失。

盡管MC利用三角面片擬合曲面的算法設(shè)計(jì)非常巧妙，在曲面渲染時(shí)也表現(xiàn)出不錯(cuò)的性能，但對于面片的生成，卻需要消耗“過多”的內(nèi)存和計(jì)算時(shí)長。這主要是因?yàn)镸C算法需要遍歷三維點(diǎn)云數(shù)據(jù)中的所有體素，而與等值面相交的體素?cái)?shù)量僅占所有體素的1/10甚至更低。雷達(dá)探測的空間范圍大（數(shù)百公里），數(shù)據(jù)密度高，一次計(jì)算需要處理數(shù)百萬個(gè)體素，而近9成的時(shí)間被用于計(jì)算這些無效體素。雷達(dá)遠(yuǎn)端的插值點(diǎn)不僅難以還原真實(shí)的探測值，還極大地增加了點(diǎn)云的密度，使得體素產(chǎn)生的面片面積過小，造成面片的計(jì)算規(guī)模陡增，網(wǎng)絡(luò)結(jié)構(gòu)復(fù)雜化，大幅增加了三維渲染的時(shí)長，影響了三維回波渲染的實(shí)時(shí)性和可交互性。

1.3 雷達(dá)三維回波實(shí)時(shí)渲染算法流程

為克服上述不足，本算法以傳統(tǒng)MC算法為基礎(chǔ)，針對雷達(dá)基數(shù)據(jù)的空間結(jié)構(gòu)特點(diǎn)，著重在三角面片頂點(diǎn)位置的計(jì)算、三角剖分構(gòu)型唯一性的判定以及包含等值面體素的加速遍歷，這3方面提出修改和優(yōu)化。雷達(dá)三維回波實(shí)時(shí)渲染的總體算法流程如圖1所示。

圖1 總體算法流程圖

2 雷達(dá)三維回波實(shí)時(shí)渲染算法

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

國內(nèi)常用的天氣雷達(dá)分為S波段、C波段和X波段3類，每個(gè)波段各有多種型號(hào)，每個(gè)型號(hào)的雷達(dá)又可采用多種方式的體掃模式。以我國東南沿海地區(qū)最常用的SA型為例，它的體掃模式主要有VCP11（Volume Coverage Pattern）、VCP21和VCP31這3種。其中，VCP21常被用于探測非強(qiáng)對流的明顯降水情況，該模式的體掃周期為6 min，方位角方向的分辨率為1°，徑向分辨率為0.25 km，探測距離為460 km［15］。1個(gè)體掃周期內(nèi)完成9個(gè)不同仰角掃描，仰角與高度、水平距離的關(guān)系如圖2所示。

圖2 SA雷達(dá)VCP21體掃模式

圖2中：hd為水平距離；vh為垂直高度，單位均為km；ev為探測仰角，單位為（°）。通常情況下，MC算法分析的數(shù)據(jù)在空間上等間隔分布結(jié)構(gòu)，即各個(gè)位置的體素密度相同。顯然，雷達(dá)探測數(shù)據(jù)結(jié)構(gòu)并不滿足這一特征，通行的做法是將極坐標(biāo)系下的數(shù)據(jù)結(jié)構(gòu)通過空間插值轉(zhuǎn)換成平面直角坐標(biāo)系下的三維結(jié)構(gòu)，但額外增加的體素信息不能準(zhǔn)確地反映大氣中的水氣狀況，除了方便計(jì)算，并不能帶來性能上的有益提升。本算法不做插值和坐標(biāo)系轉(zhuǎn)換，而直接根據(jù)雷達(dá)探測所構(gòu)成的數(shù)據(jù)空間結(jié)構(gòu)，分析各個(gè)等值面閾值與基本反射率之間的數(shù)值關(guān)系。

如圖3所示，a～h為雷達(dá)探測可識(shí)別分辨率下某處相鄰的8個(gè)頂點(diǎn)，各頂點(diǎn)的數(shù)值為該位置探測到的基本反射率，用Z（λ，θ，φ）表示，單位為dBZ。這8個(gè)頂點(diǎn)構(gòu)成一個(gè)六面體，將其視作一個(gè)體素，它是后續(xù)MC計(jì)算的基本單位。

圖3 極坐標(biāo)系下相鄰8頂點(diǎn)的空間分布

2.2 極坐標(biāo)系下三角面片頂點(diǎn)位置的計(jì)算

體素中的三角面片是擬合等值面的基本單位，要得到三角面片的空間位置，就需要準(zhǔn)確計(jì)算等值面在體素棱邊相交點(diǎn)的坐標(biāo)。通常在由X、Y、Z構(gòu)成的空間直角坐標(biāo)系下，可以利用線性插值計(jì)算各個(gè)交點(diǎn)的坐標(biāo)。雷達(dá)數(shù)據(jù)是極坐標(biāo)系下的空間三維結(jié)構(gòu)，計(jì)算方法需做相應(yīng)轉(zhuǎn)換。假設(shè)D1=（λ1，θ1，φ1）和D2=（λ2，θ2，φ2）是體素棱邊上兩個(gè)端點(diǎn)的坐標(biāo)，閾值為T1的等值面在該棱邊上有交點(diǎn)，那么該交點(diǎn)的坐標(biāo)P（λ，θ，φ）計(jì)算方法為：

2.3 基于相鄰體素狀態(tài)的快速判別算法

VCP21體掃模式類似圖3所示，由8個(gè)頂點(diǎn)構(gòu)成的六面體（體素）約有1.32×106個(gè)，采用MC算法對所有體素逐一遍歷，將消耗大量時(shí)間，且大多數(shù)體素與待求取的等值面并不相交。為提升計(jì)算效率，采用一種區(qū)域拓?fù)鋽U(kuò)展的思想，先任意選取一個(gè)與等值面相交的體素，將其作為基準(zhǔn)，篩選出與該體素相鄰且包含等值面的體素，再分別將這些體素作為基準(zhǔn)，采用相同方法擴(kuò)展篩選，直到找出所有與等值面相交的體素，再做后續(xù)處理。

上述過程的關(guān)鍵是如何判別與基準(zhǔn)體素相鄰的體素是否與等值面相交。如圖4所示，有4個(gè)相鄰的體素，每個(gè)體素包含8個(gè)頂點(diǎn)，其中CD是4個(gè)體素公有的邊，ABCD是上方2個(gè)體素公有的面。如果A、B、C、D4個(gè)頂點(diǎn)的基本反射率數(shù)值不全部小于或大于等值面的閾值，則可判定基準(zhǔn)體素的相鄰體素與等值面相交。圖4中MN表示左上角的基準(zhǔn)體素與等值面的交線，由于MN也存在于右上角的體素中，因此，右上角的體素也與等值面相交。

圖4 相鄰體素等值面與交線的示意圖

盡管A、B、C、D的取值不定，但相較于等值面的閾值，只有小于、大于或等于這3種狀態(tài)。為簡化表述，不妨假設(shè)當(dāng)頂點(diǎn)的基本反射率數(shù)值小于或等于等值面的閾值時(shí)，設(shè)定為狀態(tài)“T”，當(dāng)大于時(shí)設(shè)定為狀態(tài)“F”。那么，對于體素中任意一面的狀態(tài)可總結(jié)為如圖5所示的4種狀態(tài)。

圖5 體素中一個(gè)面與等值面相交的情況

圖5中，體素一個(gè)面的頂點(diǎn)用黑色實(shí)心圓和白色空心圓分別表示狀態(tài)“T”和“F”。其中，圖5（a）有1個(gè)頂點(diǎn)與其他頂點(diǎn)的狀態(tài)不同；圖5（b）、（c）各有2個(gè)頂點(diǎn)與其他頂點(diǎn)的狀態(tài)不同；圖5（d）中4個(gè)頂點(diǎn)為同一狀態(tài)。圖5中的虛線為體素的一個(gè)面與等值面相交的示意，相交的準(zhǔn)確位置可采用2.2節(jié)所述方法計(jì)算得到。其他情況均可通過旋轉(zhuǎn)或翻轉(zhuǎn)圖形歸并到圖5所示的4種情況之一。在遍歷計(jì)算相鄰體素時(shí)，可借助哈希表輔助算法執(zhí)行［16］。

2.4 體素中三角剖分構(gòu)型的唯一性判定算法

傳統(tǒng)MC算法中，當(dāng)狀態(tài)“T”和“F”的頂點(diǎn)如圖5（c）所示分別位于對角線的兩端，那么等值面與體素相交的方式會(huì)有兩種，因此存在歧義，如圖6所示。為保證三角面片拓?fù)浣Y(jié)構(gòu)的一致性，一些學(xué)者提出了子構(gòu)型拆分和雙曲線漸近線判別［17］等方法。無論是子構(gòu)型拆分中的遞歸迭代，還是雙曲線漸近線交點(diǎn)的三線性插值計(jì)算，都使得計(jì)算的時(shí)間復(fù)雜度增加了一個(gè)量級(jí)。

圖6 等值面與體素相交方式的二義性

為加快體素等值面的判定，并消除二義性，本算法借鑒雙曲線漸近線判別法的思路［18］，并提出一點(diǎn)改進(jìn)。通常情況下，等值面與體素相交形成的交線是一對雙曲線，該雙曲線與體素面上的位置關(guān)系如圖7所示。

圖7 雙曲線與體素面上的位置關(guān)系

可見，圖7（a）中雙曲線的兩支同時(shí)與體素的面相交，雙曲線將體素的面分割成3個(gè)區(qū)域。分別將兩組對邊上的交點(diǎn)用虛線連接，這兩條虛線形成的交點(diǎn)一定和其中一條對角線上的兩個(gè)頂點(diǎn)同處一個(gè)區(qū)域，即同為狀態(tài)“T”或“F”。通過比較該交點(diǎn)處的基本反射率值與閾值的大小關(guān)系，即可判定該二義性面與等值面的連接方式。

用圖8舉例說明，P1～P4是等值面與體素中一個(gè)面相交形成的4個(gè)交點(diǎn)，這4個(gè)交點(diǎn)的坐標(biāo)可用2.2節(jié)所述方法計(jì)算得到，再通過聯(lián)立P1、P2和P3、P4的直線方程，可計(jì)算出兩條直線相交處O的坐標(biāo)。由于A～D4個(gè)頂點(diǎn)的基本反射率值是已知的，因此可通過雙線性插值，計(jì)算出O點(diǎn)的基本反射率值，計(jì)算方法為：

圖8 二義性面的唯一性判定

式中：λA、θA和φA分別表示頂點(diǎn)A極坐標(biāo)下的3個(gè)變量；Z（A）頂點(diǎn)A處的基本反射率值；Z（O）即為兩條直接相交處O點(diǎn)的基本反射率值。同樣用狀態(tài)“T”和“F”來表示該值與等值面閾值的大小關(guān)系。當(dāng)Z（O）的狀態(tài)為“F”時(shí)，可確定圖8（a）的構(gòu)型，當(dāng)當(dāng)Z（O）的狀態(tài)為“T”時(shí)，可確定圖8（b）的構(gòu)型。

3 實(shí)驗(yàn)結(jié)果及分析

3.1 實(shí)驗(yàn)結(jié)果

為實(shí)現(xiàn)雷達(dá)三維回波的可視化，需將上述過程計(jì)算得到的一系列等值面用特定格式的文件輸出。目前主流的格式有：COLLADA是一種基于XML的3D交互文件格式，它通過對場景、節(jié)點(diǎn)、幾何形狀、著色器、力學(xué)和運(yùn)動(dòng)過程的描述，實(shí)現(xiàn)對三維矢量圖形的渲染和交互。圖形語言傳輸格式（Graphics Language Transmission Format，GLTF）是另一種常用的圖形語言交換格式，它用比XML更簡潔的JSON格式存儲(chǔ)信息，對OpenGL更友好，支持大數(shù)據(jù)量的3D Web實(shí)時(shí)渲染。

本實(shí)驗(yàn)選用GLTF格式，將一組閾值的等值面依次輸出，通過Blender軟件加載20 dBZ等值面渲染后的效果如圖9所示。

圖9 采用Blender渲染的雷達(dá)三維回波（20 dBZ）

在天氣分析等教學(xué)實(shí)踐過程中，往往需要結(jié)合地理信息對雷達(dá)回波的特征加以分析，目前諸如Three.js、Cesium和PlayCanvas等開源軟件都對GLTF格式數(shù)據(jù)有很好的支持［19］。其中，Cesium是一個(gè)基于JavaScript可創(chuàng)建三維場景和虛擬地球的地圖可視化框架，在基于Web的地圖動(dòng)態(tài)數(shù)據(jù)可視化方面具有優(yōu)異的性能。

本實(shí)驗(yàn)將上述閾值的等值面依次加載到Cesium，并用氣象上常用的色標(biāo)對不同閾值的等值面渲染相應(yīng)顏色，疊加地球影像和地圖圖像的顯示的效果如圖10所示。圖10（a）、（b）加載的是一同GLTF文件，不同之處在于觀測的視角和加載的地圖信息。

圖10 采用Cesium渲染的雷達(dá)三維回波

3.2 結(jié)果分析

改進(jìn)傳統(tǒng)MC算法的計(jì)算效能，提升雷達(dá)三維回波渲染和交互的實(shí)時(shí)性是本算法改進(jìn)的主要目標(biāo)。為得到客觀量化的評價(jià)指標(biāo)，以下測試采用相同的軟硬件環(huán)境。主要硬件配置為：Intel Xeon E5-1630 v4CPU，默認(rèn)主頻3.7 GHz，DDR3 2400 MHz 128 GB內(nèi)存。算法性能的評價(jià)主要包括兩個(gè)方面：一是計(jì)算耗時(shí)，它是評估等值面提取效率的重要指標(biāo)，直接關(guān)系到是否滿足雷達(dá)數(shù)據(jù)實(shí)時(shí)渲染的性能要求；二是所生成的三角面片的總數(shù)量，雖然更多的三角面片使得圖形更加平滑，能有效減少鋸齒，但對內(nèi)存的開銷很大。

表1中，“改進(jìn)1”是以傳統(tǒng)MC算法為基礎(chǔ)，采用2.1、2.2節(jié)所述方法，不做空間插值和坐標(biāo)轉(zhuǎn)換；“改進(jìn)2”是在“改進(jìn)1”的基礎(chǔ)上，增加2.3節(jié)所述方法進(jìn)行相鄰體素狀態(tài)的快速判別。針對傳統(tǒng)MC的等值面二義性問題，上述測試均采用了2.4節(jié)所述方法進(jìn)行處理。測試過程選用了不同時(shí)刻的雷達(dá)基數(shù)據(jù)，分別記錄不同算法下的計(jì)算耗時(shí)和三角面片數(shù)量，表1中的數(shù)值是10次測試的算法平均值取整。可見，傳統(tǒng)MC完成一次渲染需要27 s，其中大量時(shí)間消耗在坐標(biāo)轉(zhuǎn)換過程的插值計(jì)算。以及等值面提取時(shí)對大量無價(jià)值體素的分析。“改進(jìn)1”優(yōu)化了上述缺陷，既降低了計(jì)算耗時(shí)，也減少了三角面片的數(shù)量，通過Edge瀏覽器加載Cesium和GLTF數(shù)據(jù)，監(jiān)測內(nèi)存發(fā)現(xiàn)，Edge進(jìn)程占用內(nèi)存減少了56%。“改進(jìn)2”進(jìn)一步降低了體素遍歷的數(shù)量，將計(jì)算耗時(shí)再度提升34%。

表1 改進(jìn)算法性能對比

在實(shí)際操作過程中，采用“改進(jìn)2”的計(jì)算方法，每次加載新的雷達(dá)基數(shù)據(jù)仍需近7 s的等待響應(yīng)時(shí)間，如果將這些雷達(dá)基數(shù)據(jù)提前預(yù)處理好，生成GLTF文件存儲(chǔ)于服務(wù)器中，那么網(wǎng)絡(luò)加載GLTF文件并通過瀏覽器渲染的耗時(shí)將大幅降低到500 ms以內(nèi)，達(dá)到實(shí)時(shí)交互、快速響應(yīng)的目標(biāo)。通過網(wǎng)頁異步傳輸、預(yù)加載等技術(shù)，可實(shí)現(xiàn)時(shí)序三維回波動(dòng)畫的平滑播放，極大增強(qiáng)師生在天氣分析教學(xué)活動(dòng)中的體驗(yàn)。

4 結(jié) 語

利用曲面重建技術(shù)實(shí)現(xiàn)天氣雷達(dá)三維回波的實(shí)時(shí)渲染，為雷達(dá)氣象學(xué)、基于雷達(dá)資料的短臨天氣分析等教學(xué)實(shí)踐活動(dòng)提供了更加直觀、立體、可交互的應(yīng)用軟件平臺(tái)。本文設(shè)計(jì)的算法，針對雷達(dá)基數(shù)據(jù)的空間結(jié)構(gòu)特點(diǎn)，著重在三角面片頂點(diǎn)位置的計(jì)算、三角剖分構(gòu)型唯一性的判定以及包含等值面體素的加速遍歷3方面提出修改和優(yōu)化，測試結(jié)果表明，改進(jìn)后的算法計(jì)算耗時(shí)縮短至25%，內(nèi)存占用降低了56%，滿足師生教學(xué)過程中，雷達(dá)數(shù)據(jù)渲染與交互的實(shí)時(shí)性要求。該算法也可應(yīng)用于數(shù)值天氣模式、氣象衛(wèi)星等資料的三維渲染，具有一定的推廣應(yīng)用價(jià)值。