非結(jié)構(gòu)化網(wǎng)格下海洋流場(chǎng)的特征提取與種子點(diǎn)選取算法

李忠偉，宮凱旋，李永，劉格格

（1.中國(guó)石油大學(xué)（華東）海洋與空間信息學(xué)院，山東青島 266580；2.中國(guó)石油大學(xué)（華東）計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，山東 青島 266580）

0 引言

隨著計(jì)算機(jī)性能的不斷提高，研究者通過海洋數(shù)值模式對(duì)多種時(shí)空尺度下的海洋現(xiàn)象和海洋狀況進(jìn)行了分析和預(yù)報(bào)，海洋數(shù)值模式也正朝著更高分辨率和更快計(jì)算速度的方向逐步發(fā)展［1］。這些海洋數(shù)值模式產(chǎn)生了海量的數(shù)據(jù)集合，通常以網(wǎng)格形式存儲(chǔ)。網(wǎng)格可以分為結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格［2］，對(duì)于以結(jié)構(gòu)化網(wǎng)格存儲(chǔ)的海洋流場(chǎng)數(shù)據(jù)，因其規(guī)則的拓?fù)浣Y(jié)構(gòu)，很容易對(duì)其可視化以展示流場(chǎng)信息；而非結(jié)構(gòu)化網(wǎng)格沒有規(guī)則的拓?fù)浣Y(jié)構(gòu)，因此可以更好地?cái)M合復(fù)雜的邊界，這個(gè)優(yōu)點(diǎn)使其在研究島嶼和近岸岸線復(fù)雜問題時(shí)的表現(xiàn)尤為突出［3］。但也正因?yàn)榉墙Y(jié)構(gòu)化網(wǎng)格不具備規(guī)則的拓?fù)浣Y(jié)構(gòu)，在點(diǎn)定位和搜索相鄰網(wǎng)格時(shí)都比較困難，給繪制非結(jié)構(gòu)化網(wǎng)格下的流場(chǎng)數(shù)據(jù)帶來阻礙。在處理非結(jié)構(gòu)化網(wǎng)格時(shí)，有的方法會(huì)先遍歷整個(gè)流場(chǎng)區(qū)域的網(wǎng)格及格點(diǎn)信息，構(gòu)建出格點(diǎn)、網(wǎng)格的相互關(guān)系［4］后再進(jìn)行運(yùn)算，這個(gè)流場(chǎng)網(wǎng)格關(guān)系重構(gòu)過程會(huì)花費(fèi)過多的時(shí)間；此外，流場(chǎng)可視化時(shí)所需要的繪制數(shù)據(jù)（如構(gòu)成流線的軌跡點(diǎn)），基本都是以格點(diǎn)數(shù)據(jù)插值得到的［5-6］，而對(duì)于海流數(shù)據(jù)存儲(chǔ)在網(wǎng)格格心的流場(chǎng)數(shù)據(jù)，必須通過插值或外推先得到格點(diǎn)處的數(shù)據(jù)［7］，再進(jìn)行流場(chǎng)可視化，這種間接的操作會(huì)導(dǎo)致部分信息丟失以至于引發(fā)流場(chǎng)中特征區(qū)域范圍的變化［8］，造成精度損失。因此，如何快速準(zhǔn)確地使用非結(jié)構(gòu)化三角網(wǎng)格格心流場(chǎng)數(shù)據(jù)，以及構(gòu)造網(wǎng)格、格點(diǎn)及格心之間的拓?fù)潢P(guān)系，避免外推或插值所造成的精度損失，是非結(jié)構(gòu)化網(wǎng)格流場(chǎng)的特征提取和可視化工作必須要解決的難題。

基于特征的可視化技術(shù)是指通過繪制有意義的對(duì)象和特征結(jié)構(gòu)來描述流場(chǎng)，其首要任務(wù)是提取流場(chǎng)數(shù)據(jù)特征，常見方法是由HELMAN 等［9］提出的矢量場(chǎng)拓?fù)浞治龇椒?，即主要是基于矢量?chǎng)臨界點(diǎn)的檢測(cè)和分類。臨界點(diǎn)及其附近的區(qū)域蘊(yùn)含著豐富的流場(chǎng)特征，如渦旋、鞍點(diǎn)區(qū)域等，因此，臨界點(diǎn)的提取和分類作為流場(chǎng)特征提取的關(guān)鍵，國(guó)內(nèi)外學(xué)者也對(duì)其進(jìn)行了一系列研究。LAVIN 等［10］在結(jié)構(gòu)化網(wǎng)格上通過線性插值的方式定位臨界點(diǎn)，吳曉莉等［11］則將Sperner 原理引入到臨界點(diǎn)檢測(cè)中，季民等［12］綜合改進(jìn)了雙線性插值和Sperner 完全標(biāo)號(hào)法，全面準(zhǔn)確地提取了海洋流場(chǎng)中的臨界點(diǎn)，BHATIA 等［13］通過判斷單純形頂點(diǎn)處的矢量值是否構(gòu)成一個(gè)內(nèi)部包含原點(diǎn)的單純形來魯棒地判定單純形內(nèi)是否存在臨界點(diǎn)。目前，二維流場(chǎng)臨界點(diǎn)提取方法已趨近于成熟，需要對(duì)比選擇一種能夠準(zhǔn)確且適用于非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)的臨界點(diǎn)提取方法。

