利用高分一号卫星与XGBoost模型的水体总氮和总磷监测技术

赵力，卢修元，谭海，马天浩

(1.四川农业大学 水利水电学院，四川 雅安 625014；2.自然资源部国土卫星遥感应用中心，北京 100048；3.辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000)

0 引言

长期以来农业面源污染形势严峻，而化肥施用量高居不下是主因。2014—2018年，全国单位耕地化肥施用量由0.83 t·ha-1增长为0.93 t·ha-1[1]，化肥中的氮、磷营养盐进入水体，导致水体富营养化，影响水生环境，危害人体健康及动植物的生存，阻碍当地社会经济可持续发展。因此，快速监测总氮(TN)、总磷(TP)的含量，客观地反映水体富营养化程度，为污染水体治理提供依据是十分必要的[2]。传统的水体氮、磷浓度监测以点位监测为主，由人工到现场采集水样并进行化学分析，该方法耗时费力，且样点稀疏，只能反映取样点位的水质情况，对未取样位置，无法准确判定其水质。而卫星遥感技术通过经验、半经验模型或是分析模型建立卫星遥感观测信息与地面点位实测数据间的关系[3]，可以高效监测水质。潘邦龙等[4]利用环境一号卫星(HJ-1)超光谱遥感数据，建立了TN的三波段回归克里格模型；王小平等[5]采用Landsat-8 OLI数据，建立了艾比湖流域TP的4波段反演模型；钱彬杰[6]利用Landsat-8 OLI卫星遥感影像，建立了氨氮和TP波段反射率比值的回归模型；王磊等[7]基于Landsat-8 OLI影像，对丹江口水库的TP浓度进行了反演；林剑远等[8]采用芬兰机载高光谱成像系统获取的遥感数据，运用波段比值对嘉兴市水体的TP和TN进行反演，为内陆城市河网水质监测提供了参考。上述学者对水体氮、磷盐的遥感监测提供了有益探索。然而，由于不同建模方式的原理及反映复杂映射关系的能力不同，导致模型预测的效果差异较大，实测值与预测值的决定系数R2低至0.50，高至0.99[9]。此外，目前国内外遥感水质监测研究主要利用中低分辨率遥感数据，多数研究仅针对较大的内陆湖泊或海洋等水体[10]。

利用遥感进行水质监测中的建模方法，以数学统计法最为常见[11]，但因其难以反映复杂的水体光学性质，导致所得模型稳定性较差[12]。近年来，以BP神经网络为代表的机器学习因为具有自适应、自学习、高效率和容错性等优点，在遥感领域也得到了广泛使用[13-15]。而XGBoost(extreme gradient boosting)模型是提升算法(boosting)的最新研究成果之一，不仅具有较强的自我学习能力和映射能力，同时还弥补了BP神经网络无法反映输入因子的重要性的问题，已被广泛用于图像处理[16]和大数据处理[17]等领域。就数据源而言，目前常用的MODIS、MERIS、Landsat、HJ-1等中低分辨率数据限制了研究尺度。而2013年我国发射升空的高分一号卫星(GF-1)搭载有高空间分辨率的宽幅多光谱相机(WFV)，为小尺度水体的遥感监测提供了基础条件[18-19]。尽管已有基于GF-1数据并采用机器学习方法反演水质的研究，但将GF-1数据与XGBoost模型联合使用进行水质监测的研究还鲜有报道。

本研究以长江上游农业生产水汇集的白水湖水库为研究区，以水体TN、TP为研究对象，分别采用XGBoost、BP神经网络2种机器学习方法，构建水质遥感反演模型，并与数学统计模型比较，旨在选出适用于研究区域的水质监测模型，为长期、便捷、大范围的水质监测提供技术支撑。

1 材料与方法

1.1 研究区概况

白水湖水库位于长江上游，川西平原北部(31°30′55″N～31°32′10″N，104°14′45″～104°16′07″E)，属中亚热带湿润季风气候，冬干春旱，夏季旱涝交错，秋多阴雨，水域面积2 km2，库容1 672万m3，平均水深10 m，湖区四面环田，在雨水冲刷和地表径流作用下，农田里过量营养盐直接进入水体，造成水环境污染。

1.2 技术路线

