中国海洋大学学报自然科学版  2026, Vol. 56 Issue (3): 23-40  DOI: 10.16441/j.cnki.hdxb.20240185

引用本文  

巩宗晴, 李芙蓉. 全球次表层海洋热浪统计特征的反演方法研究[J]. 中国海洋大学学报(自然科学版), 2026, 56(3): 23-40.
Gong Zongqing, Li Furong. Research on Retrieving Statistical Characteristics of Global Subsurface Marine Heatwaves[J]. Periodical of Ocean University of China, 2026, 56(3): 23-40.

基金项目

国家自然科学基金项目(41906011)资助
Supported by the National Natural Science Foundation of China(41906011)

作者简介

巩宗晴(1999—),女,硕士生。E-mail:1337595538@qq.com

文章历史

收稿日期:2024-04-26
修订日期:2024-10-20
全球次表层海洋热浪统计特征的反演方法研究
巩宗晴 , 李芙蓉     
中国海洋大学数学科学学院, 山东 青岛 266100
摘要:本文旨在反演长时间尺度上的全球次表层海洋热浪(Marine heatwaves, MHW),以在更长时间尺度上把握其变化趋势。本文通过初步数据分析发现次表层各层海洋热浪统计特征和月平均海水温度存在明显的空间相关性,以此统计原理为指导,采用多观测样本的地理加权回归(Multiple-replicate geographically weighted regression, MRGWR)模型作为反演模型,并基于海洋环流在2°范围内通常可以近似认为不变的特征做出了带宽的选择。本文构建了单变量指标和多变量指标两种指标体系,并将MRGWR模型与先前研究使用的海洋热浪反演方法广义线性模型(Generalized linear model, GLM)的反演效果进行了对比,证实了利用多变量指标拟合MRGWR模型在反演次表层海洋热浪方面的显著优势。利用所选最优反演模型和指标对1940—2021年间全球次表层海洋热浪的统计特征进行反演,得到其长期变化趋势,该趋势与海水温度上升趋势存在明显的对应关系。本研究为人类活动引起的全球变暖对次表层海洋热浪的影响提供了有力证据。
关键词次表层海洋热浪    物理知识驱动的统计学习    多观测样本的地理加权回归模型    广义线性模型    

海洋热浪(Marine heatwaves, MHW)是指海洋中发生的极端暖事件,会对海洋生物产生诸多负面影响,并给渔业和生态环境带来巨大损失。在过去的一个世纪里,大量的温室气体排放使得全球气候系统发生了重大变化。获取海洋热浪统计特征在全球变暖背景下的长期变化趋势具有重要的科学意义和现实价值。

近年来,已经发生了多起严重的海洋热浪事件,如2010—2011年在西澳大利亚发生的Ningaloo Niño海洋热浪事件[1]、2013—2016年东北太平洋发生的Blob海洋热浪事件[2]、2015—2016年在塔斯曼海域[3]发生的海洋热浪事件,这些事件造成了巨大的生态和经济影响,如珊瑚白化[4]、海洋无脊椎动物大量死亡[5]、渔业减产[3]。鉴于海洋热浪事件严重的破坏性后果,准确识别并量化其变化趋势越来越成为当前亟待解决的问题。

目前对于海洋热浪事件的研究主要集中在海洋表层,海洋次表层是许多海洋生物的主要栖息地,因此对次表层海洋热浪的研究显得尤为迫切。Hu等[6]发现热带西太平洋存在高强度次表层海洋热浪,平均强度峰值在150 m深度附近,最高可达8.9 ℃,是表层海洋热浪的3~6倍,且独立于表层海洋热浪,可能影响金枪鱼渔业资源。Schaeffer等[7]对澳大利亚东南大陆架MHW的研究表明,次表层水文特征对于理解MHW至关重要。Amaya等[8]发现北美大陆架海底热浪可能比表面海洋热浪更强烈且持久。Zhang等[9]指出全球海洋MHW的垂直结构受海洋动力过程影响,在过去的二十年中浅层等4种类型的MHW面积和深度都显著增长。Sun等[10]基于海洋再分析数据集确定了0~200 m水深范围的海洋热浪事件,指出仅根据海面温度识别海洋热浪的局限性,并强调地下观测监测的必要性。

海洋热浪的定义基于高频(日平均)的温度数据。但在海洋表面以下日平均海水温度数据在1993年之后才覆盖全球,而更长时期(1940年至今)的全球海水温度观测只有月平均的数据,无法直接识别海洋热浪。如何基于月平均的温度数据,在物理知识的帮助下,建立有效的反演方法来反演更长时期海洋热浪的统计特征,在更长时间尺度上把握海洋热浪的变化趋势,这是本论文要解决的问题。

Oliver等[11]曾利用逐月海表温度数据构造衍生变量训练广义线性模型(Generalized linear model, GLM),尝试重构了1900—2016年的海洋表层全球范围内海洋热浪的发生频率、平均持续时间和强度。其构建的重构自变量都是单一指标,研究采用的广义线性模型只考虑单个观测点的时间序列,忽视了观测点之间的空间相关性,未能利用邻近空间点上的样本信息。

本文将海洋表面以下50~200 m深度取4层,即50、100、150和200 m,这样的数据划分有助于更好地了解海洋热浪在不同深度层次上的分布规律。本文首先基于对次表层海洋热浪统计特征和月平均海水温度进行的初步数据分析,选择多观测样本的地理加权回归(Multiple-replicate geographically weighted regression, MRGWR)模型作为次表层海洋热浪统计特征的反演模型,并进一步依据物理知识做出了带宽的选择。进一步地,本文构建了单变量指标和多变量指标两类预测指标体系作为回归模型的预测变量。此外,本文进一步将MRGWR模型与GLM模型的反演效果进行了对比,结果表明,MRGWR模型在反演次表层海洋热浪方面具有显著优势。最后,利用所选的最优反演模型和指标,对1940—2021年期间全球次表层海洋热浪的年总频次、年总持续时间和年总累积强度三个统计特征进行了反演,得到了其长期变化趋势。

