邵緒強，劉藝林，楊 艷，林麗娜

(1. 華北電力大學控制與計算機工程學院，河北 保定 071003；2. 海軍裝備部裝備項目管理中心，北京 100071)

向量場數(shù)據(jù)可視化大多來自于對流場的模擬和觀察，在大部分情形下，向量場數(shù)據(jù)可看為流場數(shù)據(jù)[1]。

本文綜述從龐大的流體數(shù)據(jù)集中提取形成特征所需的流體信息的過程，即流體的特征提取過程。這些流體特征是對數(shù)據(jù)的抽象，通過只提取特征、從面向問題的角度進行數(shù)據(jù)縮減[2]，從而將可視化提升到更高的抽象級別[3]。該方法經(jīng)過較為復雜的提取特征，后續(xù)的可視化步驟將會較快地實現(xiàn)[4]。隨著計算機性能的不斷提高，基于特征和拓撲提取的流體可視化技術成為了新的研究熱點[5]。

流體特征可視化在醫(yī)學上應用廣泛。在四維核磁共振成像的心臟血流方面，K?HLER等[6]發(fā)現(xiàn)部分疾病的嚴重程度與血流的某些特性(如渦流核心的位置、旋渦形狀和初始位置)有關，并用λ2方法定義旋渦，且提取主動脈中的旋渦核心。實驗中渦核呈球體，其半徑與旋渦強度成正比。文獻中通過使用四維PC-MRI采集技術實現(xiàn)可靠、實時的血流測量，從而對患者進行血流動力學定性和定量分析。目前，醫(yī)學研究人員正在研究旋渦等特征的流動模式與不同病理的關系提取；在腦動脈瘤風險評估中，OELTZE-JAFRA等[7]發(fā)現(xiàn)血流中的旋渦與動脈瘤破裂相關，通過提取渦核線來準確繪制旋渦的觀測結果，在渦核線上標注箭頭以示血流運動方向。

流體特征可視化在海洋、天氣、宇宙等研究領域也受到極大重視。CHELTON等[8]以旋渦的形式對10年來全球海洋的海面高度場進行了分析，并用Okubo-Weiss準則定義旋渦，將旋渦分為氣旋型和反氣旋型旋渦。由于噪聲和閾值的原因，無法識別弱渦；1983年，LUGT[9]在大氣中尋找颶風，將物質顆粒圍繞一個共同中心的旋轉運動定義為旋渦；MARKOWSKI等[10]定義風暴相對螺旋度(storm relative helicity, SRH)，發(fā)現(xiàn)其在時間和空間上均為高度可變的，并提出該參數(shù)用于預測風暴方面時存在的缺陷；HADJIGHASEM和HALLER[11]提供了一種參考系不變性方法，從非定常流場中提取動態(tài)拉格朗日相干結構，并對木星及其他行星大氣中的旋渦邊界進行探測。

流體的特征提取為特征可視化的首要步驟[1]，本文主要綜述流體中旋渦特征提取方法。旋渦的非正式定義反映了人的主觀直覺[12-13]，但至今仍很難正式給出旋渦的正式定義[14-16]，因此文獻中旋渦提取的依據(jù)各不相同。目前，常用的流體旋渦特征提取大致分為基于點、線、幾何曲線和基于機器學習的方法。

(1) 基于點的方法?；跇颖军c處的物理量進行計算，如根據(jù)壓力[17]或螺旋度[18]等物理量來提取旋渦，其特征通常是流體中的一個連通區(qū)域。該方法的缺陷，如難以區(qū)分較為集中的旋渦，大多數(shù)方法需要設置閾值，微小特征漏判等。近年來提出的Ω方法[19]對閾值和實驗效果方面有突破性進展。

(2) 基于線的方法。通過搜索渦核線(粒子做旋轉運動所圍繞的線)進而提取旋渦，其共性問題是渦核線通常會被過度分割，導致實驗效果降低。

(3) 基于幾何曲線的方法。通過流線或跡線的幾何特性進行旋渦提取，其中曲率中心法和纏繞角度法，分別通過曲率趨于一點和角度累計達到2π等條件來提取旋渦區(qū)域。該方法時間復雜度較高，在三維空間的應用仍需一定的改進。

(4) 基于機器學習的方法。通過預先訓練模型等方法有效提高了旋渦提取的精度，使用卷積神經(jīng)網(wǎng)絡(convolutional neural network, CNN)等方法提取流場中的旋渦，對傳統(tǒng)方法假陽性、假陰性、時間復雜度較高等缺陷均有一定的改進。

根據(jù)參考坐標系的運動類型可以將旋渦提取方法分為伽利略不變性、旋轉不變性和拉格朗日不變性。目前最理想的方法為拉格朗日不變性坐標系下的旋渦提取，但提取效果提升的同時也存在一些問題，如數(shù)據(jù)計算量過于龐大等。

與旋渦提取領域其他綜述性文章相比，本文除了闡述經(jīng)典方法外還著重介紹了近幾年來的新型旋渦提取方法，如Ω準則、Rortex法、基于機器學習的方法等，讓研究者對其研究發(fā)展也有一定的深入了解。表1對常用的方法進行總結。

表1 常用的旋渦提取方法匯總Table 1 Summary of commonly used vortex extraction methods

