基于多源遥感及气象数据的河流非光学活性水质参数反演模型研究

doi:10.11873/j.issn.1004-0323.2024.1.0120

基于多源遥感及气象数据的河流非光学活性水质参数反演模型研究

兑紫璇^,¹^,², 王卿³, 王敏³, 张璟⁴, 顾倩荣^,¹

1.中国科学院上海高等研究院低碳转化科学与工程重点实验室碳数据与碳评估中心，上海 201210

2.中国科学院大学，北京 100049

3.上海市环境科学研究院，上海 200233

4.江苏省司法警官高等职业学校，江苏南京 212008

Research on the Retrieval Model of Non-optically Active Water Quality Parameters of Rivers based on Multi-source Remote Sensing and Meteorological Data

DUI Zixuan^,¹^,², WANG Qing³, WANG Min³, ZHANG Jing⁴, GU Qianrong^,¹

1.Shanghai Carbon Data Research Center，Key Laboratory of Low - Carbon Conversion Science & Engineering，Shanghai Advanced Research Institute，Chinese Academy of Sciences，Shanghai 201210，China

2.University of Chinese Academy of Sciences，Beijing 100049，China

3.Shanghai Academy of Environmental Sciences，Shanghai 200233，China

4.Jiangsu Provincial Judicial Police Officer Higher Vocational College，Nanjing 212008，China

通讯作者: 顾倩荣，博士，高级工程师，主要从事卫星遥感数据应用研究。E⁃mail: guqr@sari.ac.cn

收稿日期: 2022-10-13 修回日期: 2023-12-18

基金资助:

上海市2020年度科技创新行动计划社会发展科技攻关项目. 20dz1204302
国家自然科学基金面上项目“基于多维数据的城市群碳达峰规划评估模型”. 52178060

Received: 2022-10-13 Revised: 2023-12-18

作者简介 About authors

兑紫璇，硕士研究生，主要从事卫星遥感数据应用研究E⁃mail：duizx@sari.ac.cn , E-mail：duizx@sari.ac.cn

摘要

针对传统河流水质监测成本高、地面监测站点稀疏等问题，基于哨兵2卫星多光谱遥感数据，结合MODIS地表温度、植被指数、气溶胶光学厚度数据产品，以及ERA5气象数据产品中的地表风速数据，以非光学活性参数溶解氧（DO）、化学需氧量（COD）和氨氮（NH₃-N）的地面水质监测站的监测数据为参照，采用支持向量回归（SVR）、随机森林（RF）和多层感知机（MLP）3种机器学习方法，通过对比实验选出每种水质参数最优反演模型及其对应的输入特征组合。模型性能测试实验结果表明，利用多源遥感及气象数据反演出的DO、COD以及NH₃-N的决定系数（R²）分别为0.896、0.781、0.529，均方根误差（RMSE）分别为0.263 mg/L、0.383 mg/L、0.061 mg/L。与仅使用哨兵2卫星多光谱遥感数据的反演结果相比，R²分别提高了7.04%、19.05%、18.34%，RMSE分别降低了34.58%、37.42%、14.08%。表明多源遥感及气象数据对提高DO、COD以及NH₃-N水质参数反演准确性有重要意义。最后，模型鲁棒性评估实验表明，当模型训练数据的代表性与全局数据较为接近时，训练的模型具有较好的时空鲁棒性。

关键词： 机器学习 ; 多源遥感数据 ; 溶解氧 ; 化学需氧量 ; 氨氮 ; 水质反演模型

Abstract

In view of the high cost of traditional river water quality monitoring and the sparse ground monitoring stations， based on Sentinel-2 satellite multispectral remote sensing data， combined with MODIS surface temperature， vegetation index， aerosol optical thickness data products， and the surface wind speed data in ERA5 meteorological data products， the monitoring data of the surface water quality monitoring stations with non-optical active parameters Dissolved Oxygen （DO）， Chemical Oxygen Demand （COD） and ammonia nitrogen （NH₃-N） are taken as reference， three machine learning methods， Support Vector Regression （SVR）， Random Forest （RF） and Multilayer Perceptron （MLP）， were used to select the optimal inversion model of each water quality parameter and its corresponding input feature combination through comparative experiments. The experimental results of the model performance test show that the determination coefficients （R²） of DO， COD and NH₃-N are 0.896，0.781 and 0.529， respectively，and the Root Mean Square Error（RMSE） are 0.263 mg/L，0.383 mg/L and 0.061 mg/L， respectively. Compared with the retrieval results using only Sentinel-2 multi-spectral remote sensing data， R² increased by 7.04%， 19.05% and 18.34% respectively， and RMSE decreased by 34.58%， 37.42% and 14.08% respectively. It shows that multi-source remote sensing and meteorological data are of great significance to improve the retrieval accuracy of DO， COD and NH₃-N water quality parameters. Finally， the model robustness evaluation experiment shows that the trained model has better space-time robustness when the representativeness of the model training data is close to the global data.

Keywords： Machine learning ; Multi-source remote sensing data ; Dissolved oxygen ; Chemical oxygen demand ; Ammonia nitrogen ; Water quality retrieval model

