基于拉格朗日拟序结构的黑潮区域海流运输结构提取和分析

田丰林，任一丹，何遒，陈戈*

( 1. 中国海洋大学 信息科学与工程学院，山东 青岛 266100；2. 青岛海洋科学与技术试点国家实验室 区域海洋动力学与数值模拟功能实验室，山东 青岛 266237)

1 引言

黑潮是西太平洋海域最大的暖流系统，它起源于吕宋岛以东的洋面，其主干流沿台湾以东进入东海，顺东海大陆坡向东北流去，流经日本列岛南部，于35°N，141°E附近海域结束，后续成为黑潮延伸体。黑潮是一种典型的西边界流，具有高温、高盐等特征，对其周边海区的气象、水文产生巨大影响。黑潮流量以及流幅存在明显的季节变化和年际变化，黑潮流量夏季大，秋季小，流幅遵循春夏宽，秋冬窄的变化规律[1]。目前对于黑潮的研究工作中，研究范围多为整个流系中的某一局部区域，如选取黑潮延伸体[2-4]或者台湾东北部地区[5]，有研究考察了流轴的形态变化[6-7]，黑潮流轴多采用局部最强速度场等属性等值线方法获取[8]，也有利用模型验证模拟黑潮的动力学过程[9]，对黑潮入侵南海路径和机理做了研究[10]。黑潮不断与周边水域发生强烈的物质交换和能量输送，目前关于黑潮区域拉格朗日拟序结构（LCSs）的研究较少，LCSs作为研究海洋流场结构的有力工具，能够获取通过流线和轨迹信息等无法显现的物质运输特征，对于研究黑潮环流结构、水体性质等有着重要意义。

LCSs最早在2000年时被提出[11]，它可以精确地定义和识别具有任意时间依赖性的动力系统中的障碍[12]，是流场表面存在的一种骨架结构。利用拉格朗日的方法提取相干结构方法众多，定义LCSs的方法可以有所不同，但均应满足4个条件：客观性、有限时间性、拉格朗日不变性、空间收敛性[13-14]。有研究对这些方法做了对比，探讨了不同方法得到的结果[15-16]，其中有限时间里雅普诺夫指数（FTLE）和有限尺寸里雅普诺夫指数（FSLE）是应用较为广泛的两种LCSs方法[17]，FTLE描述了初始时刻相邻的流体质点在有限时间范围里的平均分离率，FSLE表示在有限距离里的平均分离率。d'Ovidio等[18]利用FSLE研究地中海区域的混合结构识别，黄高龙等[19]利用FTLE分析吕宋海峡附近环流模式和水体输运，解释浮标运动轨迹的变化。更进一步的研究中Haller[20]及Farazmand和Haller[21]提出了变分理论来计算LCSs的方法，LCSs表现为一些具有不同收缩和拉伸特性的物质线，可以用来区分流场中具有不同动力学特征的区域，根据流体内部物质受力情况不同，分为椭圆型、双曲型和抛物线型LCSs，椭圆型LCSs是拉格朗日涡旋的边界，双曲型LCSs代表了物质收缩或拉伸线，抛物线型LCSs可以作为射流核心[22-23]，双曲型LCSs和抛物型LCSs在理论流场中的表现形式如图1所示。近些年LCSs在各种海洋现象的研究上取得了重大进展，LCSs在涡旋识别方面已有较为成熟的理论，在涡旋对水体输运效果评估方面也得到应用[24-25]，在海流方面多是对纬向射流进行分析[26-28]。Olascoaga等[29]将LCSs应用于墨西哥湾地区，表明了LCSs对该地区内部署的浮标的运动施加了重要的约束，发现吸引型LCSs近似于浮标运动形成的中尺度图案的中心，而这种物质运输模式对于传统的欧拉分析来说是非常困难的。

随着全球海洋观测技术日趋成熟，能够获取到更高质量的网格化海洋数据集，这些数据具有高时空分辨率、更长时间跨度等特性，将这些数据用于黑潮的研究，可以更加精确和直观的表达出海流信息。本文创新性地利用变分方法计算了黑潮区域的LCSs，并且将得到的LCSs简化合并为LCSs束，基于长时间序列的海流数据得到了清晰的气候态黑潮运输模式。

2 数据和方法

2.1 数据介绍

本次采用的数据是绝对动力地势（Absolute Dynamic Topography & Absolute Geostrophic Velocities，ADT）[30]，是利用在给定时间内可用的所有卫星计算出来的一种网格状海面高度和速度数据，计算方法为ADT=SLA+MDT，其中，SLA为海平面高度异常，MDT为平均动力地势。ADT数据质量高，能够突出海洋环流细节，反映出典型的海流特征，对于研究黑潮系统的结构十分有利。该数据的空间分辨率为0.25°×0.25°，时间分辨率为1 d。本研究采用了25年的数据，时间范围为1993-2017年。