1 相關概念及旋渦提取方法

1.1 流線

流場中，每個點均與速度矢量相切的光滑曲線稱為流線。在描述流場時，可用流線簇的疏密程度表示矢量的大小。定常流中流線的位置和形狀不隨時間變化，即流線式為

其中，u，v，w為各方向的分量。

SALZBRUNN和SCHEUERMANN[20]提出流線謂詞概念，通過為每條流線設置一組謂詞定義流結構，從而揭示三維流場中的旋渦問題。ROBINSON等[39]認為，若瞬時流線在垂直于渦核線平面上的投影大致呈圓形或螺旋狀，且該投影圖形與渦核線有相同的運動趨勢時提取到旋渦。在特征可視化研究初期，該方法及其簡單的特性被廣泛應用，但因自我引用，想找到旋渦需預先計算出渦核線的位置，因此這種方法現(xiàn)實意義較小[40]。CHEN等[41]提出一種將向量場中的周期軌道和分離線作為初始流線的新方法。BANKS和SINGER[42]使用渦度場流線，根據(jù)垂直于渦核的平面內的壓力最小值進行旋渦提取。GARTH等[43]將流線拓展到流面，提出流面計算和確定渦核邊界的新方法，在流體進行復雜運動時該方法也可以獲得精確結果。由于流場中存在大量流線，流線方法經(jīng)常出現(xiàn)漏選情況，ZHENG等[44]根據(jù)流線的幾何特征來判斷其類型并進行分類，該方法能有效篩選多余流線，提高可讀性。

1.2 跡線

流場中，流體粒子運動時描繪的軌跡曲線稱為跡線。在定常流中，跡線等同于流線；在非定常流中，與流線不同，跡線的位置和形狀隨時間變化。在此需要區(qū)分流線和跡線的性質，流線由同一時間步長下的各個流體質點組成，跡線由同一流體質點在不同時間步長下的位置組成。

大多數(shù)情況下，非定常流比定常流更接近自然現(xiàn)象，但是非定常流的計算復雜度較大。隨著計算機性能的不斷提高，人們也逐漸開始研究非定常流中的旋渦提取。

跡線可表示為

WEINKAUF等[28]將Sujudi&Haimes法[27]推廣到跡線，加入時間維度后，跡線表示為

根據(jù)非定常流中跡線的特征提取旋渦。

類似于流線謂詞，SALZBRUNN等[21]將跡線分為受旋渦影響的跡線和其余跡線2部分。文獻[21]定義了基于跡線謂詞的框架，根據(jù)跡線是否具有特定的屬性為每個跡線設置布爾值，因此每個謂詞可以將跡線分為兩類；結合布爾代數(shù)可以將跡線進行更精確的分類，從而對不同的跡線結構進行更精確的定義。

1.3 平行向量場

PEIKERT和ROTH[22]在1999年引入平行向量場的概念，通過設定約束C=V×W=0找出2個向量場V和W平行的點的位置，該方法相當于求出V和W叉積運算的零等值面。在二維中，返回孤立點；在三維中，返回渦核線，如圖1所示。

圖1 三維中的平行向量場[30]Fig. 1 Parallel vector field in three-dimensional[30]

平行向量場的提出為旋渦提取領域提供了更簡潔的方法，可以對大多數(shù)旋渦提取方法重新定義并簡化計算步驟。PEIKERT和ROTH[22]通過計算v||(?v)v找到局部曲率為零的點。其中，v為速度；?v為速度梯度；SCHINDLER等[30]計算平行向量場找出渦核線種子點。

ROTH和PEIKERT[29]提出了一種使用平行向量提取渦核線的方法，利用對速度矢量v進行二階導數(shù)計算，即

可實現(xiàn)對彎曲旋渦的定位。

1.4 矢量場拓撲

矢量場拓撲[45-47]由臨界點和由臨界點發(fā)出或指向臨界點的積分曲線構成，用于清晰明了地表示場信息。1989年，HELMAN和HESSELINK[47]將矢量場拓撲方法用于可視化領域，并將臨界點定義為矢量大小為零的點，通過計算雅可比矩陣特征值判斷臨界點的類型：2個特征值實部為正對應排斥狀，和負對應吸引狀，一正一負對應鞍狀；特征值具有虛部對應旋渦狀，實部的正負對應旋出和旋入，實部為零時旋轉曲線近似閉合。如圖2所示。

圖2 矢量場拓撲分類圖[47]Fig. 2 Topological classification of vector field[47]

KASTEN等[23]利用矢量場拓撲設計一種魯棒的旋渦邊界提取算法，根據(jù)局部加速度最小找到臨界點。由于加速度為伽利略不變量，適用于較多文獻的實驗場景，但是仍高度依賴參考坐標系。MAHROUS等[48]提出了一種基于Sperner定理的渦核區(qū)域提取方法，該方法觀察每個網(wǎng)格單元邊界處的向量，以此分析向量場。速度矢量在旋渦附近表現(xiàn)出特定的形式，通過搜索這些形式找到旋渦。BUJACK等[49]在矢量場中建立局部參考系，同一時間步長下在流場不同區(qū)域中使用不同的坐標系，根據(jù)雅可比矩陣行列式極值提取臨界點。在三維中，矢量場拓撲仍有很好的實用性，HELMAN和HESSELINK[50]證明當存在2個復特征值時，總會存在一個對應的實特征值，由此可在三維空間中得到焦點或中心點，從而提取旋渦。

