随着浙江省“五水共治”以及“剿灭劣Ⅴ类水”战略的深入实施, 养殖尾水乱排直排现象面临巨大的环境压力, 养殖尾水生态化处理迫在眉睫。截止2020年, 浙江省主要通过“沉淀池+过滤坝1+曝气池+过滤坝2+生态池(简称“三池两坝”)”技术完成了内陆约5.33万hm²养殖池塘尾水治理任务(刘梅等, 2021)。但是由于养殖池塘及尾水处理池小而散、涉及面积广等问题, 且养殖池塘对水质的监测要求具有实时、动态、快速的特点, 仍缺乏有效的水质监测体系, 阻碍了现代养殖小区数字化养殖场的创建与发展(李道亮等, 2020)。

目前, 对养殖池塘或者养殖尾水处理效果采取的仍是常规的水质监测方法, 即对特定养殖池塘或尾水处理系统定期、定点进行长年累月采样及水质监测, 这种方法受人力、物力、时间及天气的限制, 采集的数据量不可能太多, 而且成本高、速度慢; 并且对于整个养殖小区而言, 这些采样点数据只能代表局部水质状况, 难以获得大范围淡水养殖水域水质参数的分布和变化情况, 不能满足对养殖池塘水质大尺度、快速、实时的监测要求。因此, 迫切需求一种实时、快速监测养殖池塘及养殖尾水处理池水质动态变化的有效手段。

水体的光谱特征主要是由水体中的各种物质对太阳辐射吸收和散射共同作用形成的(Koponen et al, 2002), 不同的水体其光谱特征存在一定差异(黄宇业等, 2019)。国内外已有学者尝试将无人机多光谱技术应用小微水域水质要素的快速反演(Bansod et al, 2018; 杨振等, 2020), 如刘彦君等(2019)利用无人机多光谱影像技术建立了小型湖泊悬浮物浓度(SS)、浊度(TUB)与高光谱参数的反演模型, 徐福祥等(2018)利用无人机实现了对海滨景区绿潮的精准监测, 但是由于其仅有4~6个波段, 信息量少, 不能全面反演水质参数变化, 尤其是非光敏化学参数(李俊生, 2007; 张兵等, 2009), 如TN、TP、NH₄⁺-N、COD_Mn等指标, 预测精度达不到养殖要求。而高光谱摄像头, 因其包含连续的窄波段(10~20 nm), 具有最丰富的空间和光谱信息, 其精度也最好。Bansod等(2018)对恒河的水质参数进行了基于高光谱图像的水质反演评价; 杨振等(2020)利用机载高光谱技术建立了矿区水库高光谱影像和悬浮物、浊度的反演模型; 陈俊英等(2019)建立了基于地物高光谱污水处理厂水质综合评价方法; 黄宇等(2020)利用无人机高光谱成像技术对河湖水质TSS、浊度、TN、TP和叶绿素a进行了反演, 取得了较好预测效果。目前, 对我国内陆淡水养殖池塘及养殖尾水处理效果的高光谱水质监测技术还未见相关报道。

鉴于此, 本文以淡水养殖池塘及配套尾水处理池为研究养殖小区, 通过无人机搭载高光谱传感器获取光谱反射率数据, 建立水质指标的反演模型并研究其浓度空间分布, 为养殖池塘水环境实时调控、养殖尾水池处理效果评价及构建数字化养殖小区提供技术支撑。

1 材料与方法 1.1 研究区概况

本研究选择湖州市典型主要养殖品种之一的青虾养殖小区尾水处理区及周边池塘为研究区域, 位于湖州市菱湖镇杨港现代生态渔业园内(120°7′8.7″N, 30°46′32.7″E)。2018年该养殖小区采用“三池两坝”技术建立完成尾水处理区, 截止2019年湖州近4万hm²的养殖池塘均采用该技术实现尾水的生态化处理, 减轻了对周边河网水体及太湖的富营养化影响。但是由于在后期的出水水质监测均采用传统方法, 费时费力, 也不利于养殖小区数字化建设及管控。

1.2 采样点布设和水质检测

于2020年12月在研究区域开展野外水质取样工作, 采样点位置如图 1所示, 在研究过程中发现, 由于养殖池塘和尾水处理池面积一般小于1 hm², 水深在2 m范围内, 因此, 水质指标空间差异并不会太明显。综上, 每个养殖池塘仅在岸边进行采样, 采样点总共为49个, 每个点采集表层0.2 m处水样500 mL进行实验室水质分析, 分析参数包括TSS、TN、TP、COD_Mn和NH₄⁺-N, 所有水质指标的检测均参照《水和废水监测分析方法(第四版)》进行, 其中TN采用碱性过硫酸钾分光光度法(GB/T 11894), TP采用钼酸铵分光光度法(GB/T 11893), NH₄⁺-N采用纳氏试剂法(HJ 535—2009), COD_Mn采用酸性高锰酸盐指数法(GB/T 11892), TSS采用重量法(GB/T 11901)。

