抚仙湖氨氮遥感估算研究

Study on Remote Sensing Estimation of Ammonia Nitrogen in Fuxian Lake

摘要: 遥感估算湖泊水质对了解湖泊水质的空间分布特征具有重要意义，本研究利用2013~2017年的Landsat 8 OLI和GF-1 WFV卫星数据，对抚仙湖的氨氮浓度开展遥感估算。使用反向传播神经网络(BPNN)和多元线性回归(MLR)模型开展交叉对比研究。研究结果表明：在样本数较小的情况下BPNN比MLR有绝对的优势，其验证误差RMSE和rRMSE分别为0.0246 mg/L和40.74%。MLR过拟合现象严重，虽建模线性较好，误差较小，但验证误差较大。将风速作为输入参数可有效提高建模的线性，但也会增加过拟合风险。本研究的将对贫营养湖泊水质的遥感估算有一定的参考和借鉴作用。

Abstract: Understanding the spatial distribution characteristics of lake water quality requires a quantitative remote sensing estimation of lake water quality. Landsat 8 OLI and GF-1 WFV satellite data from 2013 to 2017 are presented in this study. Multiple linear regression (MLR) and back propagation neural network (BPNN) models were employed. The results indicate that BPNN has an absolute advantage over MLR despite the limited sample size, and the verification error RMSE and rRMSE are 0.0246 mg/L and 40.74% respectively. Even though the linear modeling is accurate and the error is minor, the MLR over fitting phenomenon is serious because the verification error is large. Wind velocity as an input parameter can enhance the linearity of the model, but it also increases overfitting. The results of this study can be used as a reference for the remote sensing estimation of the water quality of oligotrophic lakes.

1. 引言

抚仙湖是我国面积最大的属于I类水质的深水型淡水湖泊，其蓄水量占云南九大高原湖泊的2/3以上，占全国淡水湖泊水资源的9.2% [1]。抚仙湖以其巨大的优质水资源储量，在中国具有重要战略地位。多年以来，抚仙湖面临着流域内面源污染的威胁，目前已从Ⅰ类水体下降到II类[2]。人们对抚仙湖水质变化的详细过程还未有充分的认识，主要局限于水质监测历史数据在时空分辨率的不足，导致局部湖区出现的水质恶化不容易被发现。近年建设的水质自动监测站可高频率获取水质数据，但鉴于高昂的开支，监测点较少，在空间分布上依然有很大的不足。

采用卫星对湖泊水质进行遥感估算可实现高时间频率，且空间上连续覆盖整个湖区的监测数据，可更详细掌握湖泊水质空间分布。基于卫星的水质遥感定量估算技术为解决上述问题提供了一条新的途径。目前常用多光谱卫星的空间分辨率达到10~30 m，重访周期2~5天[3]。抚仙湖主要污染物是氮磷等非光学敏感物质(nOAC)，然而目前对水质的遥感估算研究较为成熟的是对诸如叶绿素和悬浮颗粒物等光学敏感物质(OAC)的定量遥感监测[4]-[12]。目前鲜有针对清洁水体的非光学敏感物质遥感估算的报道，清洁水体的营养盐遥感估算比浑浊水体难，主要是由于贫营养水体营养盐含量很低，纯水对光线的吸收很强，卫星传感器灵敏度不足，容易造成较大误差[13]。

目前水质遥感估算模型主要有三类，即经验模型，半经验模型和物理模型，其中经验模型比较简单，在建模数据覆盖的时空范围内效果较好，但缺少普适性；物理模型在输入参数精准的情况下结果较为理想，但需要的输入参数很多且不易获取，较难实现；而介于两者之间的半分析模型也称之为生物光学模型，其依据生物光学特征，结合易于获取的光谱参数实现运算，其普适性相对较强，是许多光学敏感物质水色遥感的理想模型选择[5]。除了上述三种模型以外，在水色遥感领域还有许多关于机器学习与半分析模型对比研究的报道，由于半分析模型需要基于生物光学特征建立，而对于非光学敏感物质的遥感估算，无法使用半分析模型开展，因此诸如神经网络、遗传算法等机器学习模型成为比较理想的选择。机器学习算法基于灰色箱子的理论，单纯从数学角度探讨事物间相互联系，采用此类模型估算氮磷营养盐有望提高估算精度，减少误差[14]。

