基于GOCI 数据的黄河口高悬浮物浓度海域的叶绿素a 提取指数构建与浓度反演

张曦元，万剑华，刘善伟，宋冬梅

（中国石油大学（华东）海洋与空间信息学院，山东 青岛 266580）

海水的固有光学特征是水体中所含各类物质对光的吸收和散射的结果，光谱特性复杂且具有区域性差异[1]。GORDON H R 等[2]依据水质成分将海水分为Ⅰ类水体和Ⅱ类水体，Ⅰ类水体的光学特性主要受叶绿素a（Chlorophyll-a，Chl-a） 及其分解产物的影响；Ⅱ类水体的光谱特征受Chl-a、总悬浮物（Total Suspended Matter，TSM）及黄色物质等多重影响，其光学特性也更为复杂。渤海西南部黄河口附近及莱州湾、渤海湾近岸海域为典型的Ⅱ类水体，黄河携带大量泥沙汇入莱州湾，在莱州湾、渤海湾海域随洋流扩散，高浓度TSM 的不透光性干扰了Chl-a 浓度的光谱量测，使得该海域水体光谱特征复杂，难以准确估算Chl-a 浓度[3]。采用GOCI（Geostationary Ocean Color Imager） 卫星遥感技术构建Chl-a 浓度反演模型，能够动态监测高TSM 浓度海域的Chl-a 浓度，打破传统方法在时空上的局限。

目前，国内外学者针对复杂浑浊水体Chl-a 浓度的反演，提出了多种方法来构建其反演模型。三波段因子（Three-band factor，TBA）模型可用来估算浑浊水体的Chl-a 浓度，该方法具有明确的物理意义且精度较高[4]，因而被研究者广泛应用[5-7]；四波段半解析模型能够提高太湖高浊度水体Chl-a 浓度的估算精度[8]；改进三波段半分析模型可用于黄河口高浑浊水体Chl-a 浓度的反演，该算法需满足λ1和λ2的波长在660～690 nm 之间，且λ3的波长要大于750 nm 的条件[9]。修正型叶绿素指数是针对实测水体光谱的近红外与红光波段的反射率进行比值修正，能够削减悬浮物对叶绿素光谱信息的影响，有效提高了太湖悬浮物主导型水体的Chl-a 浓度的反演精度[10]；综合叶绿素指数能够提高长江口浑浊水域Chl-a 浓度的估算精度，但其只能对MERIS（Medium Resolution Imaging Spectrometer Instrument）等高光谱数据使用[11]。上述方法针对高浓度悬浮物水体Chl-a 浓度反演模型的研究取得了较高的精度，但其普适性较差，对卫星数据波段区间的要求较高，对于时间分辨率较高的GOCI 数据等并不成立。

因此，本文针对黄河口近岸、莱州湾及渤海湾高TSM 浓度的浑浊海域，分析了实测水体光谱曲线的特征，提出一种新的Chl-a 提取指数（Chlorophyll Extraction Index，CEI），该指数旨在高TSM浓度的水体中，通过因子组合来削弱TSM 对Schl的影响。

1 研究区与数据

1.1 研究区

黄河是世界上泥沙含量最高的河流，每年向海洋输入巨量的悬浮泥沙以及各类营养盐等物质，因黄河下游流速的降低易造成入海口物质沉积[12]。随着经济的快速发展和排水区人为活动的强烈影响，近二十年来，污水排放的增加、农业化肥的过量使用以及工业废水的排放，使得黄河口近岸以及莱州湾、渤海湾海域的海水富营养化十分严重，藻类、浮游生物的大量繁殖导致海洋生态环境严重恶化[13]。因此，本文以高TSM 浓度的黄河口近岸及莱州湾、渤海湾海域为研究区，利用遥感技术监测该区域Chl-a 浓度的变化有利于保护渤海的生态环境。

1.2 实测数据

本文实测样本数据包括光谱反射率数据和水质参数浓度数据，其观测站点位置如图1 所示。水体光谱反射率数据是2011 年7 月4 日于黄河口近岸测得的，共18 组；水质参数浓度数据包括Chl-a 和TSM 浓度，2011—2018 年共收集有效样本数据619组，采集的样本数据如表1 所示。取样时间为当地时间9∶00 到15∶00，天气晴朗、海上无大风。

