季民，任靜*，張立國，李婷，孫勇

( 1. 山東科技大學(xué) 測繪與空間信息學(xué)院，山東 青島 266590；2. 山東省國土測繪院，山東 濟南 250102)

1 引言

隨著海洋探測技術(shù)的快速發(fā)展，人類取得了前所未有的海洋大數(shù)據(jù)集合。海洋流場作為典型的矢量場，其特征結(jié)構(gòu)復(fù)雜而多變，通過對其拓撲特征結(jié)構(gòu)進行可視化變化分析，對于理解諸多海洋現(xiàn)象的產(chǎn)生、發(fā)展和演變等都具有十分重要的意義。

矢量場的拓撲特征可視化，自20 世紀90 年代提出以來就得到了迅猛發(fā)展。Helman 和Hesselink[1-2]提出了基于雅克比矩陣特征值的臨界點分類和矢量場的拓撲分析方法，而臨界點作為矢量場拓撲結(jié)構(gòu)分析的基礎(chǔ)和重要組成部分，國內(nèi)外學(xué)者針對矢量場臨界點的提取進行了一系列研究。鑒于臨界點一般位于0 值矢量處，Lavin 等[3]通過線性插值的方式進行了0 矢量搜尋，Batra 和Hesselink[4]通過對三角網(wǎng)格上的分段線性矢量場進行線性插值來提取臨界點，Bhatia 等[5]及王文濤[6]通過判斷單純形頂點處的矢量值是否構(gòu)成一個內(nèi)部包含原點的單純形來快速判定臨界點的存在。吳曉莉等[7]則將Sperner 引理引入到臨界點檢測中，進行臨界點的檢測。以上研究均采用的是局部檢測法，另有學(xué)者基于全局視角進行了臨界點提取，Polthier 和Preuβ[8-9]提出了一種利用離散Hodge分解的2D 非結(jié)構(gòu)化三角形網(wǎng)格矢量場臨界點提取方法，Chen 等[10]基于Morse 分解的方法，通過流組合化、強聯(lián)通分量提取、商圖簡化等步驟進行了臨界點區(qū)域的識別。而海洋流場作為典型的矢量場，眾多學(xué)者也對其進行了一系列拓撲特征提取與可視化研究，管倩倩[11]基于拓撲理論實現(xiàn)了海洋水團特征提取，廖忠云[12]從提取特征點出發(fā)，應(yīng)用歐拉數(shù)值積分算法進行了海洋流線追蹤與可視化，王輝贊等[13]通過提取渦旋中心和大小來進行渦旋軌跡追蹤與表達，牛嬋等[14]通過提取臨界點來反映海洋流場空間拓撲結(jié)構(gòu)。

在前述矢量場拓撲分析法中的臨界點理論、插值求解臨界點算法和Sperner 引理檢測臨界點算法的啟發(fā)下，綜合雙線性插值和Sperner 引理檢測兩種算法，通過解決網(wǎng)格插值的二義性和零值網(wǎng)格點的影響問題，實現(xiàn)了海洋流場臨界點的提取和分類。研究結(jié)果表明，本文算法提取的臨界點更加全面、合理。

2 流場臨界點提取方法原理

2.1 臨界點理論

一個矢量場的拓撲結(jié)構(gòu)由臨界點和連接臨界點的積分曲線或曲面組成[15]。依據(jù)Helman 和Hesselink提出的臨界點理論，臨界點又被稱為奇點、不動點或平衡點，在二維流場中是指矢量場的兩個分量均為0 的點，對于非退化的臨界點（x，y），可以用臨界點位置處的偏導(dǎo)數(shù)矩陣（即雅克比矩陣）來表征矢量場及其附近曲線的行為[2-3]，其公式如下：