1.5 旋渦新定義——旋渦矢量(Rortex)

由于至今仍沒有給出旋渦的明確定義，導致研究人員在各自的旋渦提取方法中的判斷依據(jù)各不相同，為初學者造成理解的混亂。

LUGT[9]對旋渦定義是：許多粒子圍繞一個共同的中心旋轉的運動。盡管其與視覺上的直接觀測結果表述類似，但不適于設計旋渦提取算法。

目前大多數(shù)旋渦提取方法將速度梯度張量(雅可比矩陣)的特征值作為判斷條件。這些基于特征值的方法存在幾個不可避免的問題，如旋轉軸的識別過程性能較差，存在剪切力的影響等。為了解決以上問題，從數(shù)學的角度出現(xiàn)了一個基于特征向量的矢量，稱為旋渦矢量[51]，其大小與垂直于旋轉軸平面中的旋轉強度有關，方向由速度梯度矩陣的實特征向量確定，表示旋渦的旋轉軸。根據(jù)旋渦矢量的定義，將旋渦定義為旋渦矢量為正數(shù)的流通區(qū)域。

由于該方法并非基于流體動力學，因此無論是不可壓縮還是可壓縮流體，旋渦矢量均可以準確地提取旋渦。該方法還有另一優(yōu)點，Rortex法中的旋渦矢量可以同時表示旋渦的旋轉強度和旋轉方向，這是眾多標量方法無法做到的。

文中提出的快速計算Rortex算法可以準確地描述局部流體的旋轉，能清晰地顯示旋渦結構。TIAN等[52]對算法進行優(yōu)化后，Rortex算法與λ2準則的計算時間大致相同。在后期研究中，GAO和LIU[24]發(fā)現(xiàn)剪切力對復特征值的虛部有一定的影響，而基于Rortex的方法可以進行準確測量。

2 流體旋渦特征提取傳統(tǒng)方法分類

2.1 基于點的方法

2.1.1 基于局部物理量滿足閾值要求提取旋渦

該類方法中，閾值的取值極大地影響了實驗效果。在早期研究中，HUNT等[17]在定常無粘性的二維流體中定義旋渦區(qū)、匯合區(qū)和流動區(qū)。通過局部壓力最小值提取旋渦區(qū)域，設置約束P≤Pthresh來找到滿足條件的區(qū)域，其中，Pthresh為人為設置的閾值，并將該低壓區(qū)域定義為旋渦區(qū)。

后期研究證明局部壓力最小不能保證存在旋渦，反之亦然，流體其他的特性也有可能引起旋渦，如螺旋度、渦度等，因此目前一般不采用壓力法進行旋渦提取，但其為旋渦提取領域提供了一個良好的開端。

螺旋度一般由速度、渦度等物理量計算[53-54]。該方法通過局部區(qū)域螺旋度達到局部區(qū)域絕對值最大值來提取旋渦。LEVY等[18]將螺旋度密度Hd、螺旋度H和單位化螺旋度Hn應用于旋渦提取中，即

其中，ω為渦度。

由式(5)可知，螺旋度密度的符號由速度矢量和渦量矢量夾角余弦的符號決定，螺旋度密度較大則表示速度和渦度之間的夾角較小。通過單位化螺旋度在最小曲率流線上達到局部區(qū)域絕對值最大值來確定渦核區(qū)域位置，由時間上反向積分找到奇異點(速度為零的點)來確定旋渦開始位置，積分得到的曲線即旋渦旋轉軸。該方法解決了由低渦度流場區(qū)域，以及速度矢量和渦度之間的角度較大(如邊界層)的高渦度低速區(qū)域引起的誤判問題。

螺旋度法的一個優(yōu)點是螺旋度歸一化到[-1,1]的范圍，用于確定渦軸方向。但是該方法無參考系不變性，因此其只能應用于定常流或固定某一幀的旋渦提取。且該方法不考慮某些特殊情況，如Hn分母為0等。

2.1.2 基于速度梯度提取旋渦

渦度作為一個伽利略不變量，與流體粒子的角速度有關，然而渦度不能區(qū)分純剪切運動與旋渦的實際旋轉運動，因此不適合用作旋渦提取。目前基于對速度梯度?v=Ω+S及其對稱部分Ω和反對稱部分S的一類方法應用廣泛。該類方法采用由HUNT[25]提出Q準則。ZHOU等[55]提出

得出P，Q，R的關系。其中，P，Q，R由速度梯度?v

導出。其中，tr為矩陣的跡；det為行列式的絕對值，且

其中，Ω為對稱部分，表示旋轉率張量；S為反對稱部分，表示應變率張量。若Q＞0，則代表旋轉力克服應變力，另外要求旋渦區(qū)域壓力低于環(huán)境壓力，由這2個條件直接提取出旋渦結構。該方法僅適用于強渦實驗環(huán)境，提取弱渦時有一定誤差。

JEONG和HUSSAIN[26]提出λ2方法(第二大復特征值方法)。速度梯度?v可分解為