PDF (4306KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

兑紫璇, 王卿, 王敏, 张璟, 顾倩荣. 基于多源遥感及气象数据的河流非光学活性水质参数反演模型研究. 遥感技术与应用[J], 2024, 39(1): 120-133 doi:10.11873/j.issn.1004-0323.2024.1.0120

DUI Zixuan, WANG Qing, WANG Min, ZHANG Jing, GU Qianrong. Research on the Retrieval Model of Non-optically Active Water Quality Parameters of Rivers based on Multi-source Remote Sensing and Meteorological Data. Remote Sensing Technology and Application[J], 2024, 39(1): 120-133 doi:10.11873/j.issn.1004-0323.2024.1.0120

1 引言

水污染对生态环境、可利用淡水资源和人类健康构成严重威胁。因此，定期对河流水质进行监测十分必要。传统的实地采样、实验室分析的监测方法人力和物力耗费巨大，而且效率十分低下。采用地面水质自动监测站的方式，因站点数量有限，且通常稀疏地分散在河流的不同区域，所以这些站点监测到的水质数据难以真实地反映出河流较大区域范围内的水质情况。与此相对，卫星遥感具有覆盖范围广、监测效率高等优点，随着利用遥感数据反演水质研究的进展，卫星遥感反演已广泛应用于水质监测领域。

水质参数根据其是否对某些特定光谱谱段具有明显的光学响应，分为光学活性参数和非光学活性参数两大类^［1］。对于光学活性水质参数，可以根据其固有光学特性与水面反射率之间的关系建立遥感反演模型。目前对于光学活性参数的遥感反演，已有大量研究且取得较好的成果，如叶绿素^［2-3］、总悬浮物^［4-5］和有色溶解有机物^［6］等。对于非光学活性水质参数，由于水体中多种水质成分和沉积物存在的相互影响，以及不同流域水体类型、水深以及周围环境的差异，导致非光学活性参数对水面反射率的影响呈现出复杂多变的特点，这使得直接利用遥感获得的多光谱特征对其进行反演非常困难，而且传统的经验统计法、物理分析法和半经验法的泛化性较差，难以取得较好的反演结果^［7］。相反通过数据驱动，能自动提取出数据之间复杂非线性关系的机器学习方法，非常适合用于提取非光学活性水质参数和多光谱遥感数据之间的相关性，并且能在一定程度上避免大气和其它环境因素带来的干扰，从而实现对这些水质参数的反演^［7-8］。

Wang等^［9］发现SPOT5卫星的红光、绿光和近红外谱段的观测值与非光学活性水质参数中的化学需氧量（Chemical Oxygen Demand，COD）和氨氮（Ammonia Nitrogen，NH₃-N）数值之间具有较强的相关性。Kuan等^［10］利用经过超分辨率算法处理后的哨兵2多光谱数据，实现了对内陆河流的NH₃-N和COD等水质参数的反演。Guo等^［11］用机器学习的方法，从多光谱遥感数据中反演出城市小型水体的溶解氧（Dissolved Oxygen，DO）和COD等水质参数。这些研究表明，DO、COD和NH₃-N水质参数与多光谱遥感数据之间存在较强的相关性关系^［12］，揭示了利用数据驱动的机器学习方法，从多光谱遥感数据中反演出这些非光学活性参数的可能性。

目前已有的应用卫星遥感数据反演非活性水质参数的研究，往往只使用卫星的多光谱遥感数据作为反演模型的输入，缺少考虑河流周围环境因素对水质参数的影响。理论上卫星观测到水面反射率应仅受水体水质的影响，但实际上卫星的观测数据还会受到大气中气溶胶光学厚度（Aerosol Optical Depth，AOD）和因风速引起水面浪高差异的影响^［7］。此外，河流周围的地表温度和植被覆盖程度会影响水体中的营养状态，进而影响各种水质参数的数值^［13］。因此，在反演河流水质参数时，河流及附近的风速、气溶胶光学厚度、地表温度、植被指数等环境因素的影响不应被忽略。

实验以谷歌地球引擎（Google Earth Engine，GEE）云遥感处理平台^［14］为分析处理平台，利用哨兵2卫星^［15］多光谱地表反射率数据产品、中分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer，MODIS）的地表温度、植被指数，和气溶胶光学厚度数据产品，以及ERA5气象数据，作为多源数据输入，以中国环境监测总站发布的全国主要流域地表水水质监测数据为目标，使用机器学习方法训练反演模型，实现对河流的DO、COD以及NH₃-N水质参数的反演。实验结果表明，使用多源数据的DO，COD和NH₃-N反演结果的确定系数（The Coefficient of Determination，R²）分别为0.896，0.781，0.529，均方根误差（Root Mean Square Error，RMSE）分别为0.263 mg/L，0.383 mg/L，0.061 mg/L，与只使用哨兵2多光谱遥感数据的反演结果相比，R²提高分别了7.04%，19.05%，18.34%，RMSE分别降低了34.58%，37.42%，14.08%。

2 数据获取与处理

2.1　地面水质监测数据

地面水质监测数据采用中国环境监测总站公开发布的全国主要流域144个地表水水质监测站点，2017年1月至2019年12月间的实时水质监测数据集^［16］中，包含的DO、COD和NH₃-N水质参数监测数据。数据集中的144个水质监测站点，从高纬度到低纬度，覆盖我国东部与南部地区的主要河流和湖泊，具有较广泛的空间代表性。

2.2　多源遥感数据与气象数据

2.2.1　哨兵2多光谱遥感数据

哨兵2^［15］包括两颗相同的太阳同步轨道卫星，分别于2015年和2017年发射，每颗卫星上搭载一个具有13个光谱谱段的多光谱成像仪（Multispectral imager，MSI），成像的空间分辨率分别为10 m、20 m和60 m。

本研究使用经过辐射定标和大气校正的哨兵2 MSI L2A级数据产品^［17］，该数据产品代表哨兵2观测到的地表反射率。数据产品中的第10谱段数据，因其主要包含高空卷云的信息，所以被丢弃不用。另外，利用哨兵2提供的质量参数筛选出云量小于10%的数据，用于反演研究。

2.2.2　MODIS数据产品

（1）地表温度

地表温度（Land Surface Temperature，LST）来自MODIS的MOD11A1^［18］（Terra）和MYD11A1^［19］（Aqua）两个数据产品，时间分辨率为每天，空间尺度为1 000 m。

在根据地面水质监测数据的坐标筛选地表温度数据时，由于云的遮挡，有些值可能是null，所以将Terra和Aqua两颗卫星的地表温度数据产品搭配使用，以减少数据的缺失。

（2）植被指数

植被指数（Vegetation Index，VI）来自MODIS MOD13Q1^［20］（Terra）和MYD13Q1^［21］（Aqua）两个数据产品。它们相互配合，可以得到时间分辨率为8 d，空间分辨率为250 m的两种植被指数，分别为归一化差异植被指数（Normalized Difference Vegetation Index，NDVI）和增强植被指数（Enhanced Vegetation Index，EVI）。

（3）气溶胶光学厚度

气溶胶光学厚度（Aerosol Optical Depth，AOD）来自MODIS的MCD19A2^［22］数据产品，时间分辨率为每天，空间分辨率为1 km。该产品提供波长分别为47 nm和55 nm的两种光学厚度数据。

2.2.3　ERA5气象数据

ERA5是欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）生产和发布的最新一代全球气象再分析资料，（https：//www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5）。 ERA5提供空间分辨率为0.1°，时间分辨率为1 h的全球大气温度、湿度、气压、风速、降水量等气象参数。研究采用其中10 m高度的u、v风速数据，来计算湖泊河流的水面风速（Wind speed，WS）。

2.3　地面水质监测数据与多源数据的时空匹配

2.3.1　空间上的匹配

地面水质监测数据与多源遥感及气象数据在空间上的匹配，分为以下3个步骤：

步骤（1）：在GEE平台^［14］上，从哨兵2、MODIS和ERA5数据产品中，根据地面水质监测数据的坐标，筛选出2017年1月至2019年12月（水质监测数据集所覆盖的时间段）内的多源遥感及气象数据。

步骤（2）：在GEE平台^［14］上，将筛选出的MODIS和ERA5数据的投影方式，转换成哨兵2数据产品所使用的EPSG：32649投影方式。

步骤（3）：以哨兵2的10 m分辨率的谱段2数据为基准，通过重采样和双三次插值方法^［23］，将哨兵2其它低空间分辨率的谱段数据，MODIS的地表温度、植被指数，和气溶胶光学厚度，以及ERA5的气象数据，插值到哨兵2谱段2数据的坐标处，最终完成多源遥感数据和气象数据在空间上的匹配。

2.3.2　时间上的匹配

经过空间匹配之后的多源遥感及气象数据需要进一步在时间上与地面水质监测数据进行匹配。因为每种数据产品的时间分辨率不同，且对时间变化的敏感程度不同，我们根据每种数据产品的特点，如表1所示，确定其在时间上与地面水质监测数据进行匹配的时间窗口大小。

表1 水质监测数据和多源遥感及气象数据匹配的时间窗口

Table 1 Time window for matching water quality monitoring data and multi-source remote sensing and meteorological data

数据名称	哨兵2	MODIS			ERA5
数据名称	哨兵2	地表温度	植被指数	气溶胶	ERA5
时间窗口	4 h	同一天	同一天	同一天	1 h

新窗口打开| 下载CSV

在以往研究中，根据地面水质监测数据的情况，通常将与遥感数据匹配的时间窗口大小选择为4 h^［11］、1 d^［10］，或7 d^［2］。本研究为了得到较好的反演效果，将哨兵2数据与地面水质监测数据匹配的时间窗口大小设置为4 h。

2.3.3　时空匹配后数据的异常值处理

对于经过与多源遥感及气象数据进行时空匹配后筛选出的地面水质监测数据，根据其数值分布的概率，进一步使用异常检测算法EllipticEnvelope^［24］将过高或过低的异常水质监测值剔除，最后得到用于训练反演模型的地面水质监测数据样本集。

筛选出的地面水质监测数据集，异常值剔除前后的统计信息的差异如表2所示。

表2 与多源遥感和气象数据时空匹配后筛选出的地面水质监测数据样本集

Table 2 The ground water quality monitoring sample data set filtered by spatial and temporal matching with the multi-source remote sensing and meteorological data

状态	水质参数	异常值处理比例	样本数	平均值/(mg $∙$ L^-1)	标准偏差	最小值/(mg $∙$ L^-1)	最大值/(mg $∙$ L^-1)
异常值剔除前	DO	-	270	8.78	2.43	3.76	18.07
	COD	-	270	2.29	2.34	0.10	15.60
	NH₃-N	-	270	0.93	1.42	0.00	10.00
异常值剔除后	DO	0.08	249	8.29	1.54	4.70	12.10
	COD	0.10	243	1.68	1.33	0.10	5.25
	NH₃-N	0.20	216	0.55	0.65	0.00	2.80

新窗口打开| 下载CSV

3 水质参数反演模型

研究使用支持向量回归^［25］（Support Vector Regression，SVR）、随机森林^［26］（Random Forest，RF）和多层感知机^［27］（Multilayer Perceptron，MLP）3种机器学习模型，利用多源遥感和气象数据，反演DO、COD及NH₃-N的水质参数。同时，使用传统的多元线性回归（Multiple Linear Regression，MLR）模型^［28］，作为与机器学习模型反演性能对比的参照^［28］，以分析机器学习模型对这3种非光学活性水质参数反演结果准确性的提升情况。

3.1　支持向量回归

基于支持向量^［29］概念的支持向量回归，将低维特征空间中的输入通过非线性变换映射到高维特征空间中，并如图1所示，寻找确保预测值与实际观测值的差异在偏差ε范围内的回归超平面f（x）^［25］。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 支持向量回归示意图

Fig.1 Schematic diagram of support vector regression

因为仅有作为支持向量的关键训练样本用于决定回归超平面的位置，所以SVR的性能不会过度依赖训练数据量的多少^［30］，因此在可获得样本数量有限的遥感任务中SVR具有一定优势。

3.2　随机森林

随机森林由多棵不存在强依赖关系的决策树组成，是一种基于装袋法思想^［31］的集成学习方法^［26］。如图2所示，其决策树在每个非叶子节点选择一种特征属性进行测试，并将测试结果通过分支输出，随机森林最终的输出值由多棵决策树共同决定，同时可以对预测变量之间的相互作用进行模拟，因此对输入变量中的异常值具有一定的鲁棒性^［30］。随机森林可以随机选择测试使用的特征子集，可以适用于高维度数据的处理。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 随机森林示意图

Fig.2 Schematic Diagram of Random Forest

3.3　多层感知机

多层感知机如图3所示，由输入层、输出层和至少一层的隐含层构成，其中隐藏层与输入层、输出层的神经元完全连接。通过有监督的训练改变网络中神经元之间的连接权重，MLP可以成为能以任意精度^［32-33］逼近集合上任何函数^［34］的通用函数逼近器^［35］。在复杂非线性拟合问题上，MLP可以提供较为精确的近似解。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 多层感知机示意图

Fig.3 Schematic diagram of multi-layer perceptron

4 实验设计与数据准备

4.1　实验数据准备

4.1.1　模型训练及测试数据集

先将地面水质监测数据与多源遥感及气象数据进行时空匹配，筛选出时空一致的数据，再根据地面水质监测值的概率分布信息剔除掉其中的异常值，最终形成用以反演模型训练和测试的样本数据集。

从样本数据集中随机抽取80%作为模型的训练集，剩余的20%作为测试集，用以模型的训练和测试。为了防止训练过程中模型过拟合，并评估模型的泛化性能，总共生成10组训练和测试集，将模型在10组测试集上性能指标的平均值作为模型最终的性能指标。

4.1.2　模型鲁棒性评估数据集

地面水质监测数据、多源遥感和气象数据中，除去作为模型训练和测试用的样本数据集之外，剩余的数据，根据地面水质监测数据与多源遥感和气象数据时空匹配的情况，可以分为以下4类，分别为：

（1）对某个地面监测站，虽然有空间上能与之匹配的多源遥感及气象数据，但如果它的某一特定时刻的水质监测数据，没有时间上能匹配的多源遥感数据，那么它在这个时刻的监测数据，称为时间未匹配的水质监测数据。

（2）对某个地面监测站，能与其空间匹配的所有多源遥感及气象数据，可以构成一个集合。如果集合不为空，并且至少有一组多源遥感及气象数据，与这个监测站的水质监测数据在时间上也能匹配，那么对于集合中剩余的，所有与这个监测站数据，在时间上不能匹配的多源遥感及气象数据，被称为时间未匹配的多源遥感及气象数据。

（3）对某个地面监测站，如果没有任何多源遥感及气象数据，能与其在空间上匹配，这个监测站的所有水质监测数据，称为空间未匹配的水质监测数据。

（4）对某个地面监测站，能与其空间匹配的所有多源遥感及气象数据，可以构成一个集合。如果集合不为空，并且没有任何一组多源遥感及气象数据，能与这个监测站的水质监测数据，在时间上匹配上，那么这个集合中的所有多源遥感及气象数据，称为空间未匹配的遥感数据。

用上述4种时空不匹配的数据集，对反演模型的鲁棒性进行评估。

模型训练和测试数据集、模型鲁棒性评估数据集所包含的地面水质监测数据、多源遥感及气象数据的情况，如表3所示。

表3 模型鲁棒性评估实验所用的数据集

Table 3 Data set used in model robustness assessment experiments

实验名称	数据集	水质参数/(个/组）
实验名称	数据集	DO	COD	NH₃-N
	所有地面水质监测数据	812 441	811 147	804 843
时空匹配实验	训练数据	249	243	216
时空匹配实验	测试数据	75	73	65
时间未匹配实验	水质监测数据	234 470	235 918	235 598
时间未匹配实验	多源遥感及气象数据	142	142	142
空间未匹配实验	水质监测数据	575 436	580 629	574 486
空间未匹配实验	多源遥感及气象数据	219	219	219

新窗口打开| 下载CSV

4.2　实验设计

4.2.1　建模参数选取

本研究在建立水质反演模型时，选择哨兵2的多光谱数据（MSI），MODIS的地表温度（LST）、植被指数（VI）、气溶胶光学厚度（AOD），以及ERA5的10 m风速（WS）作为候选的建模参数。其中，哨兵2的多个谱段的反射率可以体现水体中不同成分对光的吸收和反射情况，解析其与DO、COD和NH₃-N参数值之间的关系，可以建立反演模型。

除了多光谱反射率之外，反映河流及周围环境因素的LST、VI、AOD和WS也被考虑进去。首先，大气状态和水面平整度与卫星观测到的反射率有关，卫星观测到的水表反射的太阳辐射需要在大气、水面和水体之间传输，在这个过程中，水面的起伏波动，以及大气中的气溶胶都会改变其原本的传输路径，导致观测到的反射率与实际的水表反射率之间产生偏差。另外，水体的富营养化程度和DO、COD和NH₃-N之间也存在着密切的关联，河流的富营养化会导致水中藻类大量繁殖，影响水体的光学特性、氧气含量和各种生物过程，进而影响DO、COD和NH₃-N参数值。综上，挑选能够影响水面平整度的风速、反映大气状态的气溶胶光学厚度，以及影响和反映水体中营养状态的地表温度和植被指数作为反演模型的辅助参数，以提供更加全面、准确的水质反演结果。

4.2.2　模型训练及性能评估实验

为了从哨兵2的MSI，MODIS的LST、VI、AOD，以及ERA5的WS等多源数据中，选择出每种水质参数反演模型最合适的输入特征组合，以哨兵2的多光谱数据为基础的输入特征，评估增加考虑环境因素和大气因素对反演模型性能的影响，总共可形成四大类共14种输入特征的组合，分别为：

（1）仅以哨兵2的多光谱为模型输入，即为MSI。

（2）增加河流周边地表温度和植被指数环境因素的影响，即可形成MSI+LST、MSI+VI、MSI+LST+VI，共3种输入特征组合。

（3）增加河流上空大气气溶胶光学厚度和10 m风速的影响，即可形成MSI+AOD、MSI+WS、MSI+AOD+WS，共3种输入特征组合。

（4）同时增加环境因素和大气因素的影响，即可形成MSI+LST+AOD、MSI+LST+WS、MSI+VI+AOD、MSI+VI+WS、MSI+LST+AOD+WS、MSI+VI+AOD +WS、MSI+LSI+VI+AOD+WS，共7种输入特征组合。

将以上14种输入特征组合，作为SVR、RF和MLP模型的输入，在训练集上分别训练DO、COD及NH₃-N水质参数的反演模型。

3种水质参数的SVR、RF和MLP反演模使用Sklearn库^［36］实现，并使用网格搜索法^［37］自动搜索各个模型的最优超参数值。

当14种输入特征组合对应的各水质反演模型训练完成之后，使用测试集对模型反演结果的准确性进行测试和评估，选出最佳输入特征组合和对应的水质反演模型，作为模型训练实验的最终结果。

4.2.3　模型鲁棒性评估实验

如果模型的时空鲁棒性较好的话，当输入为不包含在模型训练和测试样本集内的多源遥感及气象数据时，模型反演出的水质参数的统计特征，也应与对应地面监测数据的统计特征相近似。

基于以上假设，分别使用4.1.2中说明的时间未匹配和空间未匹配的数据，开展时间未匹配和空间未匹配实验，以评估反演模型的时空鲁棒性。

因时间、空间未匹配实验缺少地表水质监测数据作为反演结果比较的“真值”，所以采用对比用箱形图表示的实验数据分布的统计特征，以评估模型的时空鲁棒性。同时，将模型训练过程作为时空匹配实验，将其结果的统计特征作为评估反演模型时空鲁棒性的参照。

4.3　评价指标

使用确定系数（The Coefficient of Determination，R²）^［38］、平均绝对误差（Mean Absolute Error，MAE）^［39］和均方根误差（Root Mean Square Error，RMSE）^［40］这3种指标，评估模型反演结果的准确性。指标分别定义如下：

R^{2} = 1 - \frac{\sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}}{\sum_{i = 1}^{n} {(y_{i} - \bar{y})}^{2}}

（1）

M A E = \frac{1}{n} \sum_{i = 1}^{n} |y_{i} - {\hat{y}}_{i}|

（2）

R M S E = \sqrt[]{\frac{1}{n} \sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}}

