基于溫度場的氣象時空數據動態可視化

徐曉文，聶蕓

（華北計算技術研究所北京100083）

隨著時代的發展，計算機的技術越發的先進，同時數據觀測和采集儀器的進步，可以生成大量的氣象數據。氣象數據具有很強的時空性，表現在一個氣象信息處理系統中的數據源既有同時間不同空間的數據，也有同空間不同時間的數據[1]。同時，通過不同尺度的觀測也會有不同的比例和精度。

氣象可視化科學研究已經比較成熟，能夠對氣象觀測的要素數據、分析預報的產品進行按空間高度層次、按時次的分析和顯示，如等值線分析顯示、填色分析與顯示等，在氣象水文預報業務和氣象水文信息應用系統中廣泛使用[2]。隨著氣象水文信息對日常生活、行業工作的影響日益被重視，在應用系統中提供氣象水文環境信息的感知能力成為迫切的需求，同時也為了更好地進行氣象環境應用分析，更加充分、直觀地展現氣象水文信息，對氣象信息可視化提出了更高的要求。以往單層次、單時次的靜止數據展現不能充分體現氣象數據時空連續的特征；天氣圖等面向專業人員的方式也難以被一般用戶所接受。

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

1 系統總體設計

1.1 系統實現整體步驟

在實現氣象時空數據動態可視化過程中，需要每一個中間圖像生成的數據，通過這些圖像的連續播放，從而實現氣象時空數據的動態可視化。在很多情況下，氣象數據都是單時次，不連續的，所以需要將這些不連續的單時次的數據組合起來，如果直接生成可視化圖像進行連續播放，那么大部分的氣象數據在感官上都存在明顯的撕裂感。所以需要利用這些原本的單時次，不連續的數據生成一連串的數據來補充這個畫面減少撕裂感的產生。這時候就需要應用到插值技術。

由于是做氣象時空數據動態可視化，所以需要進行插值做出中間動畫來實現動態展示。考慮到分段低次插值函數都有一致收斂性，但光滑性較差，對于氣象時空數據動態可視化要求有一定的光滑度，即最好保證存在二階連續導數。所以本文采用最常用的三次樣條函數[3]。三次樣條插值（Cubic Spline Interpolation）簡稱Spline插值，我們選擇的是三次樣條插值[3]。這種插值技術比較容易生成平滑的曲線，并且在多個關鍵點上數據相關，在保證關鍵的準確的獲取到的數據的精度上，可以合理的生成出相關的數據[4]。同時由于數據是矩陣的形式，所以利用近似最鄰近特征匹配方法進行采樣，獲取出一連串的點，從而構建出合理的三次樣條函數[5]。

之后通過線性條插值補全數據。通過這種方法形成的新的矩陣，不僅在數據上連貫性良好，還在初始氣象數據矩陣十分精確。

1.2 系統實現流程設計

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

假定初始存在5個數據場，數據矩陣分別為M1，…，M5。這是一個固定時間間隔的序列矩陣，分別依照時間順序編號1，2，…，5。如圖2所示。

圖1 系統流程圖

圖2 初始矩陣序列圖

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

2 系統實現步驟

2.1 系統初始數據的獲取

文中研究中假定初始數據為矩陣型網格數據。這種類型的數據是實際業務中最常見的格式。實際使用中，通過氣象水文資料的數據處理，就能夠直接得到這類數據，通常按照行列順序存儲。對于氣象實況等基于觀測站點的離散數據，可以通過網格化處理形成矩陣型數據。

氣象數據的網格節點數（行列數）一般是按照在一定地理范圍內固定間隔分布來計算的，比如全球預報數據，按照4度、1度劃分，可形成90*45，360*180的數據矩陣[6]。

本章以4度劃分的矩陣為例進行算法說明。

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

2.2 提取矩陣的關鍵點

文中存在5個初始數據矩陣M1，…，M5，這是一個依照時間順序排序的序列矩陣，分別賦予編號1，2，3，4，5。現在本文需要在這些矩陣中建立一定的映射關系用來方便插值。取矩陣M1和M2作為樣本，都是同等大小的m×n矩陣，要求矩陣M1的M1（x，y）{（0＜x＜90，0＜y＜45，且x，y為整數，其中x，y對應的是矩陣中的坐標位置）}求對應到矩陣M2的最近臨近相似點[7]。

考慮到數據矩陣空間上的對應與時間上的連續，所以需要考慮空間對應關系以及時間上的變化關系[8]。

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

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

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

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

矩陣中任何的元素。之后選取I最小的點作為匹配。實現過程如下所示[9]：

