面向拓扑分析的海洋流场临界点提取算法研究

季民，任静*，张立国，李婷，孙勇

( 1. 山东科技大学 测绘与空间信息学院，山东 青岛 266590；2. 山东省国土测绘院，山东 济南 250102)

1 引言

随着海洋探测技术的快速发展，人类取得了前所未有的海洋大数据集合。海洋流场作为典型的矢量场，其特征结构复杂而多变，通过对其拓扑特征结构进行可视化变化分析，对于理解诸多海洋现象的产生、发展和演变等都具有十分重要的意义。

矢量场的拓扑特征可视化，自20 世纪90 年代提出以来就得到了迅猛发展。Helman 和Hesselink[1-2]提出了基于雅克比矩阵特征值的临界点分类和矢量场的拓扑分析方法，而临界点作为矢量场拓扑结构分析的基础和重要组成部分，国内外学者针对矢量场临界点的提取进行了一系列研究。鉴于临界点一般位于0 值矢量处，Lavin 等[3]通过线性插值的方式进行了0 矢量搜寻，Batra 和Hesselink[4]通过对三角网格上的分段线性矢量场进行线性插值来提取临界点，Bhatia 等[5]及王文涛[6]通过判断单纯形顶点处的矢量值是否构成一个内部包含原点的单纯形来快速判定临界点的存在。吴晓莉等[7]则将Sperner 引理引入到临界点检测中，进行临界点的检测。以上研究均采用的是局部检测法，另有学者基于全局视角进行了临界点提取，Polthier 和Preuβ[8-9]提出了一种利用离散Hodge分解的2D 非结构化三角形网格矢量场临界点提取方法，Chen 等[10]基于Morse 分解的方法，通过流组合化、强联通分量提取、商图简化等步骤进行了临界点区域的识别。而海洋流场作为典型的矢量场，众多学者也对其进行了一系列拓扑特征提取与可视化研究，管倩倩[11]基于拓扑理论实现了海洋水团特征提取，廖忠云[12]从提取特征点出发，应用欧拉数值积分算法进行了海洋流线追踪与可视化，王辉赞等[13]通过提取涡旋中心和大小来进行涡旋轨迹追踪与表达，牛婵等[14]通过提取临界点来反映海洋流场空间拓扑结构。

在前述矢量场拓扑分析法中的临界点理论、插值求解临界点算法和Sperner 引理检测临界点算法的启发下，综合双线性插值和Sperner 引理检测两种算法，通过解决网格插值的二义性和零值网格点的影响问题，实现了海洋流场临界点的提取和分类。研究结果表明，本文算法提取的临界点更加全面、合理。

2 流场临界点提取方法原理

2.1 临界点理论

一个矢量场的拓扑结构由临界点和连接临界点的积分曲线或曲面组成[15]。依据Helman 和Hesselink提出的临界点理论，临界点又被称为奇点、不动点或平衡点，在二维流场中是指矢量场的两个分量均为0 的点，对于非退化的临界点（x，y），可以用临界点位置处的偏导数矩阵（即雅克比矩阵）来表征矢量场及其附近曲线的行为[2-3]，其公式如下：

对于二维流场，可根据雅克比矩阵两个特征值λ1和λ2的实部Re和虚部Im的正负等情况，来进行临界点分类，具体分类情况如图1 所示（图中R表示实部值，I表示虚部值）。据此临界点主要分为交点、聚点、马鞍点和中心点，而交点和聚点又可进一步分为吸引交点、排斥交点、吸引聚点、排斥聚点。

图1 二维流场临界点分类图Fig. 1 Classification of critical points in two-dimensional flow field

2.2 Sperner 引理

现有的矢量场临界点检测方法如MC（Marching Cube）方法[16]、基于几何代数法[17]、基于More 分解方法[10]、基于物理特征的方法[18]及庞加莱指数法[19-20]等，这些方法各有优缺点，而吴晓莉等[7]将Sperner 引理引入到临界点检测中，给出了一种临界点检测的新算法，该算法定义速度分量均为0 的点为临界点，借助Sperner 引理和完全标号法进行了临界点提取算法研究。