由于S2+Ω2為對稱矩陣，因此只有實特征值。若S2+Ω2的3個特征值滿足第二大特征值小于0，則該位置存在旋渦，并且λ2越小，渦度越大，通過λ2局部最小確定旋渦位置。該方法的缺點是，在旋渦密集區(qū)域識別單個旋渦的難度較大。

目前大多數(shù)基于點的方法過度依賴閾值，對于同一個流體，閾值的選取會導致不同的實驗結果?；诖?，LIU等[19]提出了Ω準則來減小閾值對實驗效果的影響，后期被歸結為第三代旋渦提取方法[56]，如圖3所示。

圖3 Q準則不同閾值下的結果對比，t=6.8 T ((a) Q=0.001;(b) Q=0.0001; (c) Q=0.01)[19]Fig. 3 Comparison of results under different thresholds of Q criterion, t=6.8 T ((a) Q=0.001; (b) Q=0.0001; (c) Q=0.01)[19]

對于Q準則和λ2準則等方法，目前仍然沒有一個閾值能適合所有的旋渦場景。經(jīng)過ZHANG等[57]的多次驗證發(fā)現(xiàn)，對于Q準則實驗，Q=0.001最優(yōu)；對于λ2準則實驗，λ2=-0.001為最優(yōu)。相比而言，Ω準則不需要根據(jù)實驗效果來尋找最適合的閾值，后期研究發(fā)現(xiàn)，Ω=0.52時，對所有場景均表現(xiàn)出較好的適應性。

LIU等[19]將Ω準則定義如下：先將式(9)重寫為

設

其中，tr是矩陣的跡，另定義一個比值形式，其中，ε為趨于0的正數(shù)，確保除數(shù)非零。實驗表明，結果對Ω的取值并不敏感，一定程度上解決過度依賴閾值的問題。

該方法重新討論流體速度場分解并將其進一步分解為旋渦部分和非旋渦部分，將Ω的物理意義定義為渦度與渦核內部整體渦度的比值，表明當旋轉力強于應變力時將形成旋渦，Q準則和λ2準則對于等值面的物理意義描述則較為模糊。

2.1.3 基于點的方法的問題

大多數(shù)方法過度依賴閾值的選取，缺乏普適性；提取弱渦時容易出現(xiàn)假陰性；在參考系變換時效果較差，后期研究對此點進行改進，將Q準則和λ2等經(jīng)典方法嵌入?yún)⒖枷挡蛔冃灾校沪笢蕜t解決了閾值問題，但也存在假陽性和假陰性，ε的取值對實驗結果也有一定的影響。后期工作中，許多基于點的方法被賦予新的使用限制[58]。

2.2 基于線的方法

2.2.1 特征向量法

SUJUDI和HAIMES[27]提出在三維矢量場中提取旋渦的算法，認為存在渦核線需要滿足2個條件：①雅可比矩陣有一個實特征值和一對共軛復特征值；②滿足w=u-(u·n)n=0，其中，u為該點速度矢量，n為實特征值對應的單位特征向量。

該方法首先需要找到臨界點，即相對于觀察者速度為零的點。在該點計算雅可比矩陣特征值，若存在一個實特征值和一對共軛復特征值，即

其中，ui為該點速度矢量；Ci為插值所需參考點的速度矢量。此步驟相當于將速度投影到與實特征值對應的特征向量垂直的平面上。結果為零意味著此處粒子僅沿著特征向量方向移動而不旋轉，由此得到渦核位置。該方法為基于線的方法研究提供了良好的基礎，但是存在特征值和特征向量復雜運算，線性插值步驟存在誤差等問題；并且KENWRIGHT和HAIMES[59]發(fā)現(xiàn)特征向量法提取的特征易受其他向量特征影響，若2個旋渦距離較近，提取的旋渦結構會發(fā)生一定角度的偏離。

通常平行向量法中至少一個向量會涉及速度梯度，對于Sujudi & Haimes法，有幾種平行向量[60]形式。

PEIKERT和ROTH[22]使用平行向量，將Sujudi & Haimes法的特征矢量法等價表示為v||Jv，相比于求解矩陣的特征值和特征向量，矩陣相乘運算的復雜度要小得多；FUCHS等[61]將as||u作為判斷依據(jù)，其中，as=?v·u，實特征值對應的單位特征向量n||u。同理，若將速度梯度?vi(i=1,2,3)作為粒子屬性，可將Levy方法表示為

對于WEINKAUF等[28]的方法，可改寫為

THEISEL和SEIDEL[62]引入特征流場框架，該方法通過將非定常流轉換為高維定常流進行計算。WEINKAUF等[28]引入特征流場，將Sujudi &Haimes法的特征矢量法從流線推廣到跡線。由于在非定常流中需要考慮時間，因此時間維度需加入到計算中。通過對一維和二維空間的推論，用式(3)表示非定常三維流場中的跡線，雅可比矩陣具有e1，e2，e3，0等4個特征值和各自對應的特征向量，即

在旋渦區(qū)域，要求其中2個特征值e1，e2為復數(shù)，對應的2個特征向量所在面即為旋渦平面。由共面性質可得，其中，。即

再使用平行向量場搜索所有a||b的點，最終找到渦核。