（3）

其中： $y_{i}$ 和 ${\hat{y}}_{i}$ 分别为第i个地面水质监测数据和其对应的模型反演结果； $\bar{y}$ 为地面水质监测数据的平均值；n为样本总数。

5 结果分析

5.1　模型的训练及性能评估结果

5.1.1　最优超参数筛选

经实验筛选出的3种机器学习水质参数反演模型的超参数配置，如表4所示。

表4 3种机器学习水质参数反演模型的最优超参数表

Table 4 Optimal hyperparameter table of three machine learning water quality parameter retrieval models

机器学习模型	参数名称	非光学活性水质参数
机器学习模型	参数名称	DO	COD	NH₃-N
SVR	C	100	10	10
SVR	gamma	7	3	3
RF	决策树个数	15	26	15
RF	树的深度	46	15	10
MLP	隐藏层数	6	5	5
MLP	每个隐藏层神经元数	100	100	100

新窗口打开| 下载CSV

5.1.2　各模型水质参数反演的最优输入特征组合

表5为SVR、RF和MLP模型在14种输入特征组合下的DO、COD和NH₃-N水质参数反演结果。从表中可以看到，与仅使用哨兵2的MSI作为输入相比，输入数据中在叠加MODIS、ERA5的数据后，各模型的性能均有所提升。对于不同水质参数，能得到最优反演性能的模型不同，及对应的输入特征组合也不同，具体分析如下：

