地轉經驗模態方法在深海聲層析中的應用

馬 芮，拜雅潔，張 宇，岳宜宛，崔學榮

（1.自然資源部海上絲路海洋資源環境組網觀測技術創新中心，山東 青島 266580；2.中國石油大學（華東） 海洋與空間信息學院，山東 青島 266580）

0 引言

層析一詞指的是“切片”、“分層”的意思，奠定層析基礎的中心定理是投影-切片定理，以實現對一個物體或一種現象與過程進行分層成像。為了彌補經典物理海洋學與衛星遙感在測量海洋內部動力學現象的局限性，利用海洋動力學過程對聲傳播特性的影響，由MUNK 和 WUNSCH[1]于 1979年首次提出作為大尺度（100 km 量級）的海洋監測技術，用于海洋內波、中尺度渦、潮流、羅斯貝波和海流等現象的觀測。因此，深海聲層析利用聲波來研究深海地質、水文和生態系統，由于聲音在水中的速度主要取決于溫度、鹽度和深度，并受到流速的影響，因而通過將聲波發送到水下并記錄反射或折射的聲波信號可以測量聲信號在發射源和一個或多個接收器之間的傳播時延，得到沿路徑的平均溫度[2]。如果有多個發射接收對，可以通過計算聲學回波時間τ用來重建由發射和接收所限定的二維或三維場，進一步研究海洋內部環境，得到有關地質構造、海底地貌、沉積物類型、水下植被、生物多樣性等方面的重要數據，對海洋科學、地質學和生態學等領域的研究都具有重要意義。

地轉經驗模態（Gravest Empirical Mode，GEM）技術[3-5]是一種從垂直積分量確定海洋垂直剖面的方法。歷史水文數據用于計算溫度T、鹽度S和特定體積異常δ的特征關系，作為壓力p和垂直集成量的函數，例如聲學回波時間τ、位勢高度?或熱含量。當這些物理變量存在時，這些關系形成了稱為GEM 的查找表。GEM 查找表是一種表示T、鹽度S和特定體積異常δ與垂直積分量之間關系的工具，分別表示為TG(p,τ)、SG(p,τ)和δG(p,τ)[6]。這些表格允許研究人員在未直接測量T，S和δ的情況下，根據垂直積分量的觀測數據，通過查找表中的關系來估算這些變量的垂直剖面分布。在形成了GEM 的查找表的過程中，歷史水文數據可來自各種來源，包括海洋觀測站、航行數據和衛星測量等。單個逆時回聲儀（Current-Pressure Inverted Echo Sounder，CPIES）τ測量、衛星等時海面高度測量等可以在結合適當的GEM 查找表時提供溫度T、鹽度S和特定體積異常δ的完整垂直剖面的估計[7]。GEM 模態僅取決于參數化變量τ和p，和所需研究的變量溫度T、鹽度S和特定體積異常δ之間的相關性，并且τ由強烈依賴于T的聲速c決定，所以GEM 是很好的研究方法[8]。

1 歷史數據的統計與處理

1.1 數據來源

20 世紀末提出建立的自動海洋觀測剖面浮標陣列（Array for Real-time Geostrophic Oceanography，ARGO）為海洋研究提供數據支撐。ARGO 浮標的工作時長可以到 5年左右，剖面觀測時間間隔一般為 5 D 或10 D，浮標到達海面后會自動通過衛星將數據傳送回地面接收站，之后再次下沉到預定深度等待下一次的觀測。本文采用的ARGO 數據來自杭州全球海洋ARGO 系統野外科學觀測研究站杭州全球海洋ARGO 系統野外科學觀測研究站（http：// www.argo.org.cn/），時間跨度為1997年7月至2022年12月。

世界海洋數據庫（World Ocean Database，WOD）有全球海洋觀測和標準深度數據剖面。WOD 的開發始于 1982年，是一個用于海洋、氣候和環境研究的工具，經過 20 多年的協調努力，是將來自機構、個人研究人員的全球的各種剖面數據整合到一個單一的數據庫中。本文采用2018年更新的WOD18 數據（https：//www.ncei.noaa.gov/access/world-ocean-database-select/dbsearch.html）。

1.2 篩選與分析統計

本文使用數據均以全球數據形式提供，對于本文的重點研究區域墨西哥灣需要進行區域篩選，選擇的范圍如圖1所示，具體位置為西經-90.5°～85°，北緯24°～28°，所選擇的位置大于CPIES 布放站點。

黑框為研究區域。圖1 墨西哥灣地圖Fig.1 Map of the Gulf of Mexico

區域篩選后對剖面數據的數據量、最大最小測量深度、剖面分辨率進行統計分析。經統計，總的剖面數量為5 154 個，其中小于1 500 m 深度的剖面個數為 2 142 個，1 500～2 000 m 之間的剖面個數為574 個，2 000～3 000 m 的剖面個數為2 541個，大于3 000 m 的剖面個數僅為21 個。說明歷史數據測量深度主要集中在2 000 m 以上，占總剖面個數的 99.9%（圖2）。

