基于SD-SVD-Burg的玉米叶片铜铅污染甄别与程度诊断

韩倩倩，杨可明，高伟，李艳茹，张建红

（中国矿业大学（北京）地球科学与测绘工程学院，北京 100083）

近年来，由于人口急剧增长、工业迅猛发展，大量未经处理的工业固体废弃物（主要来自采掘业、有色金属冶炼、化学原料及化学制品等）在土壤表面堆放和倾倒，有害废水向土壤中渗透，农药的不科学使用以及大气中的有害气体及飘尘随雨水降落到土壤等，导致土壤重金属铜（Cu）、铅（Pb）、汞（Hg）、镉（Cd）等污染问题日趋严重[1-2]。重金属在土壤中不能被降解，其超限后会污染农作物并严重影响农产品品质，且可以通过食物链进入人体，破坏人体正常的生理功能，进而对人体造成危害[3-5]。因此，如何快速准确鉴别作物受重金属污染的元素类别以及监测其所受污染程度已成为目前的研究热点。常规的生化检测方法主要是通过采集作物样本来检测其中重金属元素含量、浓度及一些生化特征参数等，但这种方法存在费时费力、工作量大等缺点，且无法满足无损、快速、便捷等污染监测需求。高光谱遥感具有光谱信息丰富且光谱连续、波段多且范围窄、从可见光到近红外有数十甚至数百个波段、理化特性与参量反演能力强等特点，因此已成为现今环境、农作物等重金属污染监测的热点研究内容之一[6-7]。

重金属在作物叶片内积累会导致叶细胞畸形，使细胞质膜的选择透性、组成和结构均发生改变，而细胞结构改变、成分分解会使光在叶片内传输时出现不连续的折射，从而引起整个叶片到冠层的光谱变化。研究发现，重金属胁迫下作物的电磁波反射光谱变异能反映受重金属污染后植物的色素、细胞结构和含水量等生化参数的特征信息，并能响应出重金属污染元素类别[8]。因此，研究作物光谱变化特征和污染信息提取的光谱技术至关重要。基于AR 模型的功率谱估计是现代功率谱估计中最常用的方法，而利用AR模型时必须计算出AR 模型的参数，其中Burg 算法被广泛应用。近年来，Burg 算法被应用在很多领域，如罗治军等[9]将Burg算法运用到经验小波变换（Empirical wavelet transform，EWT）中，使得EWT 算法能在功率谱上找到明确的模态边界；王路晗等[10]用Burg算法对风切边雷达信号进行处理；李佳等[11]利用Burg算法获得车流量的雷达信号信息，对双波束车流量的速度进行检测；钱莹晶等[12]采用Burg算法较为准确地估计出脉冲编码雷达中目标物体的多普勒频偏；Iman等[13]将Burg 算法用于测量患者心率等；但Burg 算法在高光谱信号处理与分析方面鲜有报道。本文尝试将Burg算法用于玉米叶片的高光谱数据处理，以甄别玉米受重金属污染的元素类别和所受污染程度。

本文将以不同浓度Cu 和Pb 污染胁迫下的盆栽玉米及其植株叶片光谱作为研究对象，采用Burg 算法绘制玉米叶片光谱的功率谱曲线。考虑到该算法适应于短序列，实验将首先选取玉米的紫谷、绿峰和红边3 个特征区间，并对特征区间的高光谱数据进行一阶光谱微分（Spectral derivative，SD）和奇异值分解（Singular value decomposition，SVD）处理，以消除噪声等。然后结合Burg 算法对预处理后的光谱信号进行功率谱估计，通过曲线特征对受重金属Cu 和Pb 胁迫生长的玉米叶片进行污染元素种类辨别和污染程度诊断，同时基于不同年份同等条件的Cu、Pb污染玉米叶片光谱数据检验该模型的稳定性与可行性。

1 材料与方法

1.1 实验设计

1.1.1 玉米植株培养

实验选用有底漏的花盆培植“密糯8号”玉米种子，重金属 Cu 和 Pb 以不同浓度 CuSO4·5H2O 和 Pb（NO3）2溶液的形式翻土加入到花盆中，待玉米出苗后向盆中喷洒适量且等量的NH4NO3、KH2PO4和KNO3营养液。实验设计Cu、Pb 污染浓度为0、50、200 mg·kg-1和400 mg·kg-1，记为 CK（0）、Cu（50）、Pb（50）、Cu（200）、Pb（200）、Cu（400）、Pb（400），每个浓度设置3 组平行试验，玉米培育期间保持各盆栽生长条件一致。

