唐 燁，解利軍，桂立業(yè)，何麗莎，鄭 耀

浙江大學(xué) 航空航天學(xué)院，杭州 310027

1 引言

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

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

1.1 相關(guān)概念

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

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

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

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

圖1 流面示意圖

1.2 相關(guān)工作

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

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

Scheuermann提出了一種在四面體網(wǎng)格中生成流面的算法[7]，其基本思想為通過重心坐標(biāo)在每個四面體內(nèi)計算出規(guī)則表面，通過面上的線段組成最終流面。這種方法能夠有效地利用流場中的拓撲結(jié)構(gòu)信息，但用這種方法最終生成的流面質(zhì)量十分依賴于網(wǎng)格的密度[7]。

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

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

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

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

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

圖2 流場不同區(qū)域示意圖

2 流場的分區(qū)和種子線布局

2.1 程序流程

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

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

圖3 程序流程圖

2.2 流場中數(shù)據(jù)點間相似度計算

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

同時，渦流位置也可以通過曲率值的大小來反映，而流場的結(jié)構(gòu)也與渦流的位置有關(guān)[11]。

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

其中V為點P處的速度，

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

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

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

流場數(shù)據(jù)中空間坐標(biāo)與速度值的取值范圍并不相同，因此相比于點坐標(biāo)間的距離值，曲率值或梯度值間的距離值可能會顯得很小。此時兩點間的相似度無異于兩點坐標(biāo)間的歐氏距離，因此需要將它們?nèi)繗w一化處理。

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

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

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

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

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

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

傳統(tǒng)K-means算法中初始聚類中心的選取完全是隨機的，而經(jīng)過實驗發(fā)現(xiàn)K-means算法的效果會受到初始聚類中心的影響。選取不當(dāng)?shù)某跏季垲愔行臅咕垲愡^程更加漫長[14]。

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

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

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

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

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

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

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

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

實驗表明，采用改進版的K-means算法能夠在參數(shù)一致的前提下，減少聚類過程中的迭代次數(shù)，減少聚類所需時間，如表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 種子線生成

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

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

圖5 種子線示意圖

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

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

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

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

積分使用龍格庫塔法（Runge_Kutta）來實現(xiàn)。

四階龍格庫塔法為：

種子線效果如圖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）所示。前沿推進法將流面劃分為若干流帶，相鄰兩條流線之間構(gòu)成一條流帶，記每條流線當(dāng)前最新采樣點為前沿，每個流帶對應(yīng)一個追蹤器記錄其兩個前沿。追蹤器通過單鏈表的方式連接，如圖8（b）所示，流帶通過前沿在流場中推進的過程來生成。當(dāng)流場中前沿處速度為0，或者前沿位置超出流場空間時，該流帶結(jié)束生成。

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

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

輸入：以單鏈表形式相互連接的追蹤器Tracer，流場數(shù)據(jù)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；

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

判斷合并的標(biāo)準(zhǔn)：如圖9（a）所示，流帶Ribbon1當(dāng)前的兩前沿分別為L0、M0。其相鄰流帶Ribbon2的兩前沿為M0、R0。通過L0、M0、R0三點積分得到的采樣點為L1、M1、R1。定義四邊形L0L1R1R0的高H為L1與R1兩點到線段L0R0所屬直線的距離的較小值。當(dāng)點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，流場數(shù)據(jù)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；

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

分裂流帶的方法：如圖9（b）所示，在流帶中用適當(dāng)?shù)姆椒尤胄曼cM1，將?L0M1R0加入流帶中，生成新流帶Ribbon2，其左前沿為M1，右前沿為Ribbon1的右前沿R0，更新Ribbon1的右前沿為M1。將Ribbon2右邊的流帶設(shè)置為Ribbon1右邊的流帶，Ribbon1右邊的流帶設(shè)置為Ribbon2。其過程如下偽代碼所示：

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

輸出：無

Split_Ribbon（Tracer）

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

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

3.計算新加入點M1的坐標(biāo)

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的坐標(biāo)在取線段L1R1中點的基礎(chǔ)上配合施密特插值法使得最后生成的流面更平滑。其基本思想如圖10所示。

圖10 施密特插值法示意圖

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

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

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

圖11 流面效果圖

5 實驗結(jié)果分析

測試平臺為Win 10操作系統(tǒng)，Intel CORE i7，8 GB內(nèi)存，顯卡為NVIDIA GEFORCE GTX 860M。

圖12展現(xiàn)了通過自定義種子線與通過本文算法自動布線在Delta-wing流場數(shù)據(jù)中生成流面的情況：圖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）為種子線生成的流面?？梢钥吹剿鼈兌紝儆诒容^穩(wěn)定的流面，且值覆蓋了流場中一小部分，難以對流場整體結(jié)構(gòu)進行描繪。圖12（c）是用本文方法以 K=2，α=0.5，β=0.3生成的流面，可以看到它畫出了兩塊近似對稱的區(qū)域。圖12（d）是用本文方法以 K=4，α=0.5，β =0.8生成的流面，可以看到在兩塊對稱區(qū)域的基礎(chǔ)上增加了兩塊平滑曲面，圖12（e）是用本文方法以 K=6，α=0.8，β=0.8生成的流面，可以看到相對于圖12（d）中圖像，新增了兩塊平滑曲面，它們找出了流場中變化較劇烈的部分，并生成了有代表性的圖像，并且由于圖12（e）中的圖像相比圖12（d）中的圖像并無劇烈的變化，不再追加聚類個數(shù)。另外，用Esturo的方法在Telda-wing上所得圖像[15]，只能生成一張流面，容易給人造成流場中只有該區(qū)域有流動的錯覺，從而對流場整體的狀況產(chǎn)生誤判。

圖13展示了通過自定義種子線與本文方法在5Jet流場數(shù)據(jù)中所生成流面的情況。圖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生成的流面，保持整體結(jié)構(gòu)的同時增加了平滑曲面。圖13（e）是以 K=6，α=0.5，β =0.3生成的流面，保持圖13（d）中整體結(jié)構(gòu)的同時新增了兩張平滑曲面。至此，認為5JET數(shù)據(jù)的整體流動結(jié)構(gòu)已得出，不再增加聚類個數(shù)。

圖12 Delta-Wing實驗結(jié)果圖

圖13 5JET實驗結(jié)果圖

6 結(jié)論

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

[1]宋漢戈，劉世光.三維流場可視化綜述[J].系統(tǒng)仿真學(xué)報，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.