MSI的湖泊富营养指数遥感估算方法

本发明涉及一种msi的湖泊富营养指数遥感估算方法。

背景技术：

全球面积大于25km²的2058个湖泊中，63.1％处于富营养化状态，亚洲54％、欧洲53％、北美48％、非洲28％和南美41％的湖泊受到富营养化影响。由于陆源输入营养负荷增加和生产力急剧增长，缺氧和藻华事件频繁暴发，湖泊富营养化已经成为影响大部分淡水生态系统的全球水质问题。

湖泊水生态环境监测的常用范式是富营养化状况，它是光、热、水动力学和营养物如氮和磷的结合形成的一系列生物、化学和物理过程。

目前，国内外普遍采用“野外采样—室内试验”方法，以实测总磷、叶绿素a和透明度等浓度作为富营养化参数指标，引入一种被广泛采用的数值—卡尔森富营养化状态指数(trophicstatusindex，tsi)，定量富营养状态，表征湖泊为“寡养”、“中养”或“富营养化”等级。

传统的野外离散采样和室内分析测试方法，一方面需要大量人力、物力和财力且耗时长，水样的保存时间、保存温度、不同仪器测量精度、水样前处理操作等都会对tsi指数计算的结果造成较大影响。其次，湖泊单个样本的营养程度不能代表整个湖泊，且部分湖泊由于地处偏远，难以进行现场测量，限制区域尺度或全球尺度湖泊富营养化评估，缺乏对大区域尺度上湖泊生物地球化学循环过程的理解和富营养化治理政策综合实施。

卫星遥感技术被广泛用于监测湖泊水生环境，并基于反射率光谱和存档数据提供时空连续的动态变化，可弥补传统方法的不足。

富营养指数tsi三个指标中，叶绿素a作为光学活性物质具有独特的光学特性，富营养化遥感监测主要通过估算水体中叶绿素a的浓度来实现。叶绿素a定量遥感反演的数据源主要有oli、modis、spot、msi、olci和meris，反演模型主要有基线法、波段比值模型、三波段模型、四波段模型、qaa生物光学模型和人工智能机器学习算法(gitelsonetal.,2008；goweretal.,2005；gurlinetal.,2011；mosesetal.,2009；shenetal.,2010；yangetal.,2010)。近年来，随着人工智能新型算法发展，线性回归算法、神经网络法、随机森林和boost算法在不同程度上提高叶绿素a浓度估算精度。lisijia等(2021)基于msi数据产品验证了两波段模型、三波段模型、三波段叶绿素指数模型、mci叶绿素指数模型、sci综合叶绿素指数模型以及人工智能机器算法如线性回归、支持向量机和catboost算法在中国典型湖泊的适用性，其中支持向量机获得较高的反演精度。watanabeetal.,(2015)应用landsatoli数据反演叶绿素a浓度，进而估算水体富营养指数tsi。

富营养指数tsi三个指标中，总磷和透明度参数无显著的光学特性。湖泊透明度遥感监测是通过不同遥感数据源获取的反射率与透明度实测值建立经验算法(olmansonetal.,2008；songetal.,2020)。虽然部分研究已尝试构建水体透明度的半分析算法，但其更适用于大洋或沿岸水体(leeetal.,2016)，针对光学特性复杂的内陆湖泊水体精度较低。olmansonetal.,(2008)采用landsat系列产品反演明尼苏达州湖泊透明度，进而计算tsi和评价富营养状况。总磷遥感反演监测是通过总磷实测值与水体光学活性物质(如叶绿素a和总悬浮物浓度)的相关关系而间接反演(kutseretal.,1995；sunetal.,2013)。

但富营养化指数tsi是反映湖泊水质状况的综合指标，叶绿素a、透明度和总磷等单个指标的遥感识别、监测以及定量估算不能反映湖泊综合富营养状态。songetal.,(2012)应用高光谱数据分别用估测叶绿素a、透明度和总磷的浓度值，进一步计算印第安纳水库富营养化指数tsi。sheelaetal.,(2011)应用irslissiii数据产品反演水体透明度和叶绿素a估算tsi。shikunetal.,(2019)运用landsatoli数据产品和半分析算法估算光学活性物质的总吸收系数，进而估算tsi。

