面向拓?fù)浞治龅暮Ｑ罅鲌雠R界點(diǎn)提取算法研究

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

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

1 引言

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

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

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

2 流場臨界點(diǎn)提取方法原理

2.1 臨界點(diǎn)理論

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

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

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

2.2 Sperner 引理

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

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

（1）頂點(diǎn)Vi的標(biāo)號為i，i=1, 2, 3；

（2）在ViVj邊上的頂點(diǎn)只可以用i或者j作為標(biāo)號；

（3）不在大三角形邊上的頂點(diǎn)可以隨意以1，2，3 作為標(biāo)號。

則至少存在一個(gè)細(xì)三角形其3 個(gè)頂點(diǎn)的標(biāo)號分別為1，2，3[22]。

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

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

3 基于雙線性插值的臨界點(diǎn)提取

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

3.1 臨界點(diǎn)候選網(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，時(shí)間為2019 年8 月4 日，深度為5 000 m，分辨率為0.033°。其中數(shù)據(jù)中的正負(fù)代表了海水流動(dòng)的不同方向，東向和北向?yàn)檎?/p>

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

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

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

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

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

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

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

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

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

3.2 基于雙線性插值的臨界點(diǎn)位置計(jì)算

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

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

3.3 臨界點(diǎn)分類

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

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

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

4 基于Sperner 完全標(biāo)號的臨界點(diǎn)提取

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

4.1 臨界點(diǎn)提取與分析

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

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

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

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

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

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

4.2 臨界點(diǎn)分類

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

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

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

5 兩種算法綜合的臨界點(diǎn)提取與實(shí)驗(yàn)分析

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

5.1 算法綜合處理

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

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

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

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

5.2 算法驗(yàn)證

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

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

6 結(jié)論

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