基于聚類的流面自動布局及生成算法

唐 燁，解利軍，桂立業，何麗莎，鄭 耀

浙江大學 航空航天學院，杭州 310027

1 引言

流場可視化是科學可視化的經典分支之一。通過直觀的圖像來表示抽象的流場數據[1]，既可以加快從數據中獲取信息的速度，又能讓不具備領域專業知識的人理解這些信息。如今在航空航天、氣象、自動化等領域的科學研究過程中，流場可視化已經成為不可或缺的技術。

流面可視化是流場可視化的一種方式，近幾年來越來越流行[2]，相比于常用的流線法，流面法更適合表現復雜的流場結構[3]，并且展現更豐富的流場局部信息[4]。

1.1 相關概念

流場是一種矢量場。在空間坐標系(x,y,z)中，任意點(x0,y0,z0)都有唯一確定的向量(Vx,Vy,Vz)與之對應，就說在空間(x,y,z)范圍內存在矢量場V。在流場中，該矢量場即速度場。

流線是一個無質量的質點（種子點）在某一穩定流場中以某點為起點，在該流場中運行的軌跡。流線上的每一點都與流場中該點對應的速度相切。假設V(x)是一個三維向量場，在時刻t經過點x0的流線為L(x0,t)，其中L(x0,t)滿足公式（1）：

類似于流線的概念，流面是某條種子線C的所有種子點以C中的位置為起點在穩定流場中運行的軌跡的集合，設該流面為S，S滿足公式（2），流面S中的點處處與流場中該點處的速度相切[5]。

其中S(s,t)是種子點C中s處的種子點在流場中的運動軌跡，也就是由C中s處的種子點生成的流線，S為C中所有點生成的流線的集合，如圖1所示。

圖1 流面示意圖

1.2 相關工作

Hultquist首先在1992年提出了一種高效的流面生成算法——前沿推進算法。其基本思想為相鄰的兩條流線之間形成流帶，每條流帶需要記錄兩個前沿位置點，算法通過推進每條流帶的前沿位置的方式來構造流面，其特別之處為前沿的個數并不是固定的，當某條流帶前沿距離過大時可以通過加入新點分裂流帶來提升流面精確度，當前沿過窄時可以考慮合并相鄰流帶來減少計算與存儲資源的消耗[5]。

前沿推進法非常適用于流場中存在大量的流線分散與流線聚攏的情況，即流場中速度大小基本穩定，但方向出現偏差的情況，并且也較易于實現。但當流場中的速度大小分布不勻稱時，前沿推進法就難以適用。為此Garth在2004年提出了一個基于Hultquist算法的新算法，其基本思想為，在前沿推進時前進的距離以弧長為參數而不是固定的時間參數，同時在分裂與合并流帶時考慮曲率，使得最后得出的圖像能夠更精確[6]。

Scheuermann提出了一種在四面體網格中生成流面的算法[7]，其基本思想為通過重心坐標在每個四面體內計算出規則表面，通過面上的線段組成最終流面。這種方法能夠有效地利用流場中的拓撲結構信息，但用這種方法最終生成的流面質量十分依賴于網格的密度[7]。

與以上這些顯示計算曲面的算法不同的是，Wijk在1993年提出過一種通過定義網格邊緣的標量函數的方式來隱式地計算曲面的算法[8]，但由于它無法將所有可能的曲面都計算出來，因此應用范圍并不廣。

Schneider等人在2009年提出在構造曲面的時候用施密特插值法取代線性插值法[9]，其計算所需的額外數據都可以在流線生成過程中通過計算得出，可以將曲面的精確度提升到4階。

Schulze等人提出過一種強制流面前沿線與流場中的速度垂直的方法，它可以在湍流區生成高質量的流面[10]。

以上都是確定種子線情況下所用的流面生成算法，但是當拿到流場數據時，無法確定在哪里布線能夠得到較好的結果。因此本文介紹一種種子線自動布局以及流面生成的算法。