表1 取样数据说明

由实测数据可知，黄河口近岸海域叶绿素浓度范围为0.73～15.4 μg/L，均值为3.79 μg/L；悬浮物浓度范围为225.36～2 291.84 mg/L，均值为536.78 mg/L。渤海湾南部海域叶绿素浓度范围为0.51～18.1 μg/L，均值为3.93μg/L；悬浮物浓度范围为3.08～1 086.45mg/L，均值为116.41 mg/L。莱州湾海域叶绿素浓度范围为0.45～40.7 μg/L，均值为4.46 μg/L；悬浮物浓度范围为4.17～1 760.86 mg/L，均值为141.92 mg/L。由此可见，黄河口近岸悬浮物浓度远远高于其他海域，与顺布日等[14]研究结果较为一致，严重干扰了Chl-a浓度对光谱波段敏感性的响应。

1.3 卫星数据及预处理

GOCI 是世界上首颗专门监测海洋水色产品的静止轨道卫星[15]，其参数信息如表2 所示。GOCI卫星监测范围能够完全覆盖我国渤海海域，每日可获取8 景影像，在动态变化监测方面具有较大优势。

表2 GOCI 卫星波段参数信息表

为获取与现场取样时间匹配的影像数据，本文根据具体取样时刻选择同步的GOCI 影像。影像预处理包括大气校正、几何校正、云掩膜等操作，其中大气校正精度是影响反演模型质量的主要因素。由于渤海近岸水体的浑浊度高，其复杂的光学特性严重影响到大气校正的精度。GDPS 软件内置的算法难以降低渤海浑浊水体对大气校正带来的干扰[15]，何贤强基于紫外波长提出了一种适用于海洋水色图像的大气校正算法（UV-AC），能够解决光谱匹配算法在浑浊水域失效的问题[16]。因此，本实验采用UV-AC 法对GOCI 数据进行大气校正。

2 研究方法

本文针对黄河口及其附近海域高TSM 浓度的浑浊水体，通过分析水体实测光谱数据变换曲线与Chl-a 浓度之间的相关性，确定GOCI 卫星数据所对应的敏感波段的反射率值。根据最大叶绿素指数（Maximum Chlorophyll Index，MCI） 的波段基线内差值原理构建叶绿素敏感指数（Chlorophyll-a Sensitivity Index，Schl）、悬浮物敏感指数（Total Suspended Matter Sensitivity Index，STSM），进一步构建叶绿素提取指数，该方法能够削弱TSM 对Schl的影响。根据采样点同步的GOCI 卫星数据得到CEI，以CEI 和实测Chl-a 浓度为样本集，基于随机森林算法，采用训练样本构建Chl-a 浓度的反演模型，而后利用测试样本进行模型精度检验。

2.1 水体光谱数据的特征分析

水体的光谱反射特征受到纯水、Chl-a、TSM 等有色物质的共同影响。本研究采集的海水光谱信息位于黄河口近岸，大量泥沙的汇入使得TSM 浓度相对较高；内陆污水的排放加剧了海水的富营养化程度，浮游植物的大量繁衍使得Chl-a 浓度相对较高；此外，仪器、光照等外界因素也会对数据采集产生一定的影响，因此该区域水体的光谱特征较为复杂。现场测量的海表光谱曲线如图1 所示。

图1 黄河口近岸海域实测光谱反射率曲线

黄河口近岸海域水体的光谱曲线基本符合Ⅱ类水体的特征，但又有其独特性。400～500 nm 区间由于TSM 和Chl-a 的共同强吸收作用，光谱值较低；550 nm 附近由于TSM 散射作用形成较大的反射峰；660 nm 附近出现弱反射峰，此处受TSM 吸收作用和Chl-a 反射作用的共同影响；680 nm 附近出现吸收谷，该谷值通常被认为是判定水体是否含有Chl-a的依据；750～850 nm 范围内因TSM 对水体光谱的后向散射作用而形成微小的反射峰，当水中TSM含量增大时，其反射峰也会向长波方向移动，也称之为红移现象；900 nm 后的光谱反射率存在较大噪声，对此不做具体分析。

