基于溫度場(chǎng)的氣象時(shí)空數(shù)據(jù)動(dòng)態(tài)可視化

徐曉文，聶蕓

（華北計(jì)算技術(shù)研究所北京100083）

隨著時(shí)代的發(fā)展，計(jì)算機(jī)的技術(shù)越發(fā)的先進(jìn)，同時(shí)數(shù)據(jù)觀測(cè)和采集儀器的進(jìn)步，可以生成大量的氣象數(shù)據(jù)。氣象數(shù)據(jù)具有很強(qiáng)的時(shí)空性，表現(xiàn)在一個(gè)氣象信息處理系統(tǒng)中的數(shù)據(jù)源既有同時(shí)間不同空間的數(shù)據(jù)，也有同空間不同時(shí)間的數(shù)據(jù)[1]。同時(shí)，通過(guò)不同尺度的觀測(cè)也會(huì)有不同的比例和精度。

氣象可視化科學(xué)研究已經(jīng)比較成熟，能夠?qū)庀笥^測(cè)的要素?cái)?shù)據(jù)、分析預(yù)報(bào)的產(chǎn)品進(jìn)行按空間高度層次、按時(shí)次的分析和顯示，如等值線分析顯示、填色分析與顯示等，在氣象水文預(yù)報(bào)業(yè)務(wù)和氣象水文信息應(yīng)用系統(tǒng)中廣泛使用[2]。隨著氣象水文信息對(duì)日常生活、行業(yè)工作的影響日益被重視，在應(yīng)用系統(tǒng)中提供氣象水文環(huán)境信息的感知能力成為迫切的需求，同時(shí)也為了更好地進(jìn)行氣象環(huán)境應(yīng)用分析，更加充分、直觀地展現(xiàn)氣象水文信息，對(duì)氣象信息可視化提出了更高的要求。以往單層次、單時(shí)次的靜止數(shù)據(jù)展現(xiàn)不能充分體現(xiàn)氣象數(shù)據(jù)時(shí)空連續(xù)的特征；天氣圖等面向?qū)I(yè)人員的方式也難以被一般用戶所接受。

本論文研究基于已有的氣象可視化技術(shù)基礎(chǔ)，在不改變可視化的具體實(shí)現(xiàn)的過(guò)程情況下，利用現(xiàn)有的可視化數(shù)據(jù)，通過(guò)多種方法實(shí)現(xiàn)中間插值，從而獲得氣象時(shí)空上連續(xù)的可視化效果。

1 系統(tǒng)總體設(shè)計(jì)

1.1 系統(tǒng)實(shí)現(xiàn)整體步驟

在實(shí)現(xiàn)氣象時(shí)空數(shù)據(jù)動(dòng)態(tài)可視化過(guò)程中，需要每一個(gè)中間圖像生成的數(shù)據(jù)，通過(guò)這些圖像的連續(xù)播放，從而實(shí)現(xiàn)氣象時(shí)空數(shù)據(jù)的動(dòng)態(tài)可視化。在很多情況下，氣象數(shù)據(jù)都是單時(shí)次，不連續(xù)的，所以需要將這些不連續(xù)的單時(shí)次的數(shù)據(jù)組合起來(lái)，如果直接生成可視化圖像進(jìn)行連續(xù)播放，那么大部分的氣象數(shù)據(jù)在感官上都存在明顯的撕裂感。所以需要利用這些原本的單時(shí)次，不連續(xù)的數(shù)據(jù)生成一連串的數(shù)據(jù)來(lái)補(bǔ)充這個(gè)畫(huà)面減少撕裂感的產(chǎn)生。這時(shí)候就需要應(yīng)用到插值技術(shù)。

由于是做氣象時(shí)空數(shù)據(jù)動(dòng)態(tài)可視化，所以需要進(jìn)行插值做出中間動(dòng)畫(huà)來(lái)實(shí)現(xiàn)動(dòng)態(tài)展示?？紤]到分段低次插值函數(shù)都有一致收斂性，但光滑性較差，對(duì)于氣象時(shí)空數(shù)據(jù)動(dòng)態(tài)可視化要求有一定的光滑度，即最好保證存在二階連續(xù)導(dǎo)數(shù)。所以本文采用最常用的三次樣條函數(shù)[3]。三次樣條插值（Cubic Spline Interpolation）簡(jiǎn)稱(chēng)Spline插值，我們選擇的是三次樣條插值[3]。這種插值技術(shù)比較容易生成平滑的曲線，并且在多個(gè)關(guān)鍵點(diǎn)上數(shù)據(jù)相關(guān)，在保證關(guān)鍵的準(zhǔn)確的獲取到的數(shù)據(jù)的精度上，可以合理的生成出相關(guān)的數(shù)據(jù)[4]。同時(shí)由于數(shù)據(jù)是矩陣的形式，所以利用近似最鄰近特征匹配方法進(jìn)行采樣，獲取出一連串的點(diǎn)，從而構(gòu)建出合理的三次樣條函數(shù)[5]。

