景天忠，劉麗萍，謝雨龍

(東北林業(yè)大學林學院，哈爾濱 150040)

物候模型(phenology model)預測昆蟲發(fā)育事件的時間。物候模型已成為種群動態(tài)、物種分布和進化動態(tài)的關(guān)鍵成分(Chuine and Régnière, 2017)。在害蟲管理生產(chǎn)實踐中，需要預測出害蟲某一發(fā)育事件出現(xiàn)的時間，以便確定害蟲防治的最佳時間，稱為發(fā)生期預測，其實就是對昆蟲物候?qū)W的預測。美國國家物候?qū)W網(wǎng)絡(luò)中心(U.S. National Phenology Network，USA-NPN)在2018年開發(fā)了一套可用于科研和生產(chǎn)管理的物候預測圖(Pheno Foescast map)(Crimminsetal., 2020)，包括實時圖和在與管理相關(guān)的時間和空間分辨率上短期預測害蟲的活動。Nietschkeetal.(2007)建立了一個昆蟲發(fā)育數(shù)據(jù)庫(Insect Development Database, IDD)，包含405種昆蟲的基礎(chǔ)發(fā)育溫度閾值和310種昆蟲產(chǎn)卵和至成蟲羽化的累計日度。這些昆蟲分屬11個目和89個科。加州大學州級IPM項目(http://www.ipm.ucdavis.edu/MODELS)提供了許多昆蟲的物候信息。美國還建立了一個在線的物候和日度模型，用于農(nóng)業(yè)和害蟲管理決策(http://uspest.org/cgi-bin/ddmodel.us)。

對變溫動物而言，完成一個特定的發(fā)育階段需要累積一個種特異的溫度常數(shù)(thermal constant)，這就是最早由Réaumur在1735年提出的有效積溫法則(參見Ikemoto and Takai, 2000; Damos and Savopoulou-Soultani, 2012)，可用公式表示：

K=D(T-T0)

(1)

其中K為溫度常數(shù)，即有效積溫(effective cumulative temperature)；T為環(huán)境溫度；T0為發(fā)育起點溫度；D為發(fā)育歷期(development time)，即昆蟲完成一定發(fā)育階段所經(jīng)歷的時間，常用τ(T)來表示。而發(fā)育歷期的倒數(shù)，便是發(fā)育速率(development rate)，常用r(T)來表示。要對昆蟲在變溫條件下的發(fā)育進行建模，可采用積分的方法，即

(2)

其中，a為某一階段的生理年齡(physiological age)，表示完成某一階段發(fā)育的成數(shù)，取值在0～1。0表示發(fā)育的開始，1表示發(fā)育完成。A為由溫度響應(yīng)函數(shù)r(T,A)的參數(shù)值組成的向量。昆蟲完成發(fā)育所需的熱量也稱為生理學歷期(physiological time)。

有效積溫法則是對昆蟲物候研究的基礎(chǔ)，因而昆蟲物候模型總是圍繞溫度、發(fā)育速率或發(fā)育歷期這些變量。在物候模型的發(fā)展過程中出現(xiàn)了熱性能曲線、生物物理模型、基于概率的模型、分布時滯模型、發(fā)育進度曲線、物候匹配和物候變遷等模型。除熱性能曲線和發(fā)育進度曲線外，其余模型在中文文獻中鮮有提及。因而很有必要對這些模型做一較為系統(tǒng)的介紹，以備借鑒。

1 熱性能曲線

對式(1)兩邊取倒數(shù)，得

(3)

稱為線性日度模型(linear degree-day model)或X截距模型(x-intercept method)，是一個標準的一次函數(shù)模型(y=1/Donx=T)，可用最小二乘法求出參數(shù)。1/D為發(fā)育速率。

有學者指出，上述公式存在三個不足：一是利用上述公式很難得到最適溫度范圍；二是數(shù)據(jù)點的權(quán)重在直線的上部和下部之間不成比例，且沿X軸并不同質(zhì)。溫度越低，擔負的權(quán)重越不相稱(越少)；三是對大多數(shù)數(shù)據(jù)，變量X的誤差被忽略了(Ikemoto and Takai, 2000)。因而提出一了個新的線性模型：

(DT)=K+T0D

(4)

很明顯，利用有效積溫法則，可預測害蟲和天敵的發(fā)生期，預測害蟲和天敵的地理分布及在某地一年發(fā)生的代數(shù)。但要注意，昆蟲生活場所的氣溫和百葉箱所代表的大氣氣溫有差別。室內(nèi)所計算的發(fā)育起點和有效積溫常常是在恒溫條件下所得的數(shù)據(jù)，與自然變溫狀態(tài)下不同，變溫下昆蟲的發(fā)育常較在恒溫下的為快(Colinetetal., 2015)。積溫法則無法應(yīng)用于某些有滯育階段的昆蟲世代數(shù)的計算。昆蟲的發(fā)育不僅僅受溫度一個因素的影響，還受食物、濕度等其他因素的影響。還有，昆蟲的發(fā)育速率與溫度之間通常是非線性的關(guān)系(Rebaudo and Rabhi, 2018)。因此，也有不少學者提出了多達幾十種非線性模型。這些模型包括Bayoh、Beta、Brière1、Brière2、Damos、Gaussian、HarcourtYee、Kontodimas、Lactin-1、Lactin-2、Logan1、Logan2、Performance、Wang-Lan-Ding、Ratkowsky、Taylor、Weibull、改進的Gaussian和指數(shù)方式改進的Gaussian等，它們通常以溫度為橫軸、發(fā)育速率為縱軸，稱為熱性能曲線(thermal performance curves/functions/equations/models)(Shi and Ge, 2010; Rebaudo and Rabhi, 2018)。R語言的devRate包可擬合這些曲線。Geng and Jung(2018)建立了金紋細蛾P(guān)hyllonorycterringoniella越冬蛹發(fā)育速率的Logan-6模型。這些模型還用于研究秘夜蛾Mythimnasequax的發(fā)育速率(Jacquesetal., 2019)、預測棉鈴蟲Helicoverpaarmigera在西紅柿Lycopersiconesculentum上的物候(Dalal and Arora, 2019)，以及螽斯Ruspoliadifferens當前和未來的分布(Leonardetal., 2021)。Shi and Ge(2010)利用小菜蛾P(guān)lutellaxylostella和B型煙粉虱Bemisiatabaci為研究對象，以AIC(Akaike information criterion)、BIC(Bayesian information criterion)等方法比較了12種非線性模型，結(jié)果表明Performance，Brière1和Brière2三種模型最適合應(yīng)用于溫度依賴的發(fā)育速率。Régnièreetal.(2012)開發(fā)了一個經(jīng)驗?zāi)Ｐ?，可對近上、下臨界溫度的熱響應(yīng)更好的描述。Anarcaetal.(2018)將寄主質(zhì)量包含在模型中來提升模型。