我们获取了用于验证工作的浮标轨迹数据，数据来自NOAA的全球漂流浮标计划（Global Drifter Program）[31]，选取了全球范围内1993-2017年共计25年的浮标数据，按月份提取了浮标轨迹信息。

用于验证的温盐数据是来自CMEMS发布的ARMOR3D数据[32]，是一种结合了卫星观测数据和现场数据而得到的数据产品，数据空间分辨率为0.25°×0.25°，采用的时间范围也是1993-2017年。

2.2 研究方法

2.2.1 计算LCSs的变分方法

在二维不稳定流场中，速度场可以表示如下：

柯西格林应变张量场可以表示为[21]

图1 双曲线型 LCSs（a）和抛物线型 LCSs（b）Fig. 1 Hyperbolic LCSs（a）and parabolic LCSs（b）

应变线处处与单位应变特征向量场n1相切，可以根据下式中的常微分方程计算得到平滑的轨迹

为使此方程的数值解为一个全局平滑的应变线集，依据一般张量线的尺度变换方式，将式（5）进一步表示为[21]

利用变分方法从速度场中计算得到LCSs是一些特殊的物质线，一条物质线处处切于时为收缩线，处处切于时为拉伸线[23,30]。对于，当后向积分（T＜0）时，吸引型 LCSs定义为局部最大收缩线，排斥型LCSs定义为局部最大拉伸线，前向积分（T＞0）时，吸引型LCSs定义为局部最大拉伸线，排斥型LCSs定义为局部最大收缩线。图1所示的双曲型LCSs，中心横线代表的是拉伸线，中心竖线代表的是收缩线，虚线表示周围的流形，它们将流场划分为4个具有不同运输特征的区域，可见收缩线作为流场中的运输屏障，阻止了其两边的水体和物质交换[23,33]。

2.2.2 计算气候态LCSs

在气候态下研究LCSs具有更广泛的意义，能够反映当前海流区域的综合特征和集中趋势，为不同月份和不同季节提供对比分析。我们利用1993-2017年总共25年的速度场数据，通过年际求平均得到气候态海流速度场，用计算得到的气候态海流速度场做进一步的研究。为方便计算，取每年前360 d的数据，每个月按照30 d计算，1年共计12个月[34]。

利用密集流体粒子在流场中运动一段时间后计算得到柯西格林应变张量场，张量场描述了粒子在当前位置下的运动趋势和形变特性。我们在整个区域内均匀选取N×N个网格点，在每一个网格点周围放置4个距离较近的辅助点，如图2所示。本文中辅助点与主网格点之间的距离δx1、δx2设定为两个网格点之间距离的0.05倍，使用龙格-库塔四阶积分计算辅助点处的粒子运动后的位置，使用有限差分方法计算出每个网格点处的形变梯度矩阵，得到柯西格林应变张量，进而计算特征值和特征向量，积分得到LCSs。

图2 网格点与辅助点的位置Fig. 2 Locations of grid points and auxiliary points

本文采用后向积分的方法，将吸引型LCSs作为研究对象，从每个月里选取第30天代表计算LCSs的初始时刻，积分时间为30 d，计算张量场，基于此张量场计算得到的12个不同动力系统的LCSs，分别代表了各月的气候态LCSs。具体步骤如下：

（1）生成均匀密集网格点，求每点处的柯西格林应变张量，并计算特征值和特征向量；

（3）将所有极值点按特征值大小递减排序；

（4）以特征值极值点为种子点，在对应的特征矢量场中积分，得到吸引型LCSs。

通过调整寻找极值点的范围和积分线的密集程度等参数可以保留部分吸引型LCSs，如图3所示，得到了代表12个月份的气候态LCSs，图中的线代表了求取的吸引型LCSs，表示其法向排斥率，我们将线的颜色用ln值划分。利用FTLE计算得到的LCSs体现为标量场的脊，在有限时间[]内，FTLE定义如下：

我们将获取吸引型LCSs与FTLE场进行叠加，可以发现抽取到的吸引型LCSs与FTLE的脊基本相吻合，在这种条件下该吸引型LCSs可以作为黑潮两侧边界。

2.2.3 提取 LCSs束

由于过少的LCSs无法全面地展示海流结构，过多的LCSs容易引起视觉杂乱，我们获取整个研究区域内的全部吸引型LCSs，如图4a所示，当LCSs足够密集时，把积分得到的所有吸引型LCSs进行简化合并。保留前20%和40%的吸引型LCSs（图4b，图4c），观察其覆盖的范围，发现黑潮区域的特征值要明显大于周边区域，选取的吸引型LCSs不仅要表现出黑潮区域的整体结构，还摒弃大部分周围海流结构的干扰。

