基于衛(wèi)星測高數(shù)據(jù)的海洋中尺度渦流動態(tài)特征檢測

趙文濤，俞建成，張艾群，李　巖

(1. 中國科學(xué)院 沈陽自動化研究所, 遼寧 沈陽 110016； 2. 中國科學(xué)院大學(xué), 北京 100049)

?

基于衛(wèi)星測高數(shù)據(jù)的海洋中尺度渦流動態(tài)特征檢測

趙文濤1,2，俞建成*1，張艾群1，李巖1

(1. 中國科學(xué)院 沈陽自動化研究所, 遼寧 沈陽 110016； 2. 中國科學(xué)院大學(xué), 北京 100049)

為了最終實(shí)現(xiàn)對海洋中尺度渦流(簡稱中尺度渦)的自動采樣，首先應(yīng)該發(fā)展中尺度渦動態(tài)特征識別技術(shù)。本文基于SLA(Sea Level Anomaly)數(shù)據(jù)，實(shí)現(xiàn)了對中尺度渦動態(tài)特征的檢測算法。主要內(nèi)容是制定了一個判別相鄰兩組SLA數(shù)據(jù)中的渦流，是否為同一渦流子在不同時刻的狀態(tài)的標(biāo)準(zhǔn)，即判別下一時刻SLA數(shù)據(jù)中是否存在渦流是由上一時刻確定的被檢測渦流演化而來的。通過確定這種進(jìn)化關(guān)系，可以得到被檢測渦流的一系列動態(tài)狀態(tài)信息，例如：面積變化速率、中心移動情況以及其他情況。本算法的計算量不大，從而可以應(yīng)用到實(shí)時渦流跟蹤的環(huán)境中。值得注意的是，本文中的算法不僅僅局限于應(yīng)用SLA數(shù)據(jù)，SSH(Sea Surface Height)等大部分反映海洋高度的數(shù)據(jù)也可以使用。

中尺度渦；動態(tài)特征；演化關(guān)系；自動檢測

0　引言

海洋熱量和物質(zhì)輸送過程對全球氣候有深遠(yuǎn)影響[1]。中尺度渦對于物質(zhì)的輸送作用幾乎可以和風(fēng)生流和熱生流相媲美。中尺度渦流區(qū)域在海洋活動中的重要影響，使科學(xué)家們對其產(chǎn)生了濃厚的興趣。隨著集成電路技術(shù)的發(fā)展，利用AUV等自主平臺進(jìn)行長時間海洋特征觀測已經(jīng)成為可能。雖然已經(jīng)有相關(guān)研究人員對中尺度渦進(jìn)行采樣調(diào)查[2]。但是，到目前為止還沒有很好的方法可以對中尺度渦實(shí)現(xiàn)自主采樣，不能夠準(zhǔn)確地采集敏感區(qū)域的特征數(shù)據(jù)，反應(yīng)渦流區(qū)域的內(nèi)部特征。

掌握中尺度渦流區(qū)域的移動路徑、面積變化、形狀變化等特征可以為研究者提供必要的信息，從而為AUV采樣策略的制定提供依據(jù)，最終提高對中尺度渦進(jìn)行跟蹤觀測的準(zhǔn)確性，使采集的數(shù)據(jù)更加有效地反映渦流區(qū)域的內(nèi)部特征。鑒于衛(wèi)星數(shù)據(jù)可以有效地反映渦流區(qū)域的宏觀動態(tài)特征，利用衛(wèi)星數(shù)據(jù)自主識別中尺度渦流區(qū)域的移動路徑、面積變化、形狀變化等動態(tài)特征的方法能夠有效減輕研究者的工作強(qiáng)度，提高效率。

