基于傳輸函數(shù)的中尺度渦旋時空連續(xù)可視化

田豐林，朱新升，劉巍，韓妍嬌，陳戈*

( 1. 中國海洋大學(xué) 信息科學(xué)與工程學(xué)院，山東 青島 266100；2. 青島海洋科學(xué)與技術(shù)試點國家實驗室 區(qū)域海洋動力學(xué)與數(shù)值模擬功能實驗室，山東 青島 266237；3. 青島市計量技術(shù)研究院，山東 青島 266100；4. 哈爾濱市勘察測繪研究院，黑龍江 哈爾濱 150010)

1 引言

中尺度渦旋是旋轉(zhuǎn)運動著的并伴隨有溫度異常的水團，它的直徑一般為數(shù)十至數(shù)百千米，壽命可達數(shù)天到數(shù)月不等[1]。中尺度渦旋廣泛存在于全球海洋中，根據(jù)旋轉(zhuǎn)方向不同，可將它們分為反氣旋渦和氣旋渦，在北半球，反氣旋渦順時針旋轉(zhuǎn)，氣旋渦逆時針旋轉(zhuǎn)，而南半球則相反[2]。它們在動量、質(zhì)量、熱量、營養(yǎng)物質(zhì)以及鹽和其他海水元素的運輸中發(fā)揮著重要的作用，影響著海洋環(huán)流、水分布和海洋生物分布[3–5]。中尺度渦旋的發(fā)現(xiàn)是人類理解海洋環(huán)流的一個突破，改變了人們對洋流的傳統(tǒng)看法。因此，長時間序列全球范圍的渦旋研究對于分析和理解海洋能量和質(zhì)量輸運以及全球渦旋的時空變異性具有重要意義并且成為近幾年物理海洋學(xué)領(lǐng)域研究的重點。

Dong 等[6]根據(jù)使用數(shù)據(jù)集的不同將渦旋識別方法分為兩大類：歐拉（Eulerian）方法和拉格朗日（Lagrangian）方法。歐拉數(shù)據(jù)是指時刻快照數(shù)據(jù)或空間場數(shù)據(jù)，拉格朗日數(shù)據(jù)是指水團或物質(zhì)粒子的軌跡數(shù)據(jù)[7]。其中歐拉方法可以細分為OW（Okubo-Weiss）參數(shù)法[8–9]、纏繞角法（Winding-Angle，WA）[10]和基于渦旋幾何特征的渦旋檢測方法（Vector-Geometry，VG）[11]等方法。在歐拉方法的基礎(chǔ)上，發(fā)展出了一些渦旋追蹤算法。例如Liu 等[12]運用此方法，利用1993–2010年的高度計地轉(zhuǎn)流異常數(shù)據(jù)，對北太平洋海域的6000 多個渦旋進行了追蹤研究；Faghmous 等[13]也曾基于1993–2014年的衛(wèi)星高度計數(shù)據(jù)，共提取了4500萬個中尺度渦旋和約330萬個渦旋軌跡，并建立了1993–2014年的渦旋識別及追蹤數(shù)據(jù)集。

在流場可視化方面，Jobard 和Lefer[14]于1997年提出了一種用于二維穩(wěn)定場的均勻放置流線的算法，并且在2000年將這種算法拓展到多分辨率不穩(wěn)定場的可視化上[15]，之后基于此種方法，出現(xiàn)了多種改進算法[16–17]；Weiskopf 等[18]于2005年提出了一種結(jié)合了粒子和紋理的時空連續(xù)可視化框架，用于時變流場可視化，他們使用了Pighin 等[19]提出的徑向基函數(shù)（Radial Basis Functions，RBFs）來維持線條的動態(tài)均勻分布；何玨等[20]和Tian 等[21]將Weiskopf 等[18]的可視化框架和徑向基函數(shù)運用到海洋流場中。

在傳輸函數(shù)方面，其被廣泛應(yīng)用在空間標量場可視化中，例如Kindlmann 和Durkin[22]利用二維傳輸函數(shù)進行空間標量場可視化，他們通過計算三維標量場的一階和二階方向?qū)?shù)與標量值構(gòu)造標量場的統(tǒng)計直方圖來控制可視化效果；Sereda 等[23]通過構(gòu)建體素局部極大值和極小值的LH 統(tǒng)計直方圖來進行標量場數(shù)據(jù)邊界的可視化。矢量場方面，傳輸函數(shù)也有一些應(yīng)用，例如Helgeland 和Andreassen[24]于2004年實現(xiàn)了基于紋理方法的三維流場可視化，并且將可交互的傳輸函數(shù)應(yīng)用其中。

傳統(tǒng)的渦旋識別方法得到的識別結(jié)果多數(shù)是以統(tǒng)計圖表或者圖像的形式顯示[2,5]，這些方式顯示的結(jié)果可能會造成部分原始信息的缺失，并且不利于觀察渦旋的時空連續(xù)運動情況。傳統(tǒng)的渦旋識別與追蹤都是分開進行的，無法將二者有機的結(jié)合起來。同時其不提供交互方法，無法按用戶想法調(diào)整數(shù)據(jù)的顯示方式和效果。在這種情況下，能夠恢復(fù)甚至增強數(shù)據(jù)的整體結(jié)構(gòu)和具體細節(jié)的時空連續(xù)可視化就顯得尤為重要?，F(xiàn)有的研究中，時空連續(xù)的中尺度渦旋可視化方面的研究還很少。

