1 引言(Introductioin)

随着经济的快速发展, 我国湖泊富营养化问题日趋严重, 导致一些湖泊蓝藻大量繁殖, 水华频发.蓝藻水华暴发会严重降低湖泊的使用功能(陈云等, 2000), 藻类死亡产生的湖泛也时常发生, 藻毒素也可能由食物链损害人类身体健康, 这些都威胁着周围地区的饮水安全, 因此, 进行水华的监测和预警具有重要意义.我国湖泊水华的主要藻种是蓝藻, 其标志性色素是藻蓝蛋白, 但长期以来, 所有真核藻类的共同色素叶绿素a常作为蓝藻的量化色素(马荣华等, 2009), 并不能准确实现蓝藻水华的监测和预警, 因而开展藻蓝蛋白遥感反演用于蓝藻水华预测预警具有重要意义.

藻蓝蛋白在620 nm处具有强吸收特性, 根据这一指示性光谱特征, 国内外学者已开展了一些藻蓝蛋白浓度遥感反演研究.例如, Dekker等(1993)提出了一种基线算法, Schalles等(2000)提出了一种反射比算法, Simis等(2005；2007)建立了一种嵌套比值模型, Hunter等(2008)将三波段模型应用于藻蓝蛋白反演, 马荣华等(2009)建立了六次多项式的反演模型, 段洪涛等(2012)针对太湖在嵌套比值模型基础上重新率定各参数.由于内陆水体光学特性具有很强的区域和季节变异性, 反演模型通常具有区域和季节局限性.为了克服这个困难, 一些研究提出了先进行光谱分类, 再反演水质参数的方法(Zhang et al., 2015), 例如, 周晓宇等(2013)结合水体光学分类实现了太湖悬浮物浓度的反演, 李渊等(2013)基于光学分类建立了悬浮物反演模型.上述研究证明, 光学分类能够提高悬浮物反演的精度和稳定性, 但目前还没有采用基于光谱分类的反演方法反演藻蓝蛋白浓度的研究.

因此, 本文将构建基于光谱分类的太湖藻蓝蛋白浓度反演方法, 采用先分类再反演的策略, 基于水面反射率光谱分类进行太湖藻蓝蛋白浓度反演建模.首先采用逐步迭代的K均值聚类方法实现光谱分类, 然后分别利用每一类的训练样本光谱数据建立最适用于该类的藻蓝蛋白反演模型, 最后利用每一类的检验样本光谱数据进行反演模型精度评价.为了对比, 同时采用不分类的传统方法进行反演建模和精度评价.

2 研究区和数据(Study area and data) 2.1 研究区概况

本文以太湖为研究区.太湖位于长江三角洲南缘, 横跨江、浙两省, 水域面积为2338.1 km², 平均水深1.89 m, 是我国第三大淡水湖.近年来, 随着周边城市经济快速发展, 太湖水质恶化, 富营养化程度严重, 蓝藻水华频繁爆发, 严重影响着周边地区的生态环境和人民用水安全.太湖水质空间变化很大, 北部和西部富营养化最为严重, 东南部水质较好；太湖水质的季节变化也很大, 主要是由藻类水华的“休眠、复苏、生物量增加、上浮和积聚形成水华规律”决定的(孔繁翔等, 2005).

2.2 数据获取

本研究于2014年10月23—24日在太湖开展了水面采样测量实验, 40个采样点覆盖太湖大部分水域, 具体如图 1所示.

每个采样点测量了水面反射率光谱, 并采集了水样.采用水面以上法测量水面遥感反射率(R_rs)光谱(唐军武等, 2004), 使用野外光谱辐射计ASD Field Spec RS³, 测量光谱范围为350~1000 nm, 光谱分辨率3 nm.采集的水样放入水样箱低温保存, 带回实验室进行水质参数测量.