针对氮磷等非光学敏感物质的遥感估算研究较少，其遥感监测误差也较大，由于其对卫星光谱信息不敏感，这类物质的遥感监测的理论基础主要是基于叶绿素a与此类物质的相关关系，且相关研究主要针对内陆浑浊水体中的富营养化水体[15] [16] [17]。最近的综述文献认为，机器学习算法对nOAC物质遥感估算在普适性方面效果突出[18]，许多针对nOAC的遥感反演研究表明，与传统的经验模型对比，反向传播神经网络(BPNN)模型的效果最为突出，大部分结果的相对误差在 10%以内[8] [19] [20] [21] [22]。

本研究采用2013~2017年间65景高分一号(GF-1) WFV和17景Landsat 8 OLI影像估算抚仙湖NH₃-N含量，对比BPNN和多元线性回归方法的估算精度。本研究在数据处理前期发现部分完全无云卫星影像中湖泊水面上空有类似薄雾光谱特征，并提出因湖面大风导致表面反射率异常的科学假设，设计实验研究风速对模拟效果的影响。本研究将对贫营养湖泊非光学敏感水质参数的遥感反演有一定的参考和借鉴作用。

2. 研究方法

2.1. 数据的获取

本研究所用的水质监测数据来自中国环境监测总站公布的“全国主要流域重点断面水质自动监测周报”(http://www.cnemc.cn/sssj/szzdjczb/)中抚仙湖孤山站。气象数据来自于国家气象科学数据中心(http://data.cma.cn/)“中国地面气候资料日值数据集(V3.0)”，选取的气象站点为红塔，与孤山水质自动监测站距离为30km。遥感数据中的Landsat 8 OLI卫星数据通过美国地质调查局USGS网站(http://glovis.usgs.gov/)获取，GF-1 WFV数据通过中国资源卫星应用中心(http://www.cresda.com/CN/)获取。GF-1 卫星搭载了WFV1、WFV2、WFV3和WFV4，共4个WFV传感器，每个WFV传感器拥有蓝、绿、红、近红外4个波段，空间分辨率为16 m，而Landsat 8的OLI传感器拥有8个30 m分辨率的波段(表1)。

Table 1. Band parameters of satellite data used in this study

表1. 本研究所用的卫星数据波段参数

Landsat 8 OLI				GF-1 WFV
波段编号	波段	波长(um)	中心波长(um)	波段编号	波段	波长(um)	中心波长(um)
B1	气溶胶	0.43~0.45	0.443	−	−	−	−
B2	蓝	0.45~0.51	0.48	B1	蓝	0.45~0.52	0.485
B3	绿	0.53~0.59	0.585	B2	绿	0.52~0.59	0.555
B4	红	0.64~0.67	0.655	B3	红	0.63~0.69	0.66
B5	近红外	0.85~0.88	0.865	B4	近红外	0.77~0.89	0.83
B6	SWIR1	1.57~1.65	1.61	−	−	−	−
B7	SWIR2	2.11~2.29	2.2	−	−	−	−
B9	Cirrus	1.36~1.38	1.37	−	−	−	−

与卫星影像对应的水质监测数据选取卫星拍摄日期所在的那一周的水质监测周报数据，如果水质监测周报数据有缺失，则选取相邻的一周的监测数据，卫星影像提取的表面反射率与自动站监测数据获取的时间相差在一周以内。周报数据是一周中自动站监测数据的平均值，一周之内抚仙湖的水质波动很小，可以近似看做星地同步观测，近年来采用±7天内的卫星与实测数据开展遥感估算的研究也有报道[14]。

OLI和WFV传感器分别有17和65个观测窗口满足观测时间与氨氮实测时间对应。本研究选取分别与OLI和WFV传感器观测日期对应的水质监测周报数据中的氨氮实测数据与遥感数据以及周报日期范围对应的气象数据风速指标平均值组成数据集分别建模分析，分别命名为“OLI数据集”和“WFV数据集”，建模数据集基本参数见表2、表3。本研究所用的WFV和OLI传感器的空间分辨率分别是16m和30 m，提取的表面反射率必须为来自湖泊水面的纯净像元，不能含混合有陆地的信号，孤山自动监测站实际位置在孤山岛上，且为了让数据具代表性，本研究将与孤山站位置相邻的水域附近某个点为圆心，以500 m为半径的圆形范围水域作为表面反射率提取的区域，取该圆形范围的卫星光谱数据平均值用于后续的分析，见图1中的孤山站遥感估算区域。

Table 2. Basic information on datasets of OLI

表2. OLI数据集基本情况

参数	个案数	最小值	最大值	平均值	标准差
NH4-N (mg/L)	17	0.02	0.19	0.067	0.0414
ρs (B1)^*	17	1.93	545.71	216.008	130.308
ρs (B2)^*	17	38.6	557.16	240.911	133.291
ρs (B3)^*	17	32.9	503.1	243.420	120.107
ρs (B4)^*	17	3.26	380.46	126.706	107.163
ρs (B5)^*	17	2.09	337.79	102.482	97.391
ρs (B6)^*	17	6.46	270.97	80.740	79.144
ρs (B7)^*	17	11.82	244.48	66.833	71.034
ρs (B9)^*	17	6.62	46.26	12.185	10.753
平均风速(0.1 m/s)	17	9	39	18.177	8.248
最大风速(0.1 m/s)	17	29	79	54.294	12.648

^*ρs代表大气校正后得到的表面反射率，单位：10,000 * μW/(cm² * sr * nm)，B1~B9表示波段编号。

Table 3. Basic information on datasets of WFV

表3. WFV数据集基本情况

参数	个案数	最小值	最大值	平均值	标准差
NH4-N (mg/L)	65	0.01	0.1	0.0565	0.02503
ρs (B1)^*	65	67.24	1490.47	538.9415	331.4832
ρs (B2)^*	65	30.48	1073.79	440.954	262.7551
ρs (B3)^*	65	1	802.69	298.6928	213.5808
ρs (B4)^*	65	1	797.11	287.8076	200.1192
平均风速(0.1 m/s)	65	8	46	23.5077	7.8045
最大风速(0.1 m/s)	65	35	105	63.2615	13.88713

^*ρs代表大气校正后得到的表面反射率，单位：10,000 * μW/(cm² * sr * nm)，B1~B4表示波段编号。

Figure 1. Study area location

图1. 研究区域位置图

2.2. 氨氮遥感估算模型的训练和验证

本研究抚仙湖氨氮遥感估算模型的训练和验证过程分为遥感数据的预处理，模型的训练和验证，误差计算三个方面，本研究技术路线图如下图2所示：

Figure 2. Research technology roadmap

图2. 研究技术路线图

2.2.1. 遥感数据的预处理

在定量遥感研究中，其大气影响是造成误差的主要原因，上述研究常用卫星轨道高度通常在600~800 km区间，其遥感估算主要依赖传感器获得的离水辐射或表面反射率，而水体对光线有强烈的吸收，尤其是清洁水体吸收更加强烈，再加上卫星获取的辐射信息90%都来自于大气后向散射[23]，大气校正的难度极大。目前主流的大气校正流程化软件虽然采用了成熟的模型，例如ENVI的FLAASH模块基于MODTAN模型，但没有引入辅助数据，未发挥这些成熟模型的优势。本研究使用李国春开发的RSD软件的大气校正模块对WFV和OLI数据进行大气校正，该模块基于6S模型开发，可调用辅助数据进行大气校正，可充分发挥6S模型的优势[24]。

2.2.2. 模型的训练和验证

BP神经网络具任意复杂的模式分类能力和优良的多维函数映射能力，从结构上讲，BP网络具有输入层、隐藏层和输出层，从本质上讲，BP算法就是以网络误差平方为目标函数、用梯度下降法来计算目标函数最小值。BP神经网络计算过程由正、反向计算过程组成。正向传播过程，输入模式从输入层经隐单元层逐层处理，并转向输出层，每一层神经元的状态只影响下一层神经元的状态。如输出层不能得到期望值，则转入反向传播，将误差信号沿原来的连接通路返回，通过修改各神经元的权值，使误差信号最小。

本研究基于卫星的表面反射率数据，选用BP算法和多元线性回归方法来估算抚仙湖氨氮浓度，输入信号分两批，一批为OLI传感器8个波段的表面反射率和日平均风速、日最大风速，一共10个参数，17组样本。另一批为WFV传感器4个波段的表面反射率和日平均风速、日最大风速，共6个参数，65组样本，输出为NH₃-N浓度。

在BP神经网络分析中，因数据量较小，为保证验证样本数量，OLI数据选取10组数据作为训练样本，7组数据为验证样本，WFV数据选取44组数据为训练样本，21组数据为验证样本。因数据量级相差较大，对样本数据归一到[0.1, 0.9]。隐含层的传递函数选择logsig函数，输出层为purelin，训练方法为LM，建立BP网络，训练函数选用变学习速率梯度下降法。BP神经网络采用Matlab进行分析，多元线性回归自变量进入的方法采用输入法，在SPSS 24软件中分析。

2.2.3. 误差计算

本研究采用均方根误差(RMSE)，相对均方根误差(rRMSE)和相对误差$\left|{\epsilon }_i\right|$ 作为误差评价指标，其中RMSE是指估计值和观测值偏差的平方与观测次数n比值的平方根，能很好地反映出模型的误差，RMSE计算公式为：

$\text{RMSE}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^N{\left(X_{\text{esti},i}-X_{\text{meas},i}\right)}^2}}{N-1}}$ (1)