為解決上述問題，本文結(jié)合二維流線可視化技術(shù)與傳輸函數(shù)技術(shù)提出了3 種方法來實現(xiàn)渦旋運動軌跡及其二維表面流動情況的時空連續(xù)可視化。這3 種方法分別是：基于OW 參數(shù)的渦旋可視化方法、基于柵格模板的渦旋可視化方法和基于矢量模板的渦旋可視化方法。這3 種方法的創(chuàng)新之處在于：（1）實現(xiàn)了二維流線可視化技術(shù)與渦旋識別與追蹤技術(shù)的結(jié)合；（2）利用線性插值算法來解決不同時間數(shù)據(jù)之間連續(xù)可視化的平滑過渡的問題；（3）結(jié)合可交互的傳輸函數(shù)實現(xiàn)流場數(shù)據(jù)的交互顯示。

2 數(shù)據(jù)集和處理方法

2.1 MSLA 數(shù)據(jù)集

本文所使用的MSLA-UV 速度場紋理數(shù)據(jù)集的類型是Ssalto/Duacs 網(wǎng)格化多任務(wù)高度計產(chǎn)品，MSLA-UV數(shù)據(jù)分為two-sat-merged[25]與all-sat-merged[26]兩種類型，本文所使用的是延時的all-sat-merged 類型的數(shù)據(jù)，其融合了多顆衛(wèi)星的數(shù)據(jù)，與two-sat-merged 類型數(shù)據(jù)相比，數(shù)據(jù)質(zhì)量較高。數(shù)據(jù)格式為NetCDF，空間分辨率為0.25°×0.25°，時間分辨率為1 d。使用前將數(shù)據(jù)進行格式轉(zhuǎn)換，將NetCDF 格式轉(zhuǎn)換為TIFF（Tadded Image File Format）格式（圖1）。

圖1 TIFF 格式的MSLA-UV 數(shù)據(jù)（2016年1月1 日）Fig. 1 MSLA-UV data in TIFF format（January 1, 2016）黃色區(qū)域為無效值區(qū)域，其他不同顏色表示流速差異The yellow area is the invalid value area, and other different colors indicate the flow rate difference

2.2 渦旋柵格模板數(shù)據(jù)

本文利用Faghmous 等[13]提出的渦旋識別算法來生成渦旋柵格模板數(shù)據(jù)，此方法首先對海面高度異常（Sea Level Anomaly，SLA）柵格數(shù)據(jù)進行5×5 的鄰域搜索，尋找出鄰域中的極值作為種子點從而找到渦心，然后從渦心向外迭代SLA 等高線，直到包含另一個渦心時停止（其迭代SLA 等高線的步長為0.05 cm，渦旋邊界內(nèi)包含的像素數(shù)目不小于4），從而識別出渦旋。最終結(jié)果以灰度圖的形式呈現(xiàn)出來，圖2 所示為2016年1月1 日35.0°～47.8°N，141.5°～169.5°E 黑潮延伸體的一部分。

圖2 柵格模板數(shù)據(jù)（2016年1月1 日）Fig. 2 Grid template data (January 1, 2016)黑色為無效值區(qū)域（大陸與北極冰雪覆蓋地或非渦海域），淺灰色為冷渦，深灰色為暖渦Black represents the invalid value area (continent and Arctic snowcovered areas or non-vortex sea areas), light gray represents cold vortex,and dark gray represents warm vortex

2.3 渦旋矢量模板數(shù)據(jù)

2.3.1 識別數(shù)據(jù)

渦旋矢量模板識別數(shù)據(jù)集是基于Liu 等[12]提出的算法而得到的。該方法通過半徑、振幅、渦核、封閉SLA 等值線等來表征渦結(jié)構(gòu)。其將全球海平面異常數(shù)據(jù)分塊，用并行計算方法對這些區(qū)域進行同時識別，之后再將識別結(jié)果拼接成一張全球渦旋圖。

在可視化過程中我們需要對兩張數(shù)據(jù)紋理進行插值，來平滑由于數(shù)據(jù)變化帶來的視覺突變。但是Liu 等[12]的方法生成的渦旋邊界上的節(jié)點數(shù)目會受到渦旋大小和形狀的影響，從而導(dǎo)致每個渦旋的節(jié)點數(shù)不同，并且不均勻。這個特點并不利于我們對渦旋進行可視化，因此，我們對渦旋邊緣的節(jié)點坐標進行了重構(gòu)，重構(gòu)示意圖如圖3 所示。

圖3 渦旋頂點坐標重構(gòu)示意圖Fig. 3 Reconstruction of eddy vertex coordinates

我們以渦心為中心建立笛卡爾坐標系，從X軸方向開始順時針將渦旋分成36 個扇形，即每個扇形的頂角均為10°，從而形成渦旋邊界上的37 個頂點，其中第一個頂點與最后一個頂點重合。對所有的渦旋都進行這樣的操作，這樣每個渦旋都會被重建，并且在同一軌跡上相鄰兩天的渦旋可以建立時間連續(xù)的頂點到頂點的幾何關(guān)系。從而保證了矢量化識別數(shù)據(jù)在插值過程中的時空連續(xù)性。

