本研究在宁夏六盘山南段东坡的香水河小流域（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度以内的坡向为负。

调查林分的郁闭度、密度，然后对各样地实测全部林木的树高、胸径、冠幅及枝下高等。用相对树高表征树木优势度^[16]。优势度计算公式为：

式中$\overline H $为样地平均树高，H_i为每株解析木树高，H_max为各样地内最大树高，H_R为优势度值。

然后按Kraft树冠优势度分级标准^[17]，将各样地内的林木个体优势度分为3级：I代表优势木、Ⅱ代表平均木、Ⅲ代表被压木，每个样地至少选择优势木和平均木各2株及被压木1株，进行解析木调查（共116株）。

上外包线法可以剥离出各单一因子影响并分析各单一因子影响以确定响应函数类型^[18]。建立模型过程：（1）将林龄对应的各树高值（H）除以理论最大树高H_max（33 m），得到相对树高^[3]；把相对树高变化范围分成若干区段，在每段中选取大于本段数据平均值加一倍标准差的数据点^[19]，或在一些数据偏少的区段中直接选择最大数据点，拟合外包线f(x₁)；（2）同理，将步骤1中消除了林龄因子影响的数据与下一个影响因子做散点图，并得到响应函数f(x₂)；（3）用同样方法逐个消除其他因子的影响，直到确定相对树高对最后一个因子的响应函数f(x_n)；（4）连乘响应函数，得到受多因素影响的树高生长耦合模型：

式中，H为单木树高，$f({x_1})$、$f({x_2})$、...等分别表示单木树高生长对各单一因素${x_1}$、${x_2}$、...等的响应函数，n为考虑的影响因子数。

把116株解析木数据分成：92株解析木实测数据拟合模型参数，剩余24株实测解析木数据验证模型拟合优度，选择验证解析木数据时注意了保持解析木在林龄、密度、郁闭度、海拔、坡度、坡向、优势度范围内的均匀分布。评价模型拟合优度的指标包括决定系数${R^2}$、均方根误差RMSE、总相对误差TRE。R²越接近1，RMSE和TRE越小，说明模型拟合越好。计算公式为:

式中，${Y_i}$和${\hat Y_i}$为第i个单木树高的实测值和预测值，$\bar Y$为所有实测树高的平均值，N为用于参数率定或检验的样本数。

用统计Excel软件处理野外调查数据，用统计软件Origin进行绘图及曲线拟合，并用SPSS25.0软件进行典型相关分析等统计分析，用软件1st0pt进行耦合模型的参数拟合。

单木树高与各影响因子的相关分析表明（表2），各因子的相关系数排序为：林龄 > 海拔 > 郁闭度 > 优势度 > 林分密度 > 坡向 > 坡度。

基于所有解析木数据分析得到的单木树高对林龄的响应见图1-A中的上外包线，单木相对树高随林龄增加而逐渐升高，在21 a以前增长较快，随后增速变缓，在45 a时最高可达27.68 m。确定树高生长过程符合Richard方程：$f\left( {age} \right) = a \times (1 - \exp ( - b \times age))^{\wedge}c$，R² = 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$，R² = 0.856 7。单木树高对郁闭度的响应见图1-C中的外包线，单木树高在郁闭度0.3～0.56范围内随郁闭度增大而逐渐增大，在郁闭度 > 0.56后逐渐减小。树高对郁闭度的响应呈二次多项式：$f\left( {cd} \right) = $$ h \times c{d^2} + i \times cd + j$，R² = 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]$，R² = 0.942 9。单木树高对林分密度的响应见图1-E中的外包线，可知树高随密度增加在400～1 200株·hm⁻²范围内逐渐增加，之后逐渐下降。单木树高对林分密度的响应呈二次多项式关系：$f\left( {den} \right) = o \times $$ de{n^2} + p \times den + q$，R² = 0.835 6。单木树高对坡向的响应见图1-F中的外包线，可知最适坡向为阴坡和半阴坡，在坡向180°到0°的范围内，单木树高随坡向靠近正北方向逐渐增加。响应函数为二次多项式：$f\left( {asp} \right) = r \times as{p^2} + $$ s \times asp + t$，R² = 0.522 4。单木树高对坡度的响应见图1-G中外包线，单木树高在坡度达23.5°之前随坡度增大而增加，之后转而减小，最适坡度在20°～27°。响应函数为：$f\left( {slope} \right) = u \times slop{e^2} + v \times slope + w$，R² = 0.649 7。

图  1  单木树高对立地条件和林分特征的响应

Figure 1.  Response of single tree height to site conditions and stand structure characteristics

按相关系数的大小排序，从单木树高响应林龄的模型开始，逐个增加耦合其他因子，得出模型（式6～式12）。为便于生产应用，本研究同时建立包括林龄、林分密度、海拔、坡向、坡度这些独立变量的模型（式13），然后在此基础上增加优势度建立了模型式（14）。从表3可看出，从式（6）到式（10），伴随R²升高，RMSE和TRE降低；式（11）为加入坡向后R²降低，RMSE和TRE均上升；式（12）综合考虑了所有因子，其R²升高；式（13）只考虑了独立因子，所以R²下降；式（14）在式（13）的基础上加入了优势度，导致R²升高到和式（12）相近的水平。式（6）～（14）的TRE不超过 ± 0.3%，式（10）以及式（12）的TRE最小，这3个评价指标均证明，式（12）是最优的多因子耦合模型，其次为式（14）。

