利用海洋温度剖面与海表盐度反演盐度剖面方法研究

何子康 , 王喜冬 ,3, 陈志强, 范开桂

1. 河海大学自然资源部海洋灾害预报技术重点实验室, 江苏 南京210098;

2. 河海大学海洋学院, 江苏 南京 210098;

3. 南方海洋科学与工程广东省实验室(珠海), 广东 珠海 519000

盐度是决定海洋密度和比容海面高度的关键变量, 影响着洋流的强度、混合层的深度以及全球质量、热量、营养物质的输送。盐度的确定是全球海洋环流模型初始化的一个重要问题。由于对海表淡水通量的认识不足, 以及对混合过程和温盐环流强度的认识有限, 如果不进行盐度数据的同化, 同化模型中可能会出现强漂移(Haines et al, 2006)。早期观测资料中, 同时包含温盐观测的主要是温盐深仪(CTD)数据以及站点数据。虽然CTD 剖面数据具有很高的精度, 但由于观测代价较大, 得到的剖面数量有限。随着海洋观测技术的不断发展, 自2004 年以来, Argo 观测浮标每年可以提供近10 万个温度-盐度剖面数据。这些观测资料极大地提高了我们对全球海洋2000m 以上水文性质的认识, 这些盐度剖面的同化对改善海洋状态估计产品和海洋预报系统的技能发挥了重要作用。但是仍然存在大量的观测数据由于只观测了温度剖面而不能够直接同化到数值模式中。例如温深仪(BT)只用于测量温度剖面。因此有必要推断出其相应的盐度剖面, 以便有效利用这些温度剖面估计盐度对海洋密度和层结的贡献。

前人已提出一些使用温度剖面反演盐度的方法, 其主要方法包括: 局部多元回归(Hansen et al,1999; Korotenko, 2007; Thacker, 2007;Ballabrera-Poy et al, 2009; Goes et al, 2018),经验正交函数(EOF) (Carnes et al, 1994; Maes, 1999; Maes et al, 2000; Chang et al, 2011; Yang et al, 2015), 和高阶多项式(Fox et al, 2002; Marrero-Díaz et al,2006)。基于统计回归, Vossepoel 等(1999)在热带太平洋的研究表明, 结合温盐关系与海表盐度(sea surface salinity, SSS)可以显著提高盐度的估计, 结合温盐关系、SSS 和动力高度可以进一步降低盐度的估计误差。Hansen 等(1999)在热带东太平洋应用时发现, 在温盐关系中加入SSS 可以在盐度剖面反演的几乎所有深度上提供细微到实质性的改善。因此, SSS 在盐度剖面的估算中非常重要。Agarwal等(2007)利用经验正交函数(empirical orthogonal function, EOF)得到盐度剖面的主要模态, 并使用遗传算法从 TRMM(Tropical Rainfall Measuring Mission)观测的海表温度(sea surface temperature,SST)数据和Argo 观测的海表盐度来估计每个EOF主要模态的系数。Thacker (2007)通过构建盐度剖面与温度、经纬度之间的回归模型, 提出了一个相对简单但有效的方案。Ballabrera-Poy 等(2009)通过基于神经网络的非线性回归模型估计了Argo 观测的盐度剖面。但是, 盐度剖面和SSS 之间的关系不能通过神经网络或遗传算法直接得到, 需要在模型中加入额外的信息, 如SST 或观测温度剖面。