2.3.2 渦旋追蹤數(shù)據(jù)集

本文中使用的渦旋追蹤算法來自Sun 等[27]提出的混合追蹤算法。這種混合算法兼顧了渦旋的物理屬性，如距離、面積和振幅大小，同時又兼顧了渦旋的幾何屬性，如渦旋傳播過程中邊界形狀的變化。追蹤的結(jié)果如圖4 所示，圖中的紅色折線是渦心所經(jīng)過的路徑，也就是渦旋運動所形成的軌跡，該折線上的黑色圓點即為該條軌跡上每一天識別出的渦旋的渦心所在的位置，而不同顏色的多邊形邊界是通過渦旋識別得到的渦旋每一天的邊界。

3 二維流場可視化

3.1 時空框架

流場可視化是科學(xué)可視化領(lǐng)域的一個經(jīng)典研究方向，其方法主要可以分為4 類：圖標法、幾何法、紋理法、拓撲法。其中幾何法將離散幾何對象置于海洋流場中，其特征反映海流的基本性質(zhì)，適合我們實現(xiàn)時空連續(xù)的流場可視化，所以我們選擇幾何法來可視化流場。

圖4 渦旋追蹤結(jié)果示意圖Fig. 4 Schematic diagram of eddy tracking results紅色折線是渦心所經(jīng)過的路徑，折線上的黑色圓點為渦心位置The red line is the path through the eddy center, and the black dots on the red line are the position of the eddy center

在時空流場可視化方法中，跡線和流線是兩個重要的概念。在時變流場中，跡線用于描述流場中一個粒子在某個時間段的流動軌跡，而流線描述時變流場在某一時刻任意一點處速度的切線方向，對于流場中的特定位置，某一時刻有且僅有一條流線通過該點（除奇點外）。該點處流線的切線方向即表示該點處速度場的方向[28]。用流線表達的流場中，我們可以看到流場中每個點的速度信息，但是無法獲得粒子在流場中的運動情況。而用跡線表達的流場中，我們能表達粒子的運動情況，但又不能看到流場的速度信息。Weiskopf 等[18]提出了一種不穩(wěn)定場的時空連續(xù)可視化框架，參考其觀點，我們構(gòu)建的框架如圖5 所示。黑色實心圓點代表粒子，跡線xpath用藍色虛線表示，跡線與瞬時空間切片的交點用灰色圓點表示，紅色實線表示流線xstream，t0、t表示兩個瞬時空間切片。在此框架中，我們將流場以數(shù)據(jù)的時間分辨率為間隔切分成若干個瞬時空間切片，每一個瞬時空間切片都可以看作一個二維穩(wěn)定流場。在瞬時空間切片上進行流線積分來表示流場結(jié)構(gòu)，并且通過粒子在不同切片中的軌跡移動來表現(xiàn)時間相關(guān)性。

圖5 二維矢量場時空連續(xù)框架示意圖（修改自文獻[20]）Fig. 5 Schematic diagram of space-time continuum frame of two-dimensional vector field (modified from reference[20])

在二維歐式空間中，不穩(wěn)定場可以用映射u(x,t)來表示，每點的值代表在t時刻、x位置上的速度矢量值，跡線xpath由下面的常微分方程決定[29]：

在每一時刻的空間切片中，流場都可以看成一個穩(wěn)定場，流線可以更好地表達穩(wěn)定場的結(jié)構(gòu)，流線是下面方程的解：

τ與τ0和t不同，它只是控制曲線積分的一個參數(shù)，決定了曲線的精細程度，與真實的物理時間沒有關(guān)系。積分可以求得方程的解：

τ<τ0表示我們采用后向積分方式積分流線。這樣能確保粒子位于流線的頭部，獲得較好的可視化效果。

3.2 自適應(yīng)分辨率的流線可視化

在本文的可視化方法中，粒子除了具有位置信息之外，還具有年齡信息，其記錄了粒子初始生成到當(dāng)前時刻的時間。該年齡的初始值為0，如果初始化時粒子被播撒到數(shù)據(jù)范圍以外（大陸或者北極冰蓋上），年齡會被標記為“死亡”狀態(tài)。粒子有一個設(shè)定的年齡上限，粒子年齡達到該上限時也會變成“死亡”狀態(tài)?！八劳觥睜顟B(tài)的粒子不會參加流線積分，并會基于Simplex 噪聲算法[30]以隨機的方式重新播撒到視口中。每有一個粒子死亡，就會有一個新粒子生成，所以在視口中的粒子數(shù)量是一定的，流線的數(shù)量也是恒定的。如圖6 所示，當(dāng)相機向地球拉近時，地球表面的數(shù)據(jù)區(qū)域單位面積內(nèi)被更多的粒子可視化，可視化范圍變小但分辨率提高；當(dāng)相機拉遠時，相同大小的數(shù)據(jù)范圍被更少的粒子可視化，分辨率降低。這樣就實現(xiàn)了基于視口的自適應(yīng)分辨率的可視化表達。

3.3 密度控制