PENVEN et al[3]提出了一種渦流自動檢測的方法。但是他們將渦流區(qū)域簡單地作為一個圓形區(qū)域?qū)Υ?，這種方式忽略了渦流形狀的信息。同時在他們的算法中使用了漢寧濾波器，這使初生期的渦流很難被發(fā)現(xiàn)，因?yàn)檫@個時候的渦流區(qū)域還很小，特征表現(xiàn)還不明顯。在他們的算法中沒有考慮渦流的移動速度對動態(tài)特征識別的影響，所以當(dāng)多個渦流距離很近的時候，很容易發(fā)生進(jìn)化關(guān)系檢測的錯誤。雖然CHAIGNEAU et al[4]將PENVEN et al[3]的算法進(jìn)行了改進(jìn)，將EKE(Eddy Kinetic Energy) 加入到進(jìn)化關(guān)系判別的標(biāo)準(zhǔn)中，但是其沒有針對算法的確定進(jìn)行改進(jìn)。

基于ISERN-FORTANET et al[5-6]的算法, CHELTON et al[7-8]提出了一種渦流自主檢測的算法。 但是他們只針對這種方法進(jìn)行了簡單的說明，沒有給出具體的算法和表達(dá)式。同時其算法中也缺少對于中尺度渦進(jìn)行關(guān)系判別的說明。在 MGET 工具箱中[9], 渦流之間的進(jìn)化關(guān)系主要是依靠兩組海面測高數(shù)據(jù)中渦流區(qū)域的重復(fù)情況進(jìn)行判別的。這種方法沒有考慮渦流形狀的變化等信息。如果兩個渦流距離很近，也會引起進(jìn)化關(guān)系判別的錯誤。

對于渦流區(qū)域的判別，有很多方法可以做到。BERON-VERA et al[10]提出一種客觀的(即坐標(biāo)系獨(dú)立的)方法可以用于進(jìn)行渦流邊界的識別。ISERN-FORTANET et al[5]提出一種根據(jù)高度信息進(jìn)行渦流中心區(qū)域檢測的方法。OW(Okubo-Weiss)參數(shù)可以用來進(jìn)行SLA(Satellite Level Anomaly)數(shù)據(jù)中渦流區(qū)域的檢測，ISERN-FORTANET et al[11]給出了OW參數(shù)計算的過程，以及過程中所需的速度場的計算公式。SLA數(shù)據(jù)中滿足OW渦流主導(dǎo)區(qū)域判別標(biāo)準(zhǔn)的點(diǎn)成為OW點(diǎn)，不滿足的則稱為非OW點(diǎn)。

本研究將渦流中心的運(yùn)動速度、渦流區(qū)域的形狀變化、面積變化這3個因素作為渦流進(jìn)化關(guān)系判別方法中的主要參數(shù)，對渦流區(qū)域檢測進(jìn)行了簡單的說明，對渦流進(jìn)化關(guān)系識別進(jìn)行了具體說明，最后給出了基于SLA數(shù)據(jù)的中尺度渦動態(tài)特征識別的仿真結(jié)果。

1　渦流區(qū)域檢測

1.1OW參數(shù)

ISERN-FORTANET等人給出了如下所示的OW參數(shù)計算公式[11]。

2D流場的計算公式：

(1)

OW參數(shù)的計算公式：

(2)

公式(1)中g(shù)代表重力加速度,f為科里奧利參數(shù)，h′ 代表海面高度異常值，u和v分別為緯度方向和經(jīng)度方向的速度。公式(2)中，W代表OW 參數(shù)，根據(jù)其值的不同，可以將流場分為：渦流主導(dǎo)區(qū)域(W<-W0),變形主導(dǎo)區(qū)域 (W>W0) 和背景場(|W|≤W0)，其中W0=0.2σW，σW為所有點(diǎn)的OW參數(shù)值的標(biāo)準(zhǔn)差。有時為了簡便，研究者將W0設(shè)定為常數(shù)[12]：W0=-2×10-12s-2。這種分區(qū)方法已經(jīng)被很多研究者應(yīng)用過，證明了其可靠性[5, 13-14]。以上提到的3個區(qū)域中，自動動態(tài)特征檢測方法主要針對的是渦流主導(dǎo)區(qū)域。

1.2OW點(diǎn)聚類