Sperner 引理[21]定义为：给定一个大三角形V1V2V3，并将它三角化（把它划分成有限多个较细的三角形，且每个细三角形的边都是另一个细三角形的边或者落在大三角形的边上），若将各顶点以下述的规定标记：

（1）顶点Vi的标号为i，i=1, 2, 3；

（2）在ViVj边上的顶点只可以用i或者j作为标号；

（3）不在大三角形边上的顶点可以随意以1，2，3 作为标号。

则至少存在一个细三角形其3 个顶点的标号分别为1，2，3[22]。

Sperner 完全标号：给定二维流场内的2×2 网格，各顶点的标号为中心点处速度矢量V（u，v）落在u，v为坐标轴的象限代码。当出现速度矢量与坐标轴共线的情形时，规定速度矢量与u轴正向重合的顶点标号为1；速度矢量与v轴正向重合的顶点标号为2；速度矢量与u轴负向重合的顶点标号为3；速度矢量与v轴负向重合的顶点标号为4。若4 个网格的标号分别为1，2，3，4 时，则称其是按矢量方向Sperner 完全标号的，简称Sperner 完全标号[7]，具体完全标号形式如图2 所示。

图2 二维流场完全标号示意图[7]Fig. 2 Fully numbered illustration of a two-dimensional flow field[7]

3 基于双线性插值的临界点提取

海洋流场中的临界点与流场中有意义的形状、结构、变化和现象（如涡流、激波等）密切相关，因此，临界点提取对于海洋流场拓扑结构特征的研究具有重要意义。目前的临界点检测方法各有优缺点，而线性插值方法简洁明了，易于程序实现，为此本文基于双线性插值方法进行了二维海洋流场中临界点的提取，并基于聚合思想解决了网格插值中的二义性问题，具体算法过程包含了：候选网格的筛选、临界点提取及分类等过程。

3.1 临界点候选网格的筛选

3.1.1 数据来源

本文研究所用的海洋流场数据来源于美国国家海洋和大气管理局（NOAA）的国家环境信息中心（NCEI），为NetCDF 格式的海洋再分析数据，包括了经向流速u和纬向流速v，数据空间范围为20°～42°N，98°～116°W，时间为2019 年8 月4 日，深度为5 000 m，分辨率为0.033°。其中数据中的正负代表了海水流动的不同方向，东向和北向为正。

3.1.2 基于滑动窗口的候选网格提取

一般线性流场中临界点位于两个异向网格的中间，为了便于在经、纬两个方向同时筛选包含临界点的网格，为此，本文选择2×2 滑动窗口，分别在获取的u、v方向的两个数据层上进行遍历，筛选出临近网格流向异号的单元作为临界点的候选网格单元，筛选出的网格值正负异号的情况具体如图3 所示，主要分为一正三负、两正两负、三正一负等4 种情形，在图3a，图3b，图3d 的3 种情形中，可根据网格值的正负情形直接通过线性插值获得临界点等值线的位置，具体位置如图中加粗实线所示，而图3c 的情形，由于存在正负值交叉的情形，因而带来了线性插值中的二义性问题。

图3 候选网格情况Fig. 3 Candidate grid case

3.1.3 二义性候选网格的处理

针对图3c 情形中的临界点插值，此时可能存在如图4a 和图4b 所示的两种可能连接方式，为了进一步确定临界点等值线的连接方向，本文基于降低分辨率的聚合思想，将4 个二义性候选网格的流速均值Avg 作为降低分辨率后网格中心点处的值。Avg 存在3 种情形，若Avg=0，则直接判定该中心点为临界点；若Avg＞0，则临界点等值线的连接方向与负向对角网格点的连线方向一致，如图4a；若Avg＜0，则临界点等值线的连接方向与正向对角网格点的连线方向一致，如图4b。

图4 二义性网格Fig. 4 Ambiguous grid

3.1.4 0 值候选网格单元的处理