采用ENVI 5.3软件校正GF-1 WFV1图像，提取研究水域4个波段的反射率。基于实测TN和TP数据，通过计算水质参数与波段及波段组合的Pearson相关系数，识别各个水质参数的敏感波段或波段组合，将其作为模型的输入因子。随机选取70%的样本作为训练集，剩余30%样本作为测试集，利用Python 3.7软件构建XGBoost模型，利用MATLAB 软件构建数学统计模型和BP神经网络模型，以R2和均方根误差(RMSE)为依据，比选模型。反演TN、TP，运用ArcGIS 10.2软件制作水质专题图，并分析TN和TP空间分布特征。具体技术路线如图1所示。

图1 技术路线

1.3 数据采集及预处理

1)水质参数测量。2019年8月20日，在白水湖采集表层水样87个，样点分布如图2所示。TN采用过硫酸钾氧化紫外分光光度法测定(HJ 636-2012)，TP采用钼酸铵分光光度法(GB 11893-89)测定。

图2 研究水域样点分布

2)卫星遥感影像数据获取与处理。GF-1 WFV影像空间分辨率为16 m，重访周期为4 d，具有B1(0.45～0.52 nm，蓝)、B2(0.52～0.59 nm，绿)、B3(0.63～0.69 nm，红)、B4(0.77～0.89 nm，近红外)4个波段。本文采用8月11日GF-1 WFV准同步数据。利用ENVI软件对GF1-WFV1数据进行预处理。DEM数据为GDEMV2 30 m。采用3次卷积插值完成正射校正，以15 m空间分辨率的Landsat-8 OLI影像为基准图像完成影像的几何精校正，辐射定标采用2018年GF-1 WFV1绝对辐射定标系数，FLAASH模块进行大气校正，部分输入参数如表1所示。使用归一化水体指数(normalized difference water index，NDWI)(式(2))对研究区水体进行掩膜处理，得到水体影像。

(1)

式中：Green为GF-1数据绿光波段反射率；NIR为近红外波段反射率。

表1 FLASSH 大气校正模块输入值

2 数据分析

2.1 模型输入因子分析

卫星传感器上记录的与水体信息有关的总辐射亮度主要包括水面散射辐射、水中散射辐射以及底部反射辐射3个部分。当水体污染物、水面粗糙度或微波发生变化时，水体散射和反射状态等也会变化，总辐射亮度随之改变，导致总辐射亮度与水质参数改变。正由于水体成分和水体状态会影响水体的光谱特征[20]，可以通过波段组合消除或减少与水质无关的信息，保留或放大显著相关的信息，提高反演模型精度。如通过剔除对有机污染物浓度不敏感的近红外波段，利用敏感度高的可见光波段建立模型，可以提高水体有机污染物反演精度[21]。此外，已有研究也表明，波段组合不仅可以消除部分大气影响[22]，还可以减少水面粗糙度变化的影响[23]，同时也能消除微波变化导致的吸收系数和后向散射系数干扰，减少水体反射率二向反射等问题[24]。所以，可以通过波段组合，减少污染物、粗糙度、微波变化等的干扰，从而把目标水质参数的特性表现出来。

将校正后的波段反射率信息经过单波段、双波段和三波段间的减法(如B1-B2、B3-B4等)、比值法(如B1/B2、(B1-B3)/(B2-B4)等)、对数法(ln(B1)、ln(B1/B2) 等)，组合成254个波段反射率信息，编号为Cn(n=1～254)，利用SPSS 21软件分别计算各波段组合信息与TN、TP的Person相关系数，结果如图3所示。取极显著相关(p<0.01)且相关系数最大的波段组合作为输入因子，则TN的输入因子为波段组合C34、C33、C45、C31，TP为C156、C47、C174、C144，详见表2。由表2可知，TN对B2、B3和B4的波段组合最敏感，TP对B1和B2的波段组合最敏感。Koponen等[25]的研究也表明，采用多波段组合可以在一定程度上消除水体污染物、表面粗糙度以及微波的时空变化等干扰，与单波段相比，更能真实反映水质参数对光谱曲线的影响。本文的结论与Koponen的研究相类似。

图3 波段组合信息与TN、TP 的Person相关系数

