近几十年来, 人类活动对海洋生态环境的影响越来越大, 某些近海海域的生态系统受损, 生态功能减弱, 水质急剧恶化。因此, 监测近海海域水质状况, 评估近海海域生态环境, 实现近海海域的生态修复与可持续发展已迫在眉睫。叶绿素a的质量浓度(下文简称叶绿素a浓度)是评估水域初级生产力和营养状况的关键指标, 也代表了藻类和浮游植物的关键特征, 对近海海域叶绿素a浓度进行有效监测有助于预防近岸水质的持续恶化^[1-4]。

叶绿素a浓度监测的传统方法过程繁琐, 而在线原位监测等新型监测手段受时间和空间的局限, 导致难以快速、动态、大面积评价水体叶绿素a浓度水平^[5-6]。近几十年, 遥感技术迅猛发展, 使用遥感技术对叶绿素a浓度进行监测, 具有快速、覆盖范围广、实时动态的巨大优势^[7-9]。从遥感信息中获取叶绿素a浓度和水体光谱反射率间的关系, 对水体的叶绿素a浓度进行反演, 可以实现对水体水质高频次、大范围的实时监测。经过近几十年的研究和探索, 叶绿素a浓度遥感反演在国内外取得了长足进展。根据反演模型构建的原理不同, 叶绿素a浓度反演模型可以分为经验、半经验模型、半分析和分析模型^[10-14]。目前大洋等一类水体的蓝绿波段比值法已取得较好的效果, 实现了业务化推广, 但此类算法在光学特征复杂的二类水体或者近海海域和内陆中应用受限, 因为在这些水体中光谱特征与组分浓度之间难以用简单线性关系表示^[15]。为解决此问题, 发展出了一系列更适用于二类水体的算法, 如红光-近红外二波段算法、三波段算法和四波段算法等^{[13-14, 16-19]}; 杜耀^[20]以深度学习为模型, 通过时空谱融合算法得到高分辨率遥感影像为数据源, 成功反演了白洋淀的水体叶绿素a浓度。但是大部分反演算法普适性差, 容易受到区域和季节性影响, 少数适应性较好的算法也有待进一步验证。

径向基(Radial Basis Function, RBF)神经网络能够近似各种非线性函数, 应对系统中难以明确规律的模式, 并展现出出色的泛化性能, 其学习速度快且能迅速收敛^[21], 已在非线性函数逼近、时间序列研究和信息处理等领域取得了成功应用^[22-23]。然而传统的RBF神经网络参数优化方式在训练速度方面存在局限, 且易受局部最小值困扰。已有学者研究并运用优化的RBF神经网络模型对二类水体的悬浮物浓度进行了成功反演, 例如, 谢旭等^[15]应用了高分一号(GF-1)遥感影像数据, 对粒子群最优化算法(Particle Swarm Optimization, PSO)优化后的RBF模型进行训练, 成功反演了闽江下游水体悬浮物浓度的变化情况; 吕君伟等^[21]基于2009年3月份ETM+影像数据和香港环保署实测采样点数据, 建立了光谱反射与悬浮物浓度之间的PSO-RBF反演模型, 与线性模型相比, 预测精度有了明显的提升。

本文选取香港近海海域为研究区域(图 1), 以实测叶绿素a浓度和对应采样点的遥感反射率为数据基础, 利用动态K-means聚类算法确定RBF神经网络的隐节点中心, 使用PSO算法对RBF神经网络的隐节点宽度和连接权值进行优化求解, 从而优化RBF神经网络参数, 建立叶绿素a浓度的反演模型。运用Sentinel-2多光谱成像仪(Multispectral Imager, MSI)获取的遥感影像数据进行反演, 并与传统RBF神经网络模和传统经验模型进行比较, 为今后香港近海海域叶绿素a浓度反演提供参考。

1 数据与方法 1.1 数据获取与预处理 1.1.1 实测数据获取

