Search
E-mail
RSS

Online submission Editor-in-chief Peer review Staff Login

Journal of Tropical Oceanography >

2019 , Vol. 38 >Issue 5: 32 - 41

DOI: https://doi.org/10.11978/2018137

Marine Hydrography

Using artificial intelligence for identifying the depth of upper-ocean mixed layer

  • ZHANG Kang 1, 2 ,
  • GUO Shuangxi 1, 3 ,
  • HUANG Pengqi 1, 2 ,
  • QU Ling 1, 3 ,
  • LU Yuanzheng 1, 3 ,
  • CEN Xianrong 1, 3 ,
  • YU Lusha 1, 2 ,
  • ZHOU Weidong 1, 3 ,
  • ZHOU Shengqi , 1, 3
Expand
  • 1. State Key Laboratory of Tropical Oceanography (South China Sea Institute of Oceanology, Chinese Academy of Sciences), Guangzhou 510301, China
  • 2. University of Chinese Academy of Sciences, Beijing 100049, China
  • 3. Institution of South China Sea Ecology and Environmental Engineering, Chinese Academy of Sciences, Guangzhou 510301, China
ZHOU Shengqi. E-mail: .

Copy editor: YIN Bo

Received date: 2018-12-14

  Request revised date: 2019-04-10

  Online published: 2019-10-09

Supported by

National Natural Science Foundation of China(91752108)

National Natural Science Foundation of China(41476167)

National Natural Science Foundation of China(41706029)

National Natural Science Foundation of China(41606010)

Natural Science Foundation of Guangdong Province(2016A030311042)

Natural Science Foundation of Guangdong Province(2016A030310114)

Guangzhou Science and Technology Program Key Project(201804020056)

Strategic Priority Research Program of Chinese Academy of Sciences(XDA11030302)

Institution of South China Sea Ecology and Environmental Engineering, Chinese Academy of Sciences(ISEE2018PY05)

Copyright

Copyright reserved © 2019. Office of Acta Agronomica Sinica All articles published represent the opinions of the authors, and do not reflect the official policy of the Chinese Medical Association or the Editorial Board, unless this is clearly specified.
Fold

Abstract

An artificial intelligence (AI) method for identifying upper-ocean mixed layer depth (MLD) is proposed. In this method, a linear model, whose coefficient and variance are made into a set of statistics to characterize the profile, is established between temperature (density) and pressure (or depth). A subjective priori distribution is set in an initial window. The maximum posterior probability estimate of the mean coefficient value is obtained when the window is moved down from the sea surface, by learning the new data through the Bayesian chain rule and the minimum description length principle. The existence and depth of the mixed layer are determined when the jump of the coefficient is found by using the F-distribution. Using the Argo buoy data measured in the Pacific Ocean in February 2017, and taking the value of the quality index (QI) to estimate the accuracy of the MLD results, we find that this AI method is superior to the gradient method, the threshold method, the Hybrid method, and the relative-variant method.

Key words： upper-ocean mixed layer; artificial intelligence method; Bayesian chain rule; minimum description length principle

Cite this article

ZHANG Kang , GUO Shuangxi , HUANG Pengqi , QU Ling , LU Yuanzheng , CEN Xianrong , YU Lusha , ZHOU Weidong , ZHOU Shengqi . Using artificial intelligence for identifying the depth of upper-ocean mixed layer[J]. Journal of Tropical Oceanography, 2019 , 38(5) : 32 -41 . DOI: 10.11978/2018137

