數據驅動印度洋海域全局動力學研究1)

李自剛 嚴 旺 康佳琪 江 俊 洪 靈

* (西安科技大學力學系,西安 710054)

? (西安交通大學機械結構強度與振動國家重點實驗室,西安 710049)

引言

海洋系統天然受到多種環境因素(風速、溫度、海水含鹽量等)、人類活動和天體運動等不確定擾動的影響,導致其洋流運動復雜而多樣[1-2].然而,一些觀測和研究發現,其內部蘊含著能夠對洋流的運動方向和趨勢起限制和導向作用的復雜動力結構和特征,影響甚至主導了表面漂浮物的流動形式和演化路徑.例如,經志友等[3-4]觀測發現南海內部存在自循環環流結構,以及跨海域間的貫穿流結構.這種暗藏結構在動力學上表現為對周圍流動體的吸引或排斥性效應,并能夠引起流場的渦漩,進而可能導致表面漂浮物在某些水域的長期或短期聚集.Serra 等[5]驗證了海洋表面暫態吸引性的存在,并展示了這種內在特征結構在海洋救援中的積極作用.由此,需全面系統地揭示這類系統復雜的內在結構狀態和特性,這是深入認識并利用海洋規律的關鍵.

海洋作為典型的強非線性隨機系統,對其洋流動態演化和動力結構特性的研究,主要從數值模擬和解析分析兩方面開展.前者大多以水動力經典方程(如:Langevin 方程、Fokker?Planck 方程、Navier?Stokes 方法和對流擴散方程等)為基礎[6-7],能利用數值算法針對復雜情況開展模擬,但通常計算難度較大.而傳統對簡化模型的解析分析簡單高效,但僅能從機理上研究特定條件下的局部定性行為(如海洋吸引結構的存在條件[8]和概率躍遷機理[9]),卻無法準確地描述耦合因素的影響和實際特征,且計算結果依賴于模型簡化的合理性和準確性.遙感大數據的支撐大大提高了解析法的預測精度,使其在工程應用上具有更廣闊的應用前景[10-12].如薛紅娟和顧耀林[13]通過對數據模型的拓撲分析,研究了海洋渦旋的動力學特征及其形貌.此外,一些研究將數據的互相關信息[14]、主元分量[15]、嵌入維數[16]等作為刻畫海洋系統動力特征的重要指標,并結合實驗浮漂軌跡及其統計分布[17],利用貝葉斯估計[18]、機器學習[19]等數據驅動方法,開海洋展表面漂浮物定位和軌跡預測.海洋系統的內在動力結構十分復雜且主導著其動力學響應和趨勢,但現有的方法大多將數據看作是離散的軌跡信息,僅反映了時間序列的非線性關系或數據的概率特征,缺乏對系統動態演化行為與內在結構屬性間聯系的深入研究,亟待進一步以數據為基礎從全局動力學角度揭示其內在特征及其影響機理.

本文將數據驅動廣義胞映射方法應用于海洋動力學的分析中,通過建立浮漂運動狀態間的映射關系,形成表征海洋洋流運動特征的一步轉移概率矩陣,著重分析印度洋洋流的長期和短期全局動力結構和演化路徑,以期為該區域的研究海洋運動規律和污染物擴散、優化船舶航線、實施海洋救援等提供指導和依據,也致力于推動胞映射方法在工程領域的應用和發展.

1 數據來源和分析方法

1.1 數據來源

文中以實際測量數據為研究基礎,重點關注于印度洋附近海域(東經15°～ 120°,南緯60°～ 北緯30°)洋流的全局動態特征及其引發的響應規律.數據來源于美國國家海洋與大氣管理局過去30 年(1979—2019 年)實施的全球漂流計劃項目(Global Drifter Program).該項目在全球海洋表面共部署了超過2.4 萬個電子浮漂,并通過衛星每6 小時監測浮漂器的狀態數據.目前,初始部署在印度洋附近海域的浮漂數目有3680 個,而與該區域直接相關的浮漂則達到5006 個,如圖1 所示.

圖1 印度洋海域部署的浮漂位置(黑色“*”)和運動軌跡(彩色實線)Fig.1 The initial locations and trajectories of drifters deployed in the Indian Ocean

