基于虛擬地球的電離層電子密度時空多維可視化方法*

郭長順，范湘濤，譚 劍?，詹 勤，亢敘杰

(1 中國科學(xué)院遙感與數(shù)字地球研究所 中國科學(xué)院數(shù)字地球重點(diǎn)實(shí)驗(yàn)室， 北京 100094； 2 中國科學(xué)院大學(xué), 北京 100049； 3 海南省地球觀測重點(diǎn)實(shí)驗(yàn)室， 海南 三亞 572029)(2018年9月12日收稿； 2018年11月8日收修改稿)

電離層是地球大氣層的重要圈層之一，探索電離層電子密度的空間分布狀況對研究空間天氣狀況[1]、電離層震前異常[2]、赤道異常和威德爾海異常現(xiàn)象[3]等都具有重要意義。

電離層電子密度的可視化研究不僅可以讓用戶以直觀的方式分析數(shù)據(jù)，而且有助于發(fā)現(xiàn)數(shù)據(jù)中所隱藏的復(fù)雜規(guī)律[4]。傳統(tǒng)的電離層數(shù)據(jù)可視化主要有面繪制[5-7]、定制kml文件[8]等方法，但這些方法只表達(dá)出數(shù)據(jù)的部分信息特征，無法反映出原始數(shù)據(jù)場的全貌和細(xì)節(jié)。直接體繪制方法從原始三維數(shù)據(jù)出發(fā)，可以顯示出三維數(shù)據(jù)場中各組成分量的空間分布狀況，主要包括足跡表法、錯切變形法、光線投射法和三維紋理映射法等主要算法[9]。這些體繪制方法已被廣泛應(yīng)用于空間電磁場、氣象和地球科學(xué)領(lǐng)域，以促進(jìn)可視化分析和科學(xué)研究，如三維磁層渲染[10]、颶風(fēng)[11]和臺風(fēng)[12]模擬、大氣溫度可視化[13]、地震數(shù)據(jù)可視化[14-15]，但很少用于電離層電子密度可視化。雖然基于球體八叉樹的空間網(wǎng)格被用于電離層電子密度的三維體可視化，但在索引時間序列空間數(shù)據(jù)時存儲各個八叉樹需要大量空間，這給后端數(shù)據(jù)存儲和網(wǎng)絡(luò)數(shù)據(jù)傳輸造成了挑戰(zhàn)；并且細(xì)節(jié)程度的提高會造成存儲復(fù)雜度成幾何級數(shù)的增加，即使硬件加速，可視化效率也很低[16-17]?；诠饩€投射的體繪制方法具有較好的可視化質(zhì)量，現(xiàn)在可以實(shí)現(xiàn)基于虛擬地球電磁環(huán)境體數(shù)據(jù)的可視化[18]，但尚未被用于電離層電子密度可視化，且其可視化程序無法在Web環(huán)境中實(shí)現(xiàn)。

為解決以上電離層電子密度可視化存在的問題，本文主要完成以下工作:

首先，為解決數(shù)據(jù)在Web上傳輸和實(shí)時可視化的問題，利用視頻壓縮編碼方法，將電離層電子密度數(shù)據(jù)壓縮為視頻流，便于數(shù)據(jù)在網(wǎng)絡(luò)中傳輸和實(shí)現(xiàn)不同時序數(shù)據(jù)的可視化；其次，為解決可視化質(zhì)量差的問題，利用基于GPU加速的光線投射體繪制方法對大區(qū)域、多維、多時相的電離層電子密度進(jìn)行可視化表達(dá)，并采用自適應(yīng)步長采樣以及早期光線終止法，提高繪制效率；最后，將電離層電子密度數(shù)據(jù)及豐富的地理背景集成在開源虛擬地球平臺Cesium上，可以讓更多專業(yè)人員和公眾看到更高質(zhì)量的電離層電子密度變化，并通過交互式查詢和探索原始數(shù)據(jù)，分析隱藏的電離層變化規(guī)律和關(guān)系。

1 數(shù)據(jù)和預(yù)處理

數(shù)據(jù)編碼壓縮是提高海量數(shù)據(jù)預(yù)處理能力的一個有效手段，主要有兩方面的優(yōu)點(diǎn)：一是降低內(nèi)存需求，同時提高時域相關(guān)數(shù)據(jù)的交互性；二是可以減少這些數(shù)據(jù)的傳輸和處理時間，這對基于遠(yuǎn)程實(shí)現(xiàn)的可視化系統(tǒng)來說更為重要[19]。因此，有必要對原始電離層電子密度數(shù)據(jù)進(jìn)行預(yù)處理，將數(shù)據(jù)壓縮為視頻流進(jìn)行傳輸。下面詳細(xì)介紹整個預(yù)處理過程。