1.3 高光谱影像的获取与处理

同时利用大疆无人机M600 Pro搭载高光谱成像仪GaiaSky-mini获取尾水处理区及周边养殖池塘的高光谱影像。无人机飞行高度为100 m, 无人机高光谱相机的扫描范围为400~1 000 nm, 光谱通道为176, 光谱分辨率为3.5 nm。所选实验时间均晴朗微风(1~2级), 无云层覆盖。

无人机高光谱影像的预处理主要包括镜像变换、黑白帧校正、场地校正等, 均在高光谱系统自带软件完成。然后对拍摄的大量数据进行研究区图像拼接, 将拼接后的影像导入ENVI5.3软件, 根据经纬度坐标, 定位49个水面采样点的具体位置, 以采样点为中心的6×6像素密度作为计算区域, 再将该区域内光谱反射率求平均值, 该平均值即可作为该采样点的反射率数据, 共获得49组与地面水质数据相对应的光谱反射率, 随机选取39组数据用于模型训练样本, 其余10组用于模型检验样本。

1.4 高光谱指数分析

将水样中各水质指标与采样点在遥感影像中的反射率进行相关性分析, 获得研究区不同水质参数的敏感单波段。诸多学者研究结果显示, 将反射率进行适当的数学变换, 构建光谱任意两波段组合而成的差值光谱指数(difference spectral indices, DSI)、归一化光谱指数(normalized spectral indices, NSI)和比值光谱指数(ratio spectral indices, RSI), 可以降低背景信息的干扰, 将更有效提取光谱信息(黄宇等, 2020; 黄昕晰, 2020), 计算如式(1)、式(2)和式(3)所示。

式中, R_λ1和R_λ2分别代表400~1 000 nm波段范围内任意2个波段的反射率。

1.5 水质反演模型构建流程

近年来, 高光谱技术凭借丰富的信息量和多元的处理方法得到飞速发展(陈俊英等, 2019; 杨振等, 2020)。然而, 高光谱影像在提供丰富波谱信息的同时存在数据量大、数据冗余度高的缺点, 难免会出现多重共线性问题(黄昕晰等, 2020)。因此, 高光谱数据处理的核心步骤是选择适当的数学建模方法建立光谱反射率与样品组分之间的相关性(李俊生, 2007)。本研究将获取的光谱反射率与实测水质指标(TSS, TN, TP, COD_Mn, NH₄⁺-N)进行Pearson相关性分析, 然后选择相关性达标的敏感波段分别采用线性函数、指数函数和多项式函数进行拟合, 筛选各水质要素最佳反演模型。对于近些年发展的机器学习算法(如偏最小二乘法、人工神经网络、支持向量机)则不需要筛选相关系数最高的光谱反射率, 可一定程度上提高监测精度, 但运用的波段数多, 反演运行时间相对较长, 而养殖小区一般在100 hm²范围内, 反演运行时间基本在1 min内(运行环境为Widowns32位, 酷睿i5, 独显2G), 运行时间上可接受的。因此, 为了提高预测精度, 本研究选择全波段进行机器学习算法建模分析。综上, 本文选择敏感波段函数法和机器学习算法进行建模, 比较反演模型的优缺点。其中3种机器学习算法建模原理如下:

偏最小二乘法(partial least square, PLS)通过将自变量和因变量的高维数据空间投影到相应的低维空间, 得到相互正交的特征向量, 再建立线性回归模型, 不仅克服多重共线性问题, 即在回归的同时强调自变量对因变量的解释和预测作用, 而且能较好地辨识系统信息和噪声。

径向基神经网络(radial basis function network, RBF)属于多层前馈神经网络, 具有自组织、自学习、知识推理的特点, 对大量非结构性、非精确性规律具有自适应、信息记忆和非线性逼近能力以及优化计算等优点。RBF在Matlab中实现, 功能函数采用Sigmoid函数。

支持向量机(support vector machine, SVM)是以统计学习理论为基础、训练误差为约束条件、以结构风险最小化为优化目标的一种专门用来研究有限样本预测的高效机器学习方法。具有较强的泛化能力、较好的鲁棒性以及没有局部最小点的优势。本研究SVM采用R语言中的e1017包构建预测模型, 使用radial (径向基函数)为核函数, 成本函数cost和gamma选择模型中最佳参数。

