基于VMD-MSE的玉米铜污染信息提取与预测模型

杨可明 李 燕 程 凤 高 鹏 张 超

(中国矿业大学(北京)煤炭资源与安全开采国家重点实验室， 北京 100083)

0 引言

近年来，矿山开采、污水排放、矿物肥料使用等造成的土壤重金属污染问题日益严重，给人类生产生活带来严重影响[1]。铜(Cu)是土壤重金属污染的主要元素之一，Cu在土壤中的累积对植被生长发育产生严重影响[2-3]。同时，通过食物链进入人体，将破坏人体神经系统和免疫系统[4]。因此，植被重金属污染的实时、有效监测显得尤为重要。研究表明，受重金属污染的植被光谱曲线会发生畸变[5-7]，因此很多学者利用高光谱遥感技术获取重金属胁迫下的植被光谱信息，旨在通过监测植被光谱特征变化来诊断重金属污染程度。

电磁波在传播过程中不可避免受到大气等影响而产生噪声，噪声对光谱数据中污染信息的解读造成极大干扰[8]，寻找一种有效提取作物污染弱信息的方法成为高光谱遥感重金属污染监测领域的研究热点。近年来越来越多的特征提取方法被引入到高光谱遥感领域，主要有小波变换[9-10](Wavelet transform，WT)、分形特征提取法[11-12](Fractal dimension，FD)、经验模态分解[13](Empirical mode decomposition，EMD)等。上述方法在一定程度上均可较好地提取出光谱变异信息，但仍存有一定局限性。WT难以选择合适的小波基函数，FD受到信号长度的影响，EMD分解过程中存在模态混叠现象。变分模态分解(Variational mode decomposition，VMD)是近年来提出的一种新的自适应性信号分解方法，该方法解决了EMD方法的问题，并具有很好的噪声鲁棒性[14]，目前被广泛应用于信号分析及故障特征提取等方面[15-16]，而在高光谱遥感领域鲜见报道。

本文引入VMD理论，并结合多尺度熵(Multiscale entropy，MSE)构建一种基于变分模态分解的MSE光谱弱信息提取及定量描述模型(VMD-MSE)。利用VMD对受Cu2+污染的玉米光谱数据进行分解获得模态分量，以MSE方法对其进行定量描述获得VMD-MSE模型值(VM)，探讨VM值与玉米叶片中Cu2+含量的相关关系，并基于VM值构建玉米叶片中Cu2+含量预测模型，以实现玉米重金属污染程度的诊断，为农作物重金属污染监测研究提供参考。

1 理论与方法

1.1 变分模态分解(VMD)

VMD分解过程实质为变分问题的求解过程，主要分为变分约束问题的建立和求解两部分，最终将信号分解为数个有限带宽的固有模态分量[17-18]。该过程包含自适应滤波组，表现了较好的噪声鲁棒性。

(1)变分约束问题的建立

假设每个“模态”是具有中心频率的有限带宽，变分问题描述为寻求K个模态函数uk(t)(k=1，2，…，K)，使得每个模态的估计带宽之和最小，约束条件为各模态之和等于输入信号f。具体步骤如下：①每个模态通过Hilbert变换计算与之相关的解析信号。②对于每个模态，通过加入指数项调整各自估计的中心频率，把模态的频谱变换到基带上。③通过对解调信号进行H1高斯平滑对带宽进行估计。④得到一个变分约束问题

(1)

其中

(u1，u2，…，uK)=u

(ω1，ω2,…，ωK)=ω

式中δ(t)——单位脉冲函数

u——分解得到的K个模态分量

ω——各模态分量的中心频率

f——输入的原始信号

最后求解该问题。

(2)变分约束问题的求解

①引入二次惩罚因子α和拉格朗日算子λ(t)，将约束性变分问题变为非约束性变分问题，构成扩展的拉格朗日表达式为

L({uk},{ωk},λ)=

(2)

式中α——二次惩罚因子，在高斯噪声存在的情况下可以保证信号的重构准确度

λ(t)——拉格朗日算子，用来保持约束条件的严格性

②利用ADMM算法迭代搜索求取上述扩展的拉格朗日函数的鞍点，即式(1)约束变分模型的最优解，其中解得模态分量uk及中心频率ωk分别为

(3)

(4)

1.2 多尺度熵(MSE)

MSE是基于样本熵(Sample entropy, SpEn)的一种时间序列复杂性的度量方法，用于从不同尺度度量时间序列的复杂程度[19]。假设原始数据为X={x1，x2，…，xN}，则MSE的具体计算步骤如下[20]：