表5 SVR、RF和MLP模型在14种输入特征组合下的DO、COD和NH₃-N反演结果

Table 5 DO， COD and NH₃-N retrieval results of SVR， RF and MLP models for 14 combinations of input features

序号	输入数据	指标	DO			COD			NH₃-N
序号	输入数据	指标	SVR	RF	MLP	SVR	RF	MLP	SVR	RF	MLP
1	MSI	R²	0.837	0.845	0.832	0.649	0.699	0.656	0.242	0.447	0.485
		MAE/(mg $∙$ L^-1)	0.308	0.393	0.426	0.477	0.506	0.474	0.183	0.144	0.139
		RMSE/(mg $∙$ L^-1)	0.402	0.384	0.414	0.594	0.513	0.612	0.100	0.071	0.066
2	MSI+VI	R²	0.837	0.852	0.357	0.650	0.685	0.744	0.191	0.523	0.353
		MAE/(mg $∙$ L^-1)	0.302	0.375	0.729	0.485	0.502	0.415	0.188	0.139	0.151
		RMSE/(mg $∙$ L^-1)	0.407	0.366	1.711	0.596	0.537	0.444	0.106	0.063	0.080
3	MSI+LST	R²	0.866	0.860	0.764	0.611	0.654	0.605	0.232	0.292	0.386
		MAE/(mg $∙$ L^-1)	0.306	0.384	0.465	0.506	0.546	0.543	0.182	0.173	0.154
		RMSE/(mg $∙$ L^-1)	0.336	0.347	0.594	0.654	0.592	0.703	0.100	0.092	0.078
