近年来，国内外的林分生长研究中正逐渐由全林模型向林分径阶分布模型和单木模型过渡。单木模型包括组成林分的林木生长随机过程和枯损发生概率估计, 特别是它们寿命和枯损速率差别, 并由此产生了生态学中树种演替现象。由于林木间竞争因子构造条件不同，所以单木模型又分为有距离和无距离的单木模型。从微观和实用的角度，用无距离的竞争因子构造和估计森林的生长和演替进程构成本次天然混交林生长和演替模拟的基本内容。

1 研究资料

天然红松(Pinus koraiensis)林资料取自白云庆等(1987)上世纪50—60年代设置200块临时标准地；天然云杉(Picea asperata)、冷杉(Abies nephrolepis)及天然软阔林(natural softwood)资料取自黑龙江、吉林临时标准地；天然落叶松(Larix gmelinii)资料取自1986年在塔河设置的200块临时标准地；林木的胸径生长和枯损资料来自牡丹江等500块固定样地(有3次复测)。

2 天然混交林的随机生长与演替基础模型 2.1 标准地位指数函数 2.1.1 地位指数函数

根据对地位指数函数的分析(郎奎建，1999)，构造关于地位指数S和t的二元函数：

(1)

称为地位指数函数。当基准年龄为50年时，则称ψ(s_b, t)为标准地位指数函数。利用三次多项式对(1)式中地位指数函数曲线簇中的3个参数进行拟合，获得东北各主要林分的地位指数函数曲线簇，见表 1。

2.1.2 Zeide地位指数函数

为克服地位指数表因基准年龄不同而缺乏可比性的缺点，Zeide(1977)提出一种跨国家、跨树种的相对优势树高H′_i(t)=h_i(t)/h_i(50)的标准生长曲线簇。按基准年龄120年，将H′_i(t)从100%~260%分成1, 2, …, 18级，该相对地位指数级称为Zeide指数，简记Z。假设Zeide各相对树高生长曲线函数为H′_i(t)=b₁×[1-exp(b₂t)]^b₃。按前述办法获得Zeide指数函数，见表 2。

以上定义了3种地位指数Z、S和S_b的函数, 它们均是在同形条件下的地位指数。

2.1.3 地位指数函数解析法二点法

由某林分测得任意两个相对优势树高生长的相关点(t₁, H′₁)和(t₂, H′₂), 考察比值R=H₂^′/H₁^′(t₁ < 50, t₂>50)，由于在t₁、t₂条件下, 对应Zeide指数相对树高生长曲线H_i^′=φ(t)的比值是唯一的, 且与Z呈单调递增，用迭代法可求得该林分与R相对应Z值，显然在约束Z∈[1, 18]条件下上述迭代一定收敛，利用H′(t)和H(t)转换关系，将指定基准年龄代入H(t)中, 即可获得相应地位指数S和S_b。一点法：设某林分已测得一个树高生长的相关点(t, H), 用二分法解一元方程H=ψ_b(s_b, t)，再由已得S_b求出S和Z。

2.2 树木胸径生长方程随机模拟 2.2.1 胸径生长方程的导出

在有复测的固定样地条件下可获得胸径生长相关数列y_t, y_t-5, 利用单分子式的差分原理(白云庆等，1987)，可直接导出胸径的生长方程：将固定样地复测样木胸径化为相差1年的相关数列y_t, y_t-1，对差分方程按y_t=a×y_t-1+b回归，参数必须满足a∈(0, 1), b>0，导出生长方程为y=m×[1-exp(-kt)]，式中:k=-ln(a), m=b/(1-a)。见表 3。

2.2.2 树木胸径生长历程函数

上述方程是平均地位和平均密度的胸径生长方程。用于某一林分的胸径生长还需引进一定的噪声，从而实现胸径生长方程随机模拟。

胸径生长历程方程假设在树木生长Mitscherlish速度方程上迭加一速度噪声，解一阶线性非齐次微分方程，由此获得胸径生长历程函数(白云庆等, 1987)：