上述候选网格的筛选与二义性处理等过程，只考虑了网格值为正负的情形，当2×2 的滑动窗口中存在0 值情形时，若不加以考虑，就会造成临界点的遗漏。为此，根据0 值网格的存在情况，本文将临界点的插值分为了如图5 所示的9 种情形，并且分别给出了临界点等值线的连接方式（加粗实线所示）。其中，第9 种情形比较特殊，即4 个网格的值均为0，为了解决这个问题，本文借助聚合思想，通过降低分辨率，取临近4 个网格的速度均值为降低分辨率后网格的速度分量，并以当前降低分辨率后的大网格为左上角的起始网格，依然按照2×2 的窗口进行滑动筛选，其结果必然为图5 的9 种情形之一，若还存在第9 种情形，则进行迭代聚合。

图5 网格值含0 的情况Fig. 5 Grid value with 0

3.2 基于双线性插值的临界点位置计算

当网格单元足够小时，可以认为沿着网格的边数据场是连续线性变化的[23]，因此可以利用双线性插值的方法在网格边上求出0 值P、Q点的位置，如图6a和图6b 的计算结果。分别连接u、v向计算结果中的插值点，两条连线的交点即为u、v向均为0 值的临界点，如图6c 中的M点。

图6 双线性插值过程图Fig. 6 Bilinear interpolation process chart

3.3 临界点分类

利用本文的双线性插值算法，对美国东部某海域5 000 m 深度0.033°分辨率的再分析流场数据进行了临界点提取，结果如图7a 所示，为了进一步确定临界点的分类，通过求其位置处u、v速度分量的偏导数，进行雅克比矩阵构建，然后依据图1 所示的分类方法对临界点进行了分类，具体分类结果如图7b 所示，共22 个点。

图7 基于双线性插值的临界点提取（a）与分类（b）结果Fig. 7 Extraction (a) and classification (b) of critical points based on bilinear interpolation

本文的临界点双线性插值与分类算法，简单明了，算法的时间复杂度为O（n2），空间复杂度为O（nlog2n），提取结果全面，且临界点位置计算精确，可适用于多数二维流场的特征点提取。

4 基于Sperner 完全标号的临界点提取

前文双线性插值的临界点提取算法是基于临界点速度分量均为0 的假设，但由于所采用的海洋流场再分析数据中可能存在计算舍入误差造成的近0 值速度向量，并且在实际的流场物理场景中，也存在速度非0 的临界点。由于Sperner 完全标号理论依据流向变化进行临界点网格的筛选，相比于双线性插值方法依据u、v流速数据求解，计算简便，为此本文也基于Sperner 完全标号理论进行了临界点的提取算法研究。其算法过程包括了：网格流向计算、临界点候选网格筛选、临界点识别等。

4.1 临界点提取与分析

该算法的执行主要包括以下几个步骤：

第一步，基于海洋流场再分析数据中的u、v速度分量，进行网格点速度向量计算；

第二步，进行临界点候选网格筛选，首先以2×2的滑动窗口对整个流场进行遍历，在滑动窗口遍历过程中，依据Sperner 完全标号理论，对4 网格进行标号处理，若4 个网格完全标号，则认为4 个网格内存在临界点，从而将其作为临界点候选网格予以保留，否则继续进行窗口滑动；

第三步，进行临界点提取，临界点候选网格筛选完成后，多数研究是通过插值算法进行临界点提取，为了具有更宽泛的适用性，本文在此提出了基于最小值法的临界点提取规则，根据该规则，将4 个候选网格中速度向量模最小的网格中心作为临界点，从而完成临界点的提取。

为了检测该算法提取的临界点的位置精度，本文还基于线性插值算法对上述过程提取的临界点进行了位置计算，具体计算结果如图8 所示，其中红色圆点表示插值位置，绿色圆点表示最小值网格的中心点。通过计算分析发现，插值圆点基本上位于最小值网格内，且网格值越小，插值点距离中心越近，当网格分辨率足够小时，两类点的距离可忽略不计，但是最小值法临界点的提取，可避免插值处理过程，其算法复杂度为O（n），比插值算法更加高效，并且可忽略插值算法中的0 值影响，因而，算法的适用性更广。

