基于EWT-ARIMA 组合模型的银杏液流预测与因子关联分析*

2023-11-10王子祥李颜娥武斌徐达宇吴斌

电子技术应用 2023年10期

王子祥，李颜娥，武斌，徐达宇，吴斌

（1.浙江农林大学数学与计算机学院，浙江杭州 311300；2.浙江省林业智能监测与信息技术实验室，浙江杭州 311300；3.林业感知技术与智能装备国家林业局重点实验室，浙江杭州 311300）

0 引言

进入二十一世纪，人类正面临气候恶化这一全球性问题。碳排放造成的全球变暖会对人类和生态系统造成影响严重且不可逆转的危害[1]，利用植物生态系这一巨大天然碳库去实现固碳目标被认为是经济有效的方法[2]。同时植物也是城市公园内固碳效益的主体，在常见的城市公园树木中，银杏的固碳能力明显高于其他植物，属于第一梯队[3]。除此之外，银杏还是一种的珍贵的中药材和观赏植物，具有广泛的生态、经济和社会价值[4-5]。为了生态环境的可持续发展，寻找一种准确评估树木蒸腾耗水的稳健方法十分有必要[6]，而植物蒸腾耗水量的99%以上来自茎干液流，精确测定茎干液流量能够反映林木的蒸腾耗水量[7]。液流法也成为了近年来茎分和林分尺度蒸腾耗水研究的热点方法[8]。

基于树木液流的蒸腾评估，需要克服液流这类自然序列普遍存在的非线性与高随机性。因此本文提出一种基于经验小波变换(Empirical Wavelet Transform,EWT)和ARIMA(Auto Regressive Integrated Moving Average)组合的银杏液流预测模型。银杏液流序列先经过经验小波分解,得到平稳的多分辨分析分量(Multi Resolution Analysis,MRA)，针对各分量信号分别构建ARIMA 模型进行预测并重构集成为目标预测信号，并与单一ARIMA 模型进行对比分析。此外，有许多研究证明树木液流与环境因子有着密切的关系。但以往的研究多采用线性分析工具，在理论上有局限性，因此本文借助无需模型假设的转移熵（Transfer Entropy,TE）对环境因子和银杏液流时滞内的因果关系进行了探究。

1 相关工作

发现关联与变量选择在科学实践中至关重要，但目前针对液流与环境因子关系的模型多采用单变量分析、多元线性回归[9-10]，特定的分布假设使得这类方法在理论上存在不足。随着技术的发展，许多智能算法被广泛应用到具有复杂模式的序列中，获得了不错的效果，EWT 技术[11-13]就是其中之一。常见做法是将复杂序列视作信号，经EWT 分解为若干稳定分量，再依据各自的数据特征进行重构，从而取得一个优良的预测效果。

2 基于EWT-ARIMA 组合模型的银杏液流预测

2.1 经验小波变换EWT

2.1.1 从傅里叶变换到经验小波变换

傅里叶变换是一个非常强大而灵活的方法，其原理说明由正弦基函数组成的合成信号可以分解为有限个正弦波或余弦波的叠加，且这些分量信号的频率是原始信号频率f0的整数倍。

式中,f0为基频，A0/2 代表直流系数，An代表幅度，ϕn被称为相位。根据振幅和相位等信息可以利用反变换从而恢复信号的最初波形。

对一个正弦信号和一个非线性信号做傅里叶变换，如图1 所示。

图1 正弦信号与非线性信号傅里叶图

根据图1 可以显示出非线性信号比正弦信号具备更宽的波谷和更尖锐的波峰。

根据图2 伟尔奇周期图，正弦信号和非线性信号都具备5 Hz 上的高功率，符合二者的定义。但是在10 Hz处，非线性信号还呈现出另外一个峰值，这是因为傅里叶变换只能基于固定的线性基函数，导致其无法表达信号复杂的震荡模式。

图2 正弦信号与非线性信号伟尔奇周期图

经验小波变换（EWT）是Gilles 提出的一种分解技术[14]，它结合了小波的形态优势和经验模态分解（EMD）的自适应性。可以有效减少模态混叠和端点效应，解决原始序列模型因非线性和高随机而导致的预测精度低的问题。