隨著流場的作用，粒子可能匯聚到一起，導(dǎo)致某些地方非常稀疏而某些地方非常密集，這時就需要控制粒子密度來獲得更好的顯示效果，本文采用根據(jù)自適應(yīng)徑向基函數(shù)（RBF）[18–19]計算出來的粒子影響域來控制粒子密度，即

式中，λi、xi、x、?i分別為第i個粒子的權(quán)重、中心位置、所求點位置和徑向基函數(shù)。

圖6 粒子密度隨相機變化示意圖（視口內(nèi)播撒的粒子數(shù)目相同）Fig. 6 Schematic diagram of particle density changing with the camera（the number of particles in the viewport are the same）

使用式（6）作為影響域的徑向基函數(shù)。對于一個瞬時空間切片，影響域I(x,t)為

如圖7a 和圖7b 所示，單個粒子的影響域在徑向基函數(shù)的作用下從中心到邊緣不斷降低。粒子聚集時，聚集區(qū)某一位置會受到多個粒子的影響，影響域值相互疊加。根據(jù)3?σ 法則，將高斯函數(shù)的標準差σ 設(shè)置為R/3，R是徑向基函數(shù)的半徑，設(shè)定為概率密度圖中粒子直徑的像素距離的一半， λi設(shè)為1。將這些影響域為R的粒子渲染到紋理，就生成了一張概率密度圖，如圖7c 所示為概率密度圖的一個局部，在流場的作用下，一些地方粒子匯集產(chǎn)生高值區(qū)域，一些地方粒子稀疏，甚至產(chǎn)生了黑色的0 值區(qū)域。

圖7 單個粒子的影響域（a），兩個粒子的影響域疊加（b）和局部的概率密度圖（c）Fig. 7 The influence domain of a single particle (a), the superposition of the influence domain of two particles (b), and the local probability density diagram (c)

3.4 基于GPU 的算法實現(xiàn)

一個PASS 指的是一個渲染流水線，可以包含多種著色器和一些相關(guān)的狀態(tài)設(shè)置，它是一個基本的功能單元。流場可視化的實現(xiàn)采用多PASS 技術(shù)，步驟如圖8 所示。第一步是初始化粒子分布，在時空框架中，使用可以移動的粒子作為生成流線的種子點，在全球范圍內(nèi)均勻播撒，在1°×1°的格網(wǎng)內(nèi)部播撒一顆固定位置的粒子和一顆隨機位置的粒子（在格網(wǎng)內(nèi)隨機取位置），總共播撒360×180×2 個粒子，粒子隨后被傳入VBO（Vertex Buffer Object）中，進而傳入到PASS1中的頂點著色器。

圖8 二維流線可視化流程圖Fig. 8 2D streamline visualization flowcharta.坐標紋理；b.概率密度紋理；c.顏色表示年齡信息的密度紋理，粒子年齡為0 時，顏色為白色，年齡在[0,1]時為黃色，年齡在[2,3]時顏色為紅色；d.逐步積分生成流線， τ0為 起點， τ為終點；e.地球表面生成的流線；f.地球表面生成流線的局部放大圖a. coordinates texture; b. probability density texture; c. color represents the age information of density texture; when the particle’s age is 0, the color is white;when the age is [0,1], it is yellow; when the age is [2,3], it is red. d. step-by-step generation of streamlines, τ0 and τ are the beginning and end of integral time;e. streamlines on the earth’s surface; f. localmagnification of streamlines on the earth’s surface

此過程的同時要渲染一張坐標紋理，如圖8a 所示，該紋理的R 和G 通道分別為歸一化到[0,1]的經(jīng)度和緯度?？梢酝ㄟ^采樣這張紋理的R、G 通道的值，獲得傳入粒子的經(jīng)緯坐標，進而通過經(jīng)緯坐標采樣流場數(shù)據(jù)紋理中該粒子處的流場速度值。

之后在PASS1 中的頂點著色器中進行粒子追蹤：通過4 階龍格庫塔算法計算出粒子的位移，賦予粒子新的位置并增加年齡，追蹤后的粒子信息（位置信息、年齡信息）有兩個去向：(1)傳入PASS1 的片元著色器中，在此計算粒子影響域，生成一張概率密度圖；(2)通過Transform Feedback 技術(shù)傳入PASS2 中。

PASS1 中 通 過OpenGL 的Transform Feedback 技術(shù)傳入的粒子信息在PASS2 中繼續(xù)進行粒子追蹤，不過我們同時利用PASS1 中生成的概率密度圖進行粒子年齡的控制，進而控制粒子密度。根據(jù)粒子的位置，對概率密度圖進行采樣，得到概率密度值，并根據(jù)此值更改年齡。高密度區(qū)域的粒子年齡迅速增加，快速死亡，低密度區(qū)域的粒子年齡緩慢增長，存活時間較長。圖8c 表示年齡在自然增長和概率密度圖的雙重控制下的變化情況。新撒入的粒子年齡為0，我們將其編碼為白色，隨著時間的推移，年齡增長到1，顏色變?yōu)辄S色，年齡超過2 時，年齡為紅色，紅色的粒子在不久后將會死亡。