1 数据来源与初步数据分析 1.1 海洋热浪的定义与识别

本文采用的海洋热浪的定义是Hobday等[12]在2016年基于季节性变化的阈值给出的对海洋热浪的量化定义,将气候基准期内第90百分位数对应的海温值定义为海洋热浪阈值,而将海水温度连续超过海洋热浪阈值不少于5 d的海水异常增暖事件定义为海洋热浪事件。使用更极端的第99百分位数的阈值可以识别出更强烈的海洋热浪,但相比于第90百分位数的阈值捕捉到的观测样本量会减少。本文选择气候基准期内第90百分位数的海温作为海洋热浪阈值。

次表层高频(日平均)三维海水温度数据的时间范围为1993—2020年,本文选择该完整时段(1993—2020年)作为气候基准期, 计算各层气候平均值和季节性变化的阈值。

1.2 海洋热浪统计特征

本文主要对海洋热浪的年总频次、年总持续时间和年总累积强度三个统计特征进行反演,并基于这三个统计特征分析海洋热浪的长期变化趋势。各统计特征的具体定义见表 1

表 1 海洋热浪统计特征 Table 1 Statistical characteristics of marine heatwaves
1.3 数据及处理 1.3.1 数据介绍

(1) GLORYS12V1再分析数据集

全球海洋再分析数据集提供海面以下0~5 900 m的高频(日平均)三维海水温度、盐度等数据。每一层都覆盖80°S—90°N, 180°W—180°E范围内的网格,水平网格分辨率为0.083°×0.083°,共有50个深度层次,数据维度为4 320格(经向)×2 041格(纬向)×50层(纵向),时间跨度为1993年1月至2020年12月。采用了先进数据同化方法有效融合观测数据和模型模拟,提高数据集的准确性[13]。本研究利用该数据集识别1993—2020年间的次表层海洋热浪统计特征。

(2) IAPv4观测数据集

大气物理研究所提供的IAPv4数据集是全球范围低频(月平均)三维海水温度观测数据集,覆盖海面以下0~6 000 m,水平分辨率为1°×1°,共有119个深度层次,数据维度为360格(经向)×180格(纬向)×119层(纵向),时间跨度为1940年1月至2021年12月,本研究基于此数据集计算1940—2021年海洋热浪统计特征的反演值。

1.3.2 数据预处理

(1) 剔除缺失值

剔除陆地、海冰点以及任何月份存在温度数据缺失的网格点,不纳入模型计算。

(2) 网格插值

GLORYS12V1数据集在次表层各层的水平分辨率为0.083°×0.083°,数据维度4 320格(经向)×2 041格(纬向)×50层(纵向),数据量过大,对此数据集进行三线性插值,插值到水平分辨率为0.5°×0.5°的网格上,且深度取50、100、150和200 m,得到维度为720格(经向)×360格(纬向)×4层(纵向)的数据集。而用以反演的月度观测数据水平分辨率为1°×1°,水平面的数据维度360格(经向)×180格(纬向),在IAPv4数据集不同深度层次进行双线性插值,插值到水平分辨率为0.5°×0.5°的网格上,统一成720格(经向)×360格(纬向)×4层(纵向)的网格尺寸,实现空间位置的匹配。将插值后的月平均海水温度代入反演模型,即可得到海洋热浪统计特征反演值。

(3) 间隔取样

由于需要对比多个反演指标和模型的反演效果,在所有空间点上进行计算的时间成本是巨大的。因此,在计算误差时,在50~200 m深度范围内的四个层次上,均取网格尺寸为144格(经向)×72格(纬向)的样本,即在各层每间隔5个网格点进行一次取样。但是,在拟合空间统计模型时,保持周边观测数据的密度不变,确保采样点周边依然覆盖着高分辨率的网格数据。这种处理方式不仅大幅减少了计算时间,还确保了采样点周边的原始数据密度在模型拟合时得以保持,从而能够准确反映全体数据的拟合结果。

1.4 初步数据分析

构建反演模型之前,本文根据海洋表面以下50~200 m四个深度层次1993—2020年期间的再分析数据集,以海洋热浪的完整计算周期作为基准期(1993—2020年),并由此定义海洋次表层各层的月气候平均值和第90百分位数阈值,在此基础上,利用各层全球网格点的高频(日平均)海水温度数据,逐年识别海洋热浪事件。

通过分析海洋次表层的海洋热浪与月度海水温度的全球分布,为模型的选择提供依据。

图 13给出了海洋热浪三个统计特征的全球分布热力图,图中展示的是次表层各层1993—2020年数据的平均值。可以看出,各图中相邻区域的颜色呈渐进变化,未见明显的空间突变,这表明次表层海洋热浪的统计特征在空间上具有连续性,引入周围观测的空间影响可以为海洋热浪的反演提供更多的信息。

图 1 次表层各层海洋热浪年总频次全球分布 Fig. 1 Global distributions of annual total frequency in each subsurface layer
图 2 次表层各层海洋热浪年总持续时间全球分布 Fig. 2 Global distributions of annual total duration in each subsurface layer
图 3 次表层各层海洋热浪年总累积强度全球分布 Fig. 3 Global distributions of annual total cumulative intensity in each subsurface layer