1.1.2 玉米叶片光谱测量

使用SVC HR-1024I 型地物光谱仪采集玉米叶片反射光谱。采集时选择50 W 卤素灯作为光源，探头视场角设置为4°，并垂直于玉米叶片表面5 cm 处，采集输出的光谱反射率并用平面板对光谱进行标准化。测量时每株采摘3 片代表性叶片[玉米顶部刚长出不久的新叶片（New，N），底部完全成熟的叶片（Old，O）和位于中部的叶片（Middle，M）]分别进行光谱测量，测量过程中将玉米叶片放置于黑色塑料袋包裹的硬纸片上以排除其他因素的干扰。每盆获得3组数据，去除异常光谱值后取平均值作为研究数据。

1.1.3 Cu和Pb含量测定

收集采集光谱数据后的玉米叶片，对其进行冲洗（去除灰尘杂质）、干燥（烘干水分），用高纯硝酸、高氯酸消化（排除其他离子的干扰并使Cu、Pb以离子的形式存在）等预处理，Cu 和Pb 含量的测定依据《硅酸岩石化学分析方法》（GB/T 14506.30—2010）。叶片中Cu 含量的测定采用电感耦合等离子体发射光谱仪进行定量分析，叶片中Pb 含量的测定采用WFX-120 原子吸收分光光度计[14]。对不同浓度Cu、Pb 处理下玉米叶片中的老（O）、中（M）、新（N）3种叶片分别测量3次后取平均值作为该浓度下叶片中Cu、Pb 的含量。实验所测得的玉米叶片中Cu和Pb的含量如表1所示。

表1 不同胁迫浓度下玉米叶片中Cu和Pb含量（n=3）Table 1 Content of Cu and Pb in maize leaves stressed by different stress gradient（n=3）

1.2 一阶光谱微分

光谱微分（SD）变换是一种常用的光谱特征分析方法，光谱微分变换可在一定程度上限制低频噪声对目标光谱的影响，并能有效削弱大气效应等的干扰[15]，进而增强光谱曲线的细微变化。光谱一阶微分的计算公式为：

式中：λi-1、λi和λi+1为玉米叶片光谱数据相邻波段的波长；ρ′（λi）是波长λi的一阶微分值；ρ（λi+1）和ρ（λi-1）为波长λi+1和λi-1对应的玉米叶片光谱反射率。

1.3 奇异值分解

奇异值分解（SVD）是线性代数中一种矩阵分解技术，在信号处理方面有重要应用。SVD分解是谱分析理论在任意矩阵上的推广[16]。

假设原矩阵为Am×n（m＜n），经过奇异值分解后为：

式中：Um×m和Vn×n为单位正交矩阵，其对应的列向量和行向量为单位向量且互相正交，都是矩阵A对应的特征向量；∑m×n为对角矩阵，∑m×n=diag（σ1，σ2，σ3，…，σm），其元素σi为按降序排列矩阵A的奇异值，它的对角元素即为Am×n的奇异值，∑m×n中除对角线外，其余元素均为0。相关研究表明，矩阵∑m×n中的奇异值即为原始数据集Am×n中的全部重要特征，其他特征大多为多余的或是噪声。

就光谱信号的SVD 处理而言，设实验所测得的玉米叶片光谱数据矩阵为矩阵A，则矩阵A中存在m×m的正交矩阵U=[u1，u2，…，um]和n×n的正交矩阵V=[v1，v2，…，vn]，则：

式中：σi为玉米光谱数据A的奇异值；ui为相应于奇异值σi的左奇异向量；vi为相应于奇异值σi的右奇异向量，且满足：

式中：i=1，2，…，m；λi是方阵AAT的特征值；λi′是方阵ATA的特征值。根据式（4）和式（5）可求得矩阵U和V。

根据式（6）、式（7）、式（8）得到玉米光谱数据A奇异值σi为：

1.4 基于AR模型的Burg算法

AR模型是现代谱估计中最常用的一种功率谱估计模型，而AR 模型在功率谱估计时，必须计算出模型的参数。AR 模型参数求解以线性预测理论为基础，发展出许多算法，主要有Yule⁃Walker 法、协方差法、Burg 法和改进的协方差法，Burg 算法以其良好的分辨率和稳定性得到广泛应用[11]。

