浮叶植物是湖泊生态系统的重要组成成分, 它与湖泊生态系统的健康和功能密切相关^{[1, 2]}。浮叶植物不仅能有效降低水体营养盐浓度和叶绿素含量, 抑制藻类生长和繁殖, 改善水质^[3—5], 还能为小型动物提供栖息地以及作为食物来源^{[6, 7]}。但浮叶植被过度生长而不加清除或利用, 其死亡后大量的植物残体分解, 会向水体和沉积物中释放大量总氮、总磷等营养物质, 从而造成湖泊二次污染^[8]。因此, 开展浮叶植被覆盖度的长时序、快速监测, 对湖泊的健康评估和生态管理、修复都具有重要意义。

湖泊环境通常较为复杂, 用传统的人工野外调查方法开展浮叶植被的空间分布调查, 往往耗时耗力, 因而无法完成大范围浮叶植被的时空变化监测。遥感技术具有低成本、快速、实时性和历史追溯性等优势, 已成为湖泊水生植被监测的有效工具^[9]。

卫星影像二元分类方法众多, 决策树分类法是二元分类法中最常用的一种分类方法^[10]。在使用决策树分类法时, 需要通过不断尝试、更改, 并结合实测样本点数据进行精度验证, 最终确定分类精度最好的阈值。然而, 决策树分类法的最佳阈值会受水生植被生长状况、大气条件、观测几何、水面条件和水质等内在因素和外在因素的共同影响, 不具有广泛适用性。为了提高分类精度, 有学者开展了动态阈值方法的研究, 例如Dai等^[11]通过线性像元混合法, 将由5% 浮叶植被或挺水植被与95% 的纯水组成的像元作为水生植被与纯水的阈值像元, 然后计算其FAI (Floating Algae Index)值作为阈值, 对水生植被进行提取, 并取得了较高的精度; Luo等^[12]也通过线性像元混合法, 将20% 浮叶或挺水植被与80% 纯水组成混合像元作为阈值像元, 然后计算其AVI (Aquatic vegetation index)作为阈值, 从而对水生植被和水体进行分离, 也取得了很好的效果。然而, 该分类法在通过阈值提取浮叶植被时, 都是针对卫星影像像元进行定性判别, 只能将像元分类为纯水像元或者纯浮叶植被像元, 不能准确反映像元中浮叶植被的覆盖度大小。

近年来, 用遥感数据反演水生植被覆盖度的精度明显提高^[13]。Pu等^[14]利用Hyperion高光谱遥感数据对美国佛罗里达州西部海岸的沉水植物覆盖度进行了反演, 结果表明, 当覆盖度为3级和5级时, 反演精度分别达到95.9%和78.4%。Zou等^[15]在富营养化湖中, 利用实测数据研究了沉水植物光谱特征, 并指出沉水植物光谱反射率随覆盖度的降低而下降。杜雨春子等^[16]通过分析内蒙古乌梁素海沉水植被覆盖度对光谱特征的影响, 构建了基于波段比的沉水植物覆盖度反演模型, 反演结果取得了较高的精度。然而, 以往的研究所需要的野外实地数据, 都是直接通过人工调研获得, 且多把水生植被覆盖度简单分为若干等级。当前, 用遥感影像提取水生植被覆盖度的方法中仍存在两方面的不足: (1)通过现场目视或样点法来统计水生植被覆盖度的工作量大, 无法快速高效地识别浮叶植被覆盖度; (2)把水生植被覆盖度简单分为若干等级, 无法获取水生植被覆盖度连续变化的数据, 从而导致反演结果的精度下降。

本研究用无人机航拍数据来获取浮叶植被覆盖度的验证数据, 不仅大幅减少了样地调查的工作量, 而且也能降低植被样方的抽样误差^[17], 同时无人机航拍数据也适用于小面积的浮叶植被精确监测。卫星影像能够实现大尺度监测, 但精度上不如无人机调查数据, 将卫星影像和无人机影像进行结合, 建立卫星影像植被指数与无人机实测植被覆盖度之间的关系, 从而能有效提高卫星影像大尺度植被覆盖度反演的精度。

1.   材料与方法

1.1   研究区域

