基于矢量完全信息熵的流線并行分布方法

申麗銘，郭雨蒙，王文珂＊

1.國防科技大學，氣象海洋學院，湖南 長沙 410003

2.軍事科學院軍隊政治工作研究院，北京 100091

引言

矢量場可視化是重要的科學可視化方法，可以幫助科研工作人員從大量數(shù)據(jù)中觀測與提取關鍵特征。流線可以簡潔地表示流場中流體的運動趨勢，且計算量少、易于觀測，是矢量場可視化中最常用的可視化方法之一[1]。

流線分布方法決定了矢量場可視化表達的效果，主要分為兩類：基于均勻分布的流線分布與基于特征的流線分布方法。采用均勻分布的方法在生成流線數(shù)目較少時，容易丟失矢量場中的重要特征。而基于特征的流線分布方法，將流線集中分布在矢量場的關鍵區(qū)域附近，有利于科研人員分析矢量場的重要特征。作為一種成功的信息度量理論，信息論為矢量場特征提供了定量描述工具。基于信息論的流線分布方法使得大多數(shù)流線被放置在具有復雜信息的區(qū)域附近，突出其顯著特征。然而，現(xiàn)有的基于信息論的流線分布方法，僅考慮了矢量的方向分量，這將忽略某些矢量幅度變化劇烈的現(xiàn)象，例如CFD 流場[2]中的激波。

本文提出了一種基于矢量完全信息熵的流線分布方法，該方法建立矢量方向和長度的聯(lián)合直方圖計算完全信息熵，引導流線集中放置于矢量方向或矢量大小變化迅速的高熵值區(qū)域，從而更完整地揭示矢量場的潛在信息。同時，為了提高效率，算法采用了網(wǎng)格用于流線剪枝，消除距離過近的流線，并對部分步驟進行了并行化。實驗表明，與現(xiàn)有的基于信息熵的流線分布方法相比，本文方法得到的流線分布效果既能體現(xiàn)矢量大小變化劇烈的區(qū)域，又能保留矢量方向變化的顯著特征。

1 國內外研究現(xiàn)狀

為了有效展現(xiàn)矢量場信息，往往需要生成多條流線，為此需要考慮流線如何合理布局。流線分布方法可以分為基于均勻分布的流線種子點分布與基于特征的流線種子點分布兩類。

Turk 和Banks 首先提出了一種圖像引導（Image guided）的流線種子點選擇方法[3]，嘗試生成分布均勻的流線，但是只能用于二維矢量場，且效率較低。而后，Mao 等人[4]擴展了這種圖像引導的方法，用于在曲線網(wǎng)格中的三維參數(shù)曲面上創(chuàng)建均勻分布的流線。為了解決在三維矢量場中生成流線而不考慮其在屏幕空間投影而產生的視覺混亂問題，Li 等人[5]擴展了圖像引導方法，流線的密度和渲染樣式可以根據(jù)各種標準靈活控制，可以更有效地可視化三維矢量場。

然而，均勻分布的方法在流線數(shù)目較少時會丟失矢量場中的重要特征。為此，Verma 等人[6]首先提出了基于特征的流線分布策略，可以有效反映臨界點等矢量場的重要特征。他們使用類似于拓撲分析的概念定義特征區(qū)域，并針對特征區(qū)域不同的流動模式使用不同的種子點放置模板。Chen 等人[7]定義了兩條流線的相似距離，通過阻止相似距離較小的流線生長，避免生成影響視覺效果的冗余流線。由于特征區(qū)域的流線相似距離較大，該算法產生的流線在特征區(qū)域周圍較為密集，而非特征區(qū)域的流線比較稀疏。

為了定量度量矢量場特征，近年來有學者利用信息論解決流線的分布問題。Furuya 和Itoh[8]首先將信息理論應用于流線生成，他們根據(jù)切線方向分布度量流線的復雜性，獲得其信息熵場，用于指導流線生成。Xu 等人[9]利用信息熵與特征區(qū)域的關系引導種子點放置，同時利用條件熵度量當前流線與原矢量場的信息差異，以引導信息缺失區(qū)域的進一步布種。Lu[10]提出了一種基于紋理的可視化方法，該方法根據(jù)矢量方向熵提取不同的特征結構，以類似紋理的處理方式實現(xiàn)流線的生成?？紤]到矢量大小變化劇烈的區(qū)域也可視為一類流場特征，Guo 等人[11]在其基于信息熵的流線分布方法上提出了矢量長度熵的計算方法，以揭示矢量場更多的信息。