式中，$X_{\text{esti},i}$ 指估算值，$X_{\text{meas},i}$ 指的是实测值，N为总的数据量。

相对均方根误差(rRMSE)的计算公式为：

$\text{rRMSE}=\text{RMSE}/{\text{Mean}}_{\text{meas}}$ (2)

其中${\text{Mean}}_{\text{meas}}$ 为实测氨氮浓度平均值。$\left|{\epsilon }_i\right|$ 是测量值与估算值之差所占实测值的百分比的绝对值，它可表示每组验证数据的误差，公式如下：

$\left|{\epsilon }_i\right|=\left|\left(\frac{X_{\text{esti},i}-X_{\text{meas},i}}{X_{\text{meas,i}}}\right)\times 100\right|$ (3)

3. 结果与讨论

本研究分别对OLI数据集和WFV数据集建模和验证，对比BP神经网络和多元线性回归的建模和验证效果。表4为多元线性回归模型的输入数据和计算得到的回归方程。BP神经网络模型输入数据和隐含层节点数见表5。模型的线性和误差分析结果见图3。

Table 4. Input data and regression equation of linear regression model

表4. 线性回归模型输入数据和计算得到的回归方程

模型	输入数据	回归方程(*C为氨氮浓度)
OLI1	8个波段表面反射率	*C = 0.013421 + 0.003926 × ρB1 − 0.003542 × ρB2 + 0.000655 × ρB3 − 0.002024 × ρB4 + 0.002233 × ρB5 − 0.000803 × ρB7 − 0.003353 × ρB9
OLI2	8个波段表面反射率 + 日平均风速	*C = 0.078572 + 0.003433 × ρB1 − 0.004581 × ρB2 + 0.001808 × ρB3 + 0.00035 × ρB4 + 0.000147 × ρB5 − 0.000342 × ρB7 + 0.000252 × ρB9 − 0.006231 × 平均风速
OLI3	8个波段表面反射率 + 日平均风速 + 日最大风速	*C = 0.205786 + 0.002399 × ρB1 − 0.000688 × ρB2 − 0.002025 × ρB3 + 0.002487 × ρB4 − 0.001695 × ρB5 + 0.000553 × ρB7 + 0.004995 × ρB9 + 0.000939 × 平均风速 − 0.005332 × 最大风速
WFV1	4个波段表面反射率	*C = 0.057 + 0.000118 × ρB1 − 0.000247 × ρB2 + 0.000261 × ρB3 − 0.000137 × ρB4
WFV2	4个波段表面反射率 + 日平均风速	*C = 0.044929 + 0.000114 × ρB1 − 0.000228 × ρB2 + 0.000261 × ρB3 − 0.000163 × ρB4 + 0.000595 × 平均风速
WFV3	4个波段表面反射率 + 日平均风速 + 日最大风速	*C = 0.045718 + 0.000113 × ρB1 − 0.000227 × ρB2 + 0.000261 × ρB3 − 0.000164 × ρB4 + 0.000622 × 平均风速 − 0.000022 × 最大风速