在得到OW點(diǎn)以后，需要對搜索區(qū)域內(nèi)的OW點(diǎn)進(jìn)行聚類分析，根據(jù)這些點(diǎn)之間的連接性質(zhì)可以對其進(jìn)行聚類。搜索區(qū)域的確定方法是比較自由的，具體方法將在下一節(jié)中進(jìn)行介紹。點(diǎn)和點(diǎn)之間的連接性考慮的是4連通性質(zhì)，即只有1個點(diǎn)的上下左右4個方位出現(xiàn)同樣性質(zhì)的點(diǎn)才認(rèn)為2個點(diǎn)是連通的。在聚類之前有一些預(yù)處理操作需要進(jìn)行，首先注意到一點(diǎn)：在渦流區(qū)域內(nèi)有一些點(diǎn)不滿足OW判據(jù)，從而被判定為非OW點(diǎn)，而這樣的點(diǎn)如果完全被OW點(diǎn)包圍，可以認(rèn)為其是一些噪聲點(diǎn)，這里需要將這些點(diǎn)同樣標(biāo)記為OW點(diǎn)；同時2個OW點(diǎn)集合之間有可能存在模糊連接，例如，2個OW點(diǎn)集合之間只有1個點(diǎn)使其相互連接，這時可以將這個點(diǎn)標(biāo)記為非OW點(diǎn)從而消除這種模糊連接。

假定在消除噪聲點(diǎn)和模糊連接以后，期望搜索的區(qū)域內(nèi)存在n個OW點(diǎn)，這里使用(x1,x2,…,xp…,xn),p=1,2,…,n表示。通過對其連接性的分析，可以將其聚類為k個集合，用S={S1,S2,…,Si,…,Sk},i=1,2,…,k表示。如圖 1中左圖所示，用“1”表示的點(diǎn)即為OW點(diǎn)，非OW點(diǎn)使用“0”表示；右圖為聚類以后的示意圖，其中“1”，“2”，“3”表示聚類后的集合編號，具有同樣編號的點(diǎn)被認(rèn)為是屬于同一個點(diǎn)集合，即屬于同一個渦流區(qū)域。

通過對觀測結(jié)果的歸納分析得出，沿箱梁高度方向分布的溫度和以上十分接近，但和現(xiàn)行規(guī)范要求的“按箱梁頂板與其它部分的溫度為±5℃”有很大差異；沿橋梁長度方向分布的實(shí)際溫度大致相同。

圖1　OW點(diǎn)聚類示意圖Fig.1　Schematic diagram of OW points clustering

1.3中尺度渦邊界曲線提取

(3)

通過SVM方法，可以得到邊界曲線。有了曲線上各個點(diǎn)的坐標(biāo)，可以計算對應(yīng)渦流區(qū)域的中心位置及其他信息。渦流區(qū)域的中心定義為區(qū)域的重心，當(dāng)然，讀者也可以定義其他形式的渦流區(qū)域中心。

2　渦流演化關(guān)系識別

由于渦流的運(yùn)動本質(zhì)上是水體的運(yùn)動，所以本文中將渦流的運(yùn)動簡化為慣性體的運(yùn)動。對于慣性體運(yùn)動常用的狀態(tài)轉(zhuǎn)移方程和觀測方程如式(4)所示，其中變量如式(5)所示。Kalman濾波是對線性最小均方誤差濾波的另一種處理方法，實(shí)際是維納濾波的一種遞推算法。它的工作原理主要是利用協(xié)方差矩陣和系統(tǒng)觀測方程的相關(guān)參數(shù)，計算Kalman增益，然后通過Kalman增益值對狀態(tài)向量預(yù)測值進(jìn)行修正從而得到狀態(tài)向量的最終估計值。所以，在得到新的SLA數(shù)據(jù)之前，可以利用Kalman濾波器預(yù)測感興趣的渦流中心的運(yùn)動情況。Kalman濾波中的狀態(tài)向量是2維位置向量和2維速度向量的復(fù)合向量。Kalman遞推公式如式(6)所示，其中K為Kalman增益，X為狀態(tài)向量，P為協(xié)方差矩陣，下角標(biāo)t|t-1代表其值為預(yù)測值，下角標(biāo)t|t代表其值為估計值。為了保持符號的一致性，接下來的段落中上角標(biāo)或者下角標(biāo)中將使用字母“t”表示SLA數(shù)據(jù)對應(yīng)的時刻，使用“i”表示在同一時刻SLA數(shù)據(jù)中聚類得到的不同渦流區(qū)域的編號。