之后通過(guò)線性條插值補(bǔ)全數(shù)據(jù)。通過(guò)這種方法形成的新的矩陣，不僅在數(shù)據(jù)上連貫性良好，還在初始?xì)庀髷?shù)據(jù)矩陣十分精確。

1.2 系統(tǒng)實(shí)現(xiàn)流程設(shè)計(jì)

圖1表示了實(shí)驗(yàn)的具體流程。根據(jù)流程圖示意，首先數(shù)據(jù)格式要求為矩陣形式。那么初始數(shù)據(jù)為給定的一連串的時(shí)間上的矩陣數(shù)據(jù)。

假定初始存在5個(gè)數(shù)據(jù)場(chǎng)，數(shù)據(jù)矩陣分別為M1，…，M5。這是一個(gè)固定時(shí)間間隔的序列矩陣，分別依照時(shí)間順序編號(hào)1，2，…，5。如圖2所示。

圖1 系統(tǒng)流程圖

圖2 初始矩陣序列圖

根據(jù)圖2所示，存在一個(gè)矩陣序列，包含依照時(shí)間排序的5個(gè)數(shù)據(jù)矩陣（其中2和4沒(méi)有畫(huà)出來(lái)），矩陣大小全部都為m×n，它們依照時(shí)間順序編號(hào)為1，2，3，4，5，本文主要實(shí)現(xiàn)的就是在這一個(gè)連續(xù)的數(shù)據(jù)矩陣序列中插入中間數(shù)據(jù)并利用這些數(shù)據(jù)和原始的數(shù)據(jù)生成一個(gè)動(dòng)態(tài)的可視化效果。根據(jù)生成動(dòng)畫(huà)的原理，為了獲得流暢的動(dòng)畫(huà)效果，需要在這些數(shù)據(jù)矩陣插入中間數(shù)據(jù)。

2 系統(tǒng)實(shí)現(xiàn)步驟

2.1 系統(tǒng)初始數(shù)據(jù)的獲取

文中研究中假定初始數(shù)據(jù)為矩陣型網(wǎng)格數(shù)據(jù)。這種類(lèi)型的數(shù)據(jù)是實(shí)際業(yè)務(wù)中最常見(jiàn)的格式。實(shí)際使用中，通過(guò)氣象水文資料的數(shù)據(jù)處理，就能夠直接得到這類(lèi)數(shù)據(jù)，通常按照行列順序存儲(chǔ)。對(duì)于氣象實(shí)況等基于觀測(cè)站點(diǎn)的離散數(shù)據(jù)，可以通過(guò)網(wǎng)格化處理形成矩陣型數(shù)據(jù)。

氣象數(shù)據(jù)的網(wǎng)格節(jié)點(diǎn)數(shù)（行列數(shù)）一般是按照在一定地理范圍內(nèi)固定間隔分布來(lái)計(jì)算的，比如全球預(yù)報(bào)數(shù)據(jù)，按照4度、1度劃分，可形成90*45，360*180的數(shù)據(jù)矩陣[6]。

本章以4度劃分的矩陣為例進(jìn)行算法說(shuō)明。

同時(shí)，為了獲取到比較合理的數(shù)據(jù)，本文在實(shí)驗(yàn)中所采用的數(shù)據(jù)為NOAA（美國(guó)國(guó)家海洋和大氣管理局）的數(shù)據(jù)（https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html），下 載 Air Temperature數(shù)據(jù)并提取文件里的數(shù)據(jù)，縮小矩陣規(guī)模進(jìn)行實(shí)驗(yàn)。

2.2 提取矩陣的關(guān)鍵點(diǎn)

