密云水库位于北京市密云区北部, 面积约180 km², 地处燕山群山丘陵之中, 是北京唯一饮水源供应地, 为华北地区最大水库, 也是亚洲最大人工湖, 被称为“燕山明珠”。 密云水库有两大入库河流, 分别是白河和潮河, 两支流入库总氮含量年平均值分别为2.05和5.09 mg·L^-1, 总磷年平均值分别为0.013和0.012 mg·L^-1[11]。 南水北调通过京密引水渠进入水库, 成为第三大入库水源。 2014年底南水进京, 密云水库的水源供给减小, 但涵养城市水源、 水源战略储备意义不容忽视。 近年来, 随北京市迅速城市化及北部经济发展, 水库及周边自然环境逐渐恶化, 尽管政府已开展一系列保护水库行动, 但水库污染情况依然存在。 因此, 开展对密云水库水质参数反演研究十分迫切。

1.2.1 无人机高光谱遥感数据

研究采用DJI M600Pro无人机搭载Headwall Nano-Hyperspec超微型机载高光谱成像光谱仪, 其光谱范围为400~1 000 nm, 共640个空间通道, 270个光谱通道, 光谱分辨率(FWHM)为6 nm。 GPS/IMU模块记录成像仪的姿态信息, 在Headwall提供的后处理软件中载入姿态信息后, 可对采集的条带进行几何校正、 辐射校正操作。 无人机采集镜头焦距为8 mm, 视场角为33°, 飞行高度为290 m。 选取密云水库两座主坝为研究区, 分别位于水库西部溪翁庄村的白河大坝及水库南部南碱厂村的潮河大坝。 数据采集时间为2018年10月、 2019年4月、 6月及8月, 如图1所示。 这两个区域分别为潮河、 白河入库口, 是密云水库水质变化敏感区。 首先对高光谱数据进行传感器辐射定标、 场地辐射校正、 波段降维、 去噪等处理。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 密云水库区位及采样点分布图Fig.1 Location and sampling point distribution of Miyun Reservoir

采用低空飞行方案, 大气环境单一且飞行时间较短, 因此不考虑大气影响。 依次进行传感器辐射定标、 场地辐射校正。 传感器辐射定标采用无人机设备供应方提供的配置校正文档, 在输出信号和辐射强度值之间提供模型和转换参数。 场地辐射校正采用实验室定标的参考板进行, 在提取无人机高光谱影像上定标板相应位置的辐射值后, 将定标板标准辐射反射率值(R)与高光谱遥感影像上的辐射亮度值(L)进行线性回归即可得出每个波段定标系数。 各波段遥感反射率与辐射亮度值之间的关系为

Ri=aiLi+bi(1)

式(1)中, a_i和b_i为各波段定标系数。

1.2.1.1 基于递归特征消除法的高光谱波段降维

对无人机高光谱数据进行辐射亮度校正及几何校正后, 进行波段降维处理。 因高光谱影像270个波段数据量较大, 数据冗余, 运算速度较慢, 且深度神经网络算法复杂度较高, 需要对270个波段进行降维处理。 递归特征消除法基于模型后向搜索, 在每次迭代删除特征的过程中, 将模型内部特征重要性评价指标作为度量标准, 评估当前特征子集与当前剩余特征子集之间的关系, 得到波段间的递归特征指数^[12]。 本研究将所有波段反射率递归特征指数与总氮、 总磷两个水质参数进行相关分析, 得到水质参数与各波段的相关系数热点图(图2)。 因密云水库水体较为清洁, 大于900 nm波段光谱反射率较低, 故不考虑大于900 nm波段。 结果表明, 高光谱遥感数据的550~650和750~880 nm处为对总氮浓度敏感波段, 730~900 nm为对总磷浓度敏感波段, 以此筛选出10个主要波段, 涵盖主要高光谱信息特征, 使其既能满足数据建模要求, 也能满足运行计算的硬件环境。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 水质参数与高光谱波段反射率相关系数热点图Fig.2 Correlation heatmap of hyperspectral band and water quality parameters

1.2.1.2 基于压缩感知方法的高光谱影像去噪