圖2 剖面測量深度統計圖Fig.2 Statistical map of depths of profile measurements

經計算分析，垂直空間分辨率在0～1 m 之間的剖面數為466 個，1～2 m 之間的剖面數為3 339 個，說明歷史數據垂直空間分辨率主要集中在 0～2 m，占總剖面個數的73.3%（圖3）。

圖3 剖面垂直空間分辨率統計圖Fig.3 Statistical map of vertical spatial resolution of profiles

1.3 數據處理

海洋觀測受海流、天氣、設備等多種因素的影響，數據易丟失，質量難以保證，所以對于水文數據的處理工作十分重要，需要對歷史的水文數據進行質量控制。

首先，對剖面數據空值進行判斷篩選，空值左右相近數據若存在則進行線性插值，若不存在則舍棄該剖面；其次，GEM 構建需要選定一個參考計算層，要求深度能全覆蓋到斜壓影響深度，考慮到剖面最大測量深度，本文選取 1 000 m 為參考層，為保證參考層數據存在并且相對于1 000 m 參考層的傳播時間可計算，剔除測量深度小于1 000 m 的剖面數據；GEM 構建需要相對大的垂直空間分辨率，本文對所有符合條件的剖面數據重新以10m的垂直分辨率重采樣。由于海水幾乎不可壓縮，其密度也接近均勻，根據壓力P（dbar）：

其中可知，壓力P（dbar）與深度d具有極好的線性關系，故將海水壓力P轉化為深度索引，若原始數據中不存在所需深度d的數據，則尋找最近深度索引對應的溫鹽壓數據近似代替該深度下的溫鹽壓數據，完成均勻插值到垂向空間的粗粒度化。

最后，將不同來源剖面數據的不同參考時間格式，以1950年1月1日0 點為參考時間進行統一，最終轉化為實際測量時間。

以上為ARGO 和CTD 剖面溫度、鹽度、壓力及對應深度 數據的處理方法及過程，為計算傳播時間、對應關系擬合、構建GEM 查找表提供數據基礎。

2 墨西哥灣區域GEM 構建

為深入研究墨西哥灣區域海水的溫鹽流情況，本節將上一部分處理過的歷史溫鹽深數據整合到一起，并對每個剖面進行劃歸得到參考傳播時延，按照深度分辨率進行擬合，建立一個散點的拉格朗日矩陣，將已有的溫鹽深數據投影到此二維空間上，從而建立傳播時延與溫度、鹽度以及比容異常的經驗關系，最后針對GEM 陣進行分析[9-13]。

2.1 數據預處理和計算方法

由于歷史溫鹽壓剖面深度不一，均無法達到IES 設備深度，所以為了更好的進行GEM 查表匹配，需要統一計算得到各個剖面的參考傳播時延（取壓力面為1 000 m 時的傳播時延定義為τ1000），用該參考傳播時延進行GEM 構建。

傳播時間的計算關系式如下：

式中：P為海水壓力；g為重力加速度，取為9.8 m/s2；ρ為海水的密度；C為海水的經驗聲速。其中經驗聲速和密度都是溫度T，鹽度S和壓力P的函數。本文使用MATLAB 中的GSW 工具箱中提供的函數計算海水經驗聲速和密度。

2.2 溫度GEM 陣構建

GEM 陣的構建依賴于歷史的水溫剖面與參考傳播時間τ1000，根據 MEINEN 和 WATTS（2000）[14]提出的方法，首先需要將τ1000從小到大進行排序作為坐標點中的x向量，并將每個剖面上該深度下的所有溫度、鹽度以及比容異常作為坐標點中的y向量，共同合成坐標點，然后通過3 次樣條擬合對每個深度下的所有坐標點進行擬合，擬合時垂向的最大分辨率為10 dbar（從0～1 000 dbar），40 dbar（從1 000～2 000 dbar），這是由于剖面數量隨著深度的增加減少，因而數據的密度隨深度的減小，因此壓力網格垂直分辨率減小并允許3 次樣條中的平滑參數隨深度的減小。而對于溫度GEM 陣就是得到每一個深度上傳播時延與該深度下所有溫度的統計關系的擬合，如圖4所示（以1 000 dbar 的擬合關系為例）。

圖4 1 000 dbar 時3 次樣條擬合的溫度與傳播時間τ 的關系Fig.4 Temperature versus propagation time for 3 spline fits at 1 000 dbar

最后由各個分辨率下的所有擬合曲線建立散點的拉格朗日矩陣，從而構建溫度GEM 陣，考慮到深度2 000 m以下的變化較小，因此本文只取了2 000 m以上的溫度GEM 剖面展示，下圖5 為溫度GEM 陣。

