改进SVR的内陆水体COD高光谱遥感反演

盛 辉，池海旭，许明明*，刘善伟，万剑华，王锦锦

1.中国石油大学(华东)，海洋与空间信息学院，山东 青岛 266580 2.珠海欧比特宇航科技股份有限公司，广东 珠海 519080

引 言

河流及水库作为内陆水资源的重要组成部分，对涵养水源、保护生态具有重要的作用。但近几十年来，自然因素和人为因素的双重压力使水环境在不同程度上受到了污染，水体富营养化的趋势愈加严峻，对水质进行动态监测和精细化的管理迫在眉睫。

传统的水质监测工作主要是获得监测区域内某些点的水质浓度信息，这些有限点的信息并不能反映出整个区域的情况，而遥感水质监测具有实时高效、监测范围广和适合于长期动态监测等特点[1]。化学需氧量(chemical oxygen demand，COD)可以用来表征水中有机物的含量，是遥感水质监测的重要指标。COD遥感监测主要是研究水体反射的光谱特征和COD浓度之间的关系，从而建立反演算法。水质参数反演的方法主要有两类：半经验方法和半分析方法[2]。一些学者利用半经验的分析方法研究COD和遥感参数之间的关系，通过构建的线性及生物光学等模型反演得到了内陆湖泊COD含量的分布情况，并进一步分析得到了COD的来源及影响[3]。但半经验模型的拟合精度范围大多在0.6～0.8之间，因此模型的准确性还有待提高。国内外的部分学者研究评估机器学习进行水质遥感的可行性[4-5]，表明在充足样本的情况下，人工神经网络能够对水质参数进行有效的反演分析。由于某些机器学习算法参数选择的不确定性以及COD光谱特征的微弱性，通常需要多次重复实验来获得最优参数从而提高模型的准确性。另外，研究人员在野外实验中需要花费大量的时间和精力才能够收集到足够的水质参数数据。因此，构建结构简单、适合小样本的水质参数反演模型尤为重要。研究表明，经过优化的支持向量回归机(support vector regression，SVR)具有结构简单、全局最优的特点，已经广泛应用到叶绿素、悬浮物、溶解氧和水体富营养化的水质研究中[6-7]，但是针对COD的反演研究较少。因此，本工作采用改进的SVR方法进行COD反演。

用于COD反演的数据源主要集中于Landsat，Modis和Sentinel-2等多光谱数据[8]。高光谱数据可以捕获由不同水质浓度引起的光谱变化，并在内陆水质监测中显示出巨大的潜力[9-10]。但由于水体光谱影响机理的复杂性和其他水质参数的影响[11]，狭小区域的COD反演和变化监测需要兼有高空间分辨率和高光谱分辨率的数据提供支持。实测的水表面光谱反射率具有受大气影响较小的特征[12]，将其与高光谱图像结合可以更准确地估算COD浓度[13]。林建远等通过航空高光谱数据和实测的水表面反射率反演得到了杭嘉湖平原河道的COD含量，但是航空高光谱数据的获取难度大、成本高。珠海一号高光谱卫星(orbita hyper spectral, OHS)数据具有高空间分辨率、高时间分辨率和易于获取的特点。因此，采用现场实测光谱反射率数据结合珠海一号高光谱卫星数据用于COD浓度的反演。

针对COD光谱特征的微弱性以及SVR参数选取困难、易陷入局部极值的问题，基于高光谱数据，提出了一种模拟退火—粒子群算法(simulated annealing-particle swarm algorithm，SA-PSO)优化的SVR方法进行COD浓度反演。首先，基于实测的光谱数据建立潍河—峡山水库区域的COD高光谱反演模型，减少大气校正对反演精度的影响。然后应用于OHS高光谱数据对该区域的COD含量进行遥感估算，以期为潍河流域的水质监测与综合管理和OHS卫星数据在内陆水质的反演应用提供基础数据与科学参考。

1 研究区域与数据

潍河及峡山水库位于山东半岛，是胶东地区的战略水源地之一。潍河是流经峡山水库的主要河流之一，发源于沂水县和莒县，最终注入渤海，河长200多公里，流域面积近万平方公里；峡山水库的注入河流有潍河、渠河和浯河等河流，入库量为8×108m3左右。峡山水库周围及潍河上下游区域分布着众多村庄和工厂，该河域为周边地区的居民饮水，农业灌溉和工业发展做出了重要贡献。同时，该地区的百姓生活与经济发展也对水质造成了一定的影响，采用遥感技术对该地区的水质进行监测显得尤为必要。

