光谱二阶微分波动指数的河道水质类型识别方法研究

李 澜，田 华，季铁梅，巩彩兰，胡 勇，王歆晖，何志杰

1. 中国科学院上海技术物理研究所，上海 200083 2. 中国科学院红外探测与成像技术重点实验室，上海 200083 3. 中国科学院大学，北京 100049 4. 上海市水文总站，上海 200232

引 言

城市河道是城市生态环境的重要组成部分，改革开放以来我国城市经济快速发展，但也造成河道出现不同程度的污染[1-2]。城市河道水污染治理的前提和关键是水质状况的实时监测； 目前水质监测方法以传统采样监测为主遥感监测为辅，传统的水质监测通过实验室化验确定水质指标浓度再根据相关标准确定水质级别； 遥感监测方法是通过卫星、飞机、无人机等平台搭载的传感器获取河道区域的遥感影像，再经过一系列处理分析建立水质参数与光谱反射率之间的定量关系，进而确定河道的水质级别。传统采样监测技术已较为成熟，但存在费时、费力、只能获得断面点位水质等不足。随着遥感和高光谱技术的发展，可以准确、快速地获取各种地物的各波段的光谱信息，国内外学者提出了众多基于遥感反演水质参数进而判断水质类型的研究方法[3-6]。王丽艳等使用MODIS数据反演了呼伦湖水体COD浓度并确定水质级别[7]。Wang等提出了一种综合水质评价方法，并利用该方法反演了多个水质参数浓度[8]。但河道水质类别的划分指标如溶解氧(DO)、总磷(TP)、氨氮(NH3-N)等不具有明显的光学特性，通过遥感反演各水质指标进而判断水质类型这类方法的局限性在于最终级别的识别精度受各个水质指标反演精度共同决定，单一指标反演精度不高必然整体识别精度低下，模型应用局限性大、泛化能力不强，适用性不足。

水体反射率是诸多水质参数综合作用下的结果，水体的光谱特征与水质类型之间存在某种映射关系。提出一种基于光谱二阶微分波动指数的光谱分析方法，能有效反映各波段光谱二阶微分的波动情况，在此基础上可实现城市河道水质类型的快速识别，可推广应用于利用其他遥感平台获得的高光谱数据监测内陆水体的水质类型。

1 实验部分

基于光谱二阶微分波动特性的水质光谱特征提取及分类方法，其基本思想在于利用上下包络线来提取光谱二阶微分曲线的振动范围，具体的数据处理流程如图1。

图1 基于光谱二阶微分波动指数计算流程图Fig.1 Flow chart of second-order differentialfluctuation index calculation

(1)光谱二阶微分计算。光谱二阶微分计算包括差分计算和平滑两部分，首先利用光谱差分计算原始光谱的二阶微分，随后使用Savitzky-Golay Smoothing法[9]对原始二阶光谱进行平滑处理，消除噪点和其他干扰。

(2)局部极大值极小值提取。曲线的波峰和波谷通常对应局部极大值和极小值，而波峰波谷的值可以反映曲线局部的变化范围，使用3×3的滑动窗口逐点提取局部极值，满足式(1)或(2)的即为局部极大值或极小值点。

(1)

(2)

(3)无关极值点去除。曲线上的波峰波谷对应着局部极值点，通常光谱二阶微分曲线某一波段附近只会对应一个峰值，由于仪器探测性能限制或其他原因的影响，峰值周围会出现一系列无关极值点，这会对后续处理产生干扰，为此可设定一个距离阈值m，采用循环迭代方式去除无关极值点，下面以极大值为例介绍具体实现步骤。①寻找当前所有极大值中最大值以及对应的波长位置； ②去除最大值周围距离小于阈值m的所有极大值点； ③弹出并保存最大值和波长位置。重复①②③，直到所有的无关极大值点均消除。无关极小值的去除与上述步骤类似，只需每次弹出最小值和位置即可。

(4)双包络线生成。在上一步处理后可得到光谱二阶微分曲线上的各波峰点和波谷点，使用三次样条插值法得到光谱二阶微分曲线的上下两条包络线。图2是某光谱二阶微分及双包络线示意图。

图2 光谱二阶微分曲线与双包络线Fig.2 Second-order differential curve and double envelope

(5)光谱二阶微分波动指数计算。上下包络线能表示光谱二阶微分的振动范围，而某一波段处的对应上下包络线函数值之间的差值能表示该波段附近光谱二阶微分变化的剧烈程度。因此可定义一种光谱二阶微分波动特性指数来描述这一特征，具体计算公式如式(3)

