CN117272615A - 一种基于流场速度矢量象限角的海洋涡旋识别方法 - Google Patents

Abstract

本发明设计一种基于流场速度矢量象限角的海洋涡旋识别方法，属于海洋科学和遥感领域；首先利用MATLAB读取存有地转流分量的卫星测高产品中的相关数据，并利用获取的相关数据求出所选研究范围内网格点的速度矢量象限角方向；根据自定义矩形搜索框r*r四个顶点的速度矢量象限角值的关系确定涡旋中心；在涡旋中心附近范围内利用卫星测高产品中的相关数据求出流函数值，进而求出最外层闭合等值线即为涡旋的边界；将无边界以及半径小的涡旋中心剔除，其余的涡旋中心按日期存入一个结构体中；然后将连续的时间步长数据整合到一个数组中，从而得到涡旋在不同时间点的位置信息，通过追踪连续时间数据中涡旋位置的变化进而实现对其时序运动的跟踪。

Description

一种基于流场速度矢量象限角的海洋涡旋识别方法

技术领域

本发明属于海洋科学和遥感领域，特别涉及一种基于流场速度矢量象限角的海洋涡旋识别方法。

背景技术

海洋涡旋是一种旋转的、以封闭环流为主要特征的水体，它是由各种气象因素作用和海洋动力不稳定造成的，并广泛分布于全球海域。中尺度涡旋具有极大的动能，对海洋环流、大气环境、全球温度变化、海洋能量传输都起着重要作用。作为海洋动力的重要组成部分，涡旋的形成、移动和消散一直都是重要的研究方向，因此为了研究长期的涡旋运动变化，自动涡旋跟踪算法的研发成为近年来探测海洋涡旋的主要手段之一。

目前，关于海洋涡旋自动识别和跟踪算法的研究，主要存在以下几个代表性的方法：(1)Isern-Fontanet(2003)首次提出OW(Okubo-Weiss)算法，该方法利用地转表面的速度场，基于涡旋的形变率和相对涡度特征进行自动探测，并取得了广泛应用。但是但该方法的探测结果受阈值的影响很大，且物理参数存在噪声项，物理标准可能导致检测失败或者低估涡旋的大小。(2)Wang，(2003)和Chelton D B等(2011)利用海面高度异常数据的等高线来识别涡旋中心及其特征，即闭合等值线法，这类方法需要在海面高度差上设置一个阈值，进而识别出涡旋的边界，但是阈值微小的变化将对识别的结果产生较大的影响。(3)Nencioli(2010)提出了矢量几何法，然而，该方法在搜索区域较小的情况下存在一些问题。当搜索区域较小时，速度的极小值点会增加，用于第四约束条件检测的点增多将增加错误识别涡旋中心的概率。另一方面，如果搜索区域较大，则可能只有较大的涡旋会被检测到。

由此，上述现有算法的研究虽然为海洋涡旋的识别与跟踪技术奠定了基础，但是仍然存在一些问题需要解决。

发明内容

针对上述现有技术中的不足，本发明提出了一种基于流场速度矢量象限角的海洋涡旋识别方法，能够快速、准确的探测到海洋中的涡旋，克服了现有技术中的探测效率局限性和对海洋涡旋识别遗漏的缺陷，为跟踪涡旋运动特征和规律提供了有效的技术和数据支持。

一种基于流场速度矢量象限角的海洋涡旋识别方法，具体包括以下步骤：

步骤1：利用MATLAB读取存有地转流分量的卫星测高产品中的相关数据；所述相关数据包括：经度、纬度、海面高度异常、地转流东西速度分量U以及地转流南北速度分量V；

步骤2：利用步骤1获取的相关数据求出所选研究范围内网格点的速度矢量象限角方向；

所述步骤2是根据地转流东西速度分量U及南北速度分量V的正负关系定义速度矢量象限角的方向，具体为：

(1)U和V都大于零时为第一象限；

(2)U小于零而V大于零时为第二象限；

(3)U和V都小于零时为第三象限；

(4)U大于零而V小于零时为第四象限；

步骤3：根据自定义矩形搜索框r*r四个顶点的速度矢量象限角值的关系确定涡旋中心；

