4个竹种高生长模型的比较

凡美玲，方水元，冯俊娇，张春霞，梁大洪，时培建

(1.南京林业大学南方现代林业协同创新中心，竹类研究所，江苏 南京 210037；2.贵州省桐梓县林业局，贵州 桐梓 563200)

竹子属于禾本科竹亚科植物，竹林为我国南方重要的森林资源之一，其生长不同于木本植物，生长速度快、种群密度高使得竹子得以在我国南方广泛分布[1]。用合理的模型对植物生长数据进行拟合，能够模拟植物生长的动态变化，这将有益于探索和比较生物个体的生长模式[2]。上世纪以来，相继许多用于描述温度变化对外温生物发育速率的影响非线性模型被构建，许多模型被运用于拟合动植物的动态生长。澳大利亚学者Ratkowsky等[3]提出了平方根模型，此模型最初被用于描述细菌的种群增长速率，后被扩展到对节肢动物温发育速率的研究中。Shi等发现去平方根模型和除平方根的Ratkowsky模型对昆虫的温发育数据的拟合效果较好[4]，该模型是否能够用于植物生长发育有待研究。Yin等[5]首次提出beta模型用于描述温度对农作物的影响，Yin等[6]通过对beta模型的积分提出了一个新模型(BSE)，该模型能预测植物达到最大生物量的时间，为农作物的生长培育提供理论依据。Shi等[7]用BSE模型与其它4个模型对6种农作物的干重进行拟合对比，发现BSE模型拟合效果优于其它4个模型，该模型能适应于描述多种农作物的生长。Shi等用BSE模型[2]和LRF模型[8]对6种作物的干重、4个竹种高生长、2种树木原始年轮数据和16个动物的鲜重数据进行拟合比较，拟合随着时间推移各个动植物生物量(或者体型大小)指标变化的数据，发现2种增长方程与实际上动植物的观测数据吻合，且可以通过这2个模型预测动植物达到最大生物量和最高高度的时间，揭示了这2种模型功效的普适性[2]。

唐红等[9]2014年对湖南省植物园的圣音竹(Phyllostachysedulisf.tubaeformis)高生长研究发现：用Logistic方程能很好的拟合幼竹平均高生长数据，昼夜生长量可达16.48 cm。周芳纯[1]对毛竹高生长与出土时间关系研究发现：两者关系呈“S”形曲线，可用Logistic方程拟合。曾庆南等[10]对江西林业科学院竹种园花毛竹(Ph.edulis‘Tao Kiang’)的高生长节律研究表明Logistic方程同样能够较好地拟合竹高生长过程。竹类植物在幼竹高生长时期的生长速度迅速，不同竹种每天的高生长量不同，通过对4个幼竹高生长速度的比较，以期为4个竹种的栽培管理提供理论依据。以往用Logistic方程对竹子的高生长拟合较多，本研究旨在比较不同增长方程对竹子高生长数据的拟合效果的优劣，用模型拟合生物个体间不同的生长模式，并对不同模型的适用性和拟合效果进行比较，探讨适用于拟合竹类植物高生长的最优方程。

1 材料与方法

1.1 研究区概况

研究地点位于南京林业大学南京本部校园内(32°04′N，118°48′E)，海拔约25 m，属于亚热带季风气候，研究地年平均气温为(15.6±0.9)℃，年平均降雨量为(1 045±242.7)mm左右，该地季节分明，每年6-7月的梅雨季节约22 d(1986-2016年气候数据来源于中国气象数据网)。

1.2 研究材料

研究试验材料为4个竹种的竹笋-幼竹，从竹笋露出地面1～2 cm为基准，开始对竹笋高度进行跟踪调查，于2016年春季选择4个竹种(每个竹种4株)进行跟踪调查16株幼竹高生长。竹种分别为红壳竹(Ph.iridescens)、美竹(Ph.mannii)、唐竹(Sinobambusatootsik)、斑苦竹(Pleioblastusmaculatus)，发笋期均为春季。

1.3 研究方法

研究调查日期从2016年4月4日开始，直至6月25日结束调查。调查中，选用卷尺、可伸缩鱼竿测量幼竹高度(每次测定幼竹全高，即为每个幼竹从地面到幼竹的最高处)，调查初期为每天调查1次，后期直至幼竹基本完成高生长后，调查周期为2-3 d，调查时间均为中午12点，直至幼竹高度值稳定，结束调查。

研究共选择3个非线性模型对4个竹种高生长数据进行拟合，每个竹种均选择4个个体，每个模型都对16个个体的高生长数据进行拟合，3个模型的参数均为3个，比较3个模型对16个个体高生长数据的拟合优度及模型参数的线性近似表现，3个模型如下所示。

