上海城市水体富营养化关键参量的遥感反演与时序分析

李嘉皓　田波　曹芳　胡越凯　段元强　谢泽昊　彭亚　姜文浩　范惠芳

关键词：Sentinel-2;富营养化;叶绿素a;浊度;时空变化;上海

中图分类号：P407.1 文献标志码：A DOI：10.3969/j.issn.1000-5641.2022.01.015

0引言

城市是人口高度密集的地区，高强度的人类活动会引发水质下降、湿地退化等环境和生态问题[1].就河口城市而言，局地城市人类活动产生的富含营养物质的化肥农残等污染物排放进入水体，造成水体污染和富营养化程度不断加重[2].水体富营养化导致藻华暴发，进而引起水体含氧量下降、水质恶化，造成水生生态系统失衡，对城市人类生产生活和社会经济的可持续发展造成负面影响[3].2019年《中国环境状况公报》显示，长江口近岸海域水质多为V类水[4].因此，调查长江河口城市，例如上海城市水体的富营养化状况具有重要的科学意义和迫切的现实需求.

叶绿素a（Chlorophyll-a，Chl-a）浓度和水体透明度常作为反映水体富营养化程度的关键水质参量[5-6].研究表明，水体浊度、悬浮颗粒物浓度和透明度这3个变量密切相关，可用浊度替代透明度作为评价富营养化程度的关键参量[7].而Chl-a浓度与浮游植物生物量密切相关，其水平决定水体吸收和散射特性，是反映水体富营养化程度的另一重要参数[8].同时，水体浊度会改变水下光场和水生动植物的生存条件，从而影响水体浮游植物初级生产[1，9].此外，悬浮沉颗粒物具有较强吸附作用，被认为是污染物和重金属的重要载体，影响着水环境健康[10].

内陆水体的水质监测常通过传统野外现场采样分析获得[11]，而这一方法仅适用于小尺度相对稳定的水体环境，无法实现高频次大尺度的水体信息获取，不能很好满足对复杂和快速动态变化的城市水体实现整体评估的要求.近几十年来，卫星遥感技术逐步用于内陆水体水质的监测.Hou等[1]针对长江下游大型湖库，建立了基于MODIS卫星的波段比值算法，描绘了湖库水体的浊度分布特征;Guan等[12]联用MERIS与OLCI数据，评估了长江平原大型湖泊的Chl-a浓度与富营养化程度;黄启会等[13]使用比值模型估算了百花湖（106°27′E～106°34′E，26°35′N～26°42′N）的Chl-a浓度;徐祎凡等[14]使用环境1号卫星，选用Chl-a浓度监测太湖富营养化情况，结果表明与综合营养状态指数法相比误差较小.受到卫星分辨率的影响，上述研究主要集中在单一大型湖库或是单一类型水体.而相较于大型湖泊等内陆水体，河口城市水体单元更为微小，感潮河网影响更加强烈，更易受到城市人类活动影响，对利用传统的中低空间分辨率的卫星观测水体提出了挑战.

欧空局（European Space Agency，ESA）于2015年和2017年发射的Sentinel-2系列卫星，其数据具有高时间分辨率（3～5d/景）和高空间分辨率（≤20m）的特点，可以满足高动态变化的城市小面积水体研究，在城市内陆水体动态监测中具有重要应用价值和潜力[15].本文主要采用Sentinel-2地表反射率数据，结合地表水体实测光谱以及Chl-a浓度和浊度数据，研究适用于上海城市内陆水体的Chl-a浓度和浊度遥感算法，分析上海内陆水体Chl-a浓度和浊度的时空动态变化特征.

1研究区概况与数据获取

1.1研究区

上海位于120°52′E～122°12′E、30°40′N～31°53′N的长江流域下游三角洲地区，是中国东部特大河口城市，属亚热带湿润气候，四季分明，年平均温度为16.1℃，年降水量达1173.4mm.作为典型平原感潮河网地区，主要自然水系有长江、黄浦江、苏州河和淀山湖，河道众多[16].除自然水系外，其人工湿地主要包括库塘、养殖塘、输水河，数量庞大密集，具有河口城市典型感潮水体特征（图1）[17].

1.2 野外现场数据测量

2019年10月—2020年8月，对上海市典型内陆水域开展了4个季节的Chl-a浓度、浊度和透明度测量，并现场同步采集水体遥感反射率（Rrs）.利用定标后的Manta2多参水质仪测量Chl-a浓度和浊度;利用20cm的黑白透明度盘测量水体透明度;利用ASD-HH2手持式便携光谱仪（光谱范围为350～1075nm）实测水体光谱.在北京时间10：00—15：00共采集89组水质及光谱参考样本，采样日期分别为：2019年10月25日（5组）、2019年12月26日（4组）、2019年12月30日（9组）、2020年5月11—12日（24组）、2020年7月30日—8月5日（31组）、2020年8月15—17日（16组）;采样类别分别为：永久性河流8组、永久性淡水湖7组、草本沼泽5组、库塘25组、运河/输水河33组、水产养殖场/种植塘11组（图1）.

