利用作物生长模型和时序信号甄别水稻镉胁迫

孔 丽，刘美玲，刘湘南，邹信裕

（中国地质大学（北京）信息工程学院，北京 100083）

0 引 言

随着经济的不断发展，工业化不断推进，土壤重金属污染成为当今世界面临的重大生态环境问题之一[1]。水稻是中国的主要农作物，产量位居世界第二位。水稻重金属污染不仅会影响水稻的粮食产量，同时也会威胁人类的身体健康。传统的实地调查虽然能够准确检测水稻中的重金属浓度[2]，但费时费力，且难以满足大区域尺度监测的要求。由于遥感技术能够实时、大范围、连续动态获取目标地物对环境胁迫的光谱响应信息，因而利用遥感技术对水稻重金属污染现状进行监测和评估逐步受到国内外学者的广泛关注。

目前，许多学者对水稻重金属胁迫遥感监测开展了一系列的研究。赵利婷等[3]利用遥感同化WOFOST（World Food Studies）模型模拟水稻根重进行水稻重金属污染胁迫分析。Liu等[4]利用多时相的Sentinel-2数据融合光谱信息和时空信息监测水稻的重金属胁迫。近年来，为提高农作物重金属遥感识别精度，一些学者利用了多尺度小波变换的光谱分析方法对农作物污染胁迫状况进行识别和评价[5-6]；另一些学者将农作物生理功能特征变化（如蒸腾作用、光合作用、呼吸作用等）在光谱上的响应作为农作物重金属胁迫评估指标来揭示农作物重金属污染状况，如Jin等[7]构建了基于农作物形态和生理功能变化参数的重金属胁迫理论模型。Tian等[8]将胁迫效应的谱特征与时间特征结合去识别重金属胁迫。Tang等[9]提出了基于时空特征指标的水稻重金属胁迫识别模型，构建了融合水稻胁迫的年内特征、年际特征以及空间特征的时空指标，实现了水稻重金属胁迫的监测。上述农作物重金属污染监测研究主要是在较小的空间尺度且已知重金属污染条件下，假设光谱变异仅由重金属胁迫引起（忽略其他胁迫影响），如何剥离其他胁迫对重金属胁迫的影响是开展重金属胁迫遥感准确识别亟待解决的问题。

在复杂的农田生态系统中，农作物会受到复杂胁迫的影响，包括短期胁迫（如病虫害、水分胁迫等）和长期胁迫（如重金属胁迫等）[10]。据研究表明，重金属污染具有隐蔽性、持久性和不可逆性等特点[11-12]，因此，重金属胁迫存在于农作物整个生长周期，对农作物造成的影响是持续稳定的，且多年都呈现一个稳定相似的状态。而短期胁迫的特点是年际间差异较为明显，往往存在于农作物的某一个或某几个生长周期中。由此，农作物不同胁迫效应呈现的时间特征差异为区分农作物重金属胁迫和其它胁迫提供了可行性，探索一种适合不同胁迫特征分量的分离与提取的时序信号分解算法，是水稻重金属胁迫准确识别的关键。水稻胁迫特征参数通常包含水稻固有生长趋势，年内长期胁迫和随机噪声或年内短期胁迫。时序信号分解可以有效地获取感兴趣的信号，集合经验模态分解算法（Ensemble Empirical Mode Decomposition，EEMD）是一种自适应时序信号分解方法，可以很好地揭示信号的非平稳性和非线性信息。该方法在信号分析中得到了广泛地应用[13-15]。