2 基于完全信息熵的流線分布方法

本節(jié)對提出的基于矢量完全信息熵的流線分布方法進行了詳細描述。該方法不同于已有的忽略矢量長度信息的信息熵方法，同時考慮了矢量方向與長度信息，以更完整地揭示矢量場的信息。算法利用矢量完全信息熵引導流線集中放置于矢量方向或矢量大小變化迅速的高熵值區(qū)域，而后根據(jù)現(xiàn)有流線建立中間矢量場并計算其與原始矢量場的矢量完全條件熵。條件熵值可反映兩個矢量場之間的信息差異，因此可據(jù)此進一步迭代指導放置流線直至達到停止條件。

并且，為了提高效率，算法采用網(wǎng)格以利用已有流線限制新流線生成（見2.3 節(jié)），使得在流線生成過程中可以判斷流線間距離，提前終止距離過近的流線生成，加速流線剪枝過程。

2.1 算法流程

算法流程如圖1所示，包括以下步驟：

圖1 基于矢量完全信息熵的流線分布算法流程圖Fig.1 Flow chart of streamline distribution algorithm based on vector complete information entropy

(a)計算信息熵：遍歷每個數(shù)據(jù)網(wǎng)格點，建立矢量信息直方圖，計算信息熵場（詳見2.2 節(jié)）；

(b)選擇初始種子點集：計算信息熵場的局部極大值點，以局部極大值點為中心放置種子點（種子點模板參考文獻[9]），作為初始種子點集；

(c)生成流線：為了有效捕捉矢量場特征，模板布點會布置多個種子點，使得生成的流線中某些流線間距離過小，已有方法生成流線后再判斷流線間距進行流線刪除。為了加速此過程，本文在流線生成過程中利用網(wǎng)格修剪流線（詳見2.3 節(jié)）；

(d)重建矢量場：通過擴散法（Diffusion method）[13]由生成的流線重建矢量場；

(e)計算條件熵：建立重建矢量場與原始流場的聯(lián)合直方圖，并計算條件熵場（詳見2.2 節(jié)）；

(f)判斷是否終止：判斷原始矢量場和重建矢量場之間的條件熵是否滿足終止條件（小于某個預設閾值）。如果滿足或已達到算法迭代次數(shù)，則算法停止，輸出當前流線集；否則，執(zhí)行步驟(g)；

(g)選取新種子點集：對于條件熵值高于閾值的點，按條件熵值的大小賦予其重要性，采用重要性采樣的方法選取候選種子點，返回步驟(c)。

2.2 計算信息熵與條件熵

香農提出的信息熵用于解決信息的量化度量問題，他直觀地將信息熵定義為特定離散信息的出現(xiàn)概率，條件熵表示在已知隨機變量的條件下隨機變量的不確定性。

圖2 計算信息熵所用的直方圖。（a）矢量場數(shù)據(jù)；（b）矢量方向信息熵直方圖；(c)矢量長度信息熵直方圖；(d)矢量完全信息熵直方圖。Fig.2 Histogram for calculating information entropy.(a)Vector field data;(b)Histogram of vector direction information entropy;(c)Histogram of vector magnitude information entropy;(d)Histogram of vector complete information entropy.

結合公式(2-1)和(2-3)，即可計算得出給定矢量場區(qū)域的方向信息熵。

矢量方向條件熵計算與方向信息熵相似，即通過公式（2-3）得到與后，代入公式（2-2）可得矢量方向條件熵。

另外，由圖2 可知，完全信息熵直方圖區(qū)間的數(shù)量是分別采用矢量長度與矢量方向直方圖的區(qū)間數(shù)量的乘積。為了體現(xiàn)矢量長度變化對流線分布效果的不同影響，定義為矢量長度（軸）和矢量方向的區(qū)間數(shù)的比值。若的值較大，說明對矢量長度的分組更精細，矢量大小快速變化區(qū)域的熵值增加，流線分布更側重于體現(xiàn)矢量大小的變化信息。反之，若的值較小，說明矢量方向的分組更精細，則流線分布主要體現(xiàn)的是矢量方向的變化。特別地，若本文算法與現(xiàn)有的基于信息熵的流線分布方法相同。在本文實驗中，矢量方向直方圖的分組數(shù)默認設置為，矢量長度分組數(shù)則由控制。參數(shù)對流線分布的影響將在第4 節(jié)實驗部分進一步介紹。

