海洋中的溶解氧为海洋生物生存提供了必要的生化环境, 其时空分布同时受到物理过程与生物化学过程的影响, 物理过程包括海洋通风、海气交换、混合作用, 生物化学过程则包括初级生产、呼吸作用、有机物氧化作用。通过海水溶解氧的观测与分析, 可以在一定程度上了解该海区的生物地球化学循环过程及其对不同时空尺度物理过程的响应。特别是在全球变暖的大背景下, 上层海洋层化逐渐加强, 海洋中溶解氧的含量呈现降低的趋势, 开阔大洋的脱氧现象(deoxygenation)不断加剧, 海洋中层最低含氧区(oxygen minimum zone, OMZ)的范围不断扩大(Stramma et al, 2008; Keeling et al, 2010)。贫氧(anoxic)环境产生反硝化作用(denitrification)会生成比CO₂更强烈的温室气体——N₂O, 同时(由于氧气不足)硝酸盐作为呼吸作用的替代性受体会被快速消耗, 这对于海洋碳、氮循环以及海洋生态系统都具有重要的影响, 同时对全球变暖产生正反馈(Keeling et al, 2010; Johnson et al, 2010)。因此, 海水溶解氧观测对于海洋生态系统与气候变化的研究至关重要。

直到21世纪初, 开阔大洋的溶解氧观测数据一直非常缺乏(Stramma et al, 2008), 这极大地限制了人们对全球海洋溶解氧时空分布及其受物理过程影响的认识与理解。Argo浮标的出现为海洋学的发展做出了革命性的贡献(Roemmich et al, 2009), 借助传统温盐Argo浮标成熟的技术平台, 自2002年起, 携带溶解氧传感器的Argo浮标开始出现(Körtzinger et al, 2005)。截至2018年7月, 全球投放的可观测溶解氧的生物地球化学Argo(Biogeochemical-Argo)浮标已超过900台(图 1), 总计观测了超过13万条溶解氧剖面。其观测数据在海洋学研究中发挥着重要的作用, 包括海洋净群落生产力(Net Community Production, NCP)的季节循环(Riser et al, 2008; Yang et al, 2017; Huang et al, 2018)、海洋碳输出(Martz et al, 2008)、海气通量(Kihm et al, 2010; Bushinsky et al, 2017)、以及最低含氧区的动态变化(Prakash et al, 2012; Stanev et al, 2018; Sarma et al, 2018)等研究。

然而到目前为止, 这一海量的数据集还没有得到充分的利用, 其数据的质量控制研究仍在不断完善。影响溶解氧数据质量的主要问题来源于传感器的“存放漂移”(Johnson et al, 2015), 表现为溶解氧测量值与真实值之间存在一个固定的斜率。针对其漂移的校正目前有两种方法:空气测量校正法(Bittig et al, 2018)与气候态校正法(Takeshita et al, 2003)。虽然空气测量法被认为是目前最准确的Argo溶解氧校正方法, 但也存在一定问题(刘烨等, 2019), 主要表现为: 1)各种类型的Argo浮标目前空气测量方式与校正方法并不统一, 导致不同类型浮标的空气校正结果可能存在系统性偏差; 2)2016年之前投放的溶解氧Argo浮标由于传感器的高度不够, 可能没有有效观测到空气中的氧含量; 3)此外, 也是更为重要的是, 在目前全球Argo溶解氧的数据库中(截至2018年7月), 超过70%浮标(684台)无法进行空气测量, 其观测的剖面数超过80%(110677条), 如果仅采用空气测量校正法则会损失大量的观测数据。

另一方面, 气候态校正方法不仅可以校正所有的溶解氧数据, 而且使用的数据标准是统一的气候态数据集, 因而不同类型的浮标数据校正后不会存在系统性偏差。然而, 目前的气候态校正方法也存在一些问题, 主要表现为其校正系数会受到异常剖面和异常值的影响(刘烨等, 2019)。在此背景下, 本文提出了一种基于气候态校正方法的改进方案, 以期进一步提高校正系数的准确性以及溶解氧的数据质量。

1 材料与方法 1.1 数据