2.1.2 经验小波变换原理

经验小波变换原理大致可以总结如下：首先对序列进行傅里叶变换，得到傅里叶频谱并将频率规范至[0,π]，找到并以降序排列频谱中的N个局部极大值所对应的频率ωn（n=1，2，3，…，N-1）。将两个连续局部极大值的中点设置为支撑边界，然后执行频带分割：

其中，ωi和ωi+1为两个连续局部极大值的频率，支撑边界集合Ω={Ωi}i=1,2,…,N-1。

频带划分过程如图3 所示。其中τn代表过渡空间段一半的长度。

图3 EWT 频带分割示意图

然后基于支撑边界所形成过渡空间设计自适应小波滤波器组。最后将经滤波后的模态转换成尺度函数与小波函数，以获得信号分量。

2.2 ARIMA 预测模型

ARIMA 模型(Autoregressive Integrated Moving Average model)是一种用于时间序列预测的统计方法，它可以捕获时间序列中的不同特征，如趋势、季节性、周期性等，曾极大地促进了趋势预测研究[15-17]。ARIMA 模型由三部分组成：自回归(AR)、差分(I)和移动平均(MA)。ARIMA 模型有3 个参数(p,d,q)，分别代表自回归项数、差分次数和移动平均项数。

ARIMA 预测原理：ARIMA 模型预测先通过ADF 方法来检验时间序列是否具备平稳特征，若不平稳则进行d阶差分直至差分后转化为平稳的时间序列；然后根据自相关系数图(ACF)和偏自相关系数图(PACF)或者基于AIC 与BIC 准则选取模型参数，最后确定模型参数并进行检验和预测。

3 实验结果与分析

3.1 实验数据集

本文采用了一份来校园内银杏树的液流数据，数据由单针热脉冲树干液流监测仪分别于2021 年8 月29 日～9 月4 日期间，在一棵银杏树上进行数据的采集而成。试验地点位于浙江农林大学东湖校区(30.24°～30.26°N，119.72°～119.73°W)，在监测树干液流的同时，同步监测大气温度(Air Temperature,TA)、水汽压亏缺（Vapor Pressure Deficit,VPD）、大气湿度（Air Relative Humidity,RH）、土壤湿度（Soil Moisture,SM）等气象因子。在树木茎干的茎向直线上将热源探针插入，然后分别在其上、下部位约6 mm 和8 mm 处插入两探针,开始进行数据监测，如图4 所示。

图4 树干液流监测仪以及探针距离示意图

8 月29 日～9 月4 日监测期间，属于夏季的晴天气温较高。绘制银杏液流变化趋势如图5 所示。

图5 银杏液流序列变化趋势图

其中,图5（a)代表银杏液流完整的变化趋势，图5(b)是选取了图5(a)其中8 月30 日数据放大后的效果。通过观察图5 银杏液流曲线可知，银杏液流呈现周期性变化，日出后开始快速上升，直到中午前后达到峰值。

接下来对完整的银杏液流序列进行经验小波变换，经过EWT 分解，得到两组多分辨分析分量信号（mra1 和res）以及对应的频谱，发现经分解后mra1 分量信号大致以零轴为上下对称，序列较平稳，对应的频谱也有一致的变化趋势。同时，res 分量信号被分解得比较彻底，分解后的曲线十分平滑，为后续的建模预测建立了基础。分解效果图6 所示。

图6 银杏液流序列EWT 分解与mra 分量频谱对照图

3.2 ARIMA 模型预测

3.2.1 模型平稳性检验

ARIMA 模型要求时序数据是稳定的,接下来对银杏液流和经分解得到的两组分辨率分析分量信号进行ADF(Augmented Dickey-Fuller test)检验，ADF 通过判断序列是否存在单位根来判断序列的平稳性。ADF 单位根检验结果如表1 所示。

表1 原序列与mra 分量ADF 单位根检验

3 组序列P 值分别为0.0245、0、0.003，均小于显著性水平且T 检验值均小于3 个临界值，表明3 组序列所检验的参数落在拒绝域，拒绝原假设，认为不存在单位根，故均为平稳序列，可以通过ARIMA 建模。