2.3 基于Dsep 網(wǎng)格修剪流線

在生成種子點集時，為了有效捕捉矢量場特征，本文使用模板布點：二維時放置4 個，位于以信息熵極大值點為中心的菱形頂點上；三維時放置8 個，位于以信息熵極大值點為中點的正八面體頂點上。模板布點通常會導致流線在某一區(qū)域過于密集，因此修剪冗余流線是十分必要的?，F(xiàn)有的流線分布方法在流線生成結束后進行流線剪枝。通過比較任意兩條流線的距離（或相似度），以刪除距離相近（或形狀相似）的冗余流線，從而避免視覺混亂。此類剪枝過程要計算每兩條流線的最短距離，且每當生成新的流線后都需要執(zhí)行一次剪枝過程，時間開銷較大。

考慮到計算采樣點到流線距離的過程可能反復出現(xiàn)，本文通過建立一個Dsep網(wǎng)格結構，使得在流線生成過程中可以判斷流線間距離，提前終止距離過近的流線生成，加速流線剪枝過程。Dsep網(wǎng)格是一個完全包圍原矢量場的標準網(wǎng)格結構，每個網(wǎng)格單元主要存儲一個是否可用的標志，有“可用”和“不可用”兩種狀態(tài)。Dsep網(wǎng)格的作用是讓已生成流線的網(wǎng)格標記為“不可用”，新流線生成時只能通過“可用”的網(wǎng)格單元，如圖3所示。為了讓流線分布體現(xiàn)特征，即在特征附近允許容納更多的流線，一個點對其它流線的限制距離取決于其信息熵值。若某一點熵值很高，說明其重要度高，則該點附近應容納更多流線以突出顯示，因此對其它流線的間隔限制應該很低。反之若熵值很低，則間隔限制應設置較大。

圖3 給出了兩個流線的“不可用”網(wǎng)格的示例。一般情況下，若流線通過了網(wǎng)格，且該網(wǎng)格點信息熵值沒有高過一個閾值，則此網(wǎng)格的標志直接由“可用”變?yōu)椤安豢捎谩?。如圖3 中的紅色圓圈區(qū)域具有低熵值，將流線通過的網(wǎng)格置為“不可用”（圖中黃色網(wǎng)格）。但是，若流線通過的某些點的熵值高于閾值，為了讓更多流線突出表示特征，這些點所在的網(wǎng)格仍被設置為“可用”，如圖3 中臨界點O所在的網(wǎng)格。另外，對于流線通過的熵值較低（低于某給定閾值）的采樣點所在的網(wǎng)格，要求其它流線距其間隔更遠，那么除包含它的網(wǎng)格外，其相鄰的網(wǎng)格也被設置為“不可用”。為方便區(qū)分，這些因低熵值點而標記為“不可用”的相鄰網(wǎng)格在圖3 中用藍色顯示。

圖3 網(wǎng)格示意圖：黃色網(wǎng)格為流線通過的“不可用”網(wǎng)格，藍色網(wǎng)格為低熵值的“不可用”網(wǎng)格。Fig.3 Diagram of Dsep grid:yellow grids are “unavailable” grids with streamlines passing through,and blue grids are “unavailable”grids with low entropy.

本文算法中，由于流線間的距離由網(wǎng)格結構控制，流線生成時只需判斷下一個網(wǎng)格是否“可用”，若不可用則提前終止流線生成。雖然Dsep網(wǎng)格結構需要額外開銷用于計算每個點的限制間隔等操作，但仍可有效對流線剪枝過程進行加速。在第4 節(jié)的實驗中會對網(wǎng)格限制法與流線生成后剪枝法進行效率對比。

3 流線分布方法的并行實現(xiàn)

如圖1所示，本文算法流程中需要大量計算，如計算信息熵場與條件熵場，根據(jù)已有流線重建矢量場等；并且相比僅考慮矢量方向的信息熵計算，矢量完全信息熵計算過程增加了維度，引入了額外的計算負擔。為提高算法效率，本文采用了并行技術對算法進行加速。

3.1 信息熵與條件熵并行計算

在2.2 小節(jié)的熵值計算過程中，每個點的鄰域都需要進行一次信息熵與條件熵的計算，且任意兩個點的熵計算互不影響，彼此獨立，為此，采用如下的并行計算過程：