1.1 实测数据

2019年10月26日—10月28日，分别在潍河流域的古县大桥、峡山水库、辉村和金口大桥附近水域进行水体取样和光谱采集，采集光谱数据12组，28日有效测得COD浓度数据 22组。将采水器采集到的表面水样装在琥珀瓶中，在瓶体贴标签记录好采样地点、时间和当时的天气情况并拍照保存，最终在实验室化验得到各采样点COD的浓度值。

图1 潍河及峡山水库现场实验点位图

光谱测量采用YW-TRIOS-AWRMMS水面移动测量系统(光谱范围320～943.7 nm, 共190个波段)，该仪器内置GPS定位系统，携带方便可用于各种情况的光谱测量。测量过程中，保持伸缩杆方向和太阳入射平面135°的夹角可使3根传感器分别有效测量太阳辐照度、水面辐照度和天空辐照度三种参数，并可在显示屏上查看光谱曲线。在每个采样点测量10组光谱以最小化不确定性，通过光谱数据的处理程序及Mobley(1999)的方法对光谱数据进行处理，即可得到离水辐亮度和遥感反射率等光谱数据信息参数。

1.2 遥感数据

OHS卫星于2018年4月26日成功发射升空，弥补了我国在高光谱数据上的不足，开启了商业航天遥感的新时代。OHS卫星空间分辨率为10 m, 4颗高光谱卫星的重访周期为2.5 d，单次成像范围为150 km×320 km，光谱分辨率为3～8 nm，光谱范围为466～940 nm, 共32个波段。由于其具有幅宽大、高空间分辨率、高光谱分辨率和高时间分辨率的特点，非常有利于反演光学特性复杂多变的内陆水体、小型水库和河流的水质参数等相关研究。OHS数据是当前能够进行小区域水质参数反演的最佳数据源。

采用珠海欧比特宇航科技股份有限公司提供的2景潍河流域2019年10月28日的OHS高光谱数据。首先对影像数据进行辐射定标、大气校正、影像拼接与陆地掩膜等处理，最终得到仅有潍河流域区域在内的影像产品。大气校正是为了消除大气分子和气溶胶的影响，是准确获得水质遥感信息的前提。FLAASH模型用于遥感影像的大气纠正具有很好的应用效果；故使用FLAASH模型进行大气校正，然后根据采样点的坐标位置，在影像上提取各点对应的反射率信息。

2 改进SVR的内陆水体COD高光谱遥感反演模型

建立COD的反演模型，首先要对光谱数据进行归一化和相关性分析等预处理，然后基于对COD敏感的波段结合SA-PSO算法建立SVR反演模型。

2.1 光谱数据与COD的相关性分析

为了削弱野外环境给光谱测量造成的影响，同时也使光谱数据和OHS影像光谱范围相对应，采用归一化的方法对466～940 nm的光谱数据进行预处理。归一化的公式如式(1)

(1)

式(1)中：RN(λi)为水体遥感反射率的归一化结果，R(λi)为原始的遥感反射率，λi为i处波长，n为466～940 nm处的波段数。

Pearson相关性分析是从统计学的角度研究两个或多个随机变量间关联度强弱的方法。相关系数的大小可以描述变量间的密切程度，两个变量的相关系数表达式如式(2)

(2)

式(2)中：yi是各个采样点水质参数的浓度，xi是各个采样点的地表反射率。

2.2 基于SA-PSO算法优化的SVR模型

SVR是基于机器学习的一种算法，通过核函数将非线性数据映射到高维空间构造决策函数进行线性回归，在解决小样本、非线性和高维模式识别问题方面具有其独特的优势，经常用于小样本水质参数的反演[14]。SVR模型的数学关系为式(3)

f(xi)=ωTφ(xi)+B

(3)

式(3)中，xi为样本数据，ω为待辨识的权重，φ(xi)为非线性映射，f(xi)为特征空间中的线性函数，B为常数项。

由于SVR模型存在参数选取难的问题，如何确定最优参数直接影响到SVR模型的学习和泛化能力。根据鸟群觅食行为提出的粒子群算法PSO常用于对SVR模型参数进行优化选取，将只有位置和速度两个属性的粒子模拟为鸟，每个粒子在空间中的极值Pbest作为潜在的最优解Gbest; 所有粒子根据式(4)不断调整各自的位置和速度，直至获得最优解。同时，在粒子群更新过程中引入模拟退火算法SA可提高PSO算法的全局寻优能力，即利用其在寻优过程中的突跳能力对粒子群算法进行改进。SA算法本质是模拟高温物体退火过程中寻找全局最优解的过程。初始温度T的确定采用适应度和和接受概率的方法，由式(5)决定

