Sentinel-2/MSI深度学习超分辨率重建及河湖水质遥感反演

王世瑞,沈芳,魏小岛

(1.华东师范大学 河口海岸学国家重点实验室,上海 200062;2.上海勘测设计研究院有限公司,上海 200050;3.中国长江三峡集团有限公司 长江生态环境工程研究中心,北京 100038)

0 引言

在人类活动和气候变化的持续影响下,许多内陆水体正面临水质恶化的问题[1]。水质监测对控制水污染的延伸和维持水环境健康至关重要,是水环境管理的主要方法[2]。与传统的水质监测方法相比,卫星遥感技术不但具有周期性观测、效率高和成本低等优点,能够满足大范围水质同步监测需求,而且可以揭示常规措施难以发现的污染物迁移分布特征[3]。

Sentinel-2/MSI卫星数据由于获取免费且空间分辨率高,被广泛应用于遥感水质监测[4-5]。但该卫星多光谱波段的空间分辨率并不一致,包括10 m、20 m和60 m。通常,使用全波段数据时会将Sentinel-2/MSI的空间分辨率重采样至10 m。但传统的空间域重采样算法,如最近邻、双线性或双三次插值(bicubic)等,未利用光谱域信息,难以有效降低波段中混合像元的数量[6]。近年来,针对Sentinel-2影像结合光谱域和空间域的插值算法逐渐得到发展。传统算法方面,SupReMe(super-resolution for multispectral multiresolution estimation)算法通过对Sentinel-2高分辨率波段纹理信息进行编码并传递给低分辨率波段作为规则器,同时在光谱维上求解一个凸优化问题,以实现影像的超分辨重建[7]。深度学习算法方面,Lanaras等[8]基于波段融合思想,利用残差神经网络学习低/高分辨率图像间的纹理映射关系,融合空间及光谱信息对全局进行回归训练,实现Sentinel-2/MSI L1C大气顶层反射率影像的超分辨重建,并将此算法命名为Dsen2。Dsen2不仅能降低空间异质性,而且能出色的保留光谱特征。然而Dsen2神经网络模型结构单一,仍有一定改进空间,比如通过注意力机制对模型中的重要特征赋予权值,能够使信息提取过程更为有效,起到优化模型的作用[9]。

水质遥感监测通过探究水体的反射率、散射系数和吸收系数等光学特性与对应实测水质参数之间的关系,建立水质遥感反演模型,从而获得整个水域水质参数的空间分布情况和浓度变化趋势[10]。反映水质的参数主要有叶绿素、总磷、总氮、浊度、氨氮和高锰酸盐指数(CODMn)等。这些参数与水生态环境息息相关,其值不在正常范围时,水生动植物的生活环境会受到影响[11]。估算水质参数含量的遥感方法较多,主要有经验方法[12]、半解析方法[13]、机器学习方法和深度学习方法。近年来,由于机器学习和深度学习具有更高的预测精度,已被广泛应用于水质遥感监测领域[14-15],前者集成度较高,难以有所创新,后者更复杂,应用潜力更高,有望实现更智能的水质监测。

本研究在Dsen2深度学习模型的基础上,结合卷积注意模块(convolutional block attention module,CBAM)进行迁移学习,对低分辨率波段(>10 m)进行超分辨率重建;研究超分辨率算法在水质参数定量反演中的应用成效,分析影像重建前后对水质参数的敏感性;构建多种回归模型,对比影像重建前后水质参数的反演精度,最后选取最优模型进行制图。本研究成果可为精细化水质遥感监测提供方法借鉴与数据支撑。

1 研究区与数据

1.1 研究区域

上海市紧邻长江口,平均海拔2.19 m,属于亚热带季风气候,四季分明;水系众多,河网密布,存在约14 000条河流;最大内陆河为黄浦江,长82.5 km,发源于淀山湖,上游以农业区为主,下游则以城镇区为主,最后汇入长江。上海市除众多河流外,且有许多湖泊、水库及水塘等。

1.2 数据及预处理