用蒙特卡罗(Monte-Carlo)法产生白噪声(王梓坤，1978), 假设在某时间t由伪随机数函数RND(x)产生n个服从在[0, 1]区间上均匀分布的随机数为R₁，R₂，…R_n。根据中心极限定理可知，随着n的增大，其和趋于正态分布：。作线性变。从而由伪随机数产生在某t时树木生长速度白噪声(指定方差D)。由此生成独立的随机变数序列f(t₁), f(t₂), f(t₃), ……。则树木胸径生长历程函数:

(2)

2.3 林分枯损模型 2.3.1 树冠竞争指数

由疏开木冠幅方程，式中:c₁=0.7~1.2, c₂=14.7~24.7, c₃=0.025，随树种而异，则树种无距离竞争指数可用1 hm²林地按疏开木冠幅总量来定义：，式中：m为径阶数; n_i为1 hm²林地径阶株数; c_wi为径阶平均冠幅(疏开木)。乘100是为去小数位。当CCF=100时为疏开木林分, 当CCF=600时为近饱和林分。

2.3.2 树种耐荫性指数

在生态学中, 树种耐荫性是一种定性指标，可由冠高比C_r和CCF来定量表征, 则p₁=ln(c_r/100)/ln(CCF/100)称为树种耐荫性指数。一般阳性树种p₁=0.7、强荫性树种p₁=0.4。

2.3.3 林分的枯损模型

枯损函数是指某林龄t时(注意t用CCF来表示), 林分在单位面积上的定期正常枯损株数按径阶的分布函数, 记为M₀(t, d)。从固定样地定期复测资料中可获得M₀(t, d)的列表的表达式。显然它与CCF、p₁和径阶有关。为构建林分的枯损模型, 将M₀(t, d)分解为枯损率强度函数和径阶枯损分布函数的乘积(郎奎建，1996)。

枯损率强度是指某树种(species, 记SP)林分的枯损株数与林分总株数之比。根据枯损特性, 可假设枯损率强度是关于CCF的幂函数p₂=a₂(sp)·CCF^b₂(sp)。其参数见表 4。

由于枯损分布是生物寿命的逆分布, 概率论有一种描述寿命分布的概率模型就是指数分布：f(x)=λ·exp[-λ(x-d₀)/d_h], x≥d₀, d₀起测径阶、d_h径阶距。由固定样地资料可获得林分各径阶的活立木株数, 枯损木株数的分布, 则成活率=1-枯损率, 把它看成随机变数(径阶x)出现的频数。用矩法估计指数分布参数λ=(m×σ)^-1/2 (取它们的几何平均数)。由此获得各树种指数分布参数λ：红落樟0.3055、云冷杉0.3055、水胡黄、0.2703281、硬阔0.2020272、软阔0.32337276。

2.4 在林分初始胸径分布未知条件下，林分状态结构的Weibull分布模拟 2.4.1 林分胸径的Weibull分布概率模型

设d₀起测径阶的下限，具有概率密度函数

(3)

的连续随机变量称为服从参数b, c的Weibull分布。

根据森林分子概念, 不管林分胸径状态初始分布如何, 经过长期自然生长(包括枯损), 使得林分状态结构概率分布趋于正态分布，也就是说正态分布可以作为林分状态结构的极限分布。此时Weibull分布参数c趋于3.6。

2.4.2 林分状态结构的动态Weibull分布模拟

Weibull分布参数c的时间方程:利用大量固定样地资料，由Weibull分布对子样d₁, d₂, …, d_n的似然函数，用牛顿迭代法求出它的最大似然估计值c(白云庆等，1987)。可获得相关资料(t_i, c_i), i=1, 2, …, 采用单分子式:

(4)

在已知参数α=3.6下的迭代回归法可获得式(4)中参数β, γ的估计值, 具体见表 5。

Weibull分布参数b时间方程导出：在已知c(t)和林分平均胸径D_g的条件下, 由Weibull分布二阶原点矩m₂=Γ|2/c+1|b²和林分平均胸径定义：m₂=D_g²可得:

(5)

林分状态结构函数拟合：假设对该林分胸径分布知之不多，已知林分各树种每公顷株数N及平均胸径D_g，设:d₀为起测径阶；N₀为上径阶累计株数。则由Weibull分布函数, x径阶株数为:

(6)

由于公式(6)求出胸径分布f_N(x)与预定的D_g不相容，所以需用比值法调整f_N(x)。