现场光谱测量采用水面以上测量法[18]进行，避免水体镜面反射以及水面泡沫等干扰，测量参数包括水体上行辐射亮度（L），天空漫射光辐射亮度（L）和标准灰板辐射亮度（L），并计算水体反射率R[18].其公式为

式（1）中：ρ为标准参考灰板的反射率，约等于20%;r为水面的菲涅耳反射，根据NASA海洋光学协议，平静水面假定使用理论值0.022，在5m/s左右风速下，r可以取0.025[19].

1.3遥感数据获取与预处理

本文使用2018年12月—2019年12月的Sentinel-2地表反射率数据（https：//earthengine.google.com），数据来源于ESASentinel-2星座的S2-A和S2-B兩个卫星，卫星搭载有多光谱仪（MultiSpectralInstrument，MSI），可观测从可见光到近红外的13个光谱带，空间分辨率为10m、20m和60m[15].主要利用GoogleEarthEngine（GEE）云平台和ENVI5.3软件数据获取调用及处理.利用GEE平台，调用几何精校正（UTM/WGS84投影）和辐射定标处理的Sentinel-2（L2A）数据，该数据已进行Sen2Cor大气校正处理.虽然该大气校正方法可用于内陆水体的研究，但仍存在一定偏差[20-21].考虑到河口城市内陆水体流动性较强、变化较快的特点，本研究采用卫星过境同步现场实测反射率数据（前后时间偏差1h以内）对Sentinel-2遥感反射率进行校准.具体方法如下：

使用Sentinel-2的光谱响应函数对实测光谱进行等效光谱模拟，转换为相应波段的等效遥感反射率.本文选用10m空间分辨率的B2（443nm）、B3（560nm）、B4（665nm）、B8（842nm）和20m分辨率的B5（705nm）、B6（740nm）、B7（783nm）、B8A（865nm）共8個波段进行反演（图2（a））.结果表明，相较于实测获得的水色遥感反射率，基于Sen2Cor的大气校正后的Sentinel-2的地表反射率存在明显高估现象，其主要原因是Sen2Cor处理器的最初设计目标为陆地，水体区域会因为日照产生太阳反射.本文参考Hedley等[22]提出的方法，使用近红外波段对Sentinel-2数据进行太阳反射校准，公式为

式（2）中：R为校准后波段i的反射率，R′为校准前波段i的反射率，R为校准前近红外波段的反射率，b为常数.使用5个同步观测点位作为测试数据，对b赋值，计算校准反射率与实测反射率之间的均方根误差（RMSE），当b=0.61时，其RMSE最小，为0.0037sr.通过上述公式对其余8个同步观测点的太阳反射进行校准，并与地面实测光谱反射率进行对比分析（图2（b）），结果表明：经校正后的RMSE为0.0028sr，优于Sòria-Perpinyà等人仅用Sen2Cor处理器的校正结果（0.024sr）[21].

影像质量控制除了大气校正之外还完成了去云处理，主要利用GEE平台，根据Sentinel-2影像QA60波段标记去云.

2研究方法

2.1城市水体提取与光谱分析

城市水体空间信息的提取是遥感反演水体水质参数的基础.由于水体在近红外（NIR）波段存在强吸收现象，本文选取3种常用水体指数对水域进行提取，分别为：归一化水体指数NDWI，改进的归一化水体指数mNDWI，城镇水体指数CIWI[23].使用大津法（otsu）对水体指数运算结果进行自动阈值分割[24]，以区分水体和非水体.使用随机生成的1000个验证点计算提取精度（表1），结果表明，mNDWI水体提取精度较高.

对89组站点实测光谱、Chl-a浓度及浊度的统计分析表明（图3），Chl-a浓度的观测值位于0.76～161.69μg/L，平均值为13.75μg/L;浊度观测值位于4.1～122.30NTU，平均值为42.50NTU.

通过分析不同Chl-a浓度和浊度的光谱发现（图3），在675nm附近，由于Chl-a吸收，出现了反射率的吸收谷.710nm附近出现叶绿素a的荧光峰，其位置大小会受到Chl-a浓度的升高向长波段方向发生位移.Chl-a浓度越高，其反射峰也会相应增加[13]（图3（a））.对比图3（b）中浊度与光谱的关系可以发现，740nm附近形成小的吸收谷，810nm附近有小反射峰，其大小均会随着浊度增加而增加，这主要由悬浮物散射造成[1].