1.1 數(shù)據(jù)來源及特點(diǎn)

本文用到的電離層電子密度數(shù)據(jù)是由中國科學(xué)院國家空間科學(xué)中心提供的電離層電子密度同化數(shù)據(jù)。數(shù)據(jù)格式為NetCDF，數(shù)據(jù)經(jīng)度范圍為70°～140°，緯度范圍為15°～55°，高度范圍為100～920 km，時間跨度為4 h，經(jīng)度和緯度間隔都為2°，高度間隔為20 km，時間間隔為1 h。

1.2 數(shù)據(jù)分解

原始電離層電子密度數(shù)據(jù)具有經(jīng)度、緯度、高度、時間四維特征，本文采用文獻(xiàn)[20]中的算法將電離層電子密度數(shù)據(jù)先按照時間維度劃分為一系列空間結(jié)構(gòu)相一致的三維體數(shù)據(jù)，然后將每一時相的數(shù)據(jù)按高度分層，處理為一系列結(jié)構(gòu)相一致的二維網(wǎng)格數(shù)據(jù)，如圖1所示，每個網(wǎng)格點(diǎn)包含某個空間位置的電子密度信息。

圖1 將四維電離層電子密度數(shù)據(jù)處理為一系列二維網(wǎng)格Fig.1 Split of the four-dimensional ionospheric electron density data into a series of two-dimensional grids

1.3 色彩空間轉(zhuǎn)換

當(dāng)原始電離層電子密度數(shù)據(jù)被處理為二維網(wǎng)格數(shù)據(jù)后，將數(shù)據(jù)從標(biāo)量空間轉(zhuǎn)換到RGB空間，以便后期將數(shù)據(jù)編碼為視頻。本文針對每一高度層中的網(wǎng)格，取每個網(wǎng)格4個頂點(diǎn)數(shù)據(jù)的平均值作為該網(wǎng)格的電子密度值。然后，將每個網(wǎng)格的電子密度值通過下式歸一化到[0,1]之間，以便于轉(zhuǎn)換到RGB色彩空間：

Di,j=(GEi,j-GEmin)/(GEmax-GEmin),

(1)

式中:GEi,j是每個網(wǎng)格所代表的電子密度值，GEmax和GEmin是整個電離層電子密度數(shù)據(jù)集的最大值和最小值，Di,j稱為原始標(biāo)量值。

將原始標(biāo)量值Di,j乘以255后賦予RGB色彩空間的R色彩分量，G色彩分量和B色彩分量都賦予0，因此轉(zhuǎn)換后的RGB色彩分量可以用下式表示：

(R,G,B)=(Di,j×255,0,0).

(2)

借助MATLAB軟件，設(shè)計程序?qū)⒉煌叨葘拥腞、G、B色彩分量合成為RGB圖像。圖2(a)～2(d)分別是色彩分量合成后不同高度層的RGB圖像。最后，將不同高度層的RGB圖像合成為一張代表某一時相電離層電子密度數(shù)據(jù)的大圖像。

圖2 不同高度層的RGB圖像Fig.2 RGB images of different height layers

1.4 基于視頻編碼的數(shù)據(jù)壓縮

當(dāng)不同時相的電離層電子密度原始數(shù)據(jù)被處理為RGB圖像后，利用視頻編碼方法將一系列RGB圖像壓縮為視頻數(shù)據(jù)。視頻文件的編碼和解碼時間受到多種因素的影響，其中，最大的影響因素是視頻壓縮編解碼器和色度子采樣模型。目前適用于Web環(huán)境的視頻編解碼器主要有H.264、VP8、VP9等，色度子采樣模型主要有YUV4∶4∶4、YUV 4∶2∶0、YUV 4∶2∶2等[20]。綜合考慮視頻質(zhì)量和解碼效率兩方面，本文選擇的視頻編解碼器為VP9，色度子采樣模型為YUV4∶4∶4。

2 基于GPU加速的電離層電子密度光線投射體繪制方法

為實(shí)現(xiàn)高質(zhì)量和高效率的可視化，本文基于WebGL可編程渲染管道，設(shè)計頂點(diǎn)著色程序和片段著色程序。