最后根据模型的相关评价标准筛选各水质要素最佳反演模型, 再将最优的监测模型反演到无人机高光谱影像上, 分别制作各水质参数的空间分布图, 具体计算路线如图 2所示。

1.6 模型筛选评价指标

水质参数预测模型采用决定系数(coefficient of determination, R²)、均方根误差(root mean square error, RMSE)、相对分析误差(relative prediction deviation, RPD)以及平均绝对误差(relative absolute error, MAE)对比模型稳定性及精度。其中R²越大, 说明模型越稳定, 一般大于0.8即说明模型稳定。RMSE越小, 则模型精度越高。RPD > 2.0时表示模型可靠性高, 具有极好的预测能力, RPD值介于1.5~2.0之间, 表明模型可对样品做粗略估测, 模型有待改进, RPD < 1.5时, 表明模型预测能力差(叶勤等, 2017)。MAE < 20%可用于养殖水质的预测, 满足养殖水体水质反演误差要求。其计算公式如下:

式中, y_i表示实测值, 表示预测值, y表示实测值的平均值, n代表样本数量, S_D表示样本观测值方差。

2 结果与分析 2.1 养殖池塘和尾水处理区光谱特征

研究养殖小区共49个采样点的光谱反射率曲线如图 3所示, 其中图 3a列出了40个养殖池塘采样点的反射率曲线, 图 3b列出了9个尾水区采样点的反射率曲线。可以看出池塘水体的光谱特征在400~ 560 nm范围内, 光谱反射率呈上升趋势, 至560~ 580 nm附近形成一个反射高峰, 这主要是是由水体中叶绿素和胡萝卜素在此波段具有较弱的吸收特点, 同时加上藻类和悬浮物的散射作用而引起反射高峰(Gitelson et al, 1993)。而到580 nm后, 反射率曲线呈下降趋势, 直至680~710 nm处附近则形成一个峰谷, 主要是由于水体中的叶绿素a在此波段具有强吸收的特点(唐军武等, 2004); 而紧接着在790~810 nm范围内形成一个峰值, 主要是由于水中悬浮物的散射作用(黄宇业等, 2019)。而尾水处理区由于经过过滤、沉淀、曝气等净化过程, 水体中悬浮物和藻类去除率可达到70%以上, 因此, 尾水处理区在680~710 nm和790~810 nm范围内均未形成明星的峰值(图 3b), 而其他波段的变化规律和养殖池塘具有一致性。综合分析光谱曲线可知, 不同池塘水体的光谱曲线变化趋势总体一致, 但是由于不同采样点具有不同的水质指标浓度, 导致其峰谷值及曲线高低变换存在一定差异。

2.2 水质参数与单波段反射率的相关性

将养殖池塘和尾水处理区采样点的水质指标(TSS, TN, TP, COD_Mn, NH₄⁺-N)分别与其对应的光谱反射率值进行相关性分析, 得到如图 4所示的400~1 000 nm处的相关系数曲线。

从图 4可知, TSS和NH₄⁺-N与各波段反射率相关性变化趋势较一致, 在700~820 nm范围内保持较高的相关性, 在730 nm附近出现峰值, 其最大相关系数分别为0.69和0.75; TN和TP与各波段反射率相关性变化趋势较一致, 其相关系数不高, 同样是在730 nm处出现峰值, 其最大相关系数分别为0.33和0.35; COD_Mn与各波段的反射率呈正负相关性, 在400~920 nm范围内呈正相关, 相关系数在690 nm附近出现峰值, 其最大相关系数为0.48, 在920~1 000 nm范围内呈负相关, 相关系数绝对值范围为0~0.23。

2.3 水质参数与双波段反射率的相关性

