基于无人机高光谱技术的水质预测反演系统设计与实现

马启良，原居林，张爱华，郦 晶，刘 梅

(1.湖州师范学院 信息技术中心，浙江 湖州 313000;2.浙江省淡水水产研究所，浙江 湖州 313000)

水是生命之源，是人类生存、生产、生活的基础.充足优质的水资源是生态系统得以健康循环的首要条件.然而，水资源极易受到污染，特别是一些内陆水体，由于自然封闭其污染问题更加突出.同时，水资源作为污染物的载体，具有动态的扩散和蔓延特性，从而进一步加剧了水体的污染程度.随着经济的高速发展和人类活动的增强，河道、池塘、河流、湖泊等水质污染问题日益严重，已成为制约人类可持续发展的关键因素.因此，有必要利用高新技术手段对水质污染问题展开研究.

随着无人机技术和高光谱遥感技术的日益成熟，无人机高光谱技术在水质预测方面得到了长足的发展，科研工作者也在不断地探索和研究光谱信息与水质指标之间的关系.当前水质监测已由定性描述转向定量分析，光谱技术可监测很多水质参数，其反演精度不断提高，在水资源的保护、规划和可持续发展方面发挥着重大作用.目前，无人机搭载成像高光谱仪已被广泛用于河道、池塘、水库等小面积水体的悬浮物浓度、浑度、透明度、总磷浓度、总氮浓度、化学需氧量、高锰酸钾指数、溶解氧、叶绿素a浓度、pH值、重金属含量等水质指标的监测.如:韩亚超等基于CASI 1500H航空高光谱数据对海湾悬浮物的监测进行研究，反演总体精度优于70%[1]；Bansod等利用机载可见光/红外成像光谱仪(AVIRIS-NG)对恒河的水质参数叶绿素a、浊度和总磷进行反演评价，结果表明，地面实况数据与光谱比率指数之间的Pearson相关性强于单一光谱带，且相关性大于0.6[2]；杨振等利用机载高光谱技术建立矿区水库高光谱影像和悬浮物、浊度的反演模型，并建立了悬浮物和浊度的空间分布图[3]；黄宇等利用无人机高光谱成像技术对河湖水质悬浮物、浊度、总氮、总磷和叶绿素a进行反演，且相关性都在0.6以上[4].

本文在前人研究成果的基础上，根据实际应用需求，对多种水质指标、光谱变换方法、回归分析模型进行集成，并以可视化的软件形式呈现给使用者，实现对不同水质指标的预测反演.

1 需求分析与系统设计

1.1 需求分析

目前，常规的水质监测方法主要分为两种:一种是通过水质检测传感器获取水质指标数据，这种方法虽然能够获取实时的水质指标数据，但获取的数据只具有局部代表性；另一种是实地采样检测，实验人员需在河道、池塘、水库、河流、湖泊等定点、定剖面进行长年累月地采样及水质检测，这种方法受人力、物力、时间和天气的限制，采集的数据量不可能太多，不仅成本高、速度慢，而且对整个水体而言，其采样点数据只具有局部和典型的代表意义，不能满足对水体水质大尺度、实时的监测需求.

随着无人机高光谱遥感技术的快速发展，研究人员不断地尝试利用该新技术，开拓水质检测的新思路、新方法，充分利用无人机高光谱遥感技术监测范围广、速度快、精度高等特点来满足实时、动态、大尺度的水质监测和评价要求.但面对拍摄的高光谱数据，如何有效、快速地分析出光谱与水质指标之间的关系是迫切需要解决的问题.目前的方法是借助高光谱设备公司的技术人员，对光谱数据进行处理，或借助国外的相关遥感软件(如ENVI等)，协助研究人员进行分析处理.但这些处理数据的方法费时费力，不利于实际应用，因此迫切需要将研究人员的一些好的、快速的方法进行集成，设计一整套高光谱图像与水质指标关系分析的操作流程.

1.2 系统设计

本文根据研究和应用需求，设计一套集模型训练与预测反演于一体的智能化水质预测反演系统，其主要框架结构如图1所示.该系统集成8个常用水质指标、9个光谱数据变换方法、5个光谱数据与水质指标的回归分析模型，具有光谱数据采集和水质预测反演等功能.

图1 水质预测反演系统框架Fig.1 Framework of the water quality prediction and inversion system

2 关键技术

2.1 光谱反射率数据提取

为获取水体采样点处的光谱反射率值，系统提供了标注采样点位置和导出光谱数据的功能.选择“开始标记采样点”后，点击左键可标注采样点，而点击右键则可取消标注的前一采样点；当采样点标注完成后，可以导出所有采样点处的光谱反射率值.假设所有采样点光谱反射率值为集合D={d1,d2,…,di}，其中，i为采样点数，di为第i个采样点的光谱反射率值，di={band1,band2,…,bandk}，k为高光谱图像的光谱波段数，bandk为第k个波段的光谱值，且光谱值是以标记点为中心的6X6区域的平均值.