Burg 算法为一种递推算法，不需要计算序列的自相关，是直接通过原始序列求出模型参数，进而求得功率谱。其基本思想是根据线性预测器的前向和后向预测功率的总均方误差之和最小的准则直接从观测数据来估计反射系数，然后利用Lenvinson-Durbin 算法的递推公式求出AR 模型参数[17]。其具体算法如下：

（1）定义前向和后向误差分别为ef（n）和eb（n），预测误差功率为E。初始化前向、后向误差及预测误差功率：

式中：x（n）为玉米叶片光谱信号序列，n=0，1，…，N-1。

（2）计算反射系数：

式中：m=1，2，…，p，p为模型的阶数表示m-1阶前向预测误差为m-1阶后向预测误差。

（3）根据Lenvinson 递推算法计算模型参数及预测误差功率：

式中：k=1，2，…，m-1；am（k）表示第m阶模型的第k个系数。

（4）递推高一阶前、后向预测误差：

把m更新为m+1，重复步骤（2）至步骤（4），直到预测误差功率Em满足要求。

（5）得到参数后，通过公式（18）计算玉米叶片光谱信号的功率谱密度为：

式中：P为功率谱密度；ak为 AR 模型参数；σ2为序列的方差。

2 结果与讨论

2.1 SD-SVD-Burg算法的应用与讨论

实验所获得的Cu 和Pb 胁迫下玉米叶片光谱如图1 所示。从曲线看出，玉米叶片的反射光谱曲线整体呈现先上升后下降的趋势，在波长700～2 500 nm间健康叶片的反射波峰强度高于受重金属元素胁迫的叶片光谱，但受Cu、Pb 胁迫叶片的光谱曲线差异较小，仅从光谱曲线无法辨别玉米受重金属元素的污染类别与污染程度。根据植物光谱的光谱特性，将Cu和Pb胁迫下玉米叶片光谱划分为紫谷、绿峰和红边3个特征区间作为研究对象，各特征区间对应的范围如表2所示。波长380～450 nm出现一个吸收谷即紫谷，主要由栅栏组织中的光合色素强烈吸收所引起，由于相对地反射绿光，因此健康植物的视觉效果是绿色的。波长500～600 nm 出现一个小反射峰即绿峰，是描述色素状态和健康状况的重要指示波段。波长670～780 nm 形成陡坡即红边，是植物反射光谱中最具特征的指示波段。重金属会对植物的色素、细胞结构和含水量等生化参数产生影响，因此选择特征波段来反映健康植物和受重金属污染的植物是可行的。

表2 各光谱特征区间名称及范围Table 2 Names and ranges of each spectral characteristic interval

实验分别对Cu 与Pb 胁迫梯度下各特征区间的光谱经式（1）处理，得到一阶微分光谱，然后分别对各特征区间的微分光谱进行SVD 处理，获得各特征区间的光谱奇异值。因Burg 算法适应于短序列信号，经光谱一阶微分以及SVD 处理后数据在去除噪声的同时有效压缩了数据。然后利用Burg 算法绘制功率谱曲线，首先输入序列x（n），即各特征区间经预处理后的玉米叶片光谱以矩阵形式作为一个信号序列，接着设置模型阶数，阶数是Burg 算法的重要参数，使用不同的阶数对算法的结果影响也不同。因此，实验分别选择阶数p为 3、4、5 进行仿真，通过 Lenvinson 递推算法求得模型参数，利用式（18）求得功率谱密度。最终健康玉米叶片及Cu、Pb 胁迫下紫谷、绿峰、红边的信号处理结果如图2～图5所示。

分析图2和图3可知，在紫谷区间，健康玉米叶片的功率谱曲线在阶数p=3 时，曲线较平缓，无波峰波谷，而受Cu、Pb 污染的玉米叶片均出现1 个波峰，由此可定性分析玉米叶片是否已受Cu 或Pb 污染。从图3看出，在p=5时，受Cu污染的曲线表现为单波峰，而受Pb 污染的曲线表现为双波峰，因此通过此特征可快速定性分析玉米所受重金属污染的元素类别（Cu 或 Pb）。随着 Cu 和 Pb 胁迫浓度的增加，受 Cu 或Pb 污染的曲线波峰均出现左移现象，因此可通过曲线波峰的位置来判断玉米叶片受重金属离子的胁迫程度。在p=4 时，受Cu 污染的曲线波峰较陡且随着胁迫浓度的增加，波峰逐渐有变缓的趋势，而受Pb污染的曲线波峰较缓，随着胁迫浓度的增加波峰逐渐变缓，因此在p=4 时也可通过波峰坡度来判断玉米叶片所受重金属污染的类别及污染程度。