2.5 其他基础模型

其他基础模型包括林分的进界率模型和林分优势高、林分平均高的对偶回归模型等。

3 研究结果 3.1 模拟初始广义林分结构向量数据文件 3.1.1 初始林分的基础数据文件结构设计

初始林分的基础数据文件编辑分两个水平：不完全的初始胸径结构文件(即实际林分的组成树种、地位指数、林分年龄、平均胸径、平均树高、公顷株数、起源等)和完全的胸径结构文件。

3.1.2 现实林分地位确定

可用一点法或两点法确定该林分(即优势树种)的地位指数S、S_b、Z。

3.1.3 组成树种的广义林分结构向量拟合

林分密度和地位的阻尼因子:由于前述的基础模型都是平均密度、地位意义上的方程, 为了使其适用于某一林分, 定义密度阻尼因子: CF=(CCF/300)^0.3，式中300是CCF的平均水平, 幂函数指数0.3对CCF作用起到某种阻尼效果。同理构造地位的阻尼因子SW。

组成树种上层树高H_top拟合:由于每个林分只能有一个优势高生长过程，所以当组成树种平均年龄t已知时，可利用H(1)(优势树种导向曲线在50 a的条件平均高)和H(i)(组成树种导向曲线在50 a的条件平均高)的比值，求出某组成树种的相对地位si=S(优)H(i)/H(1)，则该组成树种的地位指数函数H_tpo(t)=a₁(si)+b₁(si)exp, 代入t即可求得H_tpo; 当组成树种平均年龄t未知时，可利用该树种的Schumacker树高曲线求该树种上层树高H_top=a₃exp(-b₃/D_g)×1.2CF，系数1.2是平均高向优势高转换的系数。

组成树种平均年龄的拟合:由该树种地位指数函数即公式(1)导出树种的平均年龄

林分胸径分布的拟合：用公式(3)、(4)、(5)进行Weibull分布拟合。

林分树高分布的拟合：由该树种的Schumacke树高曲线乘CF可求出各径阶树高。

林分各径阶冠高分布：由树种耐荫指数可导出各径阶冠高C_W=径阶平均高×CFF^p₁。

由此形成t时刻完备的各组成树种(包括优势树种)的胸径、树高、优势高、年龄、冠幅和径级株数等，称为广义林分结构向量。

3.2 林分向后生长与演替模型

在基础模型的支持下，由现实林分年龄起, 用预测方式模拟它未来若干年林分按径级(广义林分径级结构向量)的生长(演替)过程称之为林分的向后生长与演替模型。

3.2.1 现实林分的生长期调整因子

由于现实林分的年龄不是设计的龄阶, 所以引进生长期调整因子。由RA=(TI-AGE)/REPH(式中:TI为最近设计的龄阶；AGE为现实林分林龄；REPH为龄阶距)，将林分广义林分结构向量调整到指定龄阶。

3.2.2 林木枯损的随机发生

由广义林分结构向量和各树种枯损模型，可确定间隔期内各径级的枯损株数。由于有的径级枯死株数不够1株，所以按累计枯死株数为1确定。在林木空间分布为随机分布的条件下，由随机数确定枯死的林木，并由此模拟林木的分布格局。

3.2.3 树种进界株数

由于树种更新和幼树生长, 在模拟的林分中会发生进界现象。为正确控制和估计树种演替过程, 假设⑴在林分表中的第1个树种为目的树种, 由它确立林分进界演替方向；⑵进界现象仅在目的树种的林分状态向量中进行；⑶在期间进界株数全部进入指定的起测径阶。假设⑵比较弱些，因为很难获得任意树种组合的树种更新株数比重参数，故仅适用中国东北有明显目的树种的针阔或针叶、阔叶混交林地区。经过对大量固定样地进界率的因子分析, 发现它是林分竞争因子CCF减函数。则由树种平均进界率(JYJI)得出REPH期间进界株数YN=JYJI·N·REPH·RA。

3.2.4 树种各径级的胸径随机生长

由前述的胸径生长历程函数即公式(2)，通过引入林分密度和地位阻尼因子，将其调整到现实密度和地位的生长，并引入径级阻尼因子DOT=(D/D_g)^0.3，将其调整到指定径级的胸径生长。