對于二維流場，可根據(jù)雅克比矩陣兩個特征值λ1和λ2的實部Re和虛部Im的正負等情況，來進行臨界點分類，具體分類情況如圖1 所示（圖中R表示實部值，I表示虛部值）。據(jù)此臨界點主要分為交點、聚點、馬鞍點和中心點，而交點和聚點又可進一步分為吸引交點、排斥交點、吸引聚點、排斥聚點。

圖1 二維流場臨界點分類圖Fig. 1 Classification of critical points in two-dimensional flow field

2.2 Sperner 引理

現(xiàn)有的矢量場臨界點檢測方法如MC（Marching Cube）方法[16]、基于幾何代數(shù)法[17]、基于More 分解方法[10]、基于物理特征的方法[18]及龐加萊指數(shù)法[19-20]等，這些方法各有優(yōu)缺點，而吳曉莉等[7]將Sperner 引理引入到臨界點檢測中，給出了一種臨界點檢測的新算法，該算法定義速度分量均為0 的點為臨界點，借助Sperner 引理和完全標號法進行了臨界點提取算法研究。

Sperner 引理[21]定義為：給定一個大三角形V1V2V3，并將它三角化（把它劃分成有限多個較細的三角形，且每個細三角形的邊都是另一個細三角形的邊或者落在大三角形的邊上），若將各頂點以下述的規(guī)定標記：

（1）頂點Vi的標號為i，i=1, 2, 3；

（2）在ViVj邊上的頂點只可以用i或者j作為標號；

（3）不在大三角形邊上的頂點可以隨意以1，2，3 作為標號。

則至少存在一個細三角形其3 個頂點的標號分別為1，2，3[22]。

Sperner 完全標號：給定二維流場內(nèi)的2×2 網(wǎng)格，各頂點的標號為中心點處速度矢量V（u，v）落在u，v為坐標軸的象限代碼。當出現(xiàn)速度矢量與坐標軸共線的情形時，規(guī)定速度矢量與u軸正向重合的頂點標號為1；速度矢量與v軸正向重合的頂點標號為2；速度矢量與u軸負向重合的頂點標號為3；速度矢量與v軸負向重合的頂點標號為4。若4 個網(wǎng)格的標號分別為1，2，3，4 時，則稱其是按矢量方向Sperner 完全標號的，簡稱Sperner 完全標號[7]，具體完全標號形式如圖2 所示。

圖2 二維流場完全標號示意圖[7]Fig. 2 Fully numbered illustration of a two-dimensional flow field[7]

3 基于雙線性插值的臨界點提取

海洋流場中的臨界點與流場中有意義的形狀、結(jié)構(gòu)、變化和現(xiàn)象（如渦流、激波等）密切相關(guān)，因此，臨界點提取對于海洋流場拓撲結(jié)構(gòu)特征的研究具有重要意義。目前的臨界點檢測方法各有優(yōu)缺點，而線性插值方法簡潔明了，易于程序?qū)崿F(xiàn)，為此本文基于雙線性插值方法進行了二維海洋流場中臨界點的提取，并基于聚合思想解決了網(wǎng)格插值中的二義性問題，具體算法過程包含了：候選網(wǎng)格的篩選、臨界點提取及分類等過程。

3.1 臨界點候選網(wǎng)格的篩選

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

本文研究所用的海洋流場數(shù)據(jù)來源于美國國家海洋和大氣管理局（NOAA）的國家環(huán)境信息中心（NCEI），為NetCDF 格式的海洋再分析數(shù)據(jù)，包括了經(jīng)向流速u和緯向流速v，數(shù)據(jù)空間范圍為20°～42°N，98°～116°W，時間為2019 年8 月4 日，深度為5 000 m，分辨率為0.033°。其中數(shù)據(jù)中的正負代表了海水流動的不同方向，東向和北向為正。

3.1.2 基于滑動窗口的候選網(wǎng)格提取

