地转经验模态方法在深海声层析中的应用

马 芮，拜雅洁，张 宇，岳宜宛，崔学荣

（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]，实现对墨西哥湾区域海水的分析预测等工作，从而研究墨西哥湾区域的海洋现象。