综上，所述湖泊富营养状况的遥感监测多采用一个富营养化指标、多个富营养化指标或者以光学活性物质吸收系数为间接参数的遥感识别、监测以及定量估算，进而估算tsi。这种湖泊富营养指数tsi估算精度受限于多个指标独立反演不确定性的叠加，且由于内陆水体光学的复杂性和经验算法时空移植性差等缺陷，tsi估算难以在湖泊富营养大尺度应用中开展。

技术实现要素：

本发明的目的是为了提高湖泊富营养状态评价的系统性和精确程度，提供了一种msi的湖泊富营养指数遥感估算方法。

msi的湖泊富营养指数遥感估算方法按照以下步骤进行：

步骤一、通过连续三年实地采集湖泊水体样品，测定叶绿素a、总磷和透明度水质参数的浓度，根据修订的卡尔森富营养化指数公式，计算湖泊水体样品富营养化指数tsi值，并进行湖泊富营养状况评价，将湖泊富营养化划分三个等级:tsi<30,贫营养状态等级；30<tsi<50,中营养等级；50<tsi<100,富营养等级；

步骤二、结合连续三年的内陆湖泊水体星地同步数据，将msilevel1c遥感影像经过重采样和大气校正预处理，提取湖泊水体采样点对应的遥感反射率，得到的第一波段至第六波段的遥感反射率值依次记为rhown1至rhown6；

步骤三、计算rhown1和rhown4、rhown4和rhown5的比值，将三分之二的数据用于构建湖泊富营养指数tsi反演模型；

步骤四、步骤三中剩余三分之一的数据用于湖泊富营养指数tsi反演模型的精度验证；

步骤五、将步骤三所构建的湖泊富营养指数tsi反演模型作用于中国典型内陆湖泊的msi遥感影像上，得到典型湖泊对应富营养指数的时空分布。

步骤一中计算富营养化指数tsi采用的是修订的卡尔森富营养化指数公式，其计算公式表达为：

tsi＝0.54×tsim(chl-a)+0.297×tsim(sdd)+0.163×tsim(tp)(4)

式中，chl-a、tp和sdd分别为叶绿素a、总磷和透明度，叶绿素a的单位：μgl^-1，总磷的单位：mgl^-1，透明度的单位：m；tsim(chl-a)、tsim(sdd)和tsim(tp)分别表示叶绿素a、总磷和透明度相对应的tsi值。

步骤二中所述重采样预处理是基于eas的snap7.0.0软件平台，将msilevel1c影像数据重采样至10m。

步骤二中大气校正采用的是eas开发的cr2cc-nets算法，可在snap7.0.0软件平台实现，在该算法中，所需盐度、温度和海拔高度参数为实地采样湖泊水体测定的盐度值、水温值和海拔高度值；臭氧层厚度为默认值330du；气压为默认值1000pa；总悬浮物bpart因子默认值为1.72，bwit因子默认值为3.1；叶绿素指数chlexp因子默认值为1.04，chlfak因子默认值为21.0；tosa神经网络训练数据集阈值默认值为0.05；大气校正反射率神经网络训练数据集阈值默认值为0.1；云量检测下行透射比默认值为0.955。

步骤三中所述富营养指数tsi反演模型如下：

y＝-34.04×x1-1.114×x2+97.376，r²＝0.80(5)

其中y为富营养指数tsi，x1为rhown4和rhown5的比值，x2为rhown1和rhown4的比值。

步骤四中所述精度验证的表达式如下：

其中，rmse是均方根误差，mape是平均绝对百分比误差，n是样本数，yi是测量值，yi'是估算值。

依据本专利方法计算出来的tsi具有较高的可信度。湖泊富营养化指数tsi和rhown4/rhown5(rrs(665)/rrs(704))、rhown1/rhown4(rrs(443)/rrs(665))的相关性模型是本发明首先确立的。

本发明方法通过构建湖泊富营养化tsi指数反演模型，确定了一种新的、便捷的估算湖泊富营养化指数方法，即只需提取湖泊哨兵2号msilevel1c产品rhown4/rhown5(rrs(665)/rrs(704))、rhown1/rhown4(rrs(443)/rrs(665))反射率比，就可通过该模型计算出湖泊的富营养化指数tsi值，进而针对湖泊富营养状况评价。

