基于穷举法的鄱阳湖叶绿素a浓度高光谱反演模型与应用研究<br/>——以GF-5卫星AHSI数据为例

基于穷举法的鄱阳湖叶绿素a浓度高光谱反演模型与应用研究
——以GF-5卫星AHSI数据为例

2020-06-18韦安娜田礼乔陈晓玲余永明

华中师范大学学报（自然科学版） 2020年3期

韦安娜，田礼乔，陈晓玲*,余永明

(1.武汉大学测绘遥感信息工程国家重点实验室，武汉 430079;2.北京城建勘测设计研究院有限责任公司，北京 100101)

湖泊是地球水圈的重要组成单元，也是人类赖以生存及社会可持续发展的重要依托[1].随着社会经济的高速发展，在自然和人为等因素的共同作用下，我国湖泊水质状况持续恶化[2]，尤其是水体富营养化现已成为我国湖泊水环境所面临的主要问题[3].叶绿素a(chlorophyll-a,Chl-a)是衡量水体富营养化程度的一个重要指标，同时也是Ⅱ类水体水质监测和水色遥感反演的主要参数之一[4].因此，进行叶绿素a浓度监测对防治湖泊水体富营养化、保证流域生态可持续开发及利用具有重要的意义.

传统的水质监测主要依靠人工调查，可精确测定水体局部叶绿素a浓度，但成本高、耗时长，且相对于大范围水体，其空间代表性不足[5].遥感技术以其大范围、多时相、高动态及低成本等优势为水质监测拓展了新的途径[6]，逐渐成为大面积湖泊水质变化动态监测的重要手段.国内外学者开展了大量的针对内陆水体的叶绿素a浓度遥感反演研究工作，如Feng等[7]根据MERIS数据构建了针对鄱阳湖的叶绿素a浓度反演算法，Sakuno等[8]开展了基于Landsat-8 OLI数据的日本宍道湖叶绿素a浓度反演研究，朱云芳等[9]提出了针对高分一号WFV数据的太湖叶绿素a浓度反演算法.然而，湖泊、河流等内陆水体光学特性复杂，常用的宽波段多光谱数据无法捕捉其精细的光谱特征，在水质参数反演中有一定的局限性.相较于多光谱数据，高光谱遥感数据具有纳米级的光谱分辨率[10]，包含更多的光谱信息，有望提高水质参数反演的精度[11].目前，高光谱卫星数据逐渐应用到内陆水体叶绿素a浓度反演研究中，如杨煜等[12]利用环境一号卫星高光谱数据进行了巢湖的叶绿素a浓度反演，陈军等[13]基于Hyperion影像建立了太湖叶绿素a四波段半分析模型，但针对鄱阳湖的高光谱叶绿素a浓度反演研究较少.

我国GF-5卫星于2018年5月9日发射升空，是世界上首颗全谱段对大气/陆地综合观测的高光谱卫星，其上搭载的可见短波红外高光谱相机(Advanced HyperSpectral Imager,AHSI)是我国高分系列卫星中第一个高光谱传感器，其光谱分辨率高，谱段数多，精度高，信息量大，可充分捕捉目标地物光谱细节特征，有利于获取内陆水体水质参数的细微变化，在内陆水体水色遥感研究中有巨大潜力[14].因此，针对目前多光谱卫星传感器在光学特性复杂水体的叶绿素a反演精度难以进一步提高的现状[15]，开展GF-5卫星AHSI数据的内陆水体叶绿素a浓度反演研究具有重要的科学意义与现实价值.

本文以鄱阳湖为研究区域，结合2009年10月、2011年7月和2016年10月现场实测数据，利用穷举法建立了基于GF-5卫星AHSI高光谱数据的叶绿素a浓度反演模型，并将其反演结果与Sentinel-2A卫星多光谱成像仪(MultiSpectral Instrument,MSI)的输出进行了对比分析.

1 研究区域与数据

1.1 研究区概况

鄱阳湖位于江西省北部，是我国的第一大淡水湖，平均水深8.4 m.作为吞吐型通江湖泊，鄱阳湖水体边界季节性变化剧烈，丰水期(4～9月)平均面积为5 100 km2，枯水期(10月～次年3月)平均面积为510 km2，形成了“丰水一片，枯水一线”的独特景象[16].鄱阳湖光学特性异常复杂，其水位、湖泊面积、光学特性等的高动态变化为研究不同性质内陆湖泊水体定量遥感难题提供了天然实验场，对其他湖泊的研究也具有一定的借鉴意义.近年来，鄱阳湖水环境形势越来越严峻，甚至开始出现水体富营养化的趋势[17]，研究其叶绿素a浓度分布状况对于开展水体富营养化防治、缓解其水环境压力具有重要的现实意义.图1所示为研究区域及实测站点的分布图.

1.2 数据获取与处理