利用光谱响应函数将各站点实测水体遥感反射率转换为等效波段遥感反射率，与实测Chl-a浓度、浊度构建反演模型.将实测89组数据随机分为两组，其中60组数据用于建模，29组数据用于评价模型精度.最后利用卫星同步过境前后1h现场测得的Chl-a浓度与浊度数据（共8组）分别和卫星反演获得的数据进行统计分析，以评估反演模型的稳定性.

2.2Chl-a浓度反演及验证

Chl-a在蓝光波段中具有高吸收性，通常可用蓝绿波段比反演Chl-a浓度[25].河口城市内陆二类水体中，悬浮物以及黄色物质（CDOM）对较短波段反射率造成干扰[26].因此，反演二类水体中的Chl-a浓度需选择长波（大于600nm）波段或使用波段组合形式，以降低悬浮物和CDOM的影响[26].

现有的内陆水体Chl-a浓度反演算法主要有经验算法[26-28]、半经验算法[29-30]和解析算法[31].利用同步实测等效遥感反射率与实测Chl-a浓度测试了几种模型（表2），并使用确定系数（R2）和均方根误差（RMSE）检验了这些模型的性能.

结果表明，665nm（B4）、705nm（B5）和740nm（B6）这3个波段建立反演算法能较好地反演Chl-a浓度（图4）.其主要原因是，在长波段665nm与705nm附近分别为Chl-a的吸收谷和荧光峰，颗粒物和CDOM吸收较低（图3），能较好地将颗粒物和CDOM的影响去除.近红外波段740nm反射率主要是纯水吸收作用，可将水体组分中后向散射影响去除[30].

使用实测光谱数据进行验证（图4（e）），R为0.96，RMSE为4.33μg/L.通过实测Chl-a浓度与星地同步观测对比（图5（a）），具有较好的拟合性能（R=0.79，RMSE=0.86μg/L）.

2.3 浊度反演及验证

数据分析表明，水体透明度与浊度呈良好的对数关系（R2为0.9011），由于透明度的测量存在一定的主观性，本项研究采用浊度表征水体透明度.水体浊度是3种光学成分（叶绿素、悬浮物和有色可溶性有机物）的综合表现[1].基于实测光谱和同步观测浊度数据，已有多种浊度遥感反演算法，包括经验算法、半经验算法、半解析算法和物理算法[1，28，32-34].在不同研究区域，由于水体悬浮物粒径、密度、浮游植物数量或有色可溶性有机物等的差异[35]，需要选择区域内最优波段或波段组合算法（表3）.

根据光谱测量信息（图3）分析，长波段反射率与浊度相关性较强，R均达到0.8以上，本文根据Petus等[32]的算法提出了一种以Sentinel-2-B7波段（783nm）的浊度反演算法.实测光谱数据验证表明（图6（f）），该模型R为0.87，RMSE为12.04NTU.实测浊度与星地同步观测对比（图5（b））也显示了较好的反演结果（R=0.91，RMSE=8.18NTU）.

3结果与讨论

3.1上海城市水体Chl-a浓度的时空分布特征

通过计算2018年12月—2019年11月每个季度中位数值的Chl-a浓度，将4个季度的均值进一步平均来估计年度Chl-a的平均值.根據年平均Chl-a浓度（图7（a）），上海内陆水体Chl-a浓度有50%左右小于等于10μg/L，处于低Chl-a水平;15%位于10～20μg/L，处于中等Chl-a水平;7.4%位于20～30μg/L，处于高Chl-a水平;23.28%大于30μg/L，处于极高Chl-a水平.上海市不同类型水体的Chl-a浓度存在着明显差异，具体而言，养殖场/种植塘具有最高的Chl-a水平，高值区达到200μg/L，其他类型水体Chl-a浓度较低，其中永久性淡水湖和库塘水体Chl-a浓度为8μg/L左右，永久性河流和运河/输水河具有最低的Chl-a浓度，大都在6μg/L左右.

上海内陆水体Chl-a浓度变化呈现一定的季节变化特征（图8）.夏季（6—8月）呈现了最高的Chl-a浓度，春季（3—5月）和秋季（9—11月）次之，冬季（12—次年2月）的Chl-a浓度最低，这一结果与Guan等[12]在长江平原大型湖泊得出的结果相似.

3.2 城市水体浊度的时空分布特征