不同的臨界點(diǎn)附近具有不同的流場(chǎng)特征，因此，可以將臨界點(diǎn)分類定位到更具體的流場(chǎng)特征。在結(jié)構(gòu)化網(wǎng)格上，一般通過求解臨界點(diǎn)處的雅克比矩陣來將臨界點(diǎn)分類，而非結(jié)構(gòu)化網(wǎng)格由于不同區(qū)域的網(wǎng)格大小不一，dx和dy不固定，很難構(gòu)造恰當(dāng)?shù)难趴吮染仃?，更多的是采用一些拓?fù)浞椒▉矸诸?，如WANG等［14］在改進(jìn)Bhatia 的臨界點(diǎn)檢測(cè)方法以更快檢測(cè)臨界點(diǎn)的同時(shí)，通過龐加萊指數(shù)將臨界點(diǎn)分為鞍點(diǎn)和其他類型的臨界點(diǎn)，從而實(shí)現(xiàn)了鞍點(diǎn)區(qū)域的快速定位。但若想對(duì)非結(jié)構(gòu)化網(wǎng)格流場(chǎng)進(jìn)一步分析，還需要將臨界點(diǎn)類型進(jìn)一步細(xì)分，即在非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)的臨界點(diǎn)處構(gòu)造合適的雅克比矩陣。

在流場(chǎng)可視化技術(shù)中，通常用流線來表現(xiàn)流場(chǎng)，流線可以有效描述流場(chǎng)中粒子的運(yùn)動(dòng)軌跡［15］。在計(jì)算流線時(shí)，首先在流場(chǎng)中放置種子點(diǎn)，然后再考慮每個(gè)軌跡點(diǎn)處的矢量值并積分軌跡點(diǎn)的路徑以生成流線［16］。因此，種子點(diǎn)的選取在很大程度上決定了流線生成的質(zhì)量及可視化效果［17］，也一直是流場(chǎng)可視化的熱門研究領(lǐng)域。DOVEY 等［18］在規(guī)則網(wǎng)格上均勻放置種子點(diǎn)來生成流線，但這種方法生成的流線并不會(huì)均勻分布。TREINISH 等［19］將近似臨界點(diǎn)視為種子點(diǎn)并計(jì)算出一組初始流線，再識(shí)別矢量場(chǎng)中矢量值變化較大的區(qū)域，在這些區(qū)域中放置種子點(diǎn)以計(jì)算額外的流線。VERMA 等［20］根據(jù)臨界點(diǎn)類型，使用特定的種子點(diǎn)模板來捕捉臨界點(diǎn)附近的流動(dòng)行為。近年來，許多基于原始流場(chǎng)導(dǎo)出標(biāo)量場(chǎng)的方法被提出來定量地分析流場(chǎng)特征區(qū)域［21］。巴振宇等［22］將臨界點(diǎn)與信息熵有機(jī)地結(jié)合起來，先用信息熵找到流場(chǎng)變化劇烈的區(qū)域，再利用臨界點(diǎn)的拓?fù)浣Y(jié)構(gòu)信息對(duì)局部區(qū)域的種子點(diǎn)放置做精細(xì)控制。黃冬梅等［23］設(shè)計(jì)了兩種基于信息熵理論的種子點(diǎn)放置算法：貪婪策略算法和蒙特卡羅算法，這兩種算法先生成均勻分布的種子點(diǎn)，再利用流場(chǎng)信息熵篩選區(qū)域中包含較多流場(chǎng)特征的種子點(diǎn)。DU 等［24］基于四叉樹結(jié)構(gòu)和信息熵選擇種子點(diǎn)，顯著地提高了種子點(diǎn)選擇的速度。從上述種子點(diǎn)選取算法可以看出，種子點(diǎn)的選擇大多都關(guān)注于流場(chǎng)關(guān)鍵特征，因此，很多種子點(diǎn)被放置在臨界點(diǎn)附近，但這可能會(huì)遺漏一些其他位置的流場(chǎng)信息，造成整體的視覺效果較差。相比之下，通過原始流場(chǎng)導(dǎo)出的標(biāo)量場(chǎng)來定位流場(chǎng)特征可以更好地選擇種子點(diǎn)，生成的流線也可以更好地覆蓋流場(chǎng)特征區(qū)域。對(duì)于非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)，因其不規(guī)則的拓?fù)潢P(guān)系，很難由原始流場(chǎng)計(jì)算定位流場(chǎng)特征的標(biāo)量場(chǎng)，所以也無法進(jìn)一步定量地分析流場(chǎng)特征區(qū)域以選擇種子點(diǎn)。

本文為解決上述非結(jié)構(gòu)化網(wǎng)格流場(chǎng)特征提取和可視化過程中存在的問題，在流場(chǎng)數(shù)據(jù)預(yù)處理、臨界點(diǎn)的定位及分類和種子點(diǎn)的選取上分別進(jìn)行改進(jìn)創(chuàng)新。首先，通過Delaunay 三角剖分算法重構(gòu)流場(chǎng)網(wǎng)格格心的拓?fù)潢P(guān)系，從而直接利用原始流場(chǎng)的格心海流數(shù)據(jù)來提取流場(chǎng)特征；其次，基于重構(gòu)的流場(chǎng)數(shù)據(jù)，選擇適用于非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)的龐加萊指數(shù)法來提取臨界點(diǎn)并構(gòu)造臨界點(diǎn)處的雅克比矩陣，從而將臨界點(diǎn)細(xì)分；最后，改進(jìn)“最大得分”法，使其適用于非結(jié)構(gòu)化三角網(wǎng)格，在非結(jié)構(gòu)化網(wǎng)格上計(jì)算導(dǎo)出“最大得分”標(biāo)量場(chǎng)來定位流場(chǎng)特征區(qū)域，根據(jù)重構(gòu)網(wǎng)格的密度屬性設(shè)置閾值來選擇種子點(diǎn)，并通過臨界點(diǎn)來驗(yàn)證由種子點(diǎn)生成的流線能否很好地表現(xiàn)出流場(chǎng)特征。最后，通過對(duì)多組實(shí)驗(yàn)結(jié)果分析，證明本文方法在流場(chǎng)特征提取與提升可視化效果工作中的準(zhǔn)確性和有效性。

1 流場(chǎng)數(shù)據(jù)預(yù)處理