本文的主要思想：流場中不同區域的流動情況并不一致。如圖2所示，以流線法表示，可能有直線區域、曲線區域、渦流區域?？梢酝ㄟ^聚類算法將這些區域進行劃分。劃分后的每個聚類可以用一個流面表示。生成流面的方法為，首先從分區中心點出發，沿著流場的曲率場在該分區范圍內積分得到種子線。選擇在曲率場中積分是因為：（1）曲率場中任意點的方向與流場中該點處的速度方向垂直，而垂直于流場速度方向的種子線生成的流面能夠覆蓋較大的范圍。（2）沿著曲率場生成的種子線所生成的流面更能表現出流場的流動結構。得到種子線后，通過前沿推進法配合施密特插值法生成光滑流面。

圖2 流場不同區域示意圖

2 流場的分區和種子線布局

2.1 程序流程

該算法處理的對象是流場數據V(X)∈R3，其中V∈R3，X∈R3，X代表流場空間中點的坐標，V(X)代表點X處的速度。

圖3為程序的流程圖，對于流場中每個點P，首先算出它的曲率|C|以及速度模的梯度G用于之后的聚類以及種子線生成，詳見2.2節。隨后可以用改進版的K-means算法將流場中的點集劃分為K個類，詳見2.3節，其中K的值由用戶自己輸入。然后對于每一個類的聚類中心點，用四階龍格庫塔法在曲率場中進行積分得到K條種子線，詳見第3章。對每條種子線，通過帶施密特插值的前沿推進法生成并繪制流面，詳見第4章。如果需要，用戶可迭代該過程，修改參數再次進行聚類與流面生成工作。

圖3 程序流程圖

2.2 流場中數據點間相似度計算

在三維空間點集中，兩點間的相似度一般用它們的歐式距離來表示，但在流場數據集中，每個點除了坐標信息外，還有對應的速度信息。速度有兩個基本屬性：大小和方向。但由于關注的是潛在的種子線位置，因此希望把種子線布置在能夠代表周圍流線（速度）發生變化的區域。因此，并不是直接對流場中不同點間的速度進行比較，而要比較速度的變化情況。

同時，渦流位置也可以通過曲率值的大小來反映，而流場的結構也與渦流的位置有關[11]。

在流場空間中某處布置種子點形成流線S，則通過計算可以得到S中任意點P處的曲率|C|[12]，曲率也可以量化流線的彎折程度，也就是速度方向的變化程度。其中

其中V為點P處的速度，

點P處速度大小的變化情況|G|可以用速度模的梯度來表示，其中

某點處速度的變化情況可以通過|C|與|G|的值來反映，則流場中點P1、P2間的距離Dis(P1,P2)可以定義為Dis(P1,P2)=|C1-C2|+|G1-G2|。

當P1、P2兩點坐標間的歐式距離較大，而它們速度的變化量又相同的時候，Dis(P1,P2)為0，在聚類時會被歸為一類，如果這兩點的速度變化量較大，則它們應該代表兩個潛在的布線區域。因此判定相似度時需要考慮兩點間的歐式距離：Dis(P1,P2)=|C1-C2|+|G1-G2|+|X1-X2|。

流場數據中空間坐標與速度值的取值范圍并不相同，因此相比于點坐標間的距離值，曲率值或梯度值間的距離值可能會顯得很小。此時兩點間的相似度無異于兩點坐標間的歐氏距離，因此需要將它們全部歸一化處理。

在對不同的流場數據進行聚類時，三個參數所占的權重未必一樣，因此需要對三個參數進行加權，最終兩點間相似度Dis(P1,P2)表示為：

其中，α為兩點間相似度中二者的歐式距離所占的份額，β為速度變化量的兩個參數中曲率所占的份額，0<α<1，0<β<1。

2.3 基于改進的K-means算法進行流場聚類

聚類是將對象的集合按照一定的標準分割成多個相似子對象集合的過程，同一類中的對象間相似度較高，不同類間對象間相似度較低。

聚類算法可以分為兩大類：基于層次與基于劃分。在此基礎上還有三小類[13]：基于密度、基于模型、基于網格。