Table 5. Input data and number of hidden layer nodes of back propagation neural network model

表5. BP神经网络模型输入数据和隐含层节点数

模型	输入数据	隐含层节点数
OLI-BP1	8个波段表面反射率	9
OLI-BP2	8个波段表面反射率 + 日平均风速	9
OLI-BP3	8个波段表面反射率 + 日平均风速 + 日最大风速	9
WFV-BP1	4个波段表面反射率	13
WFV-BP2	4个波段表面反射率 + 日平均风速	10
WFV-BP3	4个波段表面反射率 + 日平均风速 + 日最大风速	17

根据图3，在OLI数据集建模分析中，建模阶段，基于多元线性回归的OLI3模型建模线性最好，R²为1，建模误差最小，rRMSE为0.0925%，根据图4，|ε建模|最大值仅为0.43%，但是验证阶段误差很大，rRMSE达381.645%，而基于BP神经网络的OLI-BP3模型验证误差最小，rRMSE为83.52%。这两个模型的输入数据集一致，均包含气象数据，说明多元线性回归存在过拟合现象，不适合用于样本数较少的数据模拟。过拟合现象就是指模型严重依赖训练数据，导致建模线性好，误差小，但是验证误差大，导致模型无法应用。此外，验证数据集中实测值与模拟值散点图的线性回归方程R²对比中，OLI2大于OLI-BP3，说明OLI2验证数据线性更好，但是误差却仍然是OLI-BP3最低，从图2可以看出，OLI-BP3验证点比OLI2更加集中在趋势线周围，在图3中OLI-BP3模型验证数据绝对误差1/4位数在OLI数据集中最低，为24.57%，3/4位数略高于OLI-BP1模型，为116%。