(4)

(5)

(6)

除了區(qū)域中心的位置距離以外，渦流區(qū)域的面積變化和形狀變化也作為演化判別的重要參數(shù)。得到區(qū)域邊界曲線以后，將每個渦流區(qū)域與上一時刻確定的感興趣的渦流區(qū)域的中心進(jìn)行重合。假設(shè)v(s)=(x(s),y(s)),s∈[0,1] 表示邊界曲線。在不引起歧義的情況下，以下的段落將使用v 表示v(s)。即可以使用vt,vt+1,i表示t時刻確定的感興趣的渦流區(qū)域的邊界和在t+1 時刻的SLA數(shù)據(jù)中發(fā)現(xiàn)的第i個渦流區(qū)域。通過平移操作，可以將兩個區(qū)域的中心移動到同一個位置，即坐標(biāo)原點(diǎn)。

平移vt后，得到新的邊界曲線表達(dá)式為：

(7)

(8)

圖2　不重合面積示意圖Fig.2　Schematic diagram of noncoincidence area

邊界曲線的長度和扭曲程度，這2個參數(shù)可以用來量化曲線的形狀特征。這2個參數(shù)在SNAKE(也稱：activecontour)方法中[16-17]已經(jīng)有成功的應(yīng)用，但是，對這2個參數(shù)的應(yīng)用方法需要在SNAKE方法的基礎(chǔ)上進(jìn)行改進(jìn)。因?yàn)?，長度和扭曲程度是2個不同的特征，將這2個參數(shù)分別進(jìn)行比較，可以更加全面地反映曲線形狀的差別。由這2個參數(shù)決定的消費(fèi)參數(shù)可以由以下公式計算得到，其中ω1(s)和ω2(s)為兩參數(shù)的權(quán)重參數(shù)：

(9)

最終，將渦流區(qū)域中心移動項(xiàng)、渦流面積項(xiàng)、渦流邊界曲線形狀項(xiàng)的作用進(jìn)行綜合可以得到公式(10)所示的演化關(guān)系判別公式。

(10)

通過計算演化判別公式 (10), 將搜索區(qū)域內(nèi)綜合消費(fèi)參數(shù)最小的渦流區(qū)域確定為本次SLA數(shù)據(jù)中的感興趣渦流區(qū)，即認(rèn)為上一時刻確定的感興趣的渦流區(qū)域在本次SLA數(shù)據(jù)采集時演化為此渦流區(qū)域，其過程如式(11)所示：

(11)

3　仿真結(jié)果

仿真過程中，通過對中國南海區(qū)域2014年5月2日-2014年5月31日的數(shù)據(jù)進(jìn)行分析，選取一個感興趣的渦流區(qū)域。感興趣的渦流選擇完畢后，使用公式 (10) 所示的判據(jù)進(jìn)行演化關(guān)系的判別。仿真結(jié)果如圖3所示，其中海洋中的藍(lán)色點(diǎn)代表的是OW點(diǎn)。

因?yàn)闇u流識別中最主要的消費(fèi)參數(shù)是渦流中心的匹配程度所以將公式(10)中的參數(shù)設(shè)置為：ωA=ωS=ω1(s)=ω2(s)=1，ωD=10，以此來提高中心匹配度的影響因子，而弱化其他影響因子。經(jīng)過GRBF核函數(shù)中σ 參數(shù)設(shè)定為1時可以較好地對渦流邊緣進(jìn)行平滑。狀態(tài)轉(zhuǎn)移公式中的輸入項(xiàng)U設(shè)定為0，即認(rèn)為渦流中心有保持勻速運(yùn)動的趨勢，這樣的假設(shè)也符合一般慣性體系的運(yùn)動性質(zhì)。鑒于協(xié)方差矩陣的初始值P1|0對濾波過程的影響不大，所以將其設(shè)定為單位矩陣。公式 (4) 和 (5) 中的其他項(xiàng)，可以依照式(12)所示的進(jìn)行設(shè)置，仿真結(jié)果如圖3～圖6所示。