由于在北半球(南半球)，气旋涡的旋转方向为逆时针(顺时针)，反气旋涡的旋转方向为顺时针(逆时针)；所以根据旋转方式以及象限角的关系，利用矩形搜索框四个网格顶点速度方向排序关系确定涡旋的中心，具体满足下面三个约束条件：

(1)将四个网格点中最小经度最小纬度的点为起始点，逆时针排序；无论气旋涡的四个点A1，A2，A3，A4还是反气旋涡的四个点B1，B2，B3，B4四个点的速度象限角值累加和都为10；

(2)在满足约束1的条件下，继续以气旋涡和反气旋涡的四个象限角值排序方式为约束条件；

(3)要求在满足条件1和条件2的这些网格点的附近存在一个局部速度的最小值；

步骤4：在涡旋中心附近范围内通过流函数公式，利用U、V、网格点间距离求出流函数值，通过MATLAB的contourc函数求出最外层闭合等值线即为涡旋的边界；将无边界以及半径小于b的涡旋中心剔除，其余的涡旋中心按日期存入一个结构体中，具体为：Eddy_data＝('day',{},'type',{},'lon',{},'lat',{},'radius',{})；

流函数计算公式：

具体涡旋中心附近范围确定方法如下：设置一个参数a，利用参数a确定计算流函数值的范围；首先在涡旋中心位置出发，向上、下、左、右四个方向分别移动2a个网格点，创建一个(4a+1)×(4a+1)个网格点的区域；在这个区域内，计算出最大闭合等值流线；如果这些闭合等值流线的四个端点距离该矩形区域边界的长度小于1.5个网格间距，需扩大搜索区域；通过在各个方向上各扩展一个网格点，重新计算流函数值，从而获得涡旋的新的最大闭合等值流线；重复上述过程，直至找到涡旋最大边界，同时剔除那些没有明确定界的涡旋中心；

步骤5：利用步骤4方法获取得到连续的时间步长数据整合到一个数组中，从而得到涡旋在不同时间点的位置信息，通过追踪连续时间数据中涡旋位置的变化进而实现对其时序运动的跟踪；

将结构体Eddy_data的所有数据遍历重新放到一个数组data＝[T,lon,lat,type]，这个数组包含了所有连续日期的数据，所以能够获取涡旋在不同时间点的位置信息，进而实现对其时序运动的跟踪；

具体为：通过搜索T+1时刻相同类型的涡旋中心位置，以更新涡旋的运动轨迹；如果在T+1时刻的搜索区域内无法检测到涡旋中心，将进行第二次搜索，若在T+2时刻的搜索区域内仍然无法检测到涡旋中心，将判定该涡旋的寿命已经结束，相应地终止其轨迹；在T+1时刻，如果某个涡旋中心无法与T时刻的任何涡旋相连接，则被视为新形成的涡旋，其轨迹将从T+2时刻开始更新。

本发明有益技术效果：

(1)本发明利用流场的几何特征，开发了一种基于流场速度矢量象限角的涡旋检测方法。以象限角的累加和方向连续性作为基础条件，利用流场中四个速度矢量方向象限角值从而识别涡旋位置，有效提高了涡旋识别的准确性以及提升了涡旋的探测效率。

(2)本发明利用闭合等值线法以及矢量几何法两种涡旋探测方法进行验证，可以证明利用该方法识别涡旋位置准确、遗漏少和效率高。

附图说明

图1是本发明一种基于流场速度矢量象限角的海洋涡旋识别方法流程图；

图2是本发明实施例北半球气旋涡与反气旋涡速度流向以及象限角值定义示意图；

图3为本发明实施例约束3求局部最小值的原理图；

图4是本发明实施例在某海西部海域识别出的涡旋分布情况(★为气旋涡、●为反气旋涡)；

图5是2022年冬季涡旋运动轨迹图(O为涡旋运动起点，X为涡旋结束终点)；

图6是本发明实施例和矢量几何法探测的涡旋分布图对比图；图a)为矢量几何法在某海区域探测的涡旋分布图(★为气旋涡、●为反气旋涡)；图b为本发明在某海区域探测的涡旋分布图(★为气旋涡、●为反气旋涡)。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