4	MSI+VI+LST	R²	0.890	0.879	0.841	0.671	0.673	0.717	0.251	0.232	0.325
		MAE/(mg $∙$ L^-1)	0.282	0.359	0.377	0.472	0.527	0.430	0.184	0.168	0.154
		RMSE/(mg $∙$ L^-1)	0.274	0.300	0.397	0.560	0.558	0.491	0.098	0.096	0.084
5	MSI+AOD	R²	0.872	0.863	0.856	0.676	0.671	0.722	0.307	0.507	0.432
		MAE/(mg $∙$ L^-1)	0.300	0.363	0.368	0.470	0.514	0.435	0.165	0.139	0.143
		RMSE/(mg $∙$ L^-1)	0.320	0.341	0.360	0.553	0.564	0.476	0.091	0.064	0.072
6	MSI+WS	R²	0.802	0.851	0.835	0.648	0.699	0.738	0.289	0.488	0.536
		MAE/(mg $∙$ L^-1)	0.344	0.385	0.377	0.464	0.495	0.431	0.172	0.143	0.129
		RMSE/(mg $∙$ L^-1)	0.497	0.364	0.417	0.599	0.514	0.449	0.091	0.066	0.061
7	MSI+AOD+WS	R²	0.868	0.829	0.733	0.704	0.669	0.781	0.334	0.529	0.502
		MAE/(mg $∙$ L^-1)	0.293	0.402	0.454	0.439	0.515	0.398	0.171	0.137	0.137
		RMSE/(mg $∙$ L^-1)	0.330	0.431	0.647	0.507	0.565	0.383	0.087	0.061	0.065