用24株解析木的实测数据对模型（6～14）进行检验（如图2）：式（10）、式（12）、式（14）比其他模型拟合的要好，预测点最靠近并均匀分布在45°线两侧。

图  2  单木树高生长模型（式6～14）的检验

Figure 2.  Test of tree heights of models of Eq. 6-14

林龄是影响树高生长的最重要因子。华北落叶松树高随林龄增加的生长过程表现为“S”型曲线，可划分为生长初期、快速生长期、平缓生长期3个阶段^[20]。如图1-A中下外包线表示的树高生长过程S型曲线特征比上外包线更明显，本研究表明采用Richard方程能很好表示树高生长过程^[21]。林冠郁闭度直接或间接地影响林木生长及众多服务功能^[22]，本研究中，华北落叶松单木树高在郁闭度为0.49～0.64时最好，相比多功能森林经营要求的合理郁闭度范围0.6～0.8^[23]有些偏低但也差异不大。优势度能直接反映林木对自然资源的利用能力^[24]，单木树高随优势度增大过程呈“S”型曲线特征，在优势度超过0.08后，树高增大速率趋于平缓，说明优势木和亚优势木的树高可能会接近或达到立地质量决定的潜在树高，是获得木材生产经济效益的主要贡献者。林分密度会影响单木的光照、水分、养分等资源占有量^[6]，华北落叶松单木树高生长的最适密度平均为1 100～1 300株·hm^−2[25]。海拔虽不直接影响树木生长，但会通过影响降水、温度、蒸散等环境因子间接影响树木生长^[26-27]。本研究表明六盘山区最适合华北落叶松生长的海拔范围为2000～2400 m，降水不足是六盘山区树木生长的主要限制因子，在海拔超过2200 m后，海拔升高导致的气温降低逐渐成为限制树木生长的主要因素，使得树高随海拔增加而逐渐降低^[3]。坡向与坡度会影响地面得到的太阳辐射，阳坡接受的太阳照射比阴坡更多，坡度大时接受太阳垂直照射的程度越大^[28]，从而影响潜在蒸散和水分条件^[3]。本研究表明，单木树高在阴坡和半阴坡最大，这是因光照较弱、蒸散较小导致植物可用水分较多和受干旱胁迫较轻；单木树高随坡度增大表现为先增后减，最适坡度范围为20°～25°，类似于河北塞罕坝林场的华北落叶松生长研究结果^[27]，因地势平缓时导致的林木相互遮光明显以及坡度过陡导致的土壤偏薄和水土流失均不利于林木生长。

在自然环境中，林木生长同时受多因子影响，因而时空差异很大，所以综合考虑多因子影响的模型预测精度一般较高。然而，由于缺少调查资料或追求简单易用等原因，以往的华北落叶松生长模型常仅考虑林龄影响^[8]，忽略其他因子影响，降低了预测精度和限制了应用范围。相比之下，建立多因素耦合模型，便于准确预测在不同立地及林分结构条件下的生长变化。在以往的六盘山华北落叶松林分生长指标模型中考虑了林龄、密度和海拔这3个主要因素^[3]，在祁连山青海云杉林分生长指标模型中考虑了林龄、密度、海拔和坡向^[10]。在本研究中，经过检验比较模型（6～14）得出，六盘山华北落叶松单木树高生长模型在综合考虑因子的影响时精度最佳，优点是能考虑所有立地条件和林分特征因子对树高生长的影响，会更准确地预测树高。模型在考虑林龄时能准确预测立地条件和林分结构平均情况下的树木生长；考虑海拔的影响可反映树高生长随海拔变化的空间差异；继续增加考虑了郁闭度和林分密度的影响，能在一定程度上反映林分结构对树木生长的影响，利于以林分为空间单元预测单木树高的时空变化；增加优势度的影响能精确到单株精准经营；继续增加考虑坡向和坡度，因为本研究区位于半湿润区，树高的坡向响应较弱，且本研究选取不同坡向的样地较少，无法准确体现坡向对树高生长的影响，在未来研究中应增加所缺坡向的样地，进一步研究不同坡向下的树高生长差异，使模型更精准。综上，可根据研究目的或生产需求，从这里建立的单木树高耦合模型（式6～14）中选择应用。

在宁夏六盘山半湿润区，通过华北落叶松人工纯林样地及解析木调查，研究了单木树高受不同立地因子及林分结构特征影响的规律，得到如下结论：

（1）单木树高生长受立地条件及林分结构因子影响的大小排序为：林龄 > 海拔 > 郁闭度 > 优势度 > 林分密度 > 坡向 > 坡度。

（2）单木树高生长随林龄增加和优势度增大均表现为“S”型；单木树高生长随海拔、郁闭度、林分密度、坡向、坡度的增加均呈现“先增后降”的变化。

（3）建立了能反映多因素影响的耦合模型，其中综合考虑立地因子和林分结构特征影响的模型表现最好。可根据研究需要，或考虑生产应用方便，选择利用本研究建立的单木树高生长耦合模型，从而为华北落叶松人工林的精准经营及管理决策提供科学依据。