海洋上混合层(ocean upper mixed layer)是位于大气底边界层和深层海洋之间的温度、盐度、密度等在垂向上分布近似均匀的水体。它通常形成于海洋上层的湍流混合,在海洋和大气之间的淡水、动量以及热交换过程中起着重要作用(Kantha et al, 1999; Wunsch et al, 2004), 关于它的研究对建立平均气候态以及探索全球海洋变化极其重要(Deser et al, 1996; Hanawa et al, 2001; Yu et al, 2007; Oka et al, 2015)。混合层深度(Mixed Layer Depth, MLD)主要受到海表风应力、输入海洋的净热通量、蒸发降水等因素的影响, 其变化调节着海表温度(Sea Surface Temperature, SST)以及混合层热含量的演变等与热带气旋的形成、浮游植物的丰富程度以及气候变化等现象密切相关的过程(Carton et al, 2008; Dong et al, 2008; Vanin et al, 2009; Qiu et al, 2014)。风应力的剪切作用搅动海水, 加剧混合层内的湍流混合过程, 使MLD变深; 输入海洋的热通量与蒸发降水会改变海表的浮力通量, 进而改变MLD, 正的净热通量输入到上层海洋中会削弱湍流效应, 使MLD变浅, 反之造成一个更深的MLD; 蒸发(降水)导致海洋表层盐度变高(低), 使海表浮性频率降低(增加), 加强(减弱)湍流混合过程, 进而导致MLD的深度加深(变浅) (范聪慧, 2007; Talley et al, 2011; Qiu et al, 2016)。因而准确确定MLD对于建立海气热交换机制、研究海洋表层物理、生物过程以及气候变化等有着重要意义。
为了确定混合层深度, 前人做了大量的研究工作, 提出了一系列方法(Lukas et al, 1991; Brainerd et al, 1995; Obata et al, 1996; Kara et al, 2000; Lavender et al, 2002; Thomson et al, 2003; De Boyer Montégut et al, 2004; Lorbacher et al, 2006; Holte et al, 2009; Chu et al, 2011; Holte et al, 2017; Huang et al, 2018)。在这些方法中, 阈值类方法由于具有简单易于操作的特点成为最常用的一类。它们主要包括基于温度(密度)变化的阈值法、梯度变化的梯度法, 以及集合两者的混合法(Holte et al, 2009, 2017)。关于阈值的确定, Kara和Montegut等人在比较多种阈值的取值后, 总结出全球的温度和密度的最优阈值分别为0.2℃和0.03kg·m-3 (Kara et al, 2000; De Boyer Montégut et al, 2004; Holte et al, 2009)。在南大洋混合层的研究中, 密度和温度梯度的最优阈值分别被确定为0.0005~0.05kg·m-4和0.025℃·m-1 (Dong, 2008)。但是选择一个固定值作为最优阈值是比较主观的, 从大量温度和密度数据分析结果中发现, 不同区域(季节)需要不同的阈值, 其在热带与副热带区域为0.06kg·m-2, 在副极带区域为0.09kg·m-2 (Ohno et al, 2009)。另外, 选择一个合适的最优阈值也是比较困难的, 对于具有渐变密度跃层的廓线, 过小的阈值会使识别结果较实际情况偏浅, 而较大的阈值又会造成偏深的识别结果(Lorbacher et al, 2006)。
为了改进阈值法的这些不足, 学者提出了一些相对客观的方法。Thomson等(2003)对可变长度的密度廓线片段做线性最优拟合, 提出一种split-merge (分割-合并)方法, 但是他们发现这个方法在确定混合层深度的结果上和阈值法相近。Lorbacher等(2006)基于廓线的二阶导数提出了曲率法, 用全球海洋环流实验(World Ocean Circulation Experiment, WOCE)数据测试发现, 相比较于梯度阈值法, 该方法不仅在结果上更加精确, 而且可以更好地应用于不同分辨率的数据。然而该方法被验证难以应用于存在噪声的廓线数据中(Chu et al, 2011)。对水下滑翔机(seaglider)数据的研究中, Chu 等通过对温度(密度)廓线进行线性拟合提出了最优线性拟合法(Chu et al, 2010), 并且在此基础上进行改进提出了最大角度法(Chu et al, 2011), 分别对收集自西北大西洋和墨西哥湾流海域的滑翔机数据分析发现, 最优线性拟合法和最大角度法比阈值类方法具有更高的准确性。黄鹏起等人发现, 最大角度法和最优线性拟合法在分析WOCE数据中准确性不高, 他们通过构造廓线标准差和极差的比值—相对变量, 提出了相对变化法(Huang et al, 2018), 并且对WOCE数据进行了分析, 结果表明该方法比阈值类方法、曲率法和最大角度法等有着更好的表现。然而该方法仍存在些许不足, 在真实的海洋温度(位势密度)廓线中, 相对变量取得极大值的位置常常位于混合层下边界以下的某一深度上, 从该位置向上搜索混合层深度的过程中需要设定一些与廓线相关的参数。针对不同数据, 这些限定参数的最优取值是不同的。例如对于地转海洋学实时观测阵(Array for Real-time Geostrophic Oceanography, ARGO)数据, 限定参数的最优取值不同于其在WOCE数据中的取值。这些不足为使用该方法带来了不便, 因而提出一种既可以避开阈值类方法中最优阈值选取难题, 又可以弥补非阈值类方法的不足的方法是有意义的。
贝叶斯方法是Thomas Bayes在《An essay towards solving a problem in the doctrine of chances》中首次提出(Thomas Bayes, 1763), 并且经过De Fineti、Jefferys、Savage以及Lindley等人的努力, 逐渐发展出的一套完整的统计推断方法。该方法被广泛用于机器学习、人工智能等领域。本文基于贝叶斯方法, 提出一种识别海洋混合层深度的人工智能(Artificial Intelligence, AI)方法(下称AI法)。对温度(密度)和压强(或深度)建立线性模型, 并用模型系数和误差的方差因子描述温度(密度)及其梯度和涨落信息。在垂向梯度非常小的情况下, 把模型系数首次突变发生的位置识别为混合层下边界。使用2017年2月太平洋海域的ARGO浮标数据的测试结果表明新方法能够准确识别混合层深度, 与其他方法相比准确度更高, 并且在避开阈值类方法中最优阈值的选取难题的同时, 弥补了非阈值类方法的不足。

2 方法

2.1 模型的建立

设一条廓线具有n个有效数据点, 用n维列向量$\psi ={{\left( {{\psi }_{1}},{{\psi }_{2}},...,{{\psi }_{n}} \right)}^{\prime }}$和$P={{\left( {{p}_{1}},{{p}_{2}},...,{{p}_{n}} \right)}^{\prime }}$分别表示温度(位势密度)和对应的压强(或深度)。初始时, 为廓线设置一个长度为l的位于海表的压强(或深度)窗口。它将在数据学习过程中按照步长$\text{step}=1$向下移动。将下移k-1次后的窗口命名为wk, $k=1,2,...,n+1-l$, 并且把它内部的数据点记为${{S}_{k}}=\left( {{\psi }_{k}},{{P}_{k}} \right)$, 其中, ${{\psi }_{k}}={{\left( {{\psi }_{k}},{{\psi }_{k+1}},...,{{\psi }_{k+l-1}} \right)}^{\prime }}$ ${{P}_{k}}=$ ${{\left( {{p}_{k}},{{p}_{k+1}},\cdots ,{{p}_{k+l-1}} \right)}^{\prime }}$, 分别为窗口内的温度(位势密度)和压强(或深度)数据。
在窗口wk内建立线性模型(Koch, 2007):
$\begin{align}& {{\psi }_{k}}={{\beta }_{k,0}}\text{1}+{{\beta }_{k,1}}{{P}_{k}}+{{\varepsilon }_{k}} \\& {{\varepsilon }_{k}}\tilde{\ }N\left( 0,{{\tau }_{k}}^{-1}{{I}_{l\times l}} \right) \\\end{align}$
式中: ${{\mathbf{\varepsilon }}_{k}}={{\left( {{\varepsilon }_{k,1}},{{\varepsilon }_{k,2}},...,{{\varepsilon }_{k,l}} \right)}^{\prime }}$为窗口wk内的模型误差, 在此处认为它是正态的$N\left( 0,{{\tau }_{k}}^{-1}{{I}_{l\times l}} \right)$。Il×l是l×l的单位矩阵, 表示各点模型误差相互独立; 系数βk,0和βk,1蕴含温度(密度)及其垂向梯度信息, 方差因子τ-1描述数据的涨落信息, 他们可以被用于刻画数据的特征。记${{M}_{k}}=\left( 1,{{P}_{k}} \right)$, ${{\beta }_{k}}={{\left( {{\beta }_{k,0}},{{\beta }_{k,1}} \right)}^{\prime }}$。${{\tau }_{k}}^{-1}$和${{\beta }_{k}}$作为一组描述数据特征的统计量, 其先验概率密度函数$f(\tau,\beta)$是机器根据以往数据的信息, 给温度(密度)数据特征的预先认知。它的后验概率密度函数$f\left( {{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k}} \right)$ 是机器通过学习${{\psi }_{k}}$给出的对数据特征的新认知。它的似然函数
$f\left( {{\psi }_{k}}|{{\beta }_{k}},{{\tau }_{k}} \right)={{\left( 2\pi \right)}^{\frac{1}{2}}}{{\tau }^{\frac{1}{2}}}{{e}^{-\frac{\tau }{2}{{\left( {{\psi }_{k}}-{{M}_{k}}{{\beta }_{k}} \right)}^{\prime }}\left( {{\psi }_{k}}-{{M}_{k}}{{\beta }_{k}} \right)}}$
.表示在数据特征${{\tau }_{k}},{{\beta }_{k}}$取定的条件下, 温度(密度)取实测值的概率(Koch, 2007)。