2 生物物理模型

有人根據(jù)范霍夫定律(Van’t Hoff’s law)提出了生物物理模型-Sharpe-Schoolfield模型(參見Damos and Savopoulou-Soultani, 2012)。生物物理模型將變溫動物的溫度依賴的發(fā)育看作是一個宏觀的酶促反應(yīng)。基于Arrhenius和Eyring關(guān)于化學反應(yīng)速度的理論，該模型最初由Johnson & Lewin、Hultin、Sharpe和DeMichele等人建立，后由Schoolfield、Sharpe和Magnuson改進(Schoolfieldetal., 1981)。溫度為T(絕對溫度，K)時，平均發(fā)育速率r(T)(time-1)為：

(5)

Régnièreetal.(2012)開發(fā)了一個經(jīng)驗函數(shù)，可更好的描述臨界溫度的熱響應(yīng)：

(6)

其中Tb≤T≤Tm，6個參數(shù)與式(5)類似。式(6)與(5)在大部分溫度范圍內(nèi)重合，但在臨界溫度周圍有很大的不同。

為處理卵在高溫下由于高死亡率而無法計算ΔHH的問題，Salisetal.(2016)采用縮減版的Sharpe-Schoolfield模型：

(7)

并開發(fā)了一個發(fā)育階段依賴的模型(developmental stage-dependent model, DSD)：

