基于CWT-HHT的玉米叶片铜离子污染信息探测

2023-11-09郭辉石海张全旺

农业环境科学学报 2023年10期

郭辉，石海*，张全旺

（1.安徽理工大学空间信息与测绘工程学院，安徽淮南 232001；2.安徽理工大学矿山环境与灾害协同监测煤炭行业工程研究中心，安徽淮南 232001）

近年来，由于矿山开采、污水灌溉、大气沉降等因素导致农田土壤中重金属含量普遍升高[1-2]。重金属污染下农作物的产量和质量均受到严重影响，并且积累在农作物中的重金属一旦进入食物链，将最终危害人类健康[3-4]。因此，快速、有效地监测农作物重金属污染对于农业生产、粮食安全具有重要意义。传统的光学检测法、电化学检测法和生物检测法需要对样品进行破环性采集，检测过程繁琐且容易造成二次污染，不仅费时费力且难以满足快速、实时、大面积检测的需求[5]。

高光谱遥感技术具有光谱连续精细、波段数多、光谱分辨率高等特点，为农作物重金属污染监测提供了技术保障[6]。目前，相关学者通常采用简单的光谱变换或提取植被指数的方法，开展农作物重金属污染监测研究。刘来等[7]利用光谱微分技术筛选镉污染下油菜叶片光谱敏感波段；Yu 等[8]基于光谱一、二阶微分技术选取重金属污染下冬小麦光谱特征；程凤等[9]构建一种植被指数探测玉米重金属铜污染程度；Martinez等[10]将锂胁迫下农作物光谱反射率提取的植被指数作为预测因子。以上研究取得了一定成果，但在重金属污染监测方面缺乏针对性，且易受农作物中色素、水分、有机质等生化参数的影响[11]。近年来，希尔伯特-黄变换（Hilbert-Huang transform，HHT）时频分析技术被广泛用于非线性非平稳信号的处理与分析，如张世文等[12]将HHT 算法运用于探地雷达数据分析，有效探测出覆土层厚度；王英乐等[13]利用HHT算法提取微震信号的特征参数，实现微震监测与预警；Zao 等[14]采用HHT 算法分析噪声语音信号，有效估计噪声语音的基本频率；Yin 等[15]运用HHT 算法处理电磁辐射信号，提高了煤层破坏前特征识别程度。以上研究表明HHT 算法在非线性非平稳信号的特征提取上有一定优势，但HHT 算法在高光谱数据处理与分析的应用较少。因此，本研究引入HHT 算法用于玉米叶片高光谱数据分析，以探测玉米重金属铜离子污染信息。

考虑到HHT 算法容易受到信号噪声的影响[16]，本研究结合连续小波变换（Continuous wavelet transform，CWT）与HHT 算法的优势，首先对玉米叶片原始光谱进行CWT 处理去除信号噪声，然后进行HHT时频分析，通过研究不同浓度Cu2+胁迫下玉米叶片光谱信号的瞬时能量特征进行污染信息探测。同时与常见的红边位置（REP）、红边归一化指数（NDVI705）和红边植被胁迫指数（RVSI）等植被指数监测农作物重金属污染方法作对比分析，进一步验证CWT-HHT 方法在探测玉米叶片Cu2+污染信息方面的可行性与有效性，以期为农作物重金属污染信息探测提供一种新的方法，并为农作物重金属含量反演奠定基础。

1 材料与方法

1.1 实验设计

1.1.1 玉米植株培养

以“中糯1 号”品种玉米作为研究对象，优选籽粒饱满、大小均匀的种子进行重金属Cu2+胁迫土培实验。实验采用CuSO4·5H2O分析纯溶液胁迫玉米发育生长，设置了4 种胁迫梯度，分别为0、100、300、500µg·g-1，分别标记为CK0、Cu100、Cu300、Cu500。每个浓度梯度设置3 组平行实验，共计12 盆玉米盆栽。2019年5月对玉米种子进行催芽处理，出苗后定期浇灌NH4NO3、KH2PO4和KNO3营养液。所有盆栽玉米的培育条件除了土壤中的Cu2+浓度不同，均在同一条件下进行。

1.1.2 数据采集

玉米出苗后使用美国SVC 公司生产的SVC HR-1024I 高性能地物光谱仪采集玉米叶片光谱，光谱范围为350～2 500 nm。选取位于玉米中部的叶片在室内密闭的环境下进行光谱采集。以50 W 卤素灯为光源，使用4°视场角的探头垂直于玉米叶片表面，探头距离叶片表面40 cm。为了避免因光源强度分布不均匀导致暗电流噪声影响光谱数据质量，每次测量玉米叶片光谱前，先用白板进行标准化处理。光谱仪连续测量玉米叶片反射光谱3 次，随后平均3 次测量值作为该胁迫浓度下玉米叶片光谱数据。