研究表明, 利用双波段组合因子可以更有效突出光谱特征, 其中相除因子和相差因子都是突出水质参数光谱特征波段的有效运算方法(吕航等, 2013; Yao et al, 2015)。本研究利用归一化指数、比值指数和差值指数处理400~1 000 nm的所有波段反射率, 再将这些指数分别与各水质参数进行相关分析, 得到水质参数与各波段比值的相关系数分布图(图 5), 以获取最佳的双波段组合, 用以构建水质反演模型。表 1列出了各水质参数与单波段及双波段反射率指数相关系数, 可以看出, 不同水质参数与波段数及处理方式存在较大差异, 其中除NH₄⁺-N外, TSS、TN、TP和COD_Mn和双波段相关系数均大于单波段, 其中TSS和TN与差值指数(波段168-111, 波段85-41)相关系数最大, TP与比值指数(波段65/68)相关系数最大, COD_Mn则与差值归一化指数{波段[(75-92)/ (75+92)]}相关系数最大, 而NH₄⁺-N则和单波段(729)相关系数最大。图 5为TSS、TN、TP和COD_Mn与双波段反射率指数相关系数分布图, 可以直观看出其对应的最佳波段。总体来说, 合适的预处理和数学变换对模型精度的提高有重要作用, 不同的处理方法效果不同, 应对原始数据进行不同的数学变换尽可能寻求与各水质参数相关性最好的处理方法。

2.4 最佳水质反演模型筛选与检验

通过筛选出每一种水质指标与波段反射率指数相关系数最高的单波段, 建立了各水质参数与波段因子之间的线性函数、指数函数和多项式函数, 同时通过模型训练建立了各水质参数所有波段的机器学习模型。表 2分别列出了TSS、TN、TP、COD_Mn和NH₄⁺-N各反演模型的精度评价, 从表 2可知, 波段差值指数的多项式函数用来反演TSS效果最好, R²为0.94, MAE为14.3%; 全波段RBF神经网络用来反演TN和COD_Mn效果最好, R²分别为0.89和0.85, MAE分别为15.8%和17.6%, 且偏最小二乘法、RBF神经网络和支持向量机的R²和MAE显著大于其他固定函数法; 对于TP, 则是全波段的支持向量机反演效果最好, R²为0.95, MAE为14.8%; 全波段的偏最小二乘法对NH₄⁺-N的反演效果最好, 其R²为0.88, 但是其MAE为22.1%, 大于20%, 不满足养殖水质预测精度要求, 后期需继续优化相关反演模型。

图 6列出了TSS、TN、TP、COD_Mn和NH₄⁺-N各自最佳反演模型的拟合效果。从图 6可知, TSS、TN、TP、COD_Mn和NH₄⁺-N均靠近1︰1线, 这与其决定系数R²在0.85~0.95之间相一致, 即筛选的各自模型均具有较优的稳定性, 但是各自模型在水质指数较高时, 离散度相对较高, 尤其是COD_Mn和NH₄⁺-N, 其MAE分别为17.6%和22.1%, 即模型鲁棒性相对TSS、TN和TP反演模型稍差些。TSS、TN、TP和COD_Mn其RPD均大于1.5, 说明模型达到预测精度要求, 而NH₄⁺-N其RPD为1.44, 表明模型还有待改进, 可以用来近似定量预测。综上, TSS、TN、TP、COD_Mn反演模型较可靠, 而NH₄⁺-N则须对模型进行改进以提高预测精度。

2.5 水质参数的高光谱图像反演

将无人机获取的同时期研究区的高光谱数据, 并结合构建的水质参数最佳预测模型进行反演, 计算出高光谱图像每一个像元对应的水质参数, 绘制出TSS、TN、TP、COD_Mn和NH₄⁺-N浓度空间分布图, 结果如图 7所示。从图 7可以看出, 各水质参数反演总体情况符合实际情况, 即尾水处理区的水质从沉淀池到生态池有一个浓度梯度逐渐降低的趋势, 至出水口处基本达到《淡水池塘养殖尾水排放标准SC/T 9101—2007》。尾水池北面是青虾养殖塘, 由于是12月下旬, 正是青虾捕捞上市季节, 大部分养殖池塘每天早上采用全塘地笼下网的方式进行“捕大留小”售卖, 由于人为的扰动, 引起池塘底泥上浮, 造成TSS、TN和TP含量偏高; 而尾水池南面的养殖池塘均是养殖大口黑鲈, 由于冬季大口黑鲈不摄食, 栖息在池塘底部, 故TSS、TN和TP含量相对较低; 青虾是种草养殖, 冬季时气温较低, 部分池塘水生植物未及时捞出而死亡腐烂, 释放有机细微碎屑进入养殖水体, 造成部分虾塘COD_Mn比鱼塘高, 且当天拍摄时为北风1~2级, 有机物质集中在南面下风口处, 故池塘南面COD_Mn高于北面。氨氮浓度过高, 容易引起养殖对象急性应激中毒, 是水产养殖上重点关注的水质指标之一(张武肖等, 2015)。从NH₄⁺-N反演图中可以看出个别池塘红色区域较大, 经后期实地调研及水质检测, 这些池塘NH₄⁺-N浓度确实较高, 与反演结果相一致, 并及时提醒养殖户进行水质调控, 防止鱼虾出现应激胁迫而爆发疾病。结果说明尽管NH₄⁺-N反演模型误差未达到养殖水质预测要求, 但可以近似定量反演NH₄⁺-N浓度及其趋势变化, 仍可以用来快速进行水质预警调控。