最后一步就是根據(jù)式（4）采用龍格庫塔四階積分反向積分生成流線，如圖8d 所示，在PASS3 中的幾何著色器中繪制流線，通過調(diào)整頂點數(shù)量N和積分步長

?τ，我們可以修改流線的長度和精度。

還要說明的一點是：在每個PASS 中都會進行MSLA-UV 速度場紋理數(shù)據(jù)集插值。因為數(shù)據(jù)的時間分辨率為1 d，每一天就會有一張數(shù)據(jù)紋理，繪制是連續(xù)進行的，所以如果要進行長時間序列的不穩(wěn)定場可視化，需要在每幀之間進行數(shù)據(jù)插值，來使得流場平滑變化，獲得較好的可視化效果。利用下式對流場紋理進行線性插值：

式中，Speedinter、Speedfirst、Speedsec分別為當(dāng)前幀插值后的中間紋理、前一天的紋理、后一天的紋理，timeinter是線性插值參數(shù)。最終產(chǎn)生的繪制結(jié)果如圖8e 和圖8f所示，可以看到可視化效果很好，粒子分布均勻，海洋現(xiàn)象如海流、渦旋也能夠較好的表達。

4 渦旋可視化

4.1 渦旋可視化算法流程

本文基于二維流線可視化方法實現(xiàn)了3 種渦旋可視化方法：基于OW 參數(shù)的渦旋可視化方法、基于柵格模板的渦旋可視化方法和基于矢量模板的渦旋可視化方法，以下簡稱這3 種方法分別為OW 方法、柵格方法和矢量方法。這3 種方法的基本思想就是在二維流線可視化的基礎(chǔ)上，通過3 種數(shù)據(jù)模板來控制流線的繪制范圍，繪制渦旋范圍內(nèi)的流線，拋棄渦旋范圍之外的流線，進而實現(xiàn)渦旋的可視化。算法流程如圖9 所示。

圖9 OW 方法（a）、柵格方法（b）和矢量方法（c）的可視化流程圖Fig. 9 Visual flow charts of OWmethod (a), gridmethod (b) and vectormethod (c)與流線可視化主要區(qū)別在PASS3 中的藍色部分Themain difference from streamline visualization is in the blue section of PASS3

4.1.1 基于OW 算法的渦旋可視化方法

OW 算法是一種基于物理參數(shù)的算法，被廣泛應(yīng)用到中尺度渦的研究中[31–32]，OW 算法參數(shù)W由該點的剪切變形率（ss）、正交變形率（sn）和相對渦度（ω）定義，即

各分量的計算方法為

U、V分別代表緯向和經(jīng)向上的分量，可以通過地轉(zhuǎn)關(guān)系由SLA 梯度計算得到：

式中，g是重力加速度；f是科氏力參數(shù)。

渦旋一般存在于W值為負且以旋轉(zhuǎn)為主的流場中[33]。Chelton 等[4]和Nieto 等[34]都曾使用固定的W值作為閾值來研究海洋渦旋，本文基于OW 的可視化中同樣也指定固定的W閾值來實現(xiàn)渦旋的可視化效果。

OW 方法的具體流程如圖9a 所示，前兩個PASS和二維流線可視化算法的實現(xiàn)細節(jié)相同，不同之處在第3 個PASS 中，在幾何著色器中實時計算W參數(shù)，并將其作為判定渦旋的依據(jù)，指定W閾值，在幾何著色器中繪制條件范圍內(nèi)的流線幾何，拋棄范圍之外的流線幾何，進而實現(xiàn)渦旋的可視化。

圖10 顯示了著色器中插值與W參數(shù)計算的過程。圖10a 和圖10b 分別代表當(dāng)前幀前一天和后一天的不穩(wěn)定場紋理，每一個格網(wǎng)代表一個像素，格網(wǎng)中箭頭表示的是速度矢量。根據(jù)公式（8）對流場紋理進行插值，利用公式（9）計算每個像素的W值，并設(shè)置一個閾值W0，我們只保留在W

圖10 OW 方法插值過程示意圖Fig. 10 A schematic diagram of the interpolation process in OWmethod

4.1.2 基于柵格模板的渦旋可視化方法

相似地，柵格方法也是基于二維流線的可視化方法進行的，流程圖如圖9b 所示，此方法的輸入數(shù)據(jù)是利用Faghmous 等[13]提出的渦旋識別算法生成的渦旋柵格模板數(shù)據(jù)，將TIFF 格式的模板數(shù)據(jù)作為輸入紋理引入到PASS3 中，判定依據(jù)不再是W閾值而是模板紋理值。圖11 顯示了插值的過程，格網(wǎng)表示模板紋理中的像素單元，藍色填充的地方表示渦旋存在的區(qū)域，白色填充的地方表示不存在渦旋的區(qū)域。圖11a和圖11b 分別表示第一天與第二天同一位置的兩個渦旋，由于渦旋的運動，它的形態(tài)發(fā)生了改變，插值公式與式（8）相似，不過公式中的速度值替換成模板紋理值。然后得到如圖11c 中所示的插值結(jié)果，在有效區(qū)域內(nèi)繪制流線，實現(xiàn)渦旋可視化。

4.1.3 基于矢量模板的渦旋可視化方法

