浮游植物以其细胞形态、生物地理化学功能和元素需求等一系列特征在调节海洋生物地球化学中起着关键作用^[1]。通过对浮游植物进行探测, 可以刻画出海洋水体浮游植物种群结构及分布特征, 这将为深入理解海洋中浮游植物功能多样性提供重要支撑, 为揭示海洋生物地球化学过程提供科学依据。中国近海浮游植物种群结构复杂多样, 主要的优势藻种为硅藻(Diatoms)、甲藻(Dinoflagellate)、绿藻(Chlorophytes)和定鞭金藻(Prymnesiophytes)等^[2-6]。不同浮游植物种群呈现出不同的功能性差异, 研究浮游植物种群生物量, 掌握浮游植物种群结构特征, 这对研究海洋生态动力学和碳、硅等多种元素循环等均具有重要意义^[7-10]。

针对中国近海浮游植物种群生物量, 很多学者进行过一定研究。王俊等^[11]通过分析1992年至1993年渤海生态基础调查资料, 发现渤海共有浮游植物31属70种, 其中硅藻门的种类和数量占绝对优势; 孙军等^[7]基于渤海浮游植物群落的分析和模拟资料, 研究发现硅藻和甲藻的占比最高; 赵越等^[12]利用流式细胞仪和形态学观察等方法对2015年黄海及东海的浮游植物样品进行处理分析, 发现硅藻是黄海海域的特征藻种; 陈楠生等^[13]研究分析了1957—2019年间渤海46次航次调查数据, 结果显示渤海浮游植物的生物多样性很高, 在已鉴定的140个赤潮物种中, 硅藻类占一半以上。现有研究表明, 硅藻的生物量在中国近海浮游植物种群生物量中拥有较大的占比。

硅藻也是海洋食物链的重要组成部分。在海洋富营养化的情况下, 硅藻生物量会剧增, 从而极易引起赤潮, 对海洋生态环境造成严重破坏。同时, 硅藻还会产生有害的次生代谢物, 使得有毒物质在海洋中富集, 最终通过食物链对捕食动物的发育甚至是生命造成巨大威胁^[14-15]。由于其独特的主导地位与功能优势, 硅藻已成为多个种群中的研究重点, 对硅藻生物量的反演研究成为当务之急。

目前研究硅藻生物量的方法大多基于实测采集数据。例如, 利用实验室光学显微镜、流式细胞术等获取采集水体样本中的硅藻丰富度信息; 利用高效液相色谱法(High Performance Liquid Chro­matography, HPLC)求得现场采集水体样本的色素浓度信息, 并通过CHEMTAX软件或诊断色素分析法获得硅藻浓度信息^[16-21]。这些测量方法获取的硅藻浓度较精确, 但实测数据获取难度大, 且无法长时序、大范围地对硅藻进行观测研究。随着卫星遥感的快速发展, 遥感技术在探测宏观大范围和长时序变化的目标对象上, 具有得天独厚的优势。Xi等^[22]基于诊断色素分析法获得硅藻浓度信息, 利用奇异值分解方法(Singular Value Decomposition, SVD)对星地匹配遥感反射率(Remote sensing reflectance, R_rs)数据进行处理, 探究硅藻浓度信息与R_rs之间的潜在关联, 获取全球硅藻浓度时空分布特征。现有反演方法在一定程度上可以获得硅藻浓度信息, 但主要是针对大洋水体的硅藻浓度研究。中国近海水体成分复杂, 宏观大洋水体尺度上的硅藻浓度反演模型未必适用于区域性的近海水体, 而目前专门针对中国近海水体硅藻浓度的遥感反演研究仍属空白, 估算模型的缺失阻碍了中国近海硅藻生物量的卫星反演, 更无法实现周期性对硅藻生物量的时空监测。