为揭示水稻不同胁迫呈现的时间效应差异，有必要开展水稻胁迫特征参数的时间连续模拟与计算。WOFOST模型根据气象和土壤条件以一天时间为间隔[16]模拟作物的叶面积指数（Leaf Area Index，LAI）生成密集的时间序列，准确实现环境胁迫下作物生长过程的动态模拟，从而获取时间连续的作物胁迫特征参数，为不同胁迫特征分量的分离与提取提供基础数据。陈艳玲等[17]采集冬小麦多个关键生育时期的生理生化、农田环境、气象等数据，构建基于WOFOST模型和遥感LAI数据同化的区域尺度冬小麦单产预测模型。黄健熙等[18]重点优化WOFOST模型中与品种相关的积温参数，该模型在全国尺度取得了较高的模拟精度。以上研究表明了遥感同化WOFOST模型在监测作物的长势等方面具有一定的优势。针对水稻重金属胁迫具有“年际稳定性”特征，因而水稻胁迫特征参数在不同年份存在相似性。动态时间规整（Dynamic Time Warping，DTW）方法衡量时间序列曲线的相似度，该算法已经相当成熟，应用广泛，例如语音识别、信息检索、土地覆盖分类等等[19-21]，这些研究成果为重金属胁迫研究提供了重要的参考。相对于传统的距离度量方法，DTW更适合度量不同长度以及不同节奏的时间序列的相似度，这一优势正好可以有效避免不同年份水稻物候期（如分蘖期）的差异性所引起的相似性计算误差，由此，DTW算法适合度量水稻胁迫特征参数在不同年份的相似度。基于前面的分析，本研究以湖南省株洲为试验区域，利用多年Sentinel-2数据同化作物生长模型得到不同年份的LAI时间序列，通过EEMD方法对其进行分解，剥离年内短期胁迫或噪声，提取可能包含重金属胁迫在内的年内长期胁迫；在此基础上，利用DTW计算其年内长期胁迫信号与健康水稻胁迫信号之间的距离，通过相似性度量准确识别水稻重金属胁迫。

1 材料与方法

1.1 试验区概况

试验区位于湖南省株洲市（113°05′E～113°15′ E，27°37′N～27°45′N）（图1）。该地区属于亚热带季风性湿润气候，四季分明，雨量充沛，光照充足，年平均气温为16～18 ℃，年降水量大约为1 445 mm，土壤为富含多种营养物质的红壤，有机质含量（质量分数）2%～3%，适宜水稻生长，是国家重要的商品粮食生产基地。本试验区覆盖一幅Sentinel-2影像的大小为1 000×1 000像元，其范围约为100 km2。据调查，湘江流域灌溉的土壤中重金属镉的含量较高，是主要污染物[22-23]。

1.2 数据获取与处理

本研究所采用的数据包括遥感数据、作物生长模型所需的输入数据以及实测的土壤重金属含量数据。其中作物生长模型所需的数据为气象数据和作物生长参数，气象数据来自中国气象站（http://www.cma.gov.cn），包括逐日最低温、逐日最高温以及逐日照时数。而作物参数主要有物候因子、作物初始参数、绿叶面积、CO2同化和干物质分配。土壤重金属数据针对25个采样点均选取距土壤表层10 cm深处的土壤，晒干后使用HNO3-HF-HCLO4消毒法在微波消解仪（MARS X，CEM，USA）进行消毒，并采用电感耦合等离子体质谱ICP-MS（7500a，Agilent Technologies，USA）测定Cd、Pb等元素含量。每个土壤样品平均测定3次，结果以均值表示。遥感数据是2017—2019年的Sentinel-2卫星数据（分别为2017年7月17日和24日、2017年8月6日、2017年9月15日；2018年6月17日、2018年7月17日、2018年8月8日、2018年9月10日；2019年6月14日、2019年8月16日、2019年9月15日、2019年9月22日）（https://scihub.copernicus.eu），首先对遥感影像数据进行预处理，包括辐射校正和大气校正等，然后在遥感影像中选取训练样本，利用随机森林算法提取研究区域的水稻分布范围。

1.3 研究方法

本文采用作物生长模型、EEMD与DTW方法有效地甄别并定量分析水稻重金属胁迫，具体的流程如图2所示。主要分为4个步骤：1）对数据进行预处理，基于遥感影像数据获取研究区域水稻的LAI；2）运用作物生长模型模拟水稻的LAI时间序列；3）通过EEMD方法将原始LAI时间序列分解成不同时间尺度的组分，分析合成含有重金属胁迫的信号；4）计算分解后含有重金属胁迫的信号与健康水稻之间的DTW距离并进行归一化处理，监测重金属胁迫的程度。

1.3.1 基于遥感影像反演水稻LAI

LAI是植被主要的理化参数之一，本文利用LAI作为遥感与作物生长模型之间的同化量。目前，获取LAI的方法主要有经验模型、物理模型与混合模型，经验模型主要是通过建立实地观测数据与遥感影像数据之间的数量关系来获取LAI观测值。统计模型操作简单，易于实现，精度较高，因此，本文利用实测LAI和Sentinel-2数据获取的归一化植被指数（NDVI）值建立适合本研究区的LAI经验反演模型。对比不同的拟合方程选取拟合精度最佳的拟合模型[24]。具体的LAI计算式如下：

式中NIR和R分别代表近红外波段和红波波段反射率。