一般線性流場中臨界點位于兩個異向網(wǎng)格的中間，為了便于在經(jīng)、緯兩個方向同時篩選包含臨界點的網(wǎng)格，為此，本文選擇2×2 滑動窗口，分別在獲取的u、v方向的兩個數(shù)據(jù)層上進行遍歷，篩選出臨近網(wǎng)格流向異號的單元作為臨界點的候選網(wǎng)格單元，篩選出的網(wǎng)格值正負異號的情況具體如圖3 所示，主要分為一正三負、兩正兩負、三正一負等4 種情形，在圖3a，圖3b，圖3d 的3 種情形中，可根據(jù)網(wǎng)格值的正負情形直接通過線性插值獲得臨界點等值線的位置，具體位置如圖中加粗實線所示，而圖3c 的情形，由于存在正負值交叉的情形，因而帶來了線性插值中的二義性問題。

圖3 候選網(wǎng)格情況Fig. 3 Candidate grid case

3.1.3 二義性候選網(wǎng)格的處理

針對圖3c 情形中的臨界點插值，此時可能存在如圖4a 和圖4b 所示的兩種可能連接方式，為了進一步確定臨界點等值線的連接方向，本文基于降低分辨率的聚合思想，將4 個二義性候選網(wǎng)格的流速均值A(chǔ)vg 作為降低分辨率后網(wǎng)格中心點處的值。Avg 存在3 種情形，若Avg=0，則直接判定該中心點為臨界點；若Avg＞0，則臨界點等值線的連接方向與負向?qū)蔷W(wǎng)格點的連線方向一致，如圖4a；若Avg＜0，則臨界點等值線的連接方向與正向?qū)蔷W(wǎng)格點的連線方向一致，如圖4b。

圖4 二義性網(wǎng)格Fig. 4 Ambiguous grid

3.1.4 0 值候選網(wǎng)格單元的處理

上述候選網(wǎng)格的篩選與二義性處理等過程，只考慮了網(wǎng)格值為正負的情形，當2×2 的滑動窗口中存在0 值情形時，若不加以考慮，就會造成臨界點的遺漏。為此，根據(jù)0 值網(wǎng)格的存在情況，本文將臨界點的插值分為了如圖5 所示的9 種情形，并且分別給出了臨界點等值線的連接方式（加粗實線所示）。其中，第9 種情形比較特殊，即4 個網(wǎng)格的值均為0，為了解決這個問題，本文借助聚合思想，通過降低分辨率，取臨近4 個網(wǎng)格的速度均值為降低分辨率后網(wǎng)格的速度分量，并以當前降低分辨率后的大網(wǎng)格為左上角的起始網(wǎng)格，依然按照2×2 的窗口進行滑動篩選，其結(jié)果必然為圖5 的9 種情形之一，若還存在第9 種情形，則進行迭代聚合。

圖5 網(wǎng)格值含0 的情況Fig. 5 Grid value with 0

3.2 基于雙線性插值的臨界點位置計算

當網(wǎng)格單元足夠小時，可以認為沿著網(wǎng)格的邊數(shù)據(jù)場是連續(xù)線性變化的[23]，因此可以利用雙線性插值的方法在網(wǎng)格邊上求出0 值P、Q點的位置，如圖6a和圖6b 的計算結(jié)果。分別連接u、v向計算結(jié)果中的插值點，兩條連線的交點即為u、v向均為0 值的臨界點，如圖6c 中的M點。

圖6 雙線性插值過程圖Fig. 6 Bilinear interpolation process chart

3.3 臨界點分類

利用本文的雙線性插值算法，對美國東部某海域5 000 m 深度0.033°分辨率的再分析流場數(shù)據(jù)進行了臨界點提取，結(jié)果如圖7a 所示，為了進一步確定臨界點的分類，通過求其位置處u、v速度分量的偏導(dǎo)數(shù)，進行雅克比矩陣構(gòu)建，然后依據(jù)圖1 所示的分類方法對臨界點進行了分類，具體分類結(jié)果如圖7b 所示，共22 個點。

圖7 基于雙線性插值的臨界點提?。╝）與分類（b）結(jié)果Fig. 7 Extraction (a) and classification (b) of critical points based on bilinear interpolation