采用压缩感知去噪方法对无人机高光谱数据去噪。 根据压缩感知理论, 若信号是可压缩的或在某个变换域是稀疏的, 那么可以用一个与变换基不相关的观测矩阵将高维信号 X投影到低维空间Y, 而后通过重构算法利用最优化方法将原始信号X精确或者高概率重构^[13]。 获得信号测量值的数学模型如式(2)

y=ΦΨ=ΦΨs=Θs(2)

式(2)中, Φ ∈ R^M×N为测量矩阵, Ψ 为稀疏表示基, y为X的测量值, s为稀疏向量。

压缩感知的信号重构过程是一个求解线性方程过程。 式(1)未知数的个数远大于方程的个数, 是一个不定方程, 有无穷多解。 但由于向量s是稀疏的, 如果矩阵Θ 中任意2k列线性无关, 理论证明, 同时满足上面两个要求时, 可通过求解如下问题来精确重构信号

minθ‖θ‖0s.t. y=φX=Θθ(3)

式(3)中, ‖ θ‖ ₀表示向量的范数, 即向量非零元素个数。 通过求θ的唯一解来得到信号X。

经过去噪处理获得高光谱反射率结果(图3)显示, 压缩感知算法有效剔除了异常噪声, 其光谱曲线更平滑连续, 避免了模型训练中由噪声带来的误差干扰。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 高光谱影像去噪前后反射率曲线图Fig.3 Reflectance curves of hyperspectral images before and after denoising

1.2.2 地面监测数据

在密云水库水域共计采集有效水样80个, 采样时间与无人机遥感数据采集时间相同步, 对每个样点记录其位置信息, 四期样点的分布情况如图1所示, 采样深度为表层0~50 cm水面。 监测水质指标为总磷、 总氮, 其中, 总磷是指溶解的、 悬浮的、 有机的和无机磷。 总氮是指水中各种氮化合物的总和, 不包含氮气, 检测结果如表1所示。

表1

Table 1

表1(Table 1)

表1 地面采样水质参数检测结果 Table 1 Measured values of water quality parameter of ground sampling 水质参数 范围/(mg·L-1) 均值/(mg·L-1) 水质检测描述 总氮 2.624~4.023 2.842 9 在含Pt催化剂的TOC燃烧管中, TN由NO被臭氧(O3)激发成NO2, 返回基态时, 测量其发射光(590~2 500 nm) 总磷 0.038~0.133 0.071 9 依据GB/T1893—1989水质总磷的测定, 用钼酸铵分光光度测定总磷的方法

表1 地面采样水质参数检测结果 Table 1 Measured values of water quality parameter of ground sampling

1.2.3 模型构建数据集

经过以上预处理, 根据采样点的经纬度坐标, 将80个地面有效监测数据与无人机高光谱降维后获取的10个波段光谱反射率空间信息对应匹配, 得到80行10列模型构建数据集, 以高光谱反射率信息为输入数据, 总氮、 总磷两组实测水质参数为目标输出数据, 引入基于样本的10折交叉验证, 将训练数据划分为10个子集, 其中8个用来训练, 2个用来验证。 建立深度神经网络模型进行训练, 反演得到密云水库2个敏感区域空间分布特征。

1.3.1 深度置信神经网络

深度置信神经网络(deep belief network, DBN)是深度学习的一种重要算法, 由输入层、 隐含层和输出层组成。 2006年, Hinton提出了受限玻尔兹曼机(restricted Boltzmann machine, RBM)逐层贪婪预训练方法^[14]。 DBN网络结构为本研究基础模型, 其输入数据、 输出数据分别为网络输入层及输出层, 分别对应高光谱数据10个波段光谱反射率信息及水质参数值。 其训练过程可以分为自下而上无监督学习和自上而下有监督调优两部分。

如图4所示, 首先对一个只含有一层隐含层DBN进行预训练, 作为知识源域; 随后增加目标域, 即增加固定规模隐含层, 通过预训练重构误差, 确定网络停止增长的标准。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 DBN网络结构Fig.4 The structure of DBN network