首先，將矢量場數(shù)據(jù)劃分為子塊，每個子塊關聯(lián)到一個指定的線程上，然后遍歷子塊的每個點。處理一個點時，需要先從共享內存里讀取原始矢量場或重建矢量場的鄰域數(shù)據(jù)至線程緩存中。特別的是，為了減少讀取次數(shù)，便于子塊邊界區(qū)域附近點的熵值計算，當讀取子塊數(shù)據(jù)時，邊界需要向外擴展，即多讀取一定范圍內相鄰子塊的數(shù)據(jù)。

然后，建立矢量長度與方向的信息熵直方圖或原矢量場和重建矢量場的聯(lián)合直方圖，并在每個線程獨立計算該點的信息熵或條件熵。

3.2 重建矢量場

此步驟通過對所有已生成流線使用擴散法重建矢量場，下面介紹其并行化過程。

在基于信息熵的流線分布算法[9]中，詳細描述了根據(jù)已生成流線重建矢量場的過程。該步驟以二維為例，需要最小化一個如下的能量函數(shù)：

其中，

其中：

算法 重建矢量場

指定Δt,Δx,Δy,μ的值For 第i層循環(huán)do For 每個線程e do讀取子塊數(shù)據(jù) Sub-Block(e)到線程e的緩存For Sub-Block(e)的每個點 do讀取i-1層循環(huán)的 ui-1,vi-1根據(jù)公式(3-2)、(3-3)計算ui,vi End For End For柵欄（Barrier）：等待其它線程運行結束End For

第n次迭代層中點的計算彼此獨立，而計算它們需要調用第n-1 次迭代層的各個數(shù)值，因此每一層迭代結束時必須進行同步。算法偽碼3-1 顯示了這個同步的并行過程。可以看到，該步驟通過柵欄完成數(shù)據(jù)同步，即進入下一層迭代前當前迭代層的所有線程操作必須已經完成。該步驟是對遞推公式(3-2)與(3-3)的混合并行實現(xiàn)，對同一層是并行計算，不同層則是串行，需要等待前一層所有操作完成。柵欄的存在一定程度上影響了并行加速效果，但并行方法總體上仍提高了整個算法的效率。

4 結果與分析

本節(jié)中，我們在3 個數(shù)據(jù)集上測試了基于矢量完全信息熵的流線分布方法，并與現(xiàn)有的基于信息熵的分布方法[9]進行了對比。實驗中的算法全部由C++編碼實現(xiàn)。實驗用的計算機配置6 核Intel Xeon X5675 3.06 GHz CPU、32 GB 內存、NVIDIA GeForce GTX 550 Ti 顯卡和64 位Windows 7 操作系統(tǒng)。

4.1 可視化效果分析