(12)

圖3　渦流區(qū)域跟蹤Fig.3　Tracking the eddy region

選定的渦流區(qū)域的中心運(yùn)動軌跡如圖4所示。渦流區(qū)域中心運(yùn)動的速度幅值大小如圖5所示。渦流區(qū)域的面積變化情況如圖6所示，面積單位為cells,即渦流所占網(wǎng)格的個數(shù)，網(wǎng)格的尺寸為1/4°×1/4°(以當(dāng)?shù)亟?jīng)緯度為基準(zhǔn))。

通過分析仿真結(jié)果，可以很明顯地看出，該算法可以有效地辨識渦流的動態(tài)演化過程。在觀測者確定了感興趣的渦流區(qū)域以后，本文中的算法能夠自動識別渦流區(qū)域并跟蹤其動態(tài)變化。仿真過程中從渦流的產(chǎn)生到消散，沒有出現(xiàn)目標(biāo)丟失和跟蹤錯誤的現(xiàn)象。即使在5月13日—5月18日期間，渦流中心位置、區(qū)域面積和形狀急劇變化的情況下，算法依然能夠成功跟蹤渦流區(qū)域，證明算法具有較強(qiáng)的魯棒性。

圖4　渦流區(qū)域中心運(yùn)動軌跡Fig.4　Trace of the eddy center

圖5　渦流中心運(yùn)動速度幅值Fig.5　Velocity magnitude of the eddy center movement

圖6　渦流區(qū)域面積變化Fig.6　Variation of the eddy region area

4　小結(jié)

本文提出的算法本身可以提取渦流區(qū)域的光滑邊界，并且最終根據(jù)演化判別公式的計算結(jié)果得到渦流區(qū)域的演化關(guān)系。通過仿真結(jié)果可知，相關(guān)判據(jù)和算法能夠自動識別渦流的演化關(guān)系，實(shí)現(xiàn)渦流動態(tài)信息的識別。值得說明的是，由于WA(Wind-AngleMethod) 渦流區(qū)域判別標(biāo)準(zhǔn)和基于等高線的渦流區(qū)域判別標(biāo)準(zhǔn)能夠直接給出渦流區(qū)域的邊界曲線，算法和演化判別公式可以直接進(jìn)行應(yīng)用而不用進(jìn)行聚類和邊界曲線提取的過程。

對于渦流區(qū)域的分離和融合過程，需要進(jìn)行深入的研究才能進(jìn)行自動識別。同時文中的算法和判別標(biāo)準(zhǔn)可以進(jìn)行改進(jìn)，從而進(jìn)行多渦流區(qū)域的跟蹤。未來的工作中，還將考慮AUVs的控制算法，從而實(shí)現(xiàn)渦流區(qū)域?qū)崟r跟蹤自主采樣算法。

致謝仿真所用的海面高度數(shù)據(jù)由Ssalto/Duacs提供，通過訪問Cnes支持下的Aviso網(wǎng)站進(jìn)行下載。同時感謝MGET工具箱，本文算法在制定過程中受到MGET工具箱的啟發(fā)。

[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]. 臺北市:國立臺灣大學(xué)資訊工程學(xué)系,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

國家自然科學(xué)基金項(xiàng)目資助(61233013)

趙文濤(1989-)，男，山東淄博市人，博士研究生，主要從事海洋中尺度渦流跟蹤方面的研究。E-mail:zwt2002@gmail.com

俞建成(1974-)，男，研究員，主要從事水下機(jī)器人方面的研究。E-mail:yjc@sia.cn

P731.2

A

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

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

趙文濤，俞建成，張艾群，等.基于衛(wèi)星測高數(shù)據(jù)的海洋中尺度渦流動態(tài)特征檢測[J].海洋學(xué)研究,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.