2.2 数据学习过程

数据学习是更新机器对数数据学习是更新机器对数据特征的认知的过程。用正态-伽马分布${{N}_{\gamma }}\left( \mu ,V,b,\lambda \right)$描述对数据特征的认知, 其中$\mu ,V,b,\lambda $分别为该分布的均值、协方差、尺度参数和形状参数。这种分布具有共轭性, 即先验与后验概率密度函数都是正态-伽马分布(Degroot, 1970; Bernardo et al, 1994)。以窗口wk为例, τk和βk的先验和后验分布分别被记作${{\tau }_{k}},{{\beta }_{k}}\tilde{\ }{{N}_{\gamma }}\left( {{\mu }_{k}},{{V}_{k}},{{b}_{k}},{{\lambda }_{k}} \right)$和${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k}}\tilde{\ }{{N}_{\gamma }}\left( \mu _{k}^{*},V_{k}^{*},b_{k}^{*},\lambda _{k}^{*} \right)$。学习过程分为更新和校正两个步骤。
更新步骤从前一窗口的后验分布${{\tau }_{k-1}},{{\beta }_{k-1}}$ $|{{\psi }_{k-1}}\tilde{\ }{{N}_{\gamma }}\left( \mu _{k-1}^{*},V_{k-1}^{*},b_{k-1}^{*},\lambda _{k-1}^{*} \right)$, k>1出发, 此时机器对数据的预先认知是前一窗口中的后验信息。增加数据点$\left( {{\psi }_{k+l-1}},{{p}_{k+l-1}} \right)$后, 根据贝叶斯链式法则更新认知(见附录公式A1)据特征的认知的过程。
$\begin{align} & f\left( {{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k-1}},{{\psi }_{k+l-1}} \right)= \\ & \frac{f\left( {{\tau }_{k-1}},{{\beta }_{k-1}}|{{\psi }_{k-1}} \right)f\left( {{\psi }_{k+l-1}}|{{\tau }_{k-1}},{{\beta }_{k-1}} \right)}{\iint{f\left( {{\tau }_{k-1}},{{\beta }_{k-1}}|{{\psi }_{k-1}} \right)f\left( {{\psi }_{k+l-1}}|{{\tau }_{k-1}},{{\beta }_{k-1}} \right)\text{d}{{\tau }_{k-1}}\text{d}{{\beta }_{k-1}}}} \\ \end{align}$
式中: $f\left( {{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k-1}},{{\psi }_{k+l-1}} \right)$是机器对数据信息τk和βk新的认知, 把它记为${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k-1}},{{\psi }_{k+l-1}}\tilde{\ }$ ${{N}_{\gamma }}\left( \mu _{k}^{0},V_{k}^{0},b_{k}^{0},\lambda _{k}^{0} \right)$。
然而, 更新步骤得到的分布不能使模型(1)对窗口wk内数据的拟合达到最优化。校正过程则是从新的认知${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k-1}},{{\psi }_{k+l-1}}$出发得到使拟合最优化的后验分布的过程。用g-先验(Zellner, 1986; Liang et al, 2008)对先验分布${{\tau }_{k}},{{\beta }_{k}}\tilde{\ }{{N}_{\gamma }}\left( {{\mu }_{k}},{{V}_{k}},{{b}_{k}},{{\lambda }_{k}} \right)$调整如下(附录公式A2):
$\begin{align} & {{\mu }_{k}}=\mu _{k}^{0} \\ & {{V}_{k}}=gV_{k}^{0} \\ & {{\lambda }_{k}}=\frac{1}{2} \\ \end{align}$
根据贝叶斯公式, 后验概率${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k}}\tilde{\ }$ ${{N}_{\gamma }}\left( \mu _{k}^{*},V_{k}^{*},b_{k}^{*},\lambda _{k}^{*} \right)$为:
$f\left( {{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k}} \right)=\frac{f\left( {{\psi }_{k}}|{{\tau }_{k}},{{\beta }_{k}} \right)f\left( {{\tau }_{k}},{{\beta }_{k}} \right)}{\iint{f\left( {{\psi }_{k}}|{{\tau }_{k}},{{\beta }_{k}} \right)f\left( {{\tau }_{k}},{{\beta }_{k}} \right)\text{d}{{\tau }_{k}}\text{d}{{\beta }_{k}}}}$
调整过程中g和bk的取值, 可以由最小描述长度原理(Rissanen, 1978, 2009; Hansen et al, 2003)给出(附录公式A3、A4、A5)。数据学习过程中拟合直线如图1所示。
图1 数据学习过程示意图

图中粉色方框表示窗口

Fig. 1 An example of the learning process.

The blue line is the temperature profile. The magenta box represents the current window. The green, gray, and red dotted lines represent the fitting curves in the priori, the updating, and the posterior process, respectively

2.3 识别依据和检验方法

将${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k}}\tilde{\ }{{N}_{\gamma }}\left( \mu _{k}^{*},V_{k}^{*},b_{k}^{*},\lambda _{k}^{*} \right)$沿着τk积分, 得到窗口wk内模型系数的概率密度函数为:
${{\beta }_{k}}\tilde{\ }t\left( \mu _{k}^{*},\frac{b_{k}^{*}V_{k}^{*}}{\lambda _{k}^{*}},2\lambda _{k}^{*} \right)$
这个分布的均值为$\mu _{k}^{*}$。其分布如图2所示:
图2 窗口k内拟合温度及其梯度的概率密度函数