1.3.2 基于WOFOST模型获取长时序LAI

为了将不连续的遥感信息转化为时间连续的作物信息，本文利用遥感与WOFOST作物生长模型同化方法模拟作物的生长过程。该模型基于作物生长发育的基本过程，解释了作物的生长，如光合作用以及呼吸作用。它能够模拟温度、水分对作物生长发育形成的胁迫时期和程度[25]。WOFOST模型主要有3种水平的作物生长：潜在生产水平、水分限制生产水平和营养限制生产水平。它的输入参数有作物、土壤和气候参数，本文研究综合影响水稻长势的各方面因素，包括稳定胁迫（如重金属胁迫等）、短时胁迫（病虫害等）以及特定时间发生的灾害（干旱等），依据胁迫对作物生长的影响机理，选择合适的胁迫因子嵌入到日总CO2同化的生产过程。

式中CVFf为胁迫水平下的日总CO2同化物产量，kg/hm2；CVF代表潜在生产水平下的日总CO2同化物产量，kg/hm2；f为胁迫因子，取值范围为[0,1]，f值越高，说明作物生长过程受到的胁迫越小，反之，胁迫越严重。

同化算法是模型同化的关键，同化算法的选择直接影响模型模拟的效率和精度。本文采用原理简单、易于实现的粒子群优化（Particle Swarm Optimization，PSO）算法[26-27]，其主要的思想是通过不断调整模型参数CVR（Efficiency of Conversion into Roots，干物质转化为根质量的效率），使得WOFOST模型模拟值与遥感反演实测值之间的差异逐渐减小，收敛条件是全局最优解不再变化或达到最大迭代次数，最终输出CVR的最优解即同化结果。本文使用的代价函数为

式中LAIm为实测值序列；LAIs为模拟值序列；N为实测时相数。

通过比较模型模拟的LAI与Sentinel-2数据反演的LAI，得出优化后的模拟结果更接近LAI真实值，通过WOFOST模型模拟出连续的水稻胁迫参数LAI。

1.3.3 基于EEMD分解LAI时间序列

EEMD原理较为简单，该方法是针对经验模态分解（Empirical Mode Decomposition，EMD）方法的不足向原始信号中添加白噪声，使其信号在不同尺度上具有连续性，从而避免出现模态混叠的现象，达到一个更好的分解结果。EEMD分解方法两个非常重要的参数是白噪声的幅值系数和总体平均次数，白噪声的幅值系数需要人为经验进行确定，添加的白噪声幅值系数过大，会降低分解的精度，过小，则无法解决EMD中存在的模态混叠现象。Huang等[28]提出分解以高频组分为主的原始信号时，白噪声幅值系数应相对减小，相反，白噪声幅值系数应相对增大。通常白噪声幅值系数的取值范围为0.01～0.40。总体平均次数会影响分解效果，为了同时兼顾效率与精度，总体平均次数一般选择100～300。参考前人的研究结果，结合本论文试验的需求，将总体平均次数设置为100次，白噪声幅值系数设置为0.3[10]。

时间序列经过EEMD分解后可以得到不同分量[29]。这些分量由3个波动组分和1个残差组成，其中3个组分是年内波动组分、年间波动组分和年际波动组分[30-31]。年内波动组分为时间尺度通常小于一年的高频信号。高频信号可以检测信号中存在的异常现象，其他环境胁迫和噪声信息存在于高频信号中。作物固有的生长趋势特征与季节项相似，时间尺度通常为一年，有一个固定的周期[32]。本文将波动周期为一年的年间波动组分表示为作物的固有生长趋势。去除这些因素的影响，剩余的年际波动组分和残差可能含有重金属胁迫信息。本文选取皮尔逊相关系数衡量年际信号之间的相似度。

1.3.4 动态时间规整（DTW）

DTW支持时间轴的伸缩和弯曲，能够对不同长度的时间序列进行相似性度量[33]。它是用满足一定条件的时间规整函数描述两个时间序列之间的对应关系，求出两个时间序列累加距离最小时对应的规整函数。其路径的选择满足3个条件：边界条件、连续性和单调性。DTW通过寻找到最小的扭曲路径代价来对齐不同的时间序列。

本文引入DTW方法计算时间序列曲线的距离大小，因为水稻在不同年份、不同地域条件下，其物候期也存在一定的差异，而DTW能够有效地解决这一类问题。通过计算分解后的时间序列与健康水稻时间序列之间的DTW距离，可以监测水稻重金属胁迫水平。