Vi+1=ω×Vi+c1×rand×(Pbest-Xi)+

c2×rand×(Gbest-Xi)

Xi+1=Xi+Vi+1

(4)

T=(fmax-fmin)/lnp=-|Δf|/lnp

(5)

式中：fmax和fmin和Δf为初始粒子群最大、最小适应度值及其差值，p为初始接受概率，p=0.8。

建立基于SA-PSO算法优化的SVR模型，关键是利用SA-PSO算法优化支持向量回归模型中判别函数的惩罚因子C和核函数中核的宽度g。建立SA-PSO优化SVR模型的主要步骤如图2所示。

图2 SA-PSO优化SVR流程

(1)粒子群的初始化。对粒子速度，位置，惩罚因子C和RBF核函数宽度等参数进行初始化。设置初始退火温度T，并将最大迭代次数设置为200。

(2)计算群体中每个粒子的适应度。如果新的适应性值好于原始值，则接受新位置。否则，将保留旧位置。

(3)根据式(5)更新粒子的位置和速度。

(4)收敛性的判断。如果满足终止条件，则算法停止。否则，执行退火操作，转到(2)。

(5)输出最优参数C和g。

3 实验与结果讨论

进行COD的高光谱反演，首先要对光谱进行归一化和相关性分析来确定敏感因子，然后根据敏感因子和COD浓度建立反演模型。最后通过使用与采样点相对应的OHS反射率数据来分析模型的准确性，并将SA-PSO-SVR应用于潍河流域的OHS数据进行COD的空间分析。

3.1 光谱数据与COD的相关性分析

实测的原始光谱曲线如图3(a)所示，根据等式(1)进行归一化的结果如图3(b)所示。与原始光谱曲线相比，归一化光谱曲线的形状发生了变化，反射谷和反射峰更加明显。从560，680和710 nm的反射率值可以发现，当浓度增加时，反射峰具有向短波长方向移动而反射谷向长波方向移动的趋势。

图3 原始光谱(a)和归一化后的光谱(b)

Pearson相关性分析的输入反射率形式包括单波段和波段比值等形式。研究表明, 波段比值组合可部分消除水表面光滑度和微波等其他环境因素的影响[15]，能在一定程度上有效提高水质参数反演的精度。对各个采样点的COD浓度分别与每组光谱数据的波段比值组合做Pearson相关性分析，得到COD浓度与波段比值组合的相关性分析结果如图4所示，最佳反演因子是518 nm/940.4 nm，663.6 nm/636.8 nm，729.2 nm/890.9 nm和752.3 nm/857.9 nm四种波段比值组合。

图4 466～940 nm的实测光谱比值与COD浓度的Pearson相关性分析结果

3.2 COD反演

根据选取的最佳反演因子，用12个采样点的光谱数据来训练SA-PSO-SVR模型。根据相关研究及反复实验，设置种群规模为50，进化次数为200，选取参数c1和c2均为0.5，SVR模型的C和g的取值范围为[0.1, 1 000]和[0.1, 10]。对训练数据进行归一化处理后，经过SA-PSO优化获得的惩罚因子C与核参数g分别为151.09和0.36。

利用28日测得的22个采样点数据作为检查点，在OHS影像上提取与采样点相对应的反射率，导入SA-PSO-SVR模型得到各个检查点的预测值。图5为SA-PSO-SVR精度分析结果。红色点为利用光谱数据建立模型得到的实测值和预测值的分布情况，虽然某些采样点的预测值与实测值具有一定差异，但依然可以看出SA-PSO-SVR模型的反演结果与实测值呈现较好的相关一致性，其模型决定系数为0.86。绿色点为检查点数据的预测值和实测值的分布情况，平均相对误差(MRE)和均方误差(RMSE)分别为9.04%和3.64 mg·L-1，说明该模型可以实现对COD浓度的有效反演。

图5 SA-PSO-SVR精度分析

3.3 SA-PSO-SVR精度分析

从图5可以看出，OHS数据部分反演结果与实测值差距略大，推测原因为数据质量参差不齐。由于TRIOS光谱测量仪和OHS卫星是两种截然不同的传感器，传感器的信噪比等参数并不一致，因此两种仪器获得的数据质量并不一致。研究表明，COD对水体的光学特性响应很微弱，为了减少河道两边其他非水体像元的影响，采样点的像元位置尽可能的选在距离监测点近的河道中央，但由于河道中央流速较快，可能导致河中央同采样点的COD浓度存在一定差异，进而造成部分反演结果与实测值的差异较大。