本文的臨界點雙線性插值與分類算法，簡單明了，算法的時間復(fù)雜度為O（n2），空間復(fù)雜度為O（nlog2n），提取結(jié)果全面，且臨界點位置計算精確，可適用于多數(shù)二維流場的特征點提取。

4 基于Sperner 完全標號的臨界點提取

前文雙線性插值的臨界點提取算法是基于臨界點速度分量均為0 的假設(shè)，但由于所采用的海洋流場再分析數(shù)據(jù)中可能存在計算舍入誤差造成的近0 值速度向量，并且在實際的流場物理場景中，也存在速度非0 的臨界點。由于Sperner 完全標號理論依據(jù)流向變化進行臨界點網(wǎng)格的篩選，相比于雙線性插值方法依據(jù)u、v流速數(shù)據(jù)求解，計算簡便，為此本文也基于Sperner 完全標號理論進行了臨界點的提取算法研究。其算法過程包括了：網(wǎng)格流向計算、臨界點候選網(wǎng)格篩選、臨界點識別等。

4.1 臨界點提取與分析

該算法的執(zhí)行主要包括以下幾個步驟：

第一步，基于海洋流場再分析數(shù)據(jù)中的u、v速度分量，進行網(wǎng)格點速度向量計算；

第二步，進行臨界點候選網(wǎng)格篩選，首先以2×2的滑動窗口對整個流場進行遍歷，在滑動窗口遍歷過程中，依據(jù)Sperner 完全標號理論，對4 網(wǎng)格進行標號處理，若4 個網(wǎng)格完全標號，則認為4 個網(wǎng)格內(nèi)存在臨界點，從而將其作為臨界點候選網(wǎng)格予以保留，否則繼續(xù)進行窗口滑動；

第三步，進行臨界點提取，臨界點候選網(wǎng)格篩選完成后，多數(shù)研究是通過插值算法進行臨界點提取，為了具有更寬泛的適用性，本文在此提出了基于最小值法的臨界點提取規(guī)則，根據(jù)該規(guī)則，將4 個候選網(wǎng)格中速度向量模最小的網(wǎng)格中心作為臨界點，從而完成臨界點的提取。

為了檢測該算法提取的臨界點的位置精度，本文還基于線性插值算法對上述過程提取的臨界點進行了位置計算，具體計算結(jié)果如圖8 所示，其中紅色圓點表示插值位置，綠色圓點表示最小值網(wǎng)格的中心點。通過計算分析發(fā)現(xiàn)，插值圓點基本上位于最小值網(wǎng)格內(nèi)，且網(wǎng)格值越小，插值點距離中心越近，當網(wǎng)格分辨率足夠小時，兩類點的距離可忽略不計，但是最小值法臨界點的提取，可避免插值處理過程，其算法復(fù)雜度為O（n），比插值算法更加高效，并且可忽略插值算法中的0 值影響，因而，算法的適用性更廣。

圖8 插值與最小值法提取結(jié)果對比Fig. 8 Comparison between the results of interpolation and minimum method

4.2 臨界點分類

依據(jù)Sperner 完全標號最小值法的臨界點提取過程，同樣對美國東部某海域5 000 m 深度0.033°分辨率的再分析流場數(shù)據(jù)進行了臨界點提取，結(jié)果如圖9a所示，同樣基于雅克比矩陣計算，對提取的臨界點進行了分類，具體結(jié)果如圖9b 所示，共獲得11 個臨界點。

通過對比圖7 和圖9 的提取結(jié)果，發(fā)現(xiàn)Sperner 完全標號法雖然未能檢索出更多的臨界點，但的確另外提取了幾個插值算法遺漏的臨界點，如圖9 中三角符號所示的臨界點。該算法計算效率高、速度快，可較好的彌補雙線性插值遺漏的臨界點。

