沉水植物是湖泊生态系统的重要水生植物类型之一, 可以改善生态系统的多样性, 稳定沉积物、吸收营养盐、净化水体、维持渔业生产和抑制浮游植物生长等生态功能(Blindow et al, 2002; 童昌华等, 2004; 王卫红等, 2006; Liu et al, 2018)。因此, 大范围、准确、快速获得沉水植物的时空分布对湖泊水体治理和水生态修复具有十分重要的意义。

传统监测沉水植物分布与生长状况的现场调查法和潜水测量法耗时耗力、监测覆盖范围小, 而且影响水体环境(Barillé et al, 2010)。遥感技术作为监测水生植物的有效工具, 不仅成本较低, 省时省力, 而且能够大面积同步获取地物信息, 有着传统方法不可比拟的优势(王琪等, 2016)。在多光谱应用方面, 许多学者利用Landsat系列卫星遥感影像数据(林川等, 2010)和我国环境及高分卫星(Luo et al, 2002; Chen et al, 2018)遥感数据对水生植物进行遥感监测、初步分类、提取和制图, 并获得了很好的研究结果; 其中, Brooks等(2015)利用Landsat Thematic Mapper和Multispectral Scanner Systrm影像评估了劳伦斯大湖沉水植物的分布情况, 并指出沉水植物的分布有利于净化水质, 改善水体透明度。Traganos等(2018)基于Sentinel-2 Multi Spectral Instrument数据, 利用决策树方法绘制了希腊西北部爱琴海Thermaikos海湾的沉水植物分布图。Liu等(2015)基于MODIS数据, 利用floating algal index (FAI)指数和vegetation presence frequency (VPF)方法有效地绘制了太湖水生植物分布图。Dogan等(2009)利用Quickbird卫星遥感影像数据绘制了浅水湖泊中沉水植物覆盖度的分布图, 并指出冠层水深是影响遥感大规模监测沉水植物覆盖度和丰度的主要因素。近几年来, 随着高光谱遥感技术的应用与发展(张达等, 2013), 将高光谱遥感数据用于沉水植物覆盖度等方面的监测研究, 有效的提高了反演精度。例如, Pu等(2013)利用Hyperion高光谱遥感数据对美国佛罗里达州西部海岸的沉水植物覆盖度进行了反演, 结果表明, 当覆盖度为3级和5级时, 反演精度分别达到95.9%和78.4%。Zou等(2013)在富营养化湖中, 利用实测数据研究了沉水植物光谱特征, 并指出沉水植物光谱反射率随覆盖度的减少而下降。但是当前国内外将沉水植物的实测光谱特征与遥感推广应用相结合进行沉水植物覆盖度的综合研究还较少。

本论文以内蒙古巴彦淖尔市境内乌梁素海沉水植物群落为研究对象, 基于地物光谱仪测定的沉水植物光谱数据, 开展如下研究: 1)分析不同覆盖度的沉水植物光谱特征差异; 2)探讨覆盖度对沉水植物群落光谱特征的影响; 3)建立沉水植物覆盖度的遥感反演模型; 4)利用Sentinel-2 MSI遥感数据反演得到乌梁素海沉水植物覆盖度分布图。

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

乌梁素海(108°43′~108°57′E, 40°36′~41°03′N)(图 1)位于内蒙古巴彦淖尔市乌拉特前旗境内, 是地球上同一纬度最大的自然湿地。湖泊呈南北长、东西窄的狭长形态, 湖岸线长130 km。湖泊蓄水量2.5亿~3亿m³, 水深为0.5~3.0 m, 平均水深1.5 m, 水域面积273.32 km²。所在地区的多年平均气温为7.3 ℃, 多年平均降雨量为224 mm, 湖水于每年11月初结冰, 次年4月份解冻, 冰封期为5个月。从20世纪90年代起, 由于自然气候条件干旱, 乌梁素海的自然补水量持续减少, 再加上所属地区的生活污水、农业退水和工业废水, 未经处理就直接排放, 使湖泊水质下降, 生态功能急剧退化, 湖区内沉水植物等大型水生植物疯长, 乌梁素海成为了重度富营养化草型湖泊(于瑞宏等, 2004; 张岩, 2012; 于玲红等, 2014; 李建茹等, 2017)。乌梁素海作为内蒙古“一湖两海”的重要组成部分, 其一旦消失, 不仅会造成土地沙化, 而且我国北方的生态屏障将受到重创(岳程鹏等, 2020)。因此, 利用遥感大范围、快速和周期性监测乌梁素海沉水植被, 对于保护乌梁素海生境有着极其重要的意义。