分析图2和图4可知，在绿峰区间，健康玉米叶片在p=3 时，曲线表现为单波峰，而在Cu、Pb 胁迫下曲线均出现3 个波峰，因此可通过波峰个数来判断玉米叶片是否受到重金属元素污染。从图4 可知，当p=5时，受Cu污染的曲线出现4个波峰，受Pb污染的曲线为3 个波峰，由此可辨别玉米所受重金属污染的元素类别（Cu 或 Pb）。当p=4 或 5 时，观察到 Cu 或 Pb 胁迫浓度越大，波峰越陡，因此可根据波峰的坡度来粗略诊断玉米受重金属离子的污染程度。

分析图2和图5可知，在红边区间，健康玉米叶片在p=5 时，仅出现一个波峰，而叶片受Cu、Pb 胁迫下的曲线有2～3个波峰，可由此来诊断玉米叶片是否受到重金属元素的污染。从图5 看出，受Cu 与受Pb 胁迫下的功率谱曲线未表现出明显差异，无法通过曲线特征来甄别玉米所受重金属元素的污染类别；随着重金属污染程度的增加，曲线的波峰有变缓的趋势，因此可根据波峰的坡度来判断玉米所受重金属离子的污染程度。

由图2～图5可知，在紫谷和绿峰区间通过功率谱密度曲线图，均能较好地定性分析玉米是否受Cu、Pb污染并能判断出受哪种元素污染，利用曲线的波峰变化能诊断玉米所受重金属离子的污染程度；在红边区间，可辨别出玉米是否受到重金属元素污染以及所受污染程度，但未能诊断出污染元素的类别。

为进一步从定量的角度去检测紫谷、绿峰和红边区间的功率谱曲线对重金属Cu 与Pb 胁迫程度的监测效果，分别统计了在阶数p=5 时，不同胁迫浓度下功率谱曲线的平均功率与玉米叶片中Cu、Pb 含量的相关系数，结果如表3 所示。由表可知，Cu 胁迫下紫谷、绿峰及红边3 个区间的功率谱曲线平均功率与玉米叶片中Cu含量的相关系数均在0.9以上，表现出较强的相关性，说明对重金属Cu 胁迫程度的监测达到很好的效果；Pb 胁迫下在绿峰区间功率谱密度曲线平均功率与玉米叶片中Pb 含量的相关系数达到0.995 8，效果理想。

2.2 算法稳定性验证与讨论

为验证SD-SVD-Burg 算法在监测重金属污染信息方面的稳定性和普适性，选取不同年份同等实验条件的Cu、Pb 污染玉米叶片光谱数据作为验证组对该模型进行可行性检验。根据表3，选取与Cu 和Pb 相关性最高的绿峰区间作为研究对象，验证组统计的Cu 与Pb 胁迫下绿峰功率谱曲线平均功率与玉米叶片中Cu 与Pb 含量的相关系数如表4 所示。由表4 可知，在绿峰区间，功率谱曲线的平均功率与玉米叶片中Cu 含量相关系数为0.921 3，与Pb 含量的相关系数达0.991 5，均表现出较强的相关性，因此该算法在区分作物重金属污染类别和监测污染程度方面具有一定的有效性和可行性。

3 结论

表3 Cu和Pb胁迫下紫谷、绿峰、红边功率谱曲线平均功率与玉米叶片中Cu和Pb含量的相关系数Table 3 The correlation coefficient between the average power of purple valley，green peak and red edge power spectrum curve under Cu and Pb stress and the content of Cu and Pb in maize leaves

（1）在紫谷与绿峰区间，健康玉米叶片与受Cu、Pb 污染的叶片波峰个数均不相同，因此可根据曲线的波峰个数定性分析玉米是否受到污染，并能判断出玉米所受重金属污染的元素类别（Cu或Pb）。利用曲线的波峰位置及波峰坡度变换可定性诊断玉米所受重金属元素的污染程度。

（2）通过统计不同胁迫浓度下功率谱曲线的平均功率与玉米叶片中Cu、Pb含量的相关系数，得出绿峰曲线的平均功率可较好地定量分析Cu、Pb 重金属污染程度。利用不同年份的玉米叶片光谱数据对该算法进行稳定性验证，表明SD-SVD-Burg 算法可以作为监测玉米是否受Cu、Pb污染的一种参考方法。