ROTH和PEIKERT[29]提出了一種利用二階導數(shù)提取渦核線的方法，即發(fā)現(xiàn)特征向量法實質上是找到加速度a與速度v平行的點，或者曲率為零的點。一階導數(shù)加速度a和二階導數(shù)b被定義為Jv和速度對于時間的二階導數(shù)。

文獻中為渦核線提出了2個新的屬性：旋轉強度和解的質量；將雅可比矩陣復特征值的虛部作為旋轉強度，將渦核與速度矢量之間夾角的余弦作為解的質量，對于實際旋渦而言，這個角度應該小于一定值。這2個新的屬性使得用戶可以對渦施加新的閾值，以消除弱渦。該方法還引入平行向量場，由此還可與脊線等方法結合。

2.2.2 預測校正方法

SCHINDLER等[30,63]提出一種直接在SPH數(shù)據(jù)上提取渦核線的算法，即將其分為種子點生成和預測校正2部分。

種子點生成部分需要考慮整體數(shù)據(jù)集，該部分極大限制了算法的性能。平行向量計算步驟中，該方法未進行常規(guī)SPH方法的插值計算，而是直接使用原始數(shù)據(jù)進行平行向量計算以降低復雜度，經(jīng)驗證發(fā)現(xiàn)誤差可以忽略。

預測校正部分通過計算平行向量c=(c1,c2,c3)=v×w，并求解梯度?c1，?c2，?c3找出叉積?ci×?cj最大的一對i和j，該叉積的方向即為預測方向。若給定一點x，校正點x′=x+s?ci+t?cj，其中，s和t的值由

得出，再計算校正點。

與LEVY等[18]、SUJUDI和HAIMES[27]、FUCHS等[61]、WEINKAUF等[28]4種方法優(yōu)劣勢進行對比。該方法在降低時間復雜度方面做出許多改進，如直接使用原始數(shù)據(jù)進行平行向量計算避免復雜的插值運算；設置容錯濾波和旋渦強度，在雅可比矩陣具有虛特征值并且渦核線和速度矢量角度小于45°時識別旋渦，從而過濾微弱特征等。在輸出效果上，設置最大連續(xù)次數(shù)防止渦核線過度分割(圖4)，并選用三次插值實現(xiàn)渦核線平滑處理。

2.2.3 基于線的方法的問題

該類方法要對數(shù)據(jù)集中滿足條件的各個數(shù)據(jù)點進行求解雅可比矩陣特征值等計算，并連接起來形成渦核線特征，在連接和計算過程中存在特征過度分割、時間復雜度較高等難于解決的問題。預測校正法中提到設置最大連續(xù)次數(shù)以保持連續(xù)性[30]，但仍然需要通過實驗效果來人為控制連續(xù)次數(shù)大小；CABRAL等[64-65]結合LIC方法、λ2準則來提取渦核線，并且用矢量箭頭來表示粒子旋轉方向。

2.3 基于幾何曲線的方法

2.3.1 曲率中心法

該方法由SADARJOEN等[31-32]提出。首先對二維流體進行采樣，形成能覆蓋流體的大量流線。確定每個采樣點的曲率中心(圖5(a))，對于旋渦區(qū)域，曲率中心會趨于一點；對于非旋渦區(qū)域，曲率中心無規(guī)則性(圖5(b))。由曲率中心點網(wǎng)格化形成曲率中心密度(curvature center density, CCD)標量場，再通過設置閾值CCD＞0.8CCDmax來找到旋渦區(qū)域。該方法同樣存在閾值和弱渦問題，并且真實流體中相鄰旋渦的相互作用導致流線形狀不規(guī)則，存在峰值丟失等問題。

圖4 Levy標準(黃色)的渦核線是相連的，F(xiàn)uchs方法的渦核線(紅色)被分割為2個較短部分[63]Fig. 4 The vortex core lines of levy (yellow) are connected, and the vortex core lines (red) of fuchs are divided into two shorter parts[63]

圖5 旋渦區(qū)域與非旋渦區(qū)域[32]Fig. 5 Vortex region and non-vortex region[32]

2.3.2 纏繞角度法

SADARJOEN等[31-32]提出纏繞角概念，在二維流體中選取滿足一定旋轉條件的流線來提取旋渦。設Si為二維流體中的流線，纏繞角為

其中，Si由點Pi，j和線段(Pi，j，Pi，j+1)組成，∠(A，B，C)為線段AB和BC之間的夾角。滿足如下2個條件時存在旋渦：①流線的纏繞角為2π時得到旋渦，不完全等于2π時得到弱渦[66]；②流線的起始點在終止點附近(由閾值判斷)。

該方法大致分為聚類和量化2個步驟。聚類部分將屬于同一旋渦的流線歸為一組，與屬于其他旋渦的流線進行區(qū)分；量化部分計算采樣點和旋渦的屬性，如流線中心、方向、角速度等。與曲率中心法相比，纏繞角度法時間復雜度較高，但在一定程度上可以克服弱渦問題。葉靈偉和馮威[67]提到可通過跟蹤流線上的關鍵點附近區(qū)域的流線來避免全局搜索，降低時間復雜度。江俊等[68]選取已提取出的局部對稱中心處3×3流線方向場模塊I替換式(19)中的Pi，j來計算纏繞角度值，實驗證明該方法具有較好的穩(wěn)定性。

2.3.3 基于幾何曲線的方法的問題