同时,本节选取次表层50 m处的1、4、7和10月的海水温度平均值,观察其全球分布情况。观察图 4可以发现,月平均海水温度整体上呈带状分布,从赤道向两极,颜色随纬度发生规律性变化,表明温度自赤道向两极递减,具有典型的空间结构,表现出显著的空间自相关性。

图 4 次表层50 m月平均海水温度的全球分布 Fig. 4 Global distributions of subsurface 50 m monthly mean seawater temperature

根据地理学第一定律,任何事物都存在空间相关性,距离越近,空间相关性越强。上述初步数据分析表明,海洋热浪统计特征与月平均海水温度数据均存在明显的空间相关性。据此推测,二者之间的关系在空间上也可能具有连续性,因此,引入空间统计模型进行拟合是合理的。

2 次表层海洋热浪统计特征的反演方法

基于日平均海水温度数据集识别次表层各层海洋热浪的统计特征作为响应变量,基于月平均海水温度数据构造反演指标作为解释变量,以2001—2020年共20 a的数据作为训练集去拟合各模型,1993—2000年共8 a的数据作为测试集对模型进行评估,最终,根据测试集均方根误差最小原则,为次表层每层的各个统计特征筛选出最优的反演指标和模型。

2.1 反演模型

通过初步数据分析,在海洋和大气动力过程对海温变化的重要调控作用下,次表层月度海水温度和海洋热浪统计特征在空间上存在相邻相近的特点,基于全球观测点在空间上的相关性,周围观测对于当前观测点海洋热浪事件的影响能够为海洋热浪的反演提供更多的信息。因此,本文尝试引入空间统计模型——多观测样本的地理加权回归(Multiple-replicate geographically weighted regression, MRGWR)模型作为次表层海洋热浪统计特征的反演模型。

多观测样本的地理加权回归模型是地理加权回归模型(Geographically weighted regression, GWR)在多观测条件下的推广[14]。GWR假定在每一个观测点上只有一个观测样本,而MRGWR允许每个观测点上有多个观测样本。设在空间位置((ui, vi))上具有pi个观测样本,观测值分别为

$ \begin{gathered} \left(y_{i 1} ; x_{i 11}, x_{i 12}, \cdots, x_{i 1 d}\right), \left(y_{i 2} ; x_{i 21}, x_{i 22}, \cdots, x_{i 2 d}\right), \cdots, \\ \left(y_{i p_i} ; x_{i p_i 1}, x_{i p_i 2}, \cdots, x_{i p_i d}\right) 。\end{gathered} $

其中,$\left(y_{i r} ; x_{i r 1}, x_{i r 2}, \cdots, x_{i r d}\right), r=1, 2, \cdots, p_i$是在位置$\left(u_i, v_i\right)$处各变量的第r组观测值,i=1, 2, …, n

MRGWR表达式为

$ \begin{gathered} y_{i r}=\beta_0\left(u_i, v_i\right)+\sum\limits_{k=1}^d \beta_k\left(u_i, v_i\right) x_{i r k}+\varepsilon_{i r}, \\ r=1, 2, \cdots, p_i, i=1, 2, \cdots, n, \end{gathered} $ (1)

在MRGWR模型中,在观测点(ui, vi)处的参数$\boldsymbol{\beta}_i=\left[\beta_1\left(u_i, v_i\right), \beta_2\left(u_i, v_i\right), \cdots \beta_d\left(u_i, v_i\right)\right]^{\mathrm{T}}$估计值为

$ \widehat{\boldsymbol{\beta}}_i=\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{W}_i \boldsymbol{X}\right)^{-1}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{W}_i \boldsymbol{Y}\right), $

其中,