鄱阳湖位于江西省北部、长江中下游南岸, 是我国最大的天然吞吐型、季节性淡水湖泊^[18]。受季风气候及复杂水系结构的影响, 在鄱阳湖的干湿交替带形成了我国最大的洪泛型湖泊湿地^{[19, 20]}。鄱阳湖湿地植物种类丰富, 是长江中下游重要的湿地植物种质资源库。1983年第一次科学考察显示, 鄱阳湖湿地植被总面积2262 km², 其中浮叶植物约占23.2%, 2013年第二次科学考察显示, 浮叶植被分布面积明显下降, 主要物种菱(Trapa bispinosa Roxb)仅零星分布^{[21, 22]}。根据鄱阳湖第二次科学考察的结果, 鄱阳湖常见浮叶植物有菱、荇菜[Nymphoides peltatum (Gmel.) O.Kuntze]、芡(Euryale ferox Salisb.ex Konig et Sims)、水鳖[Hydrocharis dubia (Bl.) Backer]等^[22]。

通过观察两次鄱阳湖科考的植被分布图和鄱阳湖的卫星影像, 本研究选择南矶山自然保护区内部分子湖(北深湖、南深湖、红兴湖、石湖)、康山湖、大伍湖作为调查区域, 这些区域水生植被广泛分布, 适合布设调查样点。南矶山自然保护区位于鄱阳湖主湖区南部, 为赣江北支、中支和南支汇入鄱阳湖开放水域时, 冲积形成的赣江三角洲前沿地带^[23], 保护区内碟形湖众多, 它们在鄱阳湖丰水期高水位时融入主湖体, 在枯水期低水位时又与主湖体分离, 成为独立的水域, 形成湖中湖的独特景观。南矶山自然保护区内植物资源十分丰富, 植物资源共计115科304属443种^[24]。康山湖位于鄱阳湖南部, 面积约86.6 km², 为隔断湖汊, 水位比较稳定, 第二次鄱阳湖科考发现康山湖内有大量水生植被。大伍湖位于鄱阳湖中西部, 面积约9.2 km², 它是一个与主湖体自然连通的碟形湖。

1.2   数据与预处理