Argo浮标数据来源于法国Argo数据中心(ftp://ftp.ifremer.fr/ifremer/argo), 截至2018年7月共计携带有溶解氧传感器的Argo浮标926台(图 1), 总计观测了137714条剖面数据。由于溶解氧数据的质量控制通常使用氧分压p_O2或溶解氧饱和度Sat_O2进行校正, 浮标观测的溶解氧浓度([O₂])需要转化为饱和度:

其中, 氧气的溶解度(solubility)[O_{2, sol}]根据海水温度、盐度计算(Bittig et al, 2016)。

溶解氧气候态数据集使用美国国家海洋与大气局(National Oceanic and Atmospheric Administration, NOAA)下属的国家海洋数据中心发布的世界海洋地图集(World Ocean Atlas, WOA)13v2(Garcia et al, 2013)(www.nodc.noaa.gov/OC5/woa13):气候态逐月的温度(T)、盐度(S)、溶解氧饱和度(Sat_O2)剖面数据, 水平分辨率1°×1°。

数据处理过程中, 需要首先剔除掉一些由于观测剖面较少无法有效获得校正系数的浮标, 根据Bittig等(2015)的研究表明, 至少需要20个有效的m_i值可以获得较为可信的校正系数, 因此本文通过20个有效值的标准首先对数据进行筛选。在926台浮标中, 83台浮标观测的剖面数不足20个、55台浮标虽然剖面数满足标准但观测到的表层值(深度 < 10m)不足20个、30台浮标具有足够的浮标表层值但对应的WOA气候态表层值不足20个、4台浮标虽然有足够的浮标与WOA表层值但两者匹配数据不足20个。最终满足条件的浮标共计754台(占全部浮标个数的81%), 共计观测127611条剖面(占全部剖面个数的93%)。

1.2 Takeshita气候态校正方法

Takeshita等(2013)首先提出了基于WOA气候态数据的溶解氧校正方法, 使用气候态月平均的海表面溶解氧饱和度Sat_{O2, WOA, i}(经过时间与空间内插到浮标所在时间与位置), 与浮标观测的海表面饱和度(Sat_{O2, float, i})比值, 每个观测剖面得到一个斜率系数m_i(公式(2)), 收集该浮标观测的所有斜率系数m_i取平均得到校正系数M, 本文记为M₀, (公式(3))。这里, N代表该浮标观测所得m_i的个数。

1.3 校正方法改进1:中值法替代平均值法

由于目前的校正系数易受到异常剖面和异常值的影响(刘烨等, 2019), 本文针对气候态校正方法提出了三种改进策略, 层层推进。为了比较不同改进策略的效果, 分别将改进后的校正系数命名为M₁、M₂和M₃。首先使用中值法替代平均值法:即对于一台浮标观测的所有m_i值按大小排序后取中值作为校正系数, 记为M₁。

1.4 校正方法改进2:异常溶解氧剖面剔除标准

目前的气候态校正方法(Takeshita et al, 2013)在计算每个剖面的斜率系数m_i时, 只比较每个剖面的表层溶解氧饱和度, 对于每个m_i值并没有提供任何评估方法。而实际上, 浮标现场观测的溶解氧饱和度除近表面受生物活动的影响外, 整体与气候态平均值的相似程度很高(刘烨等, 2019), 因而可以首先对浮标观测的溶解氧剖面与WOA气候态剖面进行比较, 由此剔除掉一些异常剖面, 这些异常剖面计算的m_i值则不会对最终计算的校正系数M产生影响, 从而提高M值的准确性。本文引入两个剖面(浮标观测与WOA)的平均绝对误差(mean absolute error, MAE)作为评估标准(公式(4)), 通过结果分析(见后文), 选择10%作为判别阈值, 符合这一阈值的m_i记为m_i*(即MAE_i≤10%的m_i值), 再通过中值法得到校正系数, 记为M₂。其中, Sat_{O2, floatprof, i}和Sat_{O2, WOAprof, i}分别代表某一条浮标观测的溶解氧饱和度剖面与对应的WOA溶解氧饱和度剖面(由于WOA的观测深度有限, 将浮标观测的溶解氧剖面内插到WOA剖面对应深度上进行计算), N_{obs, i}代表该剖面上点的个数。

1.5 校正方法改进3:异常斜率系数的剔除标准

通过第二步的异常剖面判别标准, 可以剔除掉一些影响最终计算校正系数的m_i值, 但在实际分析中发现, 对于某一台浮标的全部m_i值, 仍然可能存在离散性很大的情况。基于此, 本文提出改进方法的第三步, 在异常剖面判别的基础上, 引入所有m_i*值的绝对中值偏差(median absolute deviation, MAD), 它定义为所有m_i*值的中值与m_i*值差的绝对值的中值。并以三倍绝对中值偏差(3×MAD)作为偶然误差的容许值, 以此剔除掉离散性较大的m_i值(Leys et al, 2013)。在实际计算过程中, 由于每次有离散程度较大的m_i值被剔除后, 其MAD值会发生改变, 因而采用迭代方法:当有m_i值被剔除(即m_i值的个数发生变化)后, 重新计算MAD, 直到没有m_i值被剔除(即m_i值的个数不发生变化)时, 此时完成迭代, 计算剩余所有m_i值的中值, 记为M₃。

2 结果与讨论

由于“存放漂移”的存在, 理论上传感器的校正系数M值通常大于1, 但在实际计算中发现, 仍然存在少量传感器会出现M < 1的情况(Johnson et al, 2015)。整体上, 溶解氧传感器的漂移程度一般比较小, Johnson等(2015)的结果显示, M值一般出现在0.83—1.20。本文选取了m_i变化范围较大的四台Argo浮标作为范例(编号分别为2900731、3900347、2901550和2900765), 图 2显示了这四台浮标基于气候态校正方法计算的斜率系数m_i的时间序列, 其中浮标2900731(图 2a)的m_i时间序列表现为:在观测的前期(2008年)和后期(2010年)m_i较为稳定, 而在观测中期出现一个显著的高值区(2009年), 最高达到4.5, 这一高值区可能是由于现场观测数据显著低于气候态数据导致的, 但与传感器自身的稳定性无关; 浮标3900347(图 2b)则表现为:在浮标运行的大部分时间内m_i值都稳定在1.1左右, 只在后期(2008年夏季开始)突然变高, 最大值接近11.0, 增长了接近10倍, 这一异常的增高可能与传感器的退化或生物附着有关; 浮标2901550(图 2c)表现为:浮标观测前期m_i值基本保持不变, 但自2014年夏季开始, m_i值突然变化, 先突然减小而后出现较大波动, 自2015年春季之后则迅速增加, 最终达到接近7.0;浮标2900765(图 2d)同样表现为观测前期m_i值较为稳定, 之后开始增大, 与浮标3900347(图 2b)不同的是, 此台浮标后期的异常高值较多, 为准确计算校正系数带来很大的困难。

基于传统的气候态校正方法, 计算得到四台浮标的M₀值分别为1.5993、1.6035、0.88032以及1.2883, 其中前两台浮标的计算结果显著偏高, 后两台浮标虽然在数值上符合M值的正常范围, 但偏离该浮标前期稳定的m_i值, 分别出现一定程度的低估(图 2c)和高估(图 2d)。

2.1 中值法的改进

在传感器测量到的溶解氧剖面大部分为正常剖面时, 采用平均值法也可以得到较为准确的校正系数。但当异常剖面较多时, 平均值会显著地受到这些偏大(或偏小)的异常值影响, 当这些异常值的数量较少(或分布均匀), 中值则较少受到影响, 具有较强的稳健性。图 2显示出, 对于前两台浮标(图 2a和图 2b), 采用中值法之后, M₁值与大部分m_i值非常接近, 分别为1.2352和1.1881, 也符合M值的正常范围(Johnson et al, 2015); 对于后两台浮标(图 2c和图 2d), M₁值虽然与大部分m_i值仍然存在差距, 但都比M₀更接近。

通过两个校正系数校正后的溶解氧剖面与WOA气候态的比较也可以发现, 通过M₁值校正后平均绝对偏差(MAE)在10%以上的剖面数量明显减少, 说明中值法计算的M₁值受到少量异常剖面的影响较小, 更接近大部分正常剖面计算的m_i值, 因而大幅度地提高了Argo浮标观测溶解氧的整体数据质量。

2.2 异常溶解氧剖面剔除标准的改进

第二步校正方法的改进通过比较浮标观测的溶解氧剖面与WOA剖面的相似性, 剔除掉一些异常的溶解氧剖面, 首先需要确定剔除标准。图 4显示了所有Argo浮标观测的溶解氧饱和度剖面与WOA气候态剖面的MAE的分布, 以及与温度和盐度剖面的平均绝对偏差(MAE_T和MAE_S)的散点图。可以看出, 大部分(90.42%)浮标观测的溶解氧剖面与WOA的平均绝对偏差都在10%以内, 因此可以使用10%作为阈值标准剔除一些明显与WOA剖面形态不相同的剖面, 从而进一步提高M值的计算(即M₂)。此外, 图 4显示了溶解氧高偏差的剖面并不明显对应于温度偏差较高或盐度偏差较高的剖面, 说明无法使用温度或盐度剖面作为溶解氧剖面的剔除标准。

增加了异常剖面的剔除标准后, 校正系数M₂的准确性得到了进一步的提高。从图 2c中可以看出, 浮标2901550后期观测的(偏离前期稳定的m_i值范围的)大部分剖面都与WOA剖面存在较大的差别, 在M₂的计算中对这些m_i值予以有效的剔除, 计算得到的M₂值为1.0361, 已与前期稳定的m_i值一致; 与其对应的是, 浮标2900765(图 2d)显示出异常剖面剔除标准对于其M值计算的改善较小, 说明该浮标后期m_i值的持续增高是由于传感器的整体退化导致的, 因此其剖面形态与WOA剖面仍存在很高的相似性(只是传感器的灵敏度不断下降), 同样的情况也发生在浮标3900347(图 2b)中。通过异常剖面的剔除标准, 我们也可以识别出哪些异常的m_i值是由于观测异常导致的, 而哪些是由于传感器退化导致的。

从图 3中也可以看出, 使用M₂值校正后的全部溶解氧剖面与WOA气候态剖面更为接近, 大于10%平均绝对偏差的剖面数量进一步减少; 图 5a则显示出, 通过增加异常剖面的剔除标准后, 浮标所有m_i值的绝对中值偏差显著下降, 说明该标准剔除掉的大部分是一些离散较大的m_i值(偏高或偏低, 如图 2a和图 2c显示)。此外, 图 5b绘制了同一台浮标计算的M₁和M₂的散点图, 可以看出, 对于绝大多数浮标, 异常剖面的剔除标准并没有改变最终校正系数的计算, 再次说明中值法对于异常值的稳健性; 不过对于某些浮标, M₂值比M₁值更接近于1, 如前文所述, 理论上漂移校正的M值应大于1, 虽然存在一些M < 1的传感器, 但其数量较少。因此, M₂值的分布更符合预期。

2.3 异常m_i值剔除标准的改进

在剔除异常剖面的改进之后, 本文提出第三步改进方案, 基于三倍绝对中值偏差的标准, 进一步剔除掉一些离散性较大的m_i值。在图 2的四台浮标中, 仅浮标2900765(图 2d)表现为较大的改进, 计算的M₃值与该浮标前期观测的稳定的m_i值基本一致。图 5c和图 5d显示, 通过增加异常m_i值的剔除标准后, 绝对中值偏差显著下降(MAD₃ < MAD₂), 但对于校正系数M的计算影响不大, 绝大多数M₃值与M₂值基本一致, 仅对于个别浮标(如图 2d)产生作用, 这说明: 1)在使用中值法以及增加异常剖面的剔除标准之后, M₂值已比较稳定, 受异常值的影响非常小, 绝大多数浮标观测的m_i值变化不大; 2)个别浮标仍然存在M₂值偏离大部分m_i值的现象, 而M₃则可以进一步避免异常m_i值的影响(如图 2d)。从校正后的溶解氧饱和度剖面与WOA气候态剖面的平均绝对偏差(MAE)的折线统计图(图 3)中也可以看出, M₃校正后的结果与M₂基本没有差别。