表2 模型输入因子

2.2 TN反演模型的建立与评价

1)XGBoost。XGBoost是一种基于决策树的监督模型，是集成学习方法的一种，由华盛顿大学的Chen等[26]提出，其核心在于新的模型在相应损失函数梯度方向建立，修正残差的同时控制模型的复杂度。本研究以相关性最高的波段组合C34、C33、C45和C31作为输入因子，以实测TN作为输出，XGBoost模型基于树模型进行提升计算，采用线性回归目标函数，多次调参对比结果，最终设定XGBoost最大深度为3，限定孩子节点中最小的样本权重和为3，详细参数设定见表3。训练所得模型中波段组合的重要性C34>C45>C33>C31分别为32、31、23和0，包含波段B2、B3和B4，模型的决定系数R2为0.970，RMSE为1.24×10-2mg·L-1。

表3 XGBoost的TN反演模型参数

2)BP神经网络。BP神经网络由Webos在1974年提出，主要特点是信号向前传递，误差反向传递，调整层间的权值、阈值以达到预设精度要求[27]。BP神经网络构建中，隐含层神经元节点数的选择直接影响网络对复杂问题的映射能力[28]，隐含层结点数过少，则网络收敛速度减慢，精度降低，结点数过多，则导致计算量增加，模型过拟合，泛化能力降低。但目前尚无确定隐含层最佳结点数的直接方法[29]。因此，通常采用的方法是多次实验对比分析。本实验经过60个训练样本训练不同节点数的单隐含层BP神经网络模型，节点数3～10时，模型的R2和RMSE如表4所示，结果表现为9节点时R2最高，4节点次之。分别计算剩余27个测试样本在8种不同节点数的R2和RMSE，结果如表5所示，R2表现为4节点最高，7节点次之。综上，将单层4节点BP神经网络模型作为TN的BP神经网络最优模型，此时模型的R2为0.769，RMSE为3.47×10-2mg·L-1。

3)数学统计模型。取训练样本(n=60)，以波段组合反射率C34、C33、C45、C31为自变量，TN浓度为因变量，构建多元线性回归模型，模型中波段组合的权重为C31>C45>C34>C33，分别为0.518 12、0.144 13、0.138 03和0，模型的R2和RMSE分别为0.755和3.50×10-2mg·L-1，详见式(2)。

表4 BP神经网络的TN模型精度(n=60)

表5 BP神经网络的TN预测精度(n=27)

TNST=-0.138 03×C34-0.144 13×C45-0.518 12×C31+0.765 847

(2)

代入表2中因子，则反演TN的数学统计模型包含4个波段，如式(3)所示。

(3)

4)TN模型评价。以30%(27个)的采样点数据作为模型测试样本，将波段组合反射率分别输入XGBoost、单层4节点BP神经网络和数学统计模型，输出TN预测值，作预测值与实测值的散点图(图4)。由图4可知，XGBoost模型(XG)、BP神经网络模型(BP)、数学统计模型(ST)的R2分别为0.835 0、0.801 4、0.756 7；各模型的RMSE值分别为3.25×10-2、2.97×10-2、3.47×10-2mg·L-1。XGBoost模型的R2比BP神经网络和数学统计模型分别高了4.19%和10.35%，RMSE比数学统计模型降低了6.34%。由此可见，在3种模型中，XGBoost的TN反演模型具有最高的R2和较低的RMSE，对研究区TN浓度的监测效果明显优于BP神经网络和数学统计模型。

图4 TN预测值与实测值比较

2.3 TP反演模型的建立与评价

1)XGBoost。以相关性最高的波段组合C156、C47、C144和C174作为输入因子，以实测TP作为输出，XGBoost基于树模型进行提升计算，采用线性回归目标函数，对比多次调参结果，选定XGBoost最大深度为3，限定孩子节点中最小的样本权重和为1，详细参数设定如表6所示。训练所得模型中波段组合的重要性C47>C156>C174>C144，分别为96、65、45和8，包含波段B1和B2，模型的R2为0.896 6，RMSE为1.97×10-3mg·L-1。

表6 XGBoost的TP反演模型参数