1.1 數(shù)據(jù)來源

本文研究所用的海洋流場(chǎng)數(shù)據(jù)來自于中山大學(xué)的全球Global-FVCOM 海洋數(shù)值預(yù)報(bào)系統(tǒng)，模式水平分辨率在北極地區(qū)近岸約2 km，在陸架區(qū)域逐漸增長(zhǎng)到10 km 及大洋中心至25 km，包括359 191 個(gè)三角形計(jì)算網(wǎng)格點(diǎn)及689 133 個(gè)有限體積計(jì)算單元，是目前國(guó)內(nèi)分辨率最高的全球海洋預(yù)報(bào)模式。

數(shù)據(jù)格式為NetCDF，經(jīng)向流速u和緯向流速v存儲(chǔ)在三角網(wǎng)格格心位置，時(shí)間為2016 年9 月29 日，為全球海洋表面流。海洋流場(chǎng)數(shù)據(jù)解析如表1 所示，其中包含三角網(wǎng)格的頂點(diǎn)信息以及存儲(chǔ)在三角網(wǎng)格質(zhì)心的矢量信息，nv變量代表構(gòu)成三角網(wǎng)格的3 個(gè)頂點(diǎn)索引，由頂點(diǎn)索引信息可以從lon、lat變量中得到3 個(gè)頂點(diǎn)的坐標(biāo)。

表1 原始流場(chǎng)數(shù)據(jù)解析Table 1 Analysis of original flow field data 單位：個(gè)

1.2 非結(jié)構(gòu)化三角網(wǎng)格流場(chǎng)數(shù)據(jù)重構(gòu)

為避免將格心矢量數(shù)據(jù)外推或插值到格點(diǎn)再計(jì)算流場(chǎng)特征所造成的精度損失，本文使用Delaunay三角剖分算法重構(gòu)非結(jié)構(gòu)化三角網(wǎng)格，即將岸線點(diǎn)和格心點(diǎn)作為給定點(diǎn)集對(duì)海洋區(qū)域三角剖分，重新構(gòu)造了格心間的拓?fù)潢P(guān)系，新網(wǎng)格點(diǎn)由原始網(wǎng)格的岸線點(diǎn)和格心點(diǎn)構(gòu)成。

岸線點(diǎn)提取示意圖見圖1（彩色效果見《計(jì)算機(jī)工程》官網(wǎng)HTML 版，下同）。首先提取原始流場(chǎng)數(shù)據(jù)的岸線點(diǎn)，如圖1（a）所示。由于構(gòu)成岸線的網(wǎng)格邊僅會(huì)出現(xiàn)一次，因此本文通過三角網(wǎng)格頂點(diǎn)集nv來提取岸線點(diǎn)。首先將每個(gè)三角網(wǎng)格的頂點(diǎn)索引組合成3 個(gè)向量，即3 條邊，并刪除重復(fù)出現(xiàn)的向量，然后獲得構(gòu)成非公共向量的網(wǎng)格點(diǎn)索引集合，最后根據(jù)索引從原始流場(chǎng)信息lon、lat變量中得到岸線點(diǎn)的坐標(biāo)，如圖1（b）所示。

圖1 岸線點(diǎn)提取Fig.1 Extraction of shoreline points

接著由格心點(diǎn)和提取到的岸線點(diǎn)來重構(gòu)網(wǎng)格并生成網(wǎng)格拓?fù)潢P(guān)系。原始網(wǎng)格及格心點(diǎn)如圖2（a）所示，圖中包括岸線點(diǎn)、三角網(wǎng)格和格心點(diǎn)，將格心點(diǎn)和岸線點(diǎn)作為給定點(diǎn)集來構(gòu)造新的流場(chǎng)網(wǎng)格。接下來通過逐點(diǎn)插入法生成Delaunay 三角網(wǎng)格，最后由岸線點(diǎn)和格心點(diǎn)構(gòu)造出全新的三角網(wǎng)格，如圖2（b）所示?？梢钥吹剑貥?gòu)后的網(wǎng)格對(duì)岸線及島嶼的邊緣擬合并不會(huì)改變，并且將原始流場(chǎng)的格心直接相連。同時(shí)，在流場(chǎng)網(wǎng)格重構(gòu)的過程中獲得了完整的流場(chǎng)拓?fù)湫畔?，具體拓?fù)湫畔⑷绫? 所示，這些拓?fù)湫畔⒖梢源蠓群?jiǎn)化點(diǎn)定位和查詢相鄰格點(diǎn)或相鄰三角網(wǎng)格等操作，為后續(xù)流場(chǎng)特征的提取與可視化研究奠定基礎(chǔ)。

圖2 原始流場(chǎng)網(wǎng)格和重構(gòu)后的網(wǎng)格對(duì)比Fig.2 Comparison of original flow field grid and the reconstructed grid

表2 重構(gòu)后網(wǎng)格的拓?fù)湫畔able 2 Topology information of the reconstructed grid

2 本文方法

2.1 臨界點(diǎn)提取及分類

2.1.1 臨界點(diǎn)理論

臨界點(diǎn)是指流場(chǎng)中經(jīng)向流速u和緯向流速v均為0 的點(diǎn)，又被稱為奇點(diǎn)或固定點(diǎn)。臨界點(diǎn)及其附近區(qū)域有著豐富的流場(chǎng)特征，如渦旋、鞍點(diǎn)等。臨界點(diǎn)可根據(jù)該點(diǎn)處的雅克比矩陣是否等于零分為非孤立臨界點(diǎn)和孤立臨界點(diǎn)，本文主要研究的是孤立臨界點(diǎn)。對(duì)于孤立臨界點(diǎn)，可以用臨界點(diǎn)處的雅克比矩陣來表示矢量場(chǎng)及其附近流線的行為，其公式如下：