頂點(diǎn)著色程序針對電離層電子密度三維網(wǎng)格的每個頂點(diǎn)運(yùn)行一次，主要負(fù)責(zé)完成一系列坐標(biāo)轉(zhuǎn)換工作，最終確定頂點(diǎn)在屏幕上的位置。片段著色程序針對每個屏幕像素運(yùn)行一次，它主要負(fù)責(zé)計算頂點(diǎn)的紋理坐標(biāo)，進(jìn)行紋理重采樣，解析相應(yīng)的數(shù)據(jù)值，以及將數(shù)據(jù)值通過傳遞函數(shù)映射為相應(yīng)的顏色和不透明度，然后經(jīng)過一定的融合處理，把最終渲染結(jié)果顯示在屏幕上。整體流程如圖3所示。

2.1 電離層電子密度體數(shù)據(jù)坐標(biāo)轉(zhuǎn)換

在CPU中計算每個體素頂點(diǎn)在地理坐標(biāo)系下的三維坐標(biāo)，以三維數(shù)組的形式存儲。在實(shí)時渲染中，因?yàn)镚PU渲染管線不直接采用球面坐標(biāo)，因此3D幾何圖元必須以笛卡爾坐標(biāo)表示[11]。當(dāng)?shù)厍虮灰暈榍蝮w時，在頂點(diǎn)著色程序中，每個體素頂點(diǎn)的球面坐標(biāo)被轉(zhuǎn)換為原點(diǎn)位于地球中心的笛卡爾坐標(biāo)，轉(zhuǎn)換公式為：

圖3 基于GPU加速的電離層電子密度光線投射體繪制Fig.3 GPU-accelerated ray-casting-based volume rendering algorithm of ionospheric electron density

(3)

式中：α、β、r為球面坐標(biāo)，x、y、z為笛卡爾坐標(biāo)。

當(dāng)每個體素頂點(diǎn)的球面坐標(biāo)轉(zhuǎn)換為笛卡爾坐標(biāo)后，使用由Cesium提供的模型視圖投影(MVP)矩陣實(shí)現(xiàn)從笛卡爾坐標(biāo)到屏幕坐標(biāo)的最終坐標(biāo)轉(zhuǎn)換。

在光線投射體繪制方法中，不能通過笛卡爾坐標(biāo)系下的光線采樣位置來獲取體數(shù)據(jù)值，需要完成笛卡爾坐標(biāo)到紋理坐標(biāo)的轉(zhuǎn)換，這項(xiàng)工作是由片段著色程序來完成[18]。假設(shè)數(shù)據(jù)的計算范圍是A0(lat0,lon0,h0)至A1(lat1,lon1,h1)，首先把頂點(diǎn)著色程序傳遞過來的笛卡爾坐標(biāo)(x，y，z)轉(zhuǎn)換為球面坐標(biāo)(r,α,β)，轉(zhuǎn)換公式為：

(4)

將球面坐標(biāo)通過插值轉(zhuǎn)換為紋理坐標(biāo)的公式為：

(5)

只有當(dāng)紋理坐標(biāo)u，v和w的值在0～1之間時，該紋理坐標(biāo)才會對體繪制效果產(chǎn)生作用。如果紋理坐標(biāo)值超出該范圍，則不對該處的體數(shù)據(jù)進(jìn)行繪制，這樣可以減少無效數(shù)據(jù)的繪制，提高數(shù)據(jù)繪制效率。

2.2 自適應(yīng)步長采樣

本文在客戶端界面中設(shè)計了電離層電子密度層數(shù)控制模塊，用戶可以根據(jù)自身需求自由調(diào)節(jié)要顯示的層數(shù)。在重采樣過程中，為了既保證渲染質(zhì)量又提高渲染效率，根據(jù)最大和最小層之間的高度間隔自適應(yīng)計算采樣步長。當(dāng)最大和最小層之間的高度間隔由小變大時，增加采樣步長；反之，減小采樣步長。自適應(yīng)采樣步長計算公式為

Step=FixedStep×(Hmax-Hmin)/H0，

(6)

式中：FixedStep為兩層數(shù)據(jù)之間的采樣步長基數(shù)，Hmax和Hmin分別為要顯示的體數(shù)據(jù)的最大高度和最小高度，H0為兩層數(shù)據(jù)之間的高度間隔。

2.3 傳遞函數(shù)設(shè)計