为此, 针对渤海和北黄海水体, 本文基于2016年6月黄渤海和2018年7月黄渤海航次实测数据, 首先利用CHEMTAX软件, 将色素浓度数据转化为对应的藻种浓度数据, 提取硅藻浓度信息; 之后, 利用SVD, 提取标准化后的反射率光谱特征, 将SVD后得到协方差矩阵, 探究硅藻浓度与矩阵之间的潜在联系, 构建硅藻浓度反演模型; 将模型应用于2020年的6月份GOCI(Geostationary Ocean Color Imager)数据, 获取硅藻浓度月平均分布产品; 并与已有硅藻浓度反演模型进行定性地探讨比较。

1 数据与方法 1.1 研究区域概况

本文研究区域是渤海和北黄海。渤海位于37°N~ 41°N、117°E~121°E之间, 是我国最北部的半封闭内海, 仅在东部以渤海海峡与北黄海相通。受黄河、辽河与海河等河流影响, 大量泥沙冲积推挤海床, 造成海底平坦。渤海平均深度只有18 m, 丰富的营养物质经过河流进入渤海, 为浮游植物提供良好的生存条件^[23]。北黄海是一个半封闭海域, 被山东半岛、辽东半岛和朝鲜半岛围绕, 受黄海暖流、季节性黄海冷水团和河流淡水的共同作用影响, 其生态环境复杂多样^[24]。

1.2 数据获取与处理

本研究所采用的是2016年6月黄渤海和2018年7月黄渤海航次的实测数据, 主要包括遥感反射率R_rs数据以及基于HPLC方法获取的色素浓度数据, 通过质量控制后参与模型构建的实测航次数据具体信息与采样的航次站点分布分别如表 1与图 1所示。

1.2.1 遥感反射率数据获取与处理

现场测量R_rs数据是由美国ASD公司生产的FieldSpec4Hi-ResNG地物光谱仪测量获得, 它可以测量的光谱波长介于350~2 500 nm之间, 其中在700 nm以内的光谱分辨率是3 nm。在进行作业的过程中, 选择天气晴朗的环境, 使用ASD测量目标地物完毕后使用ASD自带的RS3软件对采集到的数据进行处理。本航次所采用的光谱范围在400~700 nm之间, 对数据进行重采样的采样间隔为1 nm。为保证模型的准确性, 剔除实测R_rs数据中小于零异常值, 经过处理后的R_rs光谱如图 2所示, 其中不同颜色的光谱曲线代表不同站点的R_rs光谱。

1.2.2 浮游植物种群浓度数据获取与处理

参照Aiken的质量控制方法^[25]: 1) 剔除总叶绿素a浓度和辅助色素浓度之差大于两者之和的30%; 2) 剔除12种辅助色素中浓度小于0.001 mg·m^–3的数目小于6个的样本^[26]。在航次所采集色素浓度基础上, 利用CHEMTAX (CHEMical TAXonomy, CHEMTAX)软件获取浮游植物种群浓度。CHMTAX软件可通过由每种浮游植物不同色素的比例组成的初始比率矩阵, 将色素信息转化为对应的藻种浓度信息。本研究的CHEMTAX软件计算所得藻种及其藻种包含色素信息是参考Aiken等^[25]和Zhang等^[26]的研究, CHEMTAX软件计算是使用Pan等^[27]和Sun等^[28]研究中的初始色素比率矩阵。在CHEMTAX软件计算时, 本研究定义64个矩阵, 每个矩阵各藻种和对应藻种色素之间的比率均是初始比率矩阵乘以0.65~1.35之间的随机数。在获取的64组种群浓度中, 选取6组种群浓度均方根误差(R_MSE)最低的平均值作为最终藻种浓度, 并且为保证后续模型的准确性, 只保留CHEMTAX软件计算所得大于0.005 mg·m^–3的藻种浓度。

1.2.3 GOCI卫星数据获取与处理

本研究使用的卫星数据为GOCI数据, GOCI是韩国航天局于2009年发射的COMS(Communication Ocean and Meteorological Satellite)卫星上搭载的水色遥感器, 它可以覆盖我国渤海、黄海以及东海部分区域, 空间分辨率约为500 m, 时间分辨率可达1 h, 每天可连续进行8次观测, 并设有6个可见光波段(412 nm、443 nm、490 nm、555 nm、660 nm和680 nm)。

