昆虫物候模型研究进展

景天忠，刘丽萍，谢雨龙

(东北林业大学林学院，哈尔滨 150040)

物候模型(phenology model)预测昆虫发育事件的时间。物候模型已成为种群动态、物种分布和进化动态的关键成分(Chuine and Régnière, 2017)。在害虫管理生产实践中，需要预测出害虫某一发育事件出现的时间，以便确定害虫防治的最佳时间，称为发生期预测，其实就是对昆虫物候学的预测。美国国家物候学网络中心(U.S. National Phenology Network，USA-NPN)在2018年开发了一套可用于科研和生产管理的物候预测图(Pheno Foescast map)(Crimminsetal., 2020)，包括实时图和在与管理相关的时间和空间分辨率上短期预测害虫的活动。Nietschkeetal.(2007)建立了一个昆虫发育数据库(Insect Development Database, IDD)，包含405种昆虫的基础发育温度阈值和310种昆虫产卵和至成虫羽化的累计日度。这些昆虫分属11个目和89个科。加州大学州级IPM项目(http://www.ipm.ucdavis.edu/MODELS)提供了许多昆虫的物候信息。美国还建立了一个在线的物候和日度模型，用于农业和害虫管理决策(http://uspest.org/cgi-bin/ddmodel.us)。

对变温动物而言，完成一个特定的发育阶段需要累积一个种特异的温度常数(thermal constant)，这就是最早由Réaumur在1735年提出的有效积温法则(参见Ikemoto and Takai, 2000; Damos and Savopoulou-Soultani, 2012)，可用公式表示：

K=D(T-T0)

(1)

其中K为温度常数，即有效积温(effective cumulative temperature)；T为环境温度；T0为发育起点温度；D为发育历期(development time)，即昆虫完成一定发育阶段所经历的时间，常用τ(T)来表示。而发育历期的倒数，便是发育速率(development rate)，常用r(T)来表示。要对昆虫在变温条件下的发育进行建模，可采用积分的方法，即

(2)

其中，a为某一阶段的生理年龄(physiological age)，表示完成某一阶段发育的成数，取值在0～1。0表示发育的开始，1表示发育完成。A为由温度响应函数r(T,A)的参数值组成的向量。昆虫完成发育所需的热量也称为生理学历期(physiological time)。

有效积温法则是对昆虫物候研究的基础，因而昆虫物候模型总是围绕温度、发育速率或发育历期这些变量。在物候模型的发展过程中出现了热性能曲线、生物物理模型、基于概率的模型、分布时滞模型、发育进度曲线、物候匹配和物候变迁等模型。除热性能曲线和发育进度曲线外，其余模型在中文文献中鲜有提及。因而很有必要对这些模型做一较为系统的介绍，以备借鉴。

1 热性能曲线

对式(1)两边取倒数，得

(3)

称为线性日度模型(linear degree-day model)或X截距模型(x-intercept method)，是一个标准的一次函数模型(y=1/Donx=T)，可用最小二乘法求出参数。1/D为发育速率。

有学者指出，上述公式存在三个不足：一是利用上述公式很难得到最适温度范围；二是数据点的权重在直线的上部和下部之间不成比例，且沿X轴并不同质。温度越低，担负的权重越不相称(越少)；三是对大多数数据，变量X的误差被忽略了(Ikemoto and Takai, 2000)。因而提出一了个新的线性模型：

(DT)=K+T0D

(4)

很明显，利用有效积温法则，可预测害虫和天敌的发生期，预测害虫和天敌的地理分布及在某地一年发生的代数。但要注意，昆虫生活场所的气温和百叶箱所代表的大气气温有差别。室内所计算的发育起点和有效积温常常是在恒温条件下所得的数据，与自然变温状态下不同，变温下昆虫的发育常较在恒温下的为快(Colinetetal., 2015)。积温法则无法应用于某些有滞育阶段的昆虫世代数的计算。昆虫的发育不仅仅受温度一个因素的影响，还受食物、湿度等其他因素的影响。还有，昆虫的发育速率与温度之间通常是非线性的关系(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和指数方式改进的Gaussian等，它们通常以温度为横轴、发育速率为纵轴，称为热性能曲线(thermal performance curves/functions/equations/models)(Shi and Ge, 2010; Rebaudo and Rabhi, 2018)。R语言的devRate包可拟合这些曲线。Geng and Jung(2018)建立了金纹细蛾Phyllonorycterringoniella越冬蛹发育速率的Logan-6模型。这些模型还用于研究秘夜蛾Mythimnasequax的发育速率(Jacquesetal., 2019)、预测棉铃虫Helicoverpaarmigera在西红柿Lycopersiconesculentum上的物候(Dalal and Arora, 2019)，以及螽斯Ruspoliadifferens当前和未来的分布(Leonardetal., 2021)。Shi and Ge(2010)利用小菜蛾Plutellaxylostella和B型烟粉虱Bemisiatabaci为研究对象，以AIC(Akaike information criterion)、BIC(Bayesian information criterion)等方法比较了12种非线性模型，结果表明Performance，Brière1和Brière2三种模型最适合应用于温度依赖的发育速率。Régnièreetal.(2012)开发了一个经验模型，可对近上、下临界温度的热响应更好的描述。Anarcaetal.(2018)将寄主质量包含在模型中来提升模型。

2 生物物理模型

有人根据范霍夫定律(Van’t Hoff’s law)提出了生物物理模型-Sharpe-Schoolfield模型(参见Damos and Savopoulou-Soultani, 2012)。生物物理模型将变温动物的温度依赖的发育看作是一个宏观的酶促反应。基于Arrhenius和Eyring关于化学反应速度的理论，该模型最初由Johnson & Lewin、Hultin、Sharpe和DeMichele等人建立，后由Schoolfield、Sharpe和Magnuson改进(Schoolfieldetal., 1981)。温度为T(绝对温度，K)时，平均发育速率r(T)(time-1)为：

(5)

Régnièreetal.(2012)开发了一个经验函数，可更好的描述临界温度的热响应：

(6)

其中Tb≤T≤Tm，6个参数与式(5)类似。式(6)与(5)在大部分温度范围内重合，但在临界温度周围有很大的不同。

为处理卵在高温下由于高死亡率而无法计算ΔHH的问题，Salisetal.(2016)采用缩减版的Sharpe-Schoolfield模型：

(7)

并开发了一个发育阶段依赖的模型(developmental stage-dependent model, DSD)：

(8)

其中，ST为用于校正温度对发育速率影响的缩放因子；Dth为发育速率不受发育阶段影响时的发育阈值；D为发育阶段，即从产卵(D=0)至卵孵化(D=1)的发育进程。

Crespo-Pérezetal.(2015)利用式(5)研究3种马铃薯块茎蛾Phthorimaeaoperculella的日发育速率，Jacquesetal.(2019)研究了Spodopteralitura在印度落花生上的定殖和存活。生物物理模型还用于研究个体大小与温度的关系(de Jong, 2010)。

3 基于概率的模型

将昆虫饲养在完全一样的条件下，不同个体的发育历期不尽相同。这种变异性可能来自于一个种群不同个体间的代谢和遗传方面的差异。诸多研究表明，这种变异性是不对称分布的，即右偏。根据式(5)，影响发育速率的参数中，酶是唯一的受遗传变异影响的生物学参数，且研究表明生物的一些遗传特征如头壳宽度、长度、大小、智力等都是对称的分布。因此，Sharpeetal.(1977)认为限速酶的浓度也应该是对称分布。在这些假设下，式(5)可改写为

R(T)=r(T)ε

(9)

其中ε为服从正态分布的随机变量。

这样就可利用正态分布的密度函数或者二次项分布密度函数来描述发育速率的变异(横轴为发育速率，纵轴为概率密度)，即

(10)

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

(11)

发育历期的密度函数(横轴为发育历期，纵轴为概率密度)为

(12)

(13)

这类模型称为基于群组(cohort-based)的模型，或者称为发育历期的分布模型。群组为一个特定生命阶段(life stage)的所有个体，这个阶段中的个体具有几乎一样的时间年龄(chronological age)。

在发育速率变异是正态分布的情况下，可以把昆虫种群中16%、50%和84%左右的个体进入某一发育事件的日期作为生产中划分害虫发生始盛期、高蜂期和盛末期的数量标准，对应于正态分布-σ、μ和+σ分位数。

设τi(T)为单个昆虫在恒温T下完成其第i个生命阶段所需的时间，T(t)为t时刻的温度、a(t)为在t时刻一头昆虫完成某一生命阶段的成数，a的取值范围在0到1之间，可以看作是昆虫的生理年龄，0表示起始发育，1表示发育完成。那么在温度T下，昆虫完成某一生命阶段Δa部分所需的天数为Δt，则Δt=τi(T) Δa。那么，一头昆虫从t0开始的阶段i的发育可以用下式来描述：

(14)

(15)

一头昆虫，从时刻t开始，完成阶段i的发育所需的时间为Γi(1,t)，且对所有的时刻t，Γi(0,t)=t。当Γi已知时，这头昆虫的物候学就可以确定了(Yurk and Powell, 2010)。这类模型称为基于个体的模型(individual-based models)。

假设发育历期不变，上面这个个体物候学模型可以扩展成种群模型。令p(a,t)为t时刻具有年龄a的个体的密度函数，则由完全相同个体(不存在变异)组成的种群的发育可用一个平流方程(advection equation)来描述：

(16)

在上述模型中引入一个能够缩放基础发育历期曲线的表型参数α，就可解释发育历期中的持续变异。这个表型在种群内具有相同的分布。在这种情况下，式(16)变成为一个α依赖的平流方程

(17)

如果种群的每个个体都具有相同的表型α=1，则(17)式变为式(16)。

在式(17)中加入一个扩散项(diffusion term)来将由随机效应引起的发育变异加入到表型模型中，称为Fokker-Planck发育方程(Fokker-Planck development equation)：

(18)

如果种群的每个具有表型α的个体具有相同的发育历期(即变异v=0)，则式(18)变为式(17)。

在恒温条件下，发育历期可由下式来预测：

(19)

即发育历期服从均值为μτi(T)、方差为v+σ2τi(T)2的正态分布。

在变温条件下，假设(1)描述不同恒温条件下发育的PDF积分函数可度量不同温度下的发育；(2)发育速率的均值和方差均是呈线性比例，则对每一个温度T，f(r)的均值为r0，标准差为σ=cr0，应用Sharpe模型，处于生命阶段j的个体在t时刻的历期(出现时间emergence times)分布为(Gilbertetal., 2004)：

(20)

其中，可用σ对r0的线性回归来估计c。

通常对发育的观测是以dt日为间隔的。如果定义个体与真平均值之间的差异δij服从对数正态分布，理论响应值与处理组均值之间的拟合劣度υj服从正态分布，那么在观测区间[t-dt,t]上，个体i在处理j下以恒温Tj完成某一发育阶段的概率pij为

(21)

当存在删失数据时，无法估计一个处理的平均发育历期。假定对该处理υ=1，一个个体完成其大于tij的发育的似然性为1-F(εij,A) ，则无论是否存在删失数据，可用下式来表示似然性

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

(22)

其中，d=0或1，0表示无删失数据，1表示有删失数据。然后可用最大似然法求出各参数。

当估计昆虫发育的上下临界温度时，昆虫易在近临界温度处死亡。Régnièreetal.(2012)提出了利用温度转换(temperature transfer)来解决此问题。先将昆虫暴露于近临界温度(T1)下一段固定时间(t1)，既要短到避免过多死亡，又要长到保证明显发育。然后将昆虫转移到另一个适合其发育的温度(T2)下t2,i时间，使其完成发育。在这种情况下，当式(23)成立时，一个个体完成发育。

(23)

(24)

Gilbertetal.(2004)在McKendrick-von Foerster模型基础上，开发出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))为发育速率，v(T(t))为发育变异。上式的解为

