基于MODIS数据的全国144个重点湖库营养状态监测：以2018年夏季为例*

曹 畅,王胜蕾,李俊生,赵红莉,申 维,谢 娅

(1:中国地质大学(北京)地球科学与资源学院，北京 100083)(2:中国科学院空天信息创新研究院，数字地球重点实验室，北京 100094)(3:北京大学遥感与地理信息研究所，北京 100871)(4:中国科学院大学电子电气与通信工程学院,北京 100049)(5:中国水利水电科学研究院, 水资源研究所,北京 100048)

水体富营养化是水体中营养盐增多而导致浮游植物爆发性生长的一种现象，富营养化的湖库水体透明度和溶解氧浓度下降，容易暴发蓝藻水华，对水环境质量造成严重的威胁[1-3]. 因此，动态监测和评价湖库水体的营养状态具有非常重要的意义. 水体营养状态的常规监测方法是水面采集水样并送到实验室内测量叶绿素a、总氮、总磷等水质参数，进而计算营养状态指数，然后根据营养状态指数分级可以将水体营养状态分为贫营养、中营养和富营养[4-6]. 这种方法虽然精度较高但费时、费力、费用高，而且只能获得个别采样点的结果. 与之相比，基于卫星遥感监测水体营养状态具有快速、范围广、成本低等优势，适合开展大范围动态监测.

基于卫星遥感监测水体营养状态主要是基于卫星遥感反演的叶绿素a浓度计算水体营养指数，进而得到水体的营养状态[7-8]. Thiemann等利用IRE-1C数据反演德国梅克伦堡州湖泊水体中的叶绿素a浓度，基于卡尔森模型评价其水体的富营养化状态[9]；Matthews等通过测量内陆水体叶绿素a浓度，监测非洲南部多个湖库水体10年间的富营养化及其变化情况[10]. Gilerson等利用红色和近红外波段算法评价美国布拉斯加州的叶绿素a浓度，观测该州湖泊水体的营养状态[11]；朱利等利用GF-1号卫星WFV数据对太湖叶绿素a浓度、透明度、悬浮物浓度进行遥感监测，并评价水体富营养状态和水质参数空间分布[12]；Le等在三波段算法的基础上发展四波段解析算法提高三波段算法的性能，比较两种算法在反演太湖叶绿素a浓度的精度，并对太湖营养状态进行评价[13].

然而，由于内陆水体光学特性复杂，受浮游植物、悬浮泥沙、有色可溶性有机物共同影响，而且随区域和季节变化大，因为不同季节的水体温度变化对于藻类生长的影响、地表径流和风速变化对于悬浮物的影响往往存在差异，进而导致很多湖库水体叶绿素a反演算法都具有很强的区域和季节局限性[14-16]，因此限制了基于卫星遥感的大范围水体营养状态监测. 与传统的基于叶绿素a反演结果的营养状态卫星遥感监测不同，Wang等[17]发展了基于水色指数的水体营养状态等级评价方法，并且成功应用于2012年全球大型湖库水体营养状态评价，利用文献和公报中的100个水体实地监测结果进行评价得到精度为80%，这为开展大范围湖库水体营养状态评价奠定了基础.

随着工农业的快速发展，我国湖库水体富营养化的趋势较为严重，开展全国范围内重点湖库营养状态的动态监测具有重要意义，有利于宏观掌握全国湖库营养状态的时空分布，为管理决策提供数据支持. 但是，《2018年中国生态环境状态公报》只监测了111个湖库的营养状态，前几年监测的湖库数量更少. 因此，非常有必要开展基于卫星遥感的全国范围内湖库的营养状态监测，作为公报的补充，一方面可以从空间上监测更多的湖库，另一方面可以从时间上追溯湖库营养状态的变化. 因此本文将利用基于水色指数的营养状态评价方法[17]，生产全国范围内重点湖库的营养状态遥感监测产品，进而分析全国重点湖库营养状态空间分布规律及其主要影响因素.

1 数据源和研究区

1.1 数据源

