“山水林田湖草沙”是生命共同体, 它们是相互独立而又彼此依存的自然生态系统。湖泊在生态环境中扮演着重要角色, 在防洪抗旱、水资源供给、维持生态系统平衡以及减缓气候变化等方面都发挥着不可替代的作用(马荣华等, 2011)。近年来, 随着中国湖泊污染物削减和退耕还林等生态环境保护工程的实施, 典型湖泊的富营养化状况得到改善(Liu et al, 2020)。然而, 各地区湖泊水质时空变化在过去几十年中表现出一定的复杂性, 这些变化在不同地区湖泊之间存在一定的差异, 一些湖泊富营养化得到遏制的同时, 另一些湖泊水质则有所恶化(张运林等, 2022)。生态环境修复效果具有延迟效应, 需要进行长期的动态监测和评价(Duan et al, 2022), 加深对中国湖泊水质空间分布状态的认识, 探索湖泊水质状况变化趋势, 对维护现有治理成效和进一步推动湖泊生态改善具有重要意义。

叶绿素a (chlorophyll a, chl a), 悬浮颗粒物和可溶性有机物被称为水色三要素(Morel et al, 1977), 叶绿素a是浮游植物生物体的重要组成成分之一, 是反映湖泊营养状态的关键指标, 当其浓度偏高时容易引起水华暴发, 进而破坏湖泊生态环境(Pepe et al, 2001)。通过建立湖泊水体中叶绿素a含量反演模型, 利用遥感影像对湖泊水质进行评估, 一直是湖泊水色遥感的研究重点, 目前已经有大量涉及湖泊水色参数定量反演的相关研究。但过去人们往往受制于传统方法的局限性, 侧重于对典型湖泊或特定地区湖泊的深入研究, 例如, 陶星宇等(2023)基于经验分析法, 结合叶绿素a实测数据与MODIS影像波段反射率算法模型, 对西藏地区典型湖泊的叶绿素a浓度进行反演, 分析了该区域湖泊叶绿素a浓度在不同季节的时空变化特征。杜石霞等(2018)利用遥感影像研究了近40年抚仙湖岸线和面积的时空变化特征, 并对其驱动因素进行了分析。典型湖泊或局部地区的研究不足以实现对中国各区域湖泊的全面监测和评估, 在这一背景下, 利用遥感方法对中国范围内的湖泊水质参数进行反演, 探明近几十年中国湖泊叶绿素a浓度的空间分布和变化规律变得尤为重要。

随着Google Earth Engine和PIE-Engine Studio等遥感大数据计算平台的出现和发展, 为地理空间大数据分析提供了有效工具, 近年来被广泛用于处理与环境变化相关的工作。Dai等(2023)基于Landsat系列遥感卫星影像, 依据国际照明委员会(Commission Internationale De L’Eclairage, CIE)色彩空间开发了藻华检测算法, 制作了1982~2019年间的全球水华暴发数据集。Hou等(2023)利用Sentinel数据集对2019~ 2021年间中国内陆湖泊藻华进行了识别提取, 分析了藻华暴发的纬度梯度变化。Zhao等(2024)使用Sentinel-2 MSI Level-2A反射率影像数据集反演了2019~2021年间全球湖泊叶绿素a浓度, 对其时空分布及其驱动因素进行了分析。现有研究更多侧重于对水体浮游植物生物量的提取, 例如藻华的识别, 然而, 对于内陆大范围湖泊群的叶绿素a浓度评价时间周期较短, 难以反映长时间变化规律, 并且陆地不同区域和不同流域之间湖泊叶绿素a浓度的时空分布变化状态研究较为鲜见。针对以上研究进展与现状, 本文聚焦在中国内陆面积大于10 km²湖泊叶绿素a浓度的时空分布变化特征。为了捕捉更为细致的时空变化规律, 本研究拟采用具备较高分辨率和较长时间跨度的Landsat影像数据集进行反演, 揭示中国湖泊叶绿素a浓度的长期趋势和周期性变化。同时, 流域作为湖泊生态系统的自然边界, 其内部物质与能量的产生, 输移和转化密切相关, 影响着湖泊的水文环境和水体特征, 以流域为单位进行研究有助于分析湖泊水质变化规律的形成机制(冷疏影等, 2003)。因此, 为深入理解中国湖泊与区域生态环境之间的关系, 本文将依据分区域和分流域两种方式进行分析, 探究湖泊叶绿素a浓度在不同经纬度上和不同流域内的特征。这种研究方法不仅有助于填补过去研究的空白, 而且能够促进湖泊水质反演模型的改进, 同时, 对于及时监测水质变化, 发现潜在环境问题和提高水资源管理水平具有积极意义。