藻蓝蛋白浓度测量时先使用0.7μm孔径的Waterman GF/F玻璃纤维滤膜过滤水样；再在避光条件下, 加入50 mmol·L^-1、pH=7的磷酸盐缓冲液对滤膜进行研磨和搅拌；之后放入离心管中, 在4℃、15000 r·min^-1条件下离心25 min；然后将样品放在4~9℃的冰箱内避光保存不少于24 h, 再次进行离心；最后取上清液, 使用三维荧光仪测定其荧光强度, 利用标准曲线得到样品的藻蓝蛋白浓度.其中, 标准曲线是由藻蓝蛋白标样用磷酸盐缓冲液稀释至标准系列浓度绘制而成(阎荣等, 2004).对于叶绿素a浓度, 采用热乙醇法测定.

3 基于光谱分类的反演建模(Modeling based on spectra classification) 3.1 光谱分类

太湖属浅型湖泊, 水体组分复杂, 其中的浮游植物、非藻类颗粒物和黄色物质均会不同程度影响藻蓝蛋白的反演精度.简单地利用季节、区域并不能准确反映出太湖水体光学特征差异, 本文通过光谱分类将光学性质相似的水体归为一类, 然后建立适用于每一类的反演模型.本文针对水面实测反射率光谱, 采用K均值聚类方法进行光谱分类进而实现水体光学分类.K均值算法目前作为一种聚类分析方法被广泛应用(周爱武等, 2011), 属于非监督学习.相对于欧式距离, 光谱角度距离(SAD)更加注重两个向量在方向上的差异, 本文分类时距离函数采用SAD, 具体如公式(1)所示.其中, A(a₁, a₂, …a_n)、B(b₁, b₂, ……b_n)为两条光谱向量；a_n和b_n分别为光谱A和B第n个波长的遥感反射率.

(1)

传统K均值聚类结果对初始聚类中心敏感, 为了解决这个问题, 本文在K均值方法基础上采取逐步迭代的K均值聚类方法, 即先确定一个较大的类别数N, 完成一次聚类后, 计算各类别中心间的距离, 将距离最近的两类合并后重新聚类, 重复这个过程直至聚类结束.

3.2 藻蓝蛋白反演模型

目前已有反演模型中, 波段比值、基线和三波段算法均以藻蓝蛋白在620 nm处的强吸收为理论依据, 模型精简, 已被多个研究证实可较好地应用于内陆水体藻蓝蛋白反演.

3.2.1 波段比值算法

Schalles和Yacobi(2000)提出了一种简单的波段比值模型(式(2)), 其中, PC(Pigment Phycocyanin)为藻蓝蛋白浓度(μg·L^-1), 650 nm处为水体光谱曲线的反射峰, 625 nm处为藻蓝蛋白吸收形成的反射谷.

(2)

3.2.2 基线算法

Dekker(1993)通过10个浅水富营养化湖泊得出一种算法, 利用600、648、624 nm处的遥感反射率, 在水体光谱中以600、648 nm连接的直线为基线, 计算624 nm处的光谱曲线到这条直线的距离, 进而建立与藻蓝蛋白浓度的关系(式(3)).

(3)

3.2.3 三波段算法

Gitelson等(2003)提出的适应于陆地植被叶绿素含量遥感反演的三波段模型可适用于湖泊水体中的叶绿素a浓度的反演, 该模型又被Hunter等(2008)进一步应用于藻蓝蛋白的反演(式(4)), 取得了良好效果.其中, λ₁、λ₂和λ₃一般选择620、650~710 nm、740 nm.

(4)

3.3 反演建模策略

本文采用以下两种反演建模策略：①基于光谱分类的反演建模策略, 首先采用逐步迭代的K均值聚类方法实现光谱分类, 然后分别利用每一类的训练样本光谱数据建立最适用于该类的藻蓝蛋白反演模型；最后利用每一类的检验样本光谱数据进行反演模型精度评价.②不分类的传统反演建模策略, 为了对比, 同时采用不分类的传统方法进行反演建模和精度评价.训练样本和检验样本与前面光谱分类反演建模策略是一致的.