孟加拉湾是海上丝绸之路的重要航道, 位于印度洋东北部, 受印度季风的影响显著。充沛的夏季降水以及湾北部河流输入的大量淡水导致孟加拉湾内尤其是北部表层盐度有较大的变化(宣莉莉 等,2015; 张玉红 等, 2012)。西南季风期间, 外部大洋水可通过湾口西侧进入湾内, 湾内水通过湾口东侧进入赤道印度洋; 在东北季风期间, 水团交换路径发生翻转(Schott et al, 2000, 2001)。受季风影响, 阿拉伯海水也会流入孟加拉湾赤道区域, 在5—6 月及10 月至次年1 月沿赤道向东延伸, 最远可至90°E 以东的海域。此外, 在大多数月份(6 月至次年4 月)阿拉伯海水还会出现在湾内(宣莉莉 等, 2015)。孟加拉湾水与阿拉伯海水两种水团分布的季节性变化主要是由上层环流决定。11 月至次年3 月, 孟加拉湾水通过西南向的东北季风漂流及沿孟加拉湾西边界斯里兰卡沿岸南下的东印度沿岸流从湾内输运进入阿拉伯海东南部海域, 6—10 月孟加拉湾水通过湾口东侧至赤道海域之间的连续南向流输送至苏门答腊岛外海乃至赤道海域。阿拉伯海水入侵湾内有两个路径: 11 月至次年4 月是由赤道至湾口之间的连续北向流经湾口东侧进入湾内, 而在6—9 月, 阿拉伯海水由西南季风漂流输运至湾内(Shetye et al,1996; Schott et al, 2001)。表层水团的水文特征季节变化强烈, 这使得孟加拉湾上层海洋的盐度很难准确估计。

在海洋上层, 通常将近表层温度和密度都比较均匀的水层定义为混合层。在混合层以下, 存在温度随深度迅速递减、密度迅速递增的温度跃层和密度跃层。大多数情况下, 温度跃层和密度跃层深度基本一致, 但也会出现两者不重合的情况, 这时就会出现温跃层顶界深度明显大于密度跃层顶界深度的现象, 即在密度较均匀层之下、温跃层之上存在一个密度随温度深度急剧增大而温度近于均匀的水层, 称为障碍层。在孟加拉湾, 降雨和河流径流将大量淡水注入至海洋表层(Vinayachandran et al, 2002),导致密度混合层变浅, 从而在温跃层顶部形成障碍层(Lukas et al, 1991)。盐度是决定影响障碍层的主要因素, 但现有观测不足以大范围支持这方面的研究,且大量观测数据存在着盐度剖面缺测的现象, 因此通过温度剖面反演合理的盐度剖面对于了解障碍层的变化十分重要。近年来, 随着海表盐度遥感与反演技术的发展, 卫星观测的海表盐度数据增多, 为盐度剖面的反演提供了数据基础。

本文将采用回归分析方法, 基于历史观测的温盐剖面资料(图1), 建立盐度剖面与温度剖面、表层盐度、经纬度等参数之间的回归模型, 评估不同反演模型得到最优的反演方案, 同时考虑盐度卫星数据的使用对反演结果的影响, 以便将盐度卫星更好地投入实际应用。

图1 历史温盐观测各年剖面数Fig. 1 The number of historical thermohaline observation stations in each year

1 数据与方法

1.1 数据

本文所用的历史温盐剖面观测资料来自2001—2018 年的Argo 数据(ftp://data.argo.org.cn/pub/ARGO/)。这些资料要经过严格的质量控制以剔除异常数据。月均气候态数据来源于World Ocean Atlas 2013(WOA13)用于对结果进行检验评估。本文使用的卫星海表盐度资料来自 SMAP(Soil-Moisture Active Passion)提供的0.25°×0.25°的网格化海表盐度数据。船测CTD 数据来源于CLIVAR and Carbon Hydrographic Data Office (CCHDO)(https://cchdo.ucsd.edu/search?q=GO-SHIP)并用于对结果进行检验。

本研究使用2001—2016 年的Argo 数据来完成回归模型参数库的建立, 并使用2017 年的Argo 数据作为独立数据对结果进行检验(图2)。

图2 孟加拉湾区域用于建模的所用到的Argo 数据(a)以及对结果进行检验的Argo 数据(b)红框为对反演结果进行水团分析所选区域。本图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2020)4392 号的标准地图制作Fig. 2 Argo data used for modeling (a) and Argo data used for verifying the results in the Bay of Bengal (b), the red box is the selected area for water mass analysis of the inversion results

1.2 数据预处理

