基于OPTICS聚類算法的流場結構特征分析方法

王 銳，辛大波，歐進萍

(1. 哈爾濱工業大學 土木工程學院，哈爾濱 150090；2. 東北林業大學 土木工程學院，哈爾濱 150040)

0 引言

長大跨度橋梁、遠距離大跨度的輸電塔線、高聳通訊塔桅等基礎設施，分別構成了交通網絡、能源網絡和通訊網絡中的關鍵節點，是我國乃至世界各國經濟社會發展所需要的基礎設施工程。這些工程結構往往具備長、柔的特征，對風荷載敏感，易產生顯著的風致靜力和振動響應[1-3]。因此，對長、柔結構物的風工程問題進行研究，對于提高結構抗風能力，保障基礎設施網絡的安全性、可靠性和魯棒性具有十分重要的意義。

流場結構特征分析是風工程研究中的重要內容，是揭示結構物周圍的流場結構，歸納流動規律，分析結構的風荷載成因與控制方法的必要途徑。一般結構物附近的流場可以由基于計算流體動力學（CFD）的數值模擬手段獲得[4-7]。近年來，隨著如粒子圖像測速（PIV）這樣的流場可視化技術的發展，利用風洞試驗也可以獲取結構周圍流場[8-9]。

流場數據一般包含流場區域中各點對應的速度、渦量、湍動能等數據，是一種典型的高維數據。對流場數據的分析要依賴于各種各樣的降階模型[10]。本征正交分解（POD）是一種在流場動力分析中常用的降階方法，其假設每一個流場樣本可以視為不同流場模態的線性組合，并通過對特征值問題的求解，獲取各流場模態的統計顯著程度[11-12]。動模態分解（DMD）則假設每一個流場樣本可由上一時刻的流場樣本通過線性映射矩陣獲得。相對于POD，DMD的系統矩陣中囊括了更多的時序信息，從而能夠獲取更準確的模態頻率[13-14]。無論是POD還是DMD，都需要依賴于一些動力學假設，但這些假設相對于湍流這樣復雜的動力學問題過于簡單。因此，有必要研究不依賴于動力學假設的流場結構特征分析方法。

Kaiser等提出將聚類分析方法引入到流場分析技術中[10]。聚類分析是一種由數據驅動的機器學習方法，用于將樣本數據分組為不同類，每一類中的樣本具有一定程度的相似性[10]。聚類分析完全依托于流場數據，不依賴于動力學假設，因此在流場分析的應用上更有前途。目前有一些聚類分析在PIV試驗[15]和數值模擬[10,16]結果分析的應用，在對流場結構特征的分析中取得了較好的結果。值得注意的是，上述研究中普遍采用的是聚類分析中的k-means算法。k-means算法是一種基于平均值的聚類算法，即：在所有數據中隨機選取k個作為初始聚類中心，通過對樣本和聚類中心歐氏距離的計算和比較，將樣本歸于距離最近的聚類中心；分類完成后以各類樣本平均值作為新的聚類中心替換初始聚類中心，重復上述計算直至本次迭代得到的新聚類中心與本次迭代的初始聚類中心重合[17]。這樣的算法就意味著，k值的選取不唯一，計算得到的結果不唯一，且最終得到的分類在樣本中均勻分布。對于識別流場的結構特征，k-means算法稍顯粗糙。

考慮到聚類算法在流場結構特征分析中的優越性和目前常用的k-means方法的一些缺陷，本文將關注于尋找一種算法以改善基于k-means算法的流場聚類分析方法。本文選取了一種基于密度的聚類算法OPTICS （Ordering Points To Identify the Clustering Structure），并通過引入相關距離的概念替換歐氏距離，形成了基于相關距離的OPTICS算法。利用基于大渦模擬（LES）的計算流體動力學（CFD）數值模擬，獲取了經典圓柱繞流尾流中的順流向渦渦量場樣本。以識別順流向渦脫落狀態和其A模式分布的展向間距為目標，檢驗基于相關距離的OPTICS算法的有效性，并分別與原始的OPTICS算法和k-means算法對比。結果表明，基于相關距離的OPTICS算法是一種有效的用于流場結構特征分析的聚類分析算法。

1 基于OPTICS算法的流場聚類分析

1.1 OPTICS聚類算法的基本概念