为了减少外界因素对采集光谱数据所造成的影响，本文采用二阶微分变换的方式对原始反射率进行处理。原始光谱反射率及二阶微分处理后的光谱反射率与Chl-a 和TSM 浓度的相关关系如图2 所示。图2（a）为原始光谱反射率与Chl-a 和TSM 浓度的相关性分析，图2（b）为二阶微分处理后光谱反射率与Chl-a 和TSM 浓度的相关性分析。

由图2（a）可见，原始光谱的反射率与Chl-a 和TSM 浓度最大相关系数分别为-0.32 和0.98，分别出现在737 nm 和388 nm 附近。二阶微分处理后的光谱反射率与Chl-a 和TSM 浓度的相关关系如图2（b）所示，经二阶微分处理后的光谱反射率与Chl-a 浓度的相关性较原始光谱反射率有明显提升，但其与TSM 浓度的相关性略有下降，Chl-a 和TSM 浓度最大相关系数分别为0.88 和-0.85，分别出现在745 nm 和550 nm 附近，该相关系数能够满足实验要求。

图2 光谱反射率与Chl-a 和TSM 浓度的相关性分析

此外，经二阶微分处理后的光谱反射率与Chl-a和TSM 浓度的相关性分析曲线大致呈对称分布，555 nm、660 nm、680 nm 及745 nm 处呈明显对称分布。Chl-a 浓度与二阶微分光谱反射率的相关性在这四处均呈现正相关性，555 nm 附近光谱与Chl-a浓度由负相关迅速转为强正相关，且与TSM 浓度呈现明显的负相关；660 nm 附近光谱与Chl-a 浓度正相关系数处于较低水平，与TSM 浓度的负相关系数处于较高水平，其原因可能是此处TSM 吸收作用大于Chl-a 反射效果；680 nm 附近光谱与TSM浓度的负相关系数由-0.78 骤减为-0.25，其与Chl-a浓度的正相关系数也有所下降，可能是TSM 的强反射与Chl-a 的强吸收作用的结果；745 nm 附近光谱与Chl-a 浓度的正相关系数达到最高值，与TSM 浓度的负相关系数超过-0.6。

2.2 叶绿素a 提取指数构建

MCI 利用2 个端点波段的遥感反射率，构建一条跨2 个端点波段区间的光谱基线。水体的Chl-a在信号波段之间具有光谱反射峰，其遥感反射率值与波段基线处的内差值之差为MCI[17]。

根据MCI 的波段基线内差值的原理，首先确定Chl-a 与TSM 的光谱吸收波段及多个敏感波段。通过分析海表光谱曲线特征以及Chl-a、TSM 与二阶微分处理后的光谱反射率之间的相关性可知：555 nm处的光谱反射率与TSM 浓度由强正相关突然转为强负相关；660 nm 附近TSM 对光谱的吸收效果较Chla 更为明显；680 nm 附近Chl-a 的吸收作用较强，且该波段是判断水中Chl-a 浓度高低的重要依据；745 nm 处因TSM 对水体光谱的后向散射作用而形成微小的反射峰，与Chl-a 浓度显示强正相关。在此基础上，分别计算Chl-a 与TSM 吸收波段的光谱基线内差值，本文将该值定义为Schl和STSM，其具体表达式如下。

式中，R 表示GOCI 数据在该波段的反射率；λ表示该波段对应的中心波长。

如图3 所示，Schl是从680 nm 处的反射率值到660 nm 和745 nm 之间基线的反射率值的距离；STSM是从660 nm 处的反射率值到555 nm 和745 nm 之间基线的反射率值的距离。假设Chl-a 浓度是固定的，因为660 nm 处TSM 对光谱的吸收作用明显，660 nm 的反射率值将会随着TSM 浓度的增加而减少，则Schl将增大、STSM将减小。如果仅使用Schl对来反演Chl-a 浓度，Schl的增大将导致Chl-a 浓度的反演值偏高。