通過對采集的時域數據預處理后,數據結構可表示為如下形式

式中,Yr表示第r個浮漂在整個服役周期內以采樣頻率fs(d?1)獲得的運動軌跡,代表了該浮漂在第i時刻下的運動狀態,Nr為連續采集的數據點個數,Nd則表示與該區域相關的浮漂總數目.

1.2 數據驅動模型和分析方法

胞映射方法是開展復雜系統全局動力學分析的有效手段[20-24].該方法的基本思想是將連續的D維狀態空間RD離散為一個在給定尺度下的胞狀態空間,且覆蓋感興趣區域的每個胞通過Nc個整數連續編號.通過識別數據中所攜帶的動力學信息(如周期、主維數等),并將其轉變為離散空間中各胞之間的有向映射.此時,隨機擾動下原動力系統的概率演化可通過胞空間中Markov 過程描述[25]

式中,p(k)為第k步下胞的概率向量,其分量pi(k)則代表在該步映射下第i個胞的狀態概率.P(k)表示動力系統的一步轉移概率矩陣,其元素Pij表示從第j個胞到第i個像胞的轉移概率,可通過如下公式計算

式中,τ=t?t0表示一個映射步長,Ci為在空間中RD第i個胞所在的區域,p(x,τ|xj,0)表示從第j個胞映射到像胞的一步轉移概率向量.

在數據驅動胞映射方法中[26],通過選擇一個給定的時間間隔(映射步長τ) 對浮漂軌跡重新采樣,以此實現映射過程,進而從數據中獲得轉移概率矩陣的近似估計.從式(1)中,第r個浮漂映射采樣點可表示為

此時,數據驅動的轉移概率Pij可以通過在胞空間中數值統計從j個胞到第i個胞的樣本概率Sij/S j近似獲得.

真實海洋環境中,海洋表面的洋流演化和變遷與當地海域季節性的氣候變化、大洋環流等因素密切相關[27-29],例如,在北半球冬季,當季風向西南吹時,海洋表面洋流從印尼群島附近向阿拉伯海方向向西流動.在北半球夏季,隨著季風向東北方向的變化,海洋環流逆轉,向東流從索馬里延伸到孟加拉灣.因此可見,轉移概率矩陣P將會受到海洋季節改變的影響,從而表現為P與時間直接相關,即式(2)為非自治Markov 過程.在本文中,按照四季變化規律將一年分為1—3 月(TJM),4—6 月(TAJ),7—9 月(TJS)以及10—12 月(TOD) 4 個時間區間刻畫印度洋區域洋流的典型變化模式.此時,在這些時間區間內,系統的狀態演化可形成4 個相應的分段季度轉移概率矩陣:PJM,PAJ,PJS和POD.

根據胞映射方法,映射步長τ的選取應反映系統動力學的演化特征.海洋系統中,拉格朗日不相關時間(lagrangian decorrelation time)可用于刻畫海洋表面非均勻流體的擴散程度[30-31].本文選取該時間T=3 d 作為一步映射步長τ,以此構建對應不同區間的一步映射矩陣Ps1,Ps2,Ps3和Ps4.同時,考慮每月均為31 d,每季度為3 個月,此時季度一步轉移概率矩陣可通過對式(2)的左乘,表示為

為了從有限數據中獲得轉移概率的精確估計,每個胞中應包含盡量多的樣本軌跡,并確定對應的像胞,形成廣義胞映射.根據圖論理論,Pij>0 代表了從第j個結點到第i個結點的一條有向邊,由此轉移概率矩陣實際上構成了一組以映射關系表征的有向圖.利用深度優先搜索方法對該有向圖中強連通分量的搜索和挖掘,可獲得表征系統吸引性的渦旋結構及其影響區域等全局動力學特征[26],其中渦旋中心通過映射自循環胞刻畫,而渦旋區域則由可達渦旋中心的胞構成.需要說明的是,根據大偏差理論[32],海洋系統中的強隨機擾動特征使得所有樣本響應在足夠長的演化時間下以概率1 離開所在區域而遍歷所有可能空間.因此,本文通過拓撲分析得到的渦旋結構具有動態特性,其概念和性質與傳統意義下的穩態自循環集合不完全一致,在海洋系統中表征在一定時間區間內對周圍軌道(樣本)的導向和吸引,并通過分布概率密度來描述樣本的聚集程度.