基于对中国湖泊生态环境改善, 长时间监测管理以及湖泊与流域相互关系研究等方面的需求, 本文将利用Landsat系列卫星等相关数据集, 对近40年在中国范围内面积大于10 km²的湖泊叶绿素a浓度进行反演, 并且从以下几个角度进行时空分异特征分析: (1) 湖泊叶绿素a浓度的年均值空间分布和不同季节空间分布状态; (2) 湖泊叶绿素a浓度随经纬度的空间变化特征; (3) 根据时间序列数据分析湖泊叶绿素a浓度的年际变化趋势; (4) 以流域为单位分析各区域湖泊叶绿素a浓度的演变趋势。

1 研究区域与数据源 1.1 研究区域

中国湖泊资源丰富, 广泛分布于长江流域, 黄河流域和松辽河流域等地区, 其中, 以青藏高原地区的湖泊贮水量最为丰富, 其次为东部平原地区, 拥有湖泊数量最多的地区为西藏自治区, 内蒙古自治区以及黑龙江省(刘吉峰等, 2008)。各区域湖泊所处的地貌类型和气候特征不同, 这些环境因素对湖泊水文环境和水体特征产生了多样化的影响。不同海拔和气候区的气温, 光照和降水等条件存在差异, 气温影响湖泊温度, 从而影响水体的物理和化学过程, 如水的密度、流动性、氧气溶解能力等(Wang et al, 2023)。光照强度影响水体中浮游植物的光合作用和生物的活动。降水量则会影响湖泊水位, 可能导致洪涝或枯水, 进而影响湖泊生态系统(Li et al, 2020)。中国湖泊涵盖了淡水湖泊、咸水湖泊和盐水湖泊等多种类型, 淡水湖和浅水湖泊主要分布于东部地区, 咸水湖主要分布在蒙新高原地区和青藏高原地区(Qin et al, 2023)。在人口密集、经济发达的东部平原地区, 虽然水网稠密, 但大多数为浅水湖, 环境容量有限。由于气候变化和工农业污染等原因, 这些湖泊的水质也受到了威胁。全球变暖导致的气温升高和降水模式变化等影响着湖泊生态系统的稳定性(Naithani et al, 2011)。工农业污染引起的污染物排放, 如化学物质和重金属元素, 导致部分湖泊水质恶化, 产生富营养化现象。这种富营养化导致藻类过度生长形成藻华, 对湖泊的生态平衡产生不利影响(蔡阳等, 2019)。

本文利用中国湖泊矢量数据集一共筛选了面积在10 km²以上的湖泊912个。依据中国科学院地理科学与资源研究所中国流域矢量数据, 将研究区划分为九大流域片区, 各流域内的湖泊在不同面积和海拔高度区间数量如表 1所示。

1.2 数据源