臨界點(diǎn)的類型可以通過雅克比矩陣特征值的實(shí)部和虛部來進(jìn)行判斷，實(shí)部的正負(fù)決定臨界點(diǎn)附近粒子的運(yùn)動(dòng)方向，虛部則決定臨界點(diǎn)附近是否存在旋轉(zhuǎn)行為。具體分類情況如圖3 所示（圖中R表示實(shí)部，I表示虛部）。據(jù)此，臨界點(diǎn)可以分為吸引節(jié)點(diǎn)、排斥節(jié)點(diǎn)、吸引聚點(diǎn)、排斥聚點(diǎn)、馬鞍點(diǎn)、中心點(diǎn)等6 類。

圖3 流場(chǎng)臨界點(diǎn)分類Fig.3 Classification of critical points of the flow field

2.1.2 龐加萊指數(shù)法

在提取非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)的臨界點(diǎn)時(shí)，本文對(duì)比分析了龐加萊指數(shù)法［25］和Sperner 完全標(biāo)號(hào)法的臨界點(diǎn)提取效果，并選擇前者來提取臨界點(diǎn)。

龐加萊指數(shù)是矢量場(chǎng)拓?fù)渲械闹匾拍?，根?jù)式（2）、式（3）可求得簡(jiǎn)單閉合曲線的龐加萊指數(shù)：

其中，?是矢量場(chǎng)的角坐標(biāo)。

龐加萊指數(shù)始終是一個(gè)整數(shù)，當(dāng)一個(gè)簡(jiǎn)單閉曲線中包含臨界點(diǎn)時(shí)，該閉曲線的龐加萊指數(shù)可能為-1（鞍點(diǎn)）或+1（其他類型），因此，可將三角網(wǎng)格的邊視為簡(jiǎn)單閉曲線，并且由于在流場(chǎng)中單個(gè)網(wǎng)格中僅存在一個(gè)臨界點(diǎn)，即孤立臨界點(diǎn)，可以得出，當(dāng)且僅當(dāng)三角網(wǎng)格的龐加萊指數(shù)為+1 或-1 時(shí)，該三角網(wǎng)格中存在臨界點(diǎn)。

三角網(wǎng)格的龐加萊指數(shù)可以由3 個(gè)頂點(diǎn)處的角度變化之和來計(jì)算，設(shè)?0、?1和?2分別對(duì)應(yīng)定義在三角形T頂點(diǎn)和的角坐標(biāo)（0 到2π 之間）。

三角形T的龐加萊指數(shù)可以由式（4）、式（5）計(jì)算：

2.1.3 三角網(wǎng)格內(nèi)臨界點(diǎn)定位

為提取非結(jié)構(gòu)化網(wǎng)格的臨界點(diǎn)，本文提出三角網(wǎng)格質(zhì)心迭代法以定位臨界點(diǎn)在三角網(wǎng)格中的準(zhǔn)確位置。通過求包含臨界點(diǎn)的三角網(wǎng)格的質(zhì)心，并使用三角形3 個(gè)頂點(diǎn)的矢量值重心插值得到該質(zhì)心處的矢量值，繼續(xù)計(jì)算由質(zhì)心和父三角形頂點(diǎn)所構(gòu)成的3 個(gè)子三角形的龐加萊指數(shù)，判斷臨界點(diǎn)所在的子三角形，如此迭代來逼近臨界點(diǎn)的位置。

臨界點(diǎn)定位步驟如圖4 所示，其中實(shí)線三角形表示當(dāng)前包含臨界點(diǎn)的三角形，不斷迭代質(zhì)心最后得到質(zhì)心v2作為臨界點(diǎn)。

圖4 質(zhì)心迭代臨界點(diǎn)定位Fig.4 Location of critical point by centroid iteration

臨界點(diǎn)定位算法具體步驟如下：

步驟1確定包含臨界點(diǎn)的三角網(wǎng)格，頂點(diǎn)分別為a、b、c，并設(shè)定一個(gè)閾值。

步驟2求得這個(gè)三角形的質(zhì)心點(diǎn)v0，并通過三角形重心坐標(biāo)插值出該點(diǎn)的u、v速度分量，求質(zhì)心處矢量的平方范數(shù)（u2+v2）是否小于閾值，若小于給定閾值，判定該質(zhì)心點(diǎn)為臨界點(diǎn)，算法終止，否則繼續(xù)。

步驟3使用質(zhì)心點(diǎn)依次替換三角形的3 個(gè)頂點(diǎn)，形成子三角形v0ab、v0bc和av0c。

步驟4分別計(jì)算這3 個(gè)子三角形的龐加萊指數(shù)，找到包含臨界點(diǎn)的子三角形。

步驟5當(dāng)子三角形包含臨界點(diǎn)時(shí)，使用子三角形的3 個(gè)頂點(diǎn)再次重新插值出質(zhì)心點(diǎn)的u、v速度分量。

步驟6判斷子三角形質(zhì)心處矢量的平方范數(shù)是否小于閾值，當(dāng)小于閾值時(shí)判定該質(zhì)心點(diǎn)為臨界點(diǎn)，否則跳轉(zhuǎn)到步驟3。

2.1.4 三角網(wǎng)格內(nèi)臨界點(diǎn)定位

在非結(jié)構(gòu)化網(wǎng)格的臨界點(diǎn)處構(gòu)造雅克比矩陣將其準(zhǔn)確分類是一個(gè)難題，本文設(shè)計(jì)了非結(jié)構(gòu)化三角網(wǎng)格下構(gòu)造雅克比矩陣的方法，從而將臨界點(diǎn)類型進(jìn)一步細(xì)分。在2.1.3 節(jié)臨界點(diǎn)定位算法的基礎(chǔ)上，記錄最后包含臨界點(diǎn)的子三角形的頂點(diǎn)坐標(biāo)信息，然后進(jìn)行以下2 個(gè)步驟：