1)水质数据。内陆水体的水质参数以总氮、总磷、氨氮、CODMn和浊度为例,其数据来源于国家科技基础条件平台—国家地球系统科学数据中心(https://www.geodata.cn/[2022-03-01])及青悦数据(https://www.epmap.org/[2022-03-01])。监测站点主要分布于上海黄埔江(表1),监测时间从2020年11月至2021年11月。去除离散异常值后,分别获得330个总氮、238个总磷、252个氨氮、274个CODMn和351个浊度数据。

2)遥感数据及预处理。Sentinel-2/MSI是极轨双星,有A、B两颗卫星,联合重访周期为5 d,搭载13个光谱波段的多光谱成像仪。根据水质数据采集的时间和站点匹配卫星数据,水质采样与Sentinel-2/MSI卫星过境时间差应小于4 h,据此共匹配并下载103景Sentinel-2/MSI的L1C产品(https://scihub.copernicus.eu/dhus/#/home[2022-03-01])。L1C产品需大气校正以消除大气作用,获取地表真实的反射率信息。众多学者使用Sentinel-2影像时,对比了不同大气校正算法如6S、Sen2Cor、Acolite,考虑水体、建筑物、植被等区域因素,Sen2cor均表现出不错的精度[16-17]。因此本文使用SNAP软件(http://step.esa.int/main/download/snap-download/)中的Sen2Cor算法对L1C影像进行大气校正,将大气表观反射率转换成地表反射率,以输出L2A产品。

假设实测水质参数对应影像周围50 m×50 m的水体反射率具有函数映射关系,因此建立以监测站点为中心,边长为50 m的正方形,并提取每个实测水质参数对应影像的反射率值。使用改进的归一化水体指数[18]提取水体,相比于归一化水体指数,该指数可消除建筑物的影响。

在提取水体后,由于原始天顶反射率数据的问题,大气校正后,影像地表反射率数据仍存在异常值,本研究采用表观可见波长(apparent visible wavelength,AVW)[19]指数来去除。AVW指数利用加权调和平均的方法对可见光光谱反射率进行计算,其值与每个波段反射率的权重有关,能反映光谱形状,在一定程度上也可反映水色类型。计算AVW值,其中正常水体反射率计算的AVW值占比更大,取其众数位置上的反射率以去除异常值。AVW的计算如式(1)所示。

(1)

式中:R为各波段反射率值;λi为第i个波段的波长。

2 研究方法

2.1 基于迁移学习的超分辨率重建模型

本研究采用Dsen2算法提升Sentinel-2/MSI空间分辨率,探究超分辨技术对水质遥感监测的影响程度。该方法假设具有不同地面采样距离(ground sampling distance,GSD)的波段之间的量度关系在相关标度范围内是自相似的[20],即尺度不变性,意味着分辨率20 m→10 m和40 m→20 m、60 m→10 m和360 m→60 m之间的映射基本一致,在一定比例范围内,该方法可将Sentinel-2/MSI所有低分辨率波段提升至10 m。本研究假设Dsen2算法逻辑可应用于Sentinel-2 L2A产品。因此,采用基于模型的迁移学习方法赋予预训练权值,提升训练效率。在此基础上,引入CBAM进行全局注意,CBAM模块会将输入特征依次进行一维压缩和二维压缩,分别对通道和空间推算注意力图,然后将注意力图与输入特征图相乘,传输重要特征信息,以进行自适应特征选择[21],从而优化深度学习模型。新模型命名为哨兵2残差卷积注意(Sentinel-2/MSI resblock convolutional attention,S2RESCA)模型。

为了简化符号,以s倍超分辨率重建模型为例,首先对原始影像进行高斯模糊处理,模拟其调制传递函数,通过对s×s窗口取平均来下采样,获取低分辨率真实影像。当s=6时,将具有相同分辨率的波段分为3组:A=(B2,B3,B4,B8)(GSD=60 m),B=(B5,B6,B7,B8A,B11,B12)(GSD=120 m),C=(B1,B9)(GSD=360 m)。在进行训练时,C上采样的高频信息来源于A、B、C耦合的多尺度特征。假设A中高分辨率波段的空间维度为W×H,可表示为式(2)。

(2)

式中:Ca∈R2×W×H,为生成的空间分辨率为60 m的超分辨图像。

如图1所示,模型结构的主要组成部分是卷积层(Conv)、非线性激活函数(ReLU)、残差神经网络模块(Resblock)以及CBAM。B、C分别为通过Bicubic上采样得到Yb、Yc。模型训练数据为[A,Yb,Yc](GSD=60 m),其空间分辨率一致;令原始影像波段为标签数据Cx(GSD=60 m)。

图1 深度学习超分辨率重建算法过程

将Yc与经过网络提取的高频特征信息通过跳越连接进行融合,即可得到超分辨率重建影像,如式(3)所示。

Ca=Yc+R(Cx,Yc)

(3)

式中:R(Cx,Yc)为低/高分辨率之间的非线性映射,本质上是高分辨率影像具有的高频信息,空间维度为2×W×H,模型输出Ca(GSD=60 m)即为“超分辨率重建影像”,其与C(GSD=360 m)对应,可与Cx作比较,以评估模型的性能。当s=2时,B上采样的信息来自A和B,其原理及配置均与6倍神经网络一致。

采用均方根误差(RMSE)以及信号与重构误差比(signal to reconstruction error ratio,SRE)评估超分辨率重建模型的效果。RMSE越小,SRE越大,说明模型重建结果越接近真实图像,模型精度越高。

2.2 水质参数反演模型

1) 统计回归模型。利用指数、幂指数、线性及多项式等方法(表2)对水质参数(总氮、总磷、浊度、氨氮及高锰酸盐指数)建立统计回归模型,其自变量为波段或波段组合(A-B、A+B、A/B、(A-B)/(A+B),A与B表示不同的随机波段)。