文中存在5個(gè)初始數(shù)據(jù)矩陣M1，…，M5，這是一個(gè)依照時(shí)間順序排序的序列矩陣，分別賦予編號(hào)1，2，3，4，5?，F(xiàn)在本文需要在這些矩陣中建立一定的映射關(guān)系用來(lái)方便插值。取矩陣M1和M2作為樣本，都是同等大小的m×n矩陣，要求矩陣M1的M1（x，y）{（0＜x＜90，0＜y＜45，且x，y為整數(shù)，其中x，y對(duì)應(yīng)的是矩陣中的坐標(biāo)位置）}求對(duì)應(yīng)到矩陣M2的最近臨近相似點(diǎn)[7]。

考慮到數(shù)據(jù)矩陣空間上的對(duì)應(yīng)與時(shí)間上的連續(xù)，所以需要考慮空間對(duì)應(yīng)關(guān)系以及時(shí)間上的變化關(guān)系[8]。

本文所采用的方法就是求取

這個(gè)由M2矩陣的一部分元素組成的新的矩陣與M1(x,y)最相近的，此處最相近不僅僅代表空間上的相近，也有時(shí)間上數(shù)值的相近。由于采用的是近似最臨近特征匹配，所以本文需要考慮兩個(gè)矩陣中的距離因素，即為空間要素，假設(shè)這個(gè)距離為D12。為了確定距離的具體數(shù)值，本文做了以下的定義：

在M2中求取M1(x,y)中的近似最臨近特征，依照本文設(shè)定距離D12的數(shù)值為M1(x,y)與M2(x,y)對(duì)應(yīng)，那么距離D12的數(shù)值為0，然后與M2(x-1,y)，M2(x,y-1)，M2(x+1,y)與M2(x,y+1)對(duì)應(yīng)，那么距離D12的數(shù)值為1，其他距離D12的數(shù)值為2。

之后應(yīng)該計(jì)算所有條件下對(duì)每個(gè)對(duì)應(yīng)匹配占有的影響數(shù)值I，所利用的公式I=（M2（x，y）-M1（x，y））+D12。其中M1(x,y)可以分別變換成

矩陣中任何的元素。之后選取I最小的點(diǎn)作為匹配。實(shí)現(xiàn)過(guò)程如下所示[9]：

1）首先將原始數(shù)據(jù)矩陣M1與M2進(jìn)行分割，通過(guò)多區(qū)域采樣來(lái)提高數(shù)據(jù)的精確度。如圖3所示，兩個(gè)矩陣按照紅線所示分成等分的9個(gè)部分，分別是取橫軸和縱軸的三分之一作為一個(gè)分劃位置。這些分給出來(lái)的矩陣依照先從左到右再?gòu)纳系较碌囊?guī)則 分 別 命 名 為M1p1，M1p2，…，M1p9與M2p1，M2p2，…，M2p9。之后各個(gè)相應(yīng)的區(qū)域找對(duì)應(yīng)關(guān)系，比如M1p1與M2p1對(duì)應(yīng)；

圖3 原始數(shù)據(jù)矩陣區(qū)域劃分示意圖

2）以M1p1與M2p1對(duì)應(yīng)為例，本文采用的對(duì)應(yīng)關(guān)系如圖4所示。為圓點(diǎn)到星點(diǎn)的匹配過(guò)程。分割后M1p1與M2p1均為m/3×n/3大小的數(shù)據(jù)矩陣，其中圓點(diǎn)的取值為這個(gè)區(qū)域的中心點(diǎn)M1p1（m/6，n/6），星點(diǎn)的位置對(duì)應(yīng)矩陣如下：

之后的對(duì)應(yīng)關(guān)系重新依照這個(gè)方法對(duì)應(yīng)。同理M2與M3之間的關(guān)系重新匹配。每個(gè)關(guān)鍵點(diǎn)取由前一個(gè)中心點(diǎn)確定的作為樣本數(shù)據(jù)。即第一個(gè)樣本點(diǎn)的數(shù)值是Nm1p1，根據(jù)這個(gè)關(guān)系找到的是Nm2p1，但是Nm2p1可能已經(jīng)不是矩陣M2p1(m/6,n/6)的數(shù)據(jù)了。但是尋找Nm3p1的時(shí)候仍然以矩陣M2p1的M2p1(m/6,n/6)元素開(kāi)始參考運(yùn)算。Nm2p1并不會(huì)參與到M2p1與M3p1的運(yùn)算中；

圖4 數(shù)據(jù)關(guān)鍵點(diǎn)的查找示意圖

3）后運(yùn)算影響數(shù)值數(shù)據(jù)I，當(dāng)存在多個(gè)相同的I時(shí)，并且剛好I的取值為最小，那么根據(jù)氣象時(shí)空數(shù)據(jù)時(shí)間和空間都具有變換的關(guān)系，采用匹配點(diǎn)的順時(shí)針?lè)较蝽樞蜻x取第一個(gè)I值最小的點(diǎn)作為匹配點(diǎn)。