为研究渤海和北黄海硅藻的长时序变化特征, 本研究收集的卫星数据均来自于韩国海洋卫星中心(KOSC, Korea Ocean Satellite Center, http://kosc.kiost.ac.kr/eng/p10/kosc_p11.html)。从KOSC下载2020年6月1日至2020年6月30日Level-1B级日产品, 并使用GOCI官方自带数据处理软件GDPS(GOCI Data Processing Solutions)进行预处理, 获得Level-2P级$ R_{rs}^{sat} $日产品, 并计算得到$ R_{rs}^{sat} $月平均产品。

1.3 精度指标

本研究通过决定系数(R²)、均方根误差(R_MSE)、平均绝对百分比误差(M_APE)和中值误差(M_E)来对建模效果进行精度校验, 相应的计算公式如下所示:

$ {R^2} = \frac{{\sum\limits_{i = 1}^N {({y_{{\text{o}}i}} - {{\bar y}_{\text{o}}})({y_{{\text{p}}i}} - {{\bar y}_{\text{p}}})} }}{{\sqrt {\sum\limits_{i = 1}^N {({y_{{\text{o}}i}} - {{\bar y}_{\text{o}}})} \sum\limits_{i = 1}^N {({y_{{\text{p}}i}} - {{\bar y}_{\text{p}}})} } }} , $

(1)

$ {R_\rm{MSE}} = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{({y_{{\text{o}}i}} - {y_{{\text{p}}i}})}^2}} } , $

(2)

$ {M_\rm{APE}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left| {\frac{{{y_{{\text{o}}i}} - {y_{{\text{p}}i}}}}{{{y_{{\text{o}}i}}}}} \right| \times 100\% } \text{，} $

(3)

$ {M_\rm{E}} = Median\left| {\frac{{{y_{{\text{o}}i}} - {y_{{\text{p}}i}}}}{{{y_{{\text{o}}i}}}}} \right| \times 100\% , $

(4)

式中，y_oi为第i个实测硅藻浓度, y_pi为第i个预测硅藻浓度, $ {\bar y_{\text{o}}} $实测硅藻浓度平均值, $ {\bar y_{\text{p}}} $为预测硅藻浓度平均值, N为样品数, 并利用留一交叉验证法验证模型精度。

2 硅藻浓度反演模型开发

利用SVD, 可以提取浮游植物R_rs特征, 进而获取浮游植物色素或种群生物量信息, 这种方法已经被初步应用于色素浓度反演模型的构建^[29-30]。本研究基于SVD, 提取标准化后的反射率光谱特征, 利用SVD后的协方差矩阵与CHEMTAX软件计算的硅藻浓度进行线性拟合, 建模流程如图 3所示。

首先将实测R_rs数据作标准化处理: 用实测R_rs数据减去平均实测R_rs值, 之后除以实测R_rs的标准差, 提取标准化后的R_rs数据中与GOCI相同的6个波段数据, 记为矩阵X。矩阵X根据公式(5)进行SVD后得到3个矩阵: PCs(左奇异值向量矩阵U)、EOFs(右奇异值向量矩阵V)、奇异值矩阵Σ。

$ \boldsymbol{X}=\boldsymbol{U} \sum \boldsymbol{V}^{\mathrm{T}} \text {, } $

(5)

式中, V为N行N列矩阵; U为M行N列矩阵; Σ为N行N列矩阵, 是矩阵X在对角线上的奇异值, 并按递减顺序排列, 其中M为样品数, N为波段数。将CHEMTAX计算的硅藻浓度取对数, 结合公式(6)与协方差矩阵U进行线性拟合。本研究使用留一交叉验证法(Leave-One-Out Cross Validation)进行建模, 即每次使用N–1个样品数建模, 留下1个作为验证点, 循环N次建模, 最后取模型的平均作为最终模型。

$ \ln \left(y_{\mathrm{dia}}\right)=a_0+a_1 u_1+a_1 u_1+\cdots+a_n u_n, $

(6)

式中, y_dia为硅藻浓度, a₀为截距, a₁、a₂、…、a_n为回归系数, u₁、u₂、…、u_n为矩阵U中列向量的前n个元素。这里引入赤池信息量准则(AIC)来决定n的取值, 以此来寻找最优拟合^[29-30]。