将采集完光谱的玉米叶片洗净后放到烘箱中烘烤至质量恒定，用剪刀剪成小块装入样品袋中并贴上标签。经过预处理后，在化学分析实验室中采用WFX-120 原子吸收分光光度计进行玉米叶片样品Cu2+含量测定，每份叶片样品平均分成3份，将实验测出的Cu2+含量值取平均得到该胁迫浓度下玉米叶片样品中Cu2+的含量。实验中测得0、100、300、500µg·g-1胁迫浓度下的玉米叶片Cu2+含量分别为：9.77、29.62、75.78、54.51µg·g-1。

1.2 连续小波变换

连续小波变换通过小波基函数将光谱信号分解为低频和高频两个部分，在低频部分具有较低的时间分辨率和较高的频率分辨率，而在高频部分具有较高的时间分辨率和较低的频率分辨率，适合于各类信号噪声的处理。连续小波变换公式为[17]：

式中：Wf（a，b）为小波变换结果，f（λ）为玉米叶片高光谱反射率，ψa，b（λ）为小波基函数，a为尺度因子，b为平移因子。

1.3 希尔伯特-黄变换

希尔伯特-黄变换（HHT）主要用于非线性非平稳信号的处理，包含经验模态分解（Empirical mode decomposition，EMD）和Hilbert变换两部分[18]。

1.3.1 经验模态分解

EMD 可以把任意一个复杂的非平稳信号自适应地分解为有限个窄带信号，称为本征模态函数（Intrinsic mode function，IMF），每个IMF 需要满足以下两个条件：（1）零点数目与极值点数目相同或至多相差1 个；（2）函数关于局部平均对称。最终将光谱信号x（t）分解为若干个IMF 分量ci（t）与单调残差rn（t）之和：

其具体算法如下：

（1）求取信号x（t）的极值点；（2）采用三次样条函数拟合信号的极小值点与极大值点得到上下包络线；（3）计算上下包络的平均值；（4）将原始信号减去上下包络的平均值，提取IMF分量；（5）重复以上步骤至均值趋近于零；（6）计算残差；（7）重复运算直到残差中不含IMF函数。

1.3.2 Hilbert变换

信号经过经验模态分解后，就可以对每一个IMF作Hilbert变换，公式为：

式中：H[ci（t）]为第i阶IMF 分量ci（t）的Hilbert 变换函数。

根据欧拉公式可得其解析信号为：

式中：ai（t）为解析信号的幅度或能量；θi（t）为解析信号的相位。

由此可以求出IMF分量的瞬时频率wci（t）：

则原始信号x（t）的希尔伯特谱表示为：

式中：n为IMF个数；Re为残差项。

则Hilbert瞬时能量谱E（t）为：

1.3.3 CWT-HHT铜污染信息探测流程

基于CWT 与HHT 所构建的CWT-HHT 探测方法处理具体过程如下：（1）对原始玉米叶片光谱数据进行CWT 分解，有利于光谱信号降噪；（2）提取连续小波分解后的第5 层小波系数进行EMD 分解，得到有限个单一频率的IMF 分量；（3）用Hilbert 变换对IMF分量进行处理，获取解析信号，以提取信号的瞬时频率信息；（4）将解析信号表示为希尔伯特谱形式；（5）对各IMF 分量的希尔伯特谱的平方在频率域进行积分，得到光谱信号的瞬时能量；（6）根据波长、频率和能量构建三维信息图，选取瞬时能量峰值E作为特征参数，以探测不同浓度铜污染下玉米叶片光谱重金属污染信息。

1.4 数据处理

本研究采用Matlab 2018b 编程对玉米叶片光谱数据进行连续小波变换与希尔伯特-黄变换处理分析，利用SPSS 进行相关性分析，在Excel 2013 中进行显著性检验，采用Origin 2018进行制图。

2 结果与讨论

2.1 光谱波段选择

实验采集的一组玉米叶片光谱如图1 所示。从光谱曲线可以看出，不同浓度Cu2+胁迫下叶片光谱均表现出高度的相似性，在可见光范围，550 nm 附近因光合色素对绿光的强烈反射形成绿峰，650 nm 附近因叶绿素对红光的吸收形成红谷，750 nm 附近因反射率急剧上升形成叶片光谱最明显的红边特征；在近红外波段范围，760～1 250 nm 附近受叶片细胞内部结构影响形成一个高反射平台，在1 450 nm 与1 930 nm呈现出两个明显的水吸收带，难以直接进行重金属污染探测[19]。因此，需要选取合适的光谱波段以提取玉米叶片重金属污染光谱特征。为了提高探测精度，去除噪声较大的350～400 nm 波段以及受水分强吸收影响的1 300～2 500 nm 波段光谱，选取400～1 300 nm 波段进行分析。

图1 不同Cu2+胁迫浓度的玉米叶片光谱曲线Figure 1 Spectral curves of corn leaves under different Cu2+stress gradients

