1

引言

重金属对农作物的污染主要是其所生长土壤中含有能被植株吸收的铜（Cu）、铅（Pb）、镉（Cd）等重金属元素，而矿产资源开采、废弃物排放、含有重金属元素的原料使用甚至交通运输中突发事件等都能使重金属元素渗入土壤。作物可通过根系从土壤中吸收大量的重金属元素并在作物体内富集，进而污染食物链，人体一旦摄取过多的重金属就会出现贫血、记忆力减退、免疫力下降等症状。传统的作物体内重金属离子含量测定的工序和仪器繁多，需要大量、细致的人工操作才能精确实现。因此，探索一种快速、可靠的作物重金属污染监测手段是目前研究热点，而可见光与近红外电磁波的高光谱遥感监测方法则是重要研究内容。

重金属胁迫下作物的电磁波反射光谱变异能响应重金属污染特征信息，甚至能反映重金属污染元素类别。玉米叶片是玉米所有器官中富集Cu和Pb重金属总量最高的器官，虽然受重金属污染的玉米叶片细胞结构会出现较大变化，且理化反应较为复杂，但是污染叶片光谱与健康的叶片光谱相比，其变异性仍较小。因此，研究光谱微弱变化特征和污染信息提取的光谱分析技术尤为重要。而高阶谱估计结果含有相位信息，且可以抑制白噪声，适用于非线性和非高斯系统的描述，因而其估计结果更为接近实际。现今，高阶谱估计算法常用于地震勘探、机械故障诊断、天体光变周期分析、信号处理与目标识别等众多领域，但尚未查到高阶谱用于高光谱信号处理与分析的相关文献。本文以受重金属铜离子（Cu2+）和铅离子（Pb2+）胁迫生长的玉米叶片光谱数据为研究对象，将高阶谱中的双谱估计方法应用于玉米叶片微分光谱处理，以双谱三维图方式实现Cu2+和Pb2+污染定性分析与污染元素的种类辨别。

2

数据获取

玉米叶片光谱数据采集于分别受CuSO4·5H2O和Pb(NO3)2胁迫下的盆栽“中糯1号”玉米植株，并采用其中0、250、500、750、1000μg·g-15个质量比梯度胁迫的玉米叶片光谱数据作为研究对象，为了模拟真实环境下的土壤重金属污染，计算各重金属胁迫组所需的CuSO4·5H2O、Pb(NO3)2试剂剂量后，按花盆标注添加到对应花盆中并搅拌均匀，每个胁迫梯度各设置3个平行组。实验过程中，Cu2+、Pb2+胁迫的750μg·g-1和1000μg·g-15梯度下的盆栽玉米叶片在幼苗期开始渐渐枯萎，无法用于后续数据分析，故最后数据来源于0、250、500μg·g-13个胁迫梯度，Cu2+和Pb2+的3个胁迫梯度分别标注为CK(0)、Cu(250)、Cu(500)和CK(0)、Pb(250)、Pb(500)。

光谱数据采集时使用波段为350～2500nm的美国SVC HR-1024I型高性能地物光谱仪，采集过程安排在室内进行，采集到的光谱反射系数用专用的平面白板进行标准化处理。另外，每棵玉米植株分别测量了其老（Old，O）、中（Middle，M）、新（New，N）3种叶片光谱。每片光谱各测量3次，去除异常光谱值后取平均作为最终结果。图 1显示了盆栽温室的位置、玉米的生长过程及老（O）、中（M）、新（N）叶片的位置。

图1 盆栽温室的位置和玉米的生长周期

玉米叶片光谱数据采集结束后，对各叶片样品作冲洗、干燥、粉碎等预处理，再进行高纯硝酸、高氯酸消化处理，然后用WFX-120型原子吸收分光光度计测定叶片样品中Cu2+含量，用PerkinElmer，ElanDRC-e型等离子体质谱分析仪测定Pb2+含量。

3

理论与方法

一个随机过程可用累积量从时域来描述其统计特性，也可用高阶谱（包括功率谱）从频域来描述其统计特性。实际中常用有限长度数据估计一随机过程的高阶累积量谱，一般有BT、周期图等常规法与AR、MA和ARMA等模型参数法。常规法直接且容易实现，但其应用能力受到估计的统计方差和频率分辨率限制;而模型参数法谱估计和信号模型联系紧密，它不限定信号模型的最小相位系统，且频率分辨率性能远优于常规法。在谱分析时，较多应用ARMA模型参数法的高阶谱估计，主要是ARMA模型与谱分析的发展密切，也是继最大熵谱后讨论最多的谱分析。随机过程（如信号等）高阶累积量的多维傅里叶变换称为高阶谱估计。其中，三阶谱（也称双谱）定义为三阶自相关函数的二重傅里叶变换，故三阶谱的高阶谱估计即为双谱估计。本文采用ARMA模型参数法进行光谱数据序列的双谱估计。

高阶谱能抑制高斯噪声且分辨能力强，尤其适用于盲信号、非线性或非高斯信号的分析和处理，信号经过处理后能得到相位、能量和非线性等有用信息。高阶谱研究中的“热点”是三阶谱（双谱），因其阶数最低，容易实现，含有功率谱中所没有的相位信息。三阶谱（双谱）估计能够获取高阶谱的所有特征，因此应用较为广泛。