圖5 溫度GEM 陣Fig.5 Temperature GEM array

因而從上圖基本上可以看到一個簡單的關系，對于某一個深度來說，τ1000越小，對應的溫度相對要高，因此便可以通過相應的傳播時延τ1000來對應的查出此傳播時延對應的溫度[15]。

2.3 鹽度GEM 陣構建

鹽度GEM 陣與溫度GEM 陣構建方法基本一致，但是由于海水中鹽度變化較大，所以在進行三次樣條擬合過程中出現了一些離群點，未去除離群點前的擬合曲線如圖6所示（以1 000 dbar 的擬合關系為例），顯然未去除離群點的擬合關系并不能反映傳播時延τ1000與鹽度的正確關系，所以需要去除離群點從而得到每個剖面每個深度下鹽度與傳播時延τ1000的擬合關系，如圖6所示（以1 000 dbar 的擬合關系為例）。

圖6 鹽度擬合關系Fig.6 Final salinity fit relationship

然后進行鹽度GEM 陣構建，首先由未去除離群點的各個分辨率下的所有的擬合曲線建立散點的拉格朗日矩陣，從而得到鹽度GEM 陣，如圖7所示，但顯然，該GEM 陣存在較大誤差，不能很好地進行反演，所以又構建了去除離群點之后的鹽度GEM陣。

圖7 鹽度GEM 陣Fig.7 Salinity GEM array

與溫度GEM 陣一樣，考慮到深度2 000 m 以下的變化較小，因此本文只取了2 000 m 以上的鹽度GEM 剖面展示，圖7 為最終的鹽度GEM 陣。

因而從上圖基本上可以看到一個關于傳播時延τ1000和鹽度的簡單的關系，對于某一個深度來說，τ1000越大，對應的鹽度相對要高，因此便可以通過相應的傳播時延τ1000來對應的查出此傳播時延對應的鹽度。

2.4 比容異常GEM 陣構建

比容異常GEM 陣構建相比于溫鹽GEM 陣構建多了計算比容異常的步驟，比容定義為海水密度ρ的倒數，因此某深度下的比容異常為該深度下的比容與參考比容（鹽度為35 psu，溫度0 ℃）的差值。只要求得每個剖面每個深度下的比容異常的值，擬合過程與傳播時延τ1000-鹽度擬合過程基本一致。由于比容異常是溫度T和鹽度S的函數，又因為鹽度在擬合過程中出現了離群點，所以比容異常在擬合傳播時延τ1000-比容異常的關系時也出現了離群點，未去除離群點前的擬合曲線如圖8所示（以1 000 dbar 的擬合關系為例），顯然未去除離群點的擬合關系并不能反映傳播時延τ1000與比容異常的正確關系，所以需要去除離群點從而得到每個剖面每個深度下比容異常與傳播時延τ1000的擬合關系，如下圖8所示（以1 000 dbar 的擬合關系為例）。

圖8 比容異常擬合關系Fig.8 Specific volume anomaly of spline fit versus propagation time

然后進行比容異常GEM 陣構建，首先由未去除離群點的各個分辨率下的所有的擬合曲線建立散點的拉格朗日矩陣，從而得到比容異常GEM 陣，如圖9所示，但顯然，該GEM 陣存在較大誤差，不能很好的進行反演，所以又構建了去除離群點之后的鹽度GEM 陣。

圖9 比容異常GEM 陣Fig.9 Specific volume anomaly GEM array

與溫度、鹽度 GEM 陣一樣，考慮到深度2 000 m 以下的變化較小，因此本文只取了2 000 m以上的比容異常GEM 剖面展示，圖9 為最終的比容異常GEM 陣。

3 結束語

海洋是地球氣候系統的重要組成部分，海水的溫度、鹽度和流速在調節氣候、影響海洋生物多樣性、控制全球氣候變化等方面起著關鍵作用[16]。涉及海洋循環、海洋生態系統的健康、極端天氣事件的發生與預測等方面。因此本文為了得到墨西哥灣區域海水溫度、鹽度以及流速的基本情況，在WOD18 等開源網站上下載了ARGO 和CTD 的歷史溫鹽深數據，結合這些歷史實測溫鹽深數據，建立了一個散點的拉格朗日矩陣，利用地轉經驗模態GEM 方法，將已有的溫鹽深數據投影到此二維空間上，從而構建了溫度、鹽度和比容異常的GEM陣，后續可以將校正后的傳播時延與已建立的GEM 陣結合，得到每個采樣時刻的溫度和鹽度剖面。而后根據地轉流計算公式得到地轉流場，并結合底層流速計測得的海底流速，可綜合得到絕對流速場[17-18]，實現對墨西哥灣區域海水的分析預測等工作，從而研究墨西哥灣區域的海洋現象。