(8)

其中，ST為用于校正溫度對發(fā)育速率影響的縮放因子；Dth為發(fā)育速率不受發(fā)育階段影響時的發(fā)育閾值；D為發(fā)育階段，即從產(chǎn)卵(D=0)至卵孵化(D=1)的發(fā)育進程。

Crespo-Pérezetal.(2015)利用式(5)研究3種馬鈴薯塊莖蛾P(guān)hthorimaeaoperculella的日發(fā)育速率，Jacquesetal.(2019)研究了Spodopteralitura在印度落花生上的定殖和存活。生物物理模型還用于研究個體大小與溫度的關(guān)系(de Jong, 2010)。

3 基于概率的模型

將昆蟲飼養(yǎng)在完全一樣的條件下，不同個體的發(fā)育歷期不盡相同。這種變異性可能來自于一個種群不同個體間的代謝和遺傳方面的差異。諸多研究表明，這種變異性是不對稱分布的，即右偏。根據(jù)式(5)，影響發(fā)育速率的參數(shù)中，酶是唯一的受遺傳變異影響的生物學參數(shù)，且研究表明生物的一些遺傳特征如頭殼寬度、長度、大小、智力等都是對稱的分布。因此，Sharpeetal.(1977)認為限速酶的濃度也應(yīng)該是對稱分布。在這些假設(shè)下，式(5)可改寫為

R(T)=r(T)ε

(9)

其中ε為服從正態(tài)分布的隨機變量。

這樣就可利用正態(tài)分布的密度函數(shù)或者二次項分布密度函數(shù)來描述發(fā)育速率的變異(橫軸為發(fā)育速率，縱軸為概率密度)，即

(10)

fQ(r)=αr2+βr+γ0

(11)

發(fā)育歷期的密度函數(shù)(橫軸為發(fā)育歷期，縱軸為概率密度)為

(12)

(13)

這類模型稱為基于群組(cohort-based)的模型，或者稱為發(fā)育歷期的分布模型。群組為一個特定生命階段(life stage)的所有個體，這個階段中的個體具有幾乎一樣的時間年齡(chronological age)。

在發(fā)育速率變異是正態(tài)分布的情況下，可以把昆蟲種群中16%、50%和84%左右的個體進入某一發(fā)育事件的日期作為生產(chǎn)中劃分害蟲發(fā)生始盛期、高蜂期和盛末期的數(shù)量標準，對應(yīng)于正態(tài)分布-σ、μ和+σ分位數(shù)。

設(shè)τi(T)為單個昆蟲在恒溫T下完成其第i個生命階段所需的時間，T(t)為t時刻的溫度、a(t)為在t時刻一頭昆蟲完成某一生命階段的成數(shù)，a的取值范圍在0到1之間，可以看作是昆蟲的生理年齡，0表示起始發(fā)育，1表示發(fā)育完成。那么在溫度T下，昆蟲完成某一生命階段Δa部分所需的天數(shù)為Δt，則Δt=τi(T) Δa。那么，一頭昆蟲從t0開始的階段i的發(fā)育可以用下式來描述：

(14)

(15)

一頭昆蟲，從時刻t開始，完成階段i的發(fā)育所需的時間為Γi(1,t)，且對所有的時刻t，Γi(0,t)=t。當Γi已知時，這頭昆蟲的物候?qū)W就可以確定了(Yurk and Powell, 2010)。這類模型稱為基于個體的模型(individual-based models)。

假設(shè)發(fā)育歷期不變，上面這個個體物候?qū)W模型可以擴展成種群模型。令p(a,t)為t時刻具有年齡a的個體的密度函數(shù)，則由完全相同個體(不存在變異)組成的種群的發(fā)育可用一個平流方程(advection equation)來描述：

(16)

在上述模型中引入一個能夠縮放基礎(chǔ)發(fā)育歷期曲線的表型參數(shù)α，就可解釋發(fā)育歷期中的持續(xù)變異。這個表型在種群內(nèi)具有相同的分布。在這種情況下，式(16)變成為一個α依賴的平流方程

(17)

如果種群的每個個體都具有相同的表型α=1，則(17)式變?yōu)槭?16)。