在训练过程中, 每一个RBM隐含层对上一层输入的向量进行吉布斯采样, 得到隐含层神经元激活概率, 获取得到本层RBM隐含层神经元激活概率, 从而可得到全部RBM隐含层神经元激活概率。 上述吉布斯采样及参数更新方法可以表示为

P(hj=1v, θ)=σbj+∑iΔWijvi(4)

P(vi=1h, θ)=σaj+∑jΔWijhi(5)

Δθij=Δθij+αs-∂L∂θ(6)

式(4)—式(6)中, v为输入向量, θ为网络参数的集合, 其中W_ij为可视层与隐含层之间的权值矩阵, a_i为可视层与隐含层之间的权值矩阵, b_j为可视层神经元的偏置, h为隐含层n维状态变量。

1.3.2 基于随机子空间的深度置信神经网络算法

DBN算法虽然具有较高的计算性能, 但在处理高光谱等大规模数据时, 存在训练时间过长的问题, 且随着隐层数增加容易导致过拟合。 针对这些问题, 采用多个DBN学习器进行集成, 优化为基于随机子空间的集成深度置信神经网络算法(ensemble deep belief network, EDBN)。 该算法从原始特征空间随机选择子空间构建多个基学习器, 最后结合所有基学习器输出。 如图5所示, EDBN模型对数据样本进行划分, 即将特征空间T划分为K个子空间, 每个基学习器在样本子空间上独立并行训练, 从而形成K个DBN模型。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 EDBN模型结构Fig.5 The structure of EDBN model

对于随机子空间, 每个随机选择的特征向量都可以生成原始样本视图, 因此, 对于多个随机采样特征会生成多个视图, 即从不同视图对样本进行分析和描述。 根据原始样本在不同子空间投影视图, 设计多个结构相同的基学习器, 对每个基学习器进行独立并行训练。

1.3.3 EDBN水质反演模型的构建与参数优化

EDBN水质反演模型构建主要包括网络结构初始化及超参数优化调整两方面。 以无人机高光谱降维后10个波段为样本数据, 初始化网络结构为3隐层5节点, 使每一层提取到水质数据基本特征, 向前传播输送进网络下一层, 随后逐步增加隐层数和节点数, 观察并记录损失值变化。 损失函数选取平均值误差函数, 如式(7)所示。

MSE=1n∑i=1n(xi+yi)2(7)

在深度神经网络中, EDBN网络结果训练后会捕捉一组初始权重, 反向传播计算损失值来优化网络各层权重比。 通过调节隐层数和节点数使网络训练精度明显提高, 节点数梯度设置为5、 10、 15, 隐层数梯度为3、 4、 5、 6、 7, 其对应损失值如图6所示, 当损失值较小, 且继续增大隐层数和节点数精度提高不再明显时, 认为当前参数为最佳参数。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 不同隐层数和节点数对应的损失值Fig.6 Loss values corresponding to different numbers of hidden layers and nodes

样本训练批次是影响模型训练效率、 收敛速度及整体精度的重要参数; 选取3个经验值分批次讨论(图7), 较大训练批次训练精度较低, 中等训练批次效率较低, 综合考虑收敛速度、 损失值大小及训练时间, 对总氮、 总磷设置训练批次为9。 此外, 可以发现在迭代次数达到一定数量时, 损失值减小趋势趋于稳定, 依此可选取最大迭代次数。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 不同训练批次下的损失值收敛情况Fig.7 Convergence of loss values under different training batches

网络其他相关参数设置如下: 激活函数选用稀疏性较强的Relu函数, 以增强网络非线性变换特征。 初始学习率选取为0.001, 观察测试集误差变化, 当误差持续一段时间不变, 则随着迭代次数增加, 通过乘以数值为0.5的乘数因子渐渐衰减学习速率。 另外, 采用Adam优化函数防止模型训练中产生的过拟合现象。 得到EDBN水质反演模型参数表(表2)。

表2

Table 2

表2(Table 2)

表2 模型参数表 Table 2 Paramters of EDBN model 水质参数 隐层数 节点数 迭代次数 批处理量 TN 7 15 270 9 TP 6 15 160 9

表2 模型参数表 Table 2 Paramters of EDBN model