2.2 数据变换方法

根据前人的研究经验[5-6]，本文对采样点处的光谱反射率数据与水体样本的实测数据进行相关性分析.选择最相关的一个光谱波段或多个光谱波段的组合，对常用的光谱反射率数据变换分析方法进行集成，其主要包括:对原光谱反射率值与水体指标实测值进行相关性分析；对光谱反射率值进行对数、倒数、对数的倒数等数学变换，与水体指标实测值进行相关性分析；采用双波段比值、双波段差值、双波段差值归一化等双波段组合形式，与水体指标进行相关性分析；采用PCA变换方法对高光谱数据进行降维处理，选择一个关键主成分或多个关键主成分，与水体指标进行相关性分析.

2.3 数据集划分

本文选取50组水质指标实验数据，其中40组作为校正集，10组作为验证集.为增加数据集划分的合理性，本文采用化学计量学中的经典算法SPXY(Sample set Partitioning based on jointX-Ydistance)[7]实现数据集的划分.SPXY方法是KS算法(Kennard Stone)[8]的延伸，其在KS算法的基础上，既考虑光谱空间对校正集选取的影响，又考虑浓度矩阵对校正集选取的影响，并基于X-Y联合距离对样本集进行选取.在SPXY方法中，每对样本p和q的X向量的欧氏距离为:

(1)

其中，Xp(j)和Xq(j)分别为光谱数据样本p和q在第j点的波长或波段数，J为波长数目或光谱的波段数.每对样本Y向量的距离为:

(2)

其中，yp和yq分别为样本p和q的标准分数.假设X空间与Y空间的权重相同，则X-Y间的距离为:

(3)

其中，max(dX(p,q))和max(dY(p,q))分别为X和Y向量距离的最大值.SPXY方法根据样本的X-Y联合距离dX,Y(p,q)，按顺序连续选出代表性样本分配给校正集，并重复这样的过程直至选出指定样本数目.

2.4 回归分析模型

系统在回归分析模型上，主要集成了偏最小二乘法(Partial Least Squares，PLS)线性回归模型、多项式回归模型(Polynomial Curve Fitting，PCF)、径向基(Radial Basis Function，RBF)神经网络模型、多层感知器(Multilayer Perceptron，MLP)神经网络模型，以及支持向量机(Support Vector Machine,SVM)回归模型，以满足研究人员对不同回归模型的需求.

2.5 水体提取与水质预测

2.5.1 归一化差异水体指数

1996年，Mcfeeters利用水体在绿光波段表现为强反射峰，在红光和近红外波段表现为吸收谷的光谱特性，提出了归一化差异水体指数(Normalized Difference Water Index, NDWI)[9].其计算方法为:

NDWI=(bandG-bandNIR)/(bandG+bandNIR)

(4)

其中，bandG为绿光波段光谱值，bandNIR为近红外波段光谱值.根据绿光和近红外光的波长范围及其所对应的波段，选择编号为57的波段作为绿光波段，编号为141的波段作为近红外波段，处理后的图像水体区域与非水体区域的对比度得到明显增强，如图2(b)所示.图2(c)为分类后的二值图像.

图2 FCM方法水体区域提取结果Fig.2 Water area extraction results of FCM method

2.5.2 高光谱图像分类

该系统选择模糊C-Means聚类(FCM)算法对以上方法处理的结果图像进行聚类，并将图像中每个像素作为一个样本点进行聚类分析，在确定分割阈值后获得水体区域二值图像.如图2(c)所示，白色为水体区域，黑色为非水体区域.

FCM算法最早是由Dunn[10]根据RusPini定义的模糊划分概率提出的，Bezdek对Dunn的目标函数进行了改进[11].该算法通过优化目标函数得到了每个样本点对所有类中心的隶属度，从而根据样本点的类属进行自动分类.目标函数定义为:

(5)

其中，U为隶属度矩阵；V为聚类中心集；D为数据集中数据总个数；N为设置的聚类中心数(2≤N≤D)；m为加权指数，一般取值为2；uij为数据点xi到聚类中心cj的隶属度，数据点xi到所有聚类中心的隶属度总和为1，如式(6)所示:

(6)

dij为数据点xi到聚类中心cj的欧几里得距离，如式(7)所示:

(7)

FCM算法的实现步骤:

步骤1:初始化隶属度矩阵U，使其满足式(6).

步骤2:输入灰度图像，根据式(8)计算N个聚类中心点cj，j=1,…,N.

(8)

步骤3:根据目标函数式(5)计算目标函数值，如果值小于某个确定的阈值或相对上次目标函数值的改变量小于某个阈值，则算法结束.

