林腾 陈哲锋 杨义炜
(福建省地质测绘院 福建福州 350011)
叶绿素a 作为水体中藻类植物和浮游植物最丰富的色素之一,其浓度大小可以在一定程度上反映水体的富营养化程度,是进行水体营养状态评价的重要参数之一[1-2]。对叶绿素a浓度进行监测是保护水生态系统的有效手段之一,是人类社会与自然相互协调发展的必然趋势。传统叶绿素a 浓度监测主要以地面水站监测或实地水样采集监测为主,该方法可以获得监测站点或水样采集点的较高精度的水质情况,但受限于监测范围,无法对整个水面进行监测,且该方法以人工为主,工作量大,时效性不高,难以满足实时、动态和广泛的需求[3]。遥感技术的实时性、连续性、大范围、低成本等优势,正好弥补了传统监测手段的缺点,成为目前叶绿素a 浓度监测最有效的手段之一[4]。
随着遥感传感器技术的进步和水体光谱特征研究的不断深入,叶绿素a 浓度遥感监测从定性发展到定量,模型算法不断改进,监测精度也不断提高[5]。国内外学者将遥感技术应用于叶绿素a 浓度监测的研究开始于20 世纪70 年代,GORDON等[6]利用蒙特卡罗模拟构建固有光学特性与表观光学特性之间的一般关系,为水质遥感的发展奠定基础;KEINER 等[7]基于Landsat TM 数据,利用神经网络对叶绿素和悬浮物浓度进行建模反演,取得不错的效果;郭宇龙等[8]以标准三波段算法为基础,结合实测叶绿素a 浓度和光谱数据,构建了三峡水库、巢湖、洞庭湖和太湖的叶绿素a 模型,结果表明COCI 三波段算法具有较强的应用潜力;朱利等[9]利用GF-1 WFV 和HJ-1A CCD 数据,对太湖进行叶绿素a、悬浮物、透明度监测和和富营养化评价,研究表明了2 种数据的水质反演结果的一致性,验证了GF-1 WFV 数据的水质监测潜力。
古田水库是古田县人工湖,又名“翠屏湖”,是福建省级风景名胜之一,以发电为主,兼顾旅游业,与古田县经济发展和生态环境息息相关[10]。近年来,水库周边城镇化进程日益加快,工农渔业生产迅速发展,生活污水、工业废水、农业污染的排放严重影响了水体水质,导致水库水质每况愈下。目前,针对古田水库水环境遥感监测的研究还相对较少,本文基于Sentinel-2 卫星遥感影像数据,结合同步实测叶绿素a 浓度数据对2021 年古田水库每季度的叶绿素a 参数开展遥感反演建模,并对结果进行时空特征分析,为古田水库水质监测和治理提供参考依据。
古田水库位于福建省古田县中部闽江支流古田溪上,水库控制流域面积1 325 km2,总库容6.47 亿m3,地理位置:118°46′E~118°52′E,26°34′N~26°43′N。研究区属于亚热带海洋性季风气候,气候温暖湿润,年平均气温在14~20 ℃之间,雨量充沛,3~6 月为梅雨季节,7~9 月多发台风。库区水资源丰富,水质保护目标为Ⅲ类,主要有玉源溪、柏源溪、芹溪、达才溪、横洋溪等支流流入水库。水库上游分布有平湖、新舫、桥洋、大桥、吉巷等多个乡镇,既有以食用菌生产加工为主的工业生产,又有农业耕种和渔业养殖,多方面因素对库内水体水质产生较大影响,水环境保护形势严峻。
本文采用的数据为欧空局的5 景L1C 级哨兵2 号(Sentinel-2)多光谱卫星遥感影像,覆盖光谱波段13 个,空间分辨率分别为10、20、60 m,获取时间分别为2021 年3、6、9、12 月和2022 年1 月,幅宽290 km,完整覆盖古田水库。
由于L1C 级数据是经过几何精校正的正射影像,并没有进行辐射定标和大气校正,在开展古田水库叶绿素a 浓度反演之前,要先对Sentinel-2 数据进行辐射定标、大气校正和影像裁切等预处理。辐射校正的精度直接影响影像的质量,本文借助SNAP 平台和欧空局(ESA)发布的Sen2Cor 插件对原始数据进行辐射定标和大气校正处理,并将空间分辨率统一重采样为10 m,再对处理后的数据进行裁切,最终获得研究区5 个时相的Sentinel-2 多光谱影像数据。
本文通过实地水质采样来测得古田水库叶绿素a 浓度数据,采样时间为2022 年1 月13 日—14 日的上午10 时至下午13 时,与Sentinel-2 影像(时相2022 年1 月13 日)实现准同步。研究区共采集样本20 组,剔除3 组异常数据,在剩余数据中随机选取13 组数据用于建立模型,其他4 组数据用作模型验证(图1)。
图1 古田水库采样点分布图
目前,遥感技术手段在叶绿素a 浓度监测领域的应用已经相对成熟,各种类型的反演算法相继被提出,常用的可归为3类:经验模型算法、半经验/半分析模型算法和分析模型算法。经验模型算法是在大量的样本基础上,基于多光谱遥感数据的波段反射率与实测水质参数的相关性统计分析,筛选最优的波段或波段组合来构建水质参数的回归反演模型;半经验模型算法是依托实测水面光谱和水质参数的基础上,利用数理统计等手段对两者进行反演建模;分析模型算法是在辐射传输理论的基础上,计算水质参数对光谱特定波段的吸收系数和散射系数,进而模拟水质参数据浓度[11]。考虑到研究区仅有少量实测叶绿素a 参数数据,缺少同步实测光谱数据和吸收、散射等光学数据,本文采用支持向量回归机(Support Vactor Regression,SVR)模型进行叶绿素a 参数反演。SVR 是支持向量机在函数回归中的应用,它是1 个非线性拟合函数的过程,有效解决了样本数量不足的问题,一定程度上也可以提高叶绿素a 浓度反演精度[12]。SVR 叶绿素a 定量反演是以叶绿素a 参数敏感波段为输入样本,利用交叉验证法寻找最佳模型参数来构建模型,并利用验证点对模型进行验证。
水质遥感反演的对象主要为水体水域。为了提高反演效率,避免植被、建筑、陆地等非水体数据参与反演计算,需要对古田水库进行水域信息提取。目前,常用的遥感水体信息提取方法较多,包括单波段法、多波段法、植被指数法和水体指数法等。本文采用归一化水体指数方法来提取水域边界信息。归一化水体指数(NDWI)[13]主要依据的原理是可见光至近红外波段内水体与非水体地物间存在的反射光谱差异,利用水体反射最强的近红外波段与水体反射最弱的绿波段的比值运算,进一步加大两者间的差异,最大程度地削减非水体信息而增强水体信息。NDWI 的计算如式(1)。
式中:Green 为绿波段反射率,对应Sentinel-2 影像第3 波段;NIR 为近红外波段反射率,对应Sentinel-2 影像第8 波段。
经上述运算处理后,可根据NDWI 图像直方图设置合适的阈值进行水体与非水体的划分,以提取水体信息。
在水质遥感监测中,一般通过分析单波段或波段组合与水质参数浓度的相关性,寻找各水质参数的敏感性波段。相关性分析是探究2 个变量之间的相互关系而建立的分析方法,分析结果可以反映2 个变量之间的密切程度,相关性系数的数值范围是从-1 到1 之间,绝对值越接近1,则两者的相关度越大[14]。本文采用Pearson 相关系数法进行相关性分析,其计算如式(2)。
式中:xi和分别为叶绿素a 的浓度及其算术平均值,yi和y分别为波长λ 的光谱反射率及其算术平均值;R(λ)为X 和Y之间在波长λ 处的相关系数。
SVR 模型的构建是在Matlab 软件平台中利用林智仁教授开发的libsvm 程序库,以训练样本的水质敏感波段和水质参数浓度为输入变量,选择径向基核函数为核函数,通过调节惩罚系数c、核函数系数g 和不敏感系数p 等模型参数来寻找效果最优的模型[15]。拟合度和均方根误差是判断模型是否最优的主要依据。
本文采用Pearson 相关分析法,对实测叶绿素a 浓度与2022 年1 月13 日获取的Sentinel-2 中400~900 nm 波段范围内各波段(对应Sentinel-2 的B1~B8 波段)的光谱反射率进行相关性分析,分析结果见表1。从表1 可以看出,叶绿素a 浓度与B3、B5-B8 的相关系数在0.6 以上,与B5 的相关系数最大,为0.720 2。
表1 Sentinel-2 单波段与叶绿素a 参数的相关性
根据黄灵光等[16]研究发现,遥感影像各波段之间存在信息冗余的情况,不能很好地反映其与水质参数的关系。采用波段组合可以在一定程度上消除各波段间的干扰,有效减少水中其他杂质的影响,突出水质参数与水体反射率的相关性。结合前人经验总结和波段相关性分析结果,选择Sentinel-2 中与叶绿素a 相关的B3、B4 和B5 波段进行比值、差值等组合运算,再通过波段组合与叶绿素a 浓度进行相关性分析,分析结果见表2。从表2 可以看出,叶绿素a 浓度与波段组合间的相关性明显高于单波段,叶绿素a 与(B5-B4)/B3 的相关性最高,相关系数达0.876 2。最终本文选择(B5-B4)/B3 波段组合为叶绿素a 参数的敏感波段。
表2 Sentinel-2 波段组合与叶绿素a 参数的相关性
本文根据叶绿素a 参数敏感性分析结果,选择(B5-B4)/B3作为叶绿素a 参数反演模型构建的输入变量,以径向基核函数为核函数,并采用交叉验证法来寻找拟合效果最佳的模型参数,最终确定叶绿素a 参数的最优SVR 反演模型。叶绿素a 参数的SVR 模型拟合效果如图2 所示。从图2 中可以看出,叶绿素a 反演模型的拟合效果较好,拟合度为0.874,均方根误差为2.691。
图2 叶绿素a 参数的SVR 模型拟合图
为了验证模型的准确性,本文利用SVR 反演模型对测试样本进行叶绿素a 浓度估算,采用相对误差指标对得到的估算值与实测值进行对比分析,从而对模型的精度进行检验(表3)。从表3 可以看出,模型整体精度较好,4 个测试点的相对误差均控制在30%以内,平均相对误差为13.32%,可认为满足精度要求。
表3 模型精度检验结果
为研究古田水库叶绿素a 浓度的时空变化,本文分别选取2021 年4 个不同季节的月份的Sentinel-影像数据,利用已建立的叶绿素a 参数的SVR 反演模型对古田水库进行遥感反演,并在ArcGIS 软件平台下进行制图输出,得到古田水库叶绿素a 浓度的时空变化分布图(图3)。
图3 古田水库叶绿素a 浓度2021 年不同月份时空变化分布图
从空间变化来看,古田水库叶绿素a 浓度在不同季节的4 个月份存在着显著的差异。3 月和6 月的叶绿素a 浓度整体偏高,库区内部叶绿素a 浓度略高于上游,呈现出南高北低的趋势;9 月和12 月的叶绿素a 浓度整体偏低,较3 月和6 月呈现的是下降的趋势,其中库区内叶绿素a 浓度下降明显,使得空间分布上叶绿素a 浓度从南向北递增,呈现上游明显高于下游的特征。
从季节性变化来看,古田水库叶绿素a 浓度春、夏2 季明显高于秋、冬2 季,且夏季整体达到最高,主要是由于温度对叶绿素a 的调节作用。古田水库底部的氮磷营养物质在上游城镇生活污水、畜禽农业污染和工业废水的长期影响下不断堆积,一旦天气变暖,水温升高,库内底泥中的氮、磷会被释放出来,水库中藻类大量繁殖,将导致库内叶绿素a 浓度上升。
本文选用大宽幅的欧空局Sentinel-2 多光谱卫星遥感数据,结合实测叶绿素a 参数浓度数据,基于SVR 模型构建古田水库叶绿素a 参数反演模型,并进行了精度检验。通过已建立的反演模型,对2021 年3、6、9、12 月等4 个月份的古田水库叶绿素a 浓度进行遥感监测,最终获取了古田水库4 个季节叶绿素a 浓度的时空变化分布特征,得出结论如下:
(1)古田水库叶绿素a 浓度与波段组合间的相关性明显高于单波段,叶绿素a 与波段组合(B5-B4)/B3 的相关性最高。
(2)采用测试样本对叶绿素a 参数的SVR 反演模型进行检验,结果表明,叶绿素a 参数的SVR 模型反演效果较好,拟合度高达0.8 以上。SVR 反演模型可以有效地应用于小样本的叶绿素a 浓度遥感监测。
(3)反演结果表明,春、夏2 季古田水库叶绿素a 浓度整体偏高,主要是受温度和降水的影响。空间分布上,3 月和6 月库区内部叶绿素a 浓度略高于上游,呈南高北低趋势,9 月和12月从南向北递增,呈现上游明显高于下游的特征。叶绿素a 浓度问题主要与上游乡镇生活排污、农业污染、畜禽养殖等人类活动密切相关。