基于卫星测高数据的海洋中尺度涡流动态特征检测

赵文涛，俞建成，张艾群，李　岩

(1. 中国科学院 沈阳自动化研究所, 辽宁 沈阳 110016； 2. 中国科学院大学, 北京 100049)



基于卫星测高数据的海洋中尺度涡流动态特征检测

赵文涛1,2，俞建成*1，张艾群1，李岩1

(1. 中国科学院 沈阳自动化研究所, 辽宁 沈阳 110016； 2. 中国科学院大学, 北京 100049)

为了最终实现对海洋中尺度涡流(简称中尺度涡)的自动采样，首先应该发展中尺度涡动态特征识别技术。本文基于SLA(Sea Level Anomaly)数据，实现了对中尺度涡动态特征的检测算法。主要内容是制定了一个判别相邻两组SLA数据中的涡流，是否为同一涡流子在不同时刻的状态的标准，即判别下一时刻SLA数据中是否存在涡流是由上一时刻确定的被检测涡流演化而来的。通过确定这种进化关系，可以得到被检测涡流的一系列动态状态信息，例如：面积变化速率、中心移动情况以及其他情况。本算法的计算量不大，从而可以应用到实时涡流跟踪的环境中。值得注意的是，本文中的算法不仅仅局限于应用SLA数据，SSH(Sea Surface Height)等大部分反映海洋高度的数据也可以使用。

中尺度涡；动态特征；演化关系；自动检测

0　引言

海洋热量和物质输送过程对全球气候有深远影响[1]。中尺度涡对于物质的输送作用几乎可以和风生流和热生流相媲美。中尺度涡流区域在海洋活动中的重要影响，使科学家们对其产生了浓厚的兴趣。随着集成电路技术的发展，利用AUV等自主平台进行长时间海洋特征观测已经成为可能。虽然已经有相关研究人员对中尺度涡进行采样调查[2]。但是，到目前为止还没有很好的方法可以对中尺度涡实现自主采样，不能够准确地采集敏感区域的特征数据，反应涡流区域的内部特征。

掌握中尺度涡流区域的移动路径、面积变化、形状变化等特征可以为研究者提供必要的信息，从而为AUV采样策略的制定提供依据，最终提高对中尺度涡进行跟踪观测的准确性，使采集的数据更加有效地反映涡流区域的内部特征。鉴于卫星数据可以有效地反映涡流区域的宏观动态特征，利用卫星数据自主识别中尺度涡流区域的移动路径、面积变化、形状变化等动态特征的方法能够有效减轻研究者的工作强度，提高效率。

PENVEN et al[3]提出了一种涡流自动检测的方法。但是他们将涡流区域简单地作为一个圆形区域对待，这种方式忽略了涡流形状的信息。同时在他们的算法中使用了汉宁滤波器，这使初生期的涡流很难被发现，因为这个时候的涡流区域还很小，特征表现还不明显。在他们的算法中没有考虑涡流的移动速度对动态特征识别的影响，所以当多个涡流距离很近的时候，很容易发生进化关系检测的错误。虽然CHAIGNEAU et al[4]将PENVEN et al[3]的算法进行了改进，将EKE(Eddy Kinetic Energy) 加入到进化关系判别的标准中，但是其没有针对算法的确定进行改进。

基于ISERN-FORTANET et al[5-6]的算法, CHELTON et al[7-8]提出了一种涡流自主检测的算法。 但是他们只针对这种方法进行了简单的说明，没有给出具体的算法和表达式。同时其算法中也缺少对于中尺度涡进行关系判别的说明。在 MGET 工具箱中[9], 涡流之间的进化关系主要是依靠两组海面测高数据中涡流区域的重复情况进行判别的。这种方法没有考虑涡流形状的变化等信息。如果两个涡流距离很近，也会引起进化关系判别的错误。