經(jīng)典方法同樣受限于閾值，并且時間復雜度相對較高，不能直接擴展到三維空間。對于三維空間問題，REINDERS等[69]改進纏繞角度法，在每個時間步長下的每個水平面中提取二維旋渦并逐一進行連接構造三維渦管。使用特征跟蹤方法[70-71]進行實時跟蹤，從而在連續(xù)的時間步長內提取實時三維旋渦。JIANG等[72]在三維空間中考慮渦核曲率，沿每個點處的流線計算探測向量來探測渦核在特定點的方向矢量，在測量過程中提供了必要的曲率信息，該方法適應于非定常流的渦核提取。陳麗娜和楊冠杰[73]提出角度函數(shù)法，用于快速提取渦核區(qū)域，通過設置β表示樣本點處矢量方向按一定規(guī)則變化的程度，β極小點即為渦核點。該方法只能判別大致的渦核區(qū)域，但可以很好地解決時間復雜度問題。

3 基于機器學習的流體旋渦提取方法

3.1 基于BP神經(jīng)網(wǎng)絡的方法

XU等[33]設計了一種三層BP神經(jīng)網(wǎng)絡進行流場的智能化特征提取。該方法利用臨界點理論，將流場拓撲構成分為臨界點和連接臨界點的積分曲線，因此識別速度與臨界點有密切關系。在臨界點提取部分，采用Marching Cubes算法，認為當某個網(wǎng)格單元的4個頂點的8個矢量分量同時存在符號相反時，網(wǎng)格內存在臨界點。計算三層結構中隱含層和輸出層輸出值，再計算輸入樣本與各類的誤差值，其中，dki為第k個標準類的第i個理想輸出，yi為實際輸出。取其中的最小值Ev，若小于給定閾值，則認為是第v個流場特征。

在時間復雜度方面，訓練神經(jīng)網(wǎng)絡作為預處理過程可以忽略不計；特征提取部分的復雜度為O(PLM)，其中P為候選樣本區(qū)域，L為隱含層節(jié)點數(shù)，M為輸入層節(jié)點數(shù)。

根據(jù)文獻中的實驗表明，該方法顯著減少了數(shù)據(jù)計算量，并且可以較好地突出流場特征。

3.2 基于CNN的方法

文獻[34]提出了一種基于CNN的特征可視化方法，主要研究了二維流場的特征可視化。該方法可以提取流場中用戶所需的典型特征和非典型特征并加以分類。與BP神經(jīng)網(wǎng)絡相比，CNN具有較好的識別性，可以有效處理弱渦、變形、平移等情況。

該方法主要分為數(shù)據(jù)預處理、特征區(qū)預選和區(qū)域劃分3個部分。數(shù)據(jù)預處理部分主要記錄橫縱坐標和對應的水平垂直速度分量并對矢量大小進行顏色映射，再通過

對矢量數(shù)據(jù)進行單位化。特征區(qū)預選部分在水平和垂直方向平移窗口矩陣來分割矢量矩陣，將具有相同大小的樣本矩陣作為訓練后的CNN模型的輸入。并判斷其特征類別，可將所有具備相同特征的區(qū)域歸為一類。根據(jù)矢量拓撲理論可知順時針渦、逆時針渦和鞍形特征都包含臨界點，因此候選區(qū)預選部分主要搜索矢量矩陣的臨界點。同2.4.1節(jié)中搜索臨界點的方法類似，網(wǎng)格單元4個頂點的矢量分量在水平和垂直方向至少正負轉換一次時存在臨界點，以該網(wǎng)格為中心的樣本矩陣作為后續(xù)識別的候選區(qū)域。通過掃描整個數(shù)據(jù)集的數(shù)據(jù)矩陣來獲得臨界點的大致位置，由此得到候選點列表。區(qū)域劃分部分由每個候選點得到待識別矩陣，再經(jīng)訓練后的CNN模型進行識別，并將結果存入特征列表用于后續(xù)可視化。

實驗表明該方法精度高達96.67%，效果優(yōu)于BP方法，但其需要用戶進行多項干預才能使CNN模型達到較好的識別效果，具有主觀性。

DENG等[35]提出了一種具有更高的精度和有效召回率的客觀旋渦提取方法。該方法包括預處理和渦網(wǎng)2個部分。

預處理部分將物理平面上的非均勻網(wǎng)格轉換為均勻網(wǎng)格從而得到其位置和速度信息，使用IVD[74]全局方法標記流場中的點，用0代表非旋渦點，1代表旋渦點，在歸一化速度場附近的局部區(qū)域進行采樣，將這些局部區(qū)域和全局標簽作為第2部分的數(shù)據(jù)輸入。通過該步驟，可以近似達到全局方法的效果。渦網(wǎng)部分訓練CNN模型，直接將測試區(qū)域進行輸入，通過渦網(wǎng)得到識別結果，從而將輸入數(shù)據(jù)分為旋渦點和非旋渦點。

該方法進行了多組對比實驗，如在不同尺寸和形狀的流場中驗證通用性；與局部方法和全局方法進行對比驗證識別精度和速度；與傳統(tǒng)機器學習算法相比，可以更好地反映旋渦的分離運動。

3.3 基于機器學習方法的問題