图8 插值与最小值法提取结果对比Fig. 8 Comparison between the results of interpolation and minimum method

4.2 临界点分类

依据Sperner 完全标号最小值法的临界点提取过程，同样对美国东部某海域5 000 m 深度0.033°分辨率的再分析流场数据进行了临界点提取，结果如图9a所示，同样基于雅克比矩阵计算，对提取的临界点进行了分类，具体结果如图9b 所示，共获得11 个临界点。

通过对比图7 和图9 的提取结果，发现Sperner 完全标号法虽然未能检索出更多的临界点，但的确另外提取了几个插值算法遗漏的临界点，如图9 中三角符号所示的临界点。该算法计算效率高、速度快，可较好的弥补双线性插值遗漏的临界点。

图9 基于Sperner 完全标号的临界点提取（a）与分类（b）结果Fig. 9 Extraction (a) and classification (b) of critical points based on Sperner fully labeled

5 两种算法综合的临界点提取与实验分析

从两种算法的临界点提取结果图看，双线性插值算法提取结果数量更多，覆盖范围更广，但Sperner 完全标号最小值法提取了双线性插值法未提取的临界点，并且速度更快，因此可综合两种算法临界点提取结果作为最终临界点提取结果。

5.1 算法综合处理

统计分析两种算法的分类结果，如表1 所示。

表1 两种方法分类结果统计表Table 1 Statistical table of classification results of two methods

以上两种临界点提取算法各有优缺点，通过将两种算法提取结果的合并、去重等处理，可得到如图10所示的更加全面的临界点提取与分类结果，其中，中心点7 个，交点7 个，吸引聚点5 个，排斥聚点5 个，共24 个。

图10 综合提取分类结果Fig. 10 Comprehensive extraction and classification result chart

5.2 算法验证

为了进一步验证本文算法的适用性及可行性，针对不同区域、不同深度的大量海洋流场数据进行多次验证分析，其中表2 所示，是以美国沿海区域深度5 000 m、大西洋部分海域深度2 500 m 和太平洋海域深度3 000 m 的数据为例进行的对比分析。提取结果显示：双线性插值算法临界点提取结果较为全面完整，但对非0 值临界点的提取会偶有遗漏，而Sperner完全标号法计算高效，较好地弥补了双线性插值算法遗漏的部分临界点，如表2 中美国沿海深度5 000 m数据中，弥补提取了1 个中心点、1 个交点和1 个排斥聚点；对大西洋部分海域深度2 500 m 数据中，弥补提取了1 个吸引聚点、1 个排斥聚点；对太平洋海域深度3 000 m 数据中，弥补提取了1 个吸引聚点。研究表明通过两种算法的综合处理，全面完整的提取了二维流场中的拓扑结构临界点，且算法简洁明了，可为二维流场的拓扑结构分析提供一种新的研究思路。

表2 不同数据验证结果Table 2 Different data validation results

6 结论

本文依据矢量场拓扑结构分析的临界点理论和Sperner 引理，综合改进后的双线性插值算法和Sperner 完全标号法，实现了海洋流场数据的临界点特征提取。相比于现有算法，本文算法临界点提取效果理想，且易于程序实现。改进后的双线性插值算法添加滑动窗口处理，便于临界点候选单元的筛选，涉及的聚合思想为解决插值网格二义性问题提供了新思路，而0 值网格中临界点等值线连接方式的考虑，补充了流场网格插值情形，使得提取结果更加全面，且临界点位置计算精确，可适用于多数二维流场的特征点提取。同时，基于Sperner 完全标号的临界点提取算法抛开传统插值求解临界点思想，创造性地提出了最小值法临界点提取规则，避免了插值处理过程，并且可忽略插值算法中的0 值影响，效率更高，有着更广的适用性。而两种算法作为互补，算法的综合可实现更加全面的临界点提取与分类结果，为进一步实现海洋流场拓扑分析的流线追踪奠定临界点基础。