图3 基于GOCI 遥感反射率和Chl-a 浓度构建的叶绿素提取指数示意图

因此，当水体TSM 浓度较高时，Schl难免受到STSM的干扰。通过因子组合Schl-STSM来削弱TSM 对Schl的影响，能够更好地突显Chl-a 对光谱反射率特征变化的敏感性，该因子组合即为CEI。

2.3 基于随机森林算法构建Chl-a 反演模型

随机森林（Random Forests，RF） 是一种联合多个回归决策树的算法，利用bootstrap 重采样技术，对训练集采取随机有放回的重复抽取n 个样本，并在bagging 的基础上对每棵决策树进行随机特征的选择，最后，采用简单平均法对多个分类回归树得到的回归结果进行算术平均，得到最终的模型输出[17]。随机森林因其随机性地引入，使得结果不容易过拟合，具有良好的非线性数据集处理能力，可以达到反演模型的精度要求。

本实验利用2011—2018 年多次采集的黄河口近岸、莱州湾及渤海湾水域的实测Chl-a 浓度及同期的GOCI 卫星遥感影像，以Schl和Chl-a 浓度为样本集，根据空间均匀分布的原则，采用交叉验证法对不同时期所测样本进行筛选，分别得到591 个训练样本和200 个测试样本。利用591 个训练样本基于随机森林算法构建黄河口近岸高TSM 浓度海域的Chl-a 浓度反演模型，而后采用决定系数（R2）及平均相对误差（MRE）对200 个测试样本集进行模型精度检验。精度评价指标计算公式如下。

式中，yi和y 分别为Chl-a 浓度的实测值和反演值；n 为样本数量；决定系数R2在0～1 之间，值越大，模型的反演精度越高；MRE 值越小，模型的反演精度越高。

3 结果与分析

3.1 基于不同指数法构建Chl-a 反演模型

本文将619 组样本的实测Chl-a 浓度分别与三波段因子（TBA）、归一化叶绿素指数（NDCI）、Schl以及CEI 进行相关性分析，结果如图4 所示，决定系数I（R2）为Schl＞CEI ＞NDCI ＞TBA。

图4 不同指数与实测Chl-a 浓度的相关性分析

利用419 个训练样本及同期的GOCI 卫星遥感影像，基于随机森林算法分别建立TBA、NDCI、Schl及CEI 与Chl-a 浓度的反演模型，而后采用200个测试样本依次对其模型进行精度检验并对比分析，结果如表3 所示。

目前，叶绿素浓度反演精度，对于Ⅰ类水体误差不大于30%，而对于Ⅱ类水体误差40%～50%。如表3 所示，无论是训练样本还是测试样本，基于CEI 和Schl指数反演得到的Chl-a 浓度值与实测Chl-a浓度值之间的误差均在40%～50%之间，其精度评价参数均优于TBA 与NDCI。因此，认为CEI 和Schl指数能够满足Ⅱ类水体反演精度的要求，可在渤海南部海域开展Chl-a 浓度的反演研究。

3.2 针对不同类别水体构建Chl-a 反演模型

由表3 可知，基于CEI 构建的反演模型精度优于Schl，其原因可能是本实验在黄河口及其附近海域采集的样本多位于高TSM 浓度的海域中，在高TSM 浓度的海域中，采用CEI 构建的Chl-a 反演模型的精度较高，而Schl更适合TSM 浓度低的海域。因此，将研究区分为TSM 主导型水体和Chl-a 主导型水体两类，分别采用CEI 和Schl构建Chl-a 反演模型。

表3 训练样本及测试样本的评价参数结果

3.2.1 确定分区标准

GOCI 卫星每天可获取8 景影像，云雾的干扰会使得Chl-a 浓度数据出现缺失或异常，下午（14时）海雾发生频数最低[18]，因此本实验选取的GOCI影像均为同步于实测样本数据当天的14 时的影像数据。根据GOCI 卫星光谱数据采用况润元提出的光学分类法[19]，利用实测样本数据的采集时间及同步的GOCI 卫星数据，将黄河口近岸海域、莱州湾及渤海湾海域划分为Chl-a 主导型水体和TSM 主导型水体两部分。此时公式如下。