双谱估计有利于提取光谱间变异的弱信息，故可利用双谱估计对光谱信号进行处理。双谱估计的ARMA模型参数法是将已知有限长非正态分布的数据序列，通过三阶累积量估计序列确ARMA模型参数，然后进行双谱估计。设光谱信号为x（t），则3阶累积量为

利用式（1）、（2）可得双谱估计的bisp_rts矩阵和bisp_qs矩阵。首先，基于式（1）和式（2）对高阶谱工具箱中的amarts和amaqs函数进行编译;接着依据编译的amarts和amaqs函数处理信号，以获取各函数相应的ARMA模型参数;然后利用这2种ARMA模型参数构建复数型的bisp_rts和bisp_qs矩阵;最后，采用实数化后的bisp_rts和bisp_qs矩阵分别绘制各矩阵的双谱平面图和双谱三维图。

4

污染元素识别

通过原始光谱数据曲线的对比分析，不同胁迫梯度的盆栽玉米叶片光谱间相似度极高，难以用传统的光谱分析方法识别玉米叶片重金属Cu或Pb污染元素种类以及诊断污染程度的有效信息。针对这一问题，可利用高阶谱估计的ARMA模型参数法，估计得到双谱平面图和双谱三维图;然后基于双谱三维图进行叶片中Cu2+和Pb2+污染定性分析以及污染元素的种类识别。

（1）叶片光谱数据的一阶微分处理。对不同胁迫梯度下所有的玉米叶片光谱作一阶微分处理。因叶片光谱数据的波谱分辨率高，波段间隔小，所以经微分处理后的数据波动性变大，扩大了光谱间变异的差异，从而有利于高阶谱估计和提取重金属污染的有效信息。

（2）双谱三维图的绘制。将不同胁迫梯度下玉米叶片的微分光谱数据依次进行ARMA模型参数法估计，通过多次的估计结果对比分析，确定估计的参数最佳值。由实验可得，当把ARMA模型中amarts和amaqs函数中变量设置为:p、q分别设置为2、1，累积量阶数为3，maxlag为10，samp_seg为128，overlap为默认值0，flag为默认为biased（有偏）后，可获取识别Cu2+和Pb2+污染元素种类的最佳模型参数;然后运用Matlab中的mesh函数，绘制实数化bisp_rts和bisp_qs矩阵的双谱三维图。

（1）基于bisp_qs矩阵双谱三维图的叶片Cu2+和Pb2+污染定性分析。依据实数化bisp_qs矩阵，绘制不同胁迫梯度下玉米老（O）、中（M）、新（N）3种叶片微分光谱的双谱三维图如图 2、图 3所示。图 2为健康玉米叶片微分光谱的双谱三维图，图 3为不同Pb2+和Cu2+胁迫梯度下玉米叶片微分光谱的双谱三维图。

由图 2、图 3可看出，图 2中健康叶片微分光谱双谱三维图中间会有明显竖起的“尖塔”;而图 3中受Cu2+和Pb2+污染的叶片微分光谱双谱三维图中间的“尖塔”都会全部消失;因此得出，基于bisp_qs矩阵的双谱三维图中有无“尖塔”，可定性分析玉米叶片是否已受Cu2+和Pb2+的重金属污染。

图2 健康玉米老、中、新叶片微分光谱bisp_qs矩阵的双谱三维图

图3 不同Pb2+和Cu2+胁迫梯度下老、中、新玉米叶片微分光谱bisp_qs矩阵的双谱三维图

（2）基于bisp_rts矩阵双谱三维图的叶片Cu2+或Pb2+污染元素种类识别。依据实数化bisp_rts矩阵，绘制不同胁迫梯度下玉米老（O）、中（M）、新（N）3种叶片微分光谱的双谱三维图如图 4所示。

由图 4可看出，受Pb2+污染的玉米O、M叶片微分光谱双谱三维图外圈会出现红色“圆顶”形状，中间会有蓝色的低谷;受Cu2+污染的玉米O、M叶片微分光谱双谱三维图外圈会出现蓝色“刀锋”形状，中间会有蓝色“尖塔”。因此得出，基于bisp_rts矩阵的双谱三维图中出现的“圆顶”或“刀锋、尖塔”形状，能够直观可视地辨别出Cu2+或Pb2+污染元素的种类。

但根据不同胁迫梯度下玉米N叶片微分光谱的双谱三维图，不能辨别出Cu2+或Pb2+污染元素的种类，可能与新叶片的生化特征或对重金属元素的富集能力有关，其中的原因有待进一步研究。

图4 不同Pb2+和Cu2+胁迫梯度下老、中、新玉米叶片微分光谱bisp_rts矩阵的双谱三维

5

结束语

基于不同浓度重金属Cu2+、Pb2+胁迫的污染玉米叶片微分光谱数据，利用高阶谱估计的ARMA模型参数法，对各胁迫梯度下污染的老（O）、中（M）、新（N）玉米叶片微分光谱数据序列进行双谱估计，根据构建估计结果bisp_rts或bisp_qs矩阵并绘制相应的双谱三维图，提出了一种直观可视的玉米叶片Cu2+、Pb2+污染定性分析和污染元素种类识别的快速判别方法，所有玉米叶片的污染定性分析率达到100%，O和M玉米叶片的污染元素种类识别率达到100%。