本实施采用的数据集来源于哥白尼海事服务(Copernicus Marine Service)的多源卫星融合高度计产品，即网格化SLA数据产品，该产品主要由Altika、HY2B、Jason3、Sentinel系列等多测高卫星数据融合而成。该数据集提供了全球网格的MSLA(距平值是相对于1993～2012年的平均值)、地转流速度等。本发明主要利用地转流速(U、V)和SLA数据识别中尺度涡旋。

本发明通过编制matlab程序实现，如图1所示，基于流场速度矢量象限角的海洋涡旋识别与跟踪方法包括以下步骤：

一种基于流场速度矢量象限角的海洋涡旋识别方法，具体包括以下步骤：

步骤1：利用MATLAB读取存有地转流分量的卫星测高产品中的相关数据；所述相关数据包括：经度、纬度、海面高度异常、地转流东西速度分量U以及地转流南北速度分量V；并选择某海为研究区域。

步骤2：利用步骤1获取的相关数据求出所选研究范围内网格点的速度矢量象限角方向；

所述步骤2是根据地转流东西速度分量U及南北速度分量V的正负关系定义速度矢量象限角的方向，具体为：

(1)U和V都大于零时为第一象限；

(2)U小于零而V大于零时为第二象限；

(3)U和V都小于零时为第三象限；

(4)U大于零而V小于零时为第四象限；

步骤3：根据自定义矩形搜索框r*r四个顶点的速度矢量象限角值的关系确定涡旋中心；

所述矩形搜索框为人为设置的一个搜索范围，然后遍历循环整个研究区域，假如研究区域是个50*50的范围的网格，这个搜索框就从第一排第一个点到第一排最后，然后第二排第一个到最后直至所有的排列都循环一遍；搜索框的边长就是r个网格点；

由于在北半球(南半球)，气旋涡的旋转方向为逆时针(顺时针)，反气旋涡的旋转方向为顺时针(逆时针)；如图2所示，所以根据旋转方式以及象限角的关系，通过设置参数r＝3，所述参数根据卫星测高产品数据分辨率的大小设置；利用矩形搜索框四个网格顶点速度方向排序关系确定涡旋的中心，具体满足下面三个约束条件：

(1)将四个网格点中最小经度最小纬度的点为起始点，逆时针排序；无论气旋涡的四个点A1，A2，A3，A4还是反气旋涡的四个点B1，B2，B3，B4四个点的速度象限角值累加和都为10；

(2)在满足约束1的条件下，继续以气旋涡和反气旋涡的四个象限角值排序方式为约束条件；

(3)要求在满足条件1和条件2的这些网格点的附近存在一个局部速度的最小值；

在北半球的气旋涡中，通过逆时针旋转的方式对矩形的四个顶点A1、A2、A3、A4进行象限角值的排序。具体而言，对于气旋涡，排序方式依次为第四象限、第一象限、第二象限、第三象限。针对反气旋涡，通过逆时针旋转的方式对四个顶点B1、B2、B3、B4的象限角值进行排序，依次为第二象限、第三象限、第四象限、第一象限。无论是气旋涡还是反气旋涡，都满足了一个关键约束条件：即在r×r的矩形搜索区域中，四个顶点的象限值累加和始终为10。这个条件是为了初步筛选满足要求的网格点，有助于加快识别涡旋的特征。3.2从图2可知，气旋涡和反气旋涡的矩形边框四个顶点的速度方向象限角值各不相同，且这两种涡旋的四个顶点速度象限角值分别位于相邻的象限内。接下来，需要对矩形边框四个点的象限值进行排序。若排序结果为第四象限、第一象限、第二象限、第三象限，表示该区域可能是一个气旋涡(南半球则可能是反气旋涡)的近似中心位置；若排序结果为第二象限、第三象限、第四象限、第一象限，则可能是一个反气旋涡(在南半球则可能是气旋涡)的大致中心位置。实际上，约束条件二已经包含了约束条件一所涵盖的约束范围，然而，出于计算效率的考虑，保留约束条件一仍然十分必要。

如图3所示，此时实心圆圈(·)是满足条件一和条件二的点，星号(*)代表真实的涡旋中心点。采用如下策略来确定搜索涡旋中心点。首先，选择一个尺寸为r×r的网格，该网格覆盖了可能包含涡旋中心的区域，为了保证找到的是局部最小值，首先进行第一次搜索，寻找在当前网格范围内的最小值。继而以它为中心再次进行r×r网格区域搜索，以拓宽搜索范围。需要注意的是，如果r是偶数，那么第二次搜索的矩形边长网格点将会是r+1。在两次搜索中，如果找到的最小值相同，将确认该点为局部最小值，初步假定该点为涡旋的中心位置。

