基于高光谱的水体BOD含量模拟估算

王洪伟，王 波，纪 童，徐 君，剧 锋，王彩玲

1. 武警工程大学，陕西 西安 710086 2. 盐池县草原实验站，宁夏 盐池 751506 3. 甘肃农业大学草业学院，甘肃 兰州 730070 4. 西安航空学院，陕西 西安 710077 5. 中华人民共和国银川海关，宁夏 银川 750000 6. 西安石油大学，陕西 西安 710065

引 言

随着人类物质生活水平的提高和工业化的发展，水污染已经成为当今社会普遍存在的问题，其监测与治理也备受关注。生化需氧量(biochemical oxygen demand,BOD)是水体中的好氧微生物在一定温度下将水中有机物分解成无机质，这一特定时间内的氧化过程中所需要的溶解氧量，是监测水中有机物染物的一个综合指标[1]，是地表水、生活污水及绝大多数工业废水的必测指标之一。BOD值越高表明水中溶解氧会被自身微生物消耗的数值越高，造成许多的生态问题[2]。“五日培养法”为现下普遍的接受测定BOD的方法，但测定时间长、不能及时反映水质变化，不适合现场监测。

自20世纪70年代以来，随着遥感技术的快速发展，高光谱技术已成为现代遥感技术的重要组成部分[3]，利用高光谱技术反演水质指数早有研究，刘彦君等[4]利用多光谱数据，进行线性与非线性模型反演研究，对浙江农林大学东湖水体的总磷(TP)、浊度(SS)、悬浮物浓度(TUB)进行了反演。林剑远等[5]利用水质化验数据和光谱反射率进行相关性分析，建立了浙江省嘉兴市河网化学需氧量(CODcr)、生化需氧量(BOD5)、总磷(TP)、总氮(TN)的反演模型。周亚东等[6]利用GF-1号WFV遥感影像，通过多元线性回归和RBF神经网络模型建立了武汉市周围水域综合营养状态指数模型。这些成果有效解读了水体光谱特征规律，为遥感监测水质，生产生活提供了理论支撑与技术指导。

原始光谱反射数据有着数据量大，指标彼此高度相关的特性； 原始指标高度相关的特性经常会导致多重共线性问题的产生，从而导致模型失真[7]，因此如何对大量光谱数据进行处理和挑选一直是光谱反演模型的重点。主成分分析法(PCA)与偏最小二乘法(PLS)作为常用降维方法在遥感上应用广泛[8]，许多研究结果也表明应用主成分分析与偏最小二乘法筛选的主成分参数可以更好的反演各自的指标。杨国范等[9]利用比值线性回归模型与最小二乘支持向量机，对铁岭清水河库叶绿素a浓度与Landast OLI卫星数据分析，并建立了叶绿素浓度a的反演模型。何金成等[10]利用近红外光谱数据结合偏最小二乘法回归建立了BOD预测模型。

现有文献报道中，利用光谱估测水质参数BOD指标的报道较少，基于此试验利用光谱数据进行水体指标BOD的反演，测定水体样本光谱数据的同时收集水体样本并带回实验室测定BOD指标，将采集到的光谱数据与BOD指标进行Person相关性分析，挑选敏感光谱指标； 由于光谱指标之间的高度相关，为避免模型失真，在建立反演模型之前，利用主成分分析和偏最小二乘法分别对光谱指标进行处理，消除指标之间的多重共线性问题，最终建立多元线性回归模型与偏最小二乘模型，比较两种建模方法的建模精度与预测效果，选出更加适合反演BOD指标的建模方法。探索利用高光谱技术估测水体BOD值的可行性与最优方法，为实时诊断水体状况提供理论基础和关键技术，为实现对BOD指标实时监测提供可行的途径。

1 实验部分

1.1 试验地概况

于2018年对西安市地表水环境进行取样研究，取样地点集中于渭河(林家村)、浐河(田家湾)、灞河(马渡王)，共计60处采样点，每处采样点共计10次重复。

1.2 方法

1.2.1 光谱数据获取与校正

所用仪器为美国ASD (Analytica Spectra Devices.,Inc)公司制造的适用于遥感测量、农作物监测等方面的 FieldSpec®4 Hi-ResASD便携式地物光谱仪，其光谱范围为300～2 500 nm。