(27)

假定发育过程中没有生产，且r服从正态分布，个体完成了一个发育阶段(a=1)，则式(27)变为：

(28)

其中，r0为平均发育速率。发育的概率即为发育速率r，可由平均发育速率r0来近似表示。

在变温的环境条件下，处于发育阶段j的个体在t时刻的种群分布函数为

(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)。设x(t)为一时间函数，则x(t-r)中的时滞r称为离散时滞(discrete delay)，如Leslie矩阵就属于离散模型，而形如

(30)

的式子称为分布时滞，因为它反映了时滞x(t-z)的加权平均。上式实际上是卷积的表达式。

昆虫的发育通常是不整齐的，例如在卵孵化为幼虫阶段，每天都有新幼虫加入到幼虫的行列中，这样种群就具有了一个年龄分布，可以用一个代表其居留时间(residence time，dwell time)或存活时间(survival time)的随机变量来描述。可用Erlang家族函数、正态分布或二次项分布概率密度函数来描述居留时间的分布，但Erlang分布最常用，其PDF为：

(31)

其中，k称为形状(shape)参数，r为速率(rate)参数，1/r称为尺度(scale)参数。其均值μ=k/r，方差σ2=k/r2。它是gamma分布的一个特例，当k=1是时便简化为指数分布。R的VGAM包括提供了拟合Erlang分布的函数。