基于构建的深度神经网络模型, 对2018年10月、 2019年4月、 6月、 8月四期密云水库无人机高光谱数据进行反演, 得到总氮和总磷水质参数空间分布数据。 两种水质参数反演模型精度验证如图8所示, 可以直观看出样本真实值与预测值间的离散程度。 总体而言, EDBN模型对水质参数反演取得较好拟合效果。 其中, 总氮反演模型R²达到0.835 5, 异常点较少, 离散程度较小, 总磷反演模型R²达到0.770 3, 能满足较高精度的水质监测要求。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 EDBN水质参数反演模型精度验证Fig.8 Accuracy verification of EDBN model

EDBN水质参数模型误差值如表3所示, 总磷反演模型模拟误差相对较小, MSE、 MAE分别为0.000 8、 0.001 4, 其次为总氮反演模型, 其MSE、 MAE分别为0.015 3、 0.051 4; 综合来看, EDBN模型对水质参数反演效果良好, 适宜于密云水库总氮、 总磷水质参数反演。

表3

Table 3

表3(Table 3)

表3 EDBN水质参数反演模型误差表 Table 3 Errors of EBDN model 水质参数 MSE MAE R2 TN 0.015 3 0.051 4 0.835 5 TP 0.000 8 0.001 4 0.770 3

表3 EDBN水质参数反演模型误差表 Table 3 Errors of EBDN model

2.2.1 水质参数总氮时空演变特征

密云水库水质参数总氮浓度估算范围为2.71~2.9 mg·L^-1(图9), 相对偏高, 大部分时间为Ⅲ 类水。 其中, 潮河大坝区域10月浓度为2.778 mg·L^-1, 水质整体较差; 次年4月起总氮浓度稍有所降低, 平均值为2.764 mg·L^-1; 而到6月总氮浓度降为2.752 mg·L^-1。 随夏季到来, 在8月降到最小值, 为2.733 mg·L^-1。 从时间上来看, 潮河大坝总氮浓度季节差异显著, 具有较强的规律性, 两次较大波动出现在4月和8月, 水库水质先受到汛期降雨略有好转, 而后在水库自净调节作用下恢复至较低水平; 从空间上来看, 总氮浓度整体分布均匀, 春季北部略低于南部, 这可能是由于南部地区村庄聚集, 总氮浓度偏高与该区域的人类活动有关。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 潮河大坝和白河大坝水域总氮浓度空间分布图Fig.9 Spatial distributions of TN concentrations in the water areas of Chao river dam and Bai river dam

白河大坝区域10月总氮浓度为2.754 mg·L^-1, 处于常规水平; 而8月总氮浓度较低, 平均值为2.748 mg·L^-1; 而次年4月和6月总氮浓度较高, 分别为2.765和2.776 mg·L^-1。 从时间上来看, 夏季、 秋季总氮浓度较低, 春季浓度较高, 可能与夏季降水补充水源有关。 从空间分布来看, 总氮浓度有西北至东南波动趋向, 这可能是由于入库口水源注入带来一定程度的污染。

两区域对比可发现, 潮河大坝整体总氮浓度略高于白河大坝, 且潮河大坝季节差异更为明显, 产生这种现象可能与潮河大坝南部连片分布的村庄、 果园、 耕地有关, 长期施肥或粪便、 垃圾堆放引起土壤表层总氮变化, 造成总氮浓度显著差异。 从变化范围来看, 两区域范围内水质分布较为均匀, 白河大坝浓度分布波动区间更小, 没有较大梯度差异与总氮聚集现象, 这表明在实施严格水源保护政策下, 该区域农业生产和农村生活对总氮的影响有限, 工业污染等其他来源的氮污染治理难度较大, 需采取长期的保护措施, 逐步治理流域水质。

2.2.2 水质参数总磷时空演变特征

密云水库总磷浓度估算范围为0.007 6~0.012 mg·L^-1(图10), 相对浓度偏低, 大部分时间为Ⅱ类水, 夏季可达到Ⅰ类水标准。 其中, 白河大坝区域10月总磷浓度为0.008 9 mg·L^-1, 随后逐渐升高, 次年4月平均值为0.009 1 mg·L^-1, 而到6月份总磷浓度达到最高, 为0.009 5 mg·L^-1, 到8月回落至最低值0.008 1 mg·L^-1, 水质达到最佳水平。 从空间分布来看, 密云水库总磷浓度与总氮浓度相似, 略有西北至东南波动的趋向。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 潮河大坝和白河大坝水域总磷浓度空间分布图Fig.10 Spatial distributions of TP concentrations in the water areas of Chao river dam and Bai river dam