式中，R6、R4分别为GOCI 的第6 波段和第4波段的反射率值；λ6、λ4分别对应其中心波长680nm和555nm。将（R6-R4）/λ6-λ4＞0 的海域归为TSM主导型水体，反之则为Chl-a 主导型水体。

3.2.2 构建不同类型水体的Chl-a 反演模型

根据空间均匀分布的原则，采用交叉验证法对不同时期所测样本进行筛选，符合Chl-a 主导型水体的采样点共有123 个，符合TSM 主导型水体的采样点共有496 个。

基于随机森林算法构建Chl-a 浓度反演模型时，Chl-a 主导型水体采用Schl，为进一步验证CEI 在高TSM 浓度海域的精度，TSM 主导型水体采用Schl与CEI 进行对照，而后利用R2及MRE 分别对其测试样本进行模型精度检验，结果如图5、表4 所示。

图5 不同类别水体的Chl-a 反演模型精度对比

对比表3 和表4 可以看出，基于不同类别水体构建的反演模型的MRE 明显减小，R2均得到明显提高，因此，在TSM 浓度较高的海域开展Chl-a 反演模型研究时，将样本分为Chl-a 和TSM 主导型水体后，再分别采用Schl与CEI 构建Chl-a 反演模型能够提高Chl-a 的反演精度。

表4 训练样本及测试样本的评价参数结果

3.2.3 Chl-a 反演模型时空分析

本文基于随机森林训练模型对渤海海域开展了连续8 年（2011—2018 年）的Chl-a 浓度反演，实验每月选取一幅晴朗天气14 时的影像数据进行反演，按空间像素点进行平均，获得Chl-a 浓度的年平均值。年平均反演结果如图6 和图7 所示，2011—2018 年渤海海域Chl-a 浓度的年际变化无显著规律性，但其空间分布情况具有明显的区域性特征，渤海湾沿岸及黄河口近岸的莱州湾海域的浓度值最高。

图6 渤海湾海域的Chl-a 浓度反演图

图7 莱州湾海域的Chl-a 浓度反演图

通过分析渤海海域地理位置及查阅相关资料发现，渤海湾潮差小、流速低、水体交换能力较弱，内陆废水注入的营养盐有助于浮游生物的生长繁殖，因而导致Chl-a 含量增加；莱州湾Chl-a 浓度偏高的原因可能是黄河径流为近岸海域带来丰富的营养物质，近岸水深浅且光照充足，为浮游植物提供适宜的生长环境，从而导致Chl-a 浓度偏高。

4 结 论

（1）通过分析水体光谱变换曲线与Chl-a、TSM浓度的相关性，根据MCI 原理构建了Schl和STSM，进一步构建了CEI，CEI 通过Schl-STSM的差值组合大大降低了TSM 对Schl的影响。

（2）基于随机森林算法，根据GOCI 卫星的反射率与Chl-a 浓度值，分别采用TBA、NDC、Schl及CEI 构建Chl-a 浓度的反演模型，结果表明：采用TBA 和NDCI 所构建的Chl-a 反演模型的浓度值与实测值之间的MRE 均大于50%，不能满足Ⅱ类水体反演精度的要求；基于Schl构建的反演模型精度优于CEI，在高TSM 浓度的海域中，采用CEI 能够提高Chl-a 浓度反演模型的精度。

（3）根据波谱斜率将渤海南部的莱州湾、渤海湾海域划分为Chl-a 主导型水体和TSM 主导型水体两类，分别采用Schl与CEI 构建Chl-a 主导型水体和TSM 主导型水体的Chl-a 反演模型，结果显示：针对不同类别水体构建的反演模型的精度较分类前有了明显提高，且TSM 主导型水体构建的反演模型的MRE 与Chl-a 主导型水体的结果较为接近，证实了在高TSM 浓度的海域，采用CEI 指数法能够有效削弱TSM 对Schl的影响。