(1)给定嵌入维数m、相似容限r及尺度因子τ=1，2，…，τmax，建立新的粗粒序列

(5)

对于每个τ，原序列被分为N/τ个长度为τ的粗粒序列。

(2)计算样本熵

①由步骤(1)所得的粗粒化序列组成m维向量

Y(i)=(yi(τ),yi+1(τ),…，yi+m-1(τ))
(1≤i≤N-m)

(6)

②计算Y(i)、Y(j)间的距离

(7)

(8)

(9)

④将维数增至m+1，重复①～③，得到Bm+1(r)。

⑤此序列的样本熵为

(10)

当N取有限值时，式(10)表示为

(11)

利用式(6)～(11)计算每个尺度序列的SampnEn，即可得到

MSE(X)=SampnEn(y(τ),m,r)

(12)

由以上计算过程可知，MSE由τ、m和r3个参数决定。经过多次实验分析，选取τmax=5，m=2，r=0.15δ(其中δ为X的标准差)时，光谱奇异信息规律性最显著，能够较好地反映光谱信号中的弱信息特征。

1.3 VMD-MSE模型

基于不同浓度Cu2+胁迫梯度下的玉米叶片光谱形态仍极为相似的特点，难以利用其进行有效的污染诊断，提出了一种VMD-MSE光谱特征提取及定量描述模型：将各胁迫梯度下的玉米叶片光谱经VMD分解，选取能够表征光谱奇异特征的模态函数uk，计算其多尺度熵，实现光谱变异特征的定量描述。算法流程如图1所示。

图1 VMD-MSE模型流程Fig.1 Flow chart of VMD-MSE model

2 实验数据及模型应用

2.1 实验数据获取

(1)植株培养：选用有底漏的盆钵进行“中糯1号”玉米种子培育。设置胁迫浓度为0(ck)、100、200、300、400、700、800、900 μg/g的CuSO4·5H2O溶液，将其翻土加入到玉米实验盆钵中，每梯度设置3组平行实验，共24盆。2016年5月6日对玉米种子进行催芽处理，出苗后向盆栽中添加NH4NO3、KH2PO4和KNO3营养液。玉米培育期间定期通风与浇水，保持适宜的培育温度与湿度。

(2)光谱数据采集：2016年7月17日对玉米叶片反射光谱进行测量。在50 W卤素灯光源照射下，将光谱范围为350～2 500 nm的SVC HR-1024I型地物光谱仪探头视场角设置为4°并垂直于叶片进行光谱采集，采集的光谱使用平面白板进行标准化。选取每盆玉米植株的老、中、新3种代表性叶片进行光谱测试，每盆获得3组数据。各胁迫梯度光谱由3组平行实验的9条光谱数据求均值所得。各胁迫梯度下的玉米叶片光谱曲线如图2所示。

图2 不同浓度Cu2+胁迫下玉米叶片光谱数据Fig.2 Spectral data of corn leaves stressed by different Cu2+ concentrations

(3)Cu2+含量测定：2016年9月16日对采集过光谱数据的叶片进行冲洗、干燥、粉碎等样品预处理，再经高纯硝酸和高氯酸消化处理后用WFX-120型原子吸收分光光度计对叶片中Cu2+含量进行测定，每梯度测量3次后取平均值作为该梯度叶片中的Cu2+含量，测量结果如表1所示。随着胁迫梯度的增大，叶片中Cu2+含量的总体趋势表现为先增大后减小，原因可能是随着土壤中Cu2+增加，玉米对其吸收逐渐增大；而当土壤中Cu2+浓度过高时，对植被根系产生毒害作用，造成根部受损，对Cu2+吸收逐渐降低。本文选取胁迫梯度为0、100、200、300、400 μg/g的5组Cu2+含量数据用于建立回归模型，胁迫梯度为700、800、900 μg/g的3组数据用于模型验证。

2.2 VMD-MSE光谱弱信息探测

以0、100 μg/g胁迫浓度下的玉米叶片光谱数据为例，分析VMD-MSE模型的特征提取与描述

表1 不同胁迫浓度下的玉米叶片中Cu2+含量Tab.1 Cu2+ contents in corn leaves stressed by different copper concentrations μg/g

过程：

(1)对叶片光谱进行分解模态数K为2～7的VMD处理，获得第K模态函数uk，如图3所示。由植被反射光谱特性曲线可知，在960 nm和1 900 nm附近，植被光谱受水分和大气的影响，出现较大干扰信息。由图3可知，光谱信号噪声随着分解次数的增大而增大。K=3时，960 nm和1900 nm处噪声较弱，且全波段奇异信息显著；从K=4开始，1 000 nm和1 900 nm处噪声开始出现，且奇异信息逐渐弱化。因此选择K=3既实现模态分离，又有效抑制噪声，可有效提取玉米铜污染胁迫信息。