2.2 基于CWT-HHT污染监测

2.2.1 小波系数提取

选取不同的小波基取得的结果不尽相同。经过多次对比分析，采用Daubechies 小波系中的“Db5”小波函数对重金属Cu2+胁迫下的玉米叶片光谱进行连续小波变换，提取分解至第5层的小波系数，如图2所示，可以看出，在400～1 300 nm 波段范围内小波信号较原始光谱曲线更加平滑，并且在指示植物重金属污染的红边（680～750 nm）这一敏感位置小波系数变化剧烈，可见连续小波变换既有效降低了光谱噪声又增强了玉米受重金属Cu2+污染信息。

图2 不同Cu2+胁迫浓度的玉米叶片光谱小波系数Figure 2 Spectral wavelet coefficients of corn leaves under different Cu2+stress gradients

2.2.2 CWT-HHT铜污染信息探测

将连续小波变换后的光谱数据进一步作HHT 处理得到希尔伯特谱，并将希尔伯特谱的平方在频率域上进行积分，构建叶片光谱信号的波长-频率-能量三维信息图。如图3 所示，横坐标表示波长，纵坐标表示对应波长信号的频率，能量以颜色深浅显示，颜色越深对应的波段能量就越高。

图3 不同Cu2+胁迫浓度下的希尔伯特瞬时能量Figure 3 Hilbert instantaneous energy under different Cu2+stress gradients

由图3 分析可知，随着Cu2+胁迫浓度的不断增加，玉米叶片光谱的Hilbert 瞬时能量峰值呈现先增加后降低的趋势。不同浓度Cu2+污染下，瞬时能量峰值均在700 nm 波段附近出现，并且在1 100 nm 波段附近出现第2 个能量峰，表明叶片光谱在这两个波段存在能量波动，并且在700 nm 波段附近能量波动异常明显，这与提取的小波系数剧烈变化位置一致。而信号频率集中分布在0～0.5 Hz较窄的部分，这是因为原始光谱经连续小波变换后提取的第5 层小波系数的频率单一，去除了不同频率信号的干扰。

为了定量描述玉米受重金属Cu2+污染程度，提取瞬时能量峰值E作为探测指标，图4 为能量峰值E与玉米叶片Cu2+含量变化趋势。可以看出，不同浓度Cu2+胁迫下，能量峰值E与叶片Cu2+含量变化趋势一致。因此，可以认为瞬时能量峰值E作为玉米重金属Cu2+污染监测的指标具有可行性。

图4 能量峰E与叶片Cu2+含量变化趋势Figure 4 Energy peak E and transformation trend of Cu2+content in leaves

2.3 CWT-HHT探测方法应用与验证

由于玉米受到不同浓度重金属Cu2+污染后其叶片光谱信号差异微弱，难以直接进行污染程度判别。实验首先对采集的玉米叶片光谱进行CWT 处理，获得单一频率的小波系数。叶片光谱信号经CWT 处理后，能够降低光谱信号噪声并增强玉米重金属污染信息，可以提高HHT算法探测玉米受Cu2+污染程度信息的精度。随后利用HHT 算法获取光谱信号的三维信息图并提取信号的特征值，以描述在不同浓度Cu2+污染下玉米污染程度信息。因HHT 算法能够自适应分解光谱信号，在消除人为因素干扰的同时可以进一步突出信号本身的局部特征，适合玉米在重金属Cu2+污染下光谱畸变特性的探测研究。

为了验证CWT-HHT 方法在玉米Cu2+污染信息探测方面的有效性与优越性，将其与红边位置（REP）、红边归一化指数（NDVI705）和红边植被胁迫指数（RVSI）3 个常规的植被指数监测重金属污染方法进行对比分析。其中，REP 为反射光谱一阶微分在670～780 nm 最大值对应波长，NDVI705计算公式为（R750-R705）/（R750+R705），RVSI 计算公式为[（R712+R752）/2]-R732。

不同浓度Cu2+胁迫下叶片光谱的CWT-HHT 探测方法结果E值和REP、NDVI705、RVSI 与叶片中Cu2+含量的相关系数以及显著性检验结果如表1 所示。分析发现，E值与玉米叶片中Cu2+含量的相关系数最高达到0.980 3，表明玉米受重金属Cu2+污染程度越大，叶片光谱的瞬时能量特征越显著。同时，将实验获取的玉米叶片中Cu2+含量与各光谱参数进行线性拟合分析，拟合结果如图5 所示。结果表明CWTHHT 探测方法提取的E值与叶片中Cu2+含量拟合程度最优，其拟合决定系数R2最高，达到了0.949 8，远高于常规植被指数监测方法的拟合系数。综上所述，CWT-HHT方法可以有效探测出Cu2+污染下玉米叶片光谱微弱的差异信息，也验证了其在农作物重金属污染监测方面的可行性。