图3 气候态下的吸引型LCSs和FTLE场Fig. 3 Monthly climatological attracting LCSs and FTLE fields

图4 保留 100%（a），40%（b）和 20%（c）的吸引型 LCSsFig. 4 Reserved attracting LCSs for 100% (a), 40% (b) and 20% (c)

选取前40%的吸引型LCSs进行研究，经过简化合并后的LCSs我们称之为LCSs束，具体合并流程以及示例如图5所示，对于每一条吸引型LCSs，计算它与区域内其他线的相似度，将相似度最高的两条线合并为一条线。在对任意的两条线计算相似度时，由于线是由有序排列的顶点构成，对于一条线中的每个顶点坐标，去寻找另一条线中距离最近的点进行配对并计算距离，最终得到以该条线为基准的平均距离，用同样的方法得到以另一条线为基准的平均距离，对这两个距离求取平均，得到的平均值作为这两条线相似度的数值描述[35]。通过试验设定容许两线合并的最大相似度S=1°，将两条线按照逐点对取中间点或者取较长线的方法进行合并，随着保留的线越来越少，计算复杂度逐渐降低，最终能够得到反映流场突出特征的更为简洁的LCSs束。

图5 LCSs简化合并流程图（a）和相似度表示例（b，c）Fig. 5 Flow chart of simplifing and merging LCSs (a) and an example of similarity table (b, c)

将吸引型LCSs进行简化合并后，得到的效果如图6所示，合并的线的条数越多，LCSs束在图中表现得越粗，原来密集的LCSs经过合并后能够清晰地表现出特征较强的流场区域，反映出黑潮的主轴结构，该方法能够保持大范围流场中突出的海洋特性，也有利于海流结构的二维可视化研究。

3 结果与讨论

3.1 结果分析

用气候态速度场计算得到的吸引型LCSs，在黑潮主干区域能够得到较连贯黑潮流结构，在黑潮延伸体区域，虽然不同时间的流场变化幅度大，复杂且不稳定，但是通过气候态数据的计算也可以反映出典型的运输特性。图7是提取到的黑潮系统的流场结构图，揭示了12个月份的拉格朗日运输模式。

目前，黑潮流轴和路径的获取方法多种多样，最为常用的有以下3种方法：其一是基于流线的识别方法，黑潮表面流速最大的部分看作是表面黑潮流轴；其二是基于等温线的识别方法，黑潮海水的温度显著高于周围海流，日本南部海域处的主轴采用 16℃等温线进行识别，而黑潮延伸体处的主轴则用12℃等温线作为识别依据；其三是基于高度数据的识别方法，例如对于ADT数据，ADT=240 cm处的等值线可以作为黑潮主轴[36]。图8展示了利用流速和海表面高度识别的黑潮流轴。

图6 合并前的LCSs (a)和合并后的LCSs束(b)示意图Fig. 6 LCSs before merging (a) and LCSs bundles after merging (b)

图7 12个月的气候态流场结构Fig. 7 Monthly climatological flow field structures

图8 基于多年平均流速提取的黑潮主轴[8]（a）和 170 cm海表面高度轮廓线定义的黑潮延伸体路径[4]（b）Fig. 8 Kuroshio axis extracted based on average velocity[8] (a) and paths of the Kuroshio Extension defined by the 170 cm contours in sea surface height fields[4] (b)

本文提出了一种用LCSs束识别黑潮流路的方法，LCSs束本身代表了两侧海流具有较强的排斥作用，LCSs束可以作为提取黑潮主轴的有力工具，将其应用于海洋流场中具有重要价值。从图7中可以发现黑潮主干结构稳定，但在局部地区存在月际变化。在流经吕宋海峡时，产生不连续的隔离现象，这可能是因为黑潮入侵南海所导致。在黑潮主干区域，冬季提取到的流结构较为紧密和稳健，可能是由于冬季黑潮流幅窄而平稳的缘故。黑潮在经过台湾东北部时由于海底地形的作用产生向东的转折，分出一支台湾暖流，从图7中发现春季在这一地区的转折相较于其他时期更为平缓一些，而后主干沿着大陆架的约束继续向东北方向前进，直到达到日本九州岛西南部时又出现转折，在此之前一直存在着稳定和强健的吸引型LCSs。根据图7中线的颜色可以看出在转折区域多会发现特征强度较弱的信号，可能是因为产生流系分支的原因，这些区域的流轴不是很健壮。由于黑潮延伸体区域流场的复杂多变性，提取到的LCSs束的连续性较弱，但是从图7中可以发现延伸体在143°～150°E之间，不同月份有着相似弯曲形态，均在146°E处形成谷状，且冬季弯曲幅度更大。黑潮延伸体区域存在中尺度及亚中尺度涡旋等复杂流场结构，海流在流动过程中存在较强的物质交换，使得海流边界更为曲折复杂，难以提取到完整的结构。