步骤4：在涡旋中心附近范围内通过流函数公式，利用U、V、网格点间距离(0.25°对应在地球表面的距离)求出流函数值，通过MATLAB的contourc函数求出最外层闭合等值线即为涡旋的边界；将无边界将以及半径小于b(这里面设置为30km，根据卫星测高产品分辨率设置)的涡旋中心剔除，其余的涡旋中心按日期存入一个结构体中，具体为：Eddy_data＝('day',{},'type',{},'lon',{},'lat',{},'radius',{})；

流函数计算公式：

具体涡旋中心附近范围确定方法如下：设置一个参数a＝2，利用参数a确定计算流函数值的范围；首先在涡旋中心位置出发，向上、下、左、右四个方向分别移动2a个网格点，创建一个9×9个网格点的区域；在这个区域内，计算出最大闭合等值流线；如果这些闭合等值流线的四个端点距离该矩形区域边界的长度小于1.5个网格间距，需扩大搜索区域；通过在各个方向上各扩展一个网格点，重新计算流函数值，从而获得涡旋的新的最大闭合等值流线；重复上述过程，直至找到涡旋最大边界，同时剔除那些没有明确定界的涡旋中心；

图4给出了利用本发明在某海西部海域识别出的涡旋分布情况，并呈现了这些涡旋的边界轮廓，通过将海平面异常(SLA)数据添加为底图。涡旋的半径通过边界计算确定，这一过程涉及计算涡旋边界上各点到涡旋中心点的距离，然后取这些距离的平均值作为涡旋的半径。由于本发明所用的数据的分辨率为0.25°*0.25°，所以将半径小于30km的涡旋中心剔除，最后将涡旋的中心按日期存入结构体中。Eddy_data＝('day',{},'type',{},'lon',{},'lat',{},'radius',{})。依次存入的数据为时间、涡旋类型、经度、纬度、以及涡旋的半径。

步骤5：利用步骤4方法获取得到连续的时间步长数据整合，从而得到涡旋在不同时间点的位置信息，进而实现对其时序运动的跟踪；

将连续的时间步长数据整合到一个数组中，能够获取涡旋在不同时间点的位置信息，通过追踪连续时间数据中涡旋位置的变化进而实现对其时序运动的跟踪。

将结构体Eddy_data的所有数据遍历重新放到一个数组data＝[T,lon,lat,type]，这个数组包含了所有连续日期的数据，所以能够获取涡旋在不同时间点的位置信息，进而实现对其时序运动的跟踪。通常情况下，涡旋的运动速度较缓，而洋流的平均速度约为0.2m/s^(-1)。因此，本发明以3×3网格点搜索区域足以捕捉相邻两天中相同的涡旋中心位置。在这个区域内，搜索T+1时刻相同类型的涡旋中心位置，以更新涡旋的运动轨迹。然而，由于数据源的局限性或其他原因，有时可能无法识别出T+1时刻的涡旋中心，而T和T+2时刻的位置能够被准确检测到，这将导致一个涡旋轨迹被错误地分割成两个。

因此，在下一步中，为了减少这种误差，如果在T+1时刻的搜索区域内无法检测到涡旋中心，将进行第二次搜索，将搜索区域扩展至5×5个网格点。这不仅可以防止由检测错误导致涡旋轨迹的分裂，还可以避免过早的终止涡旋的寿命。若在T+2时刻的搜索区域内仍然无法检测到涡旋中心，将判定该涡旋的寿命已经结束，相应地终止其轨迹。在T+1时刻，如果某个涡旋中心无法与T时刻的任何涡旋相连接，则被视为新形成的涡旋，其轨迹将从T+2时刻开始更新。图5展示利用本发明跟踪的2022年冬季某海区域涡旋的空间移动轨迹，其中涡旋的产生和消散分别用(○)和(X)表示。