K-means算法屬于基于劃分的聚類法，與其他方法相比，有著實現簡潔與運行快速兩大優點。因此本文選擇該方法進行流場聚類。

傳統K-means算法中初始聚類中心的選取完全是隨機的，而經過實驗發現K-means算法的效果會受到初始聚類中心的影響。選取不當的初始聚類中心會使聚類過程更加漫長[14]。

因此，針對K-means算法在初始聚類中心選取上的盲區進行一定的改進。如果K個初始聚類中心互相之間具有較高的相似度，那么它們最后很有可能屬于同一個聚類，這會影響后續聚類的時間，于是希望K個初始聚類中心間的相似度盡可能得低。其過程如下：

（1）在數據集中隨機選擇一個初始種子點。

（2）對于數據集中所有點P，記錄其與聚類中心集合中的點Ci的距離的最小值D(P)。

（3）選擇一個新點加入到聚類中心集合中，選取的原則為D(P)值越大，被選中的概率也越大。

（4）如果得到K個聚類中心，執行下一步，否則，返回到第（2）步。

（5）對于數據集中所有點，將其劃分到與它相似度最高（距離最?。┑木垲愔行乃鶎俚木垲愔?。

（6）對于每一個聚類，選擇到類中所有點的距離和最低的點作為新的聚類中心。

（7）如果本次迭代得到的K個聚類中心與上次迭代得到的K個聚類中心相同，則聚類情況已經穩定，結束程序，否則，返回到第（6）步。

實驗表明，采用改進版的K-means算法能夠在參數一致的前提下，減少聚類過程中的迭代次數，減少聚類所需時間，如表1所示。

表1 兩種K-means算法比較

聚類效果如圖4所示，其中圖4（a）為 α=0.5、β=0.3、k=2時的聚類效果。圖4（b）為 α=0.8、β =0.8、k=4時的聚類效果，圖4（c）為 α=0.5、β =0.8、k=6時的聚類效果。

圖4 聚類效果圖

3 種子線生成

得到了流場的聚類與中心數據后，就可以著手種子線的生成工作，流場的聚類中心雖然能代表流場的聚類，但單憑一個點難以確定種子線，還需確定種子線生成的方向與長度。

通過觀察可以得知，與速度方向平行的種子線無法生成流面，與速度方向垂直的種子線生成的流面能覆蓋較大范圍[15]，如圖5所示，圖中藍色帶箭頭曲線代表流線方向，x為垂直流場方向種子線，紅色曲線為不與流場垂直的種子線。二者生成同一流面，但種子線X更短。

圖5 種子線示意圖

種子線的朝向需要考慮到流場局部位置的結構信息，否則，生成流面效果將大打折扣。如果流面處處與曲率方向相切并且平行于旋轉軸，那么它可以表現出螺旋流的特性。如果該面垂直于曲率方向和旋轉軸，該流面就難以傳達出足夠的有效信息。如圖6所示，圖6（a）為以垂直于速度與曲率方向生成的種子線所形成的流面，圖6（b）為垂直于速度方向平行于曲率向量方向積分形成的種子線生成的流面，可以看到圖6（b）中的內容包含了圖6（a）的主要內容。

圖6 種子線垂直曲率與平行曲率效果圖

根據流線中的曲率定義，曲率向量C由流線的速度向量V與加速度向量a通過叉積得出，因此C向量的方向一定是與V垂直，即給定一個初始點P，在常流場生成曲率場中的運動軌跡必定處處與流場中該點處的速度垂直。因此，最終種子線生成方法是：以聚類中心點P為初始點，在曲率場中以生成流線的方法通過積分生成種子線。

由于每個聚類中心只能代表各自所屬的聚類。因此，種子線沿著曲率場積分到達分區邊界時，停止積分。

積分使用龍格庫塔法（Runge_Kutta）來實現。

四階龍格庫塔法為：

種子線效果如圖7所示，其中圖7（a）為α=0.5、β=0.3、k=2時生成的種子線，圖7（b）為 α =0.8、β =0.8、k=4時生成的種子線，圖7（c）為 α=0.5、β =0.8、k=6時生成的種子線。