另外，基于光谱数据分别建立SVR和BP神经网络和线性回归(LR)模型，将SVR模型的惩罚因子C与核参数g设置为109.85和0.000 1；根据相关研究及反复实验，设置BP神经网络模型的隐含层节点为12，最大训练次数为2 000，期望误差为0.000 1，学习速率为0.01。SVR和BP神经网络和LR模型的决定系数分别为0.712，0.598和0.494。从图6(a)可以看出，SVR模型的整体反演效果较好，但是OHS数据的部分反演结果低于实际测量值。从图6(b)可以看出BP神经网络模型的RMSE和MRE分别为1.93 mg·L-1和5.18%，与其他两种模型的结果相差不大；但是从OHS数据的反演情况看，BP神经网络模型的结果陷入了局部最优，由此也证明BP神经网络模型确实存在参数选取困难的问题。另外，从图6(c)可以看出，光谱数据的RMSE和MRE分别为3.00 mg·L-1和8.53%，与其他模型结果相近，但OHS影像的反演结果要普遍高于实测值。综合来看，SA-PSO-SVR模型的反演效果优于其他模型。

图6 SVR, BP神经网络，LR模型预测评估

3.4 COD浓度空间分布分析

将SA-PSO-SVR模型应用于OHS数据，得到的COD浓度空戒分布如图7所示。可以看出，COD的浓度呈现部分区域浓度高的特点，在韩信坝、峡山水库的东北部区域、渠河与潍河的交叉口及辉村与金口大桥之间的浓度明显偏高。

图7 COD浓度分布情况

结合潍河流域地理情况并查阅相关资料发现，潍河流域附近分布着众多工厂，大部分企业将工业污水排入河中，这些河流(潍河、渠河等)携带污染物最终汇入峡山水库，是导致峡山水库东部沿岸区域的COD浓度明显高于其他区域的原因之一。另外，库区及潍河等其他河流两岸分布着众多村庄，百姓的农业生产活动依赖于该区域的水资源，农业生产带来的化肥、农药等污染物也是造成该局部区域COD浓度过高的原因之一。同时，研究表明农作物腐烂死亡降解产生的有机物也会造成局部区域的COD浓度增高。另外，在韩信坝和山阳村附近区域有多处拦河大坝，部分区域河流径流量减少，遥感影像获取的COD反射信息不充分，造成该区域的COD反演结果偏高。

4 结 论

针对SVR模型参数选取困难和易陷入局部极值的情况，本文基于水面以上的光谱数据建立了SA-PSO-SVR模型，通过珠海一号卫星数据反演得到了潍河—峡山水库流域COD的分布情况，得到以下结论:

(1)通过对潍河—峡山水库区域的实测光谱进行分析，结果表明该区域实测的水面光谱具有典型内陆湖泊水域复杂浑浊水体的光谱特征，560和710 nm附近的光谱曲线形状呈现明显的双峰特征，反射率幅值较大。当浓度增加时，反射峰具有向短波长方向移动而反射谷向长波长方向移动的趋势。

(2)引入SA-PSO算法对SVR模型的参数进行优化选取，解决了局部最优的问题，其反演效果明显也好于其他模型。将基于实测光谱建立的SA-PSO-SVR模型应用在珠海一号卫星数据上可以进行COD遥感估算，说明了基于实测光谱和高光谱卫星影像的水质参数反演方法具有良好的应用前景和推广价值。同时，反演得到的COD浓度分布情况可以为珠海一号卫星数据的内陆水质参数反演和潍河流域的综合管理提供科学依据。

珠海一号数据的验证结果说明了COD反演的有效性，但由于星载传感器和地面光谱仪信噪比等仪器参数的不同，并受到其他复杂的大气和光学水体组分的影响，该方法还需要有足够的论据对此进行进一步的详细分析。此外，虽然基于实测光谱建立的SA-PSO-SVR模型在潍河流域表现出了相对较好的性能，但在其他复杂的水环境的应用性能还有待研究。因此，接下来的工作将多次收集潍河流域的光谱、COD浓度等其他水体成分的信息，并根据长期观测资料进一步探讨COD的时空变化；同时，在包含多个内陆水域的综合数据集的基础上，对SA-PSO-SVR方法进行更全面的评估。