在式(17)中加入一個擴散項(diffusion term)來將由隨機效應(yīng)引起的發(fā)育變異加入到表型模型中，稱為Fokker-Planck發(fā)育方程(Fokker-Planck development equation)：

(18)

如果種群的每個具有表型α的個體具有相同的發(fā)育歷期(即變異v=0)，則式(18)變?yōu)槭?17)。

在恒溫條件下，發(fā)育歷期可由下式來預測：

(19)

即發(fā)育歷期服從均值為μτi(T)、方差為v+σ2τi(T)2的正態(tài)分布。

在變溫條件下，假設(shè)(1)描述不同恒溫條件下發(fā)育的PDF積分函數(shù)可度量不同溫度下的發(fā)育；(2)發(fā)育速率的均值和方差均是呈線性比例，則對每一個溫度T，f(r)的均值為r0，標準差為σ=cr0，應(yīng)用Sharpe模型，處于生命階段j的個體在t時刻的歷期(出現(xiàn)時間emergence times)分布為(Gilbertetal., 2004)：

(20)

其中，可用σ對r0的線性回歸來估計c。

通常對發(fā)育的觀測是以dt日為間隔的。如果定義個體與真平均值之間的差異δij服從對數(shù)正態(tài)分布，理論響應(yīng)值與處理組均值之間的擬合劣度υj服從正態(tài)分布，那么在觀測區(qū)間[t-dt,t]上，個體i在處理j下以恒溫Tj完成某一發(fā)育階段的概率pij為

(21)

當存在刪失數(shù)據(jù)時，無法估計一個處理的平均發(fā)育歷期。假定對該處理υ=1，一個個體完成其大于tij的發(fā)育的似然性為1-F(εij,A) ，則無論是否存在刪失數(shù)據(jù)，可用下式來表示似然性

Lij(σε,συ,A)=(1-dij)pij(σε,συ,A)+dij[1-F(εij,A)]

(22)

其中，d=0或1，0表示無刪失數(shù)據(jù)，1表示有刪失數(shù)據(jù)。然后可用最大似然法求出各參數(shù)。

當估計昆蟲發(fā)育的上下臨界溫度時，昆蟲易在近臨界溫度處死亡。Régnièreetal.(2012)提出了利用溫度轉(zhuǎn)換(temperature transfer)來解決此問題。先將昆蟲暴露于近臨界溫度(T1)下一段固定時間(t1)，既要短到避免過多死亡，又要長到保證明顯發(fā)育。然后將昆蟲轉(zhuǎn)移到另一個適合其發(fā)育的溫度(T2)下t2,i時間，使其完成發(fā)育。在這種情況下，當式(23)成立時，一個個體完成發(fā)育。

(23)

(24)

Gilbertetal.(2004)在McKendrick-von Foerster模型基礎(chǔ)上，開發(fā)出EvF模型(Extended von Foerster model)，可處理個體之間的變異(Gilbertetal., 2004)。

McKendrick-von Foerster模型為

(25)

其中，g(a,t,p(a,t))稱為種群中個體每時刻(per time)每年齡段(per age)的總增益(total gain)，或者負損失(negative loss)。Gilbertetal.(2004)將其擴展為

(26)

稱為EvF模型。其中，r(T(t))為發(fā)育速率，v(T(t))為發(fā)育變異。上式的解為

(27)

假定發(fā)育過程中沒有生產(chǎn)，且r服從正態(tài)分布，個體完成了一個發(fā)育階段(a=1)，則式(27)變?yōu)椋?/p>

(28)

其中，r0為平均發(fā)育速率。發(fā)育的概率即為發(fā)育速率r，可由平均發(fā)育速率r0來近似表示。

在變溫的環(huán)境條件下，處于發(fā)育階段j的個體在t時刻的種群分布函數(shù)為

(29)

Gilbert and Powell(2004)比較了式(29)和式(20)。這二者有著相似的計算復雜度，但EvF不比Sharpe模型費時。以高山松小蠹Dendroctonusponderosae為實例進行研究表明，EvF擬合曲線的R2為0.93，而Sharpe擬合曲線的R2為0.59(Gilbertetal., 2004)。