表2 统计回归模型

2)深度神经网络模型。具有光学敏感信号的水色参数有悬浮物浓度、浮游植物色素浓度、有色溶解有机物(CDOM)等。CODMn的光学信号呈现弱敏感特征,总氮、总磷及氨氮属于非光敏信号水质参数,尽管其与水色参数可能存在共变关系,但是常规的单、多元线性回归模型难以构建卫星数据与水质参数数据之间的量化关系。相较而言,机器学习或深度学习算法能更好地表达卫星数据与水质参数之间潜在的复杂关系,以牺牲解释性为代价,构建精度更高的反演模型,具有更好的预测性。Peterson等[22]研究发现逐级衰减深度神经网络模型(pDNN)在水质定量反演中具有相当不错的稳健性,其模型增加了第二层神经元数量以扩大特征搜索范围,并开始逐级递减。为此,本研究通过Tensorflow和Keras构建了一个具有多个隐藏层的pDNN模型,相比于浅层神经网络,可提取更深层次的特征。该模型包含6个全连接层,神经元数量依次为[12,48,32,16,8,1],输入层的数据为与水质参数相关性较高的光谱特征,输出层为水质参数浓度。调用线性激活函数(ReLU)来传递权重,优化器为均方根传递(RMSprop)算法,损失函数为均方根误差,学习率Lr的初始值为10-3,设置为自适应调整,在训练批次中损失10次不变时,按式(4)自动降低学习率,防止陷入局部最优。模型通过绝对系数(R2)、RMSE和对称平均绝对百分比误差(symmetric mean absolute percentage error,sMAPE)进行评价。

Lrnew=Lr×factor

(4)

式中:Lrnew为新的学习率;factor为衰减系数,一般取0.1。

具体研究方法主要包括3个部分:超分辨重建方法对比、水质遥感反演回归模型对比、水质参数浓度预测和制图分析,流程如图2所示。

图2 基于Sentinel-2/MSI超分辨率重建及水质遥感反演研究的流程

3 结果与分析

3.1 超分辨重建方法与质量评价

1)超分辨重建方法的定量评价。当s=6时,利用10景影像训练S2RESCA及Dsen2模型,根据最高分辨率波段对每一景影像随机采取500个影像块,尺寸为192像素×192像素,共获取2 000个影像块,90%用于训练,剩余10%用作验证。当s=2时,影像块由A和B组成,尺寸为32像素×32像素,其余均一致。Dsen2和S2RESCA均基于迁移学习训练2 250次,训练批次为16,初始学习率为10-4,损失在10次不变时,当前学习率会乘以衰减系数0.1来自动减小。采用平均绝对误差(MAE)作为L1损失函数,有助于生成高频纹理信息,增强空间细节;采用梯度下降优化算法Nadam作为优化器,有助于增加稳定性和学习效率,采用均方根误差(MSE)测量损失。实验在GTX1660s进行,Python Tensorflow2.0编程实现。

在S2RESCA及Dsen2模型训练完成之后,按式(2),定量评价4种影像重建方法在20 m和60 m波段上的结果。结果如表3所示,加粗字体为本研究提出的方法。其中S2RESCA训练过程的损失变化情况如图3所示,随着迭代次数的增加,训练与验证损失同步下降。

图3 S2RESCA模型训练过程MSE损失变化情况

表3 4种方法重建影像后各波段的RMSE、SRE结果