由于现场观测资料存在质量问题, 在使用观测剖面计算相关参数之前, 需要对观测数据进行预处理。本研究使用的质量控制方法参考自WOD13 手册(Johnson et al, 2013)。具体步骤包括:(1)剖面位置/日期/时间检查; (2)重复剖面检查; (3)深度颠倒和深度重复检验; (4)观测数据的界限检查;(5)温盐梯度检查和密度稳定性测试。质量控制结束后, 将观测层的温度和盐度剖面数据插值到标准层数据中, 并检查了标准层数据的密度稳定性, 以消除由于插值引起的密度颠倒。其中, 温盐剖面的标准层是0、2、4、 6、 8、10、15、20、 25, 30、35、40、45、50、55、60、65、70、75、80、85、90、95、100、110、120、130、140、150、160、170、180、190、200、220、240、260、280、300、350、400、500、600、700、800、900、1000m。

本研究聚焦孟加拉湾区域, 考虑到赤道区域与湾内水文性质存在较大的差异, 在赤道与湾内各选择了一个 5°×5°的方区用于反演模型的检验评估。

1.3 用温度剖面估算盐度剖面的4 种方法

模型1: 第一种估计盐度的方法是使用气候态盐度加上温盐关系, 以此来完成对盐度剖面的反演(Fox et al, 2002), 计算公式为

模型2: 第二种反演盐度的方法是在第一种方法的基础上, 加入温度的二次项, 以此完成对盐度剖面的反演(Thacker, 2007), 计算公式为

模型2 中其余各项与模型1 一致, 仅增加了aT2(z)代表温度二次项的回归系数。

模型3: 在模型2 的基础上考虑的空间参数的影响, 在拟合式中添加经纬度信息, 具体公式如下(Thacker, 2007)。

模型3 除了考虑了距离的影响外, 其余各项与模型2 一致。a x(z)代表经向的回归系数,xo代表观测经度,xc代表分析点的经度。a y(z)代表纬向的回归系数,yo代表观测纬度,yc代表分析点的纬度。

模型4: 在模型3 的基础上考虑了海表盐度的约束, 完成对盐度剖面的反演(Hansen et al, 1999;Thacker, 2007)。具体公式为

模型4 中a s(z)代表的是不同深度下表层盐度的回归系数,S(0)代表的是观测输入的表层盐度,＜S(0)＞代表的是表层盐度的气候态。

利用加权最小二乘与距离对模型局部多项式拟合, 权重系数b计算公式为

其中x、y定义与模型中一致,t为一年内的时间差,Lx、Ly、Lt为距离尺度和时间尺度。所用的时间尺度Lt是 1000h, 距离尺度取决于纬度λ, 以 km为单位(Fox et al, 2002)。

对于上述4 种回归模型, 回归系数的求解方法为

式中,A是由回归系数组成的列向量,X是由各反演模型中需要输入的观测数据组成的矩阵,B是由各个观测数据根据式(5)计算的权重系数组成的对角矩阵,Y是由观测数据组成的列向量。

1.4 比容海面高度的计算

海水物理性质的变化是影响海平面变化的重要因子, 特别是海水密度变化尤为重要。在海水质量一定的前提下, 海水密度的变化主要靠海水体积的膨胀和收缩来体现, 即以海水比容的变化来体现。一般认为海水比容变化可由两部分组成: 海水热比容变化与盐比容变化。关于海水热比容与盐比容的计算, 传统方法是藉由对于热膨胀系数和盐收缩系数的计算来完成的, 即采用Thomson 等(1987)计算比容海面高度的方法:

其中HT指热比容海平面变化,HS指盐比容变化,HSSL指总比容变化,T代表温度,S代表盐度,α代表比容,Z代表深度,z1 和z2 分别是对深度积分时所取的上限和下限。

2 结果

本部分主要评估不同反演模型反演的盐度剖面结果, 并评估不同的反演方案对结果的影响, 以得到最佳的反演方案。

2.1 误差统计

本部分共使用5 种反演方案, 方案1—4 分别对应模型1—4。值得注意的是, 方案4 中表层盐度项的输入使用的是Argo 插值得到的表层盐度。另外,实际观测常受到观测仪器的限制, 如BT 数据只能够观测到温度剖面, 考虑到盐度卫星可以提供高分辨率的表层盐度数据, 因此方案5 以卫星观测的表层盐度作为输入项来完成盐度剖面的反演。