3.4 精度评价方法

采用公式(5)和公式(6)所示的均方根误差RMSE和平均相对误差σ对模型进行精度评价.

(5)

(6)

式中, Y表示藻蓝蛋白浓度预测值(μg·L^-1), X为藻蓝蛋白浓度实测值(μg·L^-1).

4 结果分析(Results analysis) 4.1 光谱分类结果

本次实验数据有40条光谱曲线, 初始的类别数N取15, 综合考虑光谱形状和数量, 将最终分类数目定为2, 分类结果如图 2所示, 分为第一类的有11个点, 第二类有29个点, 对各类水体的水质参数统计如表 1所示.

由光谱曲线可以判断, 第二类水体为浮游植物浓度相对较高的水体, 其在625 nm和670 nm附近有明显的反射谷, 700 nm附近有明显的反射峰.而这些特征在第一类的光谱曲线中并不明显, 第一类为浮游植物浓度相对较低的水体.为定量化确定分类效果, 计算各点705 nm与675 nm的遥感反射率波段比值(R_rs(705)/R_rs(675)), 并统计两类数据, 结果见表 1.根据表 1的统计数据, 第二类水体的藻蓝蛋白和叶绿素a浓度均高于第一类.第二类的R_rs(705)/R_rs(675)值明显高于第一类, 基本是可以利用1.0作为分类阈值, R_rs(705)/R_rs(675)>1.0的为第二类水体, 否则为第一类水体.

4.2 反演建模结果

本次实验数据共40个样点, 建模数据与验证数据比例为3:1, 模型波段选择和参数率定均利用建模数据完成, 传统不分类反演、基于光谱分类反演所用的建模数据和验证数据一致.

针对已分类水体构建反演模型, 首先分析每类水体的光谱特征.选择窗口大小为15, 多项式阶数为2的Savitaky-Golay(Ruffin et al., 2008)算法对光谱曲线做平滑处理去除噪声, 求曲线的一阶导数、二阶导数, 统计其大于或小于0的频率, 分析水体的光谱特征, 得到特征波段(申茜等, 2011).在一定范围内, 特征波段附近的波段组合其拟合效果较好, 结合拟合度最大原则, 得到的未分类数据和分类的类别对应的算法模型如表 2所示.表 2中列出了每种反演模型所使用的波段组合、模型公式、训练样本建模时的拟合度等.

各类别反演模型所用波段和参数均不相同, 波段比值模型中, 不分类数据选择650 nm和630 nm波段, 拟合度仅为0.28；第一类数据选择646 nm和629 nm波段, 拟合度为0.71；第二类数据选择654 nm和620 nm波段, 拟合度为0.61, 第二类数据为蓝藻主导的水体, 这个结果也符合藻蓝蛋白在620 nm处强吸收的理论.三波段算法中, 不分类的数据波段组合为642、698、767 nm, 拟合度为0.41；第一类数据波段组合为629、646、744 nm, 拟合度为0.61；第二类数据选择波段为610、698、759 nm, 拟合度为0.53.对于基线算法, 不分类数据的3个波段为600、630、650 nm, 拟合度为0.25；第一类数据选择波段为600、630、650 nm, 拟合度为0.27；第二类数据波段组合为600、620、665 nm, 拟合度为0.19.基线算法在3个类别中差别不大, 但拟合效果都不是很好, 波段比值和三波段算法的各类别波段组合具有较大差异性, 且分类后拟合度大幅度提高, 证明了分类的必要性.

4.3 反演结果分析

为验证本文反演方法的有效性, 利用验证数据, 根据公式(6)、(7)计算未分类和分类后反演结果的均方根误差RMSE和平均相对误差σ, 结果见表 2.