本文使用的中国区域范围矢量数据采用全国地理信息资源目录服务系统网站(https://www.webmap.cn/)1︰100万基础地理信息数据2021版, 该网站主管部门为中国自然资源部。中国省级行政区划边界数据和流域片分区数据来源于资源环境科学数据注册与出版系统(https://www.resdc.cn/), 该网站数据由中国科学院地理科学与资源研究所提供。

研究所需的中国范围内湖泊矢量数据来源于HydroLAKES全球湖泊数据库(Messager et al, 2016), 该数据库由美国宾夕法尼亚州立大学与瑞士联邦水资源和水力工程研究所联合开发, 共收集整理了全球超过18 000个湖泊的相关数据, 包括湖泊的地理位置、面积、最大深度、平均深度和水位变化等多种参数。

本研究采用Google Earth Engine平台提供的Landsat系列卫星遥感数据集, Landsat系列卫星是美国航空航天局与美国地质调查局合作开发的遥感卫星, 提供了较长的历史连续空基观测记录, 可以满足研究目标的时间跨度需求。本文选择使用Landsat TM, ETM+和OLI三种传感器的影像数据集, 并根据区域和时间范围要求进行过滤。针对数据集中存在的缺失部分进行补全, 并进行影像裁切和云掩膜等预处理流程, 以获取符合研究要求且具有较高质量的遥感影像数据。一致性是不同遥感载荷对同一目标反射率在物理意义和数值上的统一性与相似性(马灵玲等, 2023), Zhang等(2021)在基于Landsat系列卫星构建湖泊透明度遥感估算模型时, 对以上涉及到的TM, ETM+和OLI三种不同传感器数据源进行了评估, 估算结果具有高度一致性, 因此, 在分析过程中由于不同传感器带来的差异可以忽略。本研究选择的遥感影像数据均为Google Earth Engine平台Collection 2中的Level 2级Tier 1产品。这些数据产品经过几何校正, 辐射定标和大气校正处理, 相较于Collection 1数据, 精度有所提高。同时, 经过处理后的Tier 1产品具有更高的数据质量, 适合用于时间序列分析。

在ArcMap中利用中国范围边界矢量对HydroLAKES全球湖泊矢量进行裁切, 得到中国范围内所有湖泊的矢量数据, 剔除数据中面积小于10 km²的湖泊后, 对符合要求的湖泊矢量建立数据库。利用流域矢量对中国范围湖泊矢量进行划分, 得到九大流域片区湖泊矢量数据。

利用上述各流域湖泊矢量数据对Landsat系列遥感数据集进行范围过滤并掩膜, 得到研究区内所有湖泊遥感影像。根据研究目标对数据集进行1986~ 2022年时间范围筛选和年际季节划分(春季: 3~5月、夏季: 6~8月、秋季: 9~11月、冬季: 12~2月), 然后过滤出影像数据集中云层含量小于70%天气状况较好的部分, 利用影像中的质量信息波段(“QA_PIXEL”和“QA_RADSAT”波段)识别并且去除云以及阴影区域像元, 最终获取准确的地表反射率值。

2 研究方法 2.1 反演算法

在水质遥感研究领域, 建立叶绿素a浓度反演模型常使用遥感波段组合计算方法。在湖泊水体中, 可以利用归一化差值植被指数(normalized difference vegetation index, NDVI)反演湖泊水体中叶绿素a的浓度, 进而评估水质状况(Rouse et al, 1974)。NDVI的计算公式为

其中, I_NDV为归一化差值植被指数, ρ_NIR和ρ_RED分别表示近红外波段和红光波段的反射率。I_NDV值越高时, 表示该区域的叶绿素a浓度越高。

归一化差值植被指数在水体中的反演与其在陆地上类似, 在建立波段反射率与湖泊叶绿素a浓度的关系时, 将归一化差值植被指数作为自变量与实测数据进行线性回归得到湖泊叶绿素a浓度的反演算法模型。本文中使用的叶绿素a反演模型是基于王一旭(2015)在中国范围内的典型湖泊太湖进行水体采样后, 根据实测数据与Landsat遥感影像波段组合反射率进行回归分析, 建立了NDVI与叶绿素a浓度之间的关系算法, 实现了对湖泊水体中叶绿素a浓度的反演。该叶绿素a浓度反演模型为

其中, chl a表示叶绿素a的浓度, 单位为μg/L。在该计算公式中, 归一化差值植被指数与叶绿素a浓度值呈线性正相关关系, 所在区域影像波段计算的归一化插值植被指数值越大时, 相应区域的叶绿素a浓度也越高。根据该模型反演得到的数值与实测值具有良好的线性关系, 决定系数R²为0.74, 具有较好的反演精度。

2.2 变化趋势分析

一元线性回归模型可以模拟空间上湖泊叶绿素a浓度的年际变化趋势, 方程斜率系数计算方法使用最小二乘法。在Python中对叶绿素a浓度年际序列数据缺失值和异常值等进行处理, 调用Scipy库中的Linregress函数对湖泊叶绿素a浓度年际序列变化趋势进行一元线性回归Slope分析。该回归方程斜率的公式为

其中, S为线性回归方程的斜率, n为累计年数, 在本研究中n为38, 为第i年的叶绿素a浓度。叶绿素a浓度时间变化趋势显著性检验使用t检验。显著性检验方法是计算回归模型的P值, P值表示在零假设成立的情况下, 观察到的斜率或更极端情况出现的概率, 如果P值小于显著性水平(P值设为0.05), 则拒绝零假设, 认为观测到的变化显著。

2.3 稳定性评价

采用变异系数(coefficient of variation, CV)来比较不同区域湖泊叶绿素a浓度的波动程度。变异系数可以衡量数据相对于自身均值的离散程度, 在湖泊叶绿素a浓度时间变化稳定性分析中, 时间序列数据分别按照湖泊和不同流域进行分组, 然后计算其变异系数。变异系数由标准差与平均值之比得到。

标准差的计算公式为

变异系数的计算公式为

其中, c_v为变异系数, σ为标准差; 为平均值。变异系数可以消除不同数据数量级的影响, 适合用来比较不同数据集合的离散程度。变异系数越小, 表示数据越稳定。变异系数可以很好地反映空间数据在时间序列上变化的差异程度, 进而评价湖泊叶绿素a浓度时间序列数据的稳定性。

2.4 研究区分区

中国地势西高东低, 呈阶梯状分布, 不同区域的海拔高度和气候环境差异明显。为了控制海拔差异等因素带来的影响, 按照表 2中的划分原则, 将中国湖泊划分成四个分区。以第二阶梯的昆仑山脉和祁连山脉为分界线, 沿35°N进行南北方向划分, 将研究区划分为a和b两个分区, 该分界线也大致位于中国秦岭淮河南北地理分界线, 以此分析叶绿素a浓度在经度方向上的变化规律。再以横断山脉为分界线, 沿100°E进行东西方向划分, 将研究区划分为c和d两个分区。以此分析叶绿素a浓度在纬度方向上的变化规律。

3 结果与分析 3.1 湖泊叶绿素a浓度空间分布特征

利用中国地区湖泊矢量和Landsat系列遥感影像数据, 结合叶绿素a算法模型进行反演, 得到1986~ 2022年中国湖泊叶绿素a浓度的时间序列数据。图 1展示了2022年研究区内所有湖泊叶绿素a浓度年均值空间分布状态。可以观察到, 叶绿素a浓度年均值范围大致在5~30 μg/L之间, 浓度分布明显不均衡, 总体呈现东高西低的态势。

以胡焕庸线为界限, 西北侧的各湖泊叶绿素a浓度普遍偏低, 该区域北部蒙新高原湖泊分布较为分散, 叶绿素a浓度相对较高。南部青藏高原地区湖泊数量较多, 但湖泊大部分处于非淡水状态(Liu et al, 2021), 营养化程度低, 分布较为均匀。而胡焕庸线东南侧的湖泊分布密集, 叶绿素a浓度则普遍偏高, 并且自北向南呈现递增趋势, 南部沿海地区叶绿素a浓度最高。以上现象说明湖泊叶绿素a浓度受到地形地貌与气候环境的影响存在明显的地域性差异。

图 1中的饼状图统计了胡焕庸线东、西两侧叶绿素a浓度各梯度区间所占比例, 可见在胡焕庸线西侧大部分湖泊叶绿素a浓度均值处于15~20 μg/L之间, 约占所有湖泊的81.9%。在胡焕庸线以东所有湖泊中, 叶绿素a浓度在15~20 μg/L区间的湖泊依然占据大多数, 但相较于西侧有所减少, 约占总量的58.7%, 其次为20~25 μg/L区间, 约占26.8%。对应营养物质叶绿素a含量分级(表 3), 将湖泊叶绿素a浓度控制标准分为六个等级(高阳俊等, 2011)。在所有湖泊中大约69%处于轻富营养化程度, 其次为中富营养化状态, 约占17%, 在胡焕庸线东侧的中富营养化状态湖泊明显多于西侧。由此可见, 胡焕庸线两侧湖泊叶绿素a浓度空间分布差异明显, 东侧的叶绿素a浓度普遍高于西侧。

利用Loess局部加权回归平滑方法对每个分区均值进行拟合与平滑处理(经过测试每次平滑计算使用原始数据的13%进行加权平滑拟合效果最佳), 生成叶绿素a浓度随经纬度变化的趋势曲线, 如图 2所示。

图 2a为35°N以北的研究区a, 该区域大部分处于温带气候区。总体上看, 由西向东随着经度变化, 叶绿素a浓度呈现上升趋势。从局部上看, 由最西端到82°E, 叶绿素a浓度均值由20 μg/L持续下降为15 μg/L左右, 82°~90°E也是该区域叶绿素a浓度最低范围, 浓度范围为10~20 μg/L, 均值在15 μg/L附近浮动。82°E以西范围内湖泊分布较少, 并且大部分处于蒙新高原西侧, 而在82°~90°E范围湖泊数量较多, 多处于高原气候区, 叶绿素a浓度较低。90°~115°E之间, 叶绿素a浓度均值持续升高至20 μg/L, 该区域湖泊分布较少, 由西北向华北和华东地区过渡, 气候类型也由高原气候向中温带和南温带气候区转变, 人类活动加剧, 叶绿素a浓度逐渐升高。在115°~ 125°E之间湖泊叶绿素a浓度均值较高, 大约在20 μg/L上下浮动, 125°E达到峰值, 浓度范围在10~35 μg/L, 均值为23 μg/L左右。该区域湖泊分布较多, 地势较为低平, 人口分布密集, 叶绿素a浓度最高的湖泊集中在东北地区和东部沿海地区。

图 2b为35°N以南的研究区b, 总体上看, 由西向东叶绿素a浓度呈现先下降, 后大幅上升, 再下降的趋势。从局部上看, 在80°~90°E之间, 叶绿素a浓度均值最低, 大约为15 μg/L, 湖泊位置集中分布于青藏高原气候区。由90°E向东, 叶绿素a浓度均值上升, 大约在107°E到达峰值25 μg/L, 该过渡区横跨横断山脉, 海拔高度开始大幅下降, 气温也随之上升, 处于高原气候区与亚热带气候区的交界区域。在107°E以东, 叶绿素a浓度普遍较高, 存在小幅度波动下降, 最高值在30 μg/L左右, 均值大约为20 μg/L。该地区湖泊数量较多, 分布广泛, 为温度较高的亚热带地区, 同时人类活动密集, 存在高强度工业生产, 对附近的湖泊生态环境造成一定影响。

图 2c为100°E以西的研究区c, 总体上看, 随着纬度上升, 叶绿素a浓度呈现先下降后上升趋势。在28°~35°N之间, 湖泊数量较多, 且叶绿素a浓度普遍较低, 平均浓度在15 μg/L左右波动, 该地区位于高原气候区, 海拔较高, 气温较低。在35°N以北的地区, 跨过昆仑山脉和祁连山脉, 由高原气候区向中温带气候区过渡, 气温开始上升, 叶绿素a浓度也随之升高, 平均浓度最高为20 μg/L。在40°N以北, 随着纬度上升, 气温降低, 并且人类活动影响有限, 叶绿素a浓度均值有所下降。

图 2d为100°E以东的研究区d, 相较于其他区域, 该区域内湖泊数量多且分布均匀, 气候由北温带向南亚热带过渡, 叶绿素a浓度普遍偏高, 由北向南呈现波动上升趋势。50°~33°N处在北温带向南温带过渡区域, 随着纬度降低, 叶绿素a浓度均值缓慢升高, 大致由18 μg/L上升至23 μg/L。在40°N附近出现较高的浓度均值, 该地区位于吉林省与黑龙江省西部交界地区, 浓度值大约为22 μg/L。35°~25°N之间为长江中下游平原, 处于亚热带区域, 湖泊分布较多, 气温较高, 人口分布密集, 湖泊富营养化时有发生, 叶绿素a浓度均值大约在23 μg/L附近浮动。随着纬度继续降低, 叶绿素a浓度略微下降。由此可见, 向赤道热带地区靠近, 并不会使得叶绿素a浓度持续增加, 温度升高有利于湖泊浮游植物光合作用, 浮游植物大量生长繁殖, 然而温度过高也会抑制一些水生植物持续生长(赵娜等, 2020), 叶绿素a浓度相应降低。

3.2 湖泊叶绿素a浓度季节分布状态

图 3展示了研究区内所有湖泊叶绿素a浓度在不同季节空间分布状态。可以观察到, 叶绿素a浓度也存在季节性变化特征。从整体分布来看, 绝大部分地区夏季是叶绿素a浓度最高季节, 均值约为18.52 μg/L, 而冬季是叶绿素a浓度最低季节, 均值约为17.53 μg/L。空间分布状态均为东高西低。这表明不同季节的气温和降水等因素会对湖泊叶绿素a浓度产生影响。

根据各个季节叶绿素a浓度梯度占比统计(图 4)可知, 随着季节更替, 各浓度区间占比呈现出有规律的变化, 在浓度小于20 μg/L区间内, 夏季数量最少, 约占65.8%, 继而是秋季和春季, 分别约占77%和84.57%。冬季数量最多, 约占88.83%。相反, 在浓度大于20 μg/L区间内, 冬季数量最少, 约占11.18%, 继而是春季和秋季, 分别约占15.44%和23%。夏季数量最多, 约占34.2%。可以发现, 随着季节更替, 湖泊叶绿素a浓度也呈现出规律的递增或递减变化, 在浓度大于20 μg/L区间, 气温越高, 降水量越大的季节占比也越多, 浓度小于20 μg/L区间则表现相反。这可能由于夏季降雨增加了地表径流, 使得湖泊营养盐输入较多。同时, 光照较强, 温度较高为藻类的繁殖提供了有利条件, 导致叶绿素a浓度相对较高。秋季和春季的气温适宜, 光照适中的条件下, 湖泊藻类的生长相对较好。而冬季的日照时间短, 水温较低以及营养盐输入较少等特点, 导致水体中的生物代谢减缓, 生长相对较慢, 喜高温藻类大量死亡, 浮游植物丰度和生物量均降低, 因此叶绿素a浓度也较低(祁梅等, 2020)。

不同的是在浓度小于15 μg/L区间内, 占比最多的季节在秋季, 约占28.56%, 其次是夏季和春季, 分别约占22.73%和17.61%。最少的在冬季, 约占11.49%。湖泊叶绿素a浓度小于15 μg/L区间集中分布在青藏高原地区, 并且该区域湖泊叶绿素a浓度在夏秋季较低, 在春冬季较高。由此可见, 温度升高时不一定总会引起叶绿素a浓度升高, 高原气候区湖泊叶绿素a浓度变化与低海拔气候区存在差异, 这可能与青藏高原湖泊水质受控于降雨和蒸发浓缩等因素有关(朱立平等, 2019)。由于青藏高原地区夏季降雨和冰川融水较多, 径流增加, 湖泊水体中的营养物质被稀释, 限制了浮游植物的生长, 使得叶绿素a浓度降低。其次, 相关研究表明, 浮游植物较少的湖泊随着表面温度升高, 叶绿素a浓度倾向于降低, 而浮游植物较多的湖泊则随着温度升高, 叶绿素a浓度也会随之上升。高原湖泊浮游植物生物量相对较少, 因此部分湖泊叶绿素a浓度呈现出夏秋季低, 春冬季高的特点(Kraemer et al, 2017)。

3.3 湖泊叶绿素a浓度时间变化趋势

结合研究区在1986~2022年湖泊叶绿素a浓度序列构建的箱式图(图 5)进行分析, 1992年和2011年附近叶绿素a浓度均值略高于其他年份, 叶绿素a浓度均值较低的年份为1990、1999、2013和2021年。图 5中水平线为每间隔10 a计算得到的叶绿素a浓度均值, 均值最低的时间段在1996~2005年, 平均浓度值为18.15 μg/L, 均值最高的时间段在2006~2015年, 平均浓度值为18.7 μg/L, 2016~2022年叶绿素a均值为18.48 μg/L, 高于1996~2005年。可以看出整个时间序列中叶绿素a浓度变化较小, 呈现先降低后高升, 再降低的趋势, 处于缓慢波动上升状态。

对1986~2022年所有湖泊叶绿素a浓度数据进行线性回归分析, 统计时间序列变化线性回归斜率, 当S >0时意味着该湖泊在近40年间的叶绿素a浓度为上升趋势, 当S < 0时则表示为下降趋势。根据回归分析结果, 绘制了湖泊叶绿素a变化趋势图(图 6), 由统计结果可以看出, 湖泊叶绿素a浓度变化率在−0.3~0.5, 上升率最大值为0.49, 下降率最大为−0.31, 平均值为0.07。其中占比最多的变化率区间为0~0.1, 占比为43.7%。其次为−0.1~0, 占比为36.4%。变化率分布在0.1~0.2的湖泊占总数量的8.9%。变化率分布在−0.2~−0.1的湖泊占总数量的7.7%。变化率大于0.2和小于−0.2的湖泊数量较少, 仅占湖泊总数的2.2%和2.1%。由各区间比例可以看出, 叶绿素a浓度处于上升趋势的湖泊多于下降趋势的湖泊。叶绿素a浓度轻度上升或下降的湖泊广泛分布于各区域。浓度上升幅度最大的区域集中于东北地区和新疆北部的部分地区。青藏高原地区的一些湖泊叶绿素a浓度明显降低, 也存在部分湖泊叶绿素浓度轻微上升。华东地区部分湖泊叶绿素a浓度为上升趋势, 部分湖泊为下降趋势。中南地区大部分湖泊处于上升趋势。

为了判断上述湖泊叶绿素a浓度变化趋势是否处于随机变动的范围, 即验证各湖泊的变化趋势是否显著, 使用t检验方法进行验证, 设定显著性水平P为0.05。根据图 6变化率显著性检验结果分析可知, 叶绿素a浓度显著上升和显著下降的湖泊在各区域均有分布, 其中约有45.2%的湖泊叶绿素a浓度变化不显著。而显著性上升湖泊数量占比大于显著下降的湖泊数量, 浓度显著上升的湖泊占比约为30%。显著下降的湖泊占比约为24.8%。表 4为显著性检验与湖泊面积区间统计, 可以看到三类变化趋势湖泊数量占比最多的均为面积在100 km²以下的区间, 分别占各自类别的95.7%, 84.1%和81.5%, 在该面积区间内, 显著上升的湖泊数量多于显著下降的湖泊数量。面积在100~250 km²的区间范围内, 显著下降的湖泊数量多于显著上升的湖泊数量。综合各区间数量占比来看, 除了变化趋势不显著的湖泊外, 面积在100 km²以下较小型湖泊叶绿素a浓度更多倾向于上升趋势。而面积在100 km²以上较大湖泊更多倾向于下降趋势。

通过计算湖泊叶绿素a浓度年序列变异系数, 可以比较不同区域湖泊相对自身的波动程度。如果变异系数较大, 则说明该湖泊叶绿素a浓度年序列波动程度较大, 即同一湖泊不同年份叶绿素a浓度差异较大, 反之则说明该湖泊叶绿素a浓度变化较为稳定。一般情况下, 将变异系数值划分五个等级来表征数据波动性, 如表 5所示。

图 7展示了在1986~2022年间叶绿素a浓度变异系数各区间分布图, 变异系数小于5%的湖泊占总数量的11.1%, 主要分布在华东地区, 中南地区和西南地区。占比最大的变异系数区间为5%~10%, 该区间的湖泊占总数量的43.9%, 该区间范围内的湖泊分布最广, 属于较低波动级别。其次为10%~15%, 该区间湖泊占比为28.7%, 属于中等级波动级别。变异系数大于15%的湖泊均集中分布于青藏高原地区, 东北地区和长江中下游地区, 其他区域分布数量较少。变异系数区间在15%~20%的湖泊占比为11.3%, 属于较高波动等级。变异系数区间在20%~25%和25%以上的湖泊数量较少, 占比约为3.7%和1.2%, 这些湖泊的叶绿素a浓度序列变化属于高波动范围, 大部分位于青藏高原地区。

3.4 各流域湖泊叶绿素a浓度时空变化特征

图 8是各流域内湖泊叶绿素a浓度梯度区间数量分组统计柱状图, 可以看到除了珠江流域和东南流域外, 其他流域湖泊叶绿素a浓度均在15~20 μg/L区间上数量最多, 而珠江流域和东南流域湖泊数量最多的区间分别为20~25 μg/L和25~30 μg/L。另外, 内陆流域和西南流域湖泊叶绿素a浓度区间分布较多的还包括10~15 μg/L, 珠江流域和东南流域湖泊数量较多的为25~30 μg/L区间, 其他流域占比区间较多的均为20~25 μg/L。由此分组柱状图进行比较可以看出内陆流域和西南流域湖泊叶绿素a浓度普遍较低, 珠江流域和东南流域湖泊叶绿素a浓度较高, 其他流域片区的浓度区间分布状态较为相似。

图 9为不同流域湖泊叶绿素a浓度年际均值时间变化序列, 分析可知, 湖泊叶绿素a浓度年序列均值较高的区域为东南流域和珠江流域, 大约为23 μg/L。年序列均值较低的区域为内陆流域和西南流域, 浓度范围大约在15~16 μg/L。各流域湖泊叶绿素a浓度变化较为稳定, 处于上升趋势的包括东南流域、松辽河流域、海河流域、淮河流域、珠江流域、长江流域和黄河流域。西南流域和内陆流域为下降趋势, 下降趋势斜率大约分别为−0.13和−0.004。变化趋势显著性检验P值的水平设为0.05, 即当P < 0.05时, 认为该变化趋势为显著。其中湖泊叶绿素a浓度年际均值变化趋势通过显著性检验的流域有松辽河流域, 珠江流域和长江流域, 这些流域湖泊叶绿素a浓度年际变化均呈现为波动上升的趋势, 变化最为明显的流域为珠江流域, 变化率大约为0.051, 其次是松辽河流域和长江流域, 变化率大约为0.037和0.017。可以看出, 无论是上升还是下降趋势, 除了个别年份有较大波动外, 这些流域年序列均值变化幅度较小, 总体较为稳定。

结合以上流域内部湖泊叶绿素a浓度年序列变化和流域整体的年际变化分析, 结果表明, 内陆流域和西南流域湖泊叶绿素a浓度普遍较低, 珠江流域和东南流域湖泊叶绿素a浓度较高, 其他流域片区的浓度区间分布状态较为相似。大部分流域内的湖泊叶绿素a浓度呈现上升趋势, 呈现下降的流域为西南流域和内陆流域, 在珠江流域和松辽河流域的湖泊叶绿素a浓度变化最为明显, 且均为显著上升趋势。

4 讨论

在本研究中, 基于Google Earth Engine遥感大数据平台, 采用空间分辨率更高的Landsat系列数据, 实现更大空间范围反演的同时大幅提升了水质参数的计算精度。研究内容从宏观层面分析了中国湖泊叶绿素a浓度的时空分异特征, 数据分析结果直观展现了中国范围内湖泊叶绿素a浓度的分布状态以及演变趋势。

本文研究结果表明, 中国湖泊叶绿素a浓度分布存在明显的季节性与地域性特征。空间上, 以胡焕庸线为界, 湖泊叶绿素a浓度在整体上呈现东南高, 西北低的态势, 这在很多大程度受到中国地势, 气候环境和人类活动差异的影响。东南地区相对湿润, 气温较高, 降雨充沛, 有利于有机质的分解和水体中营养物质的积累, 为藻类的生长提供了适宜环境。其次, 人口密集, 农业和工业发展较快, 增加了区域水体中的养分输入, 从而促进藻类的繁殖, 导致叶绿素a浓度相对较高(Yu et al, 2019)。相比之下, 西北地区地形较为干旱, 湖泊较少, 水体相对缺乏养分, 人类活动强度较弱, 使得该地区湖泊叶绿素a浓度较低。中国湖泊叶绿素a浓度存在季节性变化现象, 这可能由于夏秋季水温相对较高, 日照时间长, 光照充足, 水体较为稳定, 营养物质充沛, 为藻类的繁殖和生长提供了良好条件, 而春冬季则相反(孙力等, 2004)。可以看出湖泊叶绿素a浓度的季节性波动是温度、光照、水体稳定性和营养盐供应等多种因素综合影响的结果。

在过去的40年里, 中国湖泊叶绿素a浓度年均值呈现出缓慢波动上升的趋势。随着人口逐年增长, 工业和城市化的快速发展以及土地利用变化, 特别是农业和城市建设的扩张, 引入了大量污染物和养分进入湖泊水体, 导致了面源污染物的增加, 区域生态调节和自我恢复能力下降, 进一步加剧了湖泊的富营养化程度(陈小华等, 2018)。另外, 气候变化也可能通过影响水温, 水体循环和光照条件促进藻类的生长, 这一趋势对湖泊生态系统的稳定性产生了潜在影响(Kosten et al, 2012)。

中国各流域内湖泊叶绿素a浓度的时空分异特征呈现出明显的差异。内陆流域地理位置海拔较高, 气候干旱, 水体养分相对较少, 叶绿素a浓度相对较低。相反, 珠江流域和东南流域气候湿润, 湖泊富含养分, 为藻类繁殖提供了有利条件, 因此叶绿素a浓度也较高。在时间变化上, 除了西南流域和内陆流域的湖泊叶绿素a浓度呈现下降趋势外, 其他流域均为上升趋势。这可能与近年来人类活动的加剧有关, 随着流域人口增长, 工农业生产的快速发展以及流域内大规模经济开发, 需水量和污水排放量增大, 来自城市, 农田和工业区的污染物通过河流带入湖泊, 水体中养分富集, 使叶绿素a浓度逐渐上升, 引起水质恶化(赵杨等, 2020)。而西南流域和内陆流域的工业化, 城市化发展较缓, 大部分地区又受到高原地区特殊的气候环境影响, 如较低的气温和较短的生长季节等因素, 使得湖泊叶绿素a浓度偏低。

本文对中国湖泊进行了分区域和分流域两种形式的分析, 分区域讨论的重要性源于中国地域辽阔, 南北气温和东西海拔存在显著差异。通过分区域划分, 有助于控制多种环境因素引起的干扰。南北方向划分使得在同一纬度上具有相对一致的气温特征, 东西方向划分则考虑了相同经度上的海拔差异较小。因此, 通过分区域讨论湖泊叶绿素a浓度分布状态, 在一定程度可以减缓多种环境因素带来的交叉干扰, 从而更准确地分析湖泊叶绿素a浓度随经纬度变化的规律。分区域讨论侧重于地理差异, 有助于深入了解湖泊在不同地理环境下的特征, 更关注地理环境对湖泊的影响, 但在解释水文过程方面不如分流域方法直观, 割裂了湖泊与流域以及各流域之间自然联系。流域是由湖泊, 河流及其支流所组成的集水区域, 是由人与自然共同构成的复合系统, 包括社会、经济、资源和环境等诸多要素。流域内的水循环系统影响着湖泊水位, 水温和水质等, 将湖泊按照所在流域边界进行划分, 有助于分析湖泊的水文过程。以往将湖泊与流域分开, 相对封闭性地研究湖泊生态系统, 一定程度上阻碍了人们对湖泊与流域关系的认识, 如今以流域为单元进行综合研究和管理逐渐得到越来越多的学者重视(冷疏影等, 2003; 傅芳敏等, 2023)。

然而, 叶绿素a浓度算法模型不具有普适性, 中国陆地范围广阔, 不同区域湖泊的光学特性存在差异, 单一的叶绿素a浓度反演模型并不能完全适用于处在不同海拔, 不同气候区以及不同类型的所有湖泊。其次, 在春冬季节, 高原地区的湖泊表面结冰也会对反演结果产生一定影响, 所得数据不可避免存在着一定程度偏差。因此, 未来根据湖泊类型, 海拔高度和气候区范围等进一步细化分类, 采用更适合的模型进行反演, 将进一步提高结果的准确性。

5 结论

本文基于Landsat系列卫星等数据集, 提取了近40年的中国湖泊表面反射率影像, 以面积在10 km²以上的湖泊作为研究对象, 对湖泊叶绿素a浓度进行反演, 得到1986~2022年间中国湖泊叶绿素a浓度的时间序列数据, 以此从多个角度进行时空分异特征分析, 结论如下。

(1) 中国湖泊叶绿素a浓度存在地域性分布差异和季节性变化特征。湖泊叶绿素a浓度受到地形地貌与气候环境影响, 以胡焕庸线为界, 整体呈现东南高, 西北低的态势。所有湖泊中大约69%处于轻富营养化程度, 其次为中富营养化状态, 约占17%。胡焕庸线东侧的中富营养化状态湖泊明显多于西侧。绝大部分地区夏季湖泊叶绿素a浓度最高, 冬季最低, 各季节空间分布均为东高西低。

(2) 湖泊叶绿素a浓度随经纬度的变化规律为: 35°N以北, 叶绿素a浓度由西向东总体呈现上升趋势。35°N以南, 由西向东, 叶绿素a浓度呈现先下降, 后大幅上升, 再下降的变化趋势。100°E以西, 随纬度升高, 叶绿素a浓度变化为先下降后上升。100°E以东的湖泊叶绿素a浓度普遍偏高, 由北向南呈现波动上升趋势。

(3) 近40年间中国湖泊叶绿素a浓度年均值处于缓慢波动上升趋势。浓度均值最低的时间段在1996~2005年, 均值最高的时间段在2006~2015年。叶绿素a浓度显著上升和显著下降的湖泊在各区域均有分布。100 km²以下的湖泊大部分处于显著上升趋势, 100 km²以上的较大型湖泊则倾向于下降趋势。整体来看湖泊叶绿素a浓度变化较为稳定, 波动较大的区域位于青藏高原, 东北地区和长江中下游的部分湖泊。

(4) 各流域内湖泊叶绿素a浓度时空分异特征表现为: 内陆流域和西南流域的湖泊叶绿素a浓度普遍较低, 珠江流域和东南流域的湖泊叶绿素a浓度较高, 其他流域片区分布状态较为相似。除了西南流域和内陆流域内的湖泊叶绿素a浓度呈现下降趋势外, 其他流域均为上升趋势, 珠江流域和松辽河流域变化最为明显, 且均为显著上升趋势。