搭载在Terra和Aqua卫星的中分辨率成像光谱仪(MODIS)具有较高的时间分辨率，一天内可以覆盖全球两次(部分赤道低纬度地区除外). MOD09数据是MODIS(Terra)的地表反射率产品，包含MODIS的1～7共7个波段. MOD09A1数据是八天合成的MOD09数据，降低了部分有云等质量不好数据的影响而且将全球范围划分为460个规则的网格分块，便于生产大范围产品并进行时序产品叠加分析，因而在业务化的大范围水质监测中具有重要优势. 因此本文选择MOD09A1数据作为全国重点湖库营养状态业务化监测的遥感数据源，利用覆盖全国重点湖库的18块MOD09A1产品，网格分块行列号包括：h23v04、h23v05、h24v04、h24v05、h25v03、h25v04、h25v05、h25v06、h26v04、h26v05、h26v06、h27v04、h27v05、h27v06、h28v05、h28v06、h28v07、h29v06.

1.2 研究区

MOD09A1数据的空间分辨率是500 m，为了降低岸边混合像元或临近像元的影响，本文主要研究面积大于5 km2的湖库水体，也就是水体的联通像元数量大于20个.

在全国面积大于5 km2的湖库中选择了144个主要湖库，包括111个湖泊和33个水库. 111个湖泊中包括洱海等5个饮用水源地，也包括巢湖等较浑浊和抚仙湖等较清洁的106个非饮用水源地湖泊；33个水库包括密云水库等19个饮用水源地，以及于桥水库等较浑浊和新安江水库等较清洁的14个非饮用水源地水库. 这些湖库的空间分布覆盖了中国的五大湖区，其中20个位于蒙新高原湖区，34个位于青藏高原湖区，7个位于云贵高原湖区，66个位于东部平原湖区，17个位于东北山地与平原湖区.

2 研究方法

本文首先对MOD09A1地表反射率产品进行离水反射率校正，然后提取湖库水体分布范围，进而计算水色指数(FUI)，并利用FUI生产水体营养状态等级，最后分析全国重点湖库营养状态的空间分布规律.

2.1 离水反射率校正

MOD09A1是地表反射率产品，已对其7个波段进行了气溶胶校正、卷积云校正和瑞利散射校正[18]，但是没有对天空光及太阳耀斑的水面反射进行校正. 本文采用基于近红外波段(NIR)和短波红外波段(SWIR)的暗目标法进行校正[19]，将MOD09A1的地表反射率(R(λ))减去这个像元的近红外和短波红外的最小值，再除以π，转换成离水反射率(Rrs(λ)).

2.2 湖库水体提取

在对MOD09A1数据进行离水反射率校正的基础上，还要对湖库的水体边界进行提取. 水体在短波红外波段具有强吸收作用，因而水体像元在短波红外波段反射率很低，相反陆地像元的反射率较高[20]，因此本文采用基于短波红外波段灰度直方图的双峰谷值阈值自动化确定的方法提取水体边界[21]. 为了降低接近陆地区域混合像元和光学浅水的影响，将得到的水体部分向内腐蚀一个像元. 最后，去除联通像元数小于20的水体，只保留面积在5 km2以上的水体.

2.3 水色指数计算

利用Wang等[22]提出的从太湖MOD09数据产品中提取FUI指数的方法对全国主要湖库的FUI指数进行求解.FUI是一种比色表，将水体颜色从深蓝色到黄褐色划分为21个级别. 利用MODIS的第1(红)、3(蓝)、4(绿)波段计算CIE颜色空间中的三刺激值[x]、[y]、[z][23]，并进一步计算CIE颜色空间色度图中找到颜色对应的坐标(x,y)，进而找到色度角α，在FUI指数查找表中找到与色度角最接近的色度值，此值对应的FUI值就是水体的FUI值[22].

2.4 基于FUI的营养状态等级评价

采用基于FUI和红波段阈值法对湖库营养状态进行评价[22]. 贫营养水体一般比较清洁，藻类浓度比较低，水体一般呈现蓝色；富营养水体较为浑浊，藻类浓度较高，水体颜色多为黄绿色；中营养水体的藻类浓度以及水体颜色一般介于贫营养和富营养之间.FUI可以表示水体颜色，因此可以使用FUI的阈值分割来评价水体的营养状态；当FUI<7时，水体为贫营养；当7≤FUI<10时，水体为中营养；当FUI≥10时，水体为富营养. 但是，也有一些中营养水体的FUI≥10，使用上述分割模型会被误判为富营养水体. 通过分析发现这些误判水体的红波段离水反射率(Rrs(645))一般较小，因此增加一个阈值分割限制条件进行二次限定：当Rrs(645)<0.00625时，即使FUI≥10，水体仍然是中营养[22].

2.5 精度评价