将一个事件X或一种生物的生命历程或生命阶段划分成k个离散的阶段，或称为闷子车、车厢(compartment)或相位(phase)，第k个闷子车完成所需的时间为k个独立同分布变量的和。用Xi来表示任一闷子车，则X1,…,Xk组成X。当Xk完成了，X就完成了。假设闷子车Xi在t时刻有xi(t)个个体，每时刻有一部分个体(ri(t)，成数)进入下一个闷子车。这里的ri称为流速(flow rate)或转移速率(transition rate)，也即所谓的分布时滞。假定没有死亡，也没有迁入迁出，所有个体都完成最后发育，则可用下面的微分方程来描述种群的变化

(32)

设a为一连续变量，代表年龄，归一化后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展开，并假定整个生命阶段的发育是一个常数(ri-1=ri=ri+1=r)，则式(32)变为

(33)

其中，rΔa代表了一个个体转移到下一个闷子车的概率。解此方程，可得

(34)

在恒温下，分布时滞模型实质上von Foerster方程特征线法解的数值近似。一个重要区别是在分布时滞模型中每一个发育阶段的相位数(number of phases)必须通过数据来估计(Gilbertetal., 2004)。分布时滞模型在变温条件下很难与恒温条件下的参数连接起来，因为为每一个温度和每一个生命阶段来指定闷子车数实际上不可能的。分布时滞模型将发育过程比做一群个体通过一系列闷子车，其输出取决于闷子车之间的流速(rΔa)和闷子车数(k=1/Δa)。