圖9 基于Sperner 完全標號的臨界點提取（a）與分類（b）結(jié)果Fig. 9 Extraction (a) and classification (b) of critical points based on Sperner fully labeled

5 兩種算法綜合的臨界點提取與實驗分析

從兩種算法的臨界點提取結(jié)果圖看，雙線性插值算法提取結(jié)果數(shù)量更多，覆蓋范圍更廣，但Sperner 完全標號最小值法提取了雙線性插值法未提取的臨界點，并且速度更快，因此可綜合兩種算法臨界點提取結(jié)果作為最終臨界點提取結(jié)果。

5.1 算法綜合處理

統(tǒng)計分析兩種算法的分類結(jié)果，如表1 所示。

表1 兩種方法分類結(jié)果統(tǒng)計表Table 1 Statistical table of classification results of two methods

以上兩種臨界點提取算法各有優(yōu)缺點，通過將兩種算法提取結(jié)果的合并、去重等處理，可得到如圖10所示的更加全面的臨界點提取與分類結(jié)果，其中，中心點7 個，交點7 個，吸引聚點5 個，排斥聚點5 個，共24 個。

圖10 綜合提取分類結(jié)果Fig. 10 Comprehensive extraction and classification result chart

5.2 算法驗證

為了進一步驗證本文算法的適用性及可行性，針對不同區(qū)域、不同深度的大量海洋流場數(shù)據(jù)進行多次驗證分析，其中表2 所示，是以美國沿海區(qū)域深度5 000 m、大西洋部分海域深度2 500 m 和太平洋海域深度3 000 m 的數(shù)據(jù)為例進行的對比分析。提取結(jié)果顯示：雙線性插值算法臨界點提取結(jié)果較為全面完整，但對非0 值臨界點的提取會偶有遺漏，而Sperner完全標號法計算高效，較好地彌補了雙線性插值算法遺漏的部分臨界點，如表2 中美國沿海深度5 000 m數(shù)據(jù)中，彌補提取了1 個中心點、1 個交點和1 個排斥聚點；對大西洋部分海域深度2 500 m 數(shù)據(jù)中，彌補提取了1 個吸引聚點、1 個排斥聚點；對太平洋海域深度3 000 m 數(shù)據(jù)中，彌補提取了1 個吸引聚點。研究表明通過兩種算法的綜合處理，全面完整的提取了二維流場中的拓撲結(jié)構(gòu)臨界點，且算法簡潔明了，可為二維流場的拓撲結(jié)構(gòu)分析提供一種新的研究思路。

表2 不同數(shù)據(jù)驗證結(jié)果Table 2 Different data validation results

6 結(jié)論

本文依據(jù)矢量場拓撲結(jié)構(gòu)分析的臨界點理論和Sperner 引理，綜合改進后的雙線性插值算法和Sperner 完全標號法，實現(xiàn)了海洋流場數(shù)據(jù)的臨界點特征提取。相比于現(xiàn)有算法，本文算法臨界點提取效果理想，且易于程序?qū)崿F(xiàn)。改進后的雙線性插值算法添加滑動窗口處理，便于臨界點候選單元的篩選，涉及的聚合思想為解決插值網(wǎng)格二義性問題提供了新思路，而0 值網(wǎng)格中臨界點等值線連接方式的考慮，補充了流場網(wǎng)格插值情形，使得提取結(jié)果更加全面，且臨界點位置計算精確，可適用于多數(shù)二維流場的特征點提取。同時，基于Sperner 完全標號的臨界點提取算法拋開傳統(tǒng)插值求解臨界點思想，創(chuàng)造性地提出了最小值法臨界點提取規(guī)則，避免了插值處理過程，并且可忽略插值算法中的0 值影響，效率更高，有著更廣的適用性。而兩種算法作為互補，算法的綜合可實現(xiàn)更加全面的臨界點提取與分類結(jié)果，為進一步實現(xiàn)海洋流場拓撲分析的流線追蹤奠定臨界點基礎(chǔ)。