步骤4:根据式(9)计算新的隶属度矩阵U，再返回步骤2继续迭代循环，直至算法结束.

(9)

2.5.3 水质预测反演

根据分类后的二值图像，提取高光谱数据中水体区域的所有像素，表示为矢量矩阵W={w1,w2,…,wn}.其中，wn为一个列向量，表示第n个像素对应所有波段的光谱反射率值.为得到预测值矩阵M={m1,m2,…,mn}，其中mn为wn经过模型预测后获得的预测值，对W中每个列向量进行相应变换，并将新向量输入预测模型，预测水质指标在水体中的浓度含量，直至W中所有像素遍历完成.

在得到预测值矩阵M后，根据水质指标标准范围的最大值Vmax和最小值Vmin，将其划分为6个区间，每个区间显示的颜色分别表示为蓝色、天蓝色、绿色、绿黄色、黄色和红色.使用者可以根据自己的需求调整区间的大小，也可将每个区间设置为想要的颜色，最终使预测值矩阵M以不同颜色的图像显示，从而让使用者更直观地了解水质情况.

3 系统实现

利用Windows10 X64操作系统，在处理器为Intel(R) Core(TM) i5-10400 CPU@2.90GHz、内存大小为16 G、固态硬盘大小为512 M的硬件基础上，基于Visual C++开发工具的MFC框架[12]、OpenCV图像处理库和GDAL遥感数据处理库，开发基于无人机高光谱技术的水质预测反演系统.

3.1 界面设计

该系统是MDI多视图界面，菜单栏的“File”菜单包括“Open”打开高光谱图像、“CLose”关闭当前高光谱图像和“Exit”退出系统；“基本操作”包括“标记采样点”和“关闭标记采样点”的开与关两个操作，以及“保存采样点光谱反射率值”“水体提取”“反演浓度分布图”等功能；“数据处理”包括“回归分析”操作.主界面如图3所示.

3.2 功能实现

在打开“标记采样点”后，鼠标在图像区域时会变为十字形，其他区域为手型.点击左键可标记实际采样点的位置，点击右键可取消已标记的上一采样点，点击“取消标记采样点”可关闭鼠标左右健功能，鼠标变为手型.“保存采样点光谱反射率值”用于将标记采样点处的光谱反射率值导出.标记结果如图4所示.

图3 水质预测反演系统主界面Fig.3 The main interface of the water quality prediction and inversion system

图4 标记采样点结果图Fig.4 The result graph of marking sample points

“高光谱图像水体区域提取”功能充分利用水体区域在绿光波段和近红外波段反射率的强反射和强吸收特性，用归一化差异水体指标增强水体与周围非水体目标的差异，利用模糊C-Means聚类方法实现水体与非水体的分类，并提取水体区域轮廓.操作界面如图5所示.

图5 高光谱图像水体区域提取操作界面Fig.5 Operation interface of water extraction in hyper-spectral image

“反演浓度分布图”功能主要利用已训练的水质模型对高光谱图像水体区域中的光谱数据进行预测反演，根据水质指标的浓度区间范围对反演预测值指定颜色，并进行可视化显示，以呈现水质指标在水体中的总体浓度分布图.PCA模型对水体TN指标的浓度反演结果如图6所示.

图6 水体区域TN指标浓度分布图Fig.6 Concentration distribution graph of TN indicator in water region

在菜单栏“数据处理”下的“回归分析”功能，支持采样点光谱数据和水质指标实测数据的导入，以及光谱数据与水质指标实测值的相关性分析；集成多种常用水质指标、光谱数据分析方法，以及多种常用回归分析模型，包括PLS、PCF、RBF、MLP和SVM等；支持模型的训练和保存，以及模型的加载和测试等.功能界面如图7所示.

图7 光谱数据分析及模型训练Fig.7 Spectral data analysis and model training

3.3 系统测试

水质反演预测系统是一款单机版应用软件系统.在软件开发过程中每开发完成一个功能都要做相应的单元测试，以确保所实现功能模块的正确性和可行性；在完成系统开发后，根据实际需求对每个功能(如标记采样点、保存采样点光谱反射率值、水体提取、水质指标浓度反演等)做功能性测试，以保证软件各项功能实现预期的效果，完成水质预测反演的整个操作流程，达到完全界面化的操作需求.

4 结 论

本文开发设计了一套基于无人机高光谱技术的水质预测反演系统，为无人机高光谱技术在水质预测方面的应用研究提供了一款有效的工具，为实时预测反演水质指标提供了可视化的操作界面.目前，该系统只能对常用的8种水质指标(总氮、总磷、叶绿素a、高锰酸钾、总悬浮物、氨氮、化学需氧量、pH值)进行监测，后期研究将对该系统的光谱数据处理、水质指标和回归分析模型等方面进行扩展和完善.