可以利用与MODIS卫星同步测量的水面实测数据计算水体营养状态等级，对MODIS监测结果进行精度评价. 于2018到2019年在太湖和于桥水库分别获取了12个和46个采样点的叶绿素a浓度实测数据，利用这些叶绿素a浓度数据可以计算综合营养状态指数(TLI)[2]. 由于这些点大部分属于富营养，为了增加贫营养水体的检验数据，从文献[24]中获取了2014-2016年青藏高原纳木错、塔若错、色林错、多尔索洞错、赤布张错5个湖泊的透明度数据，利用营养状态指数计算公式：TLI(Chl.a)=10(2.5+1.086 ln Chl.a)、TLI(TP)=10(9.436+1.624 ln TP)、TLI(TN)=10(5.453+1.694 ln TN)、TLI(SD)=10(5.118-1.94 ln SD)、TLI(COD)=10(0.109+ 2.661 ln COD)(式中, 叶绿素a(Chl.a)单位为mg/m3；透明度(SD)单位为m; 其他指标单位均为mg/L)可以计算综合营养状态指数，通过TLI数值大小和分级标准对湖库水体进行营养状态分级[2]. 利用这63个采样点实测数据计算得到的营养状态等级与同步MODIS卫星遥感监测的营养状态等级进行对比，计算得到基于MODIS监测营养状态等级的总体精度为88.9%. 其中，贫营养的判别精度比较高，少量误判主要发生在部分富营养和中营养之间. 总体上来说，营养状态评价的精度基本上可以满足大范围应用研究需求. Guan等[25]在2018年实地测量的叶绿素a浓度数据中有26个湖泊与本文研究区域一致，采用综合营养状态指数法计算26个湖泊的营养状态，有25个湖泊与本文研究结果相同，MODIS的监测精度达到96.0%，邹伟等[26]在2018年7-8月份采集了长江中下游地区19个湖库水体的叶绿素a浓度，其中与本文研究区重叠的湖库共有13个，基于叶绿素a浓度计算的营养状态等级为12个湖库为富营养化，1个湖库为中营养. 这13个湖库中的12个的营养状态调查结果与本文基于MODIS监测的结果一致，精度为92.3%. Song等[27]对2014-2018年夏季全国面积大于80000 m2的湖库进行大范围高精度水体透明度反演，利用综合营养状态指数计算得到东北山地与平原湖区和东部平原湖区水体较多呈现富营养状态，对比图1，与本文研究的结果一致.

2.6 全国重点湖库营养状态产品生产及空间分析

基于上述技术流程，编写IDL程序，实现基于MOD09A1数据的全国重点湖库营养状态产品自动化批量生产. 由于水色指数的营养状态评价方法主要是基于藻类对水体颜色的影响，因此选择藻类生长情况最好的夏季来分析全国重点湖库水体营养状态的空间分布情况，进而分析其主要影响因素.

3 结果与讨论

3.1 全国重点湖库营养状态空间分布格局

基于2018年6-8月之间的MOD09A1数据，生产了全国144个重点湖库的营养状态夏季均值产品. 为了便于统计分析，将每个湖库所有像元的营养状态求一个均值，得到每个湖库平均的营养状态等级，其中贫营养、中营养、富营养的湖库数量分别为23、34、87个，144个重点湖库营养状态空间分布如图1所示，可以看出我国湖库水体的营养状态分布不均匀，东部湖库水体以中到富营养状态为主，尤其是长江中下游和东北山地与平原湖区湖库水体富营养化比例非常高，西部湖库水体以贫到中营养状态为主，尤其是青藏高原湖区贫营养比例很高. 不同水体类型营养状态等级比例如表1，其中湖泊型饮用水源地的富营养化比例要高于水库型饮用水源地.

表1 基于MODIS监测的2018年夏季全国重点湖库中不同类型水体的营养状态等级比例

五大湖区中各营养状态等级所占的比例如图2所示，东北山地与平原湖区水体富营养状态比例最大，东部平原湖区次之，2个湖区均无贫营养状态水体. 中营养状态水体在蒙新高原湖区数量最多. 青藏高原湖区贫营养状态水体数量最多，其次是云贵高原湖区. 在五大湖区中各选1个典型湖泊，营养状态分布如图3所示. 其中，兴凯湖(东北山地与平原湖区，中俄跨界湖)、巢湖(东部平原湖区)水体富营养化程度较高；博斯腾湖(蒙新高原湖区)整体上为中营养状态；色林错(青藏高原湖区)、抚仙湖(云贵高原湖区)整体上为贫营养，水质情况较好.