步驟1如圖5（a）所示，將三角形頂點(diǎn)距離臨界點(diǎn)水平和垂直方向的最遠(yuǎn)距離分別設(shè)定為dx和dy，確保臨界點(diǎn)位于區(qū)域范圍內(nèi)，從而構(gòu)造的雅克比矩陣可以體現(xiàn)出臨界點(diǎn)周圍的運(yùn)動(dòng)趨勢(shì)。

圖5 三角網(wǎng)格內(nèi)雅克比矩陣構(gòu)造Fig.5 Construction of Jacobian matrix in triangular grid

步驟2依次插值出臨界點(diǎn)上、下、左、右4 個(gè)點(diǎn)的矢量值，如圖5（b）所示，再根據(jù)這4 個(gè)點(diǎn)的信息中心差分得到臨界點(diǎn)處的雅克比矩陣，計(jì)算雅克比矩陣的特征值，細(xì)分臨界點(diǎn)類型。

2.2 基于“最大得分”法和網(wǎng)格密度的種子點(diǎn)選取

2.2.1 “最大得分”法

由FERRARI 等［26］提出的“最大得分”法，能夠基于結(jié)構(gòu)化網(wǎng)格下的相鄰網(wǎng)格節(jié)點(diǎn)的相對(duì)標(biāo)量值識(shí)別出重要網(wǎng)格節(jié)點(diǎn)，有效地定位流場(chǎng)特征。對(duì)于給定標(biāo)量場(chǎng)，該方法根據(jù)相對(duì)于鄰居網(wǎng)格點(diǎn)的強(qiáng)度水平對(duì)每個(gè)點(diǎn)進(jìn)行評(píng)分。假設(shè)給定標(biāo)量域中的點(diǎn)p0被其相鄰節(jié)點(diǎn)p1～pn所包圍，則節(jié)點(diǎn)p0的最大得分計(jì)算為：

但由于非結(jié)構(gòu)化網(wǎng)格缺乏完整的拓?fù)湫畔?，在求某一點(diǎn)的相鄰網(wǎng)格點(diǎn)時(shí)，需要耗費(fèi)大量時(shí)間，使得該方法目前僅適用于結(jié)構(gòu)化網(wǎng)格。本文在1.2 節(jié)在重構(gòu)流場(chǎng)拓?fù)湫畔⒌倪^程中建立了新網(wǎng)格頂點(diǎn)及其相鄰網(wǎng)格頂點(diǎn)的關(guān)系，因此，選擇三角網(wǎng)格中的某一頂點(diǎn)與其雙層鄰接頂點(diǎn)，通過比較它們之間的相對(duì)標(biāo)量值來計(jì)算該頂點(diǎn)的得分。結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格格點(diǎn)及相鄰格點(diǎn)分別如圖6（a）、圖6（b）所示，其中，紅色點(diǎn)表示選定的格點(diǎn)p0，綠色點(diǎn)表示其雙層鄰接頂點(diǎn)，本文計(jì)算重構(gòu)后流場(chǎng)中每個(gè)格點(diǎn)的評(píng)分，并形成一個(gè)“最大得分”標(biāo)量場(chǎng)用于定量分析流場(chǎng)特征。

2.2.2 種子點(diǎn)選取

本文首先計(jì)算原始流場(chǎng)格心處矢量的平方范數(shù)，即u2+v2，作為給定標(biāo)量場(chǎng)計(jì)算原始流場(chǎng)的“最大得分”場(chǎng)，通過這個(gè)“最大得分”場(chǎng)定位流場(chǎng)特征。但筆者在設(shè)置“最大得分”場(chǎng)閾值選取種子點(diǎn)的過程中發(fā)現(xiàn)，非結(jié)構(gòu)化三角網(wǎng)格在近岸區(qū)域網(wǎng)格更為密集，在遠(yuǎn)岸區(qū)域網(wǎng)格則較為稀疏，導(dǎo)致部分位于稀疏區(qū)域的頂點(diǎn)距其雙層鄰接點(diǎn)過遠(yuǎn)，這些網(wǎng)格稀疏區(qū)域的格點(diǎn)得分可能較低但卻存在流場(chǎng)特征，設(shè)置過高的閾值可能會(huì)錯(cuò)過這些區(qū)域。因此，考慮將網(wǎng)格密度屬性加入到閾值的設(shè)置中，保證在非結(jié)構(gòu)化三角網(wǎng)格流場(chǎng)的任何區(qū)域，都能通過“最大得分”場(chǎng)定位到流場(chǎng)特征，選取種子點(diǎn)。

首先將網(wǎng)格密度屬性設(shè)定為格點(diǎn)到其鄰居格點(diǎn)的平均距離，網(wǎng)格密度越大，相鄰格點(diǎn)的平均距離越小。按照格點(diǎn)的相鄰平均距離從小到大劃分為n+1 個(gè)區(qū)間，分別設(shè)置選取種子點(diǎn)的最大閾值tmax和最小閾值tmin，則點(diǎn)pi的閾值t為：

其中，n（pi）代表點(diǎn)pi所處的密度區(qū)間，最小值為0，最大值為n。只有當(dāng)點(diǎn)pi的最大得分大于閾值t，才會(huì)將該點(diǎn)加入種子點(diǎn)。

3 實(shí)驗(yàn)結(jié)果及分析

本文對(duì)非結(jié)構(gòu)化網(wǎng)格格心流場(chǎng)數(shù)據(jù)預(yù)處理、臨界點(diǎn)定位及分類和種子點(diǎn)選取算法實(shí)驗(yàn)的硬件配置為：CPU 處理器Intel?CoreTMi7-10875H CPU@2.30 GHz，內(nèi)存16 GB，GPU 顯卡NVIDIA GeForce RTX 2060 6 GB。

3.1 流場(chǎng)數(shù)據(jù)預(yù)處理