8	MSI+VI+LST+AOD+WS	R²	0.889	0.868	0.853	0.670	0.660	0.742	0.221	0.278	0.383
		MAE/(mg $∙$ L^-1)	0.289	0.378	0.337	0.487	0.528	0.414	0.184	0.168	0.153
		RMSE/(mg $∙$ L^-1)	0.278	0.332	0.367	0.564	0.578	0.445	0.099	0.092	0.079
9	MSI+LST+WS	R²	0.853	0.839	0.695	0.674	0.642	0.721	0.249	0.320	0.436
		MAE/(mg $∙$ L^-1)	0.320	0.411	0.506	0.479	0.542	0.446	0.178	0.161	0.149
		RMSE/(mg $∙$ L^-1)	0.371	0.404	0.780	0.557	0.618	0.482	0.099	0.087	0.073
10	MSI+LST+AOD	R²	0.896	0.866	0.837	0.644	0.656	0.748	0.314	0.314	0.458
		MAE/(mg $∙$ L^-1)	0.276	0.369	0.368	0.506	0.543	0.438	0.175	0.165	0.143
		RMSE/(mg $∙$ L^-1)	0.263	0.333	0.409	0.611	0.589	0.428	0.088	0.089	0.070
11	MSI+LST+AOD+WS	R²	0.891	0.872	0.857	0.707	0.646	0.736	0.273	0.226	0.381
		MAE/(mg $∙$ L^-1)	0.283	0.370	0.333	0.441	0.542	0.424	0.180	0.170	0.152
		RMSE/(mg $∙$ L^-1)	0.275	0.320	0.350	0.506	0.613	0.470	0.094	0.098	0.078
12	MSI+VI+WS	R²	0.835	0.831	0.817	0.677	0.690	0.764	0.217	0.482	0.293
		MAE/(mg $∙$ L^-1)	0.305	0.382	0.358	0.464	0.502	0.404	0.185	0.145	0.162
		RMSE/(mg $∙$ L^-1)	0.412	0.416	0.466	0.553	0.531	0.412	0.103	0.066	0.092
13	MSI+VI+AOD	R²	0.869	0.856	0.867	0.692	0.667	0.714	0.313	0.459	0.382
		MAE/(mg $∙$ L^-1)	0.291	0.371	0.337	0.439	0.511	0.419	0.175	0.146	0.146
		RMSE/(mg $∙$ L^-1)	0.324	0.360	0.337	0.529	0.571	0.498	0.088	0.069	0.077
14	MSI+VI+AOD+WS	R²	0.879	0.836	0.844	0.690	0.676	0.724	0.307	0.509	0.412
		MAE/(mg $∙$ L^-1)	0.291	0.396	0.344	0.452	0.502	0.409	0.177	0.139	0.145
		RMSE/(mg $∙$ L^-1)	0.303	0.410	0.391	0.526	0.561	0.483	0.089	0.064	0.077

新窗口打开| 下载CSV

溶解氧（DO）：性能最优的模型为SVR，对应的输入特征组合为实验10中的MSI+LST+AOD，其R²为0.896，MAE为0.276 mg/L，RMSE为0.263 mg/L。比只使用哨兵2的MSI作为输入的R²提高7.04%，MAE降低10.39%，RMSE降低34.58%。

化学需氧量（COD）：性能最优的模型为MLP，对应的输入特征组合为实验7中的MSI+AOD+WS，其R²为0.781，MAE为0.398 mg/L，RMSE为0.383 mg/L。比只使用哨兵2的MSI作为输入的R²提高19.05%，MAE降低16.03%，RMSE降低37.42%。

氨氮（NH₃-N）：性能最优的模型为RF，对应的输入特征组合为实验7中的MSI+AOD+WS，其R²为0.529，MAE为0.137 mg/L，RMSE为0.061 mg/L。比只使用哨兵2的MSI作为输入的R²提高18.34%，MAE降低4.86%，RMSE降低14.08%。

5.1.3　各模型在3种水质参数上的最优性能

从表5中挑选出能达到最优反演性能的SVR、RF和MLP模型的性能指标，作为表6的内容。表6中同时列出多元线性回归（MLR）模型的反演性能指标作为对比参照。

表6 不同方法的模型能到的最佳水质参数反演性能

Table 6 The best retrieval performance of water quality parameters that can be obtained by models of different methods

方法	指标	非光学活性水质参数
方法	指标	DO	COD	NH₃-N
MLR	R²	0.383	0.066	0.120
	MAE/(mg $∙$ L^-1)	1.027	1.046	0.251
	RMSE/(mg $∙$ L^-1)	1.650	1.561	0.132
SVR	R²	0.896	0.707	0.334
	MAE/(mg $∙$ L^-1)	0.276	0.441	0.171
	RMSE/(mg $∙$ L^-1)	0.263	0.506	0.087
RF	R²	0.879	0.699	0.529
	MAE/(mg $∙$ L^-1)	0.359	0.495	0.137
	RMSE/(mg $∙$ L^-1)	0.300	0.514	0.061
MLP	R²	0.867	0.781	0.536
	MAE/(mg $∙$ L^-1)	0.337	0.398	0.129
	RMSE/(mg $∙$ L^-1)	0.337	0.383	0.061

新窗口打开| 下载CSV

从DO、COD和NH₃-N3种水质参数反演性能的平均指标来看，整体上MLP模型的R²最高、MAE和RMSE最低。MLP模型反演结果的R²均大于0.536，与作为参照的传统的MLR模型相比，机器学习方法能显著提高水质参数反演的性能。 RF模型反演结果的相关性第二高，R²均大于0.529，RF模型对DO反演结果的性能要比MLP模型稍好。

5.2　模型的鲁棒性评估结果