傳遞函數(shù)實(shí)質(zhì)上是將標(biāo)量值映射為相應(yīng)的顏色值。為了在采樣點(diǎn)處應(yīng)用傳遞函數(shù)，應(yīng)首先根據(jù)采樣點(diǎn)的紋理坐標(biāo)得到采樣點(diǎn)處的標(biāo)量值。在片段著色程序中，當(dāng)采樣點(diǎn)在第z層和z+1層紋理之間時，通過采樣點(diǎn)的紋理坐標(biāo)得到這兩層相同紋理坐標(biāo)處的顏色值c1和c2。然后將c1、c2通過值解碼轉(zhuǎn)換為公式(1)所表示的原始標(biāo)量值D1、D2，對D1、D2在z和z+1層之間進(jìn)行線性插值得到采樣點(diǎn)處的原始標(biāo)量值。

本文設(shè)計多個可以在客戶端界面中由用戶自由調(diào)節(jié)的一維傳遞函數(shù)，每個傳遞函數(shù)都有不同個數(shù)的函數(shù)分界點(diǎn)和顏色區(qū)間，將傳遞函數(shù)轉(zhuǎn)化為一維紋理傳入GPU的片段著色程序參與計算。在計算電離層電子密度數(shù)據(jù)值對應(yīng)的顏色值時，首先計算出代表電離層電子密度的原始標(biāo)量值所對應(yīng)的紋理坐標(biāo)，確定其在哪兩個分界點(diǎn)之間，然后在兩個分界點(diǎn)之間進(jìn)行線性插值求得采樣點(diǎn)的顏色值。

2.4 顏色和不透明度合成

光線投射體繪制的圖像合成算法主要有兩種，分別是從前向后和從后向前進(jìn)行圖像合成[21]。本文采用第一種圖像合成算法，合成公式為

(7)

式中：C、A代表第m個采樣點(diǎn)對應(yīng)的顏色值和不透明度，Cm-1、Am-1代表前m-1個采樣點(diǎn)累積的顏色值和不透明度，Cm、Am代表光線穿過第m個采樣點(diǎn)后累積的顏色值和不透明度。

利用這種圖像合成算法的優(yōu)勢是可以通過設(shè)定適當(dāng)?shù)牟煌该鞫乳撝担绠?dāng)Am>0.98時提前終止采樣計算，以提高繪制效率。

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

為了對本文所述方法進(jìn)行驗(yàn)證，基于開源虛擬地球平臺Cesium，對多維多時相電離層電子密度數(shù)據(jù)的可視化效果和效率進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)利用HTML、CSS、JavaScript、GLSL等編程語言，在Sublime環(huán)境下進(jìn)行開發(fā)。利用JavaScript實(shí)現(xiàn)基于光線投射的體繪制框架，借助GLSL來實(shí)現(xiàn)頂點(diǎn)和片段著色程序的核心算法。用于視頻編碼的程序?yàn)镕Fmpeg，其包含用于VP9視頻編碼的庫libvpx-vp9?？梢暬布h(huán)境為CPU Intel Core i7-6700HQ 2.60 GHz，內(nèi)存為16 GB，顯卡為Nvidia GeForce GTX 960 M(4 GB)。系統(tǒng)運(yùn)行所需的服務(wù)器為node.js，測試所用客戶端為360、Chrome、Firefox等瀏覽器。

3.1 集成地理背景后的渲染質(zhì)量分析

與傳統(tǒng)電離層電子密度可視化方法[7,16](圖4(a)、4(b))相比，本文所述方法(圖5(a)～5(d))能形象直觀地展示出電離層電子密度的空間分布狀況，用戶可以在客戶端界面中自由選擇需要顯示哪一時相的數(shù)據(jù)以及數(shù)據(jù)層數(shù)，并且可以從多角度進(jìn)行觀察。此外，虛擬地球可以提供豐富的地理背景，在虛擬地球平臺上添加天地圖中文注記、天地圖邊界、以及地形等，借助地理背景，用戶可以更好的分析理解電離層電子密度數(shù)據(jù)。圖5(a)和5(b)所使用的傳遞函數(shù)具有3個分界點(diǎn)，圖5(c)和5(d)所使用的傳遞函數(shù)具有6個分界點(diǎn)。

圖4 傳統(tǒng)電離層電子密度可視化方法繪制結(jié)果Fig.4 Results of the traditional ionospheric electron density visualization rendering method