将GOCI的$ R_{{\text{rs}}}^{{\text{sat}}} $数据标准化处理, 得到矩阵X^sat, 根据公式(7)将X^sat投影到V上来获取U^sat, 利用公式(8)进行线性合成, 即可得到硅藻反演浓度$ y_{{\text{dia}}}^{{\text{sat}}} $, 反演流程如图 3所示。

$ \boldsymbol{U}^{\text {sat }}=\boldsymbol{X}^{\text {sat }} \times \boldsymbol{V} \times \Sigma^{-1} , $

(7)

式中, X^sat为GOCI卫星标准化后$ R_{{\text{rs}}}^{{\text{sat}}} $数据, V为实测数据SVD后矩阵V^T的转置矩阵, Σ为实测数据SVD后奇异值矩阵的逆矩阵。

$ \ln (y_{{\text{dia}}}^{{\text{sat}}}) = {a_0} + {a_1}u_1^{{\text{sat}}} + {a_2}u_2^{{\text{sat}}} + \cdots + {a_n}u_n^{{\text{sat}}} , $

(8)

式中, a₀为模型截距, a₁、a₂、…、a_n为模型回归系数, $ u_1^{{\text{sat}}} $、$ u_2^{{\text{sat}}} $、…、$ u_n^{{\text{sat}}} $为矩阵U^sat中列向量的前n个元素, n与建模时一致。

3 结果与分析 3.1 硅藻浓度反演模型验证

基于硅藻浓度反演模型, 反演得到对应站点处的硅藻浓度信息。图 4(a)是实测硅藻浓度与反演硅藻浓度之间的散点对比图, 中间的红色虚线为1︰1线。从图中可以看出, 实测硅藻浓度和反演硅藻浓度相关性较高(决定系数R²=0.80, P < 0.001), 散点均匀地分散在1︰1线两侧, M_APE为58.62%, M_E为39.12%, R_MSE为0.19 mg·m^–3。基于实测R_rs数据, 通过与GOCI匹配获得卫星$ R_{{\rm{rs}}}^{{\rm{sat\_mp}}} $数据(时间窗口为±3 h、空间窗口为±0.01°×0.01°), 得到实测硅藻浓度与匹配卫星$ R_{{\rm{rs}}}^{{\rm{sat\_mp}}} $数据反演得到的硅藻浓度之间的散点对比图。如图 4(b)所示, 其中R²为0.78(P < 0.001)、M_APE为44.43%, M_E为35.55%, R_MSE为0.11 mg·m^–3。匹配的卫星$ R_{{\rm{rs}}}^{{\rm{sat\_mp}}} $数据反演硅藻浓度模型反演值与实测硅藻浓度之间拟合程度较高, 这表明本研究构建的硅藻浓度反演模型具有良好的适用性。

3.2 硅藻浓度反演模型示范应用

考虑到建模所用的实测航次为夏季, 因此将本文建立的硅藻浓度反演模型应用于2020年6月份GOCI卫星Level-2级$ R_{{\rm{rs}}}^{{\rm{sat}}} $月平均产品, 获取2020年6月渤海和北黄海的硅藻浓度月平均分布图(图中多处空白区域是由于GOCI中R_rs产品缺测造成), 如图 5a所示。整体上, 渤海硅藻浓度明显高于北黄海。对于渤海海域, 辽东湾、渤海湾、莱州湾, 硅藻浓度整体保持在同一水平, 约0.4~1.0 mg·m^–3。硅藻浓度高值主要位于辽河和海河等河流入海口及沿岸海域, 约0.9 mg·m^–3以上; 随着离岸距离的增加, 硅藻浓度呈递减趋势, 在海盆中心部分水域硅藻浓度达到最低值, 约0.3~0.4 mg·m^–3。对于北黄海海域, 近岸水域的硅藻浓度较高, 北黄海中心绝大部分海域, 硅藻浓度相对较低, 约0.1 mg·m^–3。基于2020年6月GOCI $ R_{{\rm{rs}}}^{{\rm{sat}}} $月平均数据, 通过YOC(The Yellow Sea Large Marine Ecosystem Ocean Color Work Group, YOC) 模型^[31]获得相应的叶绿素a渤海与北黄海浓度空间分布, 进一步获取硅藻浓度与叶绿素a浓度占比(如图 5b所示)。可以发现, 硅藻浓度与叶绿素a浓度百分比值在近岸较高, 可达30%, 随着离岸距离增加而降低。渤海远岸的百分比值高于北黄海远岸, 分别在12%与5%左右。