光谱采集选择干燥、无风、晴朗无云或少云的天气进行，并根据天气条件及时进行标准白板校正，采集时间尽量在10：00—14：00之间，此时光照条件良好。进行地面水质采样和水体光谱数据等实验数据获取，光谱采集参数设置时间为100 ms，测量后及时进行白板校正[11]。每块样本选择2～3个光谱采样点进行高光谱数据采集，每个样点每次重复测量10次，最后以该样点的光谱反射率均值制作光谱反射率曲线。

1.2.2 BOD指标的测定

采集水体样本时，利用聚乙烯桶采集距离水面10～12 cm的水样，不使漂浮于水面的物质混入，每处试验点共取10次样本，对水样加入保存剂，以便将样本带回实验室，利用标准稀释法[12]处理水样，并在20 ℃培养箱中培养，5 d后测出培养后的溶解氧含量，取平均值作为BOD指标原始因变量。样本BOD参数变化范围如表1所示。

表1 水质参数变化范围Table 1 Variation range of water quality parameters

1.3 数据处理

普通的多元线性回归应用中有许多限制，最典型的问题就是自变量之间的多重相关性。为此，利用主成分分析降维与消除指标间多重共线性的特性，筛选多元线性模型的自变量，已期解决多重共线性对参数估计的影响，减小模型误差。偏最小二乘回归中开辟了一种有效的技术途径，通过对系统中的数据信息进行分解和筛选，提取对因变量的解释性最强的综合变量，辨识系统中的信息与噪声，从而能够更好地克服变量多重相关性在系统建模中的不良作用。

2 结果与讨论

2.1 Person相关性

图1为原始光谱DN值与水体BOD含量的相关系数图，因图中波段1 023～2 500 nm与水体BOD含量无显著相关性，因此图中只展示了350～1 023 nm波段范围内的相关系数，由图可知BOD指标与光谱在350～900 nm呈负相关，960～100 nm为正相关，350～490与920～1 000 nm与BOD指标无显著相关性，BOD敏感波段大体分布于600～900 nm，其中758 nm处为相关系数绝对值最大值0.418，根据相关系数大小与显著性原则，筛选出了35个与BOD指标极显著相关的原始光谱指标，作为多元线性回归模型与偏最小二乘模型的自变量，筛选指标相关系数绝对值由大到小分别为： 758，759，853，809，1000，810，890，813，851，1 012，807，893，618，864，816，806，782，787，785，888，796，808，924，845，663，530，887，724，863，889，757，683，628，909和689 nm。主成分分析要求建模数据量高于变量数，偏最小二乘法允许在样本点个数少于变量个数的条件下进行回归建模，Person相关系数法共筛选出35个光谱变量，因此将60组样本数据分为建模组(40)与检验组(20)。

2.2 主成分分析

主成分分析结果如图2所示。

图1 相关系数图

图2 主成分分析碎石图注： 横坐标是主成分，纵坐标为解释程度Fig.2 Principal component analysis lithotripsy

经分析共有10个主成分。其中主成分1方差贡献率为94.9%，主成分2方差贡献率为1%，而主成分3～10累积方差贡献率不足10%，且主成分2到主成分3，斜率开始趋于平缓，因此剔除主成分3～10，只保留主成分1和2(Z1和Z2)，这2个主成分既能达到降维的目的，又能反映原始数据95.9%的信息。

图3直观展现了各植被指数在主成分1和主成分2中的分布情况。横纵坐标分别代表第一主成分与第二主成分以及各自的贡献率，Z1和Z2累计贡献率高达95.9%，可以解释原有变量中的大部分信息，4个BOD含量分组中0～0.2与0.4～0.6 mol·L-1在4组中彼此独立，可以明显区分，0.2～0.4与0.6～0.8 mol·L-1彼此交叉分组不明显。

将特征向量代入主成分公式中，得到主成分Z1和Z2的表达式

Z1=0.168x1+0.168x2+0.17x3+…+0.166x34+0.169x35

Z2=0.382x1+0.353x2+0.287x3+0.28x4+0.151x5+0.122x7-0.311x26-0.264x27-0.208x29-0.364x30-0.196x32-0.271x34-0.13x35

将2个主成分分别代入多元线性回归中，得到的方程

YBOD=-0.000 004 468z1+0.000 059 19z2+9.217

(R2=0.656，RMSE=0.007)