4 分布時滯模型

分布時滯(distributed delay)模型又稱為“悶子車”(‘boxcar’)模型(Manetsch, 1976; Vansickle, 1977; Gilbertetal., 2004)。設(shè)x(t)為一時間函數(shù)，則x(t-r)中的時滯r稱為離散時滯(discrete delay)，如Leslie矩陣就屬于離散模型，而形如

(30)

的式子稱為分布時滯，因為它反映了時滯x(t-z)的加權(quán)平均。上式實際上是卷積的表達式。

昆蟲的發(fā)育通常是不整齊的，例如在卵孵化為幼蟲階段，每天都有新幼蟲加入到幼蟲的行列中，這樣種群就具有了一個年齡分布，可以用一個代表其居留時間(residence time，dwell time)或存活時間(survival time)的隨機變量來描述?？捎肊rlang家族函數(shù)、正態(tài)分布或二次項分布概率密度函數(shù)來描述居留時間的分布，但Erlang分布最常用，其PDF為：

(31)

其中，k稱為形狀(shape)參數(shù)，r為速率(rate)參數(shù)，1/r稱為尺度(scale)參數(shù)。其均值μ=k/r，方差σ2=k/r2。它是gamma分布的一個特例，當k=1是時便簡化為指數(shù)分布。R的VGAM包括提供了擬合Erlang分布的函數(shù)。

將一個事件X或一種生物的生命歷程或生命階段劃分成k個離散的階段，或稱為悶子車、車廂(compartment)或相位(phase)，第k個悶子車完成所需的時間為k個獨立同分布變量的和。用Xi來表示任一悶子車，則X1,…,Xk組成X。當Xk完成了，X就完成了。假設(shè)悶子車Xi在t時刻有xi(t)個個體，每時刻有一部分個體(ri(t)，成數(shù))進入下一個悶子車。這里的ri稱為流速(flow rate)或轉(zhuǎn)移速率(transition rate)，也即所謂的分布時滯。假定沒有死亡，也沒有遷入遷出，所有個體都完成最后發(fā)育，則可用下面的微分方程來描述種群的變化

(32)

設(shè)a為一連續(xù)變量，代表年齡，歸一化后0≤a≤1，p為種群密度，則對具有Δa寬度的悶子車，有xi(t)?p(ai,t)Δa，這里ai為第i個悶子車的年齡，因此有ai=iΔa；xi(t)為第i個悶子車在t時刻的種群；p(ai,t)為第i個悶子車在t時刻的種群密度。對p(ai-1,t)進行Taylor展開，并假定整個生命階段的發(fā)育是一個常數(shù)(ri-1=ri=ri+1=r)，則式(32)變?yōu)?/p>

(33)

其中，rΔa代表了一個個體轉(zhuǎn)移到下一個悶子車的概率。解此方程，可得

(34)

在恒溫下，分布時滯模型實質(zhì)上von Foerster方程特征線法解的數(shù)值近似。一個重要區(qū)別是在分布時滯模型中每一個發(fā)育階段的相位數(shù)(number of phases)必須通過數(shù)據(jù)來估計(Gilbertetal., 2004)。分布時滯模型在變溫條件下很難與恒溫條件下的參數(shù)連接起來，因為為每一個溫度和每一個生命階段來指定悶子車數(shù)實際上不可能的。分布時滯模型將發(fā)育過程比做一群個體通過一系列悶子車，其輸出取決于悶子車之間的流速(rΔa)和悶子車數(shù)(k=1/Δa)。

Jianetal.(2007)開發(fā)一種用于時間可變分布時滯模型的算法，研究了不同溫度下銹赤扁谷盜Cryptolestesferrugineus老齡化的速率。

Scranton and Amarasekare(2017)使用時滯微分方程(delay differential equations, DDEs)預測了氣候變化背景下昆蟲物候的變化。該模型將昆蟲的生命周期分為幼年期(J)和成年期(A)兩個階段：

(35)