Sentinel-2卫星数据　　Sentinel-2由两颗相同的、同步运行的卫星组成, 由于两者互补, 在赤道上的重访周期仅为5d^[25]。Sentinel-2 MSI有13个光谱波段, 包括从可见光和近红外到短波红外的光谱波段^[26]。从欧空局网站(https://scihub.copernicus.eu)下载了2022年6月25日、7月31日、10月4日覆盖鄱阳湖区域的L2A级高质量卫星影像共12幅, 在SNAP(Sentinel application platform) 中完成影像重采样, 在保证研究需求的同时为减少数据的处理量, 数据重采样为20m分辨率, 提取B2、B3、B4、B8、B11、B12六个波段, 然后在ENVI 中完成数据拼接, 最终生成三景鄱阳湖全景图。

野外调查数据　　使用大疆公司的DJI Air 2S型无人机进行野外调查, 该无人机搭载了用于可见光 (RGB)成像通道的CMOS影像传感器(图 1a)。分别于2022年6月23日和28日、7月30日、10月3日在鄱阳湖开展了4次野外调查, 调查时间与所下载的卫星数据时间相同或接近。采用航点拍摄和航带拍摄两种方式对鄱阳湖浮叶植被进行调查(图 1b和1c), 航点拍摄照片用于植被类型的识别, 而航带拍摄则是让无人机对某块区域各位置连续拍照, 该照片后续完成拼接, 从而得到调查区域的全景图。

图  1  野外航拍设备及不同航拍方式布点图

a. 大疆无人机; b. 航点方式布点; c. 航带方式布点

Figure  1.  Diagram of aerial survey equipment and point layout for different aerial photography methods

a. DJI Drone; b. Waypoint Layout; c. Strip Layout

下载: 全尺寸图片 幻灯片

无人机野外调查的布点流程如下: 首先根据卫星影像的光谱特征选择调查区域, 调查区域为浮叶植被分布较多的碟形湖或隔断湖汊, 在该区域分别作两种调查点布设, 第一种是直接在图上设置散点点位, 这些点位作为后续航点点位, 第二种是在图上设置一些拐点, 将这些拐点所围成的多边形区域作为后续航带范围。将以上点位生成kml文件, 并导入rainbow无人机飞行控制软件中。

使用rainbow对无人机飞行的航点和航带路径的参数进行设置, 航点方式拍摄时飞行高度设为10 m, 所得照片清晰度高, 可准确识别植被类型。航带方式拍摄时设置航向重叠度为75%, 旁向重叠度为65%, 飞机飞行高度设为110 m, 此时影像空间分辨率约为0.037 m, 所得照片可清楚分辨水体与浮叶植被, 共获得有效航点拍摄照片200余张, 航带拍摄照片1500余张, 鄱阳湖野外调查区域位置如图 2所示。

图  2  鄱阳湖采样区分布图

Figure  2.  Distribution map of sampling areas in Poyang Lake

下载: 全尺寸图片 幻灯片

对航点照片进行目视识别, 若是浮叶植被分布区, 则将相应区域的航带照片导出, 然后在Agisoft PhotoScan Professional中, 对无人机航带照片依次进行对齐处理、建立密集点云、生成网格、生成纹理、构建数字高程模型、构建正射影像等一系列处理, 最终拼接形成8个航带的正射影像图。

1.3   浮叶植被覆盖度反演模型

浮叶植被覆盖度反演模型构建主要包括三个步骤: 卫星影像NDVI指数计算, 基于无人机影像的浮叶植被实地覆盖度计算, NDVI指数与浮叶植被覆盖度关系回归。

卫星影像NDVI指数计算　　NDVI是应用最广泛的植被指数, 它与叶面积指数LAI (Leaf Area Index)、绿色生物量及植被覆盖度等植被参数有很好的相关性^{[27, 28]}, NDVI常用于检测植被生长状态、植被覆盖度和消除部分辐射误差等。NDVI的取值范围为(–1, 1), 当地面覆盖为云、水、雪时, 它的值为负, 当地面覆盖为植被时它的值为正, 数值越大则表示植被覆盖度越大。在ArcGIS中使用Float函数计算卫星影像的NDVI指数, NDVI的计算公式 ^[28]:

式中, NIR表示近红外波段的反射率, R表示红色波段的反射率。

浮叶植被实地覆盖度计算　　本文基于无人机正射影像计算浮叶植被的实地覆盖度。无人机影像为RGB图像, 仅包含红光、绿光、蓝光三种波段, 因此使用可见光植被指数可完成植被提取。可见光波段差异植被指数(Visible Light Band Difference Vegetation Index, VDVI)对植被的反射和吸收波段十分敏感, 该植被指数演变自NDVI指数, 它用绿光波段(G)取代NDVI指数的近红外波段, 同时用红光波段(R)和蓝光波段(B)的和取代NDVI指数的红光波段, 并将绿光波段(G)乘以2使其在数量上与R+B相当^{[29, 30]}。植被与非植被的VDVI指数高度分离, 因此VDVI指数被广泛用于无人机照片的植被提取^[31]。

在ArcGIS中使用Float函数计算无人机影像的VDVI值, 然后利用阈值法区分水体和植被, 分别将水体和植被像元赋值为0和1。阈值为提取范围与原图浮叶植被区域吻合度最好时的VDVI值(图 3a), 它通过目视判别法反复尝试获得。VDVI指数的计算公式 ^[29]:

图  3  浮叶植被提取效果(a)和样方布置图(b)

Figure  3.  Extraction results of floating leaf vegetation (a) and plot layout diagram (b)

下载: 全尺寸图片 幻灯片

式中, R、G、B分别表示红光、绿光、蓝光波段的亮度。

将无人机影像和Sentinel-2卫星影像进行空间叠加, 针对二者的重叠区域, 通过Sentinel-2影像生成一系列网格, 将这些网格作为样方范围叠加到无人机影像上(图 3b), 每个20 m×20 m网格覆盖约30万个无人机影像像元, 分别统计每一个网格中VDVI值为1的比重即为该样方的植被覆盖度, 它代表相同位置的哨兵卫星影像像元内的实测浮叶植被覆盖度。

卫星影像NDVI指数与浮叶植被覆盖度的回归分析　　使用ArcGIS中的多值提取至点工具, 针对哨兵数据和无人机数据的重叠区域, 提取每个样方(20 m×20 m的格网)的NDVI值和浮叶植被覆盖度数据, 然后在R软件中建立二者的回归关系。四次无人机野外调查生成的8幅影像共覆盖2041个样方。 由于10月3日北深湖部分浮叶植被已经枯萎发黄甚至发黑, 而南深湖和红兴湖的水深都在0.2 m以下, 光谱特征受底质的影响很大, 由这3幅无人机影像形成的样方数据放弃使用。最后从其余5幅无人机影像获得有效样方1200余个, 其中80%的样方用于模型回归, 其余20%的样方用于精度评估。

1.4   模型精度评估

通过以下2种方法来评估回归模型的精度。

(1) 用RMSE和MRE 评价模型精度。用于精度评估的240个样方来自南深湖、南湖、康山湖和大伍湖, 它们未参与NDVI与实地覆盖度的关系回归。构建这240个样方的实地覆盖度与反演覆盖度的散点图, 其均方根误差RMSE和平均相对误差MRE的计算公式:

式中, X_ia和X_ib分别表示实地覆盖度和反演覆盖度, n为参与精度评估的样方总数。

(2)与像元二分法的反演精度比较。以6月23日所拍摄的有大面积浮叶植被分布的大伍湖为实验区。分别使用像元二分法和本文方法对该区域的浮叶植被的覆盖度进行计算, 并使用均方根误差和平均相对误差对两种方法的分类精度进行比较。为了尽可能提高像元二分法的分类精度, 本文采用动态阈值完成像元二分法, 而且取像元二分法分类结果精度最高的数据与本文方法的反演精度进行比较。

2.   结果

2.1   NDVI指数与浮叶植被实地覆盖度的关系

NDVI指数与浮叶植被实地覆盖度的散点图见图 4。二项式具有很好的拟合效果, 其决定系数R²达到0.9, 二者的回归关系式如下:

图  4  卫星影像NDVI指数与浮叶植被覆盖度的散点图

Figure  4.  Scatter plot of satellite imagery’s NDVI index and floating leaf vegetation coverage

下载: 全尺寸图片 幻灯片

式中, X和Y分别为卫星影像的NDVI值和浮叶植被实地覆盖度。

2.2   模型的精度

采用RMSE和MRE的评估结果　　由图 5可知, 浮叶植被覆盖度的反演结果与无人机数据计算值(实测值)之间高度正相关, 其中决定系数R²达到了0.90, 均方根误差RMSE为5.75%, 平均相对误差MRE为9.0%。通过散点图可以看出, 绝大部分点位于1﹕1线附近, 表明该回归模型获得了很高的反演精度。

图  5  浮叶植被覆盖度实测值与反演值对比图

Figure  5.  Comparison between measured and inverted values of floating leaf vegetation coverage

下载: 全尺寸图片 幻灯片

与像元二分法精度的比较结果　　如图 6所示, 本文方法不仅能通过卫星影像反演浮叶植被的覆盖度, 而且能反映植被覆盖度在空间上的连续变化(图 6B), 它是一个定量的反演模型。传统的像元二分法则是基于卫星影像判断像元是否为浮叶植被, 在像元层面是一个定性判断, 但用来判断研究区的浮叶植被分布范围时, 又可看做一个定量模型, 如图 6C。

图  6  航拍影像及两种方法反演植被覆盖度的效果图

A. 大伍湖正射影像; B. 本文方法反演结果; C. 决策树分类法提取结果

Figure  6.  Aerial imagery and effectiveness of two methods for inverting vegetation coverage

A. Orthophoto of Dawu Lake; B. Inversion results by method in this study; C. Inversion results by decision tree classification method

下载: 全尺寸图片 幻灯片

基于2022年6月25日的Sentinel-2卫星数据, 计算得到大伍湖区域NDVI的取值范围在–0.123—0.368, 因此把它作为二元分类法的NDVI阈值变化区间, 图 7和图 8分别表示使用像元二分法时, 均方根误差RMSE和平均相对误差MRE随着阈值(NDVI)变化的动态, 可以看出, 均方根误差RMSE和平均相对误差MRE均随着分类阈值的增加, 呈现先减小后增大的趋势。像元二分法反演结果的均方根误差RMSE为20.25%—65.91%, 平均相对误差MRE为53.68%—133275.70%, 当NDVI阈值为0.081时, 均方根误差RMSE取得最小值20.25%, 平均相对误差MRE取得最小值53.68%, 而本文方法反演结果与实地航拍数据相比的RMSE和MRE分别为5.52%和28.33%。这说明无论如何调整阈值, 像元二分法反演结果的误差始终大于本文方法的误差。

图  7  像元二分法的均方根误差RMSE随阈值的变化图

Figure  7.  Graph of pixel-based dichotomous method’s RMSE and its variation with threshold values

下载: 全尺寸图片 幻灯片

图  8  像元二分法的平均相对误差MRE随阈值的变化图

Figure  8.  Graph of pixel-based dichotomous method’s MRE and its variation with threshold values

下载: 全尺寸图片 幻灯片

3.   讨论

3.1   文章方法较像元二分法的优势分析

像元二分法反演植被盖度仅判断像元是否为浮叶植被, 无法体现该像元浮叶植被覆盖度的准确信息。在自然水体中, 往往存在大面积较稀疏的浮叶植被, 因此有大量由浮叶植被和水组成的混合像元, 像元二分法(如决策树分类法)仅仅根据所设定的阈值而直接将其归类为纯植被或纯水像元, 无法反映像元内浮叶植被覆盖度的真实情况, 因而导致浮叶植被分布面积与实地覆盖度之间的较大差异。本文构建的方法可根据卫星影像像元的NDVI值反演像元的真实覆盖度, 能捕捉单个像元植被覆盖度的变化信息, 适用于不同植被覆盖度的情况, 因此具有更高的精度。

通过对鄱阳湖部分碟形湖和隔断湖汊的调查发现, 在浮叶植被群落边缘过渡带, 通常会形成较稀疏的浮叶植被, 群落中间则往往形成大面积更高密度的浮叶植被群落。因此, 在使用本文所述方法和像元二分法进行浮叶植被面积提取时, 边缘过渡带会有很大差异, 而在植被覆盖度较高的区域, 像元二分法与本文所述方法相比, 结果的差异较小。

3.2   本回归模型的适用性及变异因素

本文所构建的回归模型Y=277.4X²+86.572X+7.8628, 在X分别取值–0.16和0.44时, Y分别得最小值1.1和最大值100。因此, 该回归模型自变量NDVI的适用范围为–0.16—0.44。在实际应用当中, 当卫星影像像元NDVI值低于–0.16时, 可直接认定该像元为纯水像元, 当NDVI值大于0.44时, 则直接认定该像元为纯植被像元或者调整植被100%覆盖时的最大NDVI阈值。

当水较浅时, 像元的光谱会受到底质的影响^[32], 常常会导致像元NDVI值增大, 如浅水湖泊的湖底为裸土或沉水植被, 都会导致该区域的NDVI值大于深水区纯水像元的NDVI值^{[33, 34]}。因此, 在使用该回归模型对浮叶植被覆盖度进行反演时, 浅水区的反演结果会比实际值偏大。在实际应用中, 为了减小反演结果的误差, 可把浅水和深水区分开进行计算, 针对浅水区用同样的方法再构建一个回归方程, 以消除背景因素的影响。而自然湖泊的浅水区域一般只是在湖边, 占整湖面积的比例较小, 即使浅水区域反演结果存在误差, 对模型总体精度的影响有限。

拟合方程的散点图(图 4)表明, 在特定NDVI值时, 浮叶植物的盖度仍有一些差异, 表明可能有其他因素影响NDVI值与植被覆盖度的关系。其中因素之一可能与浮叶植物所处的生长阶段有关, 植物在生长的不同阶段, 其叶绿素含量存在一定的变化, 而NDVI值和叶绿素含量存在正相关关系^[35—37], 因此造成同一NDVI指数值对应不同覆盖度的情况。另外, 查阅10月3日在南深湖、北深湖拍摄的照片发现, 已经有部分浮叶植物发黄发黑, 濒临死亡, 这部分浮叶植物的叶绿素含量由于衰老而减少, 其NDVI值比生长旺期的浮叶植被偏小, 因而根据10月3日的调查数据所做的散点图相比于6—7月的散点图整体向左偏移, 表明在植物生长末期, 直接用本文的拟合方程计算的覆盖度会偏小。在对处于生长末期的浮叶植被进行覆盖度反演时, 可先用同样的方法构建一个适用于该阶段浮叶植被的回归方程, 再进行浮叶植被覆盖度的反演。

本研究还发现, NDVI值随着植被覆盖度的增加而增大, 但后期的增长速度越来越慢, 这是因为NDVI指数在植被高覆盖区容易饱和, 从而导致植被覆盖度增大时, 植被指数无法同步增长^[34]。

4.   结论

(1)建立了鄱阳湖菱植物群落NDVI指数与植被覆盖度的回归模型, 回归模型的决定系数R²达到了0.9, 均方根误差RMSE为5.75%, 平均相对误差MRE为9.0%, 可用于鄱阳湖浮叶植物菱群落面积的精确计算。(2)相比于常用的像元二分法, 本文构建的回归模型能更准确地提取浮叶植被的面积, 特别是在浮叶植被稀疏的区域。(3)文章基于无人机航拍数据和遥感影像NDVI指数构建回归模型的方法, 可以为水生植物群落快速监测和定量反演提供借鉴。(4)本模型没经过其他浮叶植物群落的实测数据验证, 也没有应用于Sentinel-2之外的遥感数据源。但是, 仍可利用文章的模型构建方法, 构建适用于不同水生植被、不同遥感数据源和不同水环境背景下的水生植被快速监测模型。