（1）數(shù)據(jù)集1 是簡單的人工模擬矢量場，數(shù)據(jù)規(guī)模為50*50，在矢量場中心有一個鞍點。矢量長度顏色映射后分布如圖4(a)所示，暖色表示長度大、冷色表示長度小，可見左下角區(qū)域矢量長度迅速從紅色變?yōu)樗{色和綠色，矢量大小變化劇烈。為了驗證實驗效果，每個對比實驗均生成了相同數(shù)量的流線。

圖4(b)～(d)中的背景場均為信息熵場的顏色映射，顏色從暖到冷代表熵值從高到低。三張圖中均生成了21 條流線。

圖4 為設置不同的長度與方向區(qū)間比值α（取值范圍為0、0.3、0.6）后所得的可視化效果。α=0為長度不占比重的情況，即等同于已有的基于信息熵的流線分布方法。隨著α 增大，左下區(qū)域流線數(shù)目逐步增加。當α=0.6 時，具有速度大小劇烈變化的矢量場左下區(qū)域聚集了更多的流線，特征更加明顯。因此，后續(xù)兩個實驗均取α=0.6 以突出體現(xiàn)本文效果。

圖4 數(shù)據(jù)1 的可視化實驗結果。(a)中熱力圖為矢量長度場映射結果，(b)-(d)中的熱力圖為信息熵場映射結果。Fig.4 Visualization experiment results of data 1.The heat map in(a)is the mapping result of vector magnitude field,and the heat map in(b)-(d)is the mapping result of information entropy field.

由當前流線與原流場矢量信息計算得到的條件熵值是一種流線可視化效果的量化評價方式，條件熵值越低則可認為可視化結果越還原原矢量場。當α 取值分別為0、0.3、0.6 時，整個矢量場的條件熵值分別為0.976765、0.97011、0.932224，也說明了α取0.6 時確實具有更好的可視化效果。

（2）數(shù)據(jù)集2 是一個120*120 的氣象風場，結果如圖5所示，每張圖中均生成24 條流線。各圖中位置相同的紅色方框標注出了最具有代表性的區(qū)域。

由圖5(c)可知，紅色方框區(qū)域內的矢量長度場中冷暖色變化劇烈，說明該區(qū)域中的矢量大小變化很快。因此，在以完全信息熵為依據(jù)得到的流線圖5(b)中，紅色方框區(qū)域內熵值明顯增加，出現(xiàn)流線聚集現(xiàn)象。作為對比，以矢量方向信息熵為依據(jù)的圖5(a)的紅色方框區(qū)域內分布了較少的流線，未能體現(xiàn)風速大小變化。

圖5 數(shù)據(jù)2 的可視化實驗結果。紅色方框標識出了三種方法效果的對比處，(a)-(c)中熱力圖分別為矢量方向信息熵場、完全信息熵場與矢量長度場映射結果。Fig.5 Visualization experiment results of data 2.The red box indicates the comparison of the effects of the three methods.The heat maps in(a)-(c)are the mapping results of the vector direction information entropy field,the complete information entropy field and the vector magnitude field,respectively.

（3）數(shù)據(jù)集3 是尺寸為210*180*27 的三維真實氣象風場，圖6 為本文算法的三維可視化效果，共生成45 條流線。可以看到，流線大部分聚集在方向變化劇烈的經典特征區(qū)域（渦旋、鞍點等）附近。同時，右上和左下的紅框區(qū)域矢量大小變化顯著，流線更加密集。

圖6 數(shù)據(jù)3 的三維可視化實驗結果Fig.6 3D visualization experimental results of data 3.(a)Streamline distribution diagram.(b)Complete information entropy field.

4.2 效率分析

本節(jié)對Dsep網(wǎng)格剪枝效率與計算步驟并行效率進行實驗對比分析。

為了降低因流線過于密集而產生的視覺混亂，現(xiàn)有基于信息熵的流線分布方法在生成所有流線后對冗余流線進行剪枝，本文提出的Dsep網(wǎng)格剪枝方法則是通過已有流線限制新流線生成。

上述兩種方法在數(shù)據(jù)集1 與數(shù)據(jù)集2 的執(zhí)行時間如表1 和表2所示，與生成所有流線后剪枝相比，本文提出的通過已有流線間隔（Dsep網(wǎng)格）限制流線生成過程的剪枝方法運行時間減少約10%，證明了本文提出的Dsep網(wǎng)格剪枝方法效率更高。

表1 對數(shù)據(jù)集1 采用不同剪枝策略的算法運行時間Table 1 Running time of different pruning strategies

表2 對數(shù)據(jù)集2 采用不同剪枝策略的算法運行時間Table 2 Running time of different pruning strategies for dataset 1 for dataset 2

表3 列出了本文方法的串行實現(xiàn)與不同線程數(shù)目下本文算法的并行運行時間。由表可知，相比串行實現(xiàn)，使用6 線程并行時算法運行時間可降低約40%，有效提高了算法效率。

表3 不同并行線程數(shù)量的本文算法運行時間Table 3 Running time of the algorithm in this paper with different number of parallel threads

5 結論與展望

本文從可視效果與效率兩個角度對提出的基于矢量完全信息熵的流線分布方法進行了分析。該方法同時考慮了矢量的方向和長度信息，能夠有效揭示更多的矢量場信息，且不遺漏原始數(shù)據(jù)的顯著特征。此外，通過在流線生成過程中利用Dsep網(wǎng)格進行流線剪枝，以及對算法過程中的耗時計算進行并行加速，本文也有效提高了算法效率。

由結果與分析可得，本文提出的算法在效率上仍有較大的提升空間，未來將進一步采用GPU 和并行加速技術對算法效率進行優(yōu)化。此外，視點選擇在三維矢量場可視化中也是十分重要的。通過信息論對不同視點所揭示的信息進行量化并推薦信息量豐富的視點，可以有效協(xié)助專家進行科學數(shù)據(jù)分析。下一步工作中，我們也將基于完全信息熵對視點選擇進行優(yōu)化。

利益沖突聲明

所有作者聲明不存在利益沖突關系。