為驗(yàn)證本文方法在數(shù)據(jù)預(yù)處理階段所耗費(fèi)時(shí)間，分別選取不同數(shù)量的岸線點(diǎn)和格心點(diǎn)構(gòu)造新的網(wǎng)格拓?fù)潢P(guān)系，統(tǒng)計(jì)構(gòu)造新的網(wǎng)格及拓?fù)潢P(guān)系階段所耗費(fèi)的時(shí)間，如表3 所示。

表3 重構(gòu)后網(wǎng)格耗時(shí)對(duì)比Table 3 Comparison of time spent in reconstructed grid

由表3 可見，本文在數(shù)據(jù)預(yù)處理階段重構(gòu)網(wǎng)格時(shí)并不會(huì)帶來額外的計(jì)算負(fù)擔(dān)，可滿足海洋流場(chǎng)數(shù)據(jù)可視化的要求，同時(shí)在原始流場(chǎng)網(wǎng)格信息不發(fā)生更改時(shí)，新構(gòu)造的網(wǎng)格和拓?fù)潢P(guān)系可重復(fù)使用。

3.2 流場(chǎng)臨界點(diǎn)定位及分類

本文通過龐加萊指數(shù)法來提取臨界點(diǎn)，遍歷重構(gòu)的非結(jié)構(gòu)化三角網(wǎng)格并計(jì)算每個(gè)三角網(wǎng)格的龐加萊指數(shù)來判斷該網(wǎng)格是否包含臨界點(diǎn)，同時(shí)選擇Sperner完全標(biāo)號(hào)法進(jìn)行對(duì)比分析。Sperner完全標(biāo)號(hào)法通過流場(chǎng)網(wǎng)格點(diǎn)處的矢量方向變化規(guī)律來判斷網(wǎng)格內(nèi)是否包含臨界點(diǎn)，和龐加萊指數(shù)法較為類似，近年來，很多研究者結(jié)合了Sperner 的思想來提取臨界點(diǎn)。本文分別使用Sperner 完全標(biāo)號(hào)法和龐加萊指數(shù)法在非結(jié)構(gòu)化三角網(wǎng)格流場(chǎng)數(shù)據(jù)上進(jìn)行了臨界點(diǎn)提取，并借助流線對(duì)比分析2種方法所提取到的臨界點(diǎn)。

選取的實(shí)驗(yàn)范圍為4°S～26°N、99°E～126°E，該區(qū)域重構(gòu)后共包含33 992 個(gè)三角網(wǎng)格，分別使用龐加萊指數(shù)法和Sperner 完全標(biāo)號(hào)法遍歷這些三角網(wǎng)格來搜索含臨界點(diǎn)的網(wǎng)格，2 種方法所求得臨界點(diǎn)的數(shù)量及耗費(fèi)時(shí)間對(duì)比如表4 所示。該區(qū)域定位臨界點(diǎn)的位置如圖7 所示，其中，橙色點(diǎn)為Sperner 完全標(biāo)號(hào)法提取，紅色點(diǎn)為龐加萊指數(shù)法提取。

圖7 不同方法提取臨界點(diǎn)位置對(duì)比Fig.7 Comparison of critical points position extracted by different methods

表4 Sperner 完全標(biāo)號(hào)法和龐加萊指數(shù)法提取臨界點(diǎn)的結(jié)果Table 4 Results of extracting critical points by Sperner complete labeling method and Poincare index method

從表4 和圖7 可以看出，2 種方法提取的臨界點(diǎn)存在差異，因此選擇2 塊差異較大區(qū)域，并借助流線來對(duì)比分析2 種方法所提取的臨界點(diǎn)，其中，Sperner完全標(biāo)號(hào)法提取的臨界點(diǎn)位置為該方法確定包含臨界點(diǎn)的三角網(wǎng)格質(zhì)心迭代10 次后的結(jié)果，龐加萊指數(shù)法提取到的臨界點(diǎn)通過2.1.3 節(jié)中的質(zhì)心迭代法定位具體位置，閾值設(shè)置為1×10-6。2 種方法提取到的臨界點(diǎn)及周圍流線如圖8所示，其中，圖8（a）、圖8（c）為Sperner 完全標(biāo)號(hào)法提取到的臨界點(diǎn)，圖8（b）、圖8（d）為龐加萊指數(shù)法提取到的臨界點(diǎn)，圖中框選區(qū)域?yàn)? 種方法提取到的不同臨界點(diǎn)。

圖8 不同方法提取的臨界點(diǎn)對(duì)比Fig.8 Comparison of the critical points extracted by different methods

從框選區(qū)域可以看出，Sperner 完全標(biāo)號(hào)法所判定的臨界點(diǎn)會(huì)在流線方向變化較大的區(qū)域誤判，并且部分區(qū)域的臨界點(diǎn)周圍也不具有臨界點(diǎn)特征。本文也對(duì)該方法判斷存在臨界點(diǎn)的三角網(wǎng)格質(zhì)心迭代來逼近臨界點(diǎn)位置，但存在不收斂的情況；相較于Sperner 完全標(biāo)號(hào)法，雖然提取所花費(fèi)的時(shí)間稍長(zhǎng)，但龐加萊指數(shù)法不僅可以準(zhǔn)確地提取到臨界點(diǎn)，還可以通過質(zhì)心迭代精確地定位到臨界點(diǎn)位置。綜上所述，本文選擇龐加萊指數(shù)法提取非結(jié)構(gòu)化三角網(wǎng)格的臨界點(diǎn)，且若要使用Sperner 完全標(biāo)號(hào)法來準(zhǔn)確地提取臨界點(diǎn)的話，可能需要綜合其他算法，而這無疑也會(huì)增長(zhǎng)該方法所花費(fèi)的時(shí)間。

