王子祥，李顏娥，武斌，徐達(dá)宇，吳斌

（1.浙江農(nóng)林大學(xué) 數(shù)學(xué)與計(jì)算機(jī)學(xué)院，浙江 杭州 311300；2.浙江省林業(yè)智能監(jiān)測(cè)與信息技術(shù)實(shí)驗(yàn)室，浙江 杭州 311300；3.林業(yè)感知技術(shù)與智能裝備國(guó)家林業(yè)局重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300）

0 引言

進(jìn)入二十一世紀(jì)，人類正面臨氣候惡化這一全球性問(wèn)題。碳排放造成的全球變暖會(huì)對(duì)人類和生態(tài)系統(tǒng)造成影響嚴(yán)重且不可逆轉(zhuǎn)的危害[1]，利用植物生態(tài)系這一巨大天然碳庫(kù)去實(shí)現(xiàn)固碳目標(biāo)被認(rèn)為是經(jīng)濟(jì)有效的方法[2]。同時(shí)植物也是城市公園內(nèi)固碳效益的主體，在常見(jiàn)的城市公園樹(shù)木中，銀杏的固碳能力明顯高于其他植物，屬于第一梯隊(duì)[3]。除此之外，銀杏還是一種的珍貴的中藥材和觀賞植物，具有廣泛的生態(tài)、經(jīng)濟(jì)和社會(huì)價(jià)值[4-5]。為了生態(tài)環(huán)境的可持續(xù)發(fā)展，尋找一種準(zhǔn)確評(píng)估樹(shù)木蒸騰耗水的穩(wěn)健方法十分有必要[6]，而植物蒸騰耗水量的99%以上來(lái)自莖干液流，精確測(cè)定莖干液流量能夠反映林木的蒸騰耗水量[7]。液流法也成為了近年來(lái)莖分和林分尺度蒸騰耗水研究的熱點(diǎn)方法[8]。

基于樹(shù)木液流的蒸騰評(píng)估，需要克服液流這類自然序列普遍存在的非線性與高隨機(jī)性。因此本文提出一種基于經(jīng)驗(yàn)小波變換(Empirical Wavelet Transform,EWT)和ARIMA(Auto Regressive Integrated Moving Average)組合的銀杏液流預(yù)測(cè)模型。銀杏液流序列先經(jīng)過(guò)經(jīng)驗(yàn)小波分解,得到平穩(wěn)的多分辨分析分量(Multi Resolution Analysis,MRA)，針對(duì)各分量信號(hào)分別構(gòu)建ARIMA 模型進(jìn)行預(yù)測(cè)并重構(gòu)集成為目標(biāo)預(yù)測(cè)信號(hào)，并與單一ARIMA 模型進(jìn)行對(duì)比分析。此外，有許多研究證明樹(shù)木液流與環(huán)境因子有著密切的關(guān)系。但以往的研究多采用線性分析工具，在理論上有局限性，因此本文借助無(wú)需模型假設(shè)的轉(zhuǎn)移熵（Transfer Entropy,TE）對(duì)環(huán)境因子和銀杏液流時(shí)滯內(nèi)的因果關(guān)系進(jìn)行了探究。

1 相關(guān)工作

發(fā)現(xiàn)關(guān)聯(lián)與變量選擇在科學(xué)實(shí)踐中至關(guān)重要，但目前針對(duì)液流與環(huán)境因子關(guān)系的模型多采用單變量分析、多元線性回歸[9-10]，特定的分布假設(shè)使得這類方法在理論上存在不足。隨著技術(shù)的發(fā)展，許多智能算法被廣泛應(yīng)用到具有復(fù)雜模式的序列中，獲得了不錯(cuò)的效果，EWT 技術(shù)[11-13]就是其中之一。常見(jiàn)做法是將復(fù)雜序列視作信號(hào)，經(jīng)EWT 分解為若干穩(wěn)定分量，再依據(jù)各自的數(shù)據(jù)特征進(jìn)行重構(gòu)，從而取得一個(gè)優(yōu)良的預(yù)測(cè)效果。

2 基于EWT-ARIMA 組合模型的銀杏液流預(yù)測(cè)