2)BP神经网络。以随机抽取的70%的数据作为训练样本，计算不同节点数的单隐含层BP神经网络模型的R2和RMSE，结果如表7所示，R2表现为10节点>9节点>7节点>其他。以剩余30%的数据作为测试样本，分别计算8种不同节点数模型的R2和RMSE，结果如表8所示，R2表现为5节点>6节点>10节点>其他。因此，将单层10节点BP神经网络模型作为TP的BP神经网络最优模型，R2为0.911，RMSE为1.96×10-3mg·L-1。

表7 BP神经网络的TP建模精度(n=60)

表8 BP神经网络的TP预测精度(n=27)

3)统计学模型。以波段组合反射率C156、C47、C144和C174为自变量，TP浓度为因变量，抽取70%的数据作为模型训练样本，构建多元线性回归模型，模型中波段组合反射率的权重为C156>C47>C174>C144，分别为0.354 57、0.146 55、0.026 76、0.022 17，与XGBoost模型的重要性存在差异，模型的R2和RMSE分别为0.696和3.49×10-3mg·L-1，模型如式(4)所示。

TPST=-0.355 7×C156-0.146 55×
C47-0.022 17×C144+0.026 76×
C174-0.084 13

(4)

代入表2中因子，则反演TP的数学统计模型包含B1和B2波段，模型为式(5)。

(5)

4)TP模型评价。将训练样本(n=27个)的波段组合反射率分别输入XGBoost、单层10节点BP神经网络和线性数学统计模型中，输出TP预测值，作预测值与实测值的散点图(图5)。由图5可知，XGBoost模型(XG)、BP神经网络模型(BP)、数学统计模型(ST)的R2分别为：0.896 6、0.754 2、0.746 2；各模型的RMSE值分别为1.97×10-3、0.974×10-3、2.67×10-3mg·L-1。XGBoost模型与BP神经网络和数学统计的R2相比，分别高了18.88%和20.16%。与数学统计模型相比,XGBoost模型的RMSE低了26.22%。由此可见，XGBoost模型能有效提高预测精度，对研究区TP浓度具有更好的反演效果，更适用于研究水体的TP监测。

2.4 氮、磷遥感评价模型应用

以水质反演效果最好的XGBoost模型为基础，将对应反射率数据输入模型，反演得到TN、TP值。利用ArcGIS分别制作TN、TP的水质专题图，结果见图6。由图6(a)可知,TN浓度进水区域高于出水区域，近岸高于湖心，符合反硝化作用下的库区水体TN分布特征。由图6(b)可知，TP浓度分布与TN相反，呈现出水区域高于进水区域，湖心高于近岸的趋势，原因可能是水体中磷酸盐被底泥吸附，吸附量随含氧量降低而降低，出水口和湖心底泥含氧量较低，磷酸盐释放，导致水体TP增加。

图5 TP预测值与实测值比较图

图6 TN、TP空间分布

3 结束语

本研究以农业生产水汇集的白水湖水库为研究区，利用校正后的GF-1 WFV影像，采用XGBoost模型和BP神经网络模型2种机器学习方法，反演水体TN和TP，并与数学统计模型比较，得到以下结论。

1)3种模型中，XGBoost模型更适合预测TN、TP，其预测值与实际值的R2分别可达到0.835、0.897。理论上而言，提高布点均匀度还可以进一步提高XGBoost模型反演TN、TP的精度。

2)同一水质参数，采用不同建模方法建模时，其输入因子的重要性不同。TN差异较大时，XGBoost模型的C31最不重要，而数学统计模型中C31最重要；TP差异较小时，XGBoost模型中C47最重要，C156次之，而数学统计模型中C156最重要，C47次之，其他不变。对不同水质指标，敏感波段不同，TN的敏感波段为绿波、红波和近红外波段，TP的敏感波段为蓝波和绿波。前人研究表明，水深影响绿波段反射率，悬浮物浓度增加会产生“红移”现象，因此，可考虑结合水深和悬浮物浓度建立TN、TP反演模型，以进一步提高精度。

本研究结果表明，采用GF-1 WFV信息，结合XGBoost模型能准确反演水体TN、TP情况。本研究为农业生产水汇集区域的水体监测提供了方法，也为其他水体利用遥感技术进行水质的总氮、总磷监测提供了思路。