接著對(duì)龐加萊指數(shù)法提取到的臨界點(diǎn)分類，在提取臨界點(diǎn)的過程中，通過龐加萊指數(shù)可以將臨界點(diǎn)分為鞍點(diǎn)和其他類型2 類，為將臨界點(diǎn)進(jìn)一步細(xì)分，本文在非結(jié)構(gòu)化三角網(wǎng)格上構(gòu)造了臨界點(diǎn)處的雅克比矩陣并用龐加萊指數(shù)分類后的臨界點(diǎn)來驗(yàn)證本文方法。在驗(yàn)證鞍點(diǎn)時(shí)，用2.1.4 節(jié)中的方法再次構(gòu)造該鞍點(diǎn)處的雅克比矩陣，求解雅克比矩陣得到其特征值的虛部和實(shí)部，判斷該點(diǎn)的類型是否為鞍點(diǎn)。

共選擇3塊區(qū)域來驗(yàn)證本文構(gòu)造雅克比矩陣的方式對(duì)臨界點(diǎn)細(xì)分的準(zhǔn)確性，分別為區(qū)域1（4°S～26°N，99°E～126°E）、區(qū)域2（31°N～47°N，126°E～154°E）和區(qū)域3（21°S～21°N，46°E～89°E），3 塊區(qū)域在重構(gòu)后分別包含33 992、34 678、66 886 個(gè)三角網(wǎng)格，質(zhì)心迭代定位臨界點(diǎn)位置時(shí)設(shè)置的閾值為1×10-8，使用2 種方法分別對(duì)臨界點(diǎn)進(jìn)行分類，分類結(jié)果分別如表5、表6 所示。從表中可以看出，在閾值為1×10-8時(shí)，本文方法分類得到的鞍點(diǎn)類型與龐加萊指數(shù)法一致，而且可以對(duì)其他類型的臨界點(diǎn)進(jìn)一步細(xì)分，由于中心點(diǎn)在受到輕微干擾時(shí)便會(huì)轉(zhuǎn)化為2 種焦點(diǎn)，本文方法并未分類得到中心點(diǎn)，這也是需要進(jìn)一步研究的問題。

表5 不同方法對(duì)臨界點(diǎn)的分類結(jié)果（區(qū)域1）Table 5 Classification results of critical points by different methods（region 1） 單位：個(gè)

表6 不同方法對(duì)臨界點(diǎn)分類結(jié)果（區(qū)域2）Table 6 Classification results of critical points by different methods（region 2） 單位：個(gè)

筆者在分類的過程中發(fā)現(xiàn)，隨著閾值的減小，臨界點(diǎn)細(xì)分的準(zhǔn)確率會(huì)逐步提升，圖9 為不同閾值時(shí)對(duì)3 塊區(qū)域的臨界點(diǎn)分類的準(zhǔn)確率變化，可以看出，隨著閾值的減小，本文方法對(duì)臨界點(diǎn)分類的準(zhǔn)確率逐漸提升，臨界點(diǎn)細(xì)分的準(zhǔn)確率達(dá)到99%以上。由此可見，本文臨界點(diǎn)分類方法可以對(duì)非結(jié)構(gòu)化三角網(wǎng)格流場(chǎng)中的臨界點(diǎn)類型進(jìn)一步細(xì)分，且準(zhǔn)確率有一定的保證，可為非結(jié)構(gòu)化三角網(wǎng)格流場(chǎng)的拓?fù)浞治龅於己玫幕A(chǔ)。

圖9 臨界點(diǎn)分類準(zhǔn)確率與閾值間的關(guān)系Fig.9 The relationship between classification accuracy of critical point and threshold

3.3 基于網(wǎng)格密度和“最大得分”法的種子點(diǎn)選取算法

流場(chǎng)中網(wǎng)格密度由某一網(wǎng)格頂點(diǎn)距離其相鄰頂點(diǎn)的平均距離來表示，平均距離越小，即網(wǎng)格密度越大，賦予的權(quán)重也越大。由此得到了流場(chǎng)網(wǎng)格密度的標(biāo)量場(chǎng)，圖10 分別為原始網(wǎng)格和網(wǎng)格密度的標(biāo)量場(chǎng)，網(wǎng)格密度區(qū)間n=28，可以看出網(wǎng)格密度與原始網(wǎng)格的疏密相匹配。

圖10 三角網(wǎng)格密度劃分Fig.10 Density division of triangular grid

本文還通過選擇重構(gòu)后三角網(wǎng)格中的某一頂點(diǎn)及其雙層鄰接頂點(diǎn)，比較它們矢量值的平方范數(shù)的標(biāo)量值大小來計(jì)算該頂點(diǎn)的評(píng)分，由此計(jì)算得到了流場(chǎng)的“最大得分”標(biāo)量場(chǎng)，如圖11 所示。最后根據(jù)網(wǎng)格頂點(diǎn)處的密度來設(shè)置網(wǎng)格頂點(diǎn)的“最大得分”閾值，若該點(diǎn)得分大于閾值，則該網(wǎng)格頂點(diǎn)附近可能存在流場(chǎng)特征，將該點(diǎn)選為種子點(diǎn)。

圖11 “最大得分”標(biāo)量場(chǎng)Fig.11 “maximum score”scalar field