通過(guò)上述方法，可以在數(shù)據(jù)矩陣M1和M2之間建立出一個(gè)對(duì)應(yīng)的關(guān)系。本文開(kāi)始分割矩陣M1與M2，得到了多個(gè)數(shù)據(jù)矩陣，分別為M1p1，M1p2，…，M1p9與M2p1，M2p2，…，M2p9。之后以M1p1與M2p1為例，通過(guò)M1p1的中心點(diǎn)找到對(duì)應(yīng)的矩陣M2p1的所需求的數(shù)據(jù)和位置，之后再利用矩陣M2p1的中心點(diǎn)找到矩陣M3p1對(duì)應(yīng)的數(shù)據(jù)和位置，之后就可以得到一連串的5個(gè)數(shù)值，分別為Nm1p1，Nm2p1，… ，Nm5p1。N就是全部 5個(gè)序列矩陣演變關(guān)系的一個(gè)基準(zhǔn)向量，建立了5個(gè)矩陣的對(duì)應(yīng)關(guān)系。同理可以獲得其他的分塊矩陣的基準(zhǔn)向量。本文以之前的矩陣獲取到的基準(zhǔn)向量為例。

2.3 獲取三次樣條函數(shù)

根據(jù)2.2可知，通過(guò)近似最臨近特征匹配本文在初始數(shù)據(jù)上獲取到了一連串的5個(gè)數(shù)值，分別為Nm1p1，Nm2p1，… ，Nm5p1。那么接下來(lái)為了構(gòu)建出這些關(guān)鍵點(diǎn)的中間插值數(shù)據(jù)，本文考慮到氣象數(shù)據(jù)的流體性質(zhì)，采用三次樣條插值[10]。首先根據(jù)初始數(shù)據(jù)是一個(gè)存在時(shí)間順序關(guān)系的數(shù)據(jù)矩陣，那么根據(jù)這個(gè)信息，賦予一個(gè)變量T。那么假設(shè)T=[1，2，3，4，5]。同理存在隨著這個(gè)序數(shù)改變的變量，即為上一小節(jié)獲取的 5 個(gè)數(shù)值Nm=[Nm1p1，Nm2p1，…，Nm5p1]，那么就可以根據(jù)這個(gè)條件得到三次樣條函數(shù)。主要步驟如下所示[11]：

1）三次樣條插值多項(xiàng)式Sn(x)在每個(gè)小區(qū)間[xi-1,xi]上是3次多項(xiàng)式，其在此區(qū)間上的表達(dá)式如下：

其中，hi=xi+1-xi。因此，只要確定了Mi的值，就確定了整個(gè)表達(dá)式；

2）在此處x的取值為T(mén)，S（x）的取值為Nm。令：

則Mi滿足如下n-1個(gè)方程：

本文采用第一種邊界方法，所以得到以下的方程：

3）根據(jù)上文信息，存在n+1個(gè)方程，有著n+1個(gè)Mi未知數(shù)，那么根據(jù)得到的方程解出Mi，就可以得到方程的解。

2.4 得到插值數(shù)據(jù)

根據(jù)2.3得到的三次樣條插值的方程，假設(shè)本文需要求取M1p1和M2p1之間的插值數(shù)據(jù)。假設(shè)兩個(gè)矩陣之間插值數(shù)據(jù)為20個(gè)。這些數(shù)據(jù)和M1p1和M2p1是等時(shí)間間隔的，所以假設(shè)序號(hào)也是等間隔的。那么假設(shè)這些數(shù)據(jù)為N1l2p1n1，N1l2p1n2，…，N1l2p1n20。那根據(jù)時(shí)間次序的序列的劃分T=[1，1.05，1.10，…，2]。根據(jù)上一小節(jié)的方程，可以得到M1p1和M2p1之間的三次樣條函數(shù)，其中自變量為T(mén)，通過(guò)T的數(shù)值來(lái)確定對(duì)應(yīng)的插值點(diǎn)的數(shù)值。

如圖5所示，在兩個(gè)矩陣之間插入數(shù)據(jù)，不過(guò)現(xiàn)在插入的只是一些關(guān)鍵點(diǎn)的數(shù)據(jù)，所以插入矩陣的數(shù)據(jù)不完全，它插入的是星號(hào)的數(shù)據(jù)，所以需要對(duì)剩余部分進(jìn)行填充。