对于涡流区域的判别，有很多方法可以做到。BERON-VERA et al[10]提出一种客观的(即坐标系独立的)方法可以用于进行涡流边界的识别。ISERN-FORTANET et al[5]提出一种根据高度信息进行涡流中心区域检测的方法。OW(Okubo-Weiss)参数可以用来进行SLA(Satellite Level Anomaly)数据中涡流区域的检测，ISERN-FORTANET et al[11]给出了OW参数计算的过程，以及过程中所需的速度场的计算公式。SLA数据中满足OW涡流主导区域判别标准的点成为OW点，不满足的则称为非OW点。

本研究将涡流中心的运动速度、涡流区域的形状变化、面积变化这3个因素作为涡流进化关系判别方法中的主要参数，对涡流区域检测进行了简单的说明，对涡流进化关系识别进行了具体说明，最后给出了基于SLA数据的中尺度涡动态特征识别的仿真结果。

1　涡流区域检测

1.1OW参数

ISERN-FORTANET等人给出了如下所示的OW参数计算公式[11]。

2D流场的计算公式：

(1)

OW参数的计算公式：

(2)

公式(1)中g代表重力加速度,f为科里奥利参数，h′ 代表海面高度异常值，u和v分别为纬度方向和经度方向的速度。公式(2)中，W代表OW 参数，根据其值的不同，可以将流场分为：涡流主导区域(W<-W0),变形主导区域 (W>W0) 和背景场(|W|≤W0)，其中W0=0.2σW，σW为所有点的OW参数值的标准差。有时为了简便，研究者将W0设定为常数[12]：W0=-2×10-12s-2。这种分区方法已经被很多研究者应用过，证明了其可靠性[5, 13-14]。以上提到的3个区域中，自动动态特征检测方法主要针对的是涡流主导区域。

1.2OW点聚类

在得到OW点以后，需要对搜索区域内的OW点进行聚类分析，根据这些点之间的连接性质可以对其进行聚类。搜索区域的确定方法是比较自由的，具体方法将在下一节中进行介绍。点和点之间的连接性考虑的是4连通性质，即只有1个点的上下左右4个方位出现同样性质的点才认为2个点是连通的。在聚类之前有一些预处理操作需要进行，首先注意到一点：在涡流区域内有一些点不满足OW判据，从而被判定为非OW点，而这样的点如果完全被OW点包围，可以认为其是一些噪声点，这里需要将这些点同样标记为OW点；同时2个OW点集合之间有可能存在模糊连接，例如，2个OW点集合之间只有1个点使其相互连接，这时可以将这个点标记为非OW点从而消除这种模糊连接。

假定在消除噪声点和模糊连接以后，期望搜索的区域内存在n个OW点，这里使用(x1,x2,…,xp…,xn),p=1,2,…,n表示。通过对其连接性的分析，可以将其聚类为k个集合，用S={S1,S2,…,Si,…,Sk},i=1,2,…,k表示。如图 1中左图所示，用“1”表示的点即为OW点，非OW点使用“0”表示；右图为聚类以后的示意图，其中“1”，“2”，“3”表示聚类后的集合编号，具有同样编号的点被认为是属于同一个点集合，即属于同一个涡流区域。

通过对观测结果的归纳分析得出，沿箱梁高度方向分布的温度和以上十分接近，但和现行规范要求的“按箱梁顶板与其它部分的温度为±5℃”有很大差异；沿桥梁长度方向分布的实际温度大致相同。

图1　OW点聚类示意图Fig.1　Schematic diagram of OW points clustering

1.3中尺度涡边界曲线提取

(3)

通过SVM方法，可以得到边界曲线。有了曲线上各个点的坐标，可以计算对应涡流区域的中心位置及其他信息。涡流区域的中心定义为区域的重心，当然，读者也可以定义其他形式的涡流区域中心。

2　涡流演化关系识别