Fig. 2 Probability density function of temperature and its gradient

当k>1时, 模型系数发生突变的数学表达为:
${{\beta }_{k}}={{\beta }_{k-1}}$
式中, ${{\beta }_{k-1}}$和${{\beta }_{k}}$分.别为窗口k-1和k内数据拟合的模型系数。在直观上, 这种突变表现为相邻两个窗口内概率分布的不同。如图3所示, 颜色表示两个相邻窗口内所计算的温度及其梯度的概率密度。图3左侧颜色条带表示后一窗口内的概率分布, 右侧条带表示前一窗口的概率分布。当两个条带间距离比较大, 即相邻两窗口内概率分布不同时, 认为此处为混合层的下边界。
图3 相邻两窗口内数据分别拟合的温度及其梯度的概率密度图像

Fig. 3 Probability density image of temperature and its gradient of data in adjacent two windows

为了评估两条颜色带的距离大小, 我们引入一种归一化的距离。该距离的概率分布函数在数学上的意义为发生突变可能性, 刚好为自由度为2和$2\lambda _{k}^{*}$的F-分布:
$\frac{\lambda _{k-1}^{*}\left( {{\beta }_{k}}-{{\beta }_{k-1}} \right)'{{\left( V_{k-1}^{\mathbf{*}} \right)}^{-1}}\left( {{\beta }_{k}}-{{\beta }_{k-1}} \right)}{2b_{k-1}^{*}}\text{ }\!\!\tilde{\ }\!\!\text{ }F\left( 2,2\lambda _{k}^{*} \right)$
式中: $\lambda _{k-1}^{*}$、$b_{k-1}^{*}$和$V_{k-1}^{*}$为表征方差信息的参数, ${{\beta }_{k}}$和${{\beta }_{k-1}}$为相邻两窗口内数据拟合的模型系数。其图像如图4所示。
图4 F-分布的概率密度函数

该概率密度函数代表发生突变的可能性

Fig. 4 Probability density function of F-distribution

对图4所示的概率密度函数进行积分即可得到其分布函数。在检验水平0.05下,模型系数突变的条件为:
$\frac{\lambda _{k-1}^{*}\left( \mathbf{\mu }_{k-1}^{*}-\mathbf{\mu }_{k}^{*} \right)'{{\left( V_{k-1}^{\mathbf{*}} \right)}^{-1}}\left( \mathbf{\mu }_{k-1}^{*}-\mathbf{\mu }_{k}^{*} \right)}{2b_{k-1}^{*}}>{{F}_{0.95}}\left( 2,2\lambda _{k}^{*} \right)$
式中：${{F}_{0.95}}\left( 2,2\lambda _{k}^{*} \right)$ 突变发生的概率为0.95所对应的距离的值,$\mathbf{\mu }_{k-1}^{*}$和$\mathbf{\mu }_{k}^{*}$分别为${{\beta }_{k-1}}$和${{\beta }_{k}}$的均值的最大后验估计(附录公式A2.1)。此时识别突变点可信水平为0.95。突变检验过程如图5。
图5 突变检验示意图

a中实线为温度-深度廓线, 三角标注为突变点(混合层下边界); b中实线表示各点发生突变的概率, 虚线表示突变发生的临界概率, 突变发生概率首次超过临界值的点为突变点

Fig. 5 An example for checking the jump.

The solid line in (a) denotes the temperature-depth profile, and the triangle is the depth of the jump (the bottom of the mixed layer). The solid line in (b) indicates the probability of mutation at each point, and the dotted line represents the critical probability of the jump occurring. The jump will occur when the probability exceeds the critical value for the first time

最后给出识别混合层存在性的方法。计算出从参考压强(或深度10m) p0=100kPa到线性系数突变位置压强pjump的温度或位势密度差商D:
$D\text{=}\frac{\psi \left( {{p}_{\text{jump}}} \right)-\psi \left( {{p}_{0}} \right)}{{{p}_{\text{jump}}}-{{p}_{0}}}$
式中: $\psi \left( {{p}_{0}} \right)$表示参考深度的温度或位势密度, $\psi \left( {{p}_{\text{jump}}} \right)$表示突变位置的温度或位势密度。设初始的阈值δ=10-2, 如果D<3δ, 那么该廓线对应的水体中存在混合层结构。测试过程中, 记录已测试廓线中D的值, 并且用它的标准差更新δ的值。

3 ARGO数据测试结果

本节选取2017年2月太平洋的ARGO浮标数据对AI方法进行了测试。该组数据共包含6549个站位, 分布于如图6中黑色圆点所示的太平洋区域。
图6 太平洋ARGO浮标及2017年2月ARGO浮标站位图

灰色五角星位置为2017年2月1日第23号站位(28°57'43"N, 166°38'35"E)和第177号站位(57°22'55"S, 150°8'7"E)

Fig. 6 ARGO buoys and its stations in February 2017 of the Pacific

如图6中灰色五角星所示, 我们从2017年2月1日的218个站位中分别选取一个混合均匀(站位23)和一个存在弱层结的站位(站位177)。将两个站位的温度和密度廓线用AI法识别的结果分别绘制成图7和图8。图中用灰色三角表示识别的混合层深度。其中图7为混合层混合充分时的站位的识别结果, 图8为在混合层内存在层结时的识别结果。结果显示, 不论混合层内混合充分还是混合不够充分的情况下, AI法都能准确识别出混合层的深度。
图7 站位23 (28°57'43"N, 166°38'35"E)的温度(a)和密度(b)均匀混合层廓线和对应混合层深度

图中三角形表示混合层下边界位置

Fig. 7 Temperature (density) uniform mixed layer and recognition results.

The solid line represents the temperature profile (a) and potential density profile (b). The triangle is the lower boundary of the mixed layer

图8 站位177 (57°22'55"S, 150°8'7"E)的温度(a)和密度(b)存在弱层结的混合层廓线和对应混合层深度

图中三角形表示混合层下边界位置

Fig. 8 The temperature (density) profile with gradual pycnocline and recognition results.

The solid line represents the temperature profile (a) and potential density profile (b). The triangle represents the lower boundary of the mixed layer

对2017年2月太平洋的ARGO数据用AI法进行测试, 并且将所得混合层深度按照区域1°×1°范围内取中值进行汇总如图9和图10。两图中红色颜色的区域代表混合层深度较深, 蓝色区域代表该处混合层深度较浅。其中由温度廓线计算的混合层深度汇总如图9, 由位势密度廓线计算的混合层深度如图10。图中北半球混合层深度显著大于南半球混合层的深度, 这是由于2月的北半球刚好处于冬季, 混合层深度较深, 而南半球此时处于夏季, 混合层的深度较浅。与此同时, 图9和图10反映出中高纬地区的混合层深度比低纬地区的混合层深度深。
图9 AI法基于温度廓线识别的混合层深度的汇总图

Fig. 9 Mixed layer depth identified based on temperature profiles

图10 AI法基于密度廓线识别的混合层深度的汇总图

Fig. 10 Mixed layer depth identified based on density profiles

4 讨论

4.1 对比其他方法的识别结果

作为对比, 本文分别用相对变化法、阈值法、梯度法、混合方法、最优线性拟合法和最大角度法对相同的廓线计算混合层深度。QI值被广泛用作评估混合层深度识别准确程度的依据(Lorbacher et al, 2006):
$\text{QI}=1-\frac{{{S}_{\text{MLD}}}}{{{S}_{\text{1}\text{.5MLD}}}}$
式中: SMLD和S1.5MLD分别是所识别混合层深度和1.5倍深度内温度(位势密度)的标准差。一般而言, 较高的QI值代表较高的混合层深度识别准确度, 较低的QI值则代表混合层深度识别准确度较差。
对比7种识别混合层深度的方法, 我们发现使用AI法识别温度混合层深度的QI值的均值为0.936, 标准差为0.104, 1/4分位数、中位数和3/4分位数分别0.922、0.969和0.990; 识别位势密度混合层深度的QI值的均值为0.921, 标准差为0.117, 其1/4分位数, 中位数和3/4分位数分别为0.902、0.959和0.988。运用其他方法对于温度廓线识别的混合层深度的QI值的均值、标准差和四分位数详见表1, 对于密度廓线所得结果详见表2。
表1 各方法从温度廓线廓线计算的QI值的结果

Tab. 1 QI values calculated from the temperature profile by each method

方法 QImean QI0.25 QI0.5 QI0.75 QIstd
阈值法 0.795 0.708 0.899 0.962 0.239
梯度法 0.683 0.488 0.800 0.948 0.316
混合法 0.713 0.549 0.865 0.963 0.329
相对变化法 0.900 0.870 0.945 0.982 0.129
最大角度法 0.710 0.556 0.820 0.956 0.307
最优线性拟合 0.782 0.725 0.927 0.985 0.314
AI法 0.936 0.922 0.969 0.990 0.104

注: QImean为QI值的均值; QI0.25为QI值的1/4分位数; QI0.5为QI值的中值; QI0.75为QI值的3/4分位数; QIstd为QI值的标准差

表2 各方法从位势密度廓线中计算的QI值的结果

Tab. 2 QI values calculated from the potential density profile by each method

方法 QImean QI0.25 QI0.5 QI0.75 QIstd
阈值法 0.783 0.672 0.915 0.974 0.267
梯度法 0.727 0.596 0.808 0.926 0.251
混合法 0.734 0.594 0.872 0.962 0.306
相对变化法 0.893 0.862 0.945 0.983 0.144
最大角度法 0.740 0.621 0.850 0.965 0.299
最优线性拟合 0.797 0.752 0.930 0.985 0.297
AI法 0.921 0.902 0.959 0.988 0.117

注: QImean为QI值的均值; QI0.25为QI值的1/4分位数; QI0.5为QI值的中值; QI0.75为QI值的3/4分位数; QIstd为QI值的标准差

从表1和表2中可以看出AI法较其他方法识别结果的QI值具有更高的均值, 更小的标准差, 即具有更好的平均识别效果和更高的稳定性。同时, 其具有更高的0.25分位数和中位数, 即该方法在混合层廓线较为糟糕时识别效果优于其他6种方法, 而在0.75分位数上的表现发现, 在混合层混合较为充分时, AI法也有着轻微的优势。
针对温度廓线和密度廓线, 用各种方法所得到QI值的频率柱状图如图11。其中图11a表示从温度廓线中识别混合层深度的QI值的频率柱状图, 图11b表示从位势密度廓线中识别的混合层深度的QI值的频率柱状图。从图中看出, AI法较其他方法在QI值较高的区间(>0.8)占据更高的频率, 在QI值较低的区间(≤0.8)有着更低的频率。这些表明, AI法在识别混合层深度上具有更高的准确性。这是由于AI法可以识别混合层的存在性, 因而在QI的低值区与有着更小的频率。不考虑AI方法的话, 在西太平洋区域与南海北部近岸区域, 最优线性拟合法均为最优的混合层识别方法(Chu et al, 2010; Qiu et al, 2019)。本文引入AI方法之后, 发现在太平洋区域AI方法具有更高的适用性。
图11 7种方法识别混合层深度的质量因子(QI)的频率直方图

a. 从温度廓线识别混合层深度的质量因子的频率直方图; b. 从位势密度廓线识别混合层深度的质量因子的频率直方图

Fig. 11 Frequency histogram of five methods for QI values

4.2 概率分布的选取

AI方法的识别结果是不依赖于概率分布的选取的, 本文选择正态-伽马分布函数作为先验分布函数出于以下两个方面的考虑:
首先考虑的是计算的复杂性。正态-伽马分布${{N}_{\gamma }}\left( {{\mu }_{k}}\mathbf{,}{{V}_{k}},{{b}_{k}},{{\lambda }_{k}} \right)$具有一个优良的属性—共轭性, 即先验分布与后验分布具有相同的分布形式, 因此计算过程只需要待定${{\mu }_{k}},{{V}_{k}},{{b}_{k}},{{\lambda }_{k}}$的取值即可。事实上,温度(位势密度) -压强(或深度)拟合模型中, 对数据信息认知的更新过程和后验分布的校正过程都涉及多重积分运算, 计算开销比较大。选取正态-伽马分布可以避免大量的积分运算, 从而有效降低计算的时间复杂性。与此同时, 计算过程中存在两个非常消耗计算量的过程。一个是更新过程中窗口向下移动一个点, 对于新数据点概率分布的估计过程,
$\begin{align} & f\left( {{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k-1}},{{\psi }_{k+l-1}} \right)= \\ & \frac{f\left( {{\tau }_{k-1}},{{\beta }_{k-1}}|{{\psi }_{k-1}} \right)f\left( {{\psi }_{k+l-1}}|{{\psi }_{k-1}} \right)}{\iint{f\left( {{\tau }_{k-1}},{{\beta }_{k-1}}|{{\psi }_{k-1}} \right)f\left( {{\psi }_{k+l-1}}|{{\psi }_{k-1}} \right)\text{d}{{\tau }_{k-1}}\text{d}{{\beta }_{k-1}}}} \\ \end{align}$
选取正态-伽马分布族后, 无需迭代即可求出公式(12)中
$\begin{align} & f\left( {{\psi }_{k+l-1}}|{{\psi }_{k-1}} \right)= \\ & \int{f\left( \tau ',\beta '|{{\psi }_{k-1}} \right)f\left( {{\psi }_{k+l-1}}|\tau ',\beta ' \right)}\text{d}\tau '\text{d}\beta ' \\ \end{align}$
的表达式
$f\left( {{\psi }_{k+l-1}}|{{\psi }_{k-1}} \right)=f\left( {{\psi }_{k+l-1}}|{{\tau }_{k-1}},{{\beta }_{k-1}} \right)$
此时公式(12)就成了公式(3)。另一个过程是在应用最小描述长度(Minimum Description Length, MDL)准则过程中, 将公式(A3)所得的描述长度极小化。
选取正态-伽马分布将确定$f\left( {{\tau }_{k}},{{\beta }_{k}} \right)$的最优估计转化为确定${{\mu }_{k}},{{V}_{k}},{{b}_{k}},{{\lambda }_{k}}$的最优估计, 由于正态-伽马分布的共轭性, 4个参数中只有bk需要迭代求解, 且迭代过程可以迅速收敛(Rissanen, 1989), 从而极大减少计算开销。综上所述, 选择正态-伽马分布是为了降低计算量, 程序中只有一次循环与n有关, 其加法运算次数正比于nl2, 给定窗口长度l后, 其时间复杂度为O(n)。
选择正态-伽马分布的另一个理由是与模型残差的正态假设相匹配。在模型中, 拟合的误差是服从正态分布的。所以, 方差因子τk-1未知的情况下, 其无偏估计为:
${{\tau }_{k}}^{-1}\text{ }\!\!\tilde{\ }\!\!\text{ }{{\chi }^{2}}\left( l-1 \right)$
式中: l仍表示窗口长度, 此时τ, β服从正态-伽马分布, 因而本文选择正态-伽马分布是为了与模型中的正态假设匹配。

5 结论

本文基于贝叶斯方法提出一种识别海洋上混合层深度的人工智能识别方法。该方法将廓线特征的变化投影到概率空间上, 再通过真实数据学习最优的概率分布函数。相较于阈值类方法(如阈值法、梯度法、混合法)的直接使用阈值, 该方法使用统一的检验水平, 从而避开了最优阈值选取困难的问题。相较于非阈值类方法(例如相对变化法、最优线性拟合法和最大角度法), 该方法充分考虑数据的随机效应, 对未知的信息做等概率处理(最小描述长度原理), 以此取代额外设定与廓线信息相关的限定参数, 弥补了非阈值类方法的不足。通过对ARGO浮标数据测试的结果显示, 该方法在混合充分以及具有层结的情况下均可以准确地识别混合层深度。与阈值法、梯度法、混合法、相对变化法、最优线性插值法和最大角度法进行比较, 并且用QI值作为识别混合层深度的准确性的评估依据, 发现AI法较其它方法结果更加准确。与此同时, 该方法还可以对混合层的存在性进行识别, 提高了混合层深度识别的精度, 有望将来在全球区域适用来提高气候模拟的准确性。

附录

数据学习过程的迭代公式:
以窗口wk为例, 更新过程对数据的认知${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k-1}},{{\psi }_{k+l-1}}\text{ }\!\!\tilde{\ }\!\!\text{ }{{N}_{\gamma }}\left( \mu _{k}^{0},V_{_{k}}^{0},b_{_{k}}^{0},\lambda _{k}^{0} \right)$由下面的公式给出:
$\mu _{k}^{0}=V_{k}^{0}\left( \left( \begin{matrix} 1 \\ {{p}_{k+1}} \\\end{matrix} \right){{\psi }_{k+1}}+{{\left( V_{k-1}^{*} \right)}^{-1}}\mu _{k-1}^{*} \right)$
$V_{k}^{0}={{\left( {{\left( V_{k-1}^{*} \right)}^{-1}}+\left( \begin{matrix} 1 \\ {{p}_{k+1}} \\\end{matrix} \right)\left( 1,{{p}_{k+1}} \right) \right)}^{-1}}$
$\begin{align} & b_{k}^{0}=b_{k-1}^{*}+\frac{{{\left( \mu _{k}^{0}-\mu _{k-1}^{*} \right)}^{\prime }}{{\left( V_{k-1}^{*} \right)}^{-1}}\left( \mu _{k}^{0}-\mu _{k-1}^{*} \right)}{2}\text{+} \\ & \ \ \ \ \ \ \frac{{{\left( {{\psi }_{k+1}}-\left( 1,{{p}_{k+1}} \right)\mu _{k}^{0} \right)}^{\prime }}\left( {{\psi }_{k+1}}-\left( 1,{{p}_{k+1}} \right)\mu _{k}^{0} \right)}{2} \\ \end{align}$
$\lambda _{k}^{0}=\frac{1+2\lambda _{k-1}^{*}}{2}$
采用g-先验分布调整的后验分布为${{\tau }_{k}},{{\beta }_{k}}|{{\psi }_{k}}\sim {{N}_{\gamma }}\left( \mu _{k}^{*},V_{_{k}}^{*},b_{_{k}}^{*},\lambda _{k}^{*} \right)$。其中:
$\mu _{k}^{*}=V_{k}^{*}\left( {{M}_{k}}^{\prime }{{\psi }_{k}}+{{g}^{-1}}{{\left( V_{k}^{0} \right)}^{-1}}{{\mu }_{k}} \right)$
$V_{k}^{*}=\left( {{g}^{-1}}{{\left( V_{k}^{0} \right)}^{-1}}+{{M}_{k}}^{\prime }{{M}_{k}} \right)$
$\begin{align} & b_{k}^{*}={{b}_{k}}+\frac{{{\left( \mu _{k}^{*}-{{\mu }_{k}} \right)}^{\prime }}{{\left( V_{k}^{0} \right)}^{-1}}\left( \mu _{k}^{*}-{{\mu }_{k}} \right)}{2g}+ \\ & \ \frac{{{\left( {{\psi }_{k}}-{{M}_{k}}\mu _{k}^{*} \right)}^{\prime }}\left( {{\psi }_{k}}-{{M}_{k}}\mu _{k}^{*} \right)}{2} \end{align}$
$\lambda _{\text{k}}^{\text{*}}=\frac{l+2}{2}$
参数和g通过DML准则确定, Hansen (2003)给出了一种混合MDL的表达式:
${{H}_{\text{MDL}}}=\min \left( \frac{1}{2}\ln \left| gV_{k}^{0} \right|-\frac{1}{2}\ln \left| V_{k}^{*} \right|-{{\lambda }_{k}}\ln \left( 2{{b}_{k}} \right)+\lambda _{k}^{*}\ln \left( 2b_{k}^{*} \right) \right)$
参数g和bk的最优估计为:
${{\hat{b}}_{k}}=\frac{{{R}_{g}}}{l}$
$\hat{g}=\frac{{{R}_{{\hat{g}}}}\text{trace}\left( {{\left( V_{k}^{0} \right)}^{-1}}{{\left( {{{\hat{g}}}^{-1}}{{\left( V_{k}^{0} \right)}^{-1}}+{{M}_{k}}^{\prime }{{M}_{k}} \right)}^{-1}} \right)+lR}{2{{R}_{{\hat{g}}}}}$
$\begin{align} & R={{\mu }_{k}}^{\prime }{{\left( V_{k}^{0} \right)}^{-1}}{{\mu }_{k}}+{{\psi }_{k}}^{\prime }{{M}_{k}}\left( {{\left( {{{\hat{g}}}^{-1}}{{\left( V_{k}^{0} \right)}^{-1}}+{{M}_{k}}^{\prime }{{M}_{k}} \right)}^{-1}} \right) \\ & \ \ \ \ \ \ {{\left( V_{k}^{0} \right)}^{-1}}\left( {{\left( {{{\hat{g}}}^{-1}}{{\left( V_{k}^{0} \right)}^{-1}}+{{M}_{k}}^{\prime }{{M}_{k}} \right)}^{-1}} \right){{M}_{k}}^{\prime }{{\psi }_{k}}- \\ & \ \ \ \ \ \ 2{{\left( {{\left( V_{k}^{0} \right)}^{-1}}{{\mu }_{k}} \right)}^{\prime }}{{\left( \hat{g}+{{M}_{k}}^{\prime }{{M}_{k}} \right)}^{-1}}\left( \hat{g}+{{M}_{k}}^{\prime }{{\psi }_{k}} \right) \\ \end{align}$
$\begin{align} & {{R}_{{\hat{g}}}}={{\psi }_{k}}^{\prime }{{\psi }_{k}}+{{\mu }_{k}}^{\prime }{{\left( \hat{g}V_{k}^{0} \right)}^{-1}}{{\mu }_{k}}- \\ & \ \ \ \ \ \ \ {{\psi }_{k}}^{\prime }{{M}_{k}}{{\left( {{{\hat{g}}}^{-1}}{{\left( V_{k}^{0} \right)}^{-1}}+{{M}_{k}}^{\prime }{{M}_{k}} \right)}^{-1}}{{M}_{k}}^{\prime }{{\psi }_{k}} \\ \end{align}$
对于隐式方程(A4.2), 用如下收敛的迭代格式求其近似解:
$\begin{align} & {{{\hat{g}}}_{1}}=0 \\ & {{{\hat{g}}}_{j+1}}=\frac{{{R}_{{{{\hat{g}}}_{j}}}}trace\left( {{\left( V_{k}^{0} \right)}^{-1}}{{\left( {{\left( {{{\hat{g}}}_{j}}V_{k}^{0} \right)}^{-1}}+{{M}_{k}}^{\prime }{{M}_{k}} \right)}^{-1}} \right)+lR}{2{{R}_{{{{\hat{g}}}_{j}}}}} \\ \end{align}$
式中j为迭代次数, 当$\left| {{{\hat{g}}}_{j+1}}-{{{\hat{g}}}_{j}} \right|<{{10}^{-6}}$ 时, 退出迭代过程, $\hat{g}={{\hat{g}}_{j+1}}$。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
范聪慧 , 2007. 多因素对海洋上混合层深度影响的数值模拟[D]. 北京: 中国科学院研究生院( 海洋研究所):1-65.

FAN CONGHUI , 2007. Numerical simulation of the ocean surface mixed layer depth influenced by the multi-factors[D]. Beijing: The Institute of Oceanology, Chinese Academy of Sciences (in Chinese with English abstract).

[2]
BAYES T , 1763. An essay towards solving a problem in the doctrine of chances[J]. Philosophical Transactions of the Royal Society of London, 53:370-418.

[3]
BERNARDO J M, SMITH A F M , 1994. Bayesian theory[M]. New York: Wiley:265-268.

[4]
BRAINERD K E, GREGG M C , 1995. Surface mixed and mixing layer depths[J]. Deep Sea Research Part I: Oceanographic Research Papers, 42(9):1521-1543.

[5]
CARTON J A, GRODSKY S A, LIU HAILONG , 2008. Variability of the oceanic mixed layer, 1960-2004[J]. Journal of Climate, 21(5):1029-1047.

[6]
CHU P C, FAN CHENWU , 2010. Optimal linear fitting for objective determination of ocean mixed layer depth from glider profiles[J]. Journal of Atmospheric and Oceanic Technology, 27(11):1893-1898.

[7]
CHU P C, FAN CHENWU , 2011. Maximum angle method for determining mixed layer depth from seaglider data[J]. Journal of Oceanography, 67(2):219-230.

[8]
DE BOYER MONTÉGUT C, MADEC G, FISCHER A S , et al, 2004. Mixed layer depth over the global ocean: an examination of profile data and a profile-based climatology[J]. Journal of Geophysical Research: Oceans, 109(C12):C12003.

[9]
DEGROOT M H , 1970. Optimal statistical decisions[M]. New York: McGraw- Hill:159-172.

[10]
DESER C, ALEXANDER M A, TIMLIN M S , 1996. Upper-ocean thermal variations in the north pacific during 1970-1991[J]. Journal of Climate, 9(8):1840-1855.

[11]
DONG SHENFU, SPRINTALL J, GILLE S T , et al, 2008. Southern ocean mixed-layer depth from Argo float profiles[J]. Journal of Geophysical Research: Oceans, 113(C6):C06013.

[12]
HANAWA K, TALLEY L D , 2001. Chapter 5.4 Mode waters[J]. International Geophysics, 77:373-386.

[13]
HANSEN M H, YU BIN , 2003. Minimum description length model selection criteria for generalized linear models[M] //GOLDSTEIN D R. Statistics and science: a festschrift for terry speed. U.S.: Institute of Mathematical Statistics: 145-163.

[14]
HOLTE J, TALLEY L , 2009. A new algorithm for finding mixed layer depths with applications to Argo data and subantarctic mode water formation[J]. Journal of Atmospheric and Oceanic Technology, 26(9):1920-1939.

[15]
HOLTE J, TALLEY L D, GILSON J , et al, 2017. An Argo mixed layer climatology and database[J]. Geophysical Research Letters, 44(11):5618-5626.

[16]
HUANG PENGQI, LU YUANZHENG, ZHOU SHENGQI , 2018. An objective method for determining ocean mixed layer depth with applications to WOCE data[J]. Journal of Atmospheric and Oceanic Technology, 35(3):441-458.

[17]
KANTHA L, CLAYSON C A , 1999. Surface density changes and the ocean mixed layer[M] //CURRY J A, WEBSTER P J. Thermodynamics of atmospheres and oceans. San Diego, U.S: Academic Press: 291-298.

[18]
KARA A B, ROCHFORD P A, HURLBURT H E , 2000. An optimal definition for ocean mixed layer depth[J]. Journal of Geophysical Research: Oceans, 105(C7):16803-16821.

[19]
KOCH K R , 2007. Introduction to Bayesian statistics[M]. New York: Springer:85-121.

[20]
LAVENDER K L, DAVIS R E, OWENS W B , 2002. Observations of open-ocean deep convection in the Labrador Sea from subsurface floats[J]. Journal of Physical Oceanography, 32(2):511-526.

[21]
LIANG FENG, PAULO R, MOLINA G , et al, 2008. Mixtures of g priors for Bayesian variable Selection[J]. Journal of the American Statistical Association, 103(481):410-423.

[22]
LORBACHER K, DOMMENGET D, NIILER P P , et al, 2006. Ocean mixed layer depth: a subsurface proxy of ocean- atmosphere variability[J]. Journal of Geophysical Research: Oceans, 111(C7):C07010.

[23]
LUKAS R, LINDSTROM E , 1991. The mixed layer of the western equatorial pacific ocean[J]. Journal of Geophysical Research: Oceans, 96(S01):3343-3357.

[24]
OBATA A, ISHIZAKA J, ENDOH M , 1996. Global verification of critical depth theory for phytoplankton bloom with climatological in situ temperature and satellite ocean color data[J]. Journal of Geophysical Research: Oceans, 101(C9):20657-20667.

[25]
OHNO Y, IWASAKA N, KOBASHI F , et al, 2009. Mixed layer depth climatology of the north pacific based on Argo observations[J]. Journal of Oceanography, 65(1):1-16.

[26]
OKA E, QIU BO, TAKATANI Y , et al, 2015. Decadal variability of subtropical mode water subduction and its impact on biogeochemistry[J]. Journal of Oceanography, 71(4):389-400.

[27]
QIU CHUNHUA, CUI YONGSHENG, REN JIE , et al, 2016. Characteristics of the surface mixed layer depths in the northern South China Sea in spring[J]. Journal of Oceanography, 72(4):567-576.

[28]
QIU CHUNHUA, HUO DAN, LIU CHANGJIAN , et al, 2019. Upper vertical structures and mixed layer depth in the shelf of the northern South China Sea[J]. Continental Shelf Research, 174:26-34.

[29]
QIU CHUNHUA, KAWAMURA H, MAO HUABIN , et al, 2014. Mechanisms of the disappearance of sea surface temperature fronts in the subtropical North Pacific Ocean[J]. Journal of Geophysical Research: Oceans, 119(7):4389-4398.

[30]
RISSANEN J , 1978. Modeling by shortest data description[J]. Automatica, 14(5):465-471.

[31]
RISSANEN J , 1989. Stochastic complexity in statistical inquiry[M]. London: World Scientific: 122-130.

[32]
RISSANEN J , 2009. Model selection and testing by the MDL principle[M] //EMMERT-STREIB F, DEHMER M. Information theory and statistical learning. New York, U.S.: Springer:25-43.

[33]
TALLEY L D, PICKARD G L, EMERY W J , et al, 2011. Descriptive physical oceanography: an introduction[M]. 6th ed. London: Academic Press: 194-195.

[34]
THOMSON R E, FINE I V , 2003. Estimating mixed layer depth from oceanic profile data[J]. Journal of Atmospheric and Oceanic Technology, 20(2):319-329.

[35]
VANIN N S, KHEN G V , 2009. Vertical structure of water masses and silicon-phosphorus ratios in the west Bering Sea and in the Sea of Okhotsk[J]. Oceanology, 49(3):350-360.

[36]
WUNSCH C, FERRARI R , 2004. Vertical mixing, energy, and the general circulation of the oceans[J]. Annual Review of Fluid Mechanics, 36:281-314.

[37]
YU LISAN, WELLER R A , 2007. Objectively analyzed air-sea heat fluxes for the global ice-free oceans (1981-2005)[J]. Bulletin of the American Meteorological Society, 88(4):527-540.

[38]
ZELLNER A , 1986. On assessing prior distributions and Bayesian regression analysis with G-prior distributions[M] //GOEL P K, ZELLNER A. Bayesian inference and decision techniques. New York, U.S.: Elsevier Science Publishers: 233-243.

Options
Outlines

/

〈 〉