傳統(tǒng)的特征提取方法大多數(shù)只能提取單一特征，并且需要復雜的數(shù)學運算，因此基于機器學習的旋渦提取方法的優(yōu)勢逐漸展現(xiàn)出來。模型可在預處理部分進行訓練和學習，實際提取時可忽略訓練的時間；訓練后的模型可以應用到類似的其他流場中；相比于傳統(tǒng)的局部和全局提取方法，有更高的精度和速度。

但是，其局限性也需在后續(xù)的研究中加以改進。旋渦識別階段很難得到用來提取旋渦的標記數(shù)據(jù)，由于缺乏旋渦的嚴格數(shù)學定義，使得訓練網(wǎng)絡具有一定的缺陷；在提取精度方面，KIM和GüNTHER[75]在訓練數(shù)據(jù)時添加噪聲并且重采樣；模型的通用性較弱，流場的尺寸不定，而部分模型(如CNN)卻需要固定大小的輸入；數(shù)據(jù)計算量較大，相比于局部方法，在有高精度的同時仍需要更低的時間復雜度。LIU等[76]提出了一種基于CNN的沖擊波提取方法，訓練包含多個卷積層的提取網(wǎng)絡，該方法具有較好的提取效率和結果。WANG等[77]提出全卷積分割網(wǎng)絡旋渦識別方法，在保證精度的同時減少參數(shù)數(shù)量并降低計算復雜度；大多數(shù)方法的訓練集缺乏權威性，DUO等[78]對少量由海洋學專家確認的精確旋渦樣本進行增強以生成訓練集，構建自動識別定位網(wǎng)絡OEDNet。

4 參考系不變性

旋渦提取方法可以根據(jù)參考坐標系運動的類型分為伽利略不變性、旋轉不變性和拉格朗日不變性。不變性是指實驗結果(識別的旋渦結構)不隨參考系的變化而出現(xiàn)假陽性和假陰性。

大多數(shù)旋渦提取方法中，實驗結果依賴于觀察者的運動，如圖6所示。飛行員以不同的路徑在矢量場中運動時觀察到的場景不同。在標量計算中，實現(xiàn)參考系不變性較為容易，例如兩輛同向行駛的列車，以其中一輛列車為坐標系，另一輛的運動距離即為該列車原本的運動距離減去參考列車的運動距離。對于矢量而言，不同參考系的選擇會影響旋渦測量的結果，如圖7所示。在參考系向右移動時，抵消部分速度，才能顯示旋渦。

圖6 使用線積分卷積和跡線(黑色較粗)示出3種不同的參考系運動，分別為靜止、線性平移和沿正弦曲線移動[79]Fig. 6 Three different reference frame motions are shown using line integral convolution and pathlines (thicker black),which are static, linear translation, and moving along a sinusoidal curve[79]

圖7 不同參考坐標系下的對比[80]Fig. 7 Comparison under different reference frames[80]

旋渦提取方法可根據(jù)參考系的不同分別選取伽利略不變性、拉格朗日性不變性和旋轉不變性，如圖8所示。其中，線代表粒子運動軌跡，即跡線；箭頭為粒子速度分量方向。伽利略不變性適用于提取等速平移的參考坐標系下的旋渦結構(圖8(a))；旋轉不變性適于提取繞固定旋轉軸旋轉的旋渦結構(圖8(b))；拉格朗日不變性可以保證無論參考系做任何旋轉或平移變換均能得到一致的旋渦結構(圖8(c))。

圖8 參考系變換下的3種不變量[36]Fig. 8 The three invariants under the transformation of reference frame[36]

4.1 伽利略不變性

參考系以

進行變化時，其中，c0為常數(shù)點；c為常向量，若實驗結果不隨參考系變化，則具有伽利略不變性。

有部分物理量具有伽利略不變性，如加速度、雅可比矩陣等。任何只計算加速度、雅可比矩陣或其導數(shù)的方法均具有伽利略不變性。如：Q準則、λ2準則等。伽利略不變性在近年的旋渦提取中應用廣泛，JEONG和HUSSAIN[26]強調了伽利略不變性的重要性。SAHNER等[3]提出了利用平行向量提取伽利略不變渦區(qū)量的脊線或谷線方法。該方法可以應用于三維空間，通過將脊線或谷線轉換為曲面進行提取。

4.2 旋轉不變性

旋轉不變性參考系以

進行變化，其中，Q(t)為等速旋轉矩陣；x0為旋轉中心。旋轉不變性能保證在參考系等速旋轉時準確提取旋渦。該方法的缺點是需要計算旋轉中心點x0和三維旋轉軸n，如圖9所示。在旋轉不變性方法中，由Q(t)旋轉矩陣將n變化為z軸，旋轉不變性雅可比矩陣可表示為

其中，

適用于伽利略不變性參考系中的旋渦提取方法，如Q準則和λ2準則，可以通過旋轉不變性雅可比矩陣Jr轉化為旋轉不變性形式。有

則

圖9 二維和三維旋轉不變性[36]Fig. 9 Two-and three-dimensional rotation invariance[36]

