梯度提升回归树在千岛湖水体CDOM反演中的应用

金则澎，毛峰，程乾，李军，张轩豪

(1.浙江工商大学 旅游与城乡规划学院，杭州 310018；2.杭州电子科技大学 管理学院，杭州 310018；3.重庆师范大学 地理与旅游学院，重庆 400047)

0 引言

近几十年来，由于我国经济的飞速发展导致湖泊资源开发利用不断增加，在人为因素和自然因素的影响下，湖泊水环境逐年恶化[1]。淡水湖泊作为人类重要的饮用水源，也面临着农业、气候变化和人类活动带来的污染问题[2]。有色可溶解有机物(colored dissolved organic matter，CDOM)是溶解有机物(dissolved organic matter，DOM)的生物活性成分，CDOM在400～500 nm处吸收高达90%的水下太阳辐射，可以保护水下生态系统免受有害紫外线辐射[3-4]。但是，高浓度的CDOM对太阳光的吸收能力较强，会抑制水下植物的光合作用，这对淡水生态系统的初级生产产生了重大影响[5-6]。而且CDOM在水中具有吸附特性，会影响水污染物的分布与迁移[7]。此外，CDOM对饮用水处理过程也会带来许多负面影响，例如CDOM会降低氧化剂和消毒剂的效力，并在氧化过程中产生不良的副产物[8]。

了解淡水中CDOM的来源、浓度和循环对于管理水生资源和预测环境变化非常重要。监测湖泊中的CDOM，研究CDOM的时空分布能有效评价水质，并制定出保护水资源的计划。水质监测分为传统理化监测、生物监测与遥感监测三类，常规水质监测的效率不高，而且不能显示湖泊整体的水质情况[9-10]。目前人们已经开发了不同类型的检索算法来反演水质CDOM，如经验法、半解析法、矩阵反演法和机器学习法。然而，这些模型并不能很好地适用于复杂水域，特别是低浓度CDOM的内陆水域。最近，Chen等[11]用休伦湖萨吉诺湾的数据对其中15种CDOM检索方法进行了比较，发现大多数算法在CDOM浓度较高时往往会低估其浓度，在CDOM浓度较低时往往会高估其浓度。因为蓝色光谱范围内的色胺吸收与由叶绿素a(chlorophyll a，Chl-a)和总悬浮物质(total suspended matter，TSM)引起的吸收是重叠的[12]。因此，迫切需要开发一种适用于计算低浓度水域CDOM值的算法。自上个世纪以来，人工神经网络(artificial neural network，ANN)就被用于CDOM反演，并显示出处理复杂内陆水域的潜力。同时，基于机器学习的方法也被用于CDOM反演，并获得了良好的效果[13]。目前已有很多基于机器学习的方法，如随机森林回归(random forest regression，RFR)、梯度提升回归树(gradient boosting regression tree，GBRT)、支持向量机回归(support vector regression，SVR)和高斯过程回归(gaussian process regression，GPR)，但哪一种更适合淡水湖泊中的低浓度CDOM检测尚不得而知。GBRT的优点是它可以拟合复杂的非线性关系，但不需要事先进行数据转换或消除孤立点扫描[14]。因此，在本研究中，对GBRT模型进行了测试，并将其与其他几种机器学习方法进行了比较，结果表明GBRT在测试的模型中确实是最好的。

一些面向海洋的海洋水色传感器，具有粗糙的空间分辨率，不适合在大多数河流、湖泊进行遥感应用。然而面向陆地的卫星传感器，如Landsat系列，具有很好的空间分辨率(Landsat-8为30 m)，一些基于Landsat-8的波段比模型已被提出用于低浓度的CDOM(CDOM在440 nm的吸收系数在0.066～1.242 m-1范围内)监测。但是Landsat系列卫星的光谱分辨率有限，不能很精确地确定最优波段比。因此，研究人员更加关注GF-5传感器，GF-5卫星影像涵盖了可见光至短波红外(400～2 500 nm)波段，空间分辨率为 30 m，在可见光波段的光谱分辨率可达5 nm，能提供更加精确的波段，有利于确定更加精确的波段比模型。