2.1 經(jīng)驗(yàn)小波變換EWT

2.1.1 從傅里葉變換到經(jīng)驗(yàn)小波變換

傅里葉變換是一個(gè)非常強(qiáng)大而靈活的方法，其原理說(shuō)明由正弦基函數(shù)組成的合成信號(hào)可以分解為有限個(gè)正弦波或余弦波的疊加，且這些分量信號(hào)的頻率是原始信號(hào)頻率f0的整數(shù)倍。

式中,f0為基頻，A0/2 代表直流系數(shù)，An代表幅度，?n被稱為相位。根據(jù)振幅和相位等信息可以利用反變換從而恢復(fù)信號(hào)的最初波形。

對(duì)一個(gè)正弦信號(hào)和一個(gè)非線性信號(hào)做傅里葉變換，如圖1 所示。

圖1 正弦信號(hào)與非線性信號(hào)傅里葉圖

根據(jù)圖1 可以顯示出非線性信號(hào)比正弦信號(hào)具備更寬的波谷和更尖銳的波峰。

根據(jù)圖2 偉爾奇周期圖，正弦信號(hào)和非線性信號(hào)都具備5 Hz 上的高功率，符合二者的定義。但是在10 Hz處，非線性信號(hào)還呈現(xiàn)出另外一個(gè)峰值，這是因?yàn)楦道锶~變換只能基于固定的線性基函數(shù)，導(dǎo)致其無(wú)法表達(dá)信號(hào)復(fù)雜的震蕩模式。

圖2 正弦信號(hào)與非線性信號(hào)偉爾奇周期圖

經(jīng)驗(yàn)小波變換（EWT）是Gilles 提出的一種分解技術(shù)[14]，它結(jié)合了小波的形態(tài)優(yōu)勢(shì)和經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）的自適應(yīng)性?？梢杂行p少模態(tài)混疊和端點(diǎn)效應(yīng)，解決原始序列模型因非線性和高隨機(jī)而導(dǎo)致的預(yù)測(cè)精度低的問(wèn)題。

2.1.2 經(jīng)驗(yàn)小波變換原理

經(jīng)驗(yàn)小波變換原理大致可以總結(jié)如下：首先對(duì)序列進(jìn)行傅里葉變換，得到傅里葉頻譜并將頻率規(guī)范至[0,π]，找到并以降序排列頻譜中的N個(gè)局部極大值所對(duì)應(yīng)的頻率ωn（n=1，2，3，…，N-1）。將兩個(gè)連續(xù)局部極大值的中點(diǎn)設(shè)置為支撐邊界，然后執(zhí)行頻帶分割：

其中，ωi和ωi+1為兩個(gè)連續(xù)局部極大值的頻率，支撐邊界集合Ω={Ωi}i=1,2,…,N-1。

頻帶劃分過(guò)程如圖3 所示。其中τn代表過(guò)渡空間段一半的長(zhǎng)度。

圖3 EWT 頻帶分割示意圖

然后基于支撐邊界所形成過(guò)渡空間設(shè)計(jì)自適應(yīng)小波濾波器組。最后將經(jīng)濾波后的模態(tài)轉(zhuǎn)換成尺度函數(shù)與小波函數(shù)，以獲得信號(hào)分量。

2.2 ARIMA 預(yù)測(cè)模型

ARIMA 模型(Autoregressive Integrated Moving Average model)是一種用于時(shí)間序列預(yù)測(cè)的統(tǒng)計(jì)方法，它可以捕獲時(shí)間序列中的不同特征，如趨勢(shì)、季節(jié)性、周期性等，曾極大地促進(jìn)了趨勢(shì)預(yù)測(cè)研究[15-17]。ARIMA 模型由三部分組成：自回歸(AR)、差分(I)和移動(dòng)平均(MA)。ARIMA 模型有3 個(gè)參數(shù)(p,d,q)，分別代表自回歸項(xiàng)數(shù)、差分次數(shù)和移動(dòng)平均項(xiàng)數(shù)。