除基线算法在第二类数据的反演结果比未分类的反演结果差, 其他反演结果分类后精度均有明显提高.波段比值算法在未分类数据平均相对误差为156%, 在第一类数据中为50.8%, 第二类数据中为31.4%；对于三波段算法, 在未分类数据中的平均相对误差为93.5%, 在第一类数据中为31.2%, 第二类数据中为42.8%；基线算法在未分类数据中的平均相对误差为59.1%, 在第一类数据中为31.4%, 第二类数据中为66.2%, 略低于未分类数据.第一类数据中三波段算法平均相对误差为31.2%, 基线算法平均相对误差为34.9%, 第二类数据中波段比值算法的平均相对误差达到了31.4%.这一结果在藻蓝蛋白反演中较为理想, 光谱分类后反演结果明显优于未分类反演结果, 证明光谱分类能剔除不同类型水体差异所带来的影响, 提高藻蓝蛋白反演精度.

综合拟合度与验证误差, 不分类最好的模型为基线算法, 第一类最好的模型为三波段模型, 第二类最好的模型是波段比值模型.根据表 1, 第一类水体中, 叶绿素a浓度相对较低, 三波段算法算法能够较好地减少无机颗粒物及CDOM对藻蓝蛋白反演的影响.第二类水体中, 叶绿素a浓度相对较高, 通过波长相邻的两个波段的反射率比值能够抵消叶绿素a浓度的影响.

对未分类的数据采用基线模型, 第一类数据采用三波段模型, 第二类数据采用波段比值模型进行反演预测, 得到的预测结果如图 3所示.基于光谱分类反演的反演结果和实测数据拟合度达到0.87, 未分类反演模型的拟合度仅有0.46.对分类后的数据做总的精度评价得出RMSE=8.47μg·L^-1, σ=31.3%, 未分类数据的RMSE=14.14μg·L^-1, σ=59.1%.

综合上述结果, 藻蓝蛋白反演应采取先分类后反演的策略, 如图 4所示, 对于第一类数据, 光谱分类后适宜采用三波段算法反演藻蓝蛋白, 第二类数据在分类后适宜采用波段比值算法反演藻蓝蛋白.

4 结论(Conclusions)

本文利用太湖实测光谱数据和水质参数数据, 实现了基于光谱分类的太湖藻蓝蛋白反演, 建立了基于光谱分类的太湖秋季藻蓝蛋白反演模型, 并进行模型精度评价和验证, 得到以下结论：

1) 经过检验, 基于光谱分类的反演建模方法反演精度为0.87, 能够有效提高藻蓝蛋白反演精度.

2) 本文采用逐步迭代的K均值方法实现光谱分类, 将本次实验所获取的太湖水面遥感反射率光谱分为两类, 经过计算可以发现, 利用R_rs(705)/R_rs(675)的阈值分割也可以实现同样的分类效果：大于1的为第二类, 小于等于1的为第一类, 证明同一季节太湖水体具有光学差异性

3) 经过验证分析, 适用于太湖第一类水体的最优藻蓝蛋白反演模型是三波段模型, 适用于太湖第二类水体的最优藻蓝蛋白反演模型是波段比值模型, 不同类别水体反演模型不同, 有必要进行分类反演.

4) 本文采用的基于光谱分类的反演建模方法适用于秋季太湖藻蓝蛋白反演, 可为今后遥感反演太湖藻蓝蛋白浓度的时空分布提供基础, 对蓝藻监测具有重要意义.

5 展望(Prospect)

本文的方法在其他季节的适用性仍然需要更多的实验数据支持验证, 今后将在此基础上进一步完善模型.对于其他的内陆水体, 可以应用这种反演策略根据光学特性进行分类, 再建立反演模型, 今后将应用于其他内陆湖泊的藻蓝蛋白反演.在今后的研究工作中, 会将此方法进一步应用于高光谱影像, 对于位于两类水体边缘位置的水体在反演过程中可能出现“跳变”的问题, 今后将采用分类加权平均的方法进行改进.