由于涡流的运动本质上是水体的运动，所以本文中将涡流的运动简化为惯性体的运动。对于惯性体运动常用的状态转移方程和观测方程如式(4)所示，其中变量如式(5)所示。Kalman滤波是对线性最小均方误差滤波的另一种处理方法，实际是维纳滤波的一种递推算法。它的工作原理主要是利用协方差矩阵和系统观测方程的相关参数，计算Kalman增益，然后通过Kalman增益值对状态向量预测值进行修正从而得到状态向量的最终估计值。所以，在得到新的SLA数据之前，可以利用Kalman滤波器预测感兴趣的涡流中心的运动情况。Kalman滤波中的状态向量是2维位置向量和2维速度向量的复合向量。Kalman递推公式如式(6)所示，其中K为Kalman增益，X为状态向量，P为协方差矩阵，下角标t|t-1代表其值为预测值，下角标t|t代表其值为估计值。为了保持符号的一致性，接下来的段落中上角标或者下角标中将使用字母“t”表示SLA数据对应的时刻，使用“i”表示在同一时刻SLA数据中聚类得到的不同涡流区域的编号。

(4)

(5)

(6)

除了区域中心的位置距离以外，涡流区域的面积变化和形状变化也作为演化判别的重要参数。得到区域边界曲线以后，将每个涡流区域与上一时刻确定的感兴趣的涡流区域的中心进行重合。假设v(s)=(x(s),y(s)),s∈[0,1] 表示边界曲线。在不引起歧义的情况下，以下的段落将使用v 表示v(s)。即可以使用vt,vt+1,i表示t时刻确定的感兴趣的涡流区域的边界和在t+1 时刻的SLA数据中发现的第i个涡流区域。通过平移操作，可以将两个区域的中心移动到同一个位置，即坐标原点。

平移vt后，得到新的边界曲线表达式为：

(7)

(8)

图2　不重合面积示意图Fig.2　Schematic diagram of noncoincidence area

边界曲线的长度和扭曲程度，这2个参数可以用来量化曲线的形状特征。这2个参数在SNAKE(也称：activecontour)方法中[16-17]已经有成功的应用，但是，对这2个参数的应用方法需要在SNAKE方法的基础上进行改进。因为，长度和扭曲程度是2个不同的特征，将这2个参数分别进行比较，可以更加全面地反映曲线形状的差别。由这2个参数决定的消费参数可以由以下公式计算得到，其中ω1(s)和ω2(s)为两参数的权重参数：

(9)

最终，将涡流区域中心移动项、涡流面积项、涡流边界曲线形状项的作用进行综合可以得到公式(10)所示的演化关系判别公式。

(10)

通过计算演化判别公式 (10), 将搜索区域内综合消费参数最小的涡流区域确定为本次SLA数据中的感兴趣涡流区，即认为上一时刻确定的感兴趣的涡流区域在本次SLA数据采集时演化为此涡流区域，其过程如式(11)所示：

(11)

3　仿真结果

仿真过程中，通过对中国南海区域2014年5月2日-2014年5月31日的数据进行分析，选取一个感兴趣的涡流区域。感兴趣的涡流选择完毕后，使用公式 (10) 所示的判据进行演化关系的判别。仿真结果如图3所示，其中海洋中的蓝色点代表的是OW点。

因为涡流识别中最主要的消费参数是涡流中心的匹配程度所以将公式(10)中的参数设置为：ωA=ωS=ω1(s)=ω2(s)=1，ωD=10，以此来提高中心匹配度的影响因子，而弱化其他影响因子。经过GRBF核函数中σ 参数设定为1时可以较好地对涡流边缘进行平滑。状态转移公式中的输入项U设定为0，即认为涡流中心有保持匀速运动的趋势，这样的假设也符合一般惯性体系的运动性质。鉴于协方差矩阵的初始值P1|0对滤波过程的影响不大，所以将其设定为单位矩阵。公式 (4) 和 (5) 中的其他项，可以依照式(12)所示的进行设置，仿真结果如图3～图6所示。

(12)

图3　涡流区域跟踪Fig.3　Tracking the eddy region

选定的涡流区域的中心运动轨迹如图4所示。涡流区域中心运动的速度幅值大小如图5所示。涡流区域的面积变化情况如图6所示，面积单位为cells,即涡流所占网格的个数，网格的尺寸为1/4°×1/4°(以当地经纬度为基准)。