2 印度洋海域全局動力學分析

由于不確定因素的影響,海洋系統的長期和短期動力結構具有明顯的區別和聯系.長期動力結構表征了系統的穩態分布趨勢和總體特征,而短期動力結構則關注于在足夠短的時間區間內系統的瞬態運動特征和內在的影響機制,其行為更加復雜且受到噪聲影響較大.下文將從這兩個方面分別開展研究.

2.1 長期動力結構

海洋表面漂浮物在較長時間后的位置分布和動力學狀態與該區域洋流的長期全局動力學結構密切相關.本文中選取印度洋海域R={15°

圖2 印度洋長期海域渦旋結構和渦旋區Fig.2 Long-term vortex structure and its regions of influence in the Indian Ocean

圖3 顯示了由式(2)獲得的該系統隨機響應概率分布及演化特征.圖3 中可見,在該區域內系統的影響由初始的均勻分布開始演化并逐漸遠離赤道附近向西南方向運動,最終經歷1 年的演化時間后,大部分概率動態地聚集在南緯20°～ 40°范圍內,特別是在上述渦旋區范圍內分布概率達到最大,而在南緯40°以南和赤道附近的概率較低.同時,在海洋與陸地交界面附近也表現較高的響應概率,這主要是由于海岸線附近地貌特征和人類活動導致了浮漂運動的“黏滯效應”.

圖3 印度洋海域洋流流動概率分布及動態演化:(a) 演化時間T=3 個月;(b) 演化時間T=6 個月;(c) 演化時間T=9 個月;(d) 演化時間T=12 個月Fig.3 Probabilities and evolution of ocean currents in the region of the Indian Ocean:(a) evolving time T=3 mon;(b) evolving time T=6 mon;(c) evolving time T=9 mon;(d) evolving time T=12 mon

為了驗證所提方法和結果的有效性,本文提取2017—2019 年在印度洋海域的浮漂樣本作為檢驗數據,其中2017 年1 月在該區域共新投入45 個有效浮漂(圖4(a)中“△”表示),且總體成均勻分布.以此為初始狀態,圖4(b)中彩色區域顯示了本文預測的2019 年1 月浮漂響應概率分布,“*”則表示浮漂真實狀態.圖中可見,雖然近海“黏滯效應”限制3 個樣本的自由漂移,但仍有39 個浮漂(約占86.7%)處于圖2 所示的渦旋區周圍,且主要集中在預測概率較高的核心區域(南緯23°～ 37°,東經40°～ 75°范圍),這表明了所提出的數據驅動胞映射方法在海洋洋流預測方面的有效性.

圖4 2017—2019 年浮漂的位置分布:(a) 2017 年1 月浮漂位置(“△”);(b) 2019 年1 月浮漂狀態(“*”)與預測結果對比Fig.4 Distributions of drifters during 2017?2019:(a) states of drifters in January,2017;(b) Comparison of real states and predicted results in January,2019

2.2 短期動力結構

海洋短期動力結構能夠揭示該區域漂浮物在經歷短暫動態演化后的運動狀態和趨勢,在海洋搜索和救援領域具有重要的實用價值.本節將以2014 年3 月8 日MH370 航班失事(雷達消失)事件為時間節點,研究此后1 周內(3 月8 日—14 日) 該海域表面洋流的短期動力學特征.分析時,提取1979—2013 年在此時間窗口內的所有歷史數據,形成以3 d為映射步長的一步短期瞬態轉移概率矩陣Pt,該矩陣表征系統在歷史平均意義下的響應轉移概率和趨勢.同時,利用胞單元中浮漂的平均速度信息形成該海域60 × 50 分辨率的平均流場圖(如圖5).為了避免小樣本容量導致Pt的統計誤差,文中借鑒插值胞映射方法的思想[35],利用胞單元中已有樣本的映射信息,生成一個映射周期內新的短時估計樣本,并將流線的切線方向作為基本誤差準則.具體插值格式如下

圖5 印度洋海域平均流場圖Fig.5 Average stream field of Indian Ocean

式中,xi(0) 為插值點在第i方向的坐標,xi(T)則表示插值得到的短時估計樣本;aij和bij則為插值待定系數,可通過真實樣本軌線來確定.基于此方法,最終形成每個胞中包含300 個有效映射的擴充樣本集.

通過對擴充樣本集的瞬態一步轉移概率矩陣Pt拓撲分析,得到印度洋海域具有吸引性的短期渦旋中心結構及其渦旋區域(如圖6).圖6 中黑色點代表短期渦旋中心,表征了該區域洋流運動形成的瞬態聚集性和平均流動趨勢,彩色區域為各短期渦旋中心所主導的渦旋中心區域.由圖6 可見,在此時間窗口中,赤道附近以北海域的洋流由西向東逐漸向遠離赤道方向流動,而南側海域則受赤道逆流的影響較大,沿逆時針遠離赤道方向環流,并且分別在緯度 ± 10°附近短暫聚集(短期渦旋中心).在南印度洋海域,由于南半球西風漂流影響,南緯50°附近的表面洋流由西向東流動,并受短期渦旋中心影響向北運動聚集在南緯40°附近海域,其中在東經50°～ 66°形成典型的環流效應(環狀渦旋).在南回歸線附近,受西澳大利亞寒流和南赤道暖流影響,形成了在印度洋中部海域的吸引結構,并使得海洋表面洋流由東向西流動.上述這些特征構成了南印度洋在該季節的逆時針環流特點.Dong 和Qin[15]指出,MH370在該時間窗口中最可能的墜海地點為南緯17°～ 33°,東經100°～ 105° 范圍內(圖6 黑色矩形框內),位于圖6 中紅色渦旋區,并受該區域的短期渦旋中心影響將在后續的一周時間內總體向西運動,這些特征與國際太平洋研究中心(International Pacific Research Center,IPRC)在該海域的相關研究結果一致[36].

圖6 2014 年3 月8 日—14 日印度洋海域短期渦旋中心(黑色點)和渦旋域(彩色區域)Fig.6 Short-term vortex cores and their regions of influence in the Indian Ocean on March 8 to 14,2014

為了進一步證實短期渦旋結構對海洋表面洋流路徑的影響,圖7 提取了典型短期渦旋中心(黑色點)附近代表性浮漂的真實軌跡(彩色實線).圖7 中能清楚地看出浮漂樣本的軌跡演化規律受到附近短期渦旋中心的主導作用,從而在短時間范圍內有靠近該結構的趨勢.但也必須指出,由歷史數據獲得的海洋全局結構僅反映了系統在時間窗口內的典型動力學特征和平均響應行為,而海洋工況的隨機性和突發性特點,使得浮漂軌跡表現為明顯的多樣性和不確定性,因而仍有一些樣本在短期內所受吸引效應并不顯著.同時,這些內在全局結構的位置和形式也是不斷發展和變化的,進而將導致漂浮物的進一步演化和聚集.

圖7 2014 年3 月8 日—14 日不同渦旋中心周圍的浮漂樣本軌跡Fig.7 Trajectories of samples near typical vortexes on March 8 to 14,2014

3 結論

本文應用數據驅動的廣義胞映射方法,研究了印度洋海域表面洋流運動的長期和短期動力特征以及演化過程,得到以下結論:

(1)提出的基于胞空間映射的數據驅動模型能夠揭示不同周期尺度(天、季節或年度)下海洋系統的內在動態結構.通過對比在不同時間跨度下的預測結果和工程實際,從概率分布和響應趨勢兩方面驗證了所提方法和結果的有效性和正確性.

(2) 從長期動態特征來看,印度洋海域在南緯20°～ 40°,東經40°～ 96°附近范圍內形成明顯的大范圍穩定的渦旋區域.該結構主導了洋流流動趨勢,并導致大部分表面漂浮物動態聚集在該區域,而在南緯40°以南和赤道附近則較為稀疏.

(3) 從短期演化來看,該區域在2014 年3 月8 日—14 日期間共存多個具有短期吸引性的渦旋結構,它們在此期間促使洋流在南緯50°附近由西向東運動,而在南回歸線附近則導致了相反的流動方向,這些特征構成了南印度洋在該季節的逆時針環流特點.