相比于原始Dsen2神经网络模型,S2RESCA模型在各波段上均有较好的表现,B1和B9(60 m)质量提升明显,RMSE整体降低15.3%,SRE提升1.16%;其余波段(20 m)RMSE整体降低5.1%,SRE无明显差别,表明S2RESCA模型重建影像在数值上更接近真实值。SupReMe重建影像各波段的RMSE均比Bicubic低,其整体的SRE仅比Dsen2低3.6%,重建效果也不错。图4(c)～图4(f)显示了4种算法应用于同一原始影像生成空间分辨率为10 m的影像结果。对于宽度在120～240 m范围的河流,空间分辨率为60 m的B1和B9波段(图4(b)),仅2～4个像素覆盖。当水域宽度小于10个像素时,则该区域影像像元易受邻近效应的影响。相比于图4(b)、图4(c),图4(d)、图4(e)、图4(f)池塘或河流边缘更为清晰,纹理细节明显增强。综合考虑RMSE、SRE及清晰度,4种方法的重建精度是S2RESCA>Dsen2>SupReMe>Bicubic,S2RESCA方法重建的影像更接近真实情况。故本研究采用S2RESCA方法对Sentinel-2/MSI进行超分辨率重建,并研究对水质反演精度的影响。

图4 4种算法对60 m波段重建到10 m的合成显示图结果

2) S2RESCA超分辨率重建影像的质量评价。采用变异系数(coefficient of variation,CV)分析影像重建前后水域与近岸像元的空间异质性,CV为数据标准差与数据平均数的比的百分值,可检测像元的离散程度,是一个无量纲的值,可消除测量尺度和量纲的影响[23]。以水域为例,根据真彩色影像裁剪河岸边缘内较为均一的水体区域(图5),研究发现,远离岸边的水体,其超分辨率重建影像与原始影像各波段的CV变化基本一致(区域1、区域3),表明影像重建前、后空间异质性变化基本一致,邻近效应可忽略不计;而近岸水体(区域2、区域4),重建影像各波段的CV均小于原始影像,且分辨率越低的波段(B1、B9),CV降低更为明显,表明空间异质性明显降低。整体而言,原始影像超分辨率重建后,其20 m波段CV降低7.6%,60 m波段CV降低31.6%。

图5 影像重建前后各波段变异系数对比结果

3.2 水质反演及评价

1)水质参数特征波段。浊度的敏感波段一般在可见光范围内,因此需要去掉780 nm后的波段,其他参数则需去掉B11、B12远红外波段,剩余波段通过遍历两波段组合,计算影像重建前后每组最优波段组合与不同实测水质参数的相关系数绝对值,绘制二者对应波段组合的相关系数差值矩阵热力图来进行敏感性分析(图6),并去除差值小于等于0的区域。各实测水质参数的样本数量与对应的每组特征波段数量均为80个。如图6所示,浊度无明显变化,氨氮、CODMn、总氮和总磷的相关性分别在含有B1与B8、B7与B8、B1与B8A和B1的波段组合中明显增强;最高提升的相关系数分别为0.086、0.082、0.169和0.089。差值矩阵表明,相关系数明显升高的波段组合与参数的敏感性更高。浊度、氨氮、CODMn、总氮和总磷与原始影像相关系数最高分别为0.668、0.725、0.368、0.687及0.639,与重建影像相关系数最高分别为0.668、0.725、0.467、0.726及0.707。重建影像与各参数的最高相关系数均高于或等于原始影像,并达到显著水平(p<0.05)。根据相关系数降序做特征选择(表4),选取前10个特征波段,氨氮、CODMn、总氮、总磷的特征波段中主要含有的波段分别为B4与B8、B7与B8、B1、B1,与图8中相关系数明显升高的特征波段相比,二者均有共同波段(B8、B7与B8、B1、B1),表明应用深度学习超分辨算法可提高影像与不同参数的相关性,且基本上使其最佳特征波段的相关系数得到显著提高。

图6 超分辨影像与原始影像之间各水质参数与两波段组合或波段的相关系数差值矩阵

表4 各水质参数前10个特征波段

2)水质遥感反演模型及评价。分别构建单元线性、多元线性及DNN回归模型对浊度、CODMn、总氮、总磷及氨氮参数进行遥感反演。单元回归模型自变量为表4各参数首个相关性最高的特征波段。多元回归模型和DNN模型自变量为表4中所有特征波段,其中DNN模型训练5 000次。每种参数的实测样本以及对应的特征波段均为160个,打乱顺序后,60%用于构建模型,40%用于验证,计算实测值与预测值的R2、RMSE及sMAPE。