4 讨论 4.1 硅藻浓度空间分布的潜在因素分析

就上述反演得到的2020年6月渤海和北黄海水体硅藻浓度的空间分布而言, 整体上, 硅藻浓度的分布呈现出近岸浓度高、离岸浓度低的趋势特点, 与前人研究一致^[32-34]。海水温度、盐度以及营养盐等是影响浮游植物空间分布的重要因素。其中, 营养盐的作用更显著^[35-36]。在众多浮游植物种群中, 硅藻偏爱在低温、高营养盐的水环境中生存^[37], 而近岸海域水深较浅, 携带大量营养盐的地表径流与海水能够充分混合, 使得河口区营养物质较为丰富, 并且水流逐渐变缓, 水体浊度降低, 有利于硅藻等浮游植物的生长繁殖^[38], 因此可能造成硅藻浓度较高。随着离岸距离的增加, 海水的深度增大, 营养物质的丰度不及近岸水域, 所以硅藻浓度也逐渐降低。

渤海的硅藻浓度整体上高于北黄海, 这可能是由于渤海整体上是一个近封闭内海, 受海河、黄河和辽河等陆源输入的影响较大, 河流冲淡水带来了丰富的营养物质与海水进行交汇, 为硅藻的生长提供了良好条件^{[11, 35, 39]}。相较于渤海, 北黄海受陆源输入的影响比较小, 同时又具备低盐度、高营养的渤海沿岸水与高盐度、低营养的黄海水, 所以可能导致硅藻浓度相对于渤海较低^[24]。

4.2 硅藻浓度反演模型的合理性和局限性分析

近年来, 许多学者对浮游植物种群生物量有专门的探究^{[9, 40-42]}, 但目前准确计算浮游植物各种群生物量仍然是个严峻的挑战。在水色遥感中, 太阳辐射光进入水中经过浮游植物等水体组分散射后, 其中一部分离开水面被搭载在卫星的传感器所接收, 所以通过R_rs可以获取浮游植物信息^[43-45]。由于多维度数据集难以提取其特征信息, 而SVD方法可以通过降低数据维度来寻找描述数据集内的主导特征信息^[46]。因此本研究通过SVD方法来降低R_rs时间和空间维度, 获取能够描述R_rs光谱特征的主导信号^{[29, 47]}。

以往硅藻浓度反演模型研究中, Xi等^{[22, 48]}基于共享航次HPLC色素数据和卫星$ R_{rs}^{sat} $数据, 通过SVD方法构建全球硅藻浓度反演模型, 而后又加入海表温度(SST)优化模型, 得到精度较高的硅藻浓度反演模型。本研究利用Xi等模型用于区域的中国近海水体, 反演硅藻浓度的模型表现并未尽如人意(模型用于本研究区域的精度效果: R²及M_E分别为0.51与70.38%)。分析其原因可能是由于模型匹配的卫星$ R_{{\rm{rs}}}^{{\rm{sat\_mp}}} $数据, 受到时间和空间窗口、大气校正等过程的误差引入, 导致其与浮游植物信息之间关系的精度受到限制。同时中国近海水体为典型二类水体, 其含量和光学特性复杂, 浑浊近海浮游植物群落结构可能与大洋开阔一类水体浮游植物群落结构有所差异^[49-50], 导致全球宏观尺度的硅藻浓度反演模型不适用于中国近海。