香港位于南海北部, 其近海海域向南与南海相通, 西侧是珠江口, 东边是一个半封闭的海域, 地理范围22°9′~22°37′N, 113°49′~114°27′E(图 1)。由于其地理位置和周围环境的特点, 香港近海海域受到陆地和海洋综合污染的明显影响, 其水质属于二类水体^[24]。本研究中所采用的香港近海叶绿素a浓度实测数据来源于香港环保署的官方网站(http://cd.epic.epd.gov.hk/), 这些数据覆盖了10个水质控制区, 并每月在香港的76个开放水域监测站进行记录。为确保数据质量, 选择了与Sentinel-2遥感影像日期相匹配、天气晴朗的数据, 叶绿素a质量浓度的单位为μg∙L^–1。

1.1.2 遥感影像数据获取及处理

Sentinel-2是搭载了MSI的高分辨率多光谱成像卫星, 由Sentinel-2A(S2A)和Sentinel-2B(S2B)两颗卫星组成, 轨道彼此相差180°, 覆盖13个工作波段, 幅宽达290 km, 放置在地球的两端。图像采集周期为每颗卫星10 d, 双星联合观测实现了5 d的时相分辨率, 空间分辨率分别为10 m、20 m、60 m, 卫星数据如表 1所示。通过欧洲航天局哥白尼数据中心(https://scihub.copernicus.eu/)下载了2019—2022年研究区域上空云覆盖率低于30%且晴朗的Sentinel-2影像, 共204幅, 选择其中云覆盖率低于10%的影像作为建模研究数据源, 共27幅, 其中18幅为Sentinel-2A, 9幅为Sentinel-2B(表 2)。

本次研究中使用的是Sentinel-2的L2A级遥感影像数据, 已经过欧空局哥白尼数据中心预处理, 在L1C级基础上完成了大气校正、辐射定标、几何校正等步骤, 只需要进行重采样、波段融合以及波段反射率提取。首先使用由欧洲航天局开发的哨兵应用平台(Sentinels Application Platform, SNAP)中的Raster-Resampling工具包, 调整重采样的分辨率需求参数为10 m, 运用了最邻近像元法进行重采样。对重采样后的遥感影像数据, 采用ENVI遥感影像处理软件(Environment for Visualizing Images, Harris Geospatial Solutions公司)的Layer Stacking插件, 选择B1至B12共12个波段(第8个波段有B8和B8a两种, 分别对应S2A和S2B卫星, 表 1所示), 利用最邻近像元法进行波段融合。最后利用ENVI软件中Input Points from ASCII把采样点的位置信息加载到波段融合后的图像中, 选择Output ROIs to ASCII得到与实测点对应的各个波段遥感反射率。数据获取后, 剔除缺测、异常值及受云影响的数据, 最终得到52组数据, 每组数据包含对应监测站点的实测值和12个波段的遥感反射率。

1.1.3 遥感反射率与叶绿素a浓度相关性分析

根据叶绿素a的光谱特性, 如图 2所示^[25], 在蓝光区域(约430~480 nm)和红光区域(约640~680 nm)有2个主要吸收峰, 在绿色区域(约520~570 nm)和685~ 700 nm附近则有两个吸收谷^[26-27]。而Sentinel-2 MSI传感器设置的B1至B4波段波长与其对应, 其中心波长分别为443 nm、490 nm、560 nm、665 nm。进一步分析水体反射率与叶绿素a浓度的相关性, 用52组数据做不同叶绿素a浓度水体的Sentinel-2 MSI波段光谱特性, 结果如图 3所示。综上研究表明, Sentinel-2 MSI传感器提供B1至B4波段的遥感数据信息, 对叶绿素a浓度相对比较敏感, 用于叶绿素a浓度遥感反演有较高的可靠性。

结合海洋叶绿素a浓度反演研究的经验, 对Sentinel-2 MSI B1—B4波段进行组合运算^[28-29], 根据表中的个数, 应产生153种组合类别(表 3)。并对153种组合与叶绿素a浓度做Pearson相关性分析, 提取达到强相关的波段组合, 结果如表 4所示: 波段组合中B2、1/B2、1/B3、B2-B4相关性较好, 相关系数都高于0.78, 其余波段组合相关性相对较差。

1.2 研究方法 1.2.1 RBF神经网络

RBF是一种具有径向对称衰减的函数, 把RBF应用于神经网络的隐含层作为激活函数就构成了RBF神经网络, 最早由Broomhead和Lowe在1988年建立, 主要用于模式识别、函数逼近和数据聚类等任务^[30]。本次研究中, 以实测点遥感反射率为输入层, 对应实测点的叶绿素a浓度值为输出层, RBF神经网络结构如图 4所示。

输入层到隐含层需要激活函数, 由于高斯基函数具备表现形式简单、径向对称、光滑性好、易于进行理论分析等优点^[31], 以高斯基函数作为激活函数, 其公式为

$ \mathop \varphi \nolimits^i (t)\left\| {{X_t} - {\phi _i}} \right\| = \exp \left( { - \frac{{{{\left\| {{X_t} - {\phi _i}} \right\|}^2}}}{{2\sigma _i^2}}} \right) . $

(1)

输出层可与隐含层节点加权通过线性变化得到, 假设输出层节点数为p, 其公式为

$ {Y_p} = \sum\nolimits_{i = 1}^m {{W_i}\exp \left( { - \frac{{{{\left\| {{X_p} - {\phi _i}} \right\|}^2}}}{{2\sigma _i^2}}} \right)} , $

(2)

式中, ${\phi _i}$为第i个隐节点中心, σ_i为第i个隐节点宽度(即扩展常数), W_i为第i个隐节点连接权值, X_p为第p个输入样本向量, m为隐含层隐含节点数及隐单元个数, Y_p为与X_p对应的输出样本值。

从式(1)、(2)可以看出, 实际使用RBF神经网络中, 以高斯基函数作为激活函数对神经网络性能的影响并不是至关重要的, 关键是如何确定各隐节点中心$ {\phi _i} $、各隐节点宽度σ_i以及各节点连接权值W_i。因此, 本文运用了PSO算法和动态K-means聚类算法, 对这3个参数的确定进行优化。

1.2.2 动态K-means聚类算法

传统K-means聚类算法需要人为确定初始聚类中心, 但不同的初始聚类中心可能导致完全不同的聚类结果^[32], 而动态K-means聚类算法具有原理简单、模型清晰、操作方便、计算快速等特点, 可以大规模同时对多种类型的数据进行聚类, 快速挖掘出数据中隐含的关系和结构, 并且有效解决了传统K-means聚类算法的局限性问题, 具体算法步骤如下。

假定预设距离限值η, 输入的第一个样本向量作为第一个聚类中心(即隐节点中心), 之后每当输入一个新的向量, 计算该向量与各聚类中心的欧氏距离:

$ d = \left\| {{X^p} - {C_j}} \right\| = {\left[ {\sum\nolimits_{i = 1}^n {{{\left( {X_i^p - {C_{ji}}} \right)}^2}} } \right]^{\frac{1}{2}}} , $

(3)

式中, n是输入向量的维数; C_j为第j个聚类中心; X^p为第p个输入向量。

根据式(3)确定了到每个聚类中心的欧式距离, 选择欧式距离最近的聚类中心, 此时有两种情况: 当d≤η, 把新输入向量归入该聚类中心, 并计算该聚类中心总输入向量的平均值, 对中心进行更新。当d > η, 把新输入向量作为新增聚类中心。当所有向量输入后, 动态K-means聚类算法结束, 各聚类中心确定。

1.2.3 粒子群最优化算法

受到鸟群觅食行为的规律启发, Kennedy和Eberhart于1995年建立了一个简化算法模型, 经过多年改进最终形成了粒子群优化算法(Particle Swarm Optimization, PSO), 也可称为粒子群算法^[33]。假设在D维搜索空间中, 有N个粒子, 每个粒子代表一个解, 则第i(1≤i≤N)个粒子信息为:

第i个粒子空间位置: X_id = (x_i1, x_i2, $ \cdots $, x_iD);

第i个粒子速度(是矢量, 包含粒子移动的距离和方向): V_id = (v_i1, v_i2, $ \cdots $, v_iD);

第i个粒子搜索到的最优位置(个体最优解): P_{id, pbest}= (p_i1, p_i2, $ \cdots $, p_iD).

通过优化目标函数计算每个粒子的适应值, 由适应值确认个体最优解P_{id, pbest}和群体最优解P_{d, gbest} = (p_{1, gbest}, p_{2, gbest}, $ \cdots $, p_{D, gbest}), 并根据判断条件决定是否按照式(4)、(5)来更新粒子的速度和空间位置:

$ v_{id}^{k + 1} = \omega v_{id}^k + {c_1}{\gamma _1}\left( {p_{id, {\text{pbest}}}^k - x_{id}^k} \right) + {c_2}{\gamma _2}\left( {p_{d, {\text{gbest}}}^k - x_{id}^k} \right) , $

(4)

$ x_{id}^{k + 1} = x_{id}^k + v_{id}^{k + 1} , $

(5)

式中, k表示迭代次数; ω表示惯性权重; C₁表示个体学习因子; C₂表示群体学习因子; γ₁, γ₂表示在区间[0, 1]内的随机数, 增加搜索的随机性; v_id^k表示粒子i在第k次迭代中第d维的速度向量; $x_{id}^k$表示粒子i在第k次迭代中第d维的位置向量; $p_{id, \;{\text{pbest}}}^k$表示粒子i在第k次迭代中第d维的历史最优位置, 即在第k次迭代后, 第i个粒子搜索得到的最优解; $ p_{d, \;{\text{gbest}}}^k $表示群体在第k次迭代中第d维的历史最优位置, 即在第k次迭代后, 整个粒子群体中的最优解。

而ω不宜为一个固定的常数^[34], 因此提出了自适应调整的策略, 即随着迭代的进行, 线性地减小ω的值, 其公式为

$ \omega = {\omega _{\max }} - ({\omega _{\max }} - {\omega _{\min }})\frac{{iter}}{{ite{r_{\max }}}} , $

(6)

式中, ω_max、ω_min分别表示最大和最小惯性权重; iter表示当前迭代次数; iter_max表示最大迭代次数。当达到规定的最大迭代数, 或是相邻更新后, 最优解的适应值之间的误差达到指定误差范围内即停止算法优化。

1.2.4 优化RBF神经网络模型

由1.2.1节可知, 影响RBF神经网络性能的主要参数有3个: 各隐节点中心$ {\phi _i} $、各隐节点宽度σ_i以及各节点连接权值W_i。本文先采用动态K-means聚类算法获取各隐节点中心, 再利用PSO算法对RBF神经网络的隐节点宽度和连接权值进行优化求解, 由RBF神经网络均方根误差作为优化目标函数以计算适应值, 如下式所示

$ \text{fitness} (i) = \sqrt {\frac{1}{N}{{\sum\nolimits_{k = 1}^N {\left[ {{y_k}(i) - {{\hat y}_k}(i)} \right]} }^2}} , $

(7)

式中, fitness(i)表示第i个粒子自适应值, N表示输入样本数量, y_k(i)表示第i个粒子的实测值, $ {\hat y_k}(i) $表示第i个粒子预测值。当迭代停止, 自适应值最小, 训练得到的各隐节点宽度σ_i以及各节点连接权值W_i最优。优化RBF反演模型流程图如图 5所示。

2 结果与分析

将处理好的52组数据分为两组, 随机选取35组作为建模数据, 剩余17组作为验证数据。经1.1.3节分析可知, 优化RBF模型选取B2、1/B2、1/B3、B2-B4作为模型输入层; 输出层为对应采样点的叶绿素a浓度值。设置粒子群规模N为50, 迭代次数k为1 000, 学习因子C₁ = C₂ = 2, 惯性权重ω_max = 0.9, ω_min = 0.2; 隐节点中心和数目由动态K-means聚类算法自动获取, 预设距离限值η为10^–2时优化的RBF反演模型最佳。

为了深入展示优化RBF模型的优势, 使用同一组建模数据建立传统经验模型(表 5)和传统RBF神经网络模型, 并分别与17组验证数据进行误差分析。结果如图 6所示, 从决定系数(R²)来看, 优化RBF模型最高, 为0.90, 其次是自然对数模型的0.75, 传统RBF模型的决定系数最小, 为0.38; 从均方根误差(RMSE)来看, 优化RBF模型最小, 为0.23 μg∙L^–1, 自然对数模型的0.61 μg∙L^–1和传统RBF模型的0.89 μg·L^–1均大于优化RBF模型。实测值与不同模型预测值的误差分析中(表 6), 优化RBF模型的平均相对误差(MRE)为11.12%, 其中相对误差低于20%有14个, 占总验证数据的82.35%; 自然对数模型和传统RBF模型的平均相对误差为33.94%和42.82%, 其中, 相对误差低于20%的自然对数模型有9个, 占总验证数据的52.94%, 传统RBF模型有6个, 占总验证数据的32.59%。

如图 7所示, 优化RBF模型和自然对数模型的验证值趋势都与实测值大致相似, 但优化RBF模型的验证值趋势与实测值更相近, 并且没有相对过高或过低的值出现, 整体反演较为平稳。而传统RBF模型只有少许点趋势相符合。无论从建模精度分析还是误差分析, 优化RBF模型反演效果优于传统经验模型和传统RBF模型。

3 讨论

基于2019—2022年香港近海海域的204幅Sentinel-2 MSI遥感影像, 运用训练好的优化RBF模型反演香港近海海域叶绿素a浓度, 其模型输入层为B2、1/B2、1/B3、B2-B4, 模型输出层为影像反演的叶绿素a浓度值。

图 8展示了香港近海海域2019—2022年多年平均叶绿素a浓度的空间分布, 近海海域整体呈现出明显的东西分布差异的特征, 最东侧和最西侧远海海域叶绿素浓度低, 中间近岸区域叶绿素浓度高, 叶绿素浓度分布存在地域性差异, 整体的反演结果低于6 μg·L^–1。根据香港环境保护署官网(https://www.epd.gov.hk/epd/sc_chi/environmentinhk/water/hkwqrc/waterquality/marine-2.html)提供的2019—2022年逐年海水水质年报中把香港近海海域划分为西部、中部、东部和南部, 本文中为便于对反演结果的计算、分析, 把南部海域的南区水质管制区划入中部海域, 并参照海水水质年报中的水质管制区编号。依照划分的海域计算叶绿素a平均浓度值, 其结果如表 7所示。西部缓冲区叶绿素a浓度最高, 为4.0 μg∙L^–1; 其次是东部缓冲区、牛尾海、吐露港及赤门、后海湾和将军澳较高, 都位于海湾内, 叶绿素a浓度分别为3.9、3.7、3.5、3.4和3.3 μg∙L^–1; 南区和大鹏湾较低, 两个水质区内大部分海域处于远海海域, 整体叶绿素a浓度分别为3.0和2.2 μg∙L^–1; 西北部叶绿素a浓度最低, 为1.8 μg·L^–1。从最值来看, 叶绿素a浓度最高值出现在吐露港及赤门的海湾内部沿岸区域, 为5.3 μg∙L^–1; 叶绿素a浓度最低值出现在西北部的远海海域, 为0.4 μg∙L^–1。而根据前人研究, 吐露港及赤门叶绿素a浓度平均差了1.0 μg∙L^–1, 牛尾海叶绿素a浓度平均差了0.5 μg∙L^–1。这是由于PSO-RBF模型仅采用站点的实测值进行评估分析, 而这两个水质区范围狭窄, 遥感数据提取容易产生偏移, 并且没有综合其他位置区域进行验证分析, 这也会导反演结果与实际结果有一定误差^[35]。从整体上看, 香港近海海域叶绿素a浓度整体分布从西到东呈现低-高-低的变化格局。从局部来看, 叶绿素a浓度在每一个海湾中由内湾向外海方向、从近岸到远海方向逐渐减少的空间变化规律。并且在小范围近海海域, 叶绿素a的浓度变化非常显著且变化多样。

就香港近海海域叶绿素a浓度整体分布从西到东呈现低-高-低的变化格局, 可能是因为西部的西北部水质控制区虽然受到珠江入流影响, 但是该区域水体交换较为迅速且陆源污染较少, 因此, 叶绿素a浓度整体上最小。对于西部的后海湾水质管制区, 中部的西部缓冲区、东部缓冲区和将军澳, 以及东部的牛尾海和吐露港及赤门水质管制区, 由于其海湾的地理环境且靠近城市, 海水中长期存在生活污水、工业废水等, 且水体交换较深圳湾控制区弱, 因此叶绿素a浓度升高。东部海域内部的牛尾海区以及吐露港与赤门海峡区是半封闭海湾, 尤其吐露港与赤门海峡区的“瓶颈状”地貌使其水体交换非常弱^[36-37], 导致叶绿素a浓度最高值常出现于该水域内, 这与前人研究结果基本一致。从局部来看, 叶绿素a浓度在每一个海湾中由内湾向外海方向、从近岸到远海方向逐渐减少, 这种变化规律反映了陆源污染对叶绿素a浓度分布的控制作用。近岸靠近城市, 水中含有大量生活污水、工业废水等高营养物质水体, 导致叶绿素a浓度升高; 而近岸到远海方向, 水体流动性逐渐增强, 与相邻海水的水体交换率提高, 减少了陆源污染物的影响, 其叶绿素a浓度也相对降低。

由于香港上空范围内符合实验要求的Sentinel-2遥感影像数据较少, 一年中大部分时间的影像数据云覆盖率高于10%, 导致获取的实验数据量较少, 这些因素都在一定程度上影响了反演精度。因此建立高精度、大样本的水质样本数据库, 对于近海水质监控是有必要的。总的来说, 优化的RBF模型基本反演出叶绿素a浓度分布规律, 具有较高的实用性, 也说明了优化的RBF模型在叶绿素a浓度反演中的应用价值。

4 结论

本文基于Sentinel-2 MSI遥感数据和香港近海海域叶绿素a浓度实测数据, 提出一种动态K-means与粒子群算法优化的RBF神经网络叶绿素a反演模型, 并与传统经验模型和传统RBF模型进行对比分析。得到以下结论。

(1) Sentinel-2 MSI传感器提供1—4波段的遥感数据信息, 对叶绿素a浓度变化相对比较敏感, 用于叶绿素a浓度遥感反演有较高的可靠性。且4个波段组合运算中, B2、1/B2、1/B3、B2-B4与叶绿素a浓度的Pearson相关性较为出色, 相关系数都超过了0.78。

(2) 在动态K-means与粒子群算法优化的RBF叶绿素a反演模型训练中表明, 建立4个输入层节点的优化RBF反演模型效果最佳, 其建模精度R²=0.90高于传统经验模型和传统RBF模型, 不同模型的误差分析中RMSE和MRE更低[RMSE=0.23(μg·L^–1), MRE=11.12 %]。

(3) 基于优化的RBF模型反演香港近海海域叶绿素a浓度分布显示, 整体从西到东呈现低-高-低的空间变化格局, 局部上呈现近岸高到外海逐渐降低的现象。

利用动态K-means与粒子群算法优化RBF模型中的参数, 结果显示该方法简单有效, 对叶绿素a浓度反演的模型研究具有一定的参考价值, 同时可为其他水质参数进行遥感反演提供参考。