潮河大坝水域总磷浓度10月浓度均值为0.009 2 mg·L^-1, 次年4月稳定在为0.009 3 mg·L^-1, 而到6月总磷浓度略有升高, 均值为0.009 5 mg·L^-1, 到8月回落至最低值0.008 8 mg·L^-1。 从时间上来看, 潮河大坝区域总磷浓度随季节变化不明显, 夏季总磷浓度相对较低; 从空间上来看, 潮河大坝北部出现波动, 略高于南部, 总体浓度相近, 分布均匀。

总体来看, 密云水库总磷浓度分布较为稳定, 两个区域均在6月达到最高, 8月回落至最低, 夏季整体水质较好, 其污染物迁移规律较为相似。 主要考虑随季节变化, 夏季降雨增多, 南水北调来水入库量增加, 使库区水源得到了有效补充, 降低了水体总磷浓度。 从变化趋势来看, 库区总磷浓度变化幅度较小。 随水库上游水资源补给增多, 白河、 潮河两个主坝附近水域水动力条件愈发强烈, 考虑白河大坝流量高于潮河, 使得白河大坝附近总磷浓度总体略低于潮河大坝区域。

研究结果表明, 集成深度置信神经网络算法模型可以实现密云水库总氮、 总磷浓度高精度反演。 与人工神经网络、 支持向量机等其他水质反演模型相比, 本模型误差更小^[10]。 集成深度置信网络综合所有基学习器训练特征后输出, 有效避免了同类模型过拟合问题, 使得检索模型更加准确, 结果具有更好的泛化能力^[15]。 同时, 基于无人机高光谱数据的反演模型实现了密云水库小尺度水质监测, 结果显示, 与采用Landsat8、 GF-1等卫星遥感数据反演密云水库水质参数结果基本一致, 是对基于卫星遥感数据宏观水质监测的有效补充^{[9, 17]}。

分析结果表明, 密云水库潮河大坝和白河大坝水域总氮、 总磷浓度在时空分布上均有一定波动规律, 同一水质参数对应的不同区域也存在一定差异。 总氮浓度变化波动相对较大, 在秋季潮河大坝出现污染聚集, 主要考虑库区缓冲区内存在连片耕地和建设用地, 农业面源污染和生活污水在一定程度上加剧了总氮污染。 而夏季总氮浓度降低, 这主要是由于夏季降水充沛, 南水北调入库水量增加, 水源相对充足, 水库水体得以更新补给, 使污染物聚集难度大。 总磷浓度波动较为平稳, 磷污染风险较小, 这表明密云水库周边区域对磷污染控制效果良好。

密云水库水质参数浓度反演结果显示, 水体总氮浓度偏高, 主要为Ⅲ 类水。 而总磷浓度总体处于Ⅱ类水, 在夏季部分区域可达Ⅰ类水标准。 这与张雨航、 秦丽欢等人的研究结论具有相似性^{[18, 19]}。 因此, 密云水库总体水质较好, 可以满足饮用水源的标准, 但仍存在富营养化风险, 仍需加强对氮、 磷实时监测, 加强水质防控工作。 对于总氮和总磷等污染物, 建议整体优化密云水库周边区域山水林田湖草空间格局, 增加林草等生态用地对污染物层层过滤拦截作用。 此外, 还需严格监管密云水库上游地区人类季节性活动, 减少上游面源污染、 畜禽养殖等污染物入河量。

基于无人机高光谱数据, 采用深度神经网络模型, 获得密云水库敏感区域高精度水质参数反演结果, 但也存在一些不足之处。 对于EDBN网络中隐层数和节点数调节步幅基于经验设定初始值, 对反演模型构建带来一定影响。 随深度神经网络算法研究不断进步, 训练过程将更加高效、 简洁。