本文的目标是：探索GBRT模型和GF-5卫星影像反演千岛湖低浓度CDOM的能力；描绘2019年每个季度千岛湖水质CDOM浓度的时空变化；讨论影响千岛湖水质CDOM时空格局的环境因子。为实现这些目标，将2018年4月至2019年9月从千岛湖收集的实测数据用于GBRT算法。在2019年，利用四个GF-5卫星影像和千岛湖水文和气象数据研究CDOM与降雨量、气压、风速等环境因素的耦合关系。

1 数据和方法

1.1 研究区域

千岛湖(图1)位于浙江省西部和安徽省南部，水域面积为580 km2，有1 078个岛屿，面积大于2 500 m2，平均深度34 m，最大容量19 267 m3，如果储水量处于108 m的正常值，则流域面积为10 480 km2[15]。记录的年平均气温为17 ℃，年平均降水量为1 636.5 mm(1961—2014年)。千岛湖水域内的用地类型以林业用地为主，其次是农业和城镇建设用地(图2)。千岛湖水域水质优良，水体透明度高，是中国长江三角洲地区的主要饮用水源，为周围至少一千万人口提供了饮用水。然而，人们对千岛湖水域的不合理开发引发了大量的水土流失，随着降雨流入湖泊的工农业和生活废水等污染物的增多，千岛湖水域的水质情况出现下降的趋势，湖泊水环境问题越来越突出[16-17]。

根据以往研究的结果，将千岛湖划分为五个次区域，以揭示千岛湖不同水生环境中CDOM的时空变化：东北、西北、西南、湖泊的中部和东南部地区。以千岛湖水域周边5 km为缓冲区，进行了土地类型分类(图2)。

1.2 实地测量

在2018年4月19日、6月18日、2019年9月8日的三次实地取样活动中，共从千岛湖采集了40个水样(图1)。使用标准取水瓶在表面(深度约10 cm以下)收集水样，在环境水温下保存在琥珀瓶(聚丙烯250 ml)中，24小时内运往实验室分析，测量浓度等水质参数。在每个采样位置，使用ASD(analytical spectral devices)测量地表水光谱。FieldSpec地物光谱仪(波长范围为325～1 075 nm，间隔为1 nm)在观测点测量了水面辐亮度(Lt)、天空光辐亮度(Li)和水面下行辐照度(Ed)。指向水面的Lt传感器的观测角是太阳高度角40°、太阳方位角90°方向。指向天空的Li传感器采用与Lt传感器相同的太阳方位角，观测的太阳天顶角是40°。为了减小不确定度，每个取样点同时测量了十条光谱，并选择光谱平均值来计算遥感反射率(Rrs)，表达如式(1)所示。

(1)

式中：水表面反射因子ρ的值受风速影响，在风速可忽略的平静天气下取值为0.022，风速5 m/s时取值0.025；风速10 m/s时，取值0.026～0.028。本研究考虑到千岛湖的地理条件和水文特征，设置ρ=0.024 5。

在实验室低压条件下(<5ATM)，通过GF/F玻璃微纤维膜(0.45 μm)对地表水样品进行过滤，以测定表面水样的色度吸收。通过Cray-60分光辐射度计测量了200～800 nm的CDOM吸光度，测量结果为1 cm长，基线校正为Milli-Q baseline。最终，吸收系数的计算如式(2)所示。

(2)

式中：A(λ)是Cary-60分光辐射度计测量的CDOM吸光度；Path是用过的试管的路径长度(1 cm)。在这项研究中，CDOM浓度是通过其在440 nm处的吸收系数来参数化的。CDOM浓度如表1所示。

表1 千岛湖2018—2019年CDOM(440)野外实测数据统计

为了使GF-5影像导出aCDOM(440)，使用其光谱响应函数RSR(relative spectral respone)将研究区实测光谱模拟出各个波段光谱，计算如式(3)所示。

(3)

式中：Rrs(Bk)是GF-5影像的第k(k=2，3，4，…，n)个波段的遥感反射率，第k个波段的Rrs(Bk)是由λm到λn波段范围的Rrsfield积分得到;Rrsfield是现场测量的光谱。

1.3 水文和气象数据