3.2.5 树种各径级的树高生长

同理，通过引入林分密度阻尼因子和径级阻尼因子，将其调整到现实密度、径级的树高生长。

3.2.6 树种径阶年龄分布

树种径阶初始年龄分布根据地位指数曲线，由径阶树高求径阶年龄A₁，由胸径生长方程求径阶年龄A₂，取平均数。在天然林生长模拟中, 各径阶年龄分布随着林木生长和枯损在发生变化, 由此导致林分平均年龄发生变化。从模拟结果证实可能出现的林龄转移滞留现象(郎奎建，2002)。

3.2.7 混交林中树种重叠效应

上述方法的实质是按纯林模式将各树种的林木生长、枯损, 经过密度、地位、径阶阻尼因子调整, 叠加形成林分结构的动态转移。这里忽略了混交林中树种间的重叠效应，即位于上层林冠的树木要压抑下层林冠树木的生长，树种间的重叠效应主要发生在针阔混交林中, 在幼年期, 针叶树往往是被压抑的对象。

为了真正实现混交林的生长模拟, 构造树种间重叠效应系数UF=(p₁+0.3)^2.17，式中:p₁为树种耐荫指数(针叶树p₁=0.4、阔叶树p₁=0.7)。公式仅对针阔混交林中位于主林层下的树种林木起作用, 它将压抑该林木的胸径生长，一旦该树种上升到主林层, 这种效应消失。参数取自对固定样地中纯林和混交林胸径生长速度分析。

3.2.8 广义林分径级结构向量的调整

以上是根据期初广义林分结构向量和基础模型，导出期末的广义林分径级结构向量(类似于无距离的单木模型的数据结构)，有可能出现因子间的某些不相容，所以采用软件的智能化方法对广义林分径级结构向量的进行逐级控制调整。

3.2.9 林分向前生长与演替模型

从起始龄级起直至现实林分年龄止, 用模拟方式重现它过去若干年林分生长与演替过程称之为林分的向前生长与演替模型。

3.2.10 森林生长三维动态图像

按标准地面积0.01 hm²，林木株数按1%取。

3.3 天然混交林生长与演替过程实例

东北林业大学凉水实验林场第27号固定标准地, 该标准地原生林型为椴树红松林, 1956年进行过皆伐和红松人促更新, 生成当今阔叶林(1红7杨2白，林分平均年龄27年)。经林分向前、向后(间隔10年)的生长与演替模拟，得该林分现实生长过程见表 6。

3.4 区域性森林资源动态予测

常规的区域性森林资源动态预估需用标准地法估计参数。由于资源预估与该区的地位、密度分布等有关，标准地不可能有这样大的代表性，而且枯损量不好估计，特别负的纯生长量几乎估计不出来。本文利用森林面积蓄积统计表，把每个龄组看成一个林分，按上述方法模拟(包括相应的参数文件)，汇总求出区域性资源预估所需部分参数(表 7)。

4 结论

从胸径广义结构角度, 由某些林学规律概括(包括正态分布可以作为分状态结构的极限分布)，研究和构造由一簇包含林分地位、生长、枯损、密度、进界等基础模型。该基础模型55个参数作为本系统的参数文件，它决定该地区天然混交林生长与演替方向和速率。

在基础模型构造中，采用蒙特卡洛法产生胸径生长速度的白噪声和林木枯死的随机发生等，由此获得林分非完全的随机生长历程。

建立密度、地位、径级组尼因子，实现基础生长模型向现实密度、地位、径级生长模型的非线性过渡。

定义混交林中树种重叠效应系数和树种的相对地位，解决将按纯林模型有机的复合，成功地实现混交林中树种间的模型叠合。

以不完全的初始胸径结构文件(一般样地文件)分3个层次(林分结构层次、各组成树种的广义林分结构向量和广义林分径级结构向量层次)进行天然混交林的向前、向后的生长与演替模拟，并采用软件的智能化技术对广义林分径级结构向量的进行逐级控制调整，从而实现与距离无关的单木模型的模拟。在数据文件和竞争指数上稍加改进，即可完成向与距离有关的单木模型的过渡。