GüNTHER等[36]引入非等距函數(shù)gp(x,t)，逆函數(shù)hp(x,t)以及轉化后的向量場wp，實現(xiàn)笛卡爾坐標系和旋轉不變性坐標系相互轉換。但是，非等距函數(shù)會導致雅可比矩陣離散化，由此提出僅在笛卡爾坐標系下計算旋轉不變雅可比矩陣的方法以解決離散化問題。文獻中使用多組數(shù)據(jù)集對比伽利略不變性方法與旋轉不變性方法。其中一組實驗選用由WIEBEL等[81]在2011年提出的無散度珠流，在一段時間內被旋渦特征可視化領域公認為是評定旋渦提取方法優(yōu)劣的一個基準。在該對比實驗中，發(fā)現(xiàn)旋轉不變性方法能正確跟蹤渦核；伽利略不變性方法的渦核線逐漸偏離跡線，如圖10所示。

圖10 伽利略和旋轉不變性渦核線對比圖[36]Fig. 1 0 Galilean and rotation invariant vortices line comparison diagram[36]

4.3 拉格朗日不變性(客觀性)

拉格朗日不變性參考系以x′=Q(t)x+c(t)進行變化，其中，Q(t)為實正交矩陣；c為隨時間變化的平移向量。拉格朗日不變性方法能夠識別進行任何平滑旋轉和平移運動的旋渦。向量場中，v，vt，J，a，Ω等物理量均不具有拉格朗日不變性，散度?v和應變率張量S是向量場中為數(shù)不多的具有拉格朗日不變性的物理量[82]。

自拉格朗日不變性的概念出現(xiàn)后，其實驗結果基本不受參考系變化影響的特質引起了可視化領域廣泛關注，近年來提出了許多經(jīng)典的拉格朗日不變性形式。

TABOR和KLAPPER[37]將Q準則應用在拉格朗日不變性參考系中，通過

將矩陣Q轉化為與應變率張量S相關的形式，其中，Ωs是由應變率張量S的單位特征向量的時間導數(shù)組成的矩陣。

DROUOT和LUCIUS[83]定義了旋轉率張量的相對值，其中，ΩS為S的特征向量的旋轉率張量。MARTINS等[84]運用該相對值并推得拉格朗日不變參考系下的Q準則和λ2準則類似，該文獻中還對Δ準則[85]和CHAKRABORTY的局部渦識別方法[86]進行拉格朗日不變性處理。

除了對經(jīng)典方法的改進外，有許多直接應用在拉格朗日不變性的方法。SADLO和PEIKERT[87]將直接Lyapunov指數(shù)(directional Lyapunov exponent,DLE)用于提取旋渦。GüNTHER等[79]通過為域中的每個點計算局部最優(yōu)參考系得到速度場、加速度場和雅可比矩陣的拉格朗日不變性形式，并且用λ2準則和Sujudi & Haimes法來驗證提取的旋渦特征的準確性。GüNTHER和THEISEL[88]將拉格朗日不變性參考系推廣到慣性粒子的渦核提取中，能夠提取在任何平滑旋轉或平移路徑上移動的慣性粒子的旋渦，但是6D平行矢量場的引入使得計算更加復雜，需要研究更快的提取方法。HADWIGER等[38]首先提出了用于非定常流場特征提取的全新客觀參考系以計算時變的觀察者速度場和近似Killing場，讓所有的觀察者以各自不同的視角共同觀測流場。該方法結合全局和局部參考系，在全局范圍內提取的同時也可以進行局部調整，相比GüNTHER等[79]的方法，其鄰域大小固定，普適性較弱。

5 總結與挑戰(zhàn)

經(jīng)過幾十年的研究，旋渦提取技術已經(jīng)有了很大的發(fā)展，但是仍存在以下問題：

(1) 對于旋渦本身，迫切需要一個正式的定義。LIU等[51]提出的旋渦矢量(Rortex)是一個很好的開端。其可以準確描述局部流體旋轉情況，清晰顯示旋渦結構。

(2) 對于基于點的方法，閾值在各場景的適應性是難點，大多數(shù)方法仍需通過實驗效果對比來人工設置最合適的閾值。對于該問題，本文對Ω方法[19]進行了論述。后期文獻[58,89]給定了ε合適的參數(shù)值，能準確識別各場景下的旋渦。

(3) 對于基于線的方法，屬于同一旋渦的渦核線經(jīng)常被分割成多個部分，導致實驗效果下降。預測校正法[30]在一定程度上增強了渦核線的連續(xù)性，但仍需減少人為因素。

(4) 對于基于幾何曲線的方法，采樣點的選取對實驗結果有極大的影響。三維空間實時提取旋渦方法時間復雜度較大，需要進一步地研究。

(5) 對于基于機器學習的方法，該類方法意在通過大量樣本的訓練從而自主判斷是否存在特征，因此訓練集必須具備權威性；仍需加強模型的通用性，以適應不同尺寸大小的流場。

(6) 對于參考系不變性，本文先對參考系不變性給出解釋，對每種不變性給出經(jīng)典方法并進行實驗效果對比。目前，最理想的方法是拉格朗日不變性參考坐標系，該方法有助于研究更為復雜的流體場景。

對于旋渦提取領域未來的發(fā)展，研究者應該更加致力于正式定義旋渦、提高旋渦提取精度、降低時間復雜度等方面，客觀方法、全局及局部結合提取方法和基于機器學習的方法能在一定程度上解決各類問題，是未來研究的一個方向。