本发明所述的湖泊富营养指数tsi反演模型，具有较高的反演精度，均方根误差可达到5.47，平均相对误差为9.48％，相对于现有技术精度显著提高；另一方面，msi传感器是搭载在哨兵sentinel卫星，自2015年发射，较高的时间分辨率(2-5天)和空间分辨率(10,20和60m)使其在内陆湖泊水生态系统监测具有广阔应用前景；本发明所述方法的湖泊富营养指数tsi反演模型，是基于中国内陆湖泊大尺度实时监测数据建立，在内陆湖泊水生态系统健康监测和评估的通用性方面具有更好的潜力，极大地提高湖泊富营养化评估的效率和规模。

本发明将有助于评估湖泊水体富营养状态的时空变化及其发展趋势，可有效评估湖泊富营养状态和生态修复的绩效，对分析湖泊生态系统的结构和功能、预测湖泊富营养化状态和制定防控政策具有重要意义，为政府部门的水环境保护、水生态资源管理的科学决策提供科技支撑。本发明的湖泊富营养化指数遥感估算方法，提高了湖泊富营养化监测的效率和精度。

附图说明

图1是实验一中中国典型湖泊采样分布图(2017-2019年)；

图2是实验一中哨兵2号level-1c产品在snap软件处理msilevel-1c产品流程图；

图3是实验一中rhown4/rhown5(rrs(665)/rrs(704))、rhown1/rhown4(rrs(443)/rrs(665))和tsi的建模拟合图；

图4是实验一中rhown4/rhown5(rrs(665)/rrs(704))、rhown1/rhown4(rrs(443)/rrs(665))和tsi的精度验证拟合图；

图5是实验一中计算得到查干湖4月的tsi的时空分布图；

图6是实验一中计算得到查干湖8月的tsi的时空分布图；

图7是实验一中计算得到查干湖10月的tsi的时空分布图。

具体实施方式

本发明技术方案不局限于以下所列举具体实施方式，还包括各具体实施方式间的任意组合。

具体实施方式一：本实施方式中msi的湖泊富营养指数遥感估算方法按照以下步骤进行：

步骤一、通过连续三年实地采集湖泊水体样品，测定叶绿素a、总磷和透明度水质参数的浓度，根据修订的卡尔森富营养化指数公式，计算湖泊水体样品富营养化指数tsi值，并进行湖泊富营养状况评价，将湖泊富营养化划分三个等级:tsi<30,贫营养状态等级；30<tsi<50,中营养等级；50<tsi<100,富营养等级；

步骤二、结合连续三年的内陆湖泊水体星地同步数据，将msilevel1c遥感影像经过重采样和大气校正预处理，提取湖泊水体采样点对应的遥感反射率，得到的第一波段至第六波段的遥感反射率值依次记为rhown1至rhown6；

步骤三、计算rhown1和rhown4、rhown4和rhown5的比值，将三分之二的数据用于构建湖泊富营养指数tsi反演模型；

步骤四、步骤三中剩余三分之一的数据用于湖泊富营养指数tsi反演模型的精度验证；

步骤五、将步骤三所构建的湖泊富营养指数tsi反演模型作用于中国典型内陆湖泊的msi遥感影像上，得到典型湖泊对应富营养指数的时空分布。

具体实施方式二：本实施方式与具体实施方式一不同的是步骤一中计算富营养化指数tsi采用的是修订的卡尔森富营养化指数公式，其计算公式表达为：

tsi＝0.54×tsim(chl-a)+0.297×tsim(sdd)+0.163×tsim(tp)(4)

式中，chl-a、tp和sdd分别为叶绿素a、总磷和透明度，叶绿素a的单位：μgl^-1，总磷的单位：mgl^-1，透明度的单位：m；tsim(chl-a)、tsim(sdd)和tsim(tp)分别表示叶绿素a、总磷和透明度相对应的tsi值。其他与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是步骤二中所述重采样预处理是基于eas的snap7.0.0软件平台，将msilevel1c影像数据重采样至10m。其他与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是步骤二中大气校正采用的是eas开发的cr2cc-nets算法，可在snap7.0.0软件平台实现，在该算法中，所需盐度、温度和海拔高度参数为实地采样湖泊水体测定的盐度值、水温值和海拔高度值；臭氧层厚度为默认值330du；气压为默认值1000pa；总悬浮物bpart因子默认值为1.72，bwit因子默认值为3.1；叶绿素指数chlexp因子默认值为1.04，chlfak因子默认值为21.0；tosa神经网络训练数据集阈值默认值为0.05；大气校正反射率神经网络训练数据集阈值默认值为0.1；云量检测下行透射比默认值为0.955。其他与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是步骤三中所述富营养指数tsi反演模型如下：