叶绿素a可以描述浮游植物生物量变化^[51], 因此通过实测总叶绿素浓度与硅藻浓度的数据, Sun等^[28]建立了二者之间的经验关系模型, 虽然在其研究中表现出较好的模拟精度(R²与M_APE分别为0.808与54.2%), 但是经验模型归根结底对于时空的依赖性较强, 并且浮游植物种群结构成分变化复杂^[52-53], 因此经验模型是否适合于其他区域和时间仍需进行进一步比较验证。此外, 根据浮游植物吸收可加性理论, 总浮游植物吸收可分解成不同藻种吸收^[54], 基于此, 樊杰等^[55]寻找硅藻吸收光谱不同波段处的高斯峰与实测硅藻浓度之间的关系来构建硅藻浓度反演模型并且具有较高的精度。基于本研究数据, 利用樊杰等^[55]的硅藻浓度反演模型, 反演得到相应的硅藻浓度, 并与本研究所反演硅藻浓度进行比较, 如图 6所示。可以发现本研究反演的硅藻浓度与实测的硅藻浓度之间的R²(图 6a)虽然低于与基于樊杰等^[55]的硅藻反演模型(图 6b), 分别为0.80与0.83, 但在M_APE、与M_E均好于基于樊杰等^[55]硅藻反演模型, 分别为58.62%与61.15%, 39.12%与52.25。由于樊杰等^[55]是针对黄渤东海的硅藻浓度所构建的硅藻反演模型, 该模型可能具有区域性, 而本研究针对渤海和北黄海的硅藻浓度构建硅藻浓度反演模型, 因此可能造成本研究模型精度较高于樊杰等^[55]模型。再者, 卫星浮游植物吸收系数产品是由R_rs产品转化而来, 因此在经过多个经验关系的步骤过程中可能会引入累积误差, 樊杰等^[55]模型精度可能会受到限制。

相比上述研究, 本研究构建硅藻浓度反演模型中HPLC色素浓度数据和R_rs数据是于中国近海站点同步测量获得, 在空间上和时间上是同步的, 以便准确获取近海浮游植物信息, 保证建模精度。由于模型是利用实测R_rs数据和实测硅藻浓度信息建立, 避免多个经验关系带来的累积误差。此外通过GOCI卫星匹配$ R_{{\rm{rs}}}^{{\rm{sat\_mp}}} $数据并应用于构建的模型, 得到的反演硅藻浓度与实测值比较, 在一定程度上表明该模型适用于GOCI卫星$ R_{{\rm{rs}}}^{{\rm{sat}}} $产品。

然而, 本研究所开发模型的表现仍受到一些局限。本研究基于2个夏季航次实测R_rs数据和HPLC色素浓度数据中, 由于R_rs随季节变化而变化, 构建的模型是否适用于4个季节还未知, 后续在模型构建中将考虑加入更多不同季节航次数据。模型构建中所需的硅藻浓度是基于CHEMTAX软件转化HPLC色素浓度获取, 尽管CHEMTAX软件被普遍认为是一种计算浮游植物群落生物量的有效方法^[56-60], 但由于其初始色素比率矩阵受时间和空间影响较大, 需要在未来工作中进一步研究获取更准确的最优初始比率矩阵, 提高模型反演精度。

5 结论

针对中国近海硅藻浓度卫星估算的模型方法仍较为欠缺的问题, 本研究基于航次实测数据构建了硅藻浓度的遥感反演模型, 模型展现出较好的精度效果(R² = 0.80, P < 0.001; M_APE = 58.62%, M_E = 39.12%), 实测硅藻浓度与模型反演值之间表现出较好的一致性。同时将模型应用于GOCI 2020年6月份的$ R_{rs}^{sat} $月平均产品, 由于近海硅藻可能受水温度、盐度以及营养盐等影响, 其浓度呈现近岸高, 离岸低的空间分布特征。总体而言, 基于SVD技术所构建的硅藻浓度反演模型具有较好的卫星应用潜力, 一定程度上弥补了中国近海水体硅藻浓度卫星遥感估算能力上的不足。后续可以尝试将该模型方法扩展到其他浮游植物优势藻种之上, 为海洋水体浮游植物群落结构的卫星遥感研究提供方法论基础, 服务于海洋水体生态环境的精细化监测。