3种模型的结果表明(表5),基于超分辨率影像构建的水质反演模型精度均高于原始影像。相比于原始影像,3种模型基于超分辨率影像反演的浊度、CODMn、总氮、总磷及氨氮,平均RMSE分别降低了2.31%、11.03%、8.09%、13.33%及4.25%,平均sMAPE分别降低了7.03%、4.88%、11.85%、8.41%及1.38%。3种模型中,DNN模型精度最高,在水质遥感监测中具有较大优势,最差的是单元统计回归模型。基于超分辨影像,DNN模型反演水质参数实测值与预测值的散点分布情况如图7所示,各参数的分布较为均匀,模型未出现过饱和现象,具有良好的稳定性。

图7 水质参数DNN模型实测值与预测值的拟合结果

表5 3种模型对影像重建前后的水质遥感反演精度评价结果

3.3 水质参数制图与分析

为了检验DNN模型的推广适用性,使用2021年2月7日位于上海市的一景Sentinel-2/MSI影像进行水质参数反演制图,主要区域包含黄浦江、兰香湖、淀山湖。采用地表水环境质量标准(GB 3838—2002、GHZB1—1999)将水质分为5类,并进行评价。

黄浦江结果显示,黄浦江水体的浊度、CODMn、总氮、总磷及氨氮的反演制图(图8(b)～图8(f))的大部分水体处于健康状态,整体较为清洁干净,与上海市水资源公报的结果一致。如图8(c)所示,CODMn的浓度均在Ⅱ类水质范围内,属于健康范围。总磷的水质级别从Ⅰ级到Ⅳ级不等,值得注意的是,磷主要来源于人类产生的含磷城市污水或者含磷化肥,从图8(e)可以看出,城镇区域的总磷含量明显高于郊区,表明生活污水的排放会影响城镇区域河段的总磷含量。图8(d)总氮浓度与总磷浓度分布相似,图8(f)氨氮的水质级别从Ⅰ级到Ⅲ级不等,处于健康水平。

图8 水质参数浓度空间分布

图8(g)～图8(k)分别为浊度、CODMn、总氮、总磷及氨氮,带有数字的圈表示水质分类级别。淀山湖与兰香湖结果显示,淀山湖(图8(a)的标识1)河口的浊度较大(图8(g)),这主要是由于上游河口的流水携带了大量的泥沙,水向外延伸时泥沙逐渐沉淀,浊度下降,呈扇形分布。兰香湖(图8(a)的标识2)各项水质参数(图8(g)～图8(k))均达到Ⅰ类水质或Ⅱ类水质的标准。根据资料显示,兰香湖是一个人工生态湖泊,常处于清澈的状态,水质普遍较好,在此类较小型的水体反演制图结果与实测资料基本一致。其他水质参数在此也有类似的分布纹理,与浊度有一定的关联性。

4 结束语

本文提出了一种基于迁移学习结合注意力机制CBAM模块的超分辨重建模型S2RESCA,通过降尺度方法,使用SRE与RMSE指标对比了其他算法,该模型可对Sentinel-2/MSI地表反射率产品实现精度更高的超分辨率重建。所有重建波段与原始波段的RMSE均为最低,SRE均为最高,明显提高了河湖的清晰度和空间细节。此外,本文评价了超分辨率重建影像的质量及进行水质遥感反演的精度,得出以下结论。

1)影像通过S2RESCA模型重建后,近岸水体像元空间异质性明显降低。特别是B1和B9波段CV降低31.6%,其余波段CV降低7.6%。

2)水质参数特征波段分析结果表明,浊度无明显变化,CODMn、总氮、总磷和氨氮相关性最高的特征波段分别是含有B7与B8、B1、B1、B8的波段组合。影像超分辨率重建后,水质参数相关系数提升明显的均为相关性最高的特征波段,即为敏感性波段,分别提升0.082、0.169、0.089和0.086。

3)水质参数遥感反演模型精度对比结果表明,对于浊度、CODMn、总氮、总磷和氨氮的超分辨率影像反演,各模型预测精度均有提升,其中DNN模型对各水质参数的反演精度更高,R2分别为0.85、0.74、0.67、0.76和0.78。通过对上海市流域水系的水质参数浓度制图分析,黄浦江、兰香湖及淀山湖水质参数浓度分布与水环境的实际情况基本一致,验证了DNN模型在不同水域均可进行有效的水质评价。

本研究成果可实现时空连续的高分辨率(10 m)水质参数浓度制图,为地表水环境管理决策提供有用的监测工具。相比水色参数SPM、Chla及CDOM,低光学信号水质参数遥感反演影响因素更为复杂。尽管本文研究了深度学习超分辨率算法及水质参数模型的影响,但水体颜色、水体类型及水色参数等与其之间的共变关系仍有待研究。