3.2.2 模型参数定阶与检验

银杏液流预测首先要确定参数p与q，通过银杏液流的自相关系数图（ACF）和偏自相关系数图（PACF）所呈现的拖尾或截尾特征来进行定阶。

根据图7，可判断ACF 与PACF 属于拖尾，p与q可定阶在3 左右，但这种定阶方式存在一定主观性。除此之外还可以使用赤池信息准则(AIC)与贝叶斯信息准则(BIC)来迭代不同p与q的取值下的估计，综合考虑准则给出的p与q来达到最优定阶。对银杏液流进行AIC 与BIC 双重准则下的定阶，结果如表2 所示。

表2 AIC 与BIC 准则定阶

图7 银杏液流自相关系数图与偏自相关系数图

根据AIC 与BIC准则，确定模型ARIMA(2,1,3)，接下来对进行检验，在指数平滑模型下，关注ARIMA 模型的残差是否平均值为零且方差为常数的正态分布，同时还需关注连续残差是否自相关。最后进行ARIMA 模型的检验，如图8 所示。

图8 ARIMA 模型检验

从图8 可以看出，ARIMA(2,1,3)的残差接近正态分布，且相互独立。模型的DW 检验值为1.953 2 接近2，即模型随机误差项不存在一阶自回归的问题，最终认为ARIMA(2,1,3)符合建模要求。

3.2.3 模型预测

原银杏液流与分量序列信号前75%设置为训练集，后25%设置测试集，分别对mra1 与res 进行上述建模流程，再将两分量经过ARIMA 模型预测得到值进行叠加重构，与原序列在测试集上的表现进行指标对比，如表3所示。其中，MSE、MAE、MAPE和R2分别表示均方误差、平均绝对误差、平均绝对百分比误差和决定系数。

表3 ARIMA 与EWT-ARIMA 模型对比

两种模型的预测结果的可视化对比如图9 和图10所示，可以看出ARIMA 模型预测结果与真实值在明显的“滞后”，导致预测指标不理想。

图9 ARIMA 模型预测值与真实值对比图

本文提出的组合模型得益于分量信号平稳的震荡模式，其中res 余项信号几乎完美表征了原银杏序列的变化趋势，在4 个模型评价指标上的表现相较于单一ARIMA 模型均有较大提升。

3.3 银杏液流与环境因子关联分析

多种环境因子会对银杏液流造成重要影响[18-20]，为了厘清液流与环境因子的内在关联，本文对环境因子与银杏液流在0.05 显著性水平下进行了Granger 因果检验。经检验8 月29 日～9 月4 日期间的气温、湿度、水汽压亏缺的P 值无限接近0，其中土壤湿度的P 值为0.048，均通过了显著性水平下的Granger 因果检验。根据Granger 因果检验思想，可以认为这4 种环境因子对银杏液流的预测是有益的。值得注意的是Granger 因果检验仅能推断出包含环境因子对液流的预测是否有益，无法得出环境因子和液流是否存在因果关系的结论。

接下来，本文借助无需模型假设的传递熵[21]对4 组环境因子与银杏液流内在关联进行了定量的分析。

温度、湿度、土壤湿度、水汽压亏缺四者的转移熵变化趋势如图11 所示。当转移熵值大于零时，存在环境因子对银杏液流的因果关系，转移熵数值越大表明因果关系越强[22]。

图11 4 种环境因子的转移熵图

根据转移熵，认为温度和湿度是对银杏液流因果关系最强的两个因子。四者转移熵最小值都出现在最开始的时滞，在时滞整体上存在多个转折。湿度与水汽压亏缺的极大值均出现在200 min 时滞（分别为0.284 046和0.273 419），除了土壤湿度外其余因子都在时滞开始阶段有快速上升的趋势，这可能与银杏对环境因子的响应能力以及其体内的储水状况有关[23]。

其中土壤湿度波动性较强，但其时滞内的转移熵在四者当中最低，且出现了多次负值，认为土壤湿度对银杏液流因果关系较弱，这可能是由于校园内有专人负责养护，土壤水分条件保持较好所致，与马靖涵等[6]在城市园林银杏液流变化研究中所得出的结论类似。