1）首先將原始數據矩陣M1與M2進行分割，通過多區域采樣來提高數據的精確度。如圖3所示，兩個矩陣按照紅線所示分成等分的9個部分，分別是取橫軸和縱軸的三分之一作為一個分劃位置。這些分給出來的矩陣依照先從左到右再從上到下的規則 分 別 命 名 為M1p1，M1p2，…，M1p9與M2p1，M2p2，…，M2p9。之后各個相應的區域找對應關系，比如M1p1與M2p1對應；

圖3 原始數據矩陣區域劃分示意圖

2）以M1p1與M2p1對應為例，本文采用的對應關系如圖4所示。為圓點到星點的匹配過程。分割后M1p1與M2p1均為m/3×n/3大小的數據矩陣，其中圓點的取值為這個區域的中心點M1p1（m/6，n/6），星點的位置對應矩陣如下：

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

圖4 數據關鍵點的查找示意圖

3）后運算影響數值數據I，當存在多個相同的I時，并且剛好I的取值為最小，那么根據氣象時空數據時間和空間都具有變換的關系，采用匹配點的順時針方向順序選取第一個I值最小的點作為匹配點。

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

2.3 獲取三次樣條函數

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

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

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

2）在此處x的取值為T，S（x）的取值為Nm。令：

則Mi滿足如下n-1個方程：

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

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

2.4 得到插值數據

根據2.3得到的三次樣條插值的方程，假設本文需要求取M1p1和M2p1之間的插值數據。假設兩個矩陣之間插值數據為20個。這些數據和M1p1和M2p1是等時間間隔的，所以假設序號也是等間隔的。那么假設這些數據為N1l2p1n1，N1l2p1n2，…，N1l2p1n20。那根據時間次序的序列的劃分T=[1，1.05，1.10，…，2]。根據上一小節的方程，可以得到M1p1和M2p1之間的三次樣條函數，其中自變量為T，通過T的數值來確定對應的插值點的數值。

如圖5所示，在兩個矩陣之間插入數據，不過現在插入的只是一些關鍵點的數據，所以插入矩陣的數據不完全，它插入的是星號的數據，所以需要對剩余部分進行填充。

圖5 中間數據的插入過程

以N1l2p1n1為例子，假設已經取到了N1l2p1n1，現在利用這個數據進行補全成插值矩陣M1l2p1n1，這個插值矩陣對應的序號為1.05。

圖6 補全插入數據的過程

2.5 進行動態可視化

如上所示，M1l2n1已經獲取到了，同理可以得到更多的插值數據。本文采用的是現有的可視化手段，插值矩陣和原始數據構造相同，所以可以直接進行可視化展示[14]。

3 系統成果與分析

3.1 系統成果展現

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

為了方便觀察，圖8為圖7通過提取的左上角部分區域的可視化顯示圖片，可以觀測到等值線有變化過程，但是大致的區域的變動不大。說明這個數據插值進行可視化是有效果的。同時在運行過程中并沒有感受到明顯的卡頓，而且動畫的過度也比較流暢，證明了優化的效果。

3.2 系統分析

本文研究需要大量的成熟數值方法、圖形圖像處理方法支撐，為了能夠快速實驗多種方法，使用了科學計算工具matlab作為實驗平臺。

實驗中采取5個連續45×90的初始數據矩陣，每個數據之間插入300個數據，去除可視化步驟。在實驗中內存占用如下：

圖7 可視化顯示序列圖

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

物理內存：4135 MB。

交換頁面：5036 MB。

虛擬內存：6341 MB。

在使用e3處理器，16 g內存的環境下，所需要的運行時間為15.051630 s；即為每個數據插入平均耗時不到0.02 s，滿足可視化要求。

4 結論

本論文研究基于已有的氣象可視化技術基礎，在不改變可視化的具體實現的過程情況下，利用現有的可視化數據，通過多種方法實現中間插值，利用插值的數據進行可視化展現，然后連續播放這些可視化效果，從而獲得氣象時空上連續的可視化效果。論文設計的方法，分離了可視化方法和數據處理方法，使得動態展現的效果能夠針對多種可視化方法實施，便于擴展應用；數據處理方法簡單，可并行執行，便于應用實現。論文有些工作還有待進一步展開[15]，包括：

1）可視化效果可擴展。考慮到與現有系統的結合使用，本文在說明氣象時空數據可視化方法時，僅選取了等值線場和色彩場方法。對于時空數據可視化的呈現方式本身，也是值得研究的問題。

2）方法驗證和改進。論文采用的實驗數據類型有限，還需要采用不同類型不同規模的真實數據，對方法進行測試，更好地確定方法的適應性、有效性。