3 讨论

目前, 利用无人机遥感对小微水域的水质要素反演已有一些探索性应用(刘彦君等, 2019), 且绝大部分是采用多光谱传感器, 由于多光谱传感器波段通常只有5个波段, 只能对光敏参数(TSS, 浊度, 叶绿素等)进行很好的预测, 而对非光敏参数(TN, TP, COD_Mn, NH₄⁺-N等)的反演则存在很大限制(胡震天等, 2020)。本研究在无人机多光谱遥感的基础上, 选择高光谱传感器, 采用无人机在水面上空获取图像数据, 用于养殖小区水环境参数(TSS, TN, TP, COD_Mn, NH₄⁺-N)的预测, 取得了较好的反演效果, 为养殖水环境的实时快速监测及预警调控提供了实际可行的技术途径。

无人机高光谱技术在定量反演水质指标时, 由于具有较高的光谱分辨率和较强的波段连续性等特点, 可以获得更为全面的光谱特征信息, 但同时由于光谱测量中的某些人为和自然因素的干扰, 光谱数据需要进行不同的数学变换以减少系统误差和背景噪声的影响, 同时降低信息重叠度和冗余度, 进一步增强信噪比, 提高光谱数据与水质参数的相关性, 有利于提高线性模型的预测精度(赵庆展等, 2016)。本研究中3种数值变换方法明显提高了TSS和COD_Mn与光谱数据的相关性, 均达到0.81以上, 差值指数对TN与光谱数据的相关性提高最大, 比值指数则对TP与光谱数据的相关性提高最大, 而氨氮是例外, 数值变换方法甚至降低了其与光谱数据的相关性, 这可能是因为在增强敏感波段信息的同时也增强了噪声信息, 给整体数据带来了较大的反演误差, 降低了模型预测精度(谭玉敏等, 2014)。结果说明不同波段数值变换方法的处理效果不同, 应合理选择处理方法以进行后续建模, 这与相关研究结论也较为一致(董广香, 2007; 马娜等, 2010)。

另外, 研究发现选择与水质要素相关性较高的光谱波段甚至全光谱波段反射率进行建模, 再根据评价指标最优的反演模型用于区域影像的反演估测时, 依旧无法确定该模型是否效果最佳(史舟等, 2015)。如本研究中的NH₄⁺-N就出现了这一情况, 因为水质参数的差异体现在反射率光谱从可见光到近红外的多个波段, 其他非敏感波段构建的光谱参数对某些水质要素的影响也不可忽视(段洪涛等, 2006)。因此, 目前对于不同水质指标浓度的遥感反演, 并没有统一的标准模型, 不同研究区域的光谱反射率信息会因为其不同的地理位置、不同的高光谱传感器、不同的无人机飞行速度以及不同的水体类型等各种因素的变化而不同(Romero-Trigueros et al, 2017)。本研究中由NH₄⁺-N和光谱数据建立的反演模型其决定系数R²较高, 达到0.81, 但其预测误差MAE达不到20%的精度要求, 有待进一步探索合适的波段变换方法及数学模型。另外, 由于池塘水体光学特征的复杂性、水质影响因子的多样性(Guimaraes et al, 2019; 黄昕晰等, 2020), 如何更好地将水质信息与高光谱传感器相结合, 以达到通过无人机高光谱技术实时全面准确地反映养殖水体的水环境状况以及尾水处理池的处理效果, 值得进一步深入研究和探索。

4 结论

本文利用无人机高光谱遥感建立光谱反射率和采样点水质参数实测数据的反演模型, 可以快速实现养殖池塘水环境的精准预测及实时调控, 为养殖水环境水质要素的预测提供了一种高效的数据获取手段。结果表明, 将反射率进行适当的数学变换后可有效提高其与某些水质参数的相关性, 且波段差值指数的多项式函数更适合反演TSS, 全波段RBF神经网络用来反演TN和COD_Mn效果最好, 对于全波段的支持向量机用来反演TP效果最好, 全波段的偏最小二乘法对NH₄⁺-N的反演效果较好。尽管某一地区筛选的水质反演模型并不一定适合其他养殖小区, 但是其水质反演技术流程具有一致性, 为该技术成熟应用于其他养殖小区提供了技术路线。