1.2.1 现场实测数据鄱阳湖现场光谱测量及水样数据采集分别于2009年10月15～23日，2011年7月17～26日，2016年10月1～14日进行.现场光谱观测时间控制在北京时间上午10时～下午2时之间.水体光谱数据由美国Spectra Vista公司生产的便携式地物光谱仪SVC HR-1024基于水面以上法测量得到，观测几何及数据处理等严格参照美国国家航空航天局(National Aeronautics and Space Administration,NASA)海洋光学测量规范的相关要求[18].经过数据筛选剩余样本共63组，光谱曲线如图2所示.

图1 研究区域及实测站点分布图Fig.1 Study area and location of sampling stations

图2 实测遥感反射率光谱曲线Fig.2 In-situ spectrum of remote sensing reflectance

从图2可以看到，鄱阳湖不同区域的遥感反射率在可见光和近红外波段范围内的变化较为显著，并有如下明显特征：400～500 nm范围内的水体反射率较低，随着波长增加呈现缓慢上升趋势[19]；遥感反射率在580、700、810 nm附近出现明显的峰值，而在680、760 nm附近出现谷值.其中，580 nm附近处的反射峰是藻类叶绿素的弱吸收和细胞的散射作用所致[20]；680 nm附近是叶绿素的吸收峰，当藻类密度较高时，水体反射率出现谷值[21]；700 nm附近，悬浮物和叶绿素的后向散射作用增加，使水体反射率增加，出现“红移”现象[22]；760 nm附近的反射谷可能是近红外波段内叶绿素a的强吸收所致；而810 nm附近的反射峰则是由于Ⅱ类水体高悬浮物的后向散射作用形成的.

此外，本研究所利用的实测叶绿素a浓度数据样本共63组(40个用于建模，23个用于验证)，浓度范围为1.33～11.37 μg·L-1，平均浓度为5.93 μg·L-1.水样采集与水体光谱数据采集同步进行，且叶绿素a浓度实验室处理及分析步骤严格按照国家环境保护标准《水质叶绿素a的测定分光光度法》(HJ 897-2017)执行[23].

1.2.2 GF-5卫星AHSI高光谱影像 GF-5卫星AHSI数据在0.4～2.5 μm的光谱范围内细分为330个谱段，在0.4～1.0 μm(可见光近红外波段)范围内谱段数为150，光谱分辨率为5 nm；在1.0～2.5 μm(短波红外波段)范围内谱段数为180，光谱分辨率为10 nm.同时，AHSI数据空间分辨率为30 m，幅宽为60 km.本研究使用的鄱阳湖区域GF-5卫星的AHSI数据来源为中国资源卫星应用中心数据分发网站(http://data.cresda.com:90/#/home)，成像时间为北京时间2018年10月7日13时28分.

AHSI数据预处理主要包括水体范围提取、辐射定标、几何校正和大气校正.本文使用图形处理软件ENVI5.3(https://www.harrisgeospatial.com/Software-Technology/ENVI)进行人工数字化提取鄱阳湖水域，然后利用给定的辐射定标系数对该影像进行辐射定标，接着以一幅经几何精纠正的ETM+影像作为参考图对AHSI影像进行几何校正，再通过FLAASH大气校正模块实现对该高光谱遥感影像的大气校正，其主要参数设置如下：大气模式设为中纬度夏季类型；气溶胶模式设为城市类型；可见度设为40 km；气溶胶反演选用2-Band (K-T)算法；传感器高度为705 km；海拔高度为0.02 km.最后，通过水—气界面校正将结果转换成遥感反射率.

2 叶绿素a浓度反演模型的构建

2.1 模型构建

Gordon等[24]基于水体辐射传输理论提出了一个被广泛认可的生物光学模型，描述了水面以上遥感反射率与水体固有光学特性的关系，据此可建立三波段组合的叶绿素a浓度反演模型.国内外众多学者[25-27]利用三波段模型在不同内陆水域进行了叶绿素a浓度反演研究，并证明了该模型具有较高的精度，尤其是在光学特性复杂的内陆浑浊水体有一定的适用性.本研究沿用三波段模型进行叶绿素a浓度的反演.

三波段模型的具体形式如下：

(1)

其中，Rrs(λi)代表波长λi处的水面之上遥感反射率.叶绿素a浓度与Rrs(λ1)、Rrs(λ2)、Rrs(λ3)相关，选取适当的波长λ1、λ2、λ3能有效地减少水中悬浮颗粒物、CDOM等组分对叶绿素a浓度反演结果的干扰.λ1、λ2、λ3分别具有以下特征：

1)λ1为叶绿素a吸收系数αChl-a的最大敏感波长，范围在660～690 nm之间[28]；

2)λ2应同时满足叶绿素a吸收系数αChl-a(λ2)≪αChl-a(λ1)，CDOM吸收系数αCDOM(λ2)≈αCDOM(λ1)，非藻类颗粒物吸收系数αNAP(λ2)≈αNAP(λ1)这三个条件；