ARIMA 預(yù)測(cè)原理：ARIMA 模型預(yù)測(cè)先通過(guò)ADF 方法來(lái)檢驗(yàn)時(shí)間序列是否具備平穩(wěn)特征，若不平穩(wěn)則進(jìn)行d階差分直至差分后轉(zhuǎn)化為平穩(wěn)的時(shí)間序列；然后根據(jù)自相關(guān)系數(shù)圖(ACF)和偏自相關(guān)系數(shù)圖(PACF)或者基于AIC 與BIC 準(zhǔn)則選取模型參數(shù)，最后確定模型參數(shù)并進(jìn)行檢驗(yàn)和預(yù)測(cè)。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)數(shù)據(jù)集

本文采用了一份來(lái)校園內(nèi)銀杏樹(shù)的液流數(shù)據(jù)，數(shù)據(jù)由單針熱脈沖樹(shù)干液流監(jiān)測(cè)儀分別于2021 年8 月29 日～9 月4 日期間，在一棵銀杏樹(shù)上進(jìn)行數(shù)據(jù)的采集而成。試驗(yàn)地點(diǎn)位于浙江農(nóng)林大學(xué)東湖校區(qū)(30.24°～30.26°N，119.72°～119.73°W)，在監(jiān)測(cè)樹(shù)干液流的同時(shí)，同步監(jiān)測(cè)大氣溫度(Air Temperature,TA)、水汽壓虧缺（Vapor Pressure Deficit,VPD）、大氣濕度（Air Relative Humidity,RH）、土壤濕度（Soil Moisture,SM）等氣象因子。在樹(shù)木莖干的莖向直線上將熱源探針插入，然后分別在其上、下部位約6 mm 和8 mm 處插入兩探針,開(kāi)始進(jìn)行數(shù)據(jù)監(jiān)測(cè)，如圖4 所示。

圖4 樹(shù)干液流監(jiān)測(cè)儀以及探針距離示意圖

8 月29 日～9 月4 日監(jiān)測(cè)期間，屬于夏季的晴天氣溫較高。繪制銀杏液流變化趨勢(shì)如圖5 所示。

圖5 銀杏液流序列變化趨勢(shì)圖

其中,圖5（a)代表銀杏液流完整的變化趨勢(shì)，圖5(b)是選取了圖5(a)其中8 月30 日數(shù)據(jù)放大后的效果。通過(guò)觀察圖5 銀杏液流曲線可知，銀杏液流呈現(xiàn)周期性變化，日出后開(kāi)始快速上升，直到中午前后達(dá)到峰值。

接下來(lái)對(duì)完整的銀杏液流序列進(jìn)行經(jīng)驗(yàn)小波變換，經(jīng)過(guò)EWT 分解，得到兩組多分辨分析分量信號(hào)（mra1 和res）以及對(duì)應(yīng)的頻譜，發(fā)現(xiàn)經(jīng)分解后mra1 分量信號(hào)大致以零軸為上下對(duì)稱，序列較平穩(wěn)，對(duì)應(yīng)的頻譜也有一致的變化趨勢(shì)。同時(shí)，res 分量信號(hào)被分解得比較徹底，分解后的曲線十分平滑，為后續(xù)的建模預(yù)測(cè)建立了基礎(chǔ)。分解效果圖6 所示。

圖6 銀杏液流序列EWT 分解與mra 分量頻譜對(duì)照?qǐng)D

3.2 ARIMA 模型預(yù)測(cè)

3.2.1 模型平穩(wěn)性檢驗(yàn)

ARIMA 模型要求時(shí)序數(shù)據(jù)是穩(wěn)定的,接下來(lái)對(duì)銀杏液流和經(jīng)分解得到的兩組分辨率分析分量信號(hào)進(jìn)行ADF(Augmented Dickey-Fuller test)檢驗(yàn)，ADF 通過(guò)判斷序列是否存在單位根來(lái)判斷序列的平穩(wěn)性。ADF 單位根檢驗(yàn)結(jié)果如表1 所示。

表1 原序列與mra 分量ADF 單位根檢驗(yàn)