圖11 柵格方法插值過程示意圖Fig. 11 A schematic diagram of the interpolation process in gridmethod

圖12 追蹤數(shù)據(jù)與識別數(shù)據(jù)的關(guān)系示意圖Fig. 12 Schematic diagram of the relationship between tracking data and identifying data

矢量方法中，首先要對數(shù)據(jù)進行預(yù)處理，然后讀取識別和追蹤數(shù)據(jù)集，通過索引號在它們之間建立連接。如圖12 所示，渦旋識別算法將每天的渦旋識別出來并編碼相應(yīng)的日期（Day）和渦旋編號（index），Day 指的是數(shù)據(jù)日期，渦旋編號index 是Day 天第index 號渦旋。追蹤數(shù)據(jù)中存儲了很多條渦旋軌跡，每條軌跡都由一系列渦旋節(jié)點組成，節(jié)點信息中存儲著識別數(shù)據(jù)中的索引號和日期，通過每個渦旋節(jié)點的信息去識別數(shù)據(jù)中尋找相應(yīng)的渦旋數(shù)據(jù)，從而獲得渦旋的頂點坐標信息。上述這一過程在CPU 中進行。GPU中的處理過程如圖9c 所示，相較于前兩種方法，矢量方法需要一個額外的繪制流程來繪制渦旋結(jié)構(gòu)，在此流程中，通過SSBO（Shader Storage Buffer Object）傳入數(shù)據(jù)，在幾何著色器中繪制渦旋幾何，然后渲染到紋理，將紋理數(shù)據(jù)傳入PASS3 中作為判定依據(jù)，實現(xiàn)渦旋的可視化。

在繪制渦旋幾何的過程中也需要對渦旋的頂點坐標進行插值，如圖13 所示，利用下式獲得插值結(jié)果，

式中，coordinter、coordprev、coordafter分別為插值后的頂點坐標、前一天數(shù)據(jù)頂點坐標、后一天數(shù)據(jù)頂點坐標；timeinter為線性插值參數(shù)。

在實際處理過程中，因為數(shù)據(jù)量較大，CPU 數(shù)據(jù)預(yù)處理的過程和繪制渦旋幾何的過程耗時較大，所以為了減少耗費的時間，對流程進行了修改：將數(shù)據(jù)預(yù)處理和渦旋幾何紋理生成的流程預(yù)先進行，然后直接將生成的紋理圖像保存下來，之后像柵格方法那樣，將保存下來的紋理作為模板直接傳入PASS3，最后在渦旋區(qū)域內(nèi)部繪制流線，實現(xiàn)渦旋可視化。

圖13 重構(gòu)數(shù)據(jù)的頂點坐標插值示意圖Fig. 13 Interpolation diagram of vertex coordinates of reconstructed data

4.2 傳輸函數(shù)

傳輸函數(shù)是一組定義了數(shù)據(jù)值及其相關(guān)屬性與顏色、不透明度等視覺元素之間映射關(guān)系的函數(shù)[15]。不透明度決定了顯示哪些特征，顏色定義了如何顯示這些特征。

一維傳輸函數(shù)的數(shù)學(xué)表達式一般為：T:x→{c,α},x∈Rn。定義域x為自變量，代表數(shù)據(jù)的屬性值；值域為顏色 c和不透明度 α，顏色具有3 個通道，分別為R、G、B，代表紅色、綠色和藍色。

用戶可以通過設(shè)定傳輸函數(shù)，突出顯示數(shù)據(jù)的某種特定屬性或者某種屬性的某段值域進而突出顯示流場特征，并且能夠解決顯示雜亂的問題，提高可視化的效果。本文分別用速度、渦度和W值作為傳輸函數(shù)自變量來進行可交互的可視化。

以流速為自變量的傳輸函數(shù)為例說明其對顯示效果控制的基本原理：首先對將要進行渦旋可視化的流速數(shù)據(jù)進行統(tǒng)計，生成一張流速分布直方圖。以這張直方圖的流速分布作為指導(dǎo)，選取適合的屬性范圍來生成一張一維紋理，在PASS3 的幾何著色器中，以計算出的屬性值采樣這張一維紋理獲得顏色值和不透明度，然后將獲得的顏色值和不透明度值傳入片元著色器中進行染色。

如圖14 所示，通過在界面中設(shè)置Key 值點來調(diào)節(jié)生成的紋理和顏色。圖14a 中把速度為0m/s 的地方映射為淡藍色，速度大于或等于3m/s 的地方映射為紅色，中間的顏色值依據(jù)Key 值點來設(shè)置，不透明度全部設(shè)置為1，兩個Key 值點之間的顏色和透明度通過算法平滑過渡。圖14b 在圖14a 的基礎(chǔ)上調(diào)節(jié)了Key 點的不透明度，最終產(chǎn)生的一維紋理供幾何著色器中紋理采樣使用。

5 渦旋可視化效果及性能對比

5.1 渦旋可視化效果

本文的可視化效果是基于自研渲染引擎實現(xiàn)的，引擎代碼使用C++、Qt 和OpenGL 實現(xiàn)，利用著色器語言GLSL 控制顯卡進行圖形繪制。

圖14 傳輸函數(shù)交互實現(xiàn)原理示意圖Fig. 14 Transfer function interaction principle diagram