$ \begin{gathered} \boldsymbol{X}=\left(\boldsymbol{X}_1, \boldsymbol{X}_2, \cdots, \boldsymbol{X}_i, \cdots, \boldsymbol{X}_n\right)^{\mathrm{T}}, \\ \boldsymbol{X}_i=\left(\begin{array}{cccc} x_{i 11} & x_{i 12} \cdots & x_{i 1 d} \\ x_{i 21} & x_{i 22} \cdots & x_{i 2 d} \\ \vdots & \vdots & \vdots & \vdots \\ x_{i p_i 1} & x_{i p_i 2} \cdots & x_{i p_i d} \end{array}\right), \\ \boldsymbol{y}=\left(y_{11}, \cdots, y_{1 p_1}, \cdots y_{n 1}, \cdots, y_{n p_1}\right)^{\mathrm{T}}, \end{gathered} $
$ \boldsymbol{W}_i=\boldsymbol{W}\left(u_i, v_i\right)=\left(\begin{array}{cccc} \boldsymbol{W}_{i 1} & 0 & \cdots & 0 \\ 0 & \boldsymbol{W}_{i 2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \boldsymbol{W}_{i n} \end{array}\right)。$

$\boldsymbol{W}_{i j}=\operatorname{diag}\left(\omega_{i j}, \cdots, \omega_{i j}\right)_{p_j \times p_j}$是在对点(ui, vi)上的观测估计参数时,所赋予观测点(uj, vj)上pj个观测样本的空间权值矩阵,其中diag为构造对角矩阵函数。

2.1.1 核函数

MRGWR模型在求回归参数 βi时,空间权重矩阵 Wi是核心,需通过空间核函数确定空间权重。若样本点距离回归点越近,其对回归点的重要性越大,赋予较大的权重,远离回归点的样本点赋予较小的权重。样本点(uj, vj)与回归点(ui, vi)之间的权值随着空间距离$d_{i j}=\sqrt{\left(u_i-u_j\right)^2+\left(v_i-v_j\right)^2}$的增大而减小。

常见的核函数包括Gauss核函数、Exponential核函数、Bi-square核函数等。本文在对海洋热浪统计特征拟合MRGWR模型时采用指数核函数,其表达式为

$ \omega_{i j}=\exp \left(-\frac{d_{i j}}{b}\right) 。$ (3)

其中b是带宽参数,决定了回归点受周围观测影响的范围大小,对于确定每个观测权重至关重要。

2.1.2 带宽

关于最优带宽的确定,通常会考虑交叉确认法、AICC信息准则等方法,其中交叉确认法最为常用。然而本文采用次表层三维海水温度数据集数据量过大,交叉验证法在计算时间方面的成本无法承受。

因此本文尝试使用物理知识来指导带宽的选择。近期的研究表明[15],海洋中尺度涡旋引起的热输运异常是海洋热浪的主要生成机制。海洋涡旋的统计特性是由海洋环流的结构决定的。海洋环流在空间上缓慢连续变化,在2°的范围内通常可以近似视为不变。由此我们推测,月平均数据(反映的海洋环流的信息)与海洋热浪的统计特征(反映的是涡旋的统计特征)之间的关系在2°的范围内可以近似视为常数。为验证我们的推测,后续本文使拟合MRGWR模型时将带宽分别取{0.01, 0.05, 0.1, 0.5, 1, 2, 3}并比较不同带宽下测试集的误差,数据显示在带宽取到2时测试集误差较小,反演效果理想。如果选择更大的带宽参数,会小程度提高模型的准确性,模型结果对于带宽的增大并不敏感,但是计算成本却显著增加。因此,本文拟合MRGWR模型时选择带宽为2。

2.2 反演指标 2.2.1 单变量指标

单变量反演指标的构造以海洋热浪事件的定义为基础,共8个,既有离散变量,也有连续变量,可以比较全面地衡量年度汇总的月平均海水温度变化对海洋热浪事件的影响。为方便表示,将这8个反演指标以X1, X2, …, X8进行表示。表 2给出了各变量的具体含义。

表 2 单变量指标及其含义 Table 2 Univariate indexes and their meanings
2.2.2 多变量指标

多变量指标相较于单变量指标能够同时兼顾多个自变量对因变量的影响,包含更多信息。本小节使用海洋表面以下50 m深度的数据进行计算,次表层其他深度层次采用相同的方法构造多变量指标。

(1) 综合指标1

本文将8个单变量指标共同作为综合指标1,对海洋热浪统计特征进行反演。

(2) 综合指标2

若将所有变量都纳入所研究的模型,无形中增加了模型的复杂度和计算成本,同时变量之间可能存在冗余信息,造成多重共线性问题。

利用方差膨胀因子(Variance inflation factor,VIF)对单变量指标进行多重共线性检验,由表 3发现VIF值存在大于10的情况,表明各单变量指标之间存在多重共线性。

表 3 各单变量指标的VIF值 Table 3 VIF values for all univariate indexes

计算各变量之间的相关系数,得到变量X1X2之间的相关系数为0.992 9,X3X4之间的相关系数为0.932 3,以及X5X6之间的相关系数为0.938 6,皆大于0.9,这表明这些变量之间具有较高的共线性(见表 4)。结合VIF值以及相关系数对各变量进行分析,对8个变量进行人为剔除,最终保留5个变量,即X2X4X6X7X8。剔除后,自变量之间不存在共线性,如表 5所示,VIF值均小于10。本文将这5个指标作为综合指标2。

表 4 各单变量指标之间的相关系数 Table 4 The correlation coefficient between all univariate indexes
表 5 变量选择后的5个变量的VIF值 Table 5 VIF values of the five variables after the variable selection

(3) 综合指标3

此外,本文运用主成分分析这一方法对8个单变量指标组成的综合指标1实施降维处理,提取出包含原始数据核心信息的综合指标(见表 6)。

表 6 8个单变量指标主成分分析结果 Table 6 Principal component analysis results of eight univariate indexes

表 6的数据可知,前两个主成分的累积贡献率达到85.458 0%,可认为前两个主成分基本提取了原数据绝大部分信息,本文以前两个主成分PCA1PCA2作为综合指标3。

(4) 综合指标4

由于年度汇总指标无法体现出单个月的海水温度变化对于海洋热浪事件的影响,导致单个月的温度变化被掩盖,本文建立每月海水温度距平指标,以每个月的月平均海水温度与该月气候平均值的差值构成的共12个月的海水温度距平作为综合指标4,能够反映出每个月海水温度的异常程度,以T1, T2…, T12来表示。

$ T_i=T(i)-T_M(i), i=1, 2, \cdots, 12 \text { 。} $ (4)

式中:T(i)表示第i月的月平均海水温度;TM(i)表示一年中第i个月的月气候平均值。

(5) 综合指标5

本文考虑利用主成分分析对构建的综合指标4中12个变量进行降维处理。由数据可知,前5个主成分的累计贡献率达到82.642 4%(见表 7),可认为前5个主成分就基本提取了原数据绝大部分信息,本文以前5个主成分PC1PC2PC3PC4PC5为作为综合指标5。

表 7 海水温度距平指标主成分分析结果 Table 7 Principal component analysis results of mean seawater temperature anomalies

将各多变量指标汇总如表 8所示。

表 8 多变量指标及其含义 Table 8 Multivariate indexes and their meanings
2.3 反演效果

本节对海洋表面以下50 m进行数据上的详细分析,对海洋次表层另外3个深度层次采用同样的分析方法。

2.3.1 MRGWR模型反演效果 2.3.1.1 MRGWR模型全球平均RMSE

以2.2节中的单变量指标和多变量指标两个指标体系分别作为反演指标,MRGWR模型作为反演模型,对发生在次表层50 m处的海洋热浪的三个统计特征进行反演。对于每一反演指标,计算拟合MRGWR反演模型时在测试集上对响应变量全球面积加权的均方根误差(Root mean squared error, RMSE),数值为四舍五入保留四位小数后的数值,结果如表 9所示。

表 9 反演指标对次表层50 m三个统计特征拟合MRGWR模型的测试误差 Table 9 Test errors of MRGWR models fitted by retrieval indexes to three statistical characteristics of subsurface 50 m

当用MRGWR模型作为次表层50 m海洋热浪统计特征的反演模型时,多变量指标包含更多信息,优于单一指标的反演效果。由测试误差选出的三个统计特征的最优反演指标都是多变量指标:年总累积强度的最优反演指标是综合指标1,而年总频次和年总持续时间的最优反演指标均为综合指标2。

2.3.1.2 MRGWR模型回归系数

在海洋表面以下50 m深度,以综合指标2为反演指标,基于20 a训练集对全球所有有效网格点对海洋热浪年总频次、年总持续时间拟合MRGWR模型,得到6个回归系数的全球分布图。在图 5图 6中,第一个子图展示的是常数项的全球分布,第二个子图到第六个子图分别展示了综合指标2中包含的5个单变量指标的系数的全球分布情况。

图 5 用综合指标2对次表层50 m海洋热浪年总频次拟合MRGWR模型的回归系数 Fig. 5 Regression coefficients of MRGWR model fitted by comprehensive index 2 to annual total frequency of subsurface 50 m
图 6 用综合指标2对次表层50 m海洋热浪年总持续时间拟合MRGWR模型的回归系数 Fig. 6 Regression coefficients of MRGWR model fitted by comprehensive index 2 to annual total duration of subsurface 50 m

以综合指标1为反演指标对海洋热浪年总累积强度拟合MRGWR模型,分别得到9个回归系数,在图 7中第一个子图展示的是常数项数值的全球分布,第二个子图到第九个子图分别展示了综合指标1中包含的8个单变量指标的系数的全球分布情况。

图 7 用综合指标1对次表层50 m海洋热浪年总累积强度拟合MRGWR模型的回归系数 Fig. 7 Regression coefficients of MRGWR model fitted by comprehensive index 1 to annual total cumulative intensity of subsurface 50 m

各个单变量指标的系数分布呈现出多样性,但相邻空间位置之间的回归关系呈现出一定的相似性,这一发现与我们先前基于物理知识的推测相符,即海洋热浪的年度统计特征和每年的月度海水温度数据之间的关系在全球范围内分布,具有相邻相近的特点。

2.3.2 MRGWR与GLM、GWR进行对比

以2.2中的单变量指标和多变量指标两个指标体系分别作为反演指标,用GLM、GWR对发生在海洋表面以下50 m海洋热浪的三个统计特征进行反演,并将结果与MRGWR模型得到的结果进行对比(见表 1012)。对于每一反演指标,计算测试集上对响应变量的反演误差(RMSE),数值为四舍五入保留四位小数后的数值。

表 10 海洋次表层50 m各反演指标对年总频次的测试误差 Table 10 Test errors of three models fitted by retrieval indexes to annual total frequency of subsurface 50 m  
表 11 海洋次表层50 m各反演指标对年总持续时间的测试误差 Table 11 Test errors of three models fitted by retrieval indexes to annual total duration of subsurface 50 m  
表 12 海洋次表层50 m各反演指标对年总累积强度的测试误差 Table 12 Test errors of three models fitted by retrieval indexes to annual total cumulative intensity of subsurface 50 m  

综合考虑各反演指标和反演模型的反演效果,综合指标2基于MRGWR模型对海洋热浪年发生频次的反演RMSE最低。而在对海洋热浪年总持续时间的反演中,综合指标2基于MRGWR模型的反演RMSE最低。在对海洋热浪年总累积强度的反演中,综合指标1基于MRGWR模型的反演RMSE最低,选择它作为年总累积强度最优反演模型。海洋次表层50 m各统计特征的最优反演指标与对应的模型如表 13所示。

表 13 海洋次表层50 m各统计特征的最优反演指标与模型 Table 13 Optimal retrieval indexes and models for three statistical characteristics of subsurface 50 m

同样地,用各单变量指标和多变量指标拟合GLM、GWR、MRGWR模型,得到海洋表面以下其他深度层次的测试误差,各层的最优反演指标和最优模型如表 14所示。可以看到在海洋表面以下各层对于三个海洋热浪统计特征的最优反演指标都是多变量指标,同时注意到各层的三个海洋热浪统计特征的最优反演模型都是MRGWR模型,展现出利用物理知识选择模型的重要性,以及MRGWR模型在处理具有多个观测样本的空间数据方面的优越性,MRGWR是最优的反演方法。

表 14 海洋次表层各层各统计特征的最优反演指标与模型 Table 14 Optimal retrieval indexes and models for three statistical characteristics of subsurface layers

在拟合GWR模型时,每个观测点应仅包含一个观测样本。若存在多个观测样本,则需对每个样本分别估计回归系数并计算其中位数,然而这种方法未充分考虑同组多个观测样本的回归系数应一致的原则。相比之下,MRGWR模型能够处理每个观测点存在多个观测样本的情况,从而提高了计算的准确性。

在拟合GLM模型的过程中,由于未能充分利用回归关系在相邻空间点上所呈现出的相似性,使得模型的估计效果不如MRGWR。同时,引入多变量指标后,MRGWR模型的RMSE值显著下降,而GLM模型由于存在过拟合现象,RMSE反而增加,本文计算了次表层50 m处采用不同反演指标和反演模型反演年总频次得到的训练误差和测试误差进行比较,来说明过拟合现象的存在。

观察表 15发现,在运用多变量指标以及部分单变量指标对年总频次进行GLM模型拟合时,存在一种现象:模型在训练集上的表现极为出色,然而应用于测试集时却出现了异常大的误差,这种情况被视为过拟合现象的表现。相较之下,MRGWR模型在训练过程中的误差与测试过程中的误差保持了高度的一致性,且误差值均维持在较低水平。这主要得益于样本容量的增加,使得这种模型能够有效地避免过拟合现象的发生。

表 15 拟合次表层50 m海洋热浪年总频次的训练误差和测试误差 Table 15 Training errors and test errors to fit annual total frequency of subsurface 50 m  
3 次表层海洋热浪统计特征长期变化趋势反演

本章基于海洋次表层50~200 m四个深度层次的最优反演指标与反演模型,由1940—2021年的月平均数据计算反演指标,对长时间尺度上的海洋各层海洋热浪统计特征(年总频次、年总持续时间和年总累积强度)进行反演,1940—2021年共82 a间各层全球各格点及全球平均的海洋热浪统计特征变化趋势。

3.1 年总频次

图 8给出了1940—2021年次表层四层海洋热浪年总频次斜率估计的全球分布。

(“///”表示斜率显著大于0,“★★★”表示斜率显著小于0。"///" represents a significant slope greater than 0, and "★★★" represents a significant slope less than 0.) 图 8 次表层各层MHW年总频次斜率估计值 Fig. 8 The estimated slopes of annual total frequency in each subsurface layer

经过数据分析得知,在显著性水平为0.05的条件下,在海洋表面以下50~200 m的四层中,60.703%~68.853%的海域其斜率估计值显著大于0,且该比例随深度增加而减小,这些海域呈现出海洋热浪日益频繁的趋势。这些变化趋势显著上升的大部分海域的增长速度不超过0.2次/10a,占比在全球海域54.068%~63.454%。全球范围内仅有2.070%~5.391%的海域海洋热浪年总频次有显著下降的趋势,未表现出明显变化趋势的海域从海洋表面以下50 m的占比29.077%增加到海洋表面以下200 m的34.244%。

将海洋次表层各层1940—2021年的海洋热浪年总频次反演值进行面积加权平均,得到了各层全球平均的海洋热浪年总频次82 a的时间序列,如图 9所示。在海洋次表层各层,全球海洋热浪年总频次都呈明显的上升趋势。本文用泰尔-森(Theil-Sen)斜率估计海洋热浪的变化趋势。

图 9 次表层各层1940—2021年海洋热浪年总频次全球平均值 Fig. 9 The global averages of annual total frequency in each subsurface layer during the period of 1940 to 2021

表 16所示,在海洋次表层各层,全球海洋热浪年总频次呈明显的上升趋势,Theil-Sen斜率估计值显示,平均每10 a增加0.068~0.072次。同时,由各层的95%置信区间可知,显著性水平0.05的条件下,认为各层的海洋热浪年总频次上升趋势都是显著的。

表 16 次表层各层海洋热浪年发生频次的Theil-Sen斜率估计值和95%置信区间 Table 16 Theil-Sen slope estimates and 95% confidence intervals of annual MHW frequency in each subsurface layer  
3.2 年总持续时间

图 10给出了海洋次表层四层海洋热浪年总持续时间变化趋势的全球分布。对数据分析研究得到,在显著性水平为0.05的条件下,斜率估计显著大于0的海域占比从海洋表面以下50 m深度的64.319%升至100 m处的64.656%,再降至200 m处的57.910%。增长速度平均每10 a约不超过5 d的海域占全球海域的54.345%~62.102%,这一数值随着深度的加深而减小,有显著下降趋势的海域占比在1.859%~5.982%。在显著性水平0.05下,未表现出明显变化趋势的海域占比为32.842%~36.280%,这些海域的斜率估计值与0无显著差异,表明其海洋热浪年总持续时间未呈现明显的趋势性变化。

(“///”表示斜率显著大于0,“★★★”表示斜率显著小于0。"///" represents a significant slope greater than 0, and "★★★" represents a significant slope less than 0.) 图 10 次表层各层MHW年总持续时间斜率估计值 Fig. 10 The estimated slopes of annual total duration in each subsurface layer

通过对海洋次表层各层1940—2021年海洋热浪年总持续时间反演值进行面积加权平均,得到了各层全球平均的海洋热浪年总持续时间82 a的时间序列(见图 11)。

图 11 次表层各层1940—2021年海洋热浪年总持续时间全球平均值 Fig. 11 The global averages of annual total duration in each subsurface layer during the period of 1940 to 2021

对于海洋次表层各层,海洋热浪年总持续时间的Theil-Sen斜率估计值和95%置信区间如表 17所示。在海洋次表层各层,全球海洋热浪年总持续时间都呈上升趋势,Theil-Sen斜率估计值显示,平均每10 a增加1.408~1.493 d,在150 m处的斜率估计值达到最快,为1.493 d/10a。同时,由于各层的95%置信区间均不包含0,在显著性水平0.05下,可以认为各层海洋热浪年总持续时间的上升趋势均具有统计显著性。

表 17 海洋次表层各层海洋热浪年总持续时间的Theil-Sen斜率估计值和95%置信区间 Table 17 Theil-Sen slope estimates and 95% confidence intervals of annual total duration in each subsurface layer  
3.3 年总累积强度

图 12给出了海洋次表层四层海洋热浪年总累积强度变化趋势的全球分布。对海洋次表层各层进行数据分析,在显著性水平为0.05的条件下斜率估计值显著大于0的海域占比在50 m处的60.131%随着深度的加深降至200 m深度处的56.419%。增长速度不超过8 ℃/10 a的海域占比在55.030%~58.332%。海洋热浪年总累积强度斜率估计值显著小于0的海域占比在4.319%~7.353%。在0.05显著性水平下未表现出明显的趋势变化的海域占比在32.516%~38.327%。

(“///”表示斜率显著大于0,“★★★”表示斜率显著小于0。"///" represents a significant slope greater than 0, and "★★★" represents a significant slope less than 0.) 图 12 次表层各层MHW年总累积强度斜率估计值 Fig. 12 The estimated slopes of annual total cumulative intensity in each subsurface layer

将海洋次表层各层1940—2021年的海洋热浪年总累积强度反演值进行面积加权,得到了各层全球平均的海洋热浪年总累积强度82 a的时间序列,如图 13所示,各层的年总累积强度均呈现总体上升趋势。

图 13 次表层各层1940—2021年海洋热浪年总累积强度全球平均值 Fig. 13 The global averages of annual total cumulative intensity in each subsurface layer during the period of 1940 to 2021

对于海洋次表层各层,海洋热浪年总累积强度的Theil-Sen斜率估计值和95%置信区间如表 18所示。

表 18 海洋次表层各层海洋热浪年总累积强度的Theil-Sen斜率估计值和95%置信区间 Table 18 Theil-Sen slope estimates and 95% confidence intervals of annual total cumulative intensity in each subsurface layer  

在海洋次表层各层,全球海洋热浪年总累积强度都呈上升趋势,Theil-Sen斜率估计值显示平均每10 a约增加1.410~2.073 ℃,在次表层100和150 m深度处的增长速度快于50和200 m深度处。同时,由于各层的95%置信区间均不包含0,在显著性水平0.05下,可认为海洋热浪年总累积强度有显著的上升趋势。

3.4 海洋热浪统计特征变化趋势与全球变暖间的关系

为进一步探索全球变暖对于海洋热浪统计特征的影响,识别人类活动对海洋热浪长期趋势变化的影响,本节将各层全球范围内海洋热浪年总频次、年总持续时间和年总累积强度的趋势分布,与各层全球范围内海水温度的长期趋势进行对比。

图 14显示,在海洋次表层50~200 m各层中,1940—2021年82 a间年平均海水温度斜率估计值显著大于0的海域占绝大部分,表明存在全球变暖趋势。该比例从50 m层的82.740%下降至100 m层的73.240%,随后逐渐增加至200 m层的75.708%,年均海水温度每10 a上升不超过0.1 ℃的海域占比在63.202%~70.265%。斜率估计值显著小于0的海域占比从次表层50 m处的0.125%增加到100 m处的0.793%,下降到150 m处的0.451%,又增加到200 m处的0.927%。

(“///”表示斜率显著大于0,“★★★”表示斜率显著小于0。"///" represents a significant slope greater than 0, and "★★★" represents a significant slope less than 0.) 图 14 次表层各层全球海水温度斜率估计值 Fig. 14 The estimated slopes of global seawater temperature in each subsurface layer

对比发现,在海洋次表层各层中,海洋热浪的三个主要统计特征与海水温度的变化趋势存在明显的对应关系,特别是年总频次的变化趋势与海水温度的变化趋势的对应性更为显著。在海水温度显著上升的区域,海洋热浪的年总频次也呈现出明显的上升趋势,在次表层50~200 m范围内,海水温度显著上升的区域中,有80.179%~86.611%的区域也观察到了海洋热浪年总频次的显著上升。反之,在海水温度显著下降的区域,海洋热浪事件的数量和年总持续时间则显著减少。

经过综合考量与分析,从长期的时间尺度观察,海洋次表层中海洋热浪事件的增加趋势与海水温度的上升趋势呈现出紧密的关联。此现象揭示了全球变暖对海洋热浪变化趋势的影响,为人类活动导致的全球变暖影响次表层海洋热浪提供了证据。

4 结语

本文对全球次表层海洋热浪的反演方法展开了研究。通过初步数据分析发现次表层各层海洋热浪的不同统计特征和月平均海水温度均在空间上存在显著变化,同时,统计特征和月平均海水温度数据在空间上均呈现相邻相近的特点,存在明显的空间相关性。本文以物理知识为指导采用MRGWR模型作为反演模型进行反演,并应用物理知识做出了带宽的选择。MRGWR模型不仅能处理每个观测点存在多个观测样本的情况,还能充分考虑回归关系在相邻空间点上所呈现出的相似性,用MRGWR模型拟合后得到的回归系数随空间的分布验证了它们之间的回归关系也存在着相邻相近的特点。本文构建了单变量指标和多变量指标两种指标体系,由于多变量指标包含更为丰富的信息,因此对每个统计特征拟合MRGWR模型时,最优反演指标均为多变量指标。引入多变量指标后,GLM模型存在过拟合现象,而MRGWR模型则因样本容量的增加而表现出良好的性能。与GLM、GWR模型对比发现MRGWR模型将反演效果进行了提升。

对于海洋次表层各层,在显著性水平为0.05的条件下,全球大部分海域的海洋热浪的年总频次、年总持续时间与年总累积强度都显著上升,在1940—2021年共82 a的长时间尺度上,海洋次表层各层中海洋热浪的三个主要统计特征与海水温度的变化趋势之间存在明显的对应关系,特别是年总频次的变化趋势与海水温度的变化趋势之间的对应性更为显著,这为人类活动导致的全球变暖对次表层海洋热浪的显著影响提供了证据。

参考文献
[1]
Pearce A F, Feng M. The rise and fall of the "marine heat wave" off Western Australia during the summer of 2010/2011[J]. Journal of Marine Systems, 2013, 111-112: 139-156. DOI:10.1016/j.jmarsys.2012.10.009 (0)
[2]
Bond N A, Cronin M F, Freeland H, et al. Causes and impacts of the 2014 warm anomaly in the NE Pacific[J]. Geophysical Research Letters, 2015, 42(9): 3414-3420. DOI:10.1002/2015GL063306 (0)
[3]
Oliver E C J, Benthuysen J A, Bindoff N L, et al. The unprecedented 2015/16 Tasman Sea marine heatwave[J]. Nature Communications, 2017, 8: 16101. DOI:10.1038/ncomms16101 (0)
[4]
Feng Y, Bethel B J, Dong C, et al. Marine heatwave events near Weizhou Island, Beibu Gulf in 2020 and their possible relations to coral bleaching[J]. Science of the Total Environment, 2022, 823: 153414. DOI:10.1016/j.scitotenv.2022.153414 (0)
[5]
Last P R, Gledhill D C, Hobday A J, et al. Long-term shifts in abundance and distribution of a temperate fish fauna: A response to climate change and fishing practices[J]. Global Ecology and Biogeography, 2011, 20: 58-72. DOI:10.1111/j.1466-8238.2010.00575.x (0)
[6]
Hu S, Li S, Zhang Y, et al. Observed strong subsurface marine heatwaves in the tropical western Pacific Ocean[J]. Environmental Research Letters, 2021, 16(10): 104024. DOI:10.1088/1748-9326/ac26f2 (0)
[7]
Schaeffer A, Roughan M. Subsurface intensification of marine heatwaves off southeastern Australia: The role of stratification and local winds[J]. Geophysical Research Letters, 2017, 44(10): 5025-5033. DOI:10.1002/2017GL073714 (0)
[8]
Amaya D J, Jacox M G, Alexander M A, et al. Bottom marine heatwaves along the continental shelves of North America[J]. Nature Communications, 2023, 14: 13614. (0)
[9]
Zhang Y, Du Y, Feng M, et al. Vertical structures of marine heatwaves[J]. Nature Communications, 2023, 14: 7063. DOI:10.1038/s41467-023-42844-9 (0)
[10]
Sun D, Li F R, Jing Z, et al. Frequent marine heatwaves hidden below the surface of the global ocean[J]. Nature geoscience, 2023, 16(12): 1099-1104. DOI:10.1038/s41561-023-01325-w (0)
[11]
Oliver E C J, Donat M G, Burrows M T, et al. Longer and more frequent marine heatwaves over the past century[J]. Nature Communications, 2018, 9: 1324. DOI:10.1038/s41467-018-03732-9 (0)
[12]
Hobday A J, Alexander L V, Perkins S E, et al. A hierarchical approach to defining marine heatwaves[J]. Progress in Oceanography, 2016, 141: 227-238. DOI:10.1016/j.pocean.2015.12.014 (0)
[13]
Lellouche J M, Greiner E, Bourdallé-Badie R, et al. The Copernicus global 1/12° oceanic and sea ice GLORYS12 reanalysis[J]. Frontiers in Earth Science, 2021, 9: 698876. DOI:10.3389/feart.2021.698876 (0)
[14]
栗春晓, 李芙蓉. 空间多观测样本的地理加权回归模型[J]. 中国海洋大学学报(自然科学版), 2024, 54(1): 156-164.
Li C X, Li F R. Geographically weighted regression for spatial data with replicates[J]. Periodical of Ocean University of China, 2024, 54(1): 156-164. DOI:10.16441/j.cnki.hdxb.20220201 (0)
[15]
Bian C, Jing Z, Wang H, et al. Oceanic mesoscale eddies as crucial drivers of global marine heatwaves[J]. Nature Communications, 2023, 14: 2970. DOI:10.1038/s41467-023-38811-z (0)
Research on Retrieving Statistical Characteristics of Global Subsurface Marine Heatwaves
Gong Zongqing , Li Furong     
School of Mathematical Science, Ocean University of China, Qingdao 266100, China
Abstract: This paper aims to retrieve global subsurface marine heatwaves (MHW) on a long-term scale, in order to characterize their trends over an extended period. Through preliminary data analysis, this paper found that there was a significant spatial correlation between the statistical characteristics of marine heatwaves in the subsurface layers and the monthly average seawater temperature. Guided by this statistical principle, the multiple-replicate geographically weighted regression (MRGWR) model was used as the retrieval model. The bandwidth was selected based on the characteristic that ocean circulation can be approximately considered constant within a 2° range. Furthermore, two indicator systems, namely univariate indicators and multivariate indicators, were constructed. The retrieval performance of the MRGWR model was compared with that of the generalized linear model (GLM), which was previously used for retrieving marine heatwaves, confirming the significant advantage of using multivariate indicators fitting to the MRGWR model in the retrieval of subsurface marine heatwaves. Using the selected optimal retrieval model and indicators, the statistical characteristics of global subsurface marine heatwaves from 1940 to 2021 were retrieved to characterize their long-term trends, which showed a clear correspondence with the rising trend of seawater temperature. This paper provides strong evidence for the impact of global warming caused by human activities on subsurface marine heatwaves.
Key words: subsurface marine heatwaves    physics-informed statistical learning    multiple-replicate geographically weighted regression model    generalized linear model