3.2 验证

为了评估以上结果的可靠性，我们对结果进行了验证，验证方法包括拉格朗日粒子运输验证、浮标轨迹验证以及温度盐度验证。

3.2.1 拉格朗日粒子的运输

采用1月份的数据进行拉格朗日粒子运输的验证，通过示踪粒子在不稳定流场中运动一段时间后观察其是否与LCSs束吻合。选定6个不同位置且间距合理的区域，在直径为1°的圆形区域内均匀地放置示踪粒子，它们的初始位置和形态如图9a所示，让粒子在1月的不稳定流场中分别运动10 d、20 d和30 d，得到的结果如图9b至图9d所示。可以发现粒子输送总是沿着LCSs的方向拉伸，初始时刻的圆形组合随着时间的增加逐渐变为回旋镖状，最后拉伸成丝状，证明LCSs束可以代表黑潮的主轴结构，也可预测物质运输的形态。

3.2.2 浮标轨迹

我们通过1993-2017年的浮标轨迹信息对结果进行验证，提取了25年中1月和6月的全部浮标轨迹，与相应月份的LCSs束进行比较，结果如图10所示，浮标的运输模式符合LCSs束的预测情况。虽然单独的浮标运动由于受局部因素影响不确定性大，轨迹较为曲折和杂乱，但从整体上看它们的轨迹分布基本覆盖了整个黑潮流域，浮标沿着黑潮蜿蜒前行，在台湾岛东北部、九州岛南部等区域发生转折，在黑潮主干区域与LCSs束的提取结果保持了高度的一致性，在黑潮延伸体区域也存在部分相似性。

3.2.3 温盐分布

此前相关研究多采用温盐等值线来确定黑潮的位置[37]，基于这种方法我们分析了温盐分布与提取到的流场结构的相关性。采用海表面和50 m深的温度和盐度数据，得到了25年的1月份温盐数据的平均值，计算了温盐异常，将其与1月份提取到的LCSs束叠加，结果如图11所示。从温度异常图中可以看出，黑潮北侧的温度异常值变化明显，LCSs束有效区分了不同温度的区域。盐度异常图中，黑潮主干大部分地区的LCSs束存在于盐度异常变化相对较大的区域。总体上看，黑潮的北侧温盐异常梯度大，边界明显，南侧温盐边界较为模糊，与提取的LCSs束具有一致性。而在黑潮延伸体区域，温盐异常分布与LCSs束的相似性较小，这部分区域是海洋涡旋活跃地带，物质交换剧烈，难以获得精确结果。

图9 示踪粒子初始位置（a）和运动 10 d（b）、20 d（c）、30 d（d）后的 LCSs输运图Fig. 9 Initial position (a) and advected images for 10 days (b), 20 days (c) and 30 days (d) of particles

图10 1月（a）和 6 月（b）的浮标轨迹Fig. 10 Trajectories of drifters for January (a) and June (b)

图11 气候态温盐异常分布Fig. 11 Distributions of climatological temperature and salt anomalies

4 总结

本文研究了黑潮区域的拉格朗日运输模式，面对复杂多变的流场环境，提出了利用气候态LCSs束来计算流场结构的方法。整合了25年的速度场数据，通过变分方法计算得到了气候态下的吸引型LCSs，并根据特征值的强度选取了黑潮附近的LCSs，将其简化合并为LCSs束。LCSs束可以作为描述黑潮形态和流场结构的工具，用于研究黑潮的时间和空间变化规律等。我们提取并分析了12个月份的气候态流场结构图，讨论了不同月份和季节中物质运输的一致性和区别，最后用拉格朗日虚拟粒子、真实浮标轨迹以及温度盐度数据作为验证，表明了其合理性和准确性。

同时本文也存在诸多不足，首先我们只关注了黑潮及其延伸体的主干部分，未来可以考虑黑潮分支的提取，如台湾暖流、对马暖流等，将研究拓展到完整的黑潮流系结构。其次，在黑潮延伸体区域，可以结合中尺度涡旋的位置、特性以及其他海洋现象展开分析，进而获得更为连续和完整的海流结构。同时，可以将本方法用于其他大范围的海域，利用并行方案来解决大数据问题，对全球海流的给定精确描述，为物质运输状态做指导，应用于海上溢油、海上垃圾堆积和扩散区域的识别和预测。