本文基于“最大得分”法生成的標(biāo)量場(chǎng)和網(wǎng)格密度選擇種子點(diǎn)，其中，tmax設(shè)置為0.95，tmin設(shè)置為0.4，密度區(qū)間n=28，并且使用隨機(jī)種子點(diǎn)播種、均勻種子點(diǎn)播種、基于信息熵的種子點(diǎn)播種和基于特征的種子點(diǎn)放置方法在該區(qū)域放置種子點(diǎn)，使用歐拉法生成流線觀察各種子點(diǎn)選取方法對(duì)流場(chǎng)的可視化效果。不同種子點(diǎn)放置方法對(duì)同一區(qū)域可視化效果如圖12 所示，其中紅色點(diǎn)為2.1 節(jié)中的方法定位到的臨界點(diǎn)，用于驗(yàn)證不同種子點(diǎn)算法對(duì)流場(chǎng)關(guān)鍵特征區(qū)域的可視化效果。本文評(píng)價(jià)種子點(diǎn)放置方法主要有兩個(gè)標(biāo)準(zhǔn)：一個(gè)是對(duì)流場(chǎng)關(guān)鍵特征區(qū)域的可視化效果，即能否很好地表現(xiàn)出臨界點(diǎn)及其附近區(qū)域的流場(chǎng)特征；另一個(gè)是對(duì)整個(gè)流場(chǎng)區(qū)域的覆蓋效果，即能否表現(xiàn)出流場(chǎng)的全局信息。

圖12 不同種子點(diǎn)放置方法生成的流線Fig.12 Streamlines generated by different seed placement methods

圖12（a）為該區(qū)域的原始網(wǎng)格信息，可以看到該區(qū)域左側(cè)網(wǎng)格稀疏，右側(cè)網(wǎng)格密集，可以很好地代表非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)。圖12（b）為在該區(qū)域隨機(jī)播種136 個(gè)種子點(diǎn)生成的流線效果，從框選區(qū)域可以看出，該方法無法表現(xiàn)出流場(chǎng)關(guān)鍵特征，且流場(chǎng)全局信息可視化效果較差。圖12（c）為均勻播種106 個(gè)種子點(diǎn)所生成的流線，該方法生成的流線可以有效覆蓋流場(chǎng)，表現(xiàn)出流場(chǎng)的全局信息，但流場(chǎng)關(guān)鍵特征區(qū)域可視化效果較差。圖12（d）為在該區(qū)域信息熵較大的區(qū)域選取109 個(gè)種子點(diǎn)生成的流線，可以看到該方法對(duì)網(wǎng)格較密集的區(qū)域可視化效果較好，但當(dāng)網(wǎng)格過于稀疏時(shí)，容易因格點(diǎn)過少，導(dǎo)致計(jì)算得到的信息熵?zé)o法體現(xiàn)出該區(qū)域的特征信息，從而無法選擇種子點(diǎn)生成流線，致使區(qū)域左側(cè)的流場(chǎng)關(guān)鍵特征區(qū)域可視化效果很差。圖12（e）為基于特征放置145 個(gè)種子點(diǎn)所生成的流線，該方法可以很好地展現(xiàn)流場(chǎng)的關(guān)鍵特征，即臨界點(diǎn)周圍的區(qū)域，但卻專注于流場(chǎng)特征，不能有效地表達(dá)出流場(chǎng)的全局信息；圖12（f）為本文方法選取96 個(gè)種子點(diǎn)所生成的流線，可以看出，本文方法相較于上述方法播種的種子點(diǎn)較少，但生成的流線卻可以有效地表征流場(chǎng)關(guān)鍵特征，而且生成的流線也更均勻，可以覆蓋整個(gè)流場(chǎng)區(qū)域，整體可視化效果更好。

4 結(jié)束語(yǔ)

本文基于非結(jié)構(gòu)化三角網(wǎng)格流場(chǎng)數(shù)據(jù)，設(shè)計(jì)了多個(gè)算法，解決了缺乏網(wǎng)格拓?fù)潢P(guān)系、三角網(wǎng)格內(nèi)的臨界點(diǎn)難以準(zhǔn)確定位與細(xì)分、現(xiàn)有種子點(diǎn)算法不適用于非結(jié)構(gòu)化網(wǎng)格等問題。在數(shù)據(jù)預(yù)處理方面，通過Delaunay 三角剖分重構(gòu)了三角網(wǎng)格格心的拓?fù)潢P(guān)系，直接用格心矢量數(shù)據(jù)來提取流場(chǎng)特征，簡(jiǎn)化了非結(jié)構(gòu)化網(wǎng)格上的點(diǎn)定位和相鄰頂點(diǎn)的查找；在臨界點(diǎn)定位及分類方面，基于重構(gòu)的流場(chǎng)網(wǎng)格，用龐加萊指數(shù)法確定存在臨界點(diǎn)的三角網(wǎng)格并構(gòu)造了質(zhì)心迭代法準(zhǔn)確定位臨界點(diǎn)位置，與Sperner 完全標(biāo)號(hào)法相比，本文選用的龐加萊指數(shù)法可以更加準(zhǔn)確地判斷存在臨界點(diǎn)的三角網(wǎng)格；在提升流場(chǎng)可視化效果方面，為定量分析非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)，改進(jìn)了“最大得分”法使其適用于非結(jié)構(gòu)化網(wǎng)格，得到原始流場(chǎng)的“最大得分”標(biāo)量場(chǎng)，并根據(jù)網(wǎng)格密度動(dòng)態(tài)設(shè)置閾值來選取種子點(diǎn)，通過對(duì)比多種種子點(diǎn)放置方法，證明本文種子點(diǎn)放置方法可以在兼顧流場(chǎng)關(guān)鍵特征的同時(shí)，更好地表現(xiàn)流場(chǎng)全局信息，提升流場(chǎng)的可視化效果。

本文算法僅適用于二維非結(jié)構(gòu)化網(wǎng)格流場(chǎng)數(shù)據(jù)，需要繼續(xù)研究從而擴(kuò)展至三維流場(chǎng)數(shù)據(jù)；同時(shí)在提取流場(chǎng)特征的過程中，未來可以借助當(dāng)前熱門的人工智能技術(shù)對(duì)流場(chǎng)特征進(jìn)行更深層次的分析，挖掘流場(chǎng)內(nèi)數(shù)據(jù)間的聯(lián)系，進(jìn)一步揭示海洋流場(chǎng)運(yùn)動(dòng)規(guī)律。