3 組序列P 值分別為0.0245、0、0.003，均小于顯著性水平且T 檢驗(yàn)值均小于3 個(gè)臨界值，表明3 組序列所檢驗(yàn)的參數(shù)落在拒絕域，拒絕原假設(shè)，認(rèn)為不存在單位根，故均為平穩(wěn)序列，可以通過(guò)ARIMA 建模。

3.2.2 模型參數(shù)定階與檢驗(yàn)

銀杏液流預(yù)測(cè)首先要確定參數(shù)p與q，通過(guò)銀杏液流的自相關(guān)系數(shù)圖（ACF）和偏自相關(guān)系數(shù)圖（PACF）所呈現(xiàn)的拖尾或截尾特征來(lái)進(jìn)行定階。

根據(jù)圖7，可判斷ACF 與PACF 屬于拖尾，p與q可定階在3 左右，但這種定階方式存在一定主觀性。除此之外還可以使用赤池信息準(zhǔn)則(AIC)與貝葉斯信息準(zhǔn)則(BIC)來(lái)迭代不同p與q的取值下的估計(jì)，綜合考慮準(zhǔn)則給出的p與q來(lái)達(dá)到最優(yōu)定階。對(duì)銀杏液流進(jìn)行AIC 與BIC 雙重準(zhǔn)則下的定階，結(jié)果如表2 所示。

表2 AIC 與BIC 準(zhǔn)則定階

圖7 銀杏液流自相關(guān)系數(shù)圖與偏自相關(guān)系數(shù)圖

根 據(jù)AIC 與BIC準(zhǔn)則，確定模型ARIMA(2,1,3)，接下來(lái)對(duì)進(jìn)行檢驗(yàn)，在指數(shù)平滑模型下，關(guān)注ARIMA 模型的殘差是否平均值為零且方差為常數(shù)的正態(tài)分布，同時(shí)還需關(guān)注連續(xù)殘差是否自相關(guān)。最后進(jìn)行ARIMA 模型的檢驗(yàn)，如圖8 所示。

圖8 ARIMA 模型檢驗(yàn)

從圖8 可以看出，ARIMA(2,1,3)的殘差接近正態(tài)分布，且相互獨(dú)立。模型的DW 檢驗(yàn)值為1.953 2 接近2，即模型隨機(jī)誤差項(xiàng)不存在一階自回歸的問(wèn)題，最終認(rèn)為ARIMA(2,1,3)符合建模要求。

3.2.3 模型預(yù)測(cè)

原銀杏液流與分量序列信號(hào)前75%設(shè)置為訓(xùn)練集，后25%設(shè)置測(cè)試集，分別對(duì)mra1 與res 進(jìn)行上述建模流程，再將兩分量經(jīng)過(guò)ARIMA 模型預(yù)測(cè)得到值進(jìn)行疊加重構(gòu)，與原序列在測(cè)試集上的表現(xiàn)進(jìn)行指標(biāo)對(duì)比，如表3所示。其中，MSE、MAE、MAPE和R2分別表示均方誤差、平均絕對(duì)誤差、平均絕對(duì)百分比誤差和決定系數(shù)。

表3 ARIMA 與EWT-ARIMA 模型對(duì)比

兩種模型的預(yù)測(cè)結(jié)果的可視化對(duì)比如圖9 和圖10所示，可以看出ARIMA 模型預(yù)測(cè)結(jié)果與真實(shí)值在明顯的“滯后”，導(dǎo)致預(yù)測(cè)指標(biāo)不理想。

圖9 ARIMA 模型預(yù)測(cè)值與真實(shí)值對(duì)比圖

本文提出的組合模型得益于分量信號(hào)平穩(wěn)的震蕩模式，其中res 余項(xiàng)信號(hào)幾乎完美表征了原銀杏序列的變化趨勢(shì)，在4 個(gè)模型評(píng)價(jià)指標(biāo)上的表現(xiàn)相較于單一ARIMA 模型均有較大提升。

3.3 銀杏液流與環(huán)境因子關(guān)聯(lián)分析