其中變量J和A分別代表幼年期和成年期的密度，函數(shù)B(T(t),A(t))和DX(T(t),X(t))描述的是密度和溫度對平均出生率和死亡率的聯(lián)合效應(yīng)；RA(t)是溫度依賴的補充率，即在t時刻的成熟個體，代表在τ(t)時間單位前出生，歷經(jīng)發(fā)育階段τ(t)后存活為SJ(t)的個體。這里τ(t)就是發(fā)育遲滯；溫度依賴的發(fā)育速率MJ(T(t))建立初始值為1/MJ(T(0))的發(fā)育遲滯τ(t)。結(jié)果表明，多度的峰值在環(huán)境溫度變暖的年分提前，但在季節(jié)性振蕩振幅增加的年分延遲。

5 發(fā)育進度曲線

在發(fā)生期預測工作中，還可以對進入某一發(fā)育階段個體的累積成數(shù)(cummulative proportion)對相應(yīng)的累計有效溫或發(fā)育歷期(儒略日期)進行擬合。如Geng and Jung (2018)建立了金紋細蛾P(guān)hyllonorycterringoniella成蟲累計羽化速率的Weibull模型、Damos and Savopoulou-Soultani(2010)采用Boltzman和Logistic回歸兩種模型研究了桃園中主要的鱗翅目有害生物復合體的發(fā)育進度。使用4參數(shù)的Weibull模型，Milosavljevietal.(2017)擬合了兩種金針蟲Limoniuscalifornicus和L.infuscatus在土壤中的累計捕獲率對有效積溫的曲線，D’Auriaetal.(2016)擬合了3種馬鈴薯害蟲的累計羽化率對有效積溫的曲線。以上這些例子都以累計有效積溫為自變量、以累計發(fā)育進度(成數(shù))為因變量進行曲線擬合。也有以儒略日期為自變量進行曲線擬合的。但由于年際間氣溫存在變異，使用累計有效積溫的模型比使用儒略日期的模型要穩(wěn)健得多。相比之下，使用儒略日期的模型的好處是數(shù)據(jù)容易獲取，在生產(chǎn)實踐中更容易應(yīng)用。

6 物候匹配模型

物候匹配模型用來描述昆蟲的物候與寄主物候的同步性。這類模型在授粉昆蟲中研究較多。

Hindleetal.(2015)采用一個統(tǒng)計模型來表示一種蝴蝶Melanargiagalathea在飛行期的某日t(儒略日)能見到的個體期望數(shù)量μ(t)：

(36)

(37)

其中，φ為方差參數(shù)，a=μ/φ，b=1/φ。假定φ獨立于年和地點。

該蝶的主要蜜源植物Centaureascabiosa的開花時間和花期在地點間和年間是不同的，可用下面的模型來模擬在儒略日x時的期望累計開花率：

(38)

(39)

其中，a=p/φ，b=(1-p)/φ。φ定義了植株間的變異程度。由于它在地點間、年間的微棲境之間是相似的，令其在所有模型中都以常量存在。以最簡約模型(AIC最小)來計算開花期，并將5%～95%的花開日子定義為花期。

最后，作者采用線性回歸來判定蝴蝶的飛行期以及植物的花期是否隨著年顯著變化。結(jié)果表明，蜜源植物的開花時間變晚了，但蝴蝶的飛行期并無明顯的變化趨勢。

7 物候變遷模型

物候變遷模型通常采用線性模型或加性模型將環(huán)境變量納入模型中來研究昆蟲物候?qū)夂蜃兓捻憫?yīng)。這類模型非常靈活，可依研究目的采用廣義線性模型、廣義加性模型(generalized additive models, GAM)、廣義線性混合效應(yīng)模型(generalized linear mixed models, GLMM)等。下面是近年來的幾個例子。為研究環(huán)境變量對昆蟲物候變化的驅(qū)動力，Searleetal.(2013)采取了一個兩相建模途徑。首先用GAM模型來擬合每個誘捕點每周捕獲的未產(chǎn)仔癭蚊DiarthronomyiachrysanthemiAhlberg的三周移動平均數(shù)。假定原始計數(shù)數(shù)據(jù)服從泊松分布，并采用log鏈接函數(shù)。使用三周移動平均數(shù)來減少每周誘捕過程中由于氣象條件的變異引起的多度隨時間的錯誤變化。這些模型用來確定兩種癭蚊每年每地點上不同世代高峰的多度峰值時刻。然后確定以下物候性狀：第1代或第2代的峰值時刻、與第1代或第2代聯(lián)系的多度、第1代與第2代之間的多度差別(取對數(shù)后相除)。接下來第二相分析使用GLMM模型將這些物候性狀與環(huán)境因子聯(lián)系起來。