F(λ)=eA(f(λ)-g(λ))

(3)

式(3)中，F(λ)表示波长为λ处附近的变化特征，f(λ)和g(λ)分别对应该波段上下包络线的函数值，A是固定的放大系数用于增加区分度。

1.1 水体光谱采集

先后多次在上海市嘉定区采集了各水质类型水体的光谱数据，光谱数据获取采用ASD公司生产的FieldSpec3便携式地物光谱仪，波长范围选择为400～900 nm。测量和处理参照国家标准《水体可见光-短波红外光谱反射率测量GB/T 36540—2018》进行。图3是获取的水体反射光谱曲线图。

图3 反射光谱曲线图Fig.3 Reflectance spectral curves

1.2 水质采样与处理

采集光谱数据的同时，还需要获取相应的水质数据。将现场采样的水体样本密封保存带回实验室化验获得各水质指标，再根据《地表水环境质量评价办法(试行)》等标准划分类别，共得到Ⅰ和Ⅱ类样本11个(由于Ⅰ类样本稀少，将其与Ⅱ类样本合并)、Ⅲ类14个、Ⅳ类20个、Ⅴ类23，劣Ⅴ类20个。

2 结果与讨论

2.1 特征波段选择

采用上述的数据处理方法计算各类水体的光谱二阶微分波动系数，合适的距离阈值和放大系数可以有效去除虚假极值点，增加不同类别之间的区分度，经多次试验最小距离阈值m设置为10，放大系数A设为500效果较好，图4是计算的各类水体的平均光谱二阶微分波动指数曲线图。

图4 各类型水体的平均波动指数Fig.4 Average fluctuation index of each type of water body

从图4可以看出，各类水体波动指数在720～740，750～770以及820～840 nm处均有明显的峰值，且不同水质类型水体峰值高度不同，随后统计了各类样本中这三个波段内的平均波动指数的均值、标准差等特征，结果见表1—表3。

上述统计结果表明，水体在720～740，750～770和820～840 nm内的光谱二阶微分平均波动系数与水质级别具有正向相关关系，随着水质级别的上升，这三个波段内的平均波动指数的平均值和中位数也随之上升，且水质级别越高上升幅度越大，变化越剧烈，同样平均波动指数的标准差也在增加。因此可利用这一特性实现河道水质类型的自动识别分类。

表1 各类型水体720～740 nm平均波动指数统计Table 1 Average fluctuation index from 720～740 nm

表2 各类型水体750～770 nm平均波动指数统计Table 2 Average fluctuation index from 750～770 nm

表3 各类型水体820～840 nm平均波动指数统计Table 3 Average fluctuation index from 820～840 nm

2.2 水质类型识别模型建立及验证

首先将各类样本按照2∶1的比例随机划分为测试集和训练集，随后采用本方法计算训练集内各样本的二阶微分波动指数曲线，得到三个特征波段内的平均波动指数，最后以这三个特征波段的平均波动指数为输入，使用LSSVM[10]构建水质类型识别模型，训练完成后进行测试，取重复10次实验的平均值作为最终结果，具体结果见图5(a)。

图5(a)显示平均识别准确度达80.65%，而通过对误分类的样本进行分析后发现，误分类样本多集中在两类样本特征边界处，故单独测试了分类误差不超过1类的识别准确度，重复测试10次，如图5(b)所示。结果表明误差不超过1类的平均识别准确度达96.77%，在31个测试样本中，仅1个样本分类误差超过1类，因此光谱二阶微分波动系数结合机器学习模型可实现城市河道水质的快速识别，且识别误差不超过1类。

3 结 论

针对目前城市水体大面积水质类型调查的需求，提出了基于光谱二阶微分波动指数的高光谱数据处理和分析方法，在对各类型水体的波动指数分析后得到如下结论: (1)720～740，750～770和820～840 nm波段内的平均波动指数与水质级别具有正相关关系，随着水质级别的上升，平均波动指数值升高，对应波段附近的光谱二阶微分振幅越大，这表明水质状况越差的水体，其光谱二阶微分波动越大。(2)以三个特征波段的平均波动指数作为输入特征，采用LSSVM模型构建水质类型识别模型，经测试，模型的平均识别准确度达80.65%，误差不超过一类的识别准确度达96.77%，能较准确地识别城市河道水质类型。