圖7 種子線效果示意圖

4 流面生成

得到種子線后就可以開展流面的生成工作，基本的思想采用前沿推進法，如圖8（a）所示。前沿推進法將流面劃分為若干流帶，相鄰兩條流線之間構成一條流帶，記每條流線當前最新采樣點為前沿，每個流帶對應一個追蹤器記錄其兩個前沿。追蹤器通過單鏈表的方式連接，如圖8（b）所示，流帶通過前沿在流場中推進的過程來生成。當流場中前沿處速度為0，或者前沿位置超出流場空間時，該流帶結束生成。

圖8 三角化法基本思想示意圖

流帶通過三角形面片法表示，即在相鄰的兩條流線的采樣點上生成三角形條帶來近似化表現流帶，三角形的形狀需要盡量規則，角度不宜過大，這可以通過局部貪心法，即每次新增對角線較短的三角形，來實現。兩條流線的前沿在流場中各生成一個新的采樣點后，與追蹤器中記錄的前沿點一起就可以構成一個四邊形，選擇較短對角線與兩前沿連線構成三角形面片加入流帶中，并更新追蹤器中的前沿信息。其過程可用如下偽代碼表示。

輸入：以單鏈表形式相互連接的追蹤器Tracer，流場數據V，三角形面片集合Output

輸出：無

Advance_Font（Tracer）

1.Left_advanced=false；

2.while（1）

3.L0=Tracer->L0，R0=Tracer->R0；

4.在V中通過龍格庫塔法以L0、R0計算L1、R1

5.L_diag=|L1R0|，R_diag=|L0R1|

6.M_diag=min（L_diag，R_diag）；

7.if（M_diag==L_diag）

8. ad_left=true；

9.else ad_left=false；

10.if（Left_advanced&&ad_left）return；

11.if（ad_left）

12.Output->push_triangle（L0，R0，L1）；

13.Tracer->L0=L0=L1；

14.Left_advanced=true；

15.else

16.Output->push_triangle（L0，R0，R1）；

17.Advance_Font（Tracer->next）；

18.Tracer->R0=R0=R1；

隨著流場中流線的分散與聚攏，流帶的寬度也會產生變化，當流線不斷聚攏時，流帶的寬度會變窄，置之不理的話會生成多余的三角形，影響到繪制的時間，此時可以考慮將其與相鄰流帶合并。當流線不斷分散時，流帶的寬度會增大，置之不理的話，會影響到生成曲面的精確度，此時可以考慮將該流帶分裂。

判斷合并的標準：如圖9（a）所示，流帶Ribbon1當前的兩前沿分別為L0、M0。其相鄰流帶Ribbon2的兩前沿為M0、R0。通過L0、M0、R0三點積分得到的采樣點為L1、M1、R1。定義四邊形L0L1R1R0的高H為L1與R1兩點到線段L0R0所屬直線的距離的較小值。當點L1與點 R1間的距離小于 H，且 L0、L1、M0、M1、R0、R1六點近似共面時，將流帶Ribbon1與流帶Ribbon2合并。

圖9 流帶合并分裂示意圖

合并流帶的方法：如圖9（a）所示，在流帶Ribbon2中，選擇對角線M0R1，將?M0R0R1加入到流帶中，更新右前沿為 R1，在流帶Ribbon1中，選擇對角線L1M0，將?L0M0L1加入到流帶中，更新追蹤器1左前沿為L1，將?L1M0R1加入到流帶中，將追蹤器1右前沿更新為R1，并將追蹤器1右邊的追蹤器更新為追蹤器2右邊的追蹤器。其過程如下偽代碼所示：

輸入：非隊尾追蹤器Tracer，流場數據V，三角形面片集合Output

輸出：無

Merge_Ribbon（Tracer）

1.L0=Tracer->L0，M0=Tracer->R0；

2.R0=Tracer->next->R0；

3.在V中通過龍格庫塔法以L0、M0、R0計算L1，M1，R1