多種環(huán)境因子會(huì)對(duì)銀杏液流造成重要影響[18-20]，為了厘清液流與環(huán)境因子的內(nèi)在關(guān)聯(lián)，本文對(duì)環(huán)境因子與銀杏液流在0.05 顯著性水平下進(jìn)行了Granger 因果檢驗(yàn)。經(jīng)檢驗(yàn)8 月29 日～9 月4 日期間的氣溫、濕度、水汽壓虧缺的P 值無(wú)限接近0，其中土壤濕度的P 值為0.048，均通過(guò)了顯著性水平下的Granger 因果檢驗(yàn)。根據(jù)Granger 因果檢驗(yàn)思想，可以認(rèn)為這4 種環(huán)境因子對(duì)銀杏液流的預(yù)測(cè)是有益的。值得注意的是Granger 因果檢驗(yàn)僅能推斷出包含環(huán)境因子對(duì)液流的預(yù)測(cè)是否有益，無(wú)法得出環(huán)境因子和液流是否存在因果關(guān)系的結(jié)論。

接下來(lái)，本文借助無(wú)需模型假設(shè)的傳遞熵[21]對(duì)4 組環(huán)境因子與銀杏液流內(nèi)在關(guān)聯(lián)進(jìn)行了定量的分析。

溫度、濕度、土壤濕度、水汽壓虧缺四者的轉(zhuǎn)移熵變化趨勢(shì)如圖11 所示。當(dāng)轉(zhuǎn)移熵值大于零時(shí)，存在環(huán)境因子對(duì)銀杏液流的因果關(guān)系，轉(zhuǎn)移熵?cái)?shù)值越大表明因果關(guān)系越強(qiáng)[22]。

圖11 4 種環(huán)境因子的轉(zhuǎn)移熵圖

根據(jù)轉(zhuǎn)移熵，認(rèn)為溫度和濕度是對(duì)銀杏液流因果關(guān)系最強(qiáng)的兩個(gè)因子。四者轉(zhuǎn)移熵最小值都出現(xiàn)在最開(kāi)始的時(shí)滯，在時(shí)滯整體上存在多個(gè)轉(zhuǎn)折。濕度與水汽壓虧缺的極大值均出現(xiàn)在200 min 時(shí)滯（分別為0.284 046和0.273 419），除了土壤濕度外其余因子都在時(shí)滯開(kāi)始階段有快速上升的趨勢(shì)，這可能與銀杏對(duì)環(huán)境因子的響應(yīng)能力以及其體內(nèi)的儲(chǔ)水狀況有關(guān)[23]。

其中土壤濕度波動(dòng)性較強(qiáng)，但其時(shí)滯內(nèi)的轉(zhuǎn)移熵在四者當(dāng)中最低，且出現(xiàn)了多次負(fù)值，認(rèn)為土壤濕度對(duì)銀杏液流因果關(guān)系較弱，這可能是由于校園內(nèi)有專人負(fù)責(zé)養(yǎng)護(hù)，土壤水分條件保持較好所致，與馬靖涵等[6]在城市園林銀杏液流變化研究中所得出的結(jié)論類似。

4 結(jié)論

本文基于樹(shù)干液流監(jiān)測(cè)儀獲取的銀杏液流以及環(huán)境因子數(shù)據(jù)，借助經(jīng)驗(yàn)小波變換算法對(duì)銀杏液流進(jìn)行充分分解，獲得兩組平穩(wěn)的分量信號(hào)，分別使用ARIMA 統(tǒng)計(jì)模型對(duì)分量信號(hào)進(jìn)行建模，最后將預(yù)測(cè)結(jié)果疊加，取得最終良好的預(yù)測(cè)效果。并以MSE、MAE、MAPE、R2為預(yù)測(cè)精度的評(píng)價(jià)指標(biāo)。

樹(shù)木液流受到多重環(huán)境因子綜合影響，本文借助轉(zhuǎn)移熵對(duì)4 種環(huán)境因子與銀杏液流時(shí)滯內(nèi)的因果關(guān)系進(jìn)行了分析，轉(zhuǎn)移熵方法相較于Granger 因果檢驗(yàn)有無(wú)需模型假設(shè)的優(yōu)勢(shì)且能給出定量的因果關(guān)系。這對(duì)進(jìn)一步厘清樹(shù)木水分關(guān)系與響應(yīng)環(huán)境變化機(jī)制，有一定的積極意義。