OPTICS是一種基于密度的聚類分析算法，即以數據點的空間密度作為判別標準，將每個高密度區域中的數據點作為一類。相對于k-means算法需要由人工設置初始聚類數目且聚類結果不唯一，OPTICS算法能夠基于初始參數識別樣本中的高密度數據點集，從而相對客觀地找到樣本中的高密度區域，給出符合樣本分布規律的聚類。對于確定的樣本，其數據分布可以唯一確定，因此基于密度判別的OPTICS算法可以獲得相對穩定的聚類中心。針對于流場結構特征分析，當選取合理的相似度指標，OPTICS算法識別出的高密度數據集即意味著該數據集中的流場樣本高度相似，也即具備同一特征。那么通過OPTICS算法的聚類，就可以相對客觀地判別流場樣本所包含的結構特征。

假設已經獲得了n組流場數據X= [x1,x2,···,xi,··· ,xn]，其中xi代表第i個流場數據。每個流場數據為一個m維向量，代表了流場中的m個測量值，也代表了m維歐氏空間中的一個點。OPTICS算法的關鍵概念，是對每一個xi建立一個半徑為ε的鄰域Nε(xi)，且應有不少于MinPts個數據點在該領域中。據此我們可以引入以下概念[18-19]：

1）密度直達：如果X中的兩個流場xi、xj(i≠j)滿足：xi∈Nε(xj)且card(Nε(xj)) ≥MinPts，則稱xi可由xj密度直達。其中，card(Nε(xj))表示集合Nε(xj)包含的元素數，card(Nε(xj)) ≥MinPts的條件被稱作核心對象條件，xj可以被稱為核心對象，其他流場數據只能由核心對象密度直達。

2）密度可達：如果存在一系列X中的流場數據xr、xr+1、···、xr+s，當i∈[r,r+s? 1]時xi可由xi+1密度直達，則稱xr可由xr+s密度可達。

3）密度相連：如果X中的流場數據xi和xj(i≠j)均可由X中的另一流場數據xk(k≠i且k≠j)密度可達，則稱xi和xj密度相連。

4）類：類是聚類分析得到的結果，一類流場數據Y是X的一個非空子集，且滿足以下兩點：①對任意xi、xj(i≠j)∈X，如果xj∈Y且xi可由xj密度可達，則xi∈Y；②對任意xi、xj(i≠j)∈X，xi和xj密度相連。如果有些X中的流場數據沒有被分在任何一類中，則稱之為噪聲。

5）核心距離：從xi到其鄰域Nε(xi)中某一數據點，使xi成為核心對象的最小距離，可以數學表達為：