Jianetal.(2007)开发一种用于时间可变分布时滞模型的算法，研究了不同温度下锈赤扁谷盗Cryptolestesferrugineus老龄化的速率。

Scranton and Amarasekare(2017)使用时滞微分方程(delay differential equations, DDEs)预测了气候变化背景下昆虫物候的变化。该模型将昆虫的生命周期分为幼年期(J)和成年期(A)两个阶段：

(35)

其中变量J和A分别代表幼年期和成年期的密度，函数B(T(t),A(t))和DX(T(t),X(t))描述的是密度和温度对平均出生率和死亡率的联合效应；RA(t)是温度依赖的补充率，即在t时刻的成熟个体，代表在τ(t)时间单位前出生，历经发育阶段τ(t)后存活为SJ(t)的个体。这里τ(t)就是发育迟滞；温度依赖的发育速率MJ(T(t))建立初始值为1/MJ(T(0))的发育迟滞τ(t)。结果表明，多度的峰值在环境温度变暖的年分提前，但在季节性振荡振幅增加的年分延迟。

5 发育进度曲线

在发生期预测工作中，还可以对进入某一发育阶段个体的累积成数(cummulative proportion)对相应的累计有效温或发育历期(儒略日期)进行拟合。如Geng and Jung (2018)建立了金纹细蛾Phyllonorycterringoniella成虫累计羽化速率的Weibull模型、Damos and Savopoulou-Soultani(2010)采用Boltzman和Logistic回归两种模型研究了桃园中主要的鳞翅目有害生物复合体的发育进度。使用4参数的Weibull模型，Milosavljevietal.(2017)拟合了两种金针虫Limoniuscalifornicus和L.infuscatus在土壤中的累计捕获率对有效积温的曲线，D’Auriaetal.(2016)拟合了3种马铃薯害虫的累计羽化率对有效积温的曲线。以上这些例子都以累计有效积温为自变量、以累计发育进度(成数)为因变量进行曲线拟合。也有以儒略日期为自变量进行曲线拟合的。但由于年际间气温存在变异，使用累计有效积温的模型比使用儒略日期的模型要稳健得多。相比之下，使用儒略日期的模型的好处是数据容易获取，在生产实践中更容易应用。