Hollowayetal.(2018)使用熵來篩選環(huán)境變量，然后使用C5.0決策樹算法來預測氣象因子對蚜蟲揚飛驅(qū)動力。響應(yīng)變量為二元的觀測值：首飛(first flight)和未見揚飛(no flight)。由于是連續(xù)的監(jiān)測，將首飛前的7～105 d視為未見揚飛。環(huán)境變量有日氣溫、日氣壓、日北大西洋濤動、累積日度等。結(jié)果表明，與一般加性模型相比，決策樹模型對首飛預測的準確率提高了20%。通過熵篩選的變量與專家推薦的變量相比，預測的準確率也提高了3%～5%。

Stemkovskietal.(2020)使用GAM來預測蜜蜂的物候相(覓食起始、覓食高蜂和覓食終止)。對每一個物種/樣地/年度組合，以積日(day-of-year)為解釋變量，具有高斯分布族的多度為響應(yīng)變量擬合GAM模型。使用三次樣條進行平滑，使用交叉驗證來避免過擬合。建立如下混合效應(yīng)加性模型來研究物候的驅(qū)動力：

DOYphase～θclim+θtopo+θtrait+esite+esp

其中，DOYphase為每個物候相的積日估計值，θclim為氣候變量(雪融日期、夏季氣溫和夏季降水)，θtopo為拓撲變量(海拔和太陽入射角)，θtrait為物種性狀(體重、蜂巢位置、越冬蟲態(tài))，esite為樣地，esp為蜜蜂種類。θ項為固定效應(yīng)，e項為隨機效應(yīng)。

使用蜜蜂物候變化的天數(shù)除以雪融日期變化的天數(shù)得到的商來表示蜜蜂對雪融的響應(yīng)情況，正值表示提前、負值表示延遲。結(jié)果表明，蜜蜂的覓食起始物候相對氣候因素非常敏感，隨著雪融日期的提前而提前。

Gutiérrez and Wilson(2021)使用具有高斯誤差結(jié)構(gòu)的廣義最小二乘模型來為每種昆蟲的平均揚飛日期和首飛日期對氣候因子建模型。結(jié)果表明，對所有的蝴蝶種類，成蟲羽化前幾個月的氣溫與物候的聯(lián)系最強。所有的種類在較暖的年份都揚飛得較早，在飛行季節(jié)較早飛行的種類表現(xiàn)出對年度氣溫變異的敏感性最強。大多數(shù)種類的物候?qū)鉁氐拿舾行愿哂趯鉁乜臻g變異的敏感性。

Duchenneetal.(2020)建立如下模型來檢測多模態(tài)(multimodality)：

Yik=μ+ρ×latitudek+τ×longitudek+θ×altitudek+φi+Eik

(40)

其中，Yik是第i年第k個觀測的儒略日，μ為總平均數(shù)(截距)，ρ和τ分別表示緯度和經(jīng)度效應(yīng)，θ為海拔效應(yīng)，φi為年隨機效應(yīng)(年度)，Eik為誤差項。該模型的殘差就代表了去掉空間和海拔變異的采集日期。

為檢測這些殘差分布的多模態(tài)，可對分布進行平滑。起初從494種昆蟲中發(fā)現(xiàn)了些模式，然后對每一種昆蟲通過科學文獻或灰色文獻(grey literature)來驗證是否存在多模態(tài)的物候，結(jié)果表明288種昆蟲存在多模態(tài)的物候。對這些種類進行下一步混合高斯模型(Gaussian mixture-models)聚類分析，給每一條記錄分配一種物候模式。然后再一次運行與式(40)相似的線性混合效應(yīng)模型：