乌梁素海水生植被生长茂盛, 在全湖范围内广泛分布, 植被种类大致可以分为挺水植被、沉水植被和漂浮藻类三类。其中, 沉水植被有龙须眼子菜、菹草、狐尾藻、茨藻、轮藻、竹叶眼子菜等, 龙须眼子菜和穗花狐尾藻为沉水植被优势种, 5月中旬开始萌芽, 10月开始逐渐衰退, 主要分布在明水区, 但是南部明水区分布较少(岳丹等, 2015; 包菡等, 2016; 卫亚星等, 2017; Cao et al, 2020)。

1.2 野外观测

本研究分别于2019年6月17日、7月19日、8月20日、9月21日在乌梁素海开展了4次为期2 d的野外调查, 共获得有效实测光谱数据218组(图 1), 沉水植物(64组)和开阔水体(154组)的遥感反射率光谱、沉水植物覆盖度和现场照片, 且选择以上月份为沉水植物生长期, 适用于沉水植物的覆盖度遥感反演研究。

光谱测量采用美国ASD公司生产的ASD Field Spec 4地物光谱仪, 其测量波段范围为350~2 500 nm, 光谱分辨率为3 nm, 光纤视场为25°, 350~1 000 nm内采样间隔为1.4 nm, 1 001~2 500 nm内采样间隔为1.1 nm。光谱测量方法采用唐军武等(2004)提出的水面以上水体光谱测量方法。遥感反射率R_rs (remote sensing reflectance)测量在无风、天气晴朗时进行, 且每个样点需要测量标准灰板、水体和天空光等辐射亮度信息, 各测15条光谱数据(即3组), 剔除异常数据之后, 取均值作为该点的光谱数据。为了避开水体对太阳直射反射和船体阴影等影响, 需要采用一定的观测角度, 仪器观测平面与太阳入射平面夹角为90°~145°之间, 且仪器与水体之间的距离为1 m。用利1 m×1 m的样方框现场测定沉水植物样方内覆盖度C (coverage), 在每个样方中, 使用细分框架和11个覆盖度等级(覆盖率为10%, 即0 (0%)、1 (1%~10%)和11 (91%~100%)来估算覆盖度。用塞氏盘测量水体透明度, 利用手持GPS定位仪记录每个站点的位置。

1.3 遥感影像数据及其预处理

本文选用2019年6月16日和2019年8月20日, 2景无云覆盖的Sentinel-2 MSI Level-2A级地表反射率产品(MSI-2A)影像数据(https://scihub.copernicus.eu/), 反演乌梁素海沉水植物覆盖度。

Sentinel-2卫星搭载的多光谱成像仪MSI, 共有13个波段, 不仅有可见光、近红外和短波红外波段, 还包含了3个“红边”波段, 和一个增强的窄近红外波段。Sentinel-2 A和B两颗卫星在过境时间上互补, 重访周期提升为5 d (Drusch et al, 2012; Tyler et al, 2016; Pahlevan et al, 2017)。

MSI-2A是利用欧空局提供的Sen2cor软件对Sentinel-2卫星影像进行大气校正后的地表反射率产品(Main-Knorn et al, 2015)。该算法的核心是大气辐射传输模型IibRadtran (Mayer et al, 2005)。IibRadtran是用于计算大气辐射传输的软件包, 该算法用离散纵坐标法计算地球大气、太阳以及热辐射传输, 包含有不同的计算模块, 适用于各种大气条件, 例如太阳和部分地球光谱的辐射、辐照度和光化通量(霍娟等, 2006; 苏伟等, 2018)。Sen2cor是基于IibRadtran代码生成的查找表, 并对Sentinel-2影像进行大气校正, 利用Python编译的可移植代码, 为用户提供命令和集成的软件(张晓月等, 2020)。

将MSI-2A数据的各个波段灰度值乘以0.000 1, 得到地表反射率, 并通过如公式(1)转化为遥感反射率。

式中, ρ_w是地表反射率, R_rs为遥感反射率。

其次, 对影像进行陆地掩膜, 裁剪, 并重采样MSI-2A产品, 得到统一的20 m空间分辨率。

另外, 本文试图利用Sentinel-2 MSI遥感数据进行沉水植物覆盖度的反演。针对上述传感器的波段设置, 进行了乌梁素海实测遥感反射率的波段积分运算, 具体见公式(2), 光谱响应函数如图 2所示。

式中, R_rs (λ)表示实测遥感反射率; F_i表示MSI数据的第i波段的光谱响应函数。

1.4 数据处理与分析

首先, 基于沉水植物光谱特征, 进行覆盖度与沉水植物光谱之间的相关分析, 确定相关性最大的单波段和波段比。相关分析同时考虑了皮尔逊(线性)和斯皮尔曼(单调性)相关分析法。

然后, 建立单波段/波段比与覆盖度之间的回归模型, 模型如下:

式中, x=R_rsi或R_rsi/R_rsj, a、b、c…为模型参数。

采用LOOCV (leave-one-out cross validation) (Wang et al, 2010)方法确定模型参数a、b等, 具体做法如下: 首先剔除第一组样本数据, 其余数据作为训练样本, 利用最小二乘法确定模型参数; 之后剔除第二组样本数据, 其余数据作为训练样本; 该迭代过程进行到最后一组样本的剔除为止。这个过程会产生一系列的模型参数a、b等, 取其平均值作为模型参数。

最后, 采用均方根误差e (root mean square error, RMSE), 平均相对误差m (mean relative percentage error, MRPE)以及Pearson相关系数评价反演模型精度, 计算公式为:

式中, x_ipre和x_imea分别表示通过模型计算反演得到的覆盖度和实际测量的覆盖度; n为采样点数。

2 结果与讨论 2.1 不同覆盖度的沉水植物光谱特征分析

沉水植物光谱曲线表现出典型的植被光谱特征(图 3), 在可见光部分的蓝波段(480 nm为中心)和红光波段(680 nm为中心)附近有较强的吸收, 形成两个吸收谷, 在540~560 nm附近形成了明显的反射绿峰; 在670~720 nm之间反射率增高最快, 出现“红边”; 当波长大于700 nm, 沉水植物光谱反射率显著增加, 形成了水生植物在近红外波段典型的高反射区。上述结果与杜雨春子等(2020)研究结果一致。

沉水植物群落的平均光谱反射率随覆盖度降低而下降, 在530~650 nm和710~935 nm差异明显, 但是在776~850 nm波段之间差异最为明显(图 3b)。由表 1可以看出, 在770~850 nm波段之间, 沉水植物群落光谱反射率最大值为0.186, 最小值为0.000 4, 开阔水体光谱反射率最大值为0.002 8, 最小值为0.000 1。在1 028~1 110 nm波段之间, 沉水植物群落光谱反射率最大值为0.141, 最小值为0.000 1, 开阔水体光谱反射率最大值为0.001 42, 最小值趋近于零。伴随沉水植物群落覆盖度的不断降低, 水体的面积逐渐增多, 绿光波段和近红外波段的光谱反射率也随之下降, 与Beget等(2007)的研究结果一致。当覆盖度为1%~20%时, 沉水植物群落光谱反射率与水体反射率非常接近, 但是在813 nm和1 069 nm处仍存在细微差异, 这足以表明沉水植物的光谱特征除了受到群落覆盖度影响以外, 往往还受到水体环境因素的影响(Hestir et al, 2008)。

2.2 覆盖度与沉水植物光谱反射率的相关分析

通过对不同覆盖度的沉水植物光谱特征分析, 可以发现沉水植物群落的光谱反射率随覆盖度的减少而降低。所以对2019年6、7、8、9月沉水植物群落的光谱反射率与覆盖度进行相关性分析。通过相关分析可知, 在700~900 nm波段附近, 沉水植物群落光谱反射率与覆盖度呈正相关, 在710~746 nm波段附近, 沉水植物群落光谱反射率与覆盖度呈显著正相关, 是光谱反射率对覆盖度变化最敏感的波段。皮尔逊相关系数与斯皮尔曼相关系数所对应的特征波段差异很大, 其中沉水植物群落的光谱反射率与覆盖度皮尔逊最高的相关系数位于727 nm处, R=0.844, 在882 nm处斯皮尔曼相关系数最高, ρ=0.833, 这是因为在727 nm之后沉水植物光谱反射率与覆盖度可能存在非线性关系。当沉水植物覆盖度减少时, 沉水植物的绿光波段和近红外波段的光谱反射率也随之下降, 水体环境逐渐成为主导因素, 而水体通过吸收和散射造成了光衰减(Visser et al, 2013), 这可能也是影响沉水植物群落与覆盖度相关性的因素。利用迭代算法对400~1 350 nm (5 nm间隔)波段区间的实测反射率逐一两两比值, 并与相对应的沉水植物覆盖度进行相关性计算, 发现在R_rs (860)/R_rs (710)处皮尔逊相关系数最高, R=0.866, 在R_rs (910)/R_rs (690)处斯皮尔曼相关系数最高, ρ=0.831。

2.3 沉水植物覆盖度的回归模型建立

通过相关性分析挑选以上单波段/波段比值为自变量。利用LOOCV方法确定参数a和b等, 建立沉水植物群落的单波段/波段比值的沉水植物覆盖度反演模型, 模型见表 2。

图 4为单波段/波段比覆盖度反演模型的实测覆盖度与反演覆盖度对比图, 从图 4可以看出, 波段比反演模型的稳定性要优于单波段, 这是因为波段比值可在一定程度上消除不同时间和空间上水表面光滑度和大气对电磁波的影响, 并部分地减少其他物质的干扰(潘应阳等, 2018), 而单波段却不能。从图可以看出在单波段模型中, 二次多项式反演模型的精度要高于线性反演模型的精度, 二次多项式模型的决定系数R²＜0.738, 均方根误差e＜15.03%, 平均相对误差m＜38.0%。在波段比模型中, 线性模型的精度要高于二次多项式模型的精度, 其中线性模型的决定系数R²＞0.716, 均方根误差e＜14.90%, 平均相对误差m＜35.65%。通过MRPE、RMSE和散点图可以看出大部分反演值与实测值有较好的一致性, 其中单波段模型最好的覆盖度反演精度集中在80%~100%, 波段比模型最好的覆盖度反演精度集中在60%~100%。

2.4 沉水植物覆盖度反演模型的遥感应用

本文试图利用MSI的2A级地表反射率产品数据对沉水植物覆盖度进行反演。首先, 基于波段响应函数, 将实测光谱反射率积分到Sentinel-2 MSI波段上, 见图 2和公式2; 然后, 对波段积分后的MSI单波段反射率/波段比值与实测覆盖度之间进行相关性分析, 见图 5; 由图 5a和5b可以发现, MSI第5波段(中心波段为705 nm)的沉水植物光谱反射率与覆盖度皮尔逊相关系数最高, R=0.87, 在Band 8A波段(中心波段为865 nm)斯皮尔曼相关系数最高, ρ=0.86。在波段比R_rs (Narrow NIR)/R_rs (Green)处沉水植物光谱反射率与覆盖度的相关系数最高, 其中皮尔逊相关系数为R=0.83, 斯皮尔曼相关系数为ρ=0.80。因此, 将波段积分后的R_rs (Narrow NIR)/R_rs (Green)作为自变量, 利用LOOCV方法建立MSI波段比覆盖度反演模型, 其中, 波段比线性函数模型为y=12.92x+31.2, R²=0.696, 波段比二次多项式函数模型为y=-1.75x²+ 22.32x+24.8, R²=0.725。

图 6为MSI覆盖度反演值与实测值对比图。由图 6可知, MSI覆盖度反演结果与实测值之间有良好的一致性, 其中波段比二次多项式反演结果最好, 决定系数R²为0.725, 均方根误差e为14.49%, 平均相对误差m为34.56%。通过散点图可以看出大部分反演的覆盖度位于1︰1线附近, 表明Sentinel-2 MSI波段比二次多项式反演型获得了较好的反演结果。

将MSI的波段比二次多项式模型应用于MSI地表反射率数据(2019年6月16日和8月20日), 图 7对比了波段比情况下MSI反演的覆盖度和实测值。从散点图可以看出, 反演覆盖度和实测覆盖度具有较好的一致性。其中, MSI的波段比决定系数R²为0.94。然而, 由于实测点位(1 m×1 m)与遥感数据空间尺度(20 m×20 m)不一致, 因此反演沉水植物覆盖度与实测值之间存在一定的偏差(Page et al, 2019)。

因此, 以Sentinel-2 MSI地表反射率R_rs(Narrow NIR)/R_rs(Green)作为自变量建立了沉水植物覆盖度反演模型, 模型为y= -2.15x²+25.596x+11.41, R²=0.94。利用Sentinel-2 MSI地表反射率数据反演的覆盖度散点大部分都位于1︰1线附近, 反演模型的均方根误差e为4.62%, 平均相对误差m为6.22%(图 8)。Sentinel-2 MSI地表反射率波段比二次多项式覆盖度反演模型可以用于估算研究地点的沉水植物覆盖度, 而且具有较高的精度。

因此, 利用以上Sentinel-2 MSI地表反射率数据的波段比二次多项式模型得到了2019年6月16日和8月20日乌梁素海沉水植物覆盖度分布图(图 9)。6月与8月沉水植物覆盖度差异明显, 6月沉水植物在中部的覆盖度最高, 东南部覆盖度次之, 东部覆盖度最少; 但到了8月, 沉水植物在东南部覆盖度最高, 东部覆盖度次之, 中部沉水植物覆盖度几乎为零。这是因为中部区域生长的沉水植物主要为菹草, 3月下旬到6月末为菹草生长旺期, 且叶片为长条形并浮在水面上, 7~9月为衰败期, 植株衰败死亡, 沉入湖底进入休眠期; 东部和东南部主要为龙须眼子菜, 龙须眼子菜5月中旬开始萌芽, 6月开始生长, 7~9月为生长旺期, 10月开始逐渐衰退。这种沉水植物覆盖度空间分布特征与现场船舶走航观测结果大体一致。

3 结论

1) 沉水植物群落的光谱反射率随覆盖度减少而下降, 在776~850 nm波段之间差异最为明显。当覆盖度为1%~20%时, 沉水植物群落光谱反射率与水体反射率非常接近, 但是在813和1 069 nm处仍存在细微差异。

2) 在710~746 nm波段附近, 沉水植物群落光谱反射率与覆盖度呈显著正相关, 是光谱反射率对覆盖度变化最敏感的波段。在727 nm之前沉水植物群落光谱反射率与覆盖度的线性关系较好, 在727 nm之后两者存在很好的非线性关系。

3) 利用LOOCV方法, 以单波段/波段比为自变量建立了沉水植物覆盖度的反演模型, 其中波段比模型的拟合效果要优于单波段模型, 反演精度较高, 适用于60%~100%沉水植物覆盖度反演。

4) 利用波段积分运算将实测光谱反射率扩展应用到Sentinel-2 MSI波段上, 建立了MSI单波段/波段比覆盖度反演模型。获得2019年6月16日和8月20日乌梁素海沉水植物覆盖度分布图。然而, MSI遥感反演的沉水植物覆盖度具有一定的不确定性, 这主要是由于实测点位与遥感数据空间尺度不一致造成的。因此, 在以后的研究中可以结合无人机高光谱对乌梁素海沉水植物覆盖度遥感反演模型进行验证。