基于小波变换的土壤重金属铬含量定量反演

郭云开 ,谢晓峰,谢 琼,张思爱

(1．长沙理工大学 交通运输工程学院，湖南 长沙 410076；2．长沙理工大学 测绘遥感应用技术研究所，湖南 长沙 410076；3．湖南工程职业技术学院 测绘地理学院，湖南 长沙 410151;4．长沙理工大学 公路地质灾害预警空间信息技术湖南省工程实验室，湖南 长沙 410114)

耕地是人类生产生活的基本元素。随着工业化的快速发展，耕地土壤面临的重金属污染问题日益严重。在土壤中，重金属铬(Cr)主要以Cr3+、Cr6+的形式存在。Cr6+属于强致突变物质，可诱发肺癌和鼻咽癌。Cr3+主要以氢氧化物的形式附着于土壤组分中，不易迁移且毒性较低[1-2]。传统监测土壤重金属含量的方法主要是以点代面采集样本，带回室内进行化学分析。这种方法虽然精度高，但是费时费力且无法在大区域内进行实时监测。随着高光谱遥感技术的发展，高光谱遥感技术以其快速获取地面信息的能力和光谱分辨率高的特点为土壤重金属污染监测提供可能[4-5]。目前，在土壤重金属污染监测领域，很多专家已经采用高光谱遥感技术展开研究。文献[6]在对原始反射率光谱进行倒数对数、一阶微分、二阶微分和包络线去除处理等光谱变换，发现基于二阶微分处理后的数据建立的偏最小二乘模型(PLS)的预测精度高于其他方法。袁中强[7]等对比原始地表反射率、原始地表反射率的一阶微分和倒数对数等光谱变换方法，发现基于原始地表反射率的一阶微分值建立的回归模型对重金属铬的预测效果最好。龚绍琦[8]等在对土壤光谱反射率进行一阶微分、倒数对数和连续统去除等预处理的基础上建立逐步回归模型预测重金属铬含量，得出一阶微分为铬含量反演的最佳光谱变换。

通过微分处理等传统光谱变换方法，提取特征波段，建立预测模型，在一定程度上可以有效预测土壤中重金属铬含量。但是土壤中重金属的含量属于微量级，其光谱信息十分微弱。此外，微分处理反射光谱的方法在增强微弱信号的同时也会扩大噪声，降低光谱的信噪比[9]。在对微弱光谱信号的处理中，小波变换是一种较为理想的提升信噪比的方法[10]。因此，本研究采用Db4小波对原始光谱进行分解与重构，提取特征光谱，利用特征光谱与重金属含量的相关性建立PLS模型预测重金属铬的含量。

1 数据的获取与处理

1.1 研究区域概况

研究对象选取湖南省岳阳县某地区(见图1)。该地区地势平坦，土地肥沃，是湖南省粮食主产区。此外，该地区有很多中小企业。随着企业的快速发展，工业废水、废渣存在处理不当的现象，这使得该地区面临耕地污染等环境挑战。

图1 岳阳县某地区卫星图

1.2 数据采集及预处理

根据农田分布情况布设采样点，每个样本点采样范围为10 m×10 m。在采样范围内，采用5点采样法采集0～20 cm范围内的表层土壤，取约1 kg土壤装袋保存。将采集回来的土壤去除石块、杂草等明显杂质，放置在阴凉的通风室进行风干处理。风干后的土壤进行磨碎、过100目尼龙筛后将样本均匀分为两份。取其中一份进行重金属含量的测定。重金属含量测定采用电热板消化法，通过电感耦合等离子发射光谱仪测定。

使用Field-Spec AvaField-3高精度地物波普仪采集室内波谱。该仪器波段范围为：300～2 500 nm。室内光谱采集过程中，以聚四氟乙烯白板为标准白板，以1 000 W卤素灯作为光源，控制视场角为5°，光源与样本距离约30 cm，照射方向与垂直方向成15°。光谱采集前开机预热10 min，每个样品重复采集多条光谱曲线，取平均作为该样品的原始光谱。

样本采集总量为40个，剔除两个光谱异常的样本。此外，为保证数据精度，本研究去除了光谱中信噪比低的340～380 nm和2 401～2 515 nm光谱[11]。由于原始光谱分辨率高，波段之间存在大量冗余，数据量庞大。为降低数据量，在保证能充分反映土壤光谱特性的前提下，对原始光谱进行6 nm光谱重采样。

2 小波分解与重构

小波分解与重构是能够同时在时域和频域内进行局部化信号分析的数学方法。它具有多分辨率分析的特点，在时、频两域都具有表征信号细节特征的能力，适合分析非平稳信号和提取信号的局部特征[12-13]。土壤的原始光谱信号属于非平稳信号，适合采用小波变换提取光谱特征。