y＝-34.04×x1-1.114×x2+97.376，r²＝0.80(5)

其中y为富营养指数tsi，x1为rhown4和rhown5的比值，x2为rhown1和rhown4的比值。其他与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是步骤四中所述精度验证的表达式如下：

其中，rmse是均方根误差，mape是平均绝对百分比误差，n是样本数，yi是测量值，yi'是估算值。其他与具体实施方式一至五之一相同。

采用下述实验验证本发明效果：

实验一：

msi的湖泊富营养指数遥感估算方法按照以下步骤进行：

步骤一、对分布在全国的42个湖泊和水库(以下统称湖泊)进行采样(图1)，每个湖泊设置3-8个采样点，共计获取431个采样点。使用塞氏盘现场测定水体透明度，精度为0.01m，每个样点重复三次，记录并计算平均值作为该点的透明度实测值。同时使用采水器采集表层水体(0-0.5m)，每个样点采集2l，置于车载冰箱中运送至实验室避光冷藏保存，运送到实验室24h内测定其总磷和叶绿素a浓度。

步骤二、步骤一所述叶绿素a和总磷预处理和测定操作如下：将水样经0.45μm孔径的玻璃纤维滤膜(中国半岛工业有限公司生产)，将滤膜剪碎后放入玻璃离心管，加入10ml质量浓度为90％的丙酮溶液萃取，混合均匀后避光冷藏静置24h，对其研磨、离心后取上清液。以质量浓度为90％丙酮溶液为空白参考，将上清液放置到紫外可见光分光光度计(型号uv-2600pc，日本岛津公司生产)测定波长为630nm、647nm、664nm和750nm的吸光度，并参考叶绿素a-分光光度法(hj897-2017)公式计算叶绿素a浓度，如下：

chl-a—水样中叶绿素a的浓度，单位为μgl^-1；

d630—待测样品在630nm波长下的吸光度值；

d647—待测样品在647nm波长下的吸光度值；

d664—待测样品在664nm波长下的吸光度值；

d750—待测样品在750nm波长下的吸光度值；

v1—水体样品的定容体积，单位为ml；

v—取样体积，单位为l；

参照地表水环境质量标准(gb3838-2002)，总磷测试方法采用钼酸铵分光光度法(gb11893-89)。每个水样重复测定叶绿素a和总磷三次，记录并计算平均值作为该点的实测值。

步骤三、根据修订的卡尔森富营养化指数公式，富营养化指数tsi由tsim(chl-a)、tsim(sdd)和tsim(tp)计算，如下：

式中，chl-a和tp为叶绿素a(μgl^-1)和总磷(mgl^-1)浓度值，sdd为透明度深度(m)。富营养化指数tsi可计算如下：

tsi＝0.54×tsim(chl-a)+0.297×tsim(sdd)+0.163×tsim(tp)(12)

湖泊富营养化可划分三个等级:tsi<30,贫营养状态等级；30<tsi<50,中营养等级；50<tsi<100,富营养等级。

步骤四、利用步骤一获取的实地采样湖泊星地同步数据，下载经过正射校正和几何精校正的msilevel-1c级多光谱数据，累计获取253景无云level-1c级多光谱数据覆盖对应湖泊，考虑影像时间窗为±7天。

步骤五、步骤四所述level-1c级多光谱数据，在eas的snap7.0.0软件平台重采样至10m。

步骤六、步骤五所述获取重采样后的数据在snap7.0.0软件平台实现cr2cc-nets大气校正，获得大气校正后的rhown1[rrs(443)]、rhown2[rrs(492)]、rhown3[rrs(560)]、rhown4[rrs(665)]、rhown5[rrs(704)]和rhown6[rrs(740)]遥感反射率。