通过分析仿真结果，可以很明显地看出，该算法可以有效地辨识涡流的动态演化过程。在观测者确定了感兴趣的涡流区域以后，本文中的算法能够自动识别涡流区域并跟踪其动态变化。仿真过程中从涡流的产生到消散，没有出现目标丢失和跟踪错误的现象。即使在5月13日—5月18日期间，涡流中心位置、区域面积和形状急剧变化的情况下，算法依然能够成功跟踪涡流区域，证明算法具有较强的鲁棒性。

图4　涡流区域中心运动轨迹Fig.4　Trace of the eddy center

图5　涡流中心运动速度幅值Fig.5　Velocity magnitude of the eddy center movement

图6　涡流区域面积变化Fig.6　Variation of the eddy region area

4　小结

本文提出的算法本身可以提取涡流区域的光滑边界，并且最终根据演化判别公式的计算结果得到涡流区域的演化关系。通过仿真结果可知，相关判据和算法能够自动识别涡流的演化关系，实现涡流动态信息的识别。值得说明的是，由于WA(Wind-AngleMethod) 涡流区域判别标准和基于等高线的涡流区域判别标准能够直接给出涡流区域的边界曲线，算法和演化判别公式可以直接进行应用而不用进行聚类和边界曲线提取的过程。

对于涡流区域的分离和融合过程，需要进行深入的研究才能进行自动识别。同时文中的算法和判别标准可以进行改进，从而进行多涡流区域的跟踪。未来的工作中，还将考虑AUVs的控制算法，从而实现涡流区域实时跟踪自主采样算法。

致谢仿真所用的海面高度数据由Ssalto/Duacs提供，通过访问Cnes支持下的Aviso网站进行下载。同时感谢MGET工具箱，本文算法在制定过程中受到MGET工具箱的启发。

[1]ZHANGZheng-guang,WANGWei,QIUBo.Oceanicmasstransportbymesoscaleeddies[J].Science,2014,345(6194):322-324.

[2]MARTINJP,LEECM,ERIKSENCC,etal.GliderobservationsofkinematicsinaGulfofAlaskaeddy[J].JGeophysRes-Oceans,2009,114(C12):43-47.

[3]PENVENP,ECHEVINV,PASAPERAJ,etal.Averagecirculation,seasonalcycle,andmesoscaledynamicsofthePeruCurrentSystem:Amodelingapproach[J].JournalofGeophysicalResearchAtmospheres,2005,110(C10):901-902.

[4]CHAIGNEAUA,GIZOLMEA,GRADOSC.MesoscaleeddiesoffPeruinaltimeterrecords:Identificationalgorithmsandeddyspatio-temporalpatterns[J].ProgressinOceanography,2008,79(2-4):106-119.

[5]ISERN-FONTANETJ,GARCIA-LADONAE,FONTJ.Identificationofmarineeddiesfromaltimetricmaps[J].JAtmosOceanTech,2003,20(5):772-778.

[6]ISERN-FONTANETJ,GARCA-LADONAE,FONTJ.VorticesoftheMediterraneanSea:Analtimetricperspective[J].JournalofPhysicalOceanography,2006,36(1):87-103.

[7]CHELTONDB,SCHLAXMG,SAMELSONRM,etal.Globalobservationsoflargeoceaniceddies[J].GeophysicalResearchLetters,2007,34(15):87-101.

[8]CHELTONDB,SCHLAXMG,SAMELSONRM.Globalobservationsofnonlinearmesoscaleeddies[J].ProgressinOceanography,2011,91(2):167-216.

[9]ROBERTSJJ,BESTBD,DUNNDC,etal.Marinegeospatialecologytools:AnintegratedframeworkforecologicalgeoprocessingwithArcGIS,Python,R,MATLAB,andC++ [J].EnvironmentalModelling&Software,2010,25(10):1 197-1 207.