圖15 分別展示了3 種方法的可視化效果，第一行、第二行和第三行分別為使用OW 方法、柵格方法和矢量方法實現(xiàn)的效果，其中OW 方法使用的閾值W0為?2.0×10?10。第一列是使用W參數(shù)作為傳輸函數(shù)自變量的效果圖，第二列是使用渦度作為傳輸函數(shù)自變量的效果圖，第三列的傳輸函數(shù)自變量為流速。從算法角度上來講，OW 方法屬于基于物理參數(shù)的識別方法，其優(yōu)點在于理論基礎(chǔ)明確，并且參數(shù)求解簡單[35]，其弊端在于嚴重依賴閾值控制識別效果，并且不同區(qū)域的最優(yōu)閾值可能不同[34]。同時此方法對噪聲比較敏感[36]，并且局限于渦旋核心區(qū)[37]，可能對渦旋面積和半徑有所低估。另外，此方法有渦旋探測過量的趨勢[11]。從圖15 的OW 方法效果圖可以看出，OW 方法存在一些沒有形成渦旋的短線，這可能是數(shù)據(jù)中存在的噪聲，由于其附近像素的W值在所設(shè)定的閾值之內(nèi)，所以進行了后向積分，被繪制了出來，但這些短線可能并不屬于任何一個渦旋。圖15 中對比其他兩種方法，OW 方法得到的渦旋邊界更小，局限于渦旋核心區(qū)。在圖16 中展示了局部渦旋的放大效果圖，可以看到OW 方法繪制的渦旋并不完整，這可能是因為渦旋內(nèi)部W值分布不均勻和閾值設(shè)定不適當(dāng)導(dǎo)致的。同時其識別出來的渦的數(shù)量較其他兩種方法多，但是結(jié)合流場來看，有些可能不是渦旋，而是識別出的某些洋流轉(zhuǎn)彎處的殘片；柵格方法和矢量方法所用的渦旋識別數(shù)據(jù)是分別用Faghmous 等[13]和Liu 等[12]的識別算法得到的。這兩種算法都基于SLA 等值線數(shù)據(jù)，通過從渦心向外迭代等值線來確定渦旋邊界，所以這兩種方法識別出的渦旋更加飽滿，但是其可能導(dǎo)致渦旋的錯誤分類（渦旋極性與相對渦度的對應(yīng)符號不一致）或者過高估計渦旋的大小[2]。從可視化角度來講，柵格方法由于使用的數(shù)據(jù)為柵格模板，數(shù)據(jù)分辨率為0.25°×0.25°，在放大多倍后，數(shù)據(jù)的像素點相應(yīng)的被放大，所以渦旋邊緣呈鋸齒狀。而且柵格方法在進行柵格插值的時候只是簡單的柵格像素線性插值，沒有考慮相鄰兩張模板中的渦旋與渦旋之間的相關(guān)性。同時Fagmous 等[13]的算法限制最小渦的大小為4 個像素，如果前一天的某個渦旋比4 像素大（?。?，后一天比4 像素?。ù螅?，這就會導(dǎo)致有些渦旋突然消失（出現(xiàn)）。矢量方法由于我們對其識別出的渦旋進行了重構(gòu)，所以其顯示精度要比柵格模板方法高，同時Liu 等[12]的算法將最小渦的大小限制在8 個像素，所以其繪制的渦旋完整、飽滿，相對于前兩種方法效果最好，但同時可能漏掉一些小渦。從圖15d、圖15e、圖15g、圖15h 可以看出，通過自變量為W值和渦度值的傳輸函數(shù)產(chǎn)生的可視化效果并不是非常令人滿意，這主要是因為柵格方法和矢量方法采用的渦旋識別算法主要依據(jù)SLA 等值線形成的閉合幾何來確定渦旋邊界，與物理參數(shù)并不是十分匹配；但是在圖15c、圖15f、圖15i 中，速度的傳輸函數(shù)都能獲得很好的顯示效果。

圖17 展示了兩張基于傳輸函數(shù)的交互效果圖，圖17a 中，把傳輸函數(shù)設(shè)置成流速在0～0.5m/s 時透明度快速增加，0.5m/s 之后透明度緩慢增加，顏色從橘黃色緩慢過渡到深藍色。在圖17b 中，淡化了流速小于0.75m/s 的流線，流速大于0.75m/s 的部分透明度快速增加。通過交互可以把傳輸函數(shù)調(diào)節(jié)成任何形式來表現(xiàn)渦旋內(nèi)部流場，為渦旋研究提供便利。

圖15 3 種方法繪制的渦旋可視化結(jié)果（2016年1月1 日）Fig. 15 Visualization renderings of the threemethods (January 1, 2016)第一、二、三行分別為OW 方法、柵格方法、矢量方法繪制的結(jié)果，第一、二、三列的傳輸函數(shù)分別為W 值、渦度和速度。傳輸函數(shù)中W 參數(shù)和渦度的橫坐標值都被放大了1010 倍。渦度值的正號表示反氣旋渦，負號表示氣旋渦；W 參數(shù)值大于0 的部分被移除，將氣旋渦的W 參數(shù)的絕對值取代了W 參數(shù)大于0 的部分The first, second and third rows are drawn by OWmethod, gridmethod and vectormethod, respectively. The transmission functions of the first, second and third columns are W value, vorticity and velocity, respectively. The x-coordinate values of W parameter and vorticity in the transfer function aremagnified 1010 times. The positive sign of vorticity valuemeans anti-cyclonic and negative represents cyclonic. The part where the value of the W parameter is greater thanzero is removed, and the part where the W parameter is greater than zero is replaced by the absolute value of the W parameter of cyclonic