1.3.1von Bertalanffy模型 von Bertalanffy模型[11]首次由von Bertalanffy提出，2001年West等[12]提出了一个新的个体增长模型(OGM)。2014年Shi等将OGM模型用于拟合3种淡水鱼身体长度的变化，指出不同的物种应该有不同的新陈代谢率指数取值[13]。研究使用的von Bertalanffy模型，其具体形式如下：

h=H·(1-exp(-3/4·K(t-t0)))4/3

此处，h表示在时间t处生物体一维几何参数的大小(长度或者高度)。k为常数；t0亦是一个常量，表示是当竹子高度为0时的起始时间；参数H表示渐近线趋于最大值的一维几何参数(即最大长度或者最大高度)。

1.3.2BSM模型 Yin等[5]提出一个beta模型，最初是为了描述温度对玉米发育速率的影响。Yin等[6]在beta模型的基础上建立了1个3个参数的个体生长模型BSE，该模型生物量为0时对应的时间被定义为0。该模型具体形式如下：

h=copt·t·(2tmax-topt-t)/(2tmax-topt)(t/topt)topt/(tmax-topt)

此处，h表示在时间t处生物体一维几何参数的大小(长度或者高度)。模型共有3个待拟合参数，tmax表示幼竹生长到最高高度的时间，topt表示最适幼竹生长的时间(生长速率最大时)，copt表示幼竹最大生长速率的值。

1.3.3Logistic模型 Logistic模型描述动物和植物的生物量和身体长度数据运用最广泛的模型，曲线的对称性为该模型最典型的特点。该模型具体形式如下：

h=H/(1+exp(b-rt))

此处，共有3个参数，h表示随时间t变化的幼竹高度，H代表待拟合的幼竹最高高度，r代表幼竹的瞬时增长率，b代表一个常数。

1.4 参数拟合方法及模型评估的指标

研究使用R语言自编程序对幼竹高生长进行拟合，使用进化遗传算法[14]来执行对模型参数的估计，该方法是针对于解决求一个函数最大值的优化问题。研究使用决定系数(R2)和曲率评判模型的优劣，R2越大拟合优度越好[15]。在非线性回归中曲率分为均方根内禀曲率和均方根参数效果曲率(下文将省略“均方根”)，其中内禀曲率是受到模型本身制约的无法通过对模型参数的的再参数化来减小，而参数效果曲率则可以通过对模型参数的再参数化来实现变化，参数效果曲率相对于临界曲率越小，表明模型参数近似线性的表现越好，模型越合理[16-17]。

2 结果与分析

2.1 3个模型的拟合参数比较

研究共用3种模型分别对4个竹种不同个体的幼竹高生长数据进行拟合，计算结果(R2)如表1所示，从表1中的拟合结果可以看出：3个模型分别对16株幼竹个体高生长数据拟合效果都很好，R2最小也可以达到0.981 66，为von Bertalanffy模型对红壳竹4号个体的拟合结果，而拟合效果最好的出现在BSM模型对斑苦竹1号个体的拟合效果，R2最大值可达0.999 28。对拟合结果R2分析发现，Logstic模型拟合后的R2的计算值有6个数值达0.998 00，BSM模型有3个数值达0.998 00，但von Bertalanffy模型拟合后R2最大为0.996 52，没有数值达0.998 00。不同模型对不同竹种的拟合结果也不同，Logstic模型和BSM模型对斑苦竹拟合效果最佳，但von Bertalanffy模型对美竹(Ph.mannii)拟合效果较差，该模型2014年被应用于描述淡鱼体长变化，研究将其应用于植物，效果欠缺，有待改善模型结构及参数使其适应性更广泛。综上所述，Logistic模型对16株幼竹个体高生长数据拟合效果最优，BSM模型拟合效果较好，von Bertalanffy模型的拟合结果较差， Logistic模型对于诸多竹种幼竹高生长的拟合效果最好，可以被广泛运用于动植物个体增长研究。

表1 3个模型对4个竹种16株个体高生长数据集的R2拟合结果Tab.1 Coefficient of determination (R2) obtained from using the three models to fit the height growth data sets of sixteen individual shoots of four bamboo species

说明：编号Ⅰ、Ⅱ、Ⅲ、Ⅳ分别代表红壳竹、美竹、唐竹和斑苦竹，编号1、2、3、4代表每个竹种的不同个体。
The number of Ⅰ、Ⅱ、Ⅲ、Ⅳ representPh.iridescens、Ph.Mannii、Sinobambusatootsik、Pleioblastusmaculatus，and the number of 1、2、3、4 represent different individuals of each bamboo species.

2.2 3个模型参数曲率的比较