[10]BERON-VERAFJ,OLASCOAGAMJ,GONIGJ.OceanicmesoscaleeddiesasrevealedbyLagrangiancoherentstructures[J].GeophysicalResearchLetters,2008,35(12):86-109.

[11]ISERN-FONTANETJ,FONTJ,GARCIA-LADONAE,etal.SpatialstructureofanticycloniceddiesintheAlgerianbasin(MediterraneanSea)analyzedusingtheOkubo-Weissparameter[J].Deep-SeaResPtIi,2004,51(25-26):3 009-3 028.

[12]SOUZAJMAC,DEBOYERMONTGUTC,LETRAONPY.ComparisonbetweenthreeimplementationsofautomaticidentificationalgorithmsforthequantificationandcharacterizationofmesoscaleeddiesintheSouthAtlanticOcean[J].OceanScience,2011,7(3):317-334.

[13]JEONGJ,HUSSAINF.Ontheidentificationofavortex[J].JournalofFluidMechanics,2006,285:69-94.

[14]PASQUEROC,PROVENZALEA,BABIANOA.Parameterizationofdispersionintwo-dimensionalturbulence[J].JFluidMech,2001,439(5)：279-303.

[15]HSUC-W,CHANGC-C,LINC-J.Apracticalguidetosupportvectorclassification[J]. 臺北市:國立臺灣大學資訊工程學系,2003,67(5):1-12.

[16]KASSM,WITKINA,TERZOPOULOSD.Snakes:Activecontourmodels[J].InternationalJournalofComputerVision,1988,1(4):321-331.

[17]TERZOPOULOSD,SZELISKIR.Activevision[M].Cambridge：MITPressCambridge，1993:3-20.

Dynamic feature detection of mesoscale eddies based on SLA data

ZHAO Wen-tao1,2, YU Jian-cheng*1, ZHANG Ai-qun1, LI Yan1

(1.StateKeyLaboratoryofRobotics,ShenyangInstituteofAutomation,ChineseAcademyofSciences,Shenyang110016,China; 2.UniversityofChineseAcademyofSciences,Beijing100049,China)

To automatically sample mesoscale eddies by AUVs, the method to automatically recognize the dynamic features of eddies must be developed. In this study, a method for dynamic feature detection of mesoscale eddies was created based on SLA (Sea Level Anomaly) data. The main innovation is that a criterion to decide the succession relations of eddies was developed. With the succession relations of eddies being confirmed, we can calculate the area changing rate of eddy region, velocity of eddy centroid movement and some other dynamic features. Since the calculation cost of this algorithm is not enormous, it can be used in real time eddy tracking context, such as tracking eddies with AUVs. It is worth noting that SSH(Sea Surface Height) data can also be used in this algorithm. And SLA or SSH data can be acquired from the ocean numerical simulation model, satellite remote sensing or other methods. As long as SLA or SSH data are provided, our algorithm can be easily utilized for dynamic feature detection of mesoscale eddies.

mesoscale eddy; dynamic feature; succession relationship; automatically detection

2016-01-18

2016-04-03

国家自然科学基金项目资助(61233013)

赵文涛(1989-)，男，山东淄博市人，博士研究生，主要从事海洋中尺度涡流跟踪方面的研究。E-mail:zwt2002@gmail.com

俞建成(1974-)，男，研究员，主要从事水下机器人方面的研究。E-mail:yjc@sia.cn

P731.2

A

1001-909X(2016)03-0062-07

10.3969/j.issn.1001-909X.2016.03.010

赵文涛，俞建成，张艾群，等.基于卫星测高数据的海洋中尺度涡流动态特征检测[J].海洋学研究,2016,34(3):62-68,doi:10.3969/j.issn.1001-909X.2016.03.010.

ZHAO Wen-tao, YU Jian-cheng, ZHANG Ai-qun, et al. Dynamic feature detection of mesoscale eddies based on SLA data[J]. Journal of Marine Sciences, 2016,34(3):62-68, doi:10.3969/j.issn.1001-909X.2016.03.010.