5.2.1　时空匹配实验（模型训练和测试实验）结果

图4为3种水质参数时间匹配实验结果的箱形图。图中黑色、灰色，和白色箱形，分别表示实验中所有地面水质监测数据、样本中训练模型用的水质监测数据，和模型用样本中作为测试的多源遥感及气象数据反演得到的水质参数值的统计特征。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 所有地面水质监测数据，样本中的模型训练数据、模型用样本中测试数据反演的水质参数的年度对比箱形图

（黑色箱形代表所有站点的水质监测数据，灰色箱形代表样本中训练反演模型的水质监测数据；白色箱形代表模型用样本中测试数据反演出的水质参数；箱形图中，箱体内的实线表示中位数，箱体的下边缘和上边缘分别表示下四分位数（Q1）和上四分位数（Q3），箱形外的上下横线表示最大值和最小值，o表示异常值）

Fig.4 Annual boxplots of all surface water quality monitoring data， model training data in the sample set， water quality parameters retrieved by the models with test data in the sample set

如果灰色箱形与黑色箱形相似，表示在时空匹配实验中，反演模型训练所用数据的统计特征，与所有地面水质监测数据的统计特征相似，即训练数据的代表性较好。如果白色箱形与灰色箱形相似，表示模型用样本中测试数据反演出的水质参数的统计特征，与样本中训练数据的统计特征相似，即模型得到了合理的训练。如果白色箱形与黑色箱形相似，则表示模型的泛化性较好。

从图4可以看出，通过箱形图分析得到的在3种水质参数反演模型的性能，与表5和表6中的结果相同。其中，DO反演模型训练数据的代表性最好，DO模型得到了合理的训练，模型的过拟合程度最低；COD反演模型得到了合理的训练，但模型训练数据中，只有2019年数据的代表性较好，其他两年因训练数据代表性较差的原因，模型有较明显的过拟合现象；NH₃-N反演模型训练得最差，特别是对于训练数据代表性最差的2017年，用样本中测试数据反演得到NH₃-N参数值的统计特征，与样本中训练数据和所有地面监测数据的统计特征有明显的差异。

5.2.2　时间未匹配实验结果

图5为3种水质参数时间未匹配实验结果的箱形图。图中黑色、灰色，和白色箱形，分别表示实验中所有地面水质监测数据、时间未匹配的水质监测数据，和模型用时间未匹配的多源遥感及气象数据反演得到的水质参数值的统计特征。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 所有地面水质监测数据，时间不匹配的水质监测数据、模型用时间不匹配的多源遥感及气象数据反演的水质参数的年度对比箱形图

（黑色箱形代表所有站点的地面水质监测数据，灰色箱形代表时间不匹配的水质监测数据，白色箱形代表模型用时间不匹配的多源遥感及气象数据反演出的水质参数；箱形图中，箱体内的实线表示中位数，箱体的下边缘和上边缘分别表示下四分位数（Q1）和上四分位数（Q3），箱形外的上下横线表示最大值和最小值，o表示异常值）

Fig.5 Annual boxplots of all ground-based water quality monitoring data， time-mismatched water quality monitoring data， water quality parameters inverted by the models with time-mismatched multi-source remote sensing data

根据表3，时间未匹配的地面水质监测数据量约占到所有地面水质监测数据总量的29%左右，因此与图4相比，灰色箱形与黑色箱形更为接近，表明时间未匹配的水质监测数据的代表性，比训练模型所用数据的代表性更好。

从图5可以看出，DO和COD反演模型都有较好的时间鲁棒性，NH₃-N反演模型的时间鲁棒性很差，这可能是因为图4中所示的NH₃-N反演模型训练数据的代表性较差等原因导致。

5.2.3　空间未匹配实验结果

图6为3种水质参数空间未匹配实验结果的箱形图。图中黑色、灰色，和白色箱形，分别表示实验中所有地面水质监测数据、空间未匹配的水质监测数据，和模型用空间未匹配的多源遥感及气象数据反演得到的水质参数值的统计特征。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 所有地面水质监测数据，空间不匹配的水质监测数据、模型用空间不匹配的多源遥感及气象数据反演的水质参数的年度对比箱形图

（黑色箱形代表所有站点的地面水质监测数据，灰色箱形代表空间不匹配的水质监测数据，白色箱形代表模型用空间不匹配的多源遥感及气象数据反演出的水质参数；箱形图中，箱体内的实线表示中位数，箱体的下边缘和上边缘分别表示下四分位数（Q1）和上四分位数（Q3），箱形外的上下横线表示最大值和最小值，o表示异常值）

Fig.6 Annual boxplots of all ground-based water quality monitoring data， space-mismatched water quality monitoring data， water quality parameters inverted by the models with space-mismatched multi-source remote sensing data

根据表3，空间未匹配的地面水质监测数据量约占到所有地面水质监测数据总量的71%左右，因此灰色箱形与黑色箱形大致相同，表明空间未匹配的地面水质监测数据与所有地面水质监测数据具有基本相同的统计特征。

从图6可以看出，DO反演模型都有较好的空间鲁棒性，COD反演模型的空间鲁棒性较差，NH₃-N反演模型的空间鲁棒性最差，这同样可能是因为图4中所示的NH₃-N反演模型训练数据的代表性较差等原因导致。

图6所示的DO和COD反演模型的空间鲁棒性，比图5所示的两个模型的时间鲁棒性要差。图6和图5的唯一差异是，白色箱形所代表的水质参数反演所用的多源遥感及气象数据的筛选条件不同。图6所用的空间未匹配多源遥感及气象数据的筛选条件是，对应某个地面水质监测站的所有多源遥感及气象数据中，没有任何一组被选入训练和测试模型用的样本数据集。图5所用的时间未匹配多源遥感及气象数据的筛选条件则为，至少有一组被选入训练和测试模型用的样本数据集。由此可见，是否能有更多的地面站点对应的监测数据，被选入模型训练和测试的样本数据集，对模型的时空鲁棒性有直接的影响。

6 讨论