为了进一步验证涡旋探测结果的准确性，选用矢量几何法对相同区域的海洋涡旋探测结果进行对比分析，如图6所示，从本发明的算法与矢量几何法得到的涡旋分布图对比中，我们可以观察到它们之间的一致性。除此之外，本发明的算法可以识别出一些半径较小的涡旋，而这些涡旋可能无法通过矢量几何法有效识别。

Claims (7)

1.一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，具体包括以下步骤：

步骤1：利用MATLAB读取存有地转流分量的卫星测高产品中的相关数据；所述相关数据包括：经度、纬度、海面高度异常、地转流东西速度分量U以及地转流南北速度分量V；

步骤2：利用步骤1获取的相关数据求出所选研究范围内网格点的速度矢量象限角方向；

步骤3：根据自定义矩形搜索框r*r四个顶点的速度矢量象限角值的关系确定涡旋中心；

步骤4：在涡旋中心附近范围内通过流函数公式，利用U、V、网格点间距离求出流函数值，通过MATLAB的contourc函数求出最外层闭合等值线即为涡旋的边界；将无边界以及半径小于b的涡旋中心剔除，其余的涡旋中心按日期存入一个结构体中，具体为：Eddy_data＝('day',{},'type',{},'lon',{},'lat',{},'radius',{})；

步骤5：利用步骤4方法获取得到连续的时间步长数据整合到一个数组中，从而得到涡旋在不同时间点的位置信息，通过追踪连续时间数据中涡旋位置的变化进而实现对其时序运动的跟踪。

2.根据权利要求1所述的一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，所述步骤2是根据地转流东西速度分量U及南北速度分量V的正负关系定义速度矢量象限角的方向，具体为：

(1)U和V都大于零时为第一象限；

(2)U小于零而V大于零时为第二象限；

(3)U和V都小于零时为第三象限；

(4)U大于零而V小于零时为第四象限。

3.根据权利要求1所述的一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，所述由于在北半球(南半球)，气旋涡的旋转方向为逆时针(顺时针)，反气旋涡的旋转方向为顺时针(逆时针)；所以根据旋转方式以及象限角的关系，利用矩形搜索框四个网格顶点速度方向排序关系确定涡旋的中心，具体满足三个约束条件。

4.根据权利要求3所述的一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，所述约束条件具体为：

(1)将四个网格点中最小经度最小纬度的点为起始点，逆时针排序；无论气旋涡的四个点A1，A2，A3，A4还是反气旋涡的四个点B1，B2，B3，B4四个点的速度象限角值累加和都为10；

(2)在满足约束1的条件下，继续以气旋涡和反气旋涡的四个象限角值排序方式为约束条件；

(3)要求在满足条件1和条件2的这些网格点的附近存在一个局部速度的最小值。

5.根据权利要求1所述的一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，步骤4所述流函数计算公式：

6.根据权利要求1所述的一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，步骤4具体涡旋中心附近范围确定方法如下：设置一个参数a，利用参数a确定计算流函数值的范围；首先在涡旋中心位置出发，向上、下、左、右四个方向分别移动2a个网格点，创建一个(4a+1)×(4a+1)个网格点的区域；在这个区域内，计算出最大闭合等值流线；如果这些闭合等值流线的四个端点距离该矩形区域边界的长度小于1.5个网格间距，需扩大搜索区域；通过在各个方向上各扩展一个网格点，重新计算流函数值，从而获得涡旋的新的最大闭合等值流线；重复上述过程，直至找到涡旋最大边界，同时剔除那些没有明确定界的涡旋中心。

7.根据权利要求1所述的一种基于流场速度矢量象限角的海洋涡旋识别方法，其特征在于，步骤5具体为：

将结构体Eddy_data的所有数据遍历重新放到一个数组data＝[T,lon,lat,type]，这个数组包含了所有连续日期的数据，所以能够获取涡旋在不同时间点的位置信息，进而实现对其时序运动的跟踪；

具体为：通过搜索T+1时刻相同类型的涡旋中心位置，以更新涡旋的运动轨迹；如果在T+1时刻的搜索区域内无法检测到涡旋中心，将进行第二次搜索，若在T+2时刻的搜索区域内仍然无法检测到涡旋中心，将判定该涡旋的寿命已经结束，相应地终止其轨迹；在T+1时刻，如果某个涡旋中心无法与T时刻的任何涡旋相连接，则被视为新形成的涡旋，其轨迹将从T+2时刻开始更新。