3)λ3用于消除水体总反向散射系数bbt的影响，该波长处水体各组分(除纯水外)的吸收系数应接近于零，即αt(λ3)=αw(λ3)且αt(λ3)≫bbt(λ3)，因此应大于710 nm[29].其中，αt为水体总吸收系数，αw为纯水吸收系数.

较之多光谱遥感数据，GF-5卫星AHSI高光谱数据可为叶绿素a浓度三波段反演模型提供不同的波段选择与组合.为了选择最佳的波段组合，可以根据经验知识[30]进行迭代得到最佳波段，也可利用穷举法得到最优波段组合.由于采用迭代法选取的波段在特征波段范围内可能只符合局部最优，而穷举法能保证全局最优解，同时可反向验证所选取的三个波段是否符合其波段特征条件，故本文采用穷举法对λ1、λ2和λ3这三个波段进行最优选择.图3为不同三波段组合的叶绿素a模型的建模精度(本文以决定系数R2为标准).

图3 各波段组合的叶绿素a浓度反演模型决定系数R2三维分布图Fig.3 Three-dimensional distribution of the determination coefficient (R2) of the three-band Chl-a inversion model

在图3中，不同颜色表示R2值的大小.红色标注区域是R2较大的波段范围，最优波段应在此范围内选取；黄色标注区域是R2较小的波段范围；而绿色标注区域则表示不参与计算的重叠波段.由此，最终根据建模精度确定最优的λ1、λ2和λ3分别为681.3 nm、709.9 nm和829.6 nm.结合前文对鄱阳湖水体光谱特征的分析发现，三波段选取结果符合上述特征波段要求，表明该模型具有一定的物理意义.

2.2 模型验证

本文模型采用均方根误差(RMSE)、平均绝对误差(MAE)和平均相对误差(MRE)进行评价，三者公式如下：

(2)

(3)

(4)

其中，Cchl-a*和Cchl-a分别表示叶绿素a浓度实测值和反演值，N为验证集样本点数量.

(5)

叶绿素a浓度的建模精度和验证精度如图4所示.此算法的建模及验证的R2、RMSE、MAE和MRE分别为0.82、1.468 μg·L-1、1.094 μg·L-1、19.948%和0.74、2.702 μg·L-1、2.173 μg·L-1、44.015%.

图4 GF-5卫星AHSI叶绿素a浓度反演模型与验证示意图Fig.4 Modeling and validation results of the Chl-a model for GF-5 AHSI

3 结果与分析

为了评估和交叉验证本文的鄱阳湖叶绿素a浓度反演模型与结果，本研究收集了2018年10月7日同步的Sentinel-2A卫星MSI影像，并分析了两者叶绿素a浓度反演结果的一致性.首先利用ENVI5.3对MSI和AHSI影像进行手动配准，精度控制在一个像素内，然后利用Acolite软件(https://odnature.naturalsciences.be/remsem/acoliteforum/)对MSI影像进行大气校正[31-32]，参数设置为Acolite默认参数，主要包括：气溶胶校正设置为暗光谱拟合(DSF)算法[33]；气溶胶校正波段选择1 600 nm和2 200 nm；550 nm处最小气溶胶光学厚度为0.01.最后，结合李亭亭等[34]基于MSI数据建立的鄱阳湖叶绿素a反演模型进行MSI的叶绿素a浓度反演，即在北湖区采用三波段模型(R2=0.65,MRE=38.53%)，在南湖区采用差值模型(R2=0.63,MRE=39.87%).

AHSI三波段组合模型的叶绿素a浓度反演结果如图5(a)所示；MSI的叶绿素a浓度分布如图5(b)所示；两者偏差与统计结果见图5(c)和图6.

显而易见，AHSI与MSI反演的叶绿素a浓度在空间分布上是较为一致的.总的来说，叶绿素a浓度在鄱阳湖北部和中部浓度较高，在鄱阳湖南部浓度较低.AHSI与MSI的叶绿素a浓度反演结果在鄱阳湖北部偏差较大，但仍不超过20%；在鄱阳湖其它水域，两者偏差更小，不超过10%.从图6的相对偏差统计结果可以看出，在超过80%水域AHSI与MSI反演的叶绿素a浓度相对偏差小于8%，表明AHSI与MSI的叶绿素a浓度反演结果具有高度的一致性.

图5 GF-5卫星AHSI叶绿素a浓度分布(a)、Sentinel-2A卫星MSI叶绿素a浓度分布(b)和两者相对偏差(c)Fig.5 Distributions of GF-5 AHSI-derived Chl-a(a),Sentinel-2A MSI-derived Chl-a(b),and relative bias between them(c)