1.3.5 重金属胁迫程度的提取

本研究利用WOFOST模型模拟出健康水稻的LAI时间序列，经过EEMD方法分解提取的无胁迫信号作为参考序列，原始LAI时间序列提取的含有重金属胁迫信号作为查询序列，计算两者之间的DTW距离，即得到胁迫指数（Stress Index，SI）。

式中SIi表示水稻像元i的重金属胁迫指数，Si为水稻像元i当年的原始LAI时间序列分解提取的可能含重金属胁迫序列，Sh代表当年健康水稻LAI分解提取的无胁迫信号，Ddtw表示两者之间的DTW距离。

为了更加直观地了解水稻重金属胁迫的程度，本文对胁迫指数SI进行归一化处理，归一化的公式如下：

式中SIn表示归一化胁迫指数，SIi0.05和SIi0.95分别表示胁迫指数5%和95%分位数，SIn的值范围为[0,1]，越接近1代表该水稻像元i受到的重金属胁迫程度越高。

在本文中，采用标准方差分析[10]将水稻重金属胁迫程度划分为轻度、中度和重度3个等级，其SIn值分别介于0～0.40、＞0.40～0.75、＞0.75～1之间。

2 结果与分析

2.1 水稻LAI的分布

通过对Sentinel-2数据的LAI进行反演，得出试验区内所有像元的LAI值，根据随机森林方法提取水稻区域进而叠加提取水稻像元LAI值（图3）。图3中白色表示非水稻的区域，红色代表更低的LAI值，绿色代表更高的LAI值。 从图中可以看出试验区西北和东南区域的LAI值相对较低，偏东北区域的LAI值相对较大，且分布较为密集。

2.2 EEMD信号分解

基于EEMD将LAI时间序列分解为7个不同波动频率的IMF分量和1个残差项，通过不同分量的合成可以得出年内波动组分、年间波动组分以及年际波动组分（图 4）。图4中的IMF1、IMF2、IMF3和IMF4的波动周期小于一年，该4个分量之和表示年内波动组分；而IMF5分量的波动周期与原始LAI时间序列的波动周期是具有一致性的，因此将IMF5分量对应于年间波动组分；剩余的IMF6、IMF7对应于年际波动组分，波动时间周期大于一年。

年内波动组分包含的信息具有较高的时间分辨率，时间跨越尺度小于一年，信号波动比较杂乱，其中不仅包含其他短期胁迫还包含自然环境中的噪声等信息。可见IMF1～4分量的波动周期均小于一年，环境中的噪声信息虽然长期存在并对信号产生一定的波动影响，但是噪声信号没有一定的规律。

通过计算相邻年份年内波动组分的皮尔逊相关系数可知，2017年和2018年的相关系数为0.13，2018年和2019年的相关系数为0.33，表明相邻年份年内波动组分之间的相关性很弱，因此，能更好地说明年内波动组分是表示短期胁迫和噪声等不稳定的信号。

IMF5分量的波动周期为一年，与原始水稻LAI具有相同的波动规律，因此这个分量能够表示水稻本身固有的生长趋势。计算年间波动组分与健康水稻LAI之间的相关系数，2017年、2018年和2019年与健康水稻之间的相关系数分别为0.96、0.85、0.90，均表现极强的相关性。说明年间波动组分揭示了水稻固有的生长趋势。

原始LAI时间序列由持续的环境胁迫、作物的固有生长趋势和其他胁迫信号共同组成，通过对上述两种组分的剔除，重金属胁迫信息存在于剩余的分量中，即年际波动组分与残差之和（IMF6、IMF7和残差）。

2.3 水稻重金属胁迫程度

2017—2019年水稻重金属胁迫的空间分布特征如图 5所示，研究结果表明，不同年份的水稻重金属胁迫的空间分布具有很大的相似性，在空间上的分布呈现一种稳定的特征。受重金属胁迫程度较低的区域大部分位于试验区域的中部地区。而重金属胁迫程度 较高的地区主要是集中在试验区域西部、东北部以及偏东南地区。湘江东侧和西侧的工厂较为密集，重金属胁迫的空间分布与工厂空间分布具有相似的趋势，可能是由于早期的污水灌溉导致土壤重金属含量超标。