小波分解是将光谱信号分解为不同的层，各层所包含的信息与土壤不同成分有关。低频系数反映决定土壤光谱基本形状的物质，高频系数反映光谱采集过程的各种噪声及土壤微量成分在原始光谱的特征[14]。土壤中重金属含量属于微量元素，主要体现在高频系数中。小波重构是把分解的系数还原的过程，是小波分解的逆过程。小波单层重构是基于小波分解后单层小波系数进行重构，获得特定频率的重构信号[15]。

2.1 小波基选择及分解尺度的确定

目前，小波基的选择没有统一的标准，需要结合信号本身的特点和小波函数的数学特性综合考虑。本研究选择的Db4小波基函数具有正交性、双正交性、紧支性和不对称性的特点。同时，Db4小波基能够解决线性相位与正交性要求的矛盾，提取出与重金属含量具有较高相关性的局部高光谱信息[16]。

分解尺度的确立是小波分解的关键。分解尺度可以根据母小波函数和光谱信号长度共同确定[17]。确定小波分解的最佳分解尺度n的方程为：

n=1b(N).

(1)

其中N为高光谱信号的长度。

此外，最佳分解尺度n也可根据母小波的滤波器长度来计算确定[17]。本研究以Db4作为小波基，根据母小波滤波器长度计算得最佳分解尺度为4层。但是该理论计算得到的最佳分解尺度只是一个概率估计，可能存在偏差。因此，综合小波分解的特点及本研究数据特点，确定小波分解尺度为5层。

2.2 特征波段提取

以波段与重金属含量的显著相关程度作为判断特征波段的依据。在软件Matlab中编写程序实现小波分解和重构，使用重构后的信号与重金属含量进行相关性分析提取特征波段，提取的主要过程为：

1)采用Db4小波基对原始土壤光谱进行5层离散分解，获得5个高频系数和1个低频系数。

2)在Matlab中采用Wrcoef函数分别对高频系数和低频系数进行单层系数重构，获得不同的重构信号。

3)汇总各层重构信号，将1～5层高频信息分别命名为B1、B2、B3、B4、B5。由于属于痕量级的土壤重金属主要体现在高频系数中，同时为充分利用高频信息，研究将与重金属含量存在显著相关的各层高频信号进行相加，将得到的数据命名为Bn～Bm加和。

4)将得到的各组数据分别与重金属铬含量进行皮尔森相关性分析，通过0.05双侧显著性检验的波段作为特征波段。

3 实验与分析

3.1 相关性分析

如前述，有效样本总量为38个，随机选取27个样本作为建模集，11个样本作为验证集。将B1～B5各层重构信号、B1～B4加和、一阶微分处理后的土壤光谱分别与重金属铬含量进行皮尔森相关性分析，结果如图2—图8所示。

图2 B1层波段与Cr相关系数

图3 B2层波段与Cr相关系数

图4 B3层波段与Cr相关系数

图5 B4层波段与Cr相关系数

图6 B5层波段与Cr相关系数

图7 B6层波段与Cr相关系数

图8 一阶微分处理后全波段与Cr相关系数

统计各组数据的特征波段个数，结果如表1所示。单层高频系数重构信号中，B4层重构信号中特征波段的个数最多，达到22个。其次是B3层重构信号，特征波段的个数为13个。B5层重构信号不存在特征波段，说明在第5层重构信号中没有重金属铬的信息，可以视为噪声去除。基于一阶微分处理后的土壤光谱提取特征波段，共计11个。统计分析B1～B4各层重构信号与重金属铬含量都存在显著相关，因此选择B1～B4层的重构信号进行加和运算，获得B1～B4层加和数据。基于B1～B4层加和数据提取特征波段，结果存在24个显著相关波段。

表1 重金属含量与各层重构信号存在相关的波段个数

3.2 多元线性建模

选择B1～B4层加和、B4层及一阶微分处理后数据分别建立多元线性模型。建模结果如表2所示。

表2 不同的小波重构信号建立多元线性模型的结果

基于B1～B4层加和建立的多元线性预测模型的决定系数为0.916。基于B4层重构信号建立的多元线性模型的决定系数为0.339。基于一阶微分处理光谱信息建立的多元线性模型的决定系数为0.521。结合表2发现，B4层重构信号的特征波段的个数大于一阶微分处理后特征波段的个数，而基于B4层重构信号建立的多元线性模型的决定系数却小于基于一阶微分处理后的光谱信息所建立模型的决定系数。出现这样的结果，可能是因为作为自变量的特征波段之间存在严重的共线性问题。