圖5 基于GPU加速的光線投射體繪制結(jié)果Fig.5 Results of GPU-accelerated ray-casting-based volume rendering algorithm

3.2 傳輸性能分析

為驗(yàn)證視頻流在數(shù)據(jù)傳輸效率上的優(yōu)勢，在局域網(wǎng)環(huán)境下，本文對基于視頻編碼方法得到的數(shù)據(jù)與網(wǎng)絡(luò)環(huán)境中使用的其他空間壓縮技術(shù)(基于有損壓縮的JPEG圖像編碼/基于無損壓縮的PNG圖像編碼)得到的數(shù)據(jù)進(jìn)行了數(shù)據(jù)傳輸速度的對比實(shí)驗(yàn)，實(shí)驗(yàn)所用瀏覽器為360、Chrome和Firefox，實(shí)驗(yàn)結(jié)果如表1所示。

從表1可以看出，在360、Chrome、Firefox瀏覽器中，基于視頻編碼的數(shù)據(jù)傳輸速度分別是JPEG和PNG編碼方法的50和65倍、36和48倍、12和20倍。因此，將原始數(shù)據(jù)壓縮為視頻流進(jìn)行傳輸可以明顯提高數(shù)據(jù)傳輸效率。

表1 數(shù)據(jù)傳輸速度比較Table 1 Comparison of data transfer speed fps

3.3 渲染性能分析

為驗(yàn)證提出的基于GPU加速的電離層電子密度光線投射體繪制方法的效率是否有所提高，本文進(jìn)行了渲染性能驗(yàn)證測試，其方法主要是評估在各種約束條件下計算密集型可視化系統(tǒng)的每秒幀數(shù)(fps)[22]。

在圖像合成過程中，在多個采樣步長基數(shù)下，對分別應(yīng)用自適應(yīng)步長采樣與固定步長采樣算法后的體繪制速度進(jìn)行對比，實(shí)驗(yàn)結(jié)果如表2所示。

表2 體繪制速度比較Table 2 Comparison of volume rendering speed

從表2可以看出，在采樣步長基數(shù)分別為100、300、500和1 000 m時，應(yīng)用自適應(yīng)步長采樣算法后體繪制速度分別提高2.5%、3.3%、12.4%和39.9%。因此，自適應(yīng)步長算法可以提高體繪制效率，而且在規(guī)定采樣步長范圍內(nèi)，采樣步長基數(shù)越大，體繪制效率提高越明顯。

在圖像合成過程中，當(dāng)采用自適應(yīng)步長采樣時，對是否應(yīng)用早期光線終止法的體繪制速度進(jìn)行了對比，實(shí)驗(yàn)結(jié)果如表3所示。

表3 體繪制速度比較Table 3 Comparison of volume rendering speed

從表3可以看出，在采樣步長基數(shù)分別為100、300、500、1 000 m時，應(yīng)用早期光線終止法后體繪制速度分別提高397.0%、124.1%、60.6%、10.4%。因此，在最終的圖像合成過程中，當(dāng)采用從前向后的圖像合成算法時，應(yīng)用早期光線終止法可以大幅度提高體繪制效率，而且在自適應(yīng)采樣步長基數(shù)較小時尤為明顯。

4 結(jié)論

本文從電離層電子密度的時空特征和實(shí)際需求出發(fā)，詳細(xì)論述數(shù)據(jù)組織、視頻壓縮編碼方法以及基于GPU加速的電離層電子密度光線投射體繪制方法。和已有的電離層電子密度可視化方法相比，本文所述方法解決了數(shù)據(jù)在Web上實(shí)時傳輸困難和可視化效果不佳的問題，實(shí)現(xiàn)了不同時相電離層電子密度的多層、動態(tài)、交互式可視化，在采樣步長基數(shù)為100 m時，采用自適應(yīng)步長采樣和早期光線終止法體繪制效率提升近4倍。另外，基于開源虛擬地球平臺Cesium，該方法可以輕松地集成到Web環(huán)境中，便于多用戶共享。

在未來的研究中，可以將本文所述方法應(yīng)用到其他的多維標(biāo)量場中，比如大氣溫度、電離層離子溫度等，以探究該方法的適用性。此外，將地學(xué)領(lǐng)域的空間分析方法和本文所述方法進(jìn)行結(jié)合，使其不僅具有數(shù)據(jù)可視化功能，而且有助于空間分析和知識發(fā)現(xiàn)。