步骤七、在cr2cc-nets大气校正算法中，所需盐度、温度和海拔高度参数为步骤一实地采样湖泊水体测定的盐度值、水温值和海拔高度值；臭氧层厚度为默认值330du；气压为默认值1000pa；总悬浮物bpart因子默认值为1.72，bwit因子默认值为3.1；叶绿素指数chlexp因子默认值为1.04，chlfak因子默认值为21.0；tosa神经网络训练数据集阈值默认值为0.05；大气校正反射率神经网络训练数据集阈值默认值为0.1；云量检测下行透射比默认值为0.955。

步骤八、基于步骤一所述实地湖泊采样gps信息，在步骤六获取的校正后数据提取3×3像元rhown1-6反射率值，步骤五、s6和s7步骤都是在snap7.0.0软件平台中执行(图2)。

步骤九、对所述步骤三计算采样点的富营养指数tsi和步骤八提取rhown1-6反射率值星地同步数据，在matlab软件中随机分组，其中287个星地同步样本数据用于构建tsi反演模型建模数据组，144个星地同步样本数据用于tsi反演模型验证数据组。

步骤十、基于步骤九获取的建模数据组，采用microsoftexcel2017软件进行富营养化指数tsi和rhown1-6反射率相关性拟合分析，该相关性可以用以下模型进行描述：tsi＝-34.04×[rrs(665)/rrs(704)]-1.114×[rrs(443)/rrs(665)]+97.376(r²＝0.80，n＝287)(图3)，数据点在回归线两侧均匀分布。其中，均方根误差(rmse)和平均绝对百分误差(mape)分别为5.47和9.48％，计算公式如下：

其中，rmse是均方根误差，mape是平均绝对百分比误差，n是样本数，yi是测量值，yi'是估算值。

步骤十一、依据tsi反演模型，可以根据(rhown4/rhown5,rrs(665)/rrs(704)；rhown1/rhown4,rrs(443)/rrs(665))来估算湖泊富营养化指数tsi。

步骤十二、为了验证tsi反演模型的准确性，基于步骤s9获取的验证数据组，进行tsi实测值与依据本专利方法估算出来的tsi值进行拟合分析(图4)，结果表明均方根误差(rmse)和平均绝对百分误差(mape)分别为7.36和12.61％，验证精度较好。

步骤十三、所述步骤十二中，基于步骤九获取的验证数据组，56.9％的采样点tsi值为50-100富营养，40.9％的采样点tsi值为30-50中营养，2.2％的采样点tsi值为0-30贫营养。

步骤十四、基于所述步骤十一获取的富营养化指数tsi反演模型，获取典型湖泊查干湖春季(4月)、夏季(8月)和秋季(10月)的msi影像，应用本专利发明的tsi反演模型，计算得到查干湖的tsi的时空分布(图5)。

依据本专利方法计算出来的tsi具有较高的可信度。湖泊富营养化指数tsi和rhown4/rhown5(rrs(665)/rrs(704))、rhown1/rhown4(rrs(443)/rrs(665))的相关性模型是本发明首先确立的。

本发明方法通过构建湖泊富营养化tsi指数反演模型，确定了一种新的、便捷的估算湖泊富营养化指数方法，即只需提取湖泊哨兵2号msilevel1c产品rhown4/rhown5(rrs(665)/rrs(704))、rhown1/rhown4(rrs(443)/rrs(665))反射率比，就可通过该模型计算出湖泊的富营养化指数tsi值，进而针对湖泊富营养状况评价。

本发明所述的湖泊富营养指数tsi反演模型，具有较高的反演精度，均方根误差在本实验中可达到5.47，平均相对误差为9.48％，相对于现有技术精度显著提高；另一方面，msi传感器是搭载在哨兵sentinel卫星，自2015年发射，较高的时间分辨率(2-5天)和空间分辨率(10,20和60m)使其在内陆湖泊水生态系统监测具有广阔应用前景；本发明所述方法的湖泊富营养指数tsi反演模型，是基于中国内陆湖泊大尺度实时监测数据建立，在内陆湖泊水生态系统健康监测和评估的通用性方面具有更好的潜力，极大地提高湖泊富营养化评估的效率和规模。

本发明将有助于评估湖泊水体富营养状态的时空变化及其发展趋势，可有效评估湖泊富营养状态和生态修复的绩效，对分析湖泊生态系统的结构和功能、预测湖泊富营养化状态和制定防控政策具有重要意义，为政府部门的水环境保护、水生态资源管理的科学决策提供科技支撑。