进一步对试验区域重金属胁迫程度进行统计分析（表1），由表1可知，2017—2019年轻度重金属胁迫分布面积比例较大，不同年份间其所占比例相差较小，分别为38%、44%和45%。重度重金属胁迫分布面积比例约占29%，且不同年份的分布面积比例几乎一致，2017年与2018年仅相差0.17个百分点，2018年与2019年相差不到1个百分点。从3 a的重金属胁迫面积统计数据可知土壤中重金属含量存在时空稳定性。

表1 研究区不同程度重金属胁迫水稻面积比例Table 1 Percentage of rice area under different heavy metal stress in the study area %

综上所述，通过计算多年的归一化胁迫指数SIn可以准确确定不同水稻重金属胁迫水平的分布情况。

2.4 精度评价

研究利用野外实测的土壤重金属Cd的含量数据对水稻重金属胁迫监测的结果进行精度评价，通过计算实测的土壤Cd含量与归一化胁迫指数之间的皮尔逊相关系数进行验证。两者之间的相关系数高达0.851（图6），由此说明归一化胁迫指数SIn对水稻重金属胁迫反应比较敏感。

3 讨 论

本文借助EEMD方法进行胁迫信号分解，剥离其他胁迫以及噪声等因素对重金属胁迫的影响，可以更加准确地研究试验区域的水稻重金属胁迫。由于重金属胁迫在土壤中可移性小，难以降解，因而土壤重金属污染具有长期性[34-35]。有研究基于长时序的LAI结合EEMD方法提取水稻重金属胁迫特征，并进行胁迫特征的稳定性分析[8]。目前，越来越多的研究将DTW方法应用于遥感研究领域中，利用动态时间规整算法计算得到的距离作为相似性度量的指标来进行物种的识别分类[35]，也有研究通过DTW方法构建指标分析重金属胁迫[9]，但是该方法没有考虑水稻同时受到多种胁迫共同影响的情况，因此建立的模型较为简单，在实际的应用中会存在局限性。因此如何在复杂的生态环境中研究重金属胁迫是一个重要的课题。本研究利用EEMD与DTW相结合的方法先将其他胁迫与噪声等因素排除，计算与健康水稻之间的距离判别水稻受重金属污染的程度，在时间和空间上表明重金属胁迫具有稳定的特性，相邻年份之间的DTW距离体现年际变化特征，可以作为衡量重金属胁迫稳定性的一个重要指标。

本文选择利用DTW方法研究重金属胁迫的一个重要的原因是DTW作为一种相似性度量的方法，它可以对时间序列进行拉伸和缩放，将两个时间序列进行匹配比较，从而使识别效果更佳。湖南由于种植一季稻，插秧时间等的不同会导致水稻不同年份之间物候期有所差异，因此引入DTW方法可以很好地消除物候期不同导致的距离增加，从而提高水稻重金属胁迫研究的精确度。

SIn是在无其他胁迫和环境影响因素的情况下计算得出的归一化胁迫指数，根据其值的大小可以判断出水稻重金属污染的程度。试验结果表明重金属分布的范围在3 a间变化较小，在空间上表现出一定的稳定特征。因此，引入时序稳定检测方法，对水稻重金属胁迫进行遥感监测具有很高的可信度。

本研究仅采用3 a的数据对水稻重金属胁迫进行研究，今后可以通过更长时序遥感数据对其进行分析；另外，本研究利用土壤重金属Cd的含量来验证水稻重金属胁迫水平，在今后的研究中将进一步通过测定水稻中的重金属含量验证研究结果。

4 结 论

本文利用Sentinel-2数据同化作物生长模型获取了叶面积指数（Leaf Area Index，LAI）时序数据，利用集合经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）方法分解LAI时间序列，提取出含有重金属胁迫的信号，利用动态时间规整（Dynamic Time Warping，DTW）方法计算其胁迫信号与健康水稻LAI之间的距离大小，监测水稻重金属胁迫水平。结果表明：

1）原始LAI时间序列由持续的环境胁迫、作物的固有生长趋势和其他胁迫信号共同组成，经EEMD分解为7个不同的分量和1个残差项，通过分析得出重金属胁迫信号存在于IMF6、IMF7和残差的组合分量中；

2）归一化胁迫指数是水稻重金属胁迫敏感的参数，与土壤重金属含量的相关系数为0.851，水稻受到的胁迫程度越高，归一化胁迫指数值越大，反之越低；

3）在试验区中，水稻重度重金属胁迫的分布面积比例相对较低，且主要集中在西部、东北部以及偏东南地区，水稻重金属胁迫具有一定的时空稳定性。