表2为3个模型对16个幼竹高生长数据拟合后的曲率对比，内禀曲率和参数效果曲率用于评估模型的非线性表现和参数结构的合理性。此处，曲率指标共有3个：内禀曲率(RMS IC)、参数效果曲率(RMS PEC)和临界曲率(CC)，当模型拟合后的内禀曲率和参数效果曲率都小于临界曲率，表明该模型近似线性的表现越强，参数结构更为合理。由表2的拟合结果可以看出：3个模型Logstic模型拟合结果中，内禀曲率和参数效果曲率均小于临界曲率，表明该模型参数的线性近似表现较好，模型参数设置合理；对于BSM模型，内禀曲率和参数效果曲率均小于临界曲率，与Logstic模型拟合结果相比，BSE模型拟合得到的内禀曲率和参数效果曲率均大于Logstic模型，表明Logstic模型的拟合结果优于BSM模型；von Bertalanffy模型拟合结果中，内禀曲率均小于临界曲率，但有4个参数效果曲率大于临界曲率，参数效果曲率均>BSM模型>Logstic模型。综上所述，Logistic模型对16株幼竹个体高生长数据拟合优度最高，参数设置更为合理，其次为BSM模型，von Bertalanffy模型的拟合优度最低[17]。

2.3 3个模型的拟合曲线比较

图1为BSM模型、von Bertalanffy模型和Logistic模型3种不同模型对斑苦竹中1号幼竹高生长数据的拟合结果(3个模型对该个体拟合效果最好)。研究中的斑苦竹种植于南京林业大学竹种园内，1号斑苦竹于2016年4月29日破土而出，直至2016年6月25日完成其高生长，最高达12 m。开始生长速度较慢，最快速度出现于6月6日，此时增长速率达到72 cm·d-1，随后生长速度减慢，其生长曲线呈典型的“慢-快-慢”趋势。

表2 3个模型对16个个体高生长数据集拟合后得到的曲率比较Tab.2 Curvatures obtained from using three models to fit the height growth data sets of sixteen bamboo shoots

研究共选择了3个非线性模型对同一竹种不同个体高生长数据进行拟合，3个模型都设置了3个待拟合参数(BSM模型中tmin设定为0)。von Bertalanffy模型是一个通用型模型，可以用于诸多生物的生长，研究该模型新陈代谢率指数选用的为3/4，但von Bertalanffy模型对幼竹高度的拟合曲线与实际数据有很大出入(图1c和图1d)，该拟合曲线在前10 d左右零生长，随后生长速度骤增，但在曲线末端未见曲线平缓，有继续增加的趋势，这与实际观测数据不符，幼竹在完成高生长后，高度会趋于一个稳定值，其下方对幼竹生长速率的拟合曲线同样可以看出，该幼竹前期生长速率为0，随后在10 d左右达到最大生长速率，达到最大生长速率共用24 d，其后幼竹继续生长，生长速率线始终没有趋于0，且整条生长速率曲线不对称，拟合后决定系数计算结果为0.988 45，该模型对其幼竹高生长拟合效果差。BSM模型最大的特点是能够包容对称和不对称曲线，能准确估计植物到达最大生物量的时间。图1a和图1b展现了BSM模型对斑苦竹中1号幼竹高生长数据及生长速率的拟合结果，该拟合结果表现为前期生长缓慢，随后进入高生长盛期，曲线最后趋于平缓，幼竹高度达到最大高度后稍有下降(幼竹开始伸枝展叶)，最后稳定于一个值，拟合曲线与幼竹实际高度重合率高。从该模型对其生长速率拟合曲线可以发现：高生长阶段达到最大生长速率的时间为距离出笋35 d，52 d后结束高生长，该生长速率曲线左侧宽于右侧表明，竹子前期到达最大生长速率的时间略多于后期生长的时间，拟合后决定系数计算结果可达0.999 28，该模型对其拟合效果较优。Logistic模型被广泛应用于描述动植物的生物量，从图1e和图1f可以看出，幼竹高生长呈“慢-快-慢”趋势，拟合曲线与幼竹高生长数据高度重合，与BSM模型不同的是，该模型对生长速率的拟合曲线更平缓，且曲线两侧对称，幼竹最大生长速率时间为距出笋33 d，与实际数据相符，拟合后决定系数计算结果能够达到0.999 01，综上所述，Logistic模型拟合结果优于BSM模型优于von Bertalanffy模型，在本研究中，3个模型对不同幼竹个体高生长数据的拟合结果表明，Logistic模型为最优模型。

图1 用3种模型对斑苦竹(Pleioblastus maculatus)1号个体高生长数据的拟合结果Fig.1 Fitted results of the height growth for No.1 bamboo of Pleioblastus maculatus by using three models说明：a、b为BSM模型；c、d为von Bertalanffy模型；e、f为Logistic模型。a、b are the BSM model; c、d are the von Bertalanffy model; and e、f are Logistic model

3 讨论