Figure 3. Simulation results of MLR and BPNN analysis

图3. 多元线性回归与BP神经网络分析模拟结果图

在WFV数据集建模分析中，建模阶段，线性最好的是WFV-BP3模型，建模数据集的RMSE和rRMSE分别为0.0152和27.89%，|ε建模|1/4位数和3/4位数在WFV数据集中均最小，分别为4.5%和28.72%，但验证误差最大的也是该模型，而WFV-BP1的验证误差是最小的，验证数据集的RMSE和rRMSE分别为0.0246和40.74%，|ε验证|1/4位数和3/4位数在WFV数据集中均最小，分别为15.89%和48.88%，从图3可以看出，WFV-BP1的验证点分布比WFV-BP3更加集中在趋势线周围，说明模型的好坏不能只看回归方程的决定系数，也要综合考虑误差计算结果。对比前人的类似研究结果发现，其氨氮和总氮的估算误差在14%~17%之间，绝对误差为0.14~0.33 mg/L之间[8] [19] [22]，本研究的相对误差达到0%左右，但是绝对误差分布在0.01~0.018 mg/L之间，本研究相对误差高于同类研究，但绝对误差远低于同类研究，主要是因为本研究主要针对贫营养湖泊，其本底值较低，对信号噪声的敏感度增加，增加了遥感估算的难度，同时较低的真实值也客观上决定了其绝对误差会远低于浑浊水体的研究结果。