3 结论

从图 2四台浮标的示例以及图 3的折线统计图中都能明显地看到, 校正系数M₃比传统方法计算的M₀有了较大程度的改进。M₃的值域范围更合理, 基本不会受到异常剖面与异常m_i值的影响, 显著提高了Argo溶解氧数据库的整体数据质量。不仅如此, 上述三种改进方法对于空气测量校正法同样适用。特别是“异常剖面的剔除标准”, 可以使得计算中值时不受这些异常剖面的影响。

需要说明的是, 气候态溶解氧数据虽然可以作为一个重要的参考值, 但由于其代表的是一个“气候平均态”分布, 而现场观测数据则代表了溶解氧的“瞬时”特征, 包含了年际变化与特殊事件的影响, 两者不可能完全一致。因此, “异常剖面的剔除标准”可能存在“误判”——即将一些正常偏离平均态的溶解氧观测剖面定为“异常”。这种“误判”在计算校正系数M时, 对最终的结果影响不大, 但不能将所有高偏差的观测剖面都判定为异常而剔除。同样, “异常m_i值的剔除标准”也仅用于计算校正系数, 一些异常高的m_i值主要是受到传感器退化的影响, 其剖面形态与气候态数据吻合程度仍然很高(如图 2b和图 2d), 说明传感器退化产生的漂移同样是一种“斜率漂移”。

因此, 本研究虽然对校正系数m值的计算提出了一种改进方法, 但不能用于对浮标数据进行分级(flagging)。一些高偏差的“异常”剖面可能是由于一些瞬时过程(如平流过程、中尺度涡、台风等)导致的; 而一些异常m_i值则是由于传感器退化造成的, 如果采用较高的斜率系数(m_i), 这些高漂移数据在一定程度上仍然可以得到校正与补救。Argo浮标观测溶解氧数据的质量控制方法仍未完善, 在本研究的基础上, 未来还需进一步研究如何进行合理的数据分级, 剔除或采用特殊方法校正一些质量低的数据点或剖面, 从而提高数据库的整体数据质量。