我们计算了不同反演方案和WOA13 月均资料(将WOA13 月平均数据使用距离加权插值到观测剖面的位置)与实测盐度剖面资料之间的均方根误差(图3a)。为了检验反演模型的性能, 借鉴在天气和海洋数值预报领域采用的方法, 定义模型的技能得分为

图3 4 种不同反演盐度剖面方案中WOA 气候态与实测资料之间的均方根误差随水深的变化(a)以及4 种不同反演盐度剖面方案的技能得分(b)Fig. 3 RMSEs between WOA13 salinity and each of the four salinity reconstructions, showing in situ profiles (a) and skill scores of the four schemes (b)

式中SS 为技能得分, RMSE(reconstruct)为重构均方差误差, RMSE(WOA)为WOA 气候态均方根误差。当重构结果的精确度大于(或小于)气候态盐度场的精确度时, 技能得分为正值(或负值), 技能得分越接近1 说明重构的盐度场准确性越高。各反演结果的垂向平均均方根误差(MRMSE)结果如表1 所示。

表1 不同方案垂向平均均方根误差Tab. 1 Mean RMSEs of different schemes

结果显示, 在孟加拉湾表层由于受淡水和径流的影响, 盐度变化较大。当反演模型不引入表层盐度约束时, 反演结果存在较大误差。对比方案1 与方案2, MRMSE 从0.2011‰下降到0.1993‰, 表明模型中加入温度的二次项可以略微减小反演结果的误差。方案3 是在方案2 的基础上添加了经纬度的参数, 通过对比两者间的 MRMSE(方案 3 为0.1946‰, 方案2 为0.1993‰), 发现引入经纬度参数可以进一步减小反演模型的误差。对比方案4 与方案5, 当引入SSS 时, 上层的盐度结果得到了极大的改善。从整个海区的均方根误差统计结果看, SSS在反演盐度场中的作用基本限于100m 以上, 很好地改善了混合层的盐度结构, 但在跃层以下SSS 对于反演结果的影响很小。虽然卫星提供的SSS 较未引入SSS 的反演结果有较大的改善(方案5), 但是与使用现场观测的SSS 反演结果相比仍有较大差距(方案4), 这主要是由于卫星观测产生的观测误差造成的。

随着深度的增加, 均方根误差开始减小, 当深度大于200m 时, 反演误差随着深度不发生显著变化, 且反演误差与气候态接近, 约为0.04‰。但是反演结果的技能得分均要大于0, 说明反演结果要优于WOA 气候态。同时, 使用表层盐度约束反演的技能得分(方案4、5)要略大于不使用SSS 约束的反演结果。

2.2 反演盐度对比容海面高度的反演误差

本部分选择了一个在建模中未被使用的CTD数据用于检验反演结果在中小尺度上的表现情况。对比反演方案计算得到的比容海面高度与观测的比容海面高度, 发现反演结果的比容海面高度与观测结果基本吻合。对4 种反演结果进行评估, 使用WOA 气候态的盐度结果会有约3cm 的均方根误差,而不引入表层盐度的反演结果可以将误差减小到2.3cm; 当引入表层盐度作为约束时, 使用反演盐度结果计算的比容海面高度误差会进一步减小, 约为1.3cm(表2)。

表2 不同方案垂向平均均方根误差Tab. 2 RMSEs of different schemes

对比方案3 与方案5 反演的盐度结果, 我们以位势密度比4m 深度高0.15kg·m-3的深度来定义混合层深度(图4 中黑线表示混合层深度), 从混合层的深度变化来看, 方案3(图5c)只能微弱地反映盐度引起的混合层深度的变化, 不能体现表层盐度变化对于混合层深度的影响, 方案5 计算得到的混合层变化与观测更加接近。在混合层以下, 反演盐度断面的变化特征也与观测基本吻合, 两种反演结果得到的盐度场并没有明显的差异。但是在表层, 方案3 不能很好地反映盐度场的变化, SMAP 表层盐度数据的使用可以得到更好的表层盐度估计。而对比盐度相对气候态异常的结果, 发现在上层海洋, 方案5 得到的反演结果会有一个更大的盐度异常, 这一结果与观测较为接近。特别是在混合层, 表层盐度的引入可以较好地模拟4 月10 日至13 日出现的负异常, 而使用方案3 得到的反演结果并不能较好地得到上层盐度的变化特征。