图1 基于MODIS监测的2018年夏季全国重点湖库水体营养状态等级分布

自然环境是影响水体营养状态的重要因素. 青藏高原湖区的海拔较高, 多在3500 m以上，气温较低[28-29]，水中藻类进行光合作用的速率受低水温限制[30]，藻类数量一般较少，水体多为贫营养状态[31]. 云贵高原湖区深水湖和浅水湖交错分布，该湖区水体营养状态差异较大[32]. 东北地区的中部为平原, 周围三面环山，位于平原地区的湖泊面积较小，且湖水较浅[33]，降水多集中在夏季，会带入营养盐等物质，造成富营养状态水体数量较多. 东部平原湖区人口密度高，面源污染带来大量的营养盐；同时长江中下游地区主要湖泊的水深多在10 m以下[34-35]，湖泊底泥在环境条件适宜时与湖水进行交互，也成为营养盐的主要来源之一[36]. 蒙新高原湖区降水少，蒸发强，该湖区草场退化且土地沙漠化，造成地表径流减少，再加上水土流失导致入湖水体中营养盐浓度上升，水体营养状态较高[37-38].

人类活动也是影响水体营养状态的重要因素，湖泊的富营化水平与社会经济发展进程密切相关[39]. 东部平原湖区和东北山地湖区人口众多，加上经济发展速度快，工业点源污染以及农业非点源污染都会为周边湖库带来营养物质[40-41]，导致湖库水体的营养状态整体偏高. 青藏高原湖区海拔较高，自然环境较恶劣，人口密度小，城市化进程缓慢，水体受人为干扰因素较小，水体多为贫营养状态.

3.2 营养状态等级监测结果与环境公报的对比

基于FUI的营养状态评价方法是利用影响水体颜色的叶绿素a和透明度信息，因此应该使用藻类生长最充分的夏季，才能对比出来不同水体的差异. 而且很多北方和青藏高原湖库在冬季结冰，也无法使用. 综上所述，基于FUI的营养状态评价方法主要适用于夏季. 《2018中国生态环境状态公报》监测了111个湖库2018年平均的营养状态，其中有40个湖库与本文研究区重叠. 经过对比发现，其中，有24个湖库(占60%)的营养状态与本文监测结果一致，另外16个湖库(占40%)的营养状态不一致，主要是公报评价为中营养，而本文评价为富营养. 公报方法与本文方法不同，主要体现在：(1)公报是基于叶绿素a、透明度、总氮、总磷、化学需氧量5个参数计算综合营养状态指数，而本文是基于水体颜色信息评价水体营养状态，其中主要是利用了对水体颜色影响较大的叶绿素a和透明度信息；(2)公报一般是基于每月一次测量的少数几个采样点的数据的评价结果，而本文是基于时间序列的全湖遥感图像平均的结果. 本文的营养状态评价方法可以作为公报的补充，一方面可以用于监测公报没有监测的湖库，另一方面对于公报监测的湖库，还可以补充其空间分布和时间变化信息.

图2 基于MODIS监测的2018年夏季全国五大湖区湖库营养状态等级比例

图3 基于MODIS监测的2018年夏季全国五大湖区典型湖泊水体营养状态

4 结论

本文以MODIS八天合成的地表反射率产品MOD09A1为遥感数据源，利用基于水色指数的营养状态分级方法，监测全国144个重点湖库营养状态等级. 利用63个水面实测数据计算得到的营养状态等级与同步MODIS卫星遥感监测结果进行对比，得到MODIS的监测精度为88.9%，精度基本上可以满足大范围应用研究需求.

利用经过检验的方法，基于MODIS数据生产了2018年夏季全国144个重点湖库的营养状态等级，分析了全国重点湖库营养状态的空间分布情况(图2). 贫营养、中营养、富营养的湖库比例分别为16%、24%、60%，富营养化的比例很高. 营养状态空间上分布不均匀，东部地区湖库(主要是东部平原和东北山地与平原湖区)多为富营养状态，东北山地与平原湖区的富营养化湖库比例最多达到94%，东部平原湖区的富营养化湖库比例达到80%. 西部湖区湖库(特别是青藏高原湖区)多为贫营养状态，青藏高原湖区贫营养湖库比例达到62%，富营养化水体仅占9%. 海拔和地表温度等自然因素与工业点源和农业面源污染等人为因素是湖库营养状态的重要影响因素.