-
开展森林精准经营[1],是未来林业的发展方向,这需区分考虑林分内每株树的生长特征及受不同因素的影响,其中单木树高是最重要的林木特征之一。华北落叶松是我国北方暖温带湿润半湿润气候区山地的主要造林树种,具有速生及适应性强等优点[2]。华北落叶松单木树高生长过程同时受林分结构因子(如林龄、密度、郁闭度等)[3]及立地因子(气候、地形、土壤等)[4]的影响。对华北落叶松林的前期研究表明,林龄直接影响单木树高生长和其他林分结构特征变化[5];林分密度和林分郁闭度等林分结构可影响树木对环境资源(养分、水分、光照等)的占有和利用,从而影响单木树高的生长、竞争、死亡和产量[6-7]。影响林木树高生长的主要立地因子是海拔、坡向、坡度和土层厚度等[8],它们主要通过影响光照、温度、养分和水分而影响树高生长过程[9]。单木模型以树木个体自身生长特征为基础,从林木生长的竞争机制出发,模拟单木生长过程[10]。在以往研究中,采用解析木数据建立了兴安落叶松(Larix gmelinii (Rupr.)Kuzen)的单木树高模型,但并未考虑林分特征因子以及立地因子对单木树高生长的影响[11];为反映多个因素对单木树高生长的综合影响,可用外包线法首先确定生长指标对单因子变化的响应函数,然后建立多因子耦合模型,如杨文娟[12]在祁连山研究青海云杉(Picea crassifolia Kom.)林分生长响应林龄、密度、海拔和坡向的模型,田奥建立了六盘山华北落叶松林分树高与林龄、密度和海拔的耦合模型[3]。在六盘山区华北落叶松林以往的研究中都只建立林龄-树高模型[13],但还没进行过单木树高生长的多因素影响及耦合模型研究。因此,本研究在六盘山区选择华北落叶松人工林典型样地,开展每木调查及解析木分析,定量研究单木树高受立地条件和林龄及其他林分结构的影响,并建立能反映多因子影响的树高生长耦合模型,以便为选择适宜造林立地、预测单木树高生长和林分结构变化、开展森林精准经营等提供理论和技术支撑。
HTML
-
本研究在宁夏六盘山南段东坡的香水河小流域(106°12′10.6″~106°16′30.5″ E,35°27′22.5″~35°33′29.7″ N)进行,海拔变化在2070~2931 m,属温带半湿润气候,年均气温3.7 ℃,年均降水量 671 mm[14]。流域内土壤以灰褐土为主。小流域森林覆盖率高达82.91%;以华山松(Pinus armandii Franch.)、白桦(Betula platyphylla Suk.)等天然次生林为主,占小流域面积的58.51%;人工林以华北落叶松纯林为主,占小流域面积的24%[15];灌丛面积占12.01%,主要有西北栒子(Cotoneaster zabelii Schneid.)、沙棘(Hippophae rhamnoides Linn.)等;草地面积占16.25%,主要有狼针茅(Stipa baicalensis Roshev.)、早熟禾(Poa annua L.)等;草甸面积占4.64%,主要有苔草(Carex tristachya Spp.)、蕨(Pteridium aguilinum L.)等[15]。
-
经全面踏査后,沿海拔梯度(2 000~2 200、2 200~2 400、2 400~2 600、2 600~2 800和2 800~3 000 m)选择了23块面积20 m × 20 m的华北落叶松纯林样地(表1),记录样地经纬度及海拔、坡度、坡向等立地条件。调查中,将正北方向记为0度,顺时针偏离正北180度以内的坡向为正,逆时针偏离正北180度以内的坡向为负。
样地编号
Sample plot
number林分密度
Stand density/
(tree·hm−2)解析木
Analytic
trees林龄
Age/a海拔
Elevation/m坡度
Slope/°坡向
Slope
aspect/°郁闭度
Canopy
density平均胸径
Mean
DBH/cm平均树高
Mean tree
height/m枝下高
Branch
height/m平均冠幅
Mean canopy
width/m1 1 850 5 14 2 033 10 15 0.72 10.30 8.50 0.85 1.48 2 1 475 5 17 2 042 15 68 0.70 11.13 9.70 1.51 1.68 3 2 000 5 16 2 086 20 115 0.70 10.09 8.73 1.34 1.51 4 1 150 5 36 2 333 18 66 0.67 18.60 19.07 6.38 1.67 5 775 5 33 2 375 20 121 0.70 20.91 20.65 6.64 1.77 6 1 100 5 34 2 285 35 −154 0.72 17.47 17.21 6.42 1.69 7 775 6 32 2 355 19 −100 0.69 20.06 19.04 6.11 1.77 8 1 275 5 32 2 560 33 156 0.71 15.21 12.98 5.12 1.59 9 800 5 32 2 624 16 −104 0.69 18.77 13.98 5.07 1.55 10 1 050 5 34 2 570 25 −122 0.72 18.03 15.87 5.66 2.39 11 1 000 5 32 2 626 23 −160 0.62 19.12 16.11 5.48 1.77 12 1 300 5 33 2 631 20 86 0.67 16.75 15.87 5.27 1.31 13 1 425 5 22 2 693 22 −64 0.73 16.13 11.51 3.53 1.53 14 600 5 19 2 886 30 −94 0.65 10.96 6.73 0.19 1.78 15 725 5 18 2 920 31 −79 0.31 9.18 5.21 0.24 1.35 16 500 5 37 2 150 15 140 0.46 21.97 18.55 6.44 1.81 17 1 150 5 15 2 823 17 −1 0.30 7.78 4.43 2.26 1.08 18 1 150 5 19 2 811 21 −112 0.58 11.59 7.45 2.98 1.50 19 550 5 25 2 813 22 −121 0.40 18.16 9.76 4.01 1.83 20 775 5 25 2 789 19 −91 0.60 15.00 11.39 4.22 1.56 21 1 525 5 26 2 784 12 −172 0.67 13.28 11.74 4.46 1.50 22 1 625 5 28 2 778 30 −131 0.81 14.97 12.52 4.22 1.32 23 1 350 5 32 2 716 36 96 0.82 16.64 10.69 5.39 1.57 Table 1. Basic information on sample plots of pure plantation of Larix principis-rupprechtii
-
调查林分的郁闭度、密度,然后对各样地实测全部林木的树高、胸径、冠幅及枝下高等。用相对树高表征树木优势度[16]。优势度计算公式为:
式中
$\overline H $ 为样地平均树高,Hi为每株解析木树高,Hmax为各样地内最大树高,HR为优势度值。然后按Kraft树冠优势度分级标准[17],将各样地内的林木个体优势度分为3级:I代表优势木、Ⅱ代表平均木、Ⅲ代表被压木,每个样地至少选择优势木和平均木各2株及被压木1株,进行解析木调查(共116株)。
-
上外包线法可以剥离出各单一因子影响并分析各单一因子影响以确定响应函数类型[18]。建立模型过程:(1)将林龄对应的各树高值(H)除以理论最大树高Hmax(33 m),得到相对树高[3];把相对树高变化范围分成若干区段,在每段中选取大于本段数据平均值加一倍标准差的数据点[19],或在一些数据偏少的区段中直接选择最大数据点,拟合外包线f(x1);(2)同理,将步骤1中消除了林龄因子影响的数据与下一个影响因子做散点图,并得到响应函数f(x2);(3)用同样方法逐个消除其他因子的影响,直到确定相对树高对最后一个因子的响应函数f(xn);(4)连乘响应函数,得到受多因素影响的树高生长耦合模型:
式中,H为单木树高,
$f({x_1})$ 、$f({x_2})$ 、...等分别表示单木树高生长对各单一因素${x_1}$ 、${x_2}$ 、...等的响应函数,n为考虑的影响因子数。把116株解析木数据分成:92株解析木实测数据拟合模型参数,剩余24株实测解析木数据验证模型拟合优度,选择验证解析木数据时注意了保持解析木在林龄、密度、郁闭度、海拔、坡度、坡向、优势度范围内的均匀分布。评价模型拟合优度的指标包括决定系数
${R^2}$ 、均方根误差RMSE、总相对误差TRE。R2越接近1,RMSE和TRE越小,说明模型拟合越好。计算公式为:式中,
${Y_i}$ 和${\hat Y_i}$ 为第i个单木树高的实测值和预测值,$\bar Y$ 为所有实测树高的平均值,N为用于参数率定或检验的样本数。 -
用统计Excel软件处理野外调查数据,用统计软件Origin进行绘图及曲线拟合,并用SPSS25.0软件进行典型相关分析等统计分析,用软件1st0pt进行耦合模型的参数拟合。
1.1. 研究区概况
1.2. 研究方法
1.2.1. 样地布设
1.2.2. 树木调查
1.2.3. 多因素影响的树高生长耦合模型建立
1.2.4. 数据处理
-
单木树高与各影响因子的相关分析表明(表2),各因子的相关系数排序为:林龄 > 海拔 > 郁闭度 > 优势度 > 林分密度 > 坡向 > 坡度。
影响因子
Influencing factors海拔
Elevation坡向
Slope aspect坡度
Slope林分密度
Stand density郁闭度
Canopy density林龄
Tree age优势度
Dominance相关系数
Correlation coefficients−0.421** 0.201 0.101 −0.248 0.419** 0.829** 0.413** 注:** 表示在0.01水平上显著相关。
Note: ** indicate a significant correlation at the level of 0.01 .Table 2. Correlation between tree height and influencing factors
-
基于所有解析木数据分析得到的单木树高对林龄的响应见图1-A中的上外包线,单木相对树高随林龄增加而逐渐升高,在21 a以前增长较快,随后增速变缓,在45 a时最高可达27.68 m。确定树高生长过程符合Richard方程:
$f\left( {age} \right) = a \times (1 - \exp ( - b \times age))^{\wedge}c$ ,R2 = 0.997 5。单木树高对海拔的响应见图1-B,树高在海拔2 000~2 200 m随海拔升高逐渐变大,之后随海拔继续升高而逐渐下降。树高对海拔的响应呈三次多项式:$f\left( {ele} \right) = d \times el{e^3} + e \times el{e^2} + f \times ele + g$ ,R2 = 0.856 7。单木树高对郁闭度的响应见图1-C中的外包线,单木树高在郁闭度0.3~0.56范围内随郁闭度增大而逐渐增大,在郁闭度 > 0.56后逐渐减小。树高对郁闭度的响应呈二次多项式:$f\left( {cd} \right) = $ $ h \times c{d^2} + i \times cd + j$ ,R2 = 0.662 2。单木树高对优势度的响应特征见图1-D中的外包线,可知单株树高随优势度增加先平缓增加,在优势度为−0.2后增速加快,在优势度为0.08后增速渐趋平缓。树高对优势度的响应符合S型曲线:$f\left( {dom} \right) = k + $ $ (l - k)/ $ $ \left[ {1 + \exp \left( {\left( {dom - m} \right)/n} \right)} \right]$ ,R2 = 0.942 9。单木树高对林分密度的响应见图1-E中的外包线,可知树高随密度增加在400~1 200株·hm−2范围内逐渐增加,之后逐渐下降。单木树高对林分密度的响应呈二次多项式关系:$f\left( {den} \right) = o \times $ $ de{n^2} + p \times den + q$ ,R2 = 0.835 6。单木树高对坡向的响应见图1-F中的外包线,可知最适坡向为阴坡和半阴坡,在坡向180°到0°的范围内,单木树高随坡向靠近正北方向逐渐增加。响应函数为二次多项式:$f\left( {asp} \right) = r \times as{p^2} + $ $ s \times asp + t$ ,R2 = 0.522 4。单木树高对坡度的响应见图1-G中外包线,单木树高在坡度达23.5°之前随坡度增大而增加,之后转而减小,最适坡度在20°~27°。响应函数为:$f\left( {slope} \right) = u \times slop{e^2} + v \times slope + w$ ,R2 = 0.649 7。 -
按相关系数的大小排序,从单木树高响应林龄的模型开始,逐个增加耦合其他因子,得出模型(式6~式12)。为便于生产应用,本研究同时建立包括林龄、林分密度、海拔、坡向、坡度这些独立变量的模型(式13),然后在此基础上增加优势度建立了模型式(14)。从表3可看出,从式(6)到式(10),伴随R2升高,RMSE和TRE降低;式(11)为加入坡向后R2降低,RMSE和TRE均上升;式(12)综合考虑了所有因子,其R2升高;式(13)只考虑了独立因子,所以R2下降;式(14)在式(13)的基础上加入了优势度,导致R2升高到和式(12)相近的水平。式(6)~(14)的TRE不超过 ± 0.3%,式(10)以及式(12)的TRE最小,这3个评价指标均证明,式(12)是最优的多因子耦合模型,其次为式(14)。
编号 NO. 模型 Model R2 RMSE/m TRE/% 式(6) f (1)=f (age)=33×(1−exp(−0.029×age))1.41 0.871 2.036 0.228 式(7) f (2)=f (age)×f (ele)=(33×(1−exp(−0.028×age))1.293)×(3.34E−9×ele3−2.58E−5×ele2+0.065×ele−53.6) 0.908 1.717 −0.212 式(8) f (3)=f (age)×f (ele)×f (cd)=(33×(1−exp(−0.025×age))1.275)×(−3.85E−6×ele3+0.03×ele2−72.64×ele+58 516)×
(5E−4×cd2−9.6E−4×cd−6.3E−4)0.911 1.697 −0.200 式(9) f (4)=f (age)×f (ele)×f (cd)×f (dom)=(33×(1−exp(−0.041×age))1.43)×(0.07×ele3−569.14×ele2+1.45E+6×ele−
1.19E+9)×(−9.26E−5×cd2+1.45E−4×cd+0.0001)×(4.41E−5+(2.7E−4)/(1+exp((dom+0.198)/−0.245)))0.961 1.117 −0.072 式(10) f (5)=f (age)×f (ele)×f (cd)×f (dom)×f (den)=(33×(1−exp(−0.044×age))1.46)×(2.02E−16×ele3−1.57E−12×ele2+
4E−9×ele−3.32E−6)×(−1.95×cd2+3.32×cd+1.02)×(567.36+(−1025.4)/(1+exp((dom−7.67)/40.33)))×(0.13×
den2−364.3×den+k/>1.18E+60.962 1.115 0.019 式(11) f (6)=f (age)×f (ele)×f (cd)×f (dom)×f (den)×f (asp)=(33×(1−exp(−0.041×age))1.431)×(−0.65×ele3+643.6×ele2+
5.47E+6×ele+4.701E+7)×(−3.62E−18×cd2+6.57E−18×cd+1.36E−18)×(3.53+(1.22)/(1+exp((dom+0.016)/
−0.001)))×(3.67E−5×den2−0.252×den+1 414.2)×(−0.02×asp2−0.497×asp+4 918.93)0.941 1.380 0.104 式(12) f (7)=f (age)×f (ele)×f (cd)×f (dom)×f (den)×f (asp)×f (slope)=(33×(1−exp(−0.043×age))1.462)×(0.049×ele3−
362.01×ele2+8.86E+5×ele−6.99E+8)×(−0.101×cd2+0.198×cd+0.038)×(−314.96+578/(1+exp((dom+
61.81)/−338)))×(2.57E−16×den2−5.56E−13×den−6.64E−10)×(−0.0034×asp2−0.177×asp−1927)×
(−0.001×slope2+0.026×slope+0.396)0.975 0.906 −0.028 式(13) f (8)=f (age)×f (ele)×f (den)×f (asp)×f (slope)=(33×(1−exp(−0.022×age))1.27)×(0.039×ele3−293.6×ele2+
7.22E+5×ele−5.7E+8)×(−3E−19×den2+7.67E−16×den+6.82E−13)×(−0.007×asp2−0.3×asp+6128.7)×
(−0.009×slope2+0.37×slope+10.31)0.916 1.650 −0.156 式(14) f (9)=f (age)×f (ele)×f (dom)×f (den)×f (asp)×f (slope)=(33×(1−exp(−0.0417×age))1.46)×(1.51E−8×ele3−
0.0001×ele2+0.29×ele−237)×(3.07+(−2.48)/(1+exp((dom+0.11)/0.24)))×(0.28×den2−767.14×den−716549)×
(1.09E−15×asp2+3.36E−14×asp+5.98E−10)×(0.111×slope2−4.186×slope−97.64)0.971 0.960 −0.034 Table 3. Forecasting and calculating models of tree height growth per plant
用24株解析木的实测数据对模型(6~14)进行检验(如图2):式(10)、式(12)、式(14)比其他模型拟合的要好,预测点最靠近并均匀分布在45°线两侧。