4.Output->push_triangle（L0，L1，M0）；

5.Output->push_triangle（M0，R0，R1）；

6.Output->push_triangle（L1，M0，R1）；

7.Tracer->L0=L1，Tracer->R0=R0；

8.Tmp=Tracer->next

9.Tracer->next=Tracer->next->next；

10.Delete tmp；

判斷分裂的標準：當新采樣點L1與點R1的距離|L1R1|大于四邊形 L0L1R1R0的高 H的兩倍，即|L1R1|≥2H時，流帶需要分裂。

分裂流帶的方法：如圖9（b）所示，在流帶中用適當的方法加入新點M1，將?L0M1R0加入流帶中，生成新流帶Ribbon2，其左前沿為M1，右前沿為Ribbon1的右前沿R0，更新Ribbon1的右前沿為M1。將Ribbon2右邊的流帶設置為Ribbon1右邊的流帶，Ribbon1右邊的流帶設置為Ribbon2。其過程如下偽代碼所示：

輸入：追蹤器Tracer，流場數據V，三角形面片集合Output

輸出：無

Split_Ribbon（Tracer）

1.L0=Tracer->L0，R0=Tracer->R0

2.通過龍格庫塔法在V中以L0、R0計算L1、R1

3.計算新加入點M1的坐標

4.Output->push_triangle（L0，M1，R0）；

5.Tracer->R0=M1；

6.New Tracer2；

7.Tracer2->L0=M1，Tracer2->R0=R0；

8.Tracer2->next=Tracer->next；

9.Tracer->next=Tracer2；

加入新點M1的坐標在取線段L1R1中點的基礎上配合施密特插值法使得最后生成的流面更平滑。其基本思想如圖10所示。

圖10 施密特插值法示意圖

在插入新點時考慮到流場沿著種子線方向的變化，決定點M1是位于中點偏上還是偏下，其計算方法為：

其中，，由于是基于中值法改進，通常r取0.5。

流面生成效果如圖11所示，直接三角化法得出的圖像由于無法適應流線變化會產生失真，如圖11（a）所示。前沿推進法產生的圖像不夠平滑，如圖11（b）所示。配合施密特插值法則能得到平滑又精確的流面，如圖11（c）所示。

圖11 流面效果圖

5 實驗結果分析

測試平臺為Win 10操作系統，Intel CORE i7，8 GB內存，顯卡為NVIDIA GEFORCE GTX 860M。

圖12展現了通過自定義種子線與通過本文算法自動布線在Delta-wing流場數據中生成流面的情況：圖12（a）為連接端點（0.334 5，－200，40.404）與（0.334 5，200，40.404）為種子線生成的流面，圖12（b）為連接端點（201.003，－200，90.909）與（201.003，200，90.909）為種子線生成的流面。可以看到它們都屬于比較穩定的流面，且值覆蓋了流場中一小部分，難以對流場整體結構進行描繪。圖12（c）是用本文方法以 K=2，α=0.5，β=0.3生成的流面，可以看到它畫出了兩塊近似對稱的區域。圖12（d）是用本文方法以 K=4，α=0.5，β =0.8生成的流面，可以看到在兩塊對稱區域的基礎上增加了兩塊平滑曲面，圖12（e）是用本文方法以 K=6，α=0.8，β=0.8生成的流面，可以看到相對于圖12（d）中圖像，新增了兩塊平滑曲面，它們找出了流場中變化較劇烈的部分，并生成了有代表性的圖像，并且由于圖12（e）中的圖像相比圖12（d）中的圖像并無劇烈的變化，不再追加聚類個數。另外，用Esturo的方法在Telda-wing上所得圖像[15]，只能生成一張流面，容易給人造成流場中只有該區域有流動的錯覺，從而對流場整體的狀況產生誤判。