Yijk=μ+ρ×latitudek+τ×longitudek+θ×altitudek+βj+φi+Eik

(41)

其中βj稱為物候模式效應(yīng)。最后從2 473種昆蟲發(fā)現(xiàn)了2 473個單模態(tài)(unimodal)的物候模式。接下來估計平均揚飛日期(mean flight date, MFD)和揚飛期(flight period length, FPL)的變化。對每個物候模式，使用下面的模型來預測采集日期：

Yk=μ+(π+α×latitudek+δ×longitudek)×yeark+(ρ1+γ1×longitudek)×latitudek+

(42)

其中Yk為觀測點k的儒略日，π為時間效應(yīng)，其余同前。最后使用系統(tǒng)發(fā)育廣義最小二乘模型(phylogenetic generalized least squares model, PGLS)來檢驗季節(jié)性早熟以及物種的空間分布是否與物候的改變相聯(lián)系：

PSz=μ+α×MFDz+β×latitudez+δ×longitudez+Ez

(43)

其中，PSz為z種昆蟲在物候上的變化(即MFD或FPL的變化)，μ為總平均數(shù)(截距)，α為由近期記錄計算的平均揚飛日期，β為記錄的平均緯度的效應(yīng)，δ為記錄的平均經(jīng)度的效應(yīng)，Ez～N(0,σ2)為誤差。結(jié)果表明近60年來，歐洲2 000多種傳粉昆蟲的MFD提前了6 d，而FPL減少了2 d。由于傳粉者物候交錯期的縮短，導致了傳粉的功能和服務(wù)的季節(jié)性分布發(fā)生了改變。

8 小結(jié)與展望

本文回顧了昆蟲物候?qū)W研究中的熱性能曲線、生物物理模型、基于概率的模型、分布時滯模型、發(fā)育進度曲線、物候匹配模型和物候變遷模型。這些模型各有特點和適用場景。物候變遷模型研究的是物候與多種環(huán)境因子的關(guān)系，其余模型則研究的是物候與積溫的關(guān)系，或者說是發(fā)育速率與溫度的關(guān)系。但這些模型都可用于對昆蟲物候的預測。熱性能曲線經(jīng)歷了線性模型到非線性的發(fā)展過程。相比其它模型，生物物理模型是描述發(fā)育速率與溫度關(guān)系的最適模型(Wagneretal., 1984)。發(fā)育進度曲線將進入某一發(fā)育階段個體的累積成數(shù)與對相應(yīng)的累計有效溫或發(fā)育歷期(儒略日期)建立聯(lián)系，它與熱性能曲線相比，在數(shù)學表達式上無異，只是縱橫坐標所代表的變量不同，因而模型的意義不同。物候匹配模型用來描述昆蟲的物候與寄主物候的同步性。物候變遷模型通常采用線性模型將環(huán)境變量納入模型中來研究昆蟲物候?qū)Νh(huán)境因素的響應(yīng)。分布時滯模型是基于生理年齡的模型，其余的模型都是基于實際年齡的。基于生理年齡的模式要比基于實際年齡的模式更可靠一些(Rossinietal., 2020)。大多數(shù)的昆蟲物候模型是基于群組的。這類模型基于發(fā)育歷期的變異性來解釋不同群組的發(fā)育。而基于個體的模型容易引入諸如兼性滯育的發(fā)育分支途徑；一年中多重世代、單一世代和部分世代共存；以及發(fā)育響應(yīng)的變異等(參見Chuine and Régnière, 2017)。物候匹配模型和物候變遷模型并非新的數(shù)學模型，而是將現(xiàn)有模型用于相關(guān)研究，解決具體問題。種內(nèi)熱響應(yīng)參數(shù)的地理變異模式、這些模式的形成機制以及對氣候變化的響應(yīng)速度等是物候模型面臨的巨大挑戰(zhàn)，是目前的研究熱點(Chmuraetal., 2019)。在全球變化背景下人們聚焦于昆蟲在大尺度時空范圍內(nèi)的物候變遷模式(pattern of phenological shift)(Chmuraetal., 2019)。利用線性混合效應(yīng)模型，將空間和時間特征納入模型，也是目前研究的熱點之一。