浊度反演结果表明，在2019年（图7（b）），4.56%的内陆水体浊度小于等于20NTU，45.05%位于20～40NTU，处于较低浊度水平;32.8%位于40～60NTU，处于中等浊度水平;12.51%位于60～80NTU，处于较高浊度水平;其余5.08%的水体大于80NTU，处于高浊度水平.从水体类型比较看，养殖场/种植塘具有最高的浊度水平，高值区达到130NTU;永久性河流及运河/输水河由于地处河口区，受感潮河网影响，浊度较高，约为45NTU;永久性淡水湖水体的浊度为30NTU;库塘最低，为20NTU.

季节尺度上（图9），不同水体均在夏季呈现了最高浊度，春季次之，秋冬季浊度最低.主要由于冬季水体温度较低，藻类代谢程度低，同时降水量/径流量低，水体面源污染少.而春夏季节降水增多/径流增大，面源污染增大，水体温度升高，导致藻类生长旺盛，浊度增大[36].

3.3城市水体Chl-a浓度与浊度相关性研究

从上述分析可知，Chl-a浓度和浊度均具有强烈的空间和时间异质性，但具有相似的时空变化规律.对5种类型水体分别采样，影像去云处理后计算采样点的中位数浊度和Chl-a浓度值，使用移动平均法排除极端数据，拟合时间序列曲线分布（2019年），计算浊度和Chl-a浓度的相关系数R（图10），分析不同类型水体Chl-a浓度和浊度要素的时间序列相关性.

对于受人类活动影响小的永久性河流、库塘和运河/输水河3类水体，其浊度和Chl-a浓度具有较强的相关性，R均在0.75左右（p<0.05）.永久性淡水湖（淀山湖）的相关性达到了0.67（p<0.05）.养殖塘的浊度和Chl-a浓度的相关性较弱，主要是养殖过程中的增氧、饲料投喂、消毒等人类活动过程对水体光学特性及水质产生影响[37].

3.4不确定性分析

在本研究中，使用经校正后的Sentinel-2地表反射率数据来估计上海市水体的Chl-a浓度和浊度这两个富营养化关键参量，但研究也存在一些不确定性：

（1）尽管Sentinel-2的重访周期为3～5d，但由于云、雾覆盖的限制，导致部分数据不可用.经统计，2019年覆盖研究区影像共有803景，云量小于30%的仅有187景，占比23.29%，从而导致时序趋势分析数据的不充分性和不确定性.因此，对于Chl-a浓度和浊度的年平均值，是根据相应的4个季度平均值而不是一年中的每日平均值估算的[12].

（2）Sentinel-2的水色大气校正方法常用的有Sen2Cor方法和ACOLITE方法等[21，38]，上海市水体的实测数据表明，两种方法都存在反射率偏高的现象.因此，本文提出了一种简单的大气校正后处理的方法，目的是在上海城市水体中获得较为可靠的遥感反射率.

（3）城市水体的破碎化导致了严重的水陆邻接效应，混合像元较多，从而导致反演的不确定性，遥感断面监测显示，应去除水边界的一个像元，排除土地的影响，尽可能降低反演的不确定性.除此之外，城市水体受人类活动影响强烈，水上的船只通行也会造成像元污染，故本研究使用各像元中位数值进行季度计算，使用均值进行年际计算，以最大化减少像元污染的概率.

4结论

本文建立了基于Sentinel-2的上海城市内陆水体富营养化关键参量的具有较高精度的遥感反演模型，利用这一模型，反演了从2018年12月到2019年12月的不同水体中Chl-a浓度和浊度，并对这两个参量的时空变化特征和相互关系进行了分析，结果表明：

（1）利用实测数据，基于Sentinel-2卫星数据，本文构建了水体Chl-a浓度和浊度的反演模型.其中，Chl-a浓度模型的R为0.96，RMSE为4.33μg/L;浊度模型的R为0.87，RMSE为12.04NTU.反演精度较高.

（2）上海城市不同类别水体的Chl-a浓度具有较大时空变异特征，由高到低排序大致为：养殖场/种植塘>永久性淡水湖≈库塘>永久性河流≈运河/输水河.

（3）类似地，上海城市水体浊度呈现出较大的空间异质性，由高到低排序大致为：养殖场/种植塘>永久性河流≈运河/输水河>永久性淡水湖>库塘.

（4）水体Chl-a浓度和浊度存在较强的季节变化，表现出夏季高冬季低的现象.

（5）在受人类影响较小的水体类型中（如运河，库塘，永久性河流），Chl-a浓度和浊度存在较强的显著正相关关系，而在受人类活动影响较剧烈的水体类型中，如养殖水体中，两者相关性较弱.

河口城市水体受流域、海洋和局地人类活动影响强烈，其水体富营养化变化相对快速复杂，时空季节性分异性强，与人类活动类型、方法和强度密切相关，基于时序遥感反演和分析Chl-a浓度和浊度变化需进一步深入研究.

（责任编辑：李万会）