图4 不同方案反演盐度剖面计算的比容海面高度(a)和CTD 航次轨迹(b)b 中红点代表船的起始位置。该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2020)4392 号的标准地图制作。Fig. 4 Steric height of salinity profile reconstruction calculated using different schemes (a) and the cruise track of CTD (b),and the red star represents the start of the cruise

图5 观测盐度断面及其重构结果a. 实测盐度断面; b.卫星观测作为表层盐度引入的方案5 反演盐度断面; c.方案3 反演盐度断面; d. 实测盐度异常断面; e. 方案5 反演盐度断面; f. 方案3 反演盐度异常断面。图中黑线表示混合层深度Fig. 5 Salinity profile results. (a) Observed salinity profiles; (b) scheme 5 inversion salinity profiles using satellite SSS; (c)scheme 3 inversion of salinity profiles; (d) measured salinity abnormal profiles; (e) inversion of salinity anomaly profiles using scheme 4; and (f) inversion of salinity anomaly profiles using scheme 3. The black curve indicates mixed-layer depth

2.3 反演盐度对障碍层及季节内振荡信号的恢复

本部分结果检验选择了2016—2017 年期间序列号为6901731 的Argo 浮标用于对结果进行检验(图6)。由于数据建模时只是用了2016 年之前的数据, 该浮标在建立反演模型时并未使用, 因此该检验也具有严格的独立性。

图6 观测Argo 的运动轨迹红点为浮标起始位置。该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2020)4392 号的标准地图制作Fig. 6 Trajectory of Argo profile

为探究孟加拉湾出现的季节内振荡现象, 选择盐度为34.2‰的等值线用来代表盐跃层的变化。从结果上看, 反演盐度能够较好地反映出盐度在跃层出现的季节内振荡, 且两种反演方案得到的盐跃层变化基本一致, 说明SSS 的引入对盐跃层上的盐度变化没有明显的影响。对比盐度反演结果在混合层以下的盐度异常(图 7d—f), 也发现两种方案的结果基本一致。但是在混合层, SSS 的引入较好地改善了上层盐度的估计结果。对比混合层以上的观测结果与反演结果, 发现SSS 的引入对结果有较大改善。

图7 2016—2017 年间观测盐度断面及其重构结果a. 实测盐度断面; b. 卫星观测作为表层盐度引入的方案4 反演盐度断面; c. 方案3 反演盐度断面; d. 实测盐度异常断面; e. 卫星观测作为表层盐度引入的方案4 反演盐度断面; f. 方案3 反演盐度异常断面;图中黑线表示混合层深度; a—c 图中的蓝色实心曲线代表盐度34.2‰的等值线Fig. 7 Results of observed salinity profiles and reconstruction during 2016-2017. (a) Observed salinity profiles; (b)reconstructed salinity profiles by scheme 4, using satellite SSS; (c) reconstructed salinity profiles by scheme 3; (d) observed salinity anomaly profiles; (e) reconstructed salinity anomaly profiles using scheme 4 and satellite SSS; (f) WOA13 climatology salinity anomaly profiles. The black curve is the mixed-layer depth, and the blue curve in (a-c) represents the depth of salinity isoline 34.2‰

为检验反演盐度剖面在障碍层厚度变化方面的反演能力, 分别使用Argo 观测剖面和反演盐度剖面计算其轨迹上障碍层厚度随时间的变化, 结果如图8 所示。通过对比反演结果与观测计算得到障碍层厚度的时间序列, 可以发现反演剖面计算的障碍层厚度与观测的结果变化趋势较为一致。特别是在冬季, 埃克曼抽吸和平流引起的淡水再分配共同导致的障碍层变厚, 该现象在反演结果被较好地体现出来。对比方案3 与方案5 的反演结果, 发现SSS 的引入能够进一步改善对障碍层的模拟性能。

图8 由Argo 温度剖面反演的盐度计算的障碍层厚度随时间的变化Fig. 8 Variation of barrier layer thickness with time, calculated by salinity inversion of Argo temperature profiles

2.4 水团分析