分别从屯溪和渔梁两个水文站以及中国气象数据网(http：//data.cma.cn/)收集了2019年千岛湖降雨量、气压、风速等气象数据。

1.4 GF-5影像数据获取和处理

GF-5卫星以高空间分辨率(30 m)，高光谱分辨率(5～10 nm)的方式，对地球表面进行观测。2019年，共采集了六张无云GF-5卫星影像，开发了CDOM反演算法，并对千岛湖的CDOM时空变化进行了分析。使用了一张与野外实测日期相同的2019年9月8日的GF-5卫星影像验证算法。利用ENVI5.3中的大气校正模块(FLAASH)进行大气校正，得到表面反射率(Rt)。在此基础上，采用5*5均值滤波，降低了图像的不确定性。ENVI中执行FLAASH大气校正如式(4)所示。

(4)

式中：L是传感器接收的像素辐亮度值；ρe是像素及其周边像素的平均表面反射率；S是天体反照率；La是被大气散射进入传感器的辐亮度；系数A和C分别取决于大气和地理位置条件，与表面无关；大气校正后输出的是表面反射率影像Rt。最后计算遥感反射率(Rrs)如式(5)所示。

(5)

式中：Lr(θ，φ)是天顶角θ和方位角φ观测角下的水表面辐亮度；Ed是向下辐照度。这两个参数使用Hydrolight估算。计算这两个物理量的过程中，需要在Hydrolight软件中输入一些参数，如太阳天顶角、影像拍摄时间、地理位置、风速和云覆盖百分比等。这些参数可从GF-5影像数据的头文件和国家气象中心提供的数据中获取。其他需要输入的参数(例如水质参数浓度和水深值)对Lr(θ，φ)和Ed的大小无影响而被设置成缺省值或是任意值。

2 算法实现

2.1 算法背景

集成学习算法是一种强大的机器学习技术，已成功地用于各种分类和回归问题，以提供更准确的响应变量估计[18]。在本研究中选择GBRT模型，因为它能够适应复杂的非线性关系，而且不需要事先的数据转换或消除孤立点扫描。GBRT模型使用CART树作为弱分类器，需要多次迭代。新生成的回归树在每次迭代中都会拟合上一棵树的误差，采用梯度下降法在每次迭代中向损失函数的负梯度移动，使损失函数下降。一般来说，GBRT模型的每一次迭代都会产生一个弱分类器，其精度不高，但集成弱分类器可以达到更高的精度。每个回归树的预测结果的加权和就是预测值(图3)。

图3 GBRT模型构建过程

GBRT模型表达如式(6)所示。

(6)

式中：n是弱学习器的数量；θi是系数(减少过拟合)；fi是弱学习器；Fm是最终的一般模型。GBRT模型有几个重要的参数：每个弱学习器的最大深度(一般不超过五个)、弱学习器的最大数量，以及学习率(更高的学习率意味着更强的修正，使模型更加复杂)。所以需要选择适合的参数来提高模型的性能。

2.2 算法架构

涉及蓝、绿和红域遥感的波段比算法已被广泛用于淡水湖的CDOM遥感[19]。此外，与使用单波段反射率相比，波段比模型可以减少大气校正的更多不确定性[20]。本研究以B25(Rrs493)、B41(Rrs561)和B56(Rrs625)组成的GF-5影像波段比作为GBRT模型的输入，通过决定系数(R2)、均方根误差(RMSE)、偏差和平均绝对百分比误差(MAPE)四个指标确定最佳组合。

现场测量的40个样本被随机分成两个独立的数据集，它们的样本比例约为3∶1，这是机器学习中常用的比例。因此，GBRT的训练和验证数据集分别包含30个和10个现场测量样本。基础学习者的数量为40人，学习率为0.1，每个学习者的最大深度为5，其他参数默认配置。

2.3 与其他CDOM检索算法的比较

为了验证GBRT在估算淡水环境中低浓度CDOM方面的优势，还对以往的两种传统模型进行了比较。对数模型，参见式(7)，其中aCDOM(440)的范围为0.51～25.1 m-1[21]。另一个模型是一元二次多项式模型，见式(8)，其中aCDOM(440)的范围为0.644～1.413 m-1[22]。第二个模型在水体类型和CDOM范围方面与构建模型更具可比性。

(7)

(8)

式中：B25、B41和B56是GF-5影像数据的波段。对于不同的研究地点，参数D、E和F需要重新校准。

3 结果分析

3.1 野外测量的CDOM和地表光谱

在以前的许多研究地点，内陆水域中的CDOM浓度范围很大。例如，在明尼苏达州的15个湖泊中，aCDOM(440)的范围为0.6～19.4 m-1[23]。与这些CDOM浓度变化较大的湖泊相比，2019年千岛湖的CDOM浓度变化相对较低。野外实测的aCDOM(440)范围为0.103～1.366 m-1，平均值为0.53 m-1。野外实测数据显示出明显的季节变化，4月份的aCDOM(440)在0.103～0.507 m-1(平均0.253 m-1)范围内，而9月份的变化范围为0.196～1.366 m-1(平均0.637 m-1)，说明秋季的CDOM浓度高于早春。

千岛湖实测的遥感反射率如图4所示，显示了复杂内陆淡水的典型光谱特征。浮游植物对叶绿素和胡萝卜素的低吸收和低散射产生了570 nm的峰值[24]。670 nm处的小反射谷可能是由叶绿素a对红光波段强烈吸收引起的[25]。在千岛湖，由藻类色素和水体自身吸收较小引起的700 nm附近的叶绿素荧光并不显著。

图4 千岛湖实测地表光谱(Rrs)

3.2 模型评估和比较

为了确定GBRT模型的最佳输入，对所有可能的波段比组合进行了计算。考虑四个指标(R2、MAPE、RMSE、偏差)，最终选择了B25/B41和B56/B25来反演CDOM浓度。最佳的GBRT模型的R2为0.90和0.93，MAPE为0.34%和15.53%，RMSE为0.35 m-1和0.16 m-1，培训和验证数据集的偏差为-0.004和-0.032，测量到的训练和验证数据集最佳模型的CDOM浓度散点图如图5(a)和图5(b)所示。结果表明，该算法的性能可以用于反演千岛湖的CDOM浓度。

图5 模型验证结果

本研究利用在千岛湖的实测资料，对前人建立的两个模型进行了检验，研究发现对数模型(R2=0.357 4，RMSE=0.394 5 m-1)(图6(a))和多项式模型(R2=0.048 4，RMSE=0.414 7 m-1)(图6(b))。在千岛湖不能很好地拟合。将GBRT与RFR、SVR、GPR进行了比较，所有模型均选择波段比B25/B41和B56/B25作为输入。表2显示了统计分析的结果。结果表明，支持向量回归对训练数据集(R2=0.44，MAPE=64.5%)和验证数据集(R2=0.27，MAPE=68.3%)都表现不佳，而高斯过程回归对验证数据集产生了过拟合(R2=0.02，MAPE=86.3%)。与随机森林回归相比，对于验证数据集，梯度提升回归树的R2增加了20%，而MAPE降低了26.17%。

图6 检验千岛湖低浓度CDOM反演的两个经验模型

表2 四种典型机器学习算法的比较

3.3 CDOM的时空变化

将训练好的GBRT模型应用于2019年的四幅GF-5影像，以估计千岛湖的CDOM浓度，如图7所示。为了进一步分析CDOM的时空变化，根据每个区域的有效像素估计了各个子区域每个季度的平均aCDOM(440) 浓度(图8)。结果表明，2019年千岛湖CDOM的季平均浓度较低，为春季0.61 m-1、夏季0.65 m-1、秋季0.75 m-1和冬季0.54 m-1，CDOM浓度季节性差异较为显著，其中秋季明显高于其他三季。四幅反演图像中，CDOM浓度高值区与千岛湖流域用地分类图(图2)中的建筑/道路、农田用地相契合，表明了人类活动对于水体CDOM浓度影响较大。千岛湖西北区和西南区的CDOM含量明显高于其他地区。西北区的新安江占千岛湖总径流量的60%，西南区的武强溪占千岛湖总径流量的25%。众所周知，径流会将陆地腐殖质带入湖泊，使得大量的CDOM通过武强溪和新安江汇入千岛湖，西北和西南区河口的水流流速比中部地区的水流流速大，河口区域的水体充分混合，使得水体分层不明显，因为紫外光无法轻易穿透表层水体，因此光漂白速率会比那些分层明显的中部区域低，从而导致西北和西南区的CDOM浓度相对较高。

图8 千岛湖季平均aCDOM(440)浓度

3.4 CDOM与水文气象因子的关系

以前的一些研究表明，降水量对水中CDOM浓度有重大影响[26]。CDOM浓度与降水量的R2为0.620 8，表明降水量对千岛湖CDOM浓度的变化影响显著。千岛湖地区年间降水量不均匀，水位有明显落差，其平均水深为37 m，最大水深为90 m。丰水期的降水量比枯水期的降水量高得多，高度动态变化的水位会使得湖泊底质变化，而水体底部有机质的释放是水中CDOM的来源之一。丰水期水流速度大，河流中的CDOM会流入主湖区；而枯水期水域的水位低，水流速度不大，从河流流入的CDOM就会减少。夏季和秋季正是千岛湖的丰水期，因此其CDOM浓度高于春季和冬季。然而，在本研究中，千岛湖的CDOM浓度峰值出现在秋季，而不是在降雨量最大的夏季(图9)。

图9 千岛湖季平均降雨量aCDOM(440)之间的相关关系

众所周知，光化学降解和光漂白是去除CDOM的重要机制。例如，在太湖进行的一项为期12 d的实验发现，当CDOM暴露在太阳辐射下时，CDOM的吸收量aCDOM(355)和aCDOM(280)分别下降了29.8%和20.8%[27]。2019年，千岛湖在夏季(5—7月)经历了高强度的UV-B辐射，这显著增加地表水的光漂白强度，从而消除部分降雨和径流带来的CDOM。因此，2019年夏季千岛湖反演的CDOM浓度并不太高。相比之下，秋季的CDOM浓度较高可能是由于以下三个原因。①经过连续多个月的降雨后，地表径流和地下水逐渐将丰富的土壤有机质释放到水中。②秋季太阳辐射减弱，减缓了CDOM的光漂白。③秋季越来越多的落叶腐烂，直接增加了土壤和水中的有机质含量；冬季是一年中降雨量最少的季节，加上温度下降导致的水面结冰，径流量大大减少，CDOM浓度处于一年中最低；春季由于降雨量的回升，以及冰雪融化加大了地表径流，导致外源有机物的输入加大，CDOM浓度相比较冬季有所上升。

分析卫星获取的CDOM浓度与水文和气象因子(气压和风速)之间的相关性发现，CDOM浓度与气压和风速的R2分别为0.082 5和0.096 4(图10)，表明这些因素对千岛湖CDOM浓度的变化没有显著影响。

图10 千岛湖季平均气压、风速以及与aCDOM(440)之间的相关关系

4 结束语

利用梯度提升回归树算法，结合GF-5高光谱卫星影像建立了千岛湖低浓度CDOM遥感模型。经比较验证，模型性能良好，RMSE=0.1 m-1，MAPE=19.0%。然后，应用GBRT模型预测了2019年千岛湖CDOM浓度的时空变化。千岛湖整体aCDOM(440)较低，最低为0.005 m-1，最高为1.472 m-1，CDOM浓度季节差异较为显著，秋季>夏季>春季>冬季。各个季节高CDOM浓度部分一般位于湖的边缘区，主要是入水口以及湖周围与人类活动密切相关的地方，如西北和西南区的入水口附近的CDOM浓度较高。最后结合水文数据，发现CDOM浓度与气压和风速的相关性很低(R2分别为0.082 5和0.096 4)，表明这些气象因素对千岛湖CDOM浓度的变化影响不大。

由于现场测量和观测仅在一年内进行，更多的环境因素或事件可能会改变CDOM和水质，因此GBRT模型仍需进一步改进。本研究表明，GBRT模型和GF-5影像能准确监测低浓度湖泊中CDOM的时空变化，可以更好地了解河流、湖泊系统中可溶解有机质与其环境因子之间的耦合关系。