在OLI数据集建模分析中，将风速因子作为输入参数可有效提高建模和验证的线性，减小误差，而在WFV数据集建模分析中，建模效果最好的模型为输入风速因子的WFV-BP3，而验证效果最好的模型为WFV-BP1。在OLI数据集中，日平均风速平均值小于2 m/s，最大值小于4 m/s，而WFV数据集中，风速要明显大于OLI数据集，本研究认为风速因子引入建模数据集可以有效提高建模的线性，但也存在增加过拟合现象的风险。

WFV和OLI数据集建模效果对比分析表明，根据图3，虽然WFV数据集建模绝对误差最大值要远高于OLI数据集，但是其1/4和3/4位数与OLI数据集没有明显差距，验证效果对比分析表明，WFV数据集验证绝对误差明显低于OLI，其四分位数明显低于OLI。WFV-BP1模型总体误差最低。该模型验证数据集的绝对误差最小值虽然不是所有模型中最低的，但是其1/4位数和3/4位数都是所有模型中最低的，说明该模型是本研究所对比的模型中最稳定的模型。由于GF-1 WFV具有比Landsat 8 OLI更高的访问频率，更短的重访周期，同一时期所获得的数据更多，模型的训练效果更好，所以在波段数比OLI少的劣势下，WFV的建模效果却优于OLI。

Figure 4. Absolute error analysis map

图4. 绝对误差分析图

4. 结论

(1) 多元线性回归模型过拟合现象严重，建模线性很好，误差很小，但验证误差较大。在样本数较小的情况下BP神经网络可以有效抑制过拟合现象，建模和验证误差最低，误差最小的模型是WFV-BP1模型，该模型验证数据集RMSE和rRMSE分别为0.0246和40.74%，|ε验证|1/4位数和3/4位数分别为15.89%和48.88%。验证数据决定系数R²较高的模型，其误差反而更大，判断模型的好坏不能只看R²。

(2) 风速对水体反射率有一定的影响，将其作为参数输入模型，可有效提高建模的线性，但也存在增加过拟合现象的风险。

(3) OLI和WFV两种不同传感器对比验证的结果表明，由于GF-1 WFV具有比Landsat 8 OLI更高的访问频率，更短的重访周期，同一时期所获得的数据更多，模型的训练效果更好，所以在波段数比OLI少的劣势下，WFV的建模效果却优于OLI。

(4) 贫营养湖泊污染物遥感估算中，卫星获得的表面反射率信号受到大气的干扰以及水中的有色物质的干扰对估算结果的影响比富营养湖泊要大得多，对卫星传感器的灵敏程度也提出了更高的要求。本研究采用了70个不同日期的卫星影像开展研究，不同日期影像的大气干扰不同，增加了大气校正难度，也增加了表面反射率误差。

(5) 与浑浊水体的类似研究相比较，本研究的相对误差较高，但绝对误差却更低，原因是本研究针对本底值较低的贫营养湖泊，其对信号噪声敏感度较高，而较低的真实值决定了其绝对误差会远低于浑浊水体。

最后，本研究的经验将对贫营养湖泊非光学敏感物质参数的遥感反演有一定的参考和借鉴作用。

基金项目

云南省地方本科高校(部分)基础研究联合专项项目(202001BA070001-061，202101BA070001-081)，云南省大学生创新创业训练计划项目(202111390023)联合资助。

NOTES

^*通讯作者。