在孟加拉湾选择图2b 中所示的两个方区来研究重构结果的水团特征。在孟加拉湾内, 大量的河流径流产生的低盐度水团称为孟加拉湾水(Bay of Bengal Water, BBW)。其主要特征是温度变化不大,但上层的盐度变化比较明显。与WOA13 相比, 更多的低盐的BBW 被较好的反演出来。由图8 可知, 在孟加拉湾内的表层盐度变异性较大, 观测盐度变化从28‰至35‰, 与WOA 气候态结果相比, 更多盐度低于29‰的低盐水被反演出来。在孟加拉湾赤道区域, 盐度变化范围没有孟加拉湾内这么大, 观测盐度主要集中在33‰～35.5‰, 反演盐度的范围与观测盐度基本吻合。同时发现表层水团中不仅包含了孟加拉湾水团的信息, 还包含了部分阿拉伯海水团的信息, 且观测到的阿拉伯海水团均出现在2016年10 月至2017 年2 月。BBW 的形成主要是由于类似东印度流(East Indian Current, EIC)南下的洋流将BBW 从孟加拉湾带到了赤道。阿拉伯海水团出现在孟加拉湾近赤道区域的主要原因是西南季风形成的赤道流将阿拉伯海的水带到了孟加拉湾近赤道区域。而WOA 气候态的盐度范围主要集中33.7‰至35.1‰, 并没有捕捉到阿拉伯海水团的信息。

在200～1000m 的中层海洋, 存在的水团主要为印度中央水(Indian Central Water, ICW), 盐度为34.8‰～35.2‰(图9、10 和11)。WOA13 气候态和重建结果均与观测结果较为一致。对比反演结果与WOA 气候态, 二者都较好地反映了这一水团的特征, 没有明显的区别。

图9 0m 处孟加拉湾WOA 气候态结果和反演结果与观测结果的对比a. WOA 气候态结果与观测结果对比; b.反演结果与观测结果对比; c. 孟加拉湾赤道WOA 气候态结果与观测结果对比; d. 孟加拉湾赤道反演结果与观测结果对比。图中等值线代表位势密度的超量Fig. 9 Comparison between WOA climatology and reconstruction results with observations in the Bay of Bengal at 0m. (a)Comparison of WOA climatology and observations; (b) comparison of reconstruction and observations; (c) comparison of observations and WOA climatology on the equator in the Bay of Bengal; and (d) comparison of observations and reconstruction on the equator in the Bay of Bengal.

图10 200m 处孟加拉湾WOA 气候态结果和反演结果与观测结果的对比Fig. 10 Comparison between WOA climatology and reconstruction results with observations in the Bay of Bengal at 200 m

图11 1000m 处孟加拉湾WOA 气候态结果和反演结果与观测结果的对比Fig. 11 Comparison between WOA climatology and reconstruction results with observations in the Bay of Bengal at 1000 m

3 总结

本文采用加权最小二乘多项式拟合方法, 利用历史上Argo 观测的温盐剖面资料建立了垂向盐度与温度、经纬度和表层盐度这些变量之间的映射关系, 实现了垂向盐度剖面的反演。结果显示4 种反演方案均能得到较为合理的盐度剖面。对比4 种方案, 当不引入SSS 时, 无论是引入温度的二次项还是加入空间参数来估计盐度, 反演结果误差基本一致。而SSS 的引入, 可以优化上层盐度的反演结果。

将反演结果与观测进行对比, 结果显示使用方案1—3 反演得到的盐度剖面计算的比容海面高度会产生超过2cm 的误差, 而引入了SSS 可以将比容高度的计算误差降低到1.5cm 以下。从盐度的变化看, 表层盐度信息的引入能够较为真实地反映海洋盐度场的垂直结构和内部变化特征, 它既能够捕捉到对上混合层有重要影响的SSS 信号, 又能够反映盐度在跃层上的季节内变化。另外, 使用反演盐度剖面计算障碍层厚度并与实际观测结果进行对比,发现反演结果对障碍层具有较好的模拟性能。通过水团分析发现, 与气候态相比, 反演结果可以更好地表征上层海洋水团的变化特征。本文所构建的盐度场既可以作为观测数据的补充, 也可以与温度一起同化至模式中, 进而改善三维温、盐、流的数值再分析和预报。