多元回归模型中BOD与主成分拟合方程R2较大，RMSE值较小，说明利用主成分Z1和Z2通过多元线性回归，可以很好的拟合水体BOD指标。

图3 主成分分析效果图

2.3 偏最小二乘模型的构建

由于自变量与因变量之间的量纲与数值都是不同的，现将BOD值与筛选的光谱指标进行标准化处理，利用R语言PLS偏最小二乘函数包建立水质BOD含量的估测模型，各主成分贡献率结果见表2。

表2 主成分贡献率Table 2 Contribution rate of principal component

如表2可知，当主成分为3时，解释率逐渐趋于平稳，因此选取comps=3时建立模型。y=0.015 703x1+0.124 092x2+0.423 545x3-0.181 04x4-0.255 47x5+…+0.331 165x34-0.189x35(R2=0.896，RMSEP=0.7469)。

使用函数包中jack.test函数对回归系数进行显著性检验结果见表3。

表3 jack.test函数显著性检验Table 3 Significance test of jack.test function

通过jack.test函数进行显著性检验，表3中“*”代表极显著影响，“**”代表显著影响，x1—x35代表波长按由小到大排列的原始光谱变量，由表3可知对水体BOD含量有显著影响的光谱指标有628，889和893 nm，其中对BOD有正向影响的光谱指标为628与889 nm，对BOD有负向影响的光谱指标为893 nm。

2.4 最优模型筛选

比较多元线性回归模型与偏最小二乘法模型，依据R2最大RMSE最小原则，最终采用偏最小二乘法模型y=0.015 703x1+0.124 092x2+0.423 545x3-0.181 04x4-0.255 47x5+…+0.331 165x34-0.189x35(R2=0.896，RMSEP=0.746 9)。

2.5 模型精度检验

将检验组的20组BOD与光谱数据代入模型进行模型精度检验，检验结果见表4。

表4 模型精度检验Table 4 Model accuracy test

从表4可以看出，偏最小二乘模型，其均方根误差较低为0.12，且估测精度R2较高。说明利用偏最小二乘法可以建立精度较好的BOD反演模型。

3 结 论

利用多元线性回归与偏最小二乘法建立水质BOD指标的反演模型。在进行光谱单波段与BOD指标相关性分析时发现在350～500 nm波段相关系数偏低，350～500 nm原始光谱反射曲线杂乱，说明该波段可能受其他水质参数影响，波段敏感性较差，不能作为模型预测波段。随着波长增加相关系数于758 nm达到最高值，且光谱最优反演波段大多分布在600～900 nm处，与林剑远等[5]得到的高光谱遥感数据与BOD指标敏感波段750～900 nm有一定不同但也有相似之处，主要原因有以下几点：

(1)水体光谱受自然条件与人为干扰，使边缘波段噪声很大，导致350～500 nm波段敏感性较差，与BOD指标的相关性较低。

(2)水体光谱易受时间空间等影响，导致光谱区别较大，但光谱趋势整体相似，且受其他水质指标影响，其光谱也会随之变化。

高光谱具有分辨率高，波段连续性强的特点，但光谱信息冗杂，数据的筛选与模型的简化一直是光谱模型研究的重点[13]，通过主成分分析和偏最小二乘法综合筛选的光谱指标，建立了多元线性回归模型与偏最小二乘法，结果表明主成分分析与偏最小二乘法可以有效降低数据维度，综合筛选指标特性，提高光谱数据与BOD参数的相关性与模型精度，其中偏最小二乘模型模型精度远高于多元线性回归模型，因为偏最小二乘法是分别从因变量与自变量中提取成分因子，保证成分因子尽可能多的反应变量的变异信息，同时也保证了两者之间相关性最大[14]，且试验中样本个数与变量个数大致一致，适用于偏最小二乘法模型。

在拟合偏最小二乘模型时，利用jack.test函数得出对水体BOD含量有显著影响的光谱指标有628，889与893 nm，说明628，889与893 nm可以作为反演BOD指标的敏感波长。林剑远等[5]以高光谱数据研究是城市河网BOD指标中发现波段565 nm为单波段与BOD指标相关系数(0.44)最佳波段，689/667 nm为组合波段与BOD指标相关系数(0.84)最佳波段，与本工作筛选的敏感波段有所不同，但有所相近。

以上试验结果为水质BOD指标的快速估算提供了依据，也为水体质量评估提供更便利的方案。