竹子每年的发笋时间主要受当年的温度和降雨量等因素影响，按发笋期发笋数量可分为初期、中期和末期，末期笋常因营养不足或虫害退笋[18-20]。实验选择是4个竹种发笋期差异主要因为竹种差异造成，选择发笋前期到中期的笋为研究对象，调查发现红壳竹与美竹发笋日期与高生长期相近，于4月初开始发笋，6月伊始结束高生长，唐竹和斑苦竹发笋期略晚，分别为4月中期和末期，高生长期约为40-55 d。

研究使用的von Bertalanffy模型，参数=3/4，但总体拟合效果差，同时也尝试参数=2/3，拟合优度低于前者，2014年Shi等[2]用该模型和OGM模型同时对3种淡鱼体长数据进行拟合，他们发现当OGM模型中的参数α越接近1，模型的拟合效果越好，与本研究的结果一致。本研究中，使用的BSM模型中限制δ=1、tmin为0，但并不会影响该模型描述数据的有效性，结果表明BSM模型对4种幼竹高生长数据拟合效果较好，对幼竹生长到最大生物量的时间估计与实际相符，Shi等[2]用2种生长模型分别对动植物生长进行拟合，作者指出：如果可以预先确定动植物增长的起始时间，那对于BSM模型的参数δ可以不限制为1，可以放宽δ作为一个参数，这样将帮助对比动植物种间增长模式的差异。对于常用的生长模型Logistic模型，拟合结果同前人研究一致，对本研究中的4种竹种个体高生长数据拟合效果好，Shi等13用4种模型同时对相同数据拟合比较发现，Logistic模型对于描述动植物生物量和身体长度都是通用的，且拟合优度都优于其它模型，是最适合描述动植物的生长模型。Richards[21]提出一个广义的Logstic模型，即Richards模型，该模型更加灵活，不再限制生长速率曲线的对称性，且该模型对于在其最小值到最大值之间的任何值都能有其最大斜率，但该模型也存在一些争议，Ratkowsky指出该模型对系数估计值不准确[21]，Birch[23]也指出该模型不能反映一些蔬菜和牧草的生长，如黑麦草(LoliumperenneL.)在低密度时就能达到其最大净产量，与实际不符。研究表明：许多非线性模型不仅适用于温度对动物生长速度和内禀增长率的影响，对于植物的生长也是适用的。

最初一些非线性模型被提出用于描述温度对昆虫发育速率的影响，如表现型模型[24]、Taylor模型[25]和Ratkowsky模型[4]等等，这些模型可以帮助预测一些重要参数，对农业生产及科学研究有着重要意义，但因为模型结构设置及参数选择不同，对不同数据的拟合优度不同，选择一个通用性的最优模型具有重要意义。评价模型拟合优度的好坏可以从几个方面评估，模型结构的复杂程度、模型适用性以及模型预测曲线与实际是否相符等[2]。随后，王如松等将数学模型用于描述温度变化对昆虫内禀增长率的影响，拟合结果显示相关系数高，拟合结果较优，指出可以利用数学模型研究一些与温度相关的现象[26]。2016年Shi等提出BSM的通用表达式(即生长发育初始时间不为0时的BSM模型)能够拟合植物的动态生长，并能帮助预测植物生长到最大生长量的时间，这为森林资源的合理利用提供了参考价值[7]。Shi等[13]将1种个体增长模型运用于作物的生长，结果表明Logistic模型能够更好地拟合动植物生长的数据，这与本研究结果一致。本研究中，BSM模型与Logistic模型对幼竹生长到最大生物量时间的预测相近，与实际数据吻合，而von Bertalanffy模型预测幼竹最大生物量的时间提前10 d，正如图1d所展示的生长速率曲线左偏。同样，Logistic模型拟合的生长曲线同前人研究一样，曲线两侧关于对称轴对称，对数据缺乏灵活性[27]，但总体Logistic模型适合描述植物生长过程中生物量的变化趋势。

研究使用3个非线性模型对不同竹种幼竹高生长数据进行拟合，将以往描述温度对昆虫发育速率影响的模型用于拟合植物的生长，将自变量温度调换为时间，结果表明，生长模型也能很好地描述植物的生长[2]。调查发现，不同竹种前期与盛期笋生长为成竹的概率较大，末期笋易发生退笋现象，所以在选择留笋养竹时，应选择初期和盛期萌发的笋为宜。利用3个模型对不同幼竹高生长数据进行拟合，4个竹种高生长均符合“S”型生长，符合“慢-快-慢”是生长规律，这为4个竹种的抚育管理提供了重要的理论依据。本研究发现，Logistic模型对动植物生长的拟合效果最优，BSM模型与von Bertalanffy模型有待改进，这些生长模型的适用性有待进一步探索，生物有机体生长模式的比较也将成为一个研究热点。