圖5 中間數(shù)據(jù)的插入過(guò)程

以N1l2p1n1為例子，假設(shè)已經(jīng)取到了N1l2p1n1，現(xiàn)在利用這個(gè)數(shù)據(jù)進(jìn)行補(bǔ)全成插值矩陣M1l2p1n1，這個(gè)插值矩陣對(duì)應(yīng)的序號(hào)為1.05。

圖6 補(bǔ)全插入數(shù)據(jù)的過(guò)程

2.5 進(jìn)行動(dòng)態(tài)可視化

如上所示，M1l2n1已經(jīng)獲取到了，同理可以得到更多的插值數(shù)據(jù)。本文采用的是現(xiàn)有的可視化手段，插值矩陣和原始數(shù)據(jù)構(gòu)造相同，所以可以直接進(jìn)行可視化展示[14]。

3 系統(tǒng)成果與分析

3.1 系統(tǒng)成果展現(xiàn)

本文采用最常用的等值線方法進(jìn)行可視化。可視化效果如圖7所示。

為了方便觀察，圖8為圖7通過(guò)提取的左上角部分區(qū)域的可視化顯示圖片，可以觀測(cè)到等值線有變化過(guò)程，但是大致的區(qū)域的變動(dòng)不大。說(shuō)明這個(gè)數(shù)據(jù)插值進(jìn)行可視化是有效果的。同時(shí)在運(yùn)行過(guò)程中并沒(méi)有感受到明顯的卡頓，而且動(dòng)畫(huà)的過(guò)度也比較流暢，證明了優(yōu)化的效果。

3.2 系統(tǒng)分析

本文研究需要大量的成熟數(shù)值方法、圖形圖像處理方法支撐，為了能夠快速實(shí)驗(yàn)多種方法，使用了科學(xué)計(jì)算工具matlab作為實(shí)驗(yàn)平臺(tái)。

實(shí)驗(yàn)中采取5個(gè)連續(xù)45×90的初始數(shù)據(jù)矩陣，每個(gè)數(shù)據(jù)之間插入300個(gè)數(shù)據(jù)，去除可視化步驟。在實(shí)驗(yàn)中內(nèi)存占用如下：

圖7 可視化顯示序列圖

圖8 可視化顯示部分序列圖

物理內(nèi)存：4135 MB。

交換頁(yè)面：5036 MB。

虛擬內(nèi)存：6341 MB。

在使用e3處理器，16 g內(nèi)存的環(huán)境下，所需要的運(yùn)行時(shí)間為15.051630 s；即為每個(gè)數(shù)據(jù)插入平均耗時(shí)不到0.02 s，滿足可視化要求。

4 結(jié)論

本論文研究基于已有的氣象可視化技術(shù)基礎(chǔ)，在不改變可視化的具體實(shí)現(xiàn)的過(guò)程情況下，利用現(xiàn)有的可視化數(shù)據(jù)，通過(guò)多種方法實(shí)現(xiàn)中間插值，利用插值的數(shù)據(jù)進(jìn)行可視化展現(xiàn)，然后連續(xù)播放這些可視化效果，從而獲得氣象時(shí)空上連續(xù)的可視化效果。論文設(shè)計(jì)的方法，分離了可視化方法和數(shù)據(jù)處理方法，使得動(dòng)態(tài)展現(xiàn)的效果能夠針對(duì)多種可視化方法實(shí)施，便于擴(kuò)展應(yīng)用；數(shù)據(jù)處理方法簡(jiǎn)單，可并行執(zhí)行，便于應(yīng)用實(shí)現(xiàn)。論文有些工作還有待進(jìn)一步展開(kāi)[15]，包括：

1）可視化效果可擴(kuò)展?？紤]到與現(xiàn)有系統(tǒng)的結(jié)合使用，本文在說(shuō)明氣象時(shí)空數(shù)據(jù)可視化方法時(shí)，僅選取了等值線場(chǎng)和色彩場(chǎng)方法。對(duì)于時(shí)空數(shù)據(jù)可視化的呈現(xiàn)方式本身，也是值得研究的問(wèn)題。

2）方法驗(yàn)證和改進(jìn)。論文采用的實(shí)驗(yàn)數(shù)據(jù)類(lèi)型有限，還需要采用不同類(lèi)型不同規(guī)模的真實(shí)數(shù)據(jù)，對(duì)方法進(jìn)行測(cè)試，更好地確定方法的適應(yīng)性、有效性。