圖13展示了通過自定義種子線與本文方法在5Jet流場數據中所生成流面的情況。圖13（a）為連接端點（1，64，90）與端點（128，64，90）為種子線所生成的流面，圖13（b）為連接端點（1，64，64）與端點（128，64，64）為種子線所生成的流面?？梢钥闯觯@兩張流面都只能反映出流場局部的流動情況，無法對該流場整體的流動情況作出判斷。圖13（c）是本文方法以 K=2，α=0.3，β=0.3生成的流面。圖13（d）是本文方法以K=4，α=0.5，β=0.8生成的流面，保持整體結構的同時增加了平滑曲面。圖13（e）是以 K=6，α=0.5，β =0.3生成的流面，保持圖13（d）中整體結構的同時新增了兩張平滑曲面。至此，認為5JET數據的整體流動結構已得出，不再增加聚類個數。

圖12 Delta-Wing實驗結果圖

圖13 5JET實驗結果圖

6 結論

本文設計了一個流場種子線自動布局以及流面自動生成的算法，通過坐標、曲率以及梯度參數對數據點進行聚類，得到中心點，在曲率場中通過四階龍格庫塔法得到種子線，然后基于前沿推進法得到流面，能夠有效地在未知流場數據中生成表現力較強的流面，也能為人工布線位置提供較好的參考。目前，該算法還需要用戶自己輸入聚類個數K，希望將來能夠開發出自動確定聚類個數的算法。該算法依然是串行執行的版本，未來將實現并行計算的版本，以提高程序運行的速度。此外，該算法所生成的流面依然存在流面間遮擋的問題，將來希望開發出能夠根據觀察者視角與流面遮擋情況自動調節局部透明度的版本。

[1]宋漢戈，劉世光.三維流場可視化綜述[J].系統仿真學報，2016（9）：1929-1936.

[2]Edmunds M，Laramee R S，Chen G，et al.Surface-based flow visualization[J].Computers&Graphics，2012，36（8）：974-990.

[3]Mcloughlin T，Laramee R S，Peikert R，et al.Over two decades of integration-based，geometric flow visualization[J].Computer Graphics Forum，2010，29（6）：1807-1829.

[4]Laramee R S，Garth C，Doleisch H，et al.Visual analysis and exploration of fluid flow in a cooling jacket[C]//Proceedings IEEE Visualization 2005，2005：623-630.

[5]Hultquist J P M.Constructing stream surfaces in steady 3D vector fields[C]//IEEE Conference on Visualization（Visualization’92），1992：171-178.

[6]Garth C，Tricoche X，Salzbrunn T，et al.Surface techniques for vortex visualization[C]//Symposium on Visualization，Konstanz（Vissym 2004），Germany，May 2004：155-164.

[7]Scheuermann G，Bobach T，Hagen H，et al.A tetrahedrabased stream surface algorithm[C]//IEEE Visualization 2001，San Diego，CA，USA，2001：151-553.

[8]Wijk J J V.Implicit stream surfaces[C]//IEEE Conference on Visualization（Visualization’93），1993：245-252.

[9]Schneider D，Wiebel A，Scheuermann G.Smooth stream surfaces of fourth order precision[J].Computer Graphics Forum，2009，28（3）：871-878.

[10]Schulze M，Germer T，R?ssl C，et al.Stream surface parametrization by flow-orthogonal front lines[C]//Computer Graphics Forum，2012：1725-1734.

[11]Banks D C，Singer B.Vortex tubes in turbulent flows：Identification，representation，reconstruction[C]//IEEE Conference on Visualization（Visualization’94），1994：132-139.

[12]Roth M.Automatic extraction of vortex core lines and other line type features for scientific visualization[microform][D].Swiss Federal Institute of Technology Zurich，2000.

[13]Rokach L.A survey of clustering algorithms[M]//Data Mining and Knowledge Discovery Handbook.Boston，MA：Springer，2009：269-298.

[14]Arthur D，Vassilvitskii S.k-means++：The advantages of careful seeding[C]//Eighteenth ACM-SIAM Symposium on Discrete Algorithms（SODA 2007），New Orleans，Louisiana，USA，January 2007：1027-1035.

[15]Esturo J M，Schulze M，R?ssl C，et al.Global selection of stream surfaces[J].Computer Graphics Forum，2013，32（2pt1）：113-122.