5.2 性能測試

性能測試基于以下計算機硬件配置：Windows 1064 位操作系統(tǒng)，12 G 內(nèi)存，Intel Core i5-4430 處理器（3.00 GHz），NVDIA GeForce GTX1050Ti 顯卡（4 G 顯存），視口大小為1200×600，概率密度圖、柵格模板數(shù)據(jù)、矢量模板數(shù)據(jù)的分辨率都為4096×2160，單位為像素。利用每秒繪制的幀數(shù)（Frames Per Second，F(xiàn)PS）作為性能測試指標。使用2016年1月1 日的MSLAUV 數(shù)據(jù)、柵格模板數(shù)據(jù)、矢量模板數(shù)據(jù)作為測試數(shù)據(jù)，OW 方法閾值設(shè)定為?2×10?10。

圖16 OW 方法（a）、柵格方法（b）、矢量方法（c）繪制的渦旋可視化效果局部放大圖Fig. 16 Local visualization renderings of OWmethod (a), gridmethod (b) and vectormethod (c)

圖17 傳輸函數(shù)交互效果圖Fig. 17 Transfer function interaction diagram

圖18 表示當(dāng) ?τ為1 時，每條流線上頂點數(shù)目（PointNum）對性能的影響，從圖中可以看出隨著頂點數(shù)目的增加，3 種方法的繪制效率都有較大幅度的降低。其主要原因是由于流線在幾何著色器中逐點生成，隨著點數(shù)增加，計算量會變大，繪制效率就會逐漸降低。

圖18 頂點數(shù)目對繪制效率的影響Fig. 18 The effect of the number of vertices on rendering efficiency

圖19 表示積分步長（?τ）對性能的影響，頂點數(shù)目為0 時，隨著積分步長的增大，繪制效率也有一定的降低，但是同頂點數(shù)目對性能的影響相比，積分步長對性能的影響更小。

圖19 積分步長對繪制效率的影響Fig. 19 The influence of integral step on drawing efficiency

總體來講，OW 方法具有最好的繪制效率，在3 種方法平均FPS 最高。主要的原因可能為設(shè)定的閾值使得渦旋范圍較小，需要繪制的渦旋范圍內(nèi)的流線較少，計算量較小，所以繪制效率會較高。柵格方法和矢量方法繪制效率差距大概在4～5 幀每秒，但是由圖18 和圖19 可見兩種方法在頂點數(shù)目和 ?τ不斷增加的情況下，有趨近于相同的繪制效率的趨勢。

6 結(jié)論

為了解決渦旋連續(xù)可視化的時空連續(xù)性問題，本文在前人的基礎(chǔ)上提出了3 種策略實現(xiàn)時空連續(xù)可視化，分別采用了3 種不同的渦旋識別手段：OW 參數(shù)、柵格渦旋模板數(shù)據(jù)和矢量渦旋模板數(shù)據(jù)來實現(xiàn)渦旋的判定，并結(jié)合二維流場可視化技術(shù)實現(xiàn)長周期的中尺度渦旋連續(xù)可視化?？傮w上講，基于OW 參數(shù)的渦旋可視化方法具有最高的繪制效率，但是其顯示結(jié)果不盡如人意，閾值的設(shè)置需要經(jīng)驗判定，并且很多渦旋無法完整顯示，某些局部還存在雜亂的短線擾亂視覺?；跂鸥衲０宓臏u旋可視化方法中規(guī)中矩，繪制效率低于OW 方法，但高于矢量方法，顯示效果上相對于OW 方法，渦旋顯示更加飽滿，雜亂少，但是相對于矢量方法顯示精度不高，放大后渦旋邊緣會出現(xiàn)鋸齒狀，同時由于最小渦旋限定為4 個像素，所以面對小渦時可能出現(xiàn)視覺上的突然出現(xiàn)或消失現(xiàn)象?；谑噶磕０宓臏u旋可視化方法有最好的繪制效果，渦旋完整，由于對渦旋進行了矢量化，所以邊界更加平滑，其缺點就是繪制前我們需要大量時間來進行數(shù)據(jù)預(yù)處理。通過引入傳輸函數(shù)，使得用戶可以通過與界面交互來控制流場屬性信息的顯示，使得研究者可以直觀的看到渦旋周圍的流場屬性信息，為研究渦旋提供一種新的工具。

當(dāng)前的工作中還有許多不盡如人意的地方，比如顯示效果和繪制性能不能兼得的問題。我們在未來的工作中將不斷提升渦旋的繪制性能和顯示效果，并且當(dāng)前方法是基于二維數(shù)據(jù)進行的，在未來的工作中，三維數(shù)據(jù)的渦旋可視化是一個巨大的挑戰(zhàn)。