图3 VMD分解结果Fig.3 Results of decomposition of VMD

(2)计算由VMD分解得到的第3模态函数u3的多尺度熵，即为VMD-MSE的模型值(记为VM)，不同尺度下的模型值记为VMi(i=1，2，…，τmax)。

3 结果与分析

3.1 VMD-MSE污染监测

对0、100、200、300、400 μg/g胁迫浓度下的玉米叶片光谱应用VMD-MSE模型进行污染信息提取与测度。各胁迫浓度下的玉米叶片光谱的VM值如图4所示。不同胁迫浓度下，VM曲线变化趋势有所差异，但总体而言，曲线变化趋势相似，即随着τ的增大而增大。在所有尺度因子上，VM值与作物受污染程度呈负相关关系，即受污染越严重，VM值越小。显然，各胁迫浓度下的光谱细节特征通过不同尺度VM(即VMi)得到体现，能更好地揭示光谱曲线畸变特征。

图4 不同胁迫浓度下玉米叶片光谱的VM值Fig.4 VM values of corn leaf spectrum under different stress concentrations

图5 不同胁迫梯度下各尺度VM值与玉米叶片中Cu2+含量相关关系Fig.5 Relationships between VM on five scales and Cu2+ contents in corn leaves stressed by different Cu2+ concentrations

不同胁迫浓度下，VMi与叶片中Cu2+含量的相关变化如图5所示，两者呈现相反的变化趋势。随着胁迫浓度的增加，叶片中Cu2+含量先增加后减小，相应地VMi表现为先减小后增大的趋势，两者相关系数绝对值均大于0.90，其中VM1与Cu2+相关性最强，相关系数为-0.968。据此可认为，各尺度VM均可在一定程度上诊断玉米受铜胁迫程度，其中VM1效果最优。

3.2 VM值与玉米铜含量的回归分析

为了进一步确定VMD-MSE模型值与玉米铜胁迫的定量关系，以不同胁迫浓度下玉米叶片中Cu2+含量为因变量y，以VMi为自变量x进行回归分析，VMi值与铜含量的线性回归方程、决定系数(R2)和均方根误差(RMSE)如表2所示。VM1～VM5的回归决定系数均大于0.88，且显著性水平P<0.05。整体上，从VM1到VM4，R2逐渐减小，RMSE逐渐增大，这可能是由于随时间序列复杂度的增加，样本熵估计准确率下降所致。

为了验证上述模型的可靠性与准确度，选用700、800、900 μg/g 3个胁迫浓度下的样本数据对模型进行检验，图6描述了玉米叶片中Cu2+含量实测值与模型预测值的相关变化。各预测模型的Cu2+含量预测值与实测数据具有很好一致性。VM1和VM2回归模型的预测结果较优，相关系数大于0.70，其中VM1模型预测结果最优，相关系数为0.992；VM3～VM5回归模型预测效果较差，相关系数小于0.650，可能是由于随着污染程度的增加，光谱信号的奇异信息随机性增加，且随度量尺度的增大，熵估计准确率下降，使VM3～VM5模型稳定性及适用性变差。因此，通过VMD-MSE模型值VM1可以建立玉米Cu污染程度的预测模型，从而实现监测玉米Cu污染状况的目的。

表2 基于不同尺度VM值的玉米叶片中Cu2+含量回归模型Tab.2 Regression models for VM on different scales with Cu2+ content in corn leaves

图6 玉米叶片中Cu2+含量实测值与预测值的相关关系Fig.6 Relationships between measured and predicted values of Cu2+ contents in corn leaves

4 结论

(1)将VMD运用到高光谱信息提取中，从相似光谱信号中挖掘有利于诊断重金属Cu污染的畸变信息。结果表明，对原始光谱数据进行3次VMD分解后，可有效抑制噪声，实现模态分离，从而提取出隐藏在光谱噪声中的污染弱信息。

(2)将MSE方法应用于VMD分解结果，揭示了污染信息不同尺度的内在差异。

(3)基于VMD和MSE构建的VMD-MSE模型可有效诊断作物铜污染程度，模型值VM与玉米叶片中Cu2+含量呈现显著负相关，两者回归精度达到0.88以上，其中VM1模型预测效果最优。本研究构建的预测模型能在一定程度上预测玉米受污染程度，但具有一定局限性，即预测精度随污染程度增大而降低，这种情况可通过增加样本数量对模型进行优化而得到改善。