式中：dc(xi)是xi的核心距離；為在xi的鄰域Nε(xi)中，距離xi第MinPts近的流場數據；d(xi,表示xi和之間的距離。

6）可達距離：對xi、xj(i≠j)∈X，xi關于xj的可達距離可以定義為：

其含義是，使xi可由xj密度直達和xj為一個核心對象同時成立的最小距離。

1.2 OPTICS聚類算法用于流場聚類的步驟

執行OPTICS聚類算法，需要首先輸入流場數據序列X，其次是鄰域半徑ε和鄰域集合最小元素數MinPts，然后執行以下步驟[18]：

1）在X中隨機選擇一組數據xi作為當前對象（執行中往往選擇第一組數據），將其標記為已處理，并描繪在決策圖順序軸的第一個位置，其對應的可達距離是無意義的，一般在程序中賦予一個較大值。

2）計算X中所有其他流場數據關于當前對象的可達距離。

3）找到關于當前對象可達距離最小的一組數據，用該組數據替換當前對象并標記為已處理，將該組數據按處理順序和可達距離描繪在決策圖對應位置。

4）依次計算未處理數據關于當前對象的可達距離，如果這些可達距離中存在小于步驟3）中計算出的可達距離，則將其對應的數據替換為當前對象；否則保持原有當前對象不變。

5）重復步驟3）、4）直至X中所有數據處理完畢。

OPTICS算法得到的結果用決策圖表示。如圖1所示，其橫軸表示流場數據的處理順序，其縱軸表示該流場數據關于處理時當前對象的可達距離。通過與已經設置好的鄰域半徑ε對比，若可達距離dr<ε，則該可達距離有意義，其對應的流場數據被聚為一類。圖1的決策圖表示該組數據被聚為三類。

圖1 OPTICS算法的決策圖示意Fig. 1 An example of the reachability-plot of OPTICS

2 基于相關距離的OPTICS聚類算法

2.1 基于歐氏距離OPTICS聚類算法的缺陷

OPTICS聚類算法依賴于鄰域半徑、核心距離和可達距離等距離概念，基于這些距離來判別數據的分布密度。傳統的OPTICS算法往往采用歐氏距離[18-19]，即將流場數據X中的xi和xj視為m維歐式空間中的兩個點，將xi和xj中的元素作為其在m維歐式空間中各個維度上的坐標，計算兩點間的幾何間距。采用歐氏距離定義清晰，計算方便，但卻忽略了流場數據的空間分布特性，割裂了一個流場數據內各個元素間的聯系，這對于識別流場的模式和結構特征是不利的。

2.2 常用相似度指標在流場結構特征識別中的適用性

常見的衡量高維數組間距離的相似度指標包括歐氏距離、切比雪夫距離、閔可夫斯基距離、馬氏距離、相關系數、夾角余弦等，其定義如表1所示。根據這些指標的定義，可以將其簡單分為兩類，其一是以空間坐標間距（ai -bi）為核心的距離概念，包括歐氏距離、曼哈頓距離、切比雪夫距離、閔可夫斯基距離和馬氏距離；其二則是以統計相關性為核心的相關系數概念，包括夾角余弦、Pearson相關系數和Spearman相關系數。

表1 衡量高維數組間相似程度的常見相似度指標Table 1 Common indexes for evaluating the similarity among high-dimensional arrays

流場樣本往往是一套以空間坐標標記的物理量，該物理量可以是流速、壓力、渦量、湍流度、或湍動能等。在對流場結構特征進行識別和分析時，需要將流場樣本的空間坐標按一定規律映射于一個一維數組的元素角標，而該一維數組中的元素為各角標對應的物理量，從而將原始的流場樣本組織為高維數組。流場的結構特征主要體現在物理量的空間分布，其含義既包括物理量本身的大小，也包括特定大小的物理量出現的空間位置。例如根據流速和渦量的空間分布，可以識別流場中旋渦的形狀、位置與強度；根據鈍體邊界層中壓力的空間分布，可以識別到鈍體繞流的分離點。因此，識別流場的結構特征，對流場樣本進行相似度的度量，既需要關注物理量本身的大小，即流場樣本數組中的元素大小；也需要關注物理量的空間坐標，即流場樣本數組中的元素角標。

為了說明基于統計相關性的相似度指標在流場分析中的優勢，這里舉例如下。假設存在三組簡單的流 場 樣 本：a= [1,0,1]T，b= [0,2,0]T，c= [3,0,3]T。a和c的元素大小不同，但其中較大值的元素出現在相同的位置，即元素的分布規律相似，這代表著這兩個樣本具有相似的流場結構；a和b中較大值的元素出現在不同位置，這代表著這兩個樣本的流場結構顯然不同。按上述指標計算a、b之間和a、c之間的距離，結果如表2所示。顯然，基于空間坐標間距（ai-bi）的距離指標普遍沒能識別出a和c之間具有更加相似的物理量分布，而基于統計相關性的夾角余弦、Pearson相關系數和Spearman相關系數均體現出了這一點。

表2 基于不同相似度指標計算的a、b、c之間的距離Table 2 Distances among a,b and c calculated by different similarity indexes

當對流場樣本這樣的高維數組進行相似度度量時，以空間坐標間距（ai-bi）為核心的距離指標，更加注重元素大小的衡量，而忽略各個維度間的聯系；以統計相關性為核心的相關系數指標則偏重于關注對元素分布的衡量，這一特性對于識別流場物理量的空間分布更加有利。具體分析各個指標，其中歐氏距離、曼哈頓距離和切比雪夫距離分別是閔可夫斯基距離在p= 2、p= 1和p趨于正無窮時的特殊情況，這一類的指標顯然割裂了不同維度的元素間的聯系。馬氏距離的定義式中包含了樣本的協方差矩陣的逆矩陣，這就要求樣本的協方差矩陣滿秩，對于流場樣本這種維度往往在104以上數量級而樣本數相對少的樣本矩陣，其協方差矩陣通常不能滿足滿秩的條件。如果采用POD方法對流場樣本降維處理，又會引入線性系統的假設。所以馬氏距離在流場結構特征識別中只能有限應用于樣本量不少于樣本維度的情況。夾角余弦和Pearson相關系數形式接近，當樣本A和B的均值均為0時，Pearson相關系數就可以寫作夾角余弦的形式，也即Pearson相關系數是排除了樣本元素均值影響的夾角余弦。當樣本元素的均值存在較大波動，基于夾角余弦的相似度衡量可能受到均值的影響。Spearman相關系數基于樣本A和B的秩次獲得，也就是忽略樣本中元素的具體數值，而利用其在所有元素中的大小排序進行計算。Spearman相關系數完全忽略了樣本元素的具體數值，對于流場這種包含同一物理量的樣本集并不適用。綜上所述，Pearson相關系數相對更適合于作為衡量樣本相似度指標，應用于基于OPTICS聚類算法的流場結構特征識別。

2.3 基于相關距離的OPTICS聚類算法

基于對流場樣本數據特征和各類距離指標的分析，本文引入基于Pearson相關系數的相關距離的概念，并用相關距離替換了傳統OPTICS算法采用的歐氏距離，以使OPTICS聚類算法更加適用于識別流場結構特征。首先，Pearson相關系數的概念可以按下式定義：

式中：A、B表示兩組相同維度的向量；ρAB為向量A、B的相關系數；Cov(A,B)為A,B的協方差；E(A)和E(B)表示A和B的期望；σA和σB表示A和B的標準差。根據Pearson相關系數可以定義向量A、B的相關距離如下式：

式中：dAB表示向量A、B的相關距離。相關距離越小，表明向量A和B越相似；反之則表明向量A和B差異越大。

3 典型圓柱繞流場結構分析實例

3.1 流場結構特征識別對象

圓柱繞流是一種經典的鈍體繞流現象，實際工程結構中具有大量的圓柱形彈性結構，例如大跨橋梁的斜拉索、吊索、長大跨度輸電線、海洋管道等。這些結構極易受到來流作用而發生流致響應，導致結構的疲勞損傷甚至破壞。對圓柱形結構的三維繞流場進行流場結構分析，對于研究其流致響應機制，開發流動控制方法具有重要的作用。因此，本文以識別圓柱尾流中順流向渦的A模式[20]為目標，通過CFD數值計算手段獲取流場分析對象，用來說明基于相關距離的OPTICS聚類算法的優越性。數值模擬對象為一個沉浸于均勻流中且與來流方向正交的靜止圓柱體，長l= 0.2 m，直徑d= 0.01 m，長細比20∶1，其軸線與z軸重合。計算域為一個高20d，半徑為35d，軸線與圓柱模型軸線重合的圓柱體，如圖2所示。計算域各邊界條件如圖2所示，來流風是U= 0.35 m/s的均勻流，從入口沿x軸方向進入，從出口自由流出。計算的雷諾數Re約為240。

圖2 計算域幾何特征與邊界條件Fig. 2 Geometry features and boundary conditions of the calculation domain

計算域采用結構化網格劃分，共計包含2187416個六面體網格。計算域被沿圓柱軸線等分為67層，每層網格如圖3所示。其中靠近圓柱壁面的網格約為0.1 mm×0.1 mm的正方形，沿周向等分，沿徑向按擴散率1.05向外側擴大。圓柱模型表面的壁面y+如圖4所示，模型表面各點的y+均小于0.5，說明第一層網格足夠精細以滿足計算需求。計算利用商業軟件Fluent執行瞬態模擬，采用三維大渦模擬（LES），亞尺度模型選用WALE。動量按二階迎風格式離散，壓力和瞬態方程采用二階格式離散。計算時間步長為5×10?5s。

圖3 網格劃分示意Fig. 3 A schematic of the mesh

圖4 圓柱壁面y+Fig. 4 y+ on the surface of the circular cylinder

為了驗證數值計算結果的準確性，本文對圓柱模型氣動力系數的統計值和斯托拉哈數（St）與相關參考文獻提供的其他試驗結果的擬合值進行了對比，對比結果如表3所示。數值計算獲得的圓柱阻力系數時均值與試驗擬合結果的偏差均小于5%，數值計算的St與試驗擬合結果的偏差約為5%，數值計算結果可信。

表3 氣動力系數驗證Table 3 The comparison of the numerical and experimental aerodynamic forces

流場分析的識別目標是在x= 2.5d切面上，識別圓柱尾流中順流向渦的A模式旋渦脫落的兩種狀態與樣本中包含的兩種展向間距（如圖5所示）。順流向渦的A模式是指在雷諾數超過190左右時，圓柱尾流中除卡門渦之外，出現了旋轉軸沿來流方向的旋渦。其出現方式為一系列沿展向按間距約4d排列的方向相反的旋渦對，與卡門渦呈相同頻率脫落，并且每半個卡門渦脫周期發生一次旋轉方向的變化。由于數值計算存在的數值誤差和尾流中順流向渦本身的不穩定性，相當于流場本身存在不屬于標準模式的噪聲。有效的流場分析方法應當能識別到x= 2.5d切面上，順流向渦的展向分布間距和旋渦脫落的兩個狀態，并排除噪聲流場的干擾。流場樣本共計按時間步長0.005 s采樣2000份，從流動的第4 s采樣至第14 s。樣本數據為x= 2.5d切面上的x方向渦量場。

圖5 圓柱尾流中順流向渦的旋渦脫落與A模式示意Fig. 5 The vortex shedding and the A-mode of the streamwise vortices in a circular cylinder wake flow

3.2 基于k-means聚類算法的流場結構特征識別

k-means聚類算法中，聚類數目N作為初始參數輸入，需要利用總方差的概念對初始聚類數目N進行判別。對于每一類，可以獲得各樣本圍繞該聚類中心的統計方差，各個類的方差之和即為總方差。理論上當N= 1時，總方差達到最大；當N與樣本總數n相同時，總方差為0，最優的聚類數目N應當在較小的總方差和較少的聚類數目中取得平衡。一般地，最佳聚類數目可以經由總方差的“肘形判據”獲取[10]。本例中嘗試將初始聚類數目N設置在2～12范圍內，獲得總方差如圖6所示。總方差隨著聚類數目N的增長單調下降，斜率逐漸減小；但是曲線上并不能明顯找到“肘部”，“肘形判據”并不能明確給出最合理的聚類數目。

圖6 k-means算法的方差與初始聚類設置的關系Fig. 6 The number of initial clusters as a function of the variances of the k-means algorithm

圖7展示了聚類數目分別為2、4、6和8時的流場分析結果，每一類流場的分析結果為各類流場的聚類中心。當N= 2時，兩類流場分別對應順流向渦脫落的兩個狀態，可以認為捕捉到了順流向渦脫落；但是A模式應當具有的約4d展向間距的順流向旋渦對在這兩個結果中體現不明顯，兩種展向間距沒能得到區分。當N> 2時，順流向渦的A模式在部分類別中得到體現，但由于分類增多，出現了兩個問題：其一，順流向渦的脫落狀態不能被這些類明確；其二，當N=8時，出現了類似的多個聚類，即過度分類。

圖7 k-means算法的流場結構特征識別結果Fig. 7 The identified flow structures by the k-means algorithm

總的來說，k-means算法難以在識別順流向渦脫落和順流向渦展向間距之間取得平衡。另外k-means算法的輸出結果嚴重依賴輸入的初始聚類數目N，且難以找到相對最優的N。k-means算法在計算中隨機選取初始聚類中心，所以其聚類結果存在隨機性，同一組樣本數據多次聚類的結果并不一致。

3.3 基于原始OPTICS算法的流場結構特征識別

OPTICS算法不需要像k-means算法一樣設置初始聚類數目，其初始參數是ε鄰域Nε(xi)中最小點數MinPts和鄰域半徑ε，而最終的聚類結果由決策圖獲得。基于OPTICS算法的流場分析結果如圖8所示，共計采用60～200間的六組MinPts，Ni表示第i類中的樣本數量，每個流場的聚類結果由各類中流場樣本的平均場表示。圖8（a）顯示，當MinPts= 60，令ε= 1010，可以在決策圖中找到兩個聚類，分別對應順流向渦脫落的兩個狀態，且順流向渦的A模式比較明顯。但是僅識別到順流向渦對的4d展向間距，5d的展向間距沒能在聚類結果中體現；且由于MinPts取值較小，造成聚類中樣本數量偏小，大量的樣本被識別為了噪聲。圖8（b～f）給出了MinPts≥80的5種情況，其決策圖中均只能找到一個聚類，這一類的聚類中心既不能體現順流向渦的脫落，也不能體現A模式渦對的展向分布。也就是說，MinPts≥80時，OPTICS算法不能夠識別出順流向渦渦量場的狀態。

圖8 原始OPTICS算法的決策圖與聚類中心Fig. 8 The reachability-plots and cluster centers of the original OPTICS algorithm

3.4 基于閔可夫斯基距離的OPTICS算法的流場結構特征識別

閔可夫斯基距離是一種典型的以空間坐標間距（ai-bi）為核心的距離指標，原始OPTICS算法采用的歐氏距離即對應于p= 2時的閔可夫斯基距離。基于閔可夫斯基距離的OPTICS算法的流場結構特征識別結果與原始的基于歐氏距離的OPTICS算法類似，圖9列出了p= 5時的流場分析的決策圖與聚類中心，可見MinPts取值在60～120之間時，選取適宜的鄰域半徑ε，都可以得到分別對應于順流向渦脫落的兩個狀態的兩個聚類中心。相對于p= 2時的聚類結果，能夠識別旋渦脫落的MinPts取值范圍大幅增大，但仍然沒能識別出5d的順流向渦展向間距。

圖9 基于閔可夫斯基距離(p = 5) OPTICS算法的決策圖與聚類中心Fig. 9 The reachability-plots and cluster centers of the OPTICS algorithm based on the Minkowski distance (p = 5)

參數p取值對閔可夫斯基距離的定義十分重要，對基于閔可夫斯基距離的OPTICS算法的流場結構識別結果也有影響。表4列出了參數p的取值對聚類中心數目的影響。當利用基于閔可夫斯基距離的OPTICS算法對順流向渦結構進行識別時，若識別出兩個聚類中心，其結果往往對應于順流向渦脫落的兩個狀態；若識別出一個聚類中心，則意味著沒能區分出任何結構特征。隨著參數p取值的增加，能夠識別旋渦脫落特征的MinPts取值范圍逐漸增大。但當采用切比雪夫距離，當MinPts= 60時出現了過度分類，在識別出順流向渦脫落特征的同時，將順流向渦脫落的其中一個狀態識別成兩種狀態；而MinPts≥80時無法區分出任何流場結構特征。需要指出的是，當聚類中心個數相同時，基于閔可夫斯基距離的OPTICS算法獲取的決策圖圖形類似，聚類中樣本數量普遍偏少，大量的樣本被識別成了噪聲。

表4 p的取值對聚類數目的影響Table 4 The effects of p on the number of clustering centers

3.5 基于相關距離的OPTICS算法的流場結構特征識別

基于相關距離的OPTICS算法的初始參數與原始OPTICS算法一致，區別僅在于其鄰域半徑ε用相關距離表達而非歐氏距離。相關距離可以根據夾角余弦、Pearson相關系數或者Spearman相關系數構造，本節將比較這些不同的相似度指標對基于相關距離的OPTICS算法的影響。

圖10列出了基于Pearson相關距離的OPTICS算法獲得的決策圖和對應的聚類中心。對于全部的MinPts的取值，均可以找到適宜的鄰域半徑ε，對渦量場樣本的特征實施聚類。當MinPts取60和80時，能夠找到6個明顯的聚類中心，其中1～4號聚類對應于間距約為5d的順流向渦對的脫落狀態，由于z= 0.08 m處渦核形狀存在差異，出現了過度分類的情況；5～6號聚類分別對應于間距約為4d的順流向渦對的兩個脫落狀態。兩種間距的順流向渦均可以被視為A模式，其間距差別源于流動的動力不穩定。值得注意的是，前四種聚類中心與后兩種聚類中心需要分別通過兩種鄰域半徑ε1和ε2在決策圖中獲得，這就表明，獲取前四種聚類與后兩種聚類的密度評估標準不同，也即意味著聚類結果的判別標準不統一。當MinPts取100和120時，能夠利用同一鄰域半徑ε找到4個明顯的聚類中心，分別對應于間距約為5d和4d的順流向渦對的兩個脫落狀態。此時相對于MinPts< 100的情況，在保證特征識別準確的同時，聚類數目也更加合理。當MinPts進一步增加，聚類中心縮減為兩個，分別對應于順流向渦脫落的兩個狀態，但沒有能夠區分兩種渦對展向間距的情況，這也導致聚類中心呈現的渦對出現了不清晰的情況。可見，當MinPts取100～120時，基于Pearson相關距離的OPTICS算法能夠通過合適的鄰域半徑，識別到順流向渦結構的展向分布間距和旋渦脫落狀態。

圖11列出了基于夾角余弦相關距離的OPTICS算法獲得的決策圖和對應的聚類中心。對比圖10和圖11，會發現基于夾角余弦的流場結構識別結果與基于Pearson相關距離的流場結構識別結果高度相似，參數規律也完全相同。區別僅在于部分聚類的樣本數目有微小差距。由表1中兩者的定義，可以知道，Pearson相關系數相當于排除均值后樣本的夾角余弦，所以采用兩種指標計算的核心距離與可達距離應當具備相同的規律，僅在具體數值上有所區別。同時，對于本文選取的渦量場樣本，每個樣本中除了渦出現的位置存在顯著的渦量值，其他大范圍的背景渦量值接近于0，這導致各個渦量場樣本的均值也接近于0，從而夾角余弦和Pearson相關系數對于渦量場樣本的計算結果也差別極小。采用兩種指標進行基于OPTICS算法的聚類，其決策圖和聚類結果也應當十分接近。

圖10 基于Pearson相關距離的OPTICS算法的決策圖與聚類中心Fig. 10 The reachability-plots and cluster centers of the OPTICS algorithm based on the Pearson-correlation distance

圖11 基于夾角余弦相關距離的OPTICS算法的決策圖與聚類中心Fig. 11 The reachability-plots and cluster centers of the OPTICS algorithm based on the cosine-correlation distance

圖12列出了基于Spearman相關距離的OPTICS算法獲得的決策圖和對應的聚類中心。Spearman相關系數僅考慮了流場樣本中渦量值的排序關系，而忽略了其具體的數值的影響，這對OPTICS算法的聚類結果產生了明顯影響。當MinPts取60和80時，能夠找到8個明顯的聚類中心，但其中僅能清晰識別到間距約為4d的順流向渦對的脫落狀態，而其他的聚類結果并不能清晰地體現渦量場的結構特征。當MinPts= 100時，可以識別出四個聚類中心，分別對應于間距約為5d和4d的順流向渦對的兩個脫落狀態；其中間距為5d的順流向渦的脫落狀態不如采用Pearson相關系數或夾角余弦時的結果清晰。而MinPts> 100時，可以找到對應于間距約為5d和4d的兩種順流向渦分布，但無法區分旋渦脫落的狀態。這樣的結果證明，基于Spearman相關距離的OPTICS算法并不適于流場結構特征的分析。

圖12 基于Spearman相關距離的OPTICS算法的決策圖與聚類中心Fig. 12 The reachability-plots and cluster centers of the OPTICS algorithm based on the Spearman-correlation distance

表5列出了本節討論的各種聚類算法用于流場結構特征分析算例的結果對比，結合具體的聚類中心分析結果，基于Pearson相關距離和基于夾角余弦的OPTICS算法在本算例的流場結構分析中獲得了準確的分析結果。

表5 不同流場結構特征分析聚類算法的對比Table 5 A comparison of among the different clustering algorithms for analyzing flow structure characteristics

4 基于相關距離的OPTICS聚類算法進行流場結構特征識別的魯棒性分析

本節考慮基于前述數值模擬獲取的流場樣本實例，通過在樣本中加入噪聲的方式，分析樣本噪聲強度和噪聲樣本的數量對聚類分析結果魯棒性的影響。

4.1 測試樣本的構造方法與聚類結果的評判指標

前述數值模擬共計獲得2000個流場樣本，在這些樣本中隨機抽取一定數量的樣本，按下式添加高斯噪聲：

式中：xi′表示流場樣本X中已添加噪聲的樣本，xi表示流場樣本X中被隨機選中的待添加噪聲的原始樣本，e表示均值為0標準差為1的高斯隨機序列，δ則為高斯隨機序列的縮放系數，可以表征噪聲的強度。測試樣本由未被抽取的原始樣本和隨機抽取后添加噪聲的樣本混合構造，每種參數工況按3種不同的隨機序列分別構造3個測試樣本。

考慮到數值模擬獲得的渦量場樣本，其渦量為101～102量級，所以本節探討以下7種δ取值：0.1、0.3、1、3、10、30、100，它們對應的噪聲樣本如圖13所示。隨機抽取用以添加噪聲的樣本數目則利用樣本代換率Rs表示：

圖13 添加不同強度噪聲的流場樣本示例Fig. 13 The flow fields with different intensities of the Gaussian-noise

本節討論的Rs取值情況為：1%、2%、5%、15%和50%。

由于流場本身動力系統的復雜性，較難提取到標準的流場結構特征聚類。所以聚類分析結果的評判，要基于測試樣本聚類結果相對于原始樣本聚類結果的偏移，即衡量測試樣本聚類中心和原始樣本聚類中心的相關距離。這一衡量指標可以按下式定義為各個聚類中心偏移量的算術平均：

式中：dcenter為基于相關距離的聚類中心偏移指標，N為聚類數目，dc′,c,i表示測試樣本C′與原始樣本C的第i個聚類中心的相關距離。在原始樣本的聚類結果中（圖10），MinPts= 100且ε= 0.644的結果可以作為用于衡量聚類結果魯棒性的基準聚類中心。

4.2 噪聲強度δ對聚類結果的影響規律

圖14列出了各種樣本代換率下，dcenter隨噪聲強度δ的變化規律。可見針對所有的Rs設置，dcenter均隨δ的增長而增加；在對數坐標下，dcenter的增長隨著δ的增長而趨緩。這表明，當噪聲強度較大時，噪聲強度對聚類中心的偏移影響減小。當樣本代換率Rs≤5%，噪聲強度造成的聚類中心偏移的相關距離普遍不超過0.05，也即意味著此時測試樣本的聚類中心與原始樣本的聚類中心高度相關，聚類結果穩定。此時，流場樣本的噪聲對流場結構識別結果幾乎沒有影響。當樣本代換率Rs≥15%，在引入較小噪聲（δ≤10）時，測試樣本的聚類中心偏移量最大約為0.10，此時測試樣本的聚類中心與原始樣本的聚類中心仍然高度相關。但當引入較大噪聲（δ≥30）時，部分測試樣本已經無法被識別到對應于原始樣本聚類中心的四種聚類，因此在圖14（d）和圖14（e）中缺失了部分樣本δ≥30的數據點，此時基于相關距離的OPTICS聚類算法無法在原有的初始參數設置下實現流場結構特征分析識別。

圖14 噪聲強度δ對聚類中心偏移量dcenter的影響Fig. 14 Effects of the noise intensity (δ) on the displacements of the clustering centers (dcenter)

4.3 樣本代換率Rs對聚類結果的影響規律

圖15列出了噪聲強度0.1≤δ≤10時，dcenter隨樣本代換率Rs的變化規律。針對所有的δ設置，dcenter均隨Rs的增長而呈現大致增長的趨勢；在對數坐標下，dcenter的增長率隨著δ的增加而下降。這表明，當噪聲強度增加時，樣本代換率對聚類中心的偏移影響減小。根據圖15，當噪聲強度有限時，樣本代換率幾乎不影響基于相關距離的OPTICS算法的流場結構特征識別的結果。

圖15 樣本代換率Rs對聚類中心偏移量dcenter的影響Fig. 15 Effects of the sample-replacement ratio (Rs) on the displacements of the clustering centers (dcenter)

5 結論

本文基于OPTICS聚類算法，在分析流場數據特征與多種相似度評判指標的基礎上，引入相關距離的概念加以改進，提出利用基于相關距離的OPTICS算法進行流場結構特征分析。該方法依托基于Pearson相關系數的相關距離指標，不需要動力學假設，僅依靠數據驅動。本文通過基于LES的CFD數值模擬，對經典的圓柱繞流問題進行瞬態計算，獲取了2000個流場樣本。以識別順流向渦的旋渦脫落狀態和A模式分布為目標，對比了k-means、原始的OPTICS、基于閔可夫斯基距離的OPTICS和基于相關距離的OPTICS等聚類算法的有效性，并對基于Pearson相關距離的OPTICS識別流場結構特征的魯棒性進行了分析。結果表明：

1）相對于k-means算法、原始的OPTICS算法和基于閔可夫斯基距離的OPTICS算法，基于Pearson相關距離的OPTICS算法可以有效識別出圓柱尾流中順流向渦脫落的兩種狀態，以及A模式順流向渦的兩種不同的展向間距。針對本文的算例，基于夾角余弦相關距離的OPTICS算法也可以獲得與基于Pearson相關距離的OPTICS算法高度相似的結果，其原因在于算例樣本的元素均值幾乎為0，樣本元素均值對于夾角余弦的影響很小；而基于Spearman相關距離的OPTICS算法的識別結果相對不夠清晰，適宜的初始參數范圍較小。

2）基于Pearson相關距離的OPTICS算法，可以在MinPts取100～120時，通過設置合適的鄰域半徑，識別到順流向渦結構的展向分布間距和旋渦脫落狀態。MinPts取值小于100會導致流場結構聚類結果的過度分類，而MinPts取值大于120時會失去識別順流向渦展向間距的能力。

3）相對于k-means算法，基于各種相似度指標的OPTICS算法均提升了聚類結果的穩定性，降低了聚類結果對初始參數的敏感性，保證了在確定的初始參數下，聚類結果唯一。相關距離概念的引入則改善了原始OPTICS算法對流場結構特征的識別能力，其中以Pearson相關系數和夾角余弦為指標計算的相關距離效果最佳。

4）基于Pearson相關距離的OPTICS算法，當初始參數設置為MinPts= 100時，在噪聲樣本量少于總樣本量的5%時，可以在高強度高斯噪聲（δ≤100）樣本的干擾下，依然識別到清晰的流場結構特征，噪聲樣本造成的聚類中心偏移的相關距離不超過0.05；而當噪聲樣本量超過樣本總量15%時，基于相關距離的OPTICS算法可以在引入δ≤10的噪聲強度下實現較穩定的流場結構特征識別，噪聲造成的聚類中心偏移的相關距離不超過0.10。