6.1　地表温度和植被指数对反演性能的影响

在反演DO、COD和NH₃-N时，温度和植被指数是两个重要因素。首先，温度会通过直接或间接的方式影响DO、COD和NH₃-N的浓度。对于所研究的3种水质参数，DO主要来源于从大气中转移到水体里的氧气，以及藻类和其他水生植物光合作用产生的氧气^［41］。COD作为衡量有机污染物对水污染程度的指标^［42］，会受到微生物消耗氧气以降解有机物这一过程的影响^［11］。而NH₃-N会在氧气的参与下被藻类吸收^［43］。所以水温会通过影响氧气溶解度的方式^［44］，进而影响DO、COD和NH₃-N的浓度。同时水温也会影响微生物的呼吸和氧化过程^{［45， 46］}，导致DO、COD和NH₃-N浓度变化。因此在反演水质参数时考虑温度的影响，能够在一定程度上提高模型的精度。

而NH₃-N作为营养物质，会导致藻类和其他水生植物过度繁殖^［47］，大量消耗水中的氧气，DO与COD的浓度也会随之变化。叶绿素a作为浮游植物中最重要的色素之一，可以衡量浮游植物生物量和水体富营养化程度^［13］。而植被指数与水体中的叶绿素a浓度有较高的相关性^{［45，48-49］}，另外河流周围植被的活跃程度（如森林的覆盖程度），也和水中的叶绿素a浓度具有一定关系^［13］。因此，在建立反演模型时，增加反映水体营养状态的植被指数等参数，可以提高模型的反演性能。

6.2　气溶胶和风速对反演性能的影响

大气中的各种由固态和液态颗粒物组成的气溶胶，会通过直接和间接的吸收和散射方式影响光在大气中的传输过程。卫星载荷的传感器接收到的水体表面出水辐射，不仅受水体成分和水底吸光度的影响，也受到大气中气溶胶散射导致的衰减影响^{［7， 50］}。在遥感数据产品生产过程中，通常用均值或某一大气剖面状态预先定义的气溶胶光学厚度来进行大气辐射校正^{［7， 51］}。这种方式并不能反映高频的气溶胶光学厚度的实际变化情况。因此在反演模型中，使用MODIS的气溶胶数据产品，可以一定程度上弥补这个问题，从而提高模型反演性能。

河流上的风速会影响水面的平整度，进而对入、出水辐射造成影响。在反演模型的输入数据中加入风速数据，可以使模型在一定程度上考虑水面平整度的影响，从而改善模型的反演性能。

从表5中可以看到，将气溶胶光学厚度和风速信息加入到模型输入中时，DO、COD和NH₃-N的反演精度会有所提高。

6.3　模型的应用

以南通附近的长江段和东莞附近的珠江段为研究区域，在谷歌地球引擎平台上收集研究区域在特定时间的哨兵2、MODIS地表温度、植被指数和气溶胶光学厚度数据以及ERA5风速数据。基于目视制作的水体掩膜栅格数据，提取出遥感及气象数据中对应于水体部分的数据，作为多源遥感及气象数据。利用表5中能获得最佳反演性能的模型及对应的输入特征组合，反演出两个区域的DO、COD和NH₃-N水质参数，并如图7所示分别绘制它们的空间分布图。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 用多源遥感及气象数据反演的一部分长江河和珠江江段的DO，COD，NH3-H水质参数空间分布

Fig.7 Spatial distribution of DO， COD and NH3-H water quality parameters in a part of Yangtze and Pearl River retrieved from multi-source remote sensing data

图7（a）显示的南通附近长江段的DO空间分布具有较为明显的变化，水质由北向南逐渐变差。相比之下，COD和NH₃-N空间分布较为均匀。在该区域西南部的近岸水域水质较差，推测是因为此处为一个半封闭的湾区，地形的原因使水域的流通性变弱，污染物会在此聚集^［52］。

图7（b）显示的南通西侧长江段，DO和COD的浓度总体上从东向西逐渐增加，西部的双山岛附近COD浓度较高，推测可能是因为双山岛对岸分布着工业区和码头，工业活动和船舶活动导致附近水质较差。另外，COD和NH₃-N在东北部岛屿附近的长江支流处，均出现了较高值。该长江支流北岸分散着几处农田，而南岸的岛屿又分布着居民区、度假村、工业区等，因此推测该狭长地带两岸较为活跃的人类活动，影响了流经此处的长江支流的水质。

从图7（c）显示的东莞附近的珠江段可以看出，该流域东北部穿插在城市或者农田之间的一些支流中，COD和NH₃-N浓度偏高。可能是因附近生活排水和农业灌溉和地表径流导致的氮和磷化合物进入水体，过量的营养元素导致水体富营养化^［11］的影响。人类的日常活动、工业废水的排放和农业化肥的使用等，都会对DO、COD和NH₃-N的浓度造成影响。

7 结论

本文基于多源遥感卫星数据，利用全国主要流域地表水水质监测站观测数据，建立了基于机器学习的反演模型。其中选择对DO、COD和NH₃-N 3种非光学活性水质参数进行研究，并为每种水质参数选择最佳的输入特征组合。结果显示基于机器学习的模型比传统的多元线性回归模型更能捕捉到反演非光学活性参数所需要的多源遥感以及气象数据的线性和非线性关系。通过实验证明，当样本的代表性与全局水质数据较为接近时，反演模型的时间和空间鲁棒性都比较好。另外对输入特征组合进行分析，发现对于溶解氧的反演，温度和反映植被动态的植被指数较为重要；对于化学需氧量和氨氮，影响反射率的风速和气溶胶更能够提高模型的性能。在未来的工作中，计划将风速数据的降尺度操作进行更精细的模拟，以在本研究的基础上进一步提高水质参数的反演效果，为遥感水质监测和水质分析提供有用的借鉴信息。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

ZHOU

Guanhua

， TANG

Junwu

， TIAN

Guoliang

， et al.

Uncertainty analysis of inland water quality remote sensing： A renew

［J］. Advances in Earth Science，2009，24（2）：150-158.