方差膨胀因子(VIF)是衡量多元线性回归模型中共线性严重程度的一种度量。一般情况下方差膨胀因子大于10，表明模型中自变量之间存在很强的共线性问题[18]。分别对B1～B4层加和、B4层、一阶微分处理土壤光谱提取的特征光谱进行共线性诊断，结果如表3所示。分析表3可得，B4层、B1～B4层加和的特征波段，方差膨胀因子(VIF)远大于10，基于一阶微分处理提取的特征波段存在多个波段的VIF大于10。因此，基于B4层、B1～B4层加和及一阶微分处理后光谱数据提取的特征波段之间都存在共线性问题，其中，基于B4层、B1～B4层加和提取的特征波段之间共线性尤为严重。

表3 方差膨胀因子汇总

3.3 偏最小二乘模型建立

偏最小二乘回归模型是一种多因变量对多自变量的回归建模方法，同时包含多元线性回归分析、典型相关分析和主成分分析的特点，能够较好的解决变量之间存在多重相关性或样本点数过少的回归难题的方法[19-20]。为消除特征波段之间的共线性影响，建立可靠的预测模型，选择PLS作为建模方法。基于B1～B4层加和、B4层和一阶微分处理后的三组数据分别建立PLS预测模型，建模结果如表4所示。

分析表4可知，基于一阶导数处理后的数据建立PLS预测模型，在8个潜在因子、自变量累积解释能力为99.4%的前提下，PLS模型的决定系数R2为0.519。基于B4层重构信号进行PLS建模，在潜在因子为7个、自变量累积解释能力100%的前提下，模型的决定系数R2为0.481。基于B1～B4层加和建模结果为：在潜在因子为9个、自变量累积解释能力99.5%的前提下，模型的决定系数R2为0.605。基于B1～B4层加和建模结果比基于B4层重构信号建模的R2提升了0.124，比基于一阶导数处理后的数据建模R2提升0.086。结果说明，直接在原始光谱的基础上进行小波分解与重构，提取特征波段，利用特征波段与重金属含量之间的相关性建立PLS模型，可以有效预测重金属铬的含量。同时，与常用的一阶微分变换的方法相比，基于小波分解与重构提取的特征波段建立PLS模型的预测精度有所提升。

表4 不同小波重构信号建立PLS模型结果

3.4 模型验证与分析

根据建模集所建立的PLS模型系数，得到采用一阶微分光谱变换所提取的特征波段与重金属铬含量相关性建立的PLS模型为：

y=-635.527×χ1+3 546.074×χ2-

1 147.016×χ3-52.217×χ4+259.237×χ5-

3 093.044×χ6+1 036.373×χ7+

3 599.764×χ8-1 426.175×χ9+

1 923.961×χ10-41.321×χ11+11.406.

(2)

利用B1～B4层加和提取的特征波段与重金属铬含量的相关性建立的PLS模型为：

y=-654.654×χ1+646.323×χ2+

2 615.525×χ3-3 940.077×χ4+

4 670.818×χ5-1 957.401×χ6-

1 797.477×χ7-2 956.227×χ8+

3 810.497×χ9+1 645.746×χ10-

220.139×χ11+626.171×χ12-

1 107.56×χ13-982.036×χ14-

3 342.157×χ15-921.935×χ16

4 046.777×χ17+465.223×χ18-

3 928.37×χ19+1 141.567×χ20-

5 354.198×χ21-710.444×χ22+

3 132.64×χ23-1 518.24×χ24+10.416.

(3)

基于验证集数据分别验证以上两个模型，计算结果如表5所示。

分析表5可知，在验证集中，基于B1～B4层加和建立PLS模型的决定系数R2为0.511，比基于一阶微分建立的PLS模型的决定系数(R2)0.351提高0.16。结果表明：利用小波变换处理原始光谱建立偏最小二乘模型具有较高的精度。

表5 预测集建模精度汇总

4 结束语

本研究采用基于土壤重金属原始光谱进行小波分解与重构的方法，提取特征波段，建立PLS预测建模，对土壤重金属铬的含量进行定量反演,并与基于常用的一阶微分处理土壤光谱建立PLS模型进行比较。结果表明：

1)利用小波分解与重构提取的有用信息可以有效预测土壤重金属铬的含量；

2)对比基于常用的一阶微分处理提取的特征波段和基于小波变换处理所提取的特征波段，发现利用小波变换处理所提取的特征波段与重金属铬的含量相关性建立的PLS模型具有较高精度。

虽然小波分解与重构的方法可以有效提取光谱中有用信息，但本研究是基于室内光谱开展的，对于复杂的室外环境，小波分解与重构提取特征波段的方法是否适用仍需在今后研究中进行探讨。