6 物候匹配模型

物候匹配模型用来描述昆虫的物候与寄主物候的同步性。这类模型在授粉昆虫中研究较多。

Hindleetal.(2015)采用一个统计模型来表示一种蝴蝶Melanargiagalathea在飞行期的某日t(儒略日)能见到的个体期望数量μ(t)：

(36)

(37)

其中，φ为方差参数，a=μ/φ，b=1/φ。假定φ独立于年和地点。

该蝶的主要蜜源植物Centaureascabiosa的开花时间和花期在地点间和年间是不同的，可用下面的模型来模拟在儒略日x时的期望累计开花率：

(38)

(39)

其中，a=p/φ，b=(1-p)/φ。φ定义了植株间的变异程度。由于它在地点间、年间的微栖境之间是相似的，令其在所有模型中都以常量存在。以最简约模型(AIC最小)来计算开花期，并将5%～95%的花开日子定义为花期。

最后，作者采用线性回归来判定蝴蝶的飞行期以及植物的花期是否随着年显著变化。结果表明，蜜源植物的开花时间变晚了，但蝴蝶的飞行期并无明显的变化趋势。

7 物候变迁模型

物候变迁模型通常采用线性模型或加性模型将环境变量纳入模型中来研究昆虫物候对气候变化的响应。这类模型非常灵活，可依研究目的采用广义线性模型、广义加性模型(generalized additive models, GAM)、广义线性混合效应模型(generalized linear mixed models, GLMM)等。下面是近年来的几个例子。为研究环境变量对昆虫物候变化的驱动力，Searleetal.(2013)采取了一个两相建模途径。首先用GAM模型来拟合每个诱捕点每周捕获的未产仔瘿蚊DiarthronomyiachrysanthemiAhlberg的三周移动平均数。假定原始计数数据服从泊松分布，并采用log链接函数。使用三周移动平均数来减少每周诱捕过程中由于气象条件的变异引起的多度随时间的错误变化。这些模型用来确定两种瘿蚊每年每地点上不同世代高峰的多度峰值时刻。然后确定以下物候性状：第1代或第2代的峰值时刻、与第1代或第2代联系的多度、第1代与第2代之间的多度差别(取对数后相除)。接下来第二相分析使用GLMM模型将这些物候性状与环境因子联系起来。

Hollowayetal.(2018)使用熵来筛选环境变量，然后使用C5.0决策树算法来预测气象因子对蚜虫扬飞驱动力。响应变量为二元的观测值：首飞(first flight)和未见扬飞(no flight)。由于是连续的监测，将首飞前的7～105 d视为未见扬飞。环境变量有日气温、日气压、日北大西洋涛动、累积日度等。结果表明，与一般加性模型相比，决策树模型对首飞预测的准确率提高了20%。通过熵筛选的变量与专家推荐的变量相比，预测的准确率也提高了3%～5%。

Stemkovskietal.(2020)使用GAM来预测蜜蜂的物候相(觅食起始、觅食高蜂和觅食终止)。对每一个物种/样地/年度组合，以积日(day-of-year)为解释变量，具有高斯分布族的多度为响应变量拟合GAM模型。使用三次样条进行平滑，使用交叉验证来避免过拟合。建立如下混合效应加性模型来研究物候的驱动力：

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

其中，DOYphase为每个物候相的积日估计值，θclim为气候变量(雪融日期、夏季气温和夏季降水)，θtopo为拓扑变量(海拔和太阳入射角)，θtrait为物种性状(体重、蜂巢位置、越冬虫态)，esite为样地，esp为蜜蜂种类。θ项为固定效应，e项为随机效应。

使用蜜蜂物候变化的天数除以雪融日期变化的天数得到的商来表示蜜蜂对雪融的响应情况，正值表示提前、负值表示延迟。结果表明，蜜蜂的觅食起始物候相对气候因素非常敏感，随着雪融日期的提前而提前。

Gutiérrez and Wilson(2021)使用具有高斯误差结构的广义最小二乘模型来为每种昆虫的平均扬飞日期和首飞日期对气候因子建模型。结果表明，对所有的蝴蝶种类，成虫羽化前几个月的气温与物候的联系最强。所有的种类在较暖的年份都扬飞得较早，在飞行季节较早飞行的种类表现出对年度气温变异的敏感性最强。大多数种类的物候对气温的敏感性高于对气温空间变异的敏感性。

Duchenneetal.(2020)建立如下模型来检测多模态(multimodality)：

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

(40)

其中，Yik是第i年第k个观测的儒略日，μ为总平均数(截距)，ρ和τ分别表示纬度和经度效应，θ为海拔效应，φi为年随机效应(年度)，Eik为误差项。该模型的残差就代表了去掉空间和海拔变异的采集日期。

为检测这些残差分布的多模态，可对分布进行平滑。起初从494种昆虫中发现了些模式，然后对每一种昆虫通过科学文献或灰色文献(grey literature)来验证是否存在多模态的物候，结果表明288种昆虫存在多模态的物候。对这些种类进行下一步混合高斯模型(Gaussian mixture-models)聚类分析，给每一条记录分配一种物候模式。然后再一次运行与式(40)相似的线性混合效应模型：

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

(41)

其中βj称为物候模式效应。最后从2 473种昆虫发现了2 473个单模态(unimodal)的物候模式。接下来估计平均扬飞日期(mean flight date, MFD)和扬飞期(flight period length, FPL)的变化。对每个物候模式，使用下面的模型来预测采集日期：

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

(42)

其中Yk为观测点k的儒略日，π为时间效应，其余同前。最后使用系统发育广义最小二乘模型(phylogenetic generalized least squares model, PGLS)来检验季节性早熟以及物种的空间分布是否与物候的改变相联系：

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

(43)

其中，PSz为z种昆虫在物候上的变化(即MFD或FPL的变化)，μ为总平均数(截距)，α为由近期记录计算的平均扬飞日期，β为记录的平均纬度的效应，δ为记录的平均经度的效应，Ez～N(0,σ2)为误差。结果表明近60年来，欧洲2 000多种传粉昆虫的MFD提前了6 d，而FPL减少了2 d。由于传粉者物候交错期的缩短，导致了传粉的功能和服务的季节性分布发生了改变。

8 小结与展望

本文回顾了昆虫物候学研究中的热性能曲线、生物物理模型、基于概率的模型、分布时滞模型、发育进度曲线、物候匹配模型和物候变迁模型。这些模型各有特点和适用场景。物候变迁模型研究的是物候与多种环境因子的关系，其余模型则研究的是物候与积温的关系，或者说是发育速率与温度的关系。但这些模型都可用于对昆虫物候的预测。热性能曲线经历了线性模型到非线性的发展过程。相比其它模型，生物物理模型是描述发育速率与温度关系的最适模型(Wagneretal., 1984)。发育进度曲线将进入某一发育阶段个体的累积成数与对相应的累计有效温或发育历期(儒略日期)建立联系，它与热性能曲线相比，在数学表达式上无异，只是纵横坐标所代表的变量不同，因而模型的意义不同。物候匹配模型用来描述昆虫的物候与寄主物候的同步性。物候变迁模型通常采用线性模型将环境变量纳入模型中来研究昆虫物候对环境因素的响应。分布时滞模型是基于生理年龄的模型，其余的模型都是基于实际年龄的。基于生理年龄的模式要比基于实际年龄的模式更可靠一些(Rossinietal., 2020)。大多数的昆虫物候模型是基于群组的。这类模型基于发育历期的变异性来解释不同群组的发育。而基于个体的模型容易引入诸如兼性滞育的发育分支途径；一年中多重世代、单一世代和部分世代共存；以及发育响应的变异等(参见Chuine and Régnière, 2017)。物候匹配模型和物候变迁模型并非新的数学模型，而是将现有模型用于相关研究，解决具体问题。种内热响应参数的地理变异模式、这些模式的形成机制以及对气候变化的响应速度等是物候模型面临的巨大挑战，是目前的研究热点(Chmuraetal., 2019)。在全球变化背景下人们聚焦于昆虫在大尺度时空范围内的物候变迁模式(pattern of phenological shift)(Chmuraetal., 2019)。利用线性混合效应模型，将空间和时间特征纳入模型，也是目前研究的热点之一。