传统的林分择伐优化模型是功能优化，如总收益最多(谢哲根等，1994)、纯收益最多(Hof et al., 2000)或净现值最大(Buongiorno et al., 1995)等，很少把林分结构作为优化目标。事实上，根据系统论结构决定功能的观点(杨春时等，1987)，只有保持结构优良，系统功能才能很好地发挥。研究表明，增加林分结构多样性，将提高物种多样性和生态系统稳定性(Pretzsch, 1999)。因此，维持和提高林分结构复杂性已成为生态系统和景观保护的经营目标(Garcia et al., 1999)。

林分结构包括空间结构和非空间结构。空间结构包括混交、大小分化和林木空间分布格局3个方面(惠刚盈等，2001)。非空间结构包括径级结构、生长量和树种多样性等。空间结构直接影响非空间结构，因此空间结构比非空间结构更重要。择伐是对实施经营林地最重要的干扰(冯佳多等，1998)。合理的择伐是调整林分空间结构的手段，如过于稠密将导致许多林分滋生昆虫和疾病(Bartuska，1999)，而通过抚育间伐可调整林分空间结构，改善林分卫生状况，防止病虫害发生与蔓延。研究还证实，择伐能提高树种多样性(Ammer et al., 1999)。因此，基于空间分析的单株树木择伐研究呈日益增长趋势(Tome et al., 1999)。

现有建立在空间结构分析基础上的林分优化经营模型(Pukkala，1999)主要是功能优化，甚至采用模拟林分(Hof et al., 2000)，空间结构仅作为约束且多属单方面的。这些模型的目标不是空间结构，也未考虑生物多样性保护的需要。本研究把多方面空间结构引入林分择伐规划，并以林分空间结构为目标，以非空间结构为主要约束，建立林分择伐空间优化模型，旨在通过空间结构优化，最大限度地保持林分结构完整性，持续发挥森林的多种功能。

1 林分择伐空间优化模型的建立 1.1 目标确定

林分择伐空间优化的实质是合理确定采伐木，以便在获取木材并保持非空间结构的同时，导向理想的空间结构。因此，林分择伐后保持理想的林分空间结构是模型的目标。什么样的空间结构是理想的，这个问题尚难回答。结合已有的研究成果，作一些假设是必要的。本研究从混交、竞争与分布格局3个方面描述林分空间结构，相应地，有3个空间结构子目标。首先对子目标作最优状态分析，然后综合得到模型目标。

择伐的基本单元是林分，在森林景观中就是一个斑块。景观生态学强调斑块与周围环境在外貌或性质上不同，但又具有空间可分辨的内部均质性(郭晋平，2001)。在Weintraub(1991)所提出的森林规划等级方法中，也强调单元同质性是经营单元划分的依据。显然，在森林景观中，要求斑块内部保持均质且与邻近斑块有着本质不同的特征。根据森林斑块的乔木层，不妨把斑块均质性理解为林分空间结构的均质性，包括混交、竞争和林木空间分布格局的均质性。

林木空间分布格局的均质性是林木均匀分布。Moeur等(1993)研究表明，林木最初与它们直接近邻者竞争，自稀疏增加了林木之间的距离，使林木分布格局从聚集变为均匀，结果大树趋于均匀分布，而年幼的林木呈聚集分布。Hanus等(1998)的研究证实了这一点。大树和林分整体均匀分布可减少林木之间冠层重叠，有效阻止病虫害蔓延和传播，还可以使地表连续覆盖，符合近自然林业要求。因此，林分择伐后，林木均匀分布，包括大树均匀分布和林分整体均匀分布，可作为林分空间结构的第1个子目标，即聚集指数越大越好。本研究把大于起测直径的树木分为两类，凡胸径大于或等于林分平均胸径的树木称为大树，否则为小树。

林分在混交方面的均质性是不同树种之间充分隔离。一般认为，混交林中树种相互隔离程度越高，林分越稳定。因此，林分择伐后保持较高的混交度是林分空间结构的第2个子目标，即混交度越大越好。林分在竞争方面的均质性是林分择伐后保持较低的竞争水平，使各保留木在满足其生态位需求上更有利。所以，确定林分采伐后保持较低的竞争指数为第3个子目标，即竞争指数越小越好。

由于林分空间结构的各个方面既相互依赖又可能相互排斥，要求各子目标同时达到最优是困难的，采用多目标规划可取得林分空间结构整体最优。基于以上分析，林分空间结构在混交度和聚集指数上都是以取大为优，而竞争指数是取小为优。同时，考虑保持空间结构在林分整体上的稳定性，避免出现变动过大。再采用乘除法综合各子目标就得到模型目标(钱颂迪，1990)。

1.2 约束设置

主要根据林分非空间结构设置约束条件，包括林分结构多样性、生态系统进展演替和采伐量不超过生长量约束。林分结构多样性是生物多样性的重要组成部分，林分结构多样性主要指林木大小多样性和树种多样性(Buongiorno et al., 1995)。生态系统进展演替是维持优势种或建群种的优势程度。采伐量不超过生长量是木材可持续利用的基本前提。同时，要求伐后空间结构质量不降低。这些约束都符合近自然林业和可持续经营的要求。

1.2.1 林木大小多样性约束

林木大小多样性用径级多样性和株数按径级倒J形分布描述。以采伐后不减少径级个数作为径级多样性约束条件；采伐后保持株数按径级的倒J形分布作为林木分布形式约束条件。径级多样性约束条件比较容易建立，倒J形分布约束稍复杂一些。

de Liocourt(1899)研究认为，理想的异龄林株数按径级依常量q值递减。此后，Meyer(1952)发现，异龄林株数按径级的分布可用负指数分布表示，公式如下

(1)

式中，N株数；e自然对数的底；D胸径；a、k常数。

Husch(1982)把q值与负指数分布联系起来，得到

(2)

式中，q相邻径级株数之比；a负指数分布的结构常数；h径级距；e自然对数的底。

显然，如果已知现实异龄林株数按径级的分布，通过对(1)式作回归分析，求出常数k和a，再把a和径级宽度h代入(2)式可求得q。de Liocourt(1899)认为，q值一般在1.2~1.5之间。也有研究认为，q值在1.3~1.7之间(Garcia et al., 1999)。如果异龄林的q值落在这个区间内，认为该异龄林的株数分布是合理的，否则是不合理的。

1.2.2 树种多样性约束

设置树种多样性约束是为满足生物多样性保护的需要。树种多样性约束包括两个方面：树种多样性Shannon-Wiener指数不降低；树种个数不减少，确保不人为造成树种消失。

1.2.3 生态系统进展演替约束

优势种群特别是建群种与森林群落是共存亡的。如果建群种被破坏，那么由它所造成的环境也会改变，不仅增加物种消失的可能性，还会使森林群落类型发生改变(中国森林编委会，1997)。因此，在林分择伐空间优化模型中引入优势种或建群种优势度不降低约束，以维持森林生态系统的进展演替。常用重要值表示一个树种的优势程度，重要值可以用某个种的多度、盖度和频度的平均值或生物量表示(金明仕，1993；张金屯，1995)。本研究用优势种和建群种的蓄积比例表示优势度。

1.2.4 采伐量不超过生长量约束

检查法(control method)是Gurnaud(1886)和Biolley(1920)创立的用于异龄林收获调整的方法，在欧洲应用较广。检查法的核心思想是采伐量不超过生长量。此约束有利于长期维持林地生产力，以便持续获得木材。根据于政中(1993)和Garcia(1999)等修改后的检查法计算生长量公式为

(3)

式中，Z林分定期平均生长量；Δ林分定期总生长量；M₂本次调查全林蓄积；M₁上次调查全林蓄积；C调查间隔期内采伐量；D调查间隔期内枯损量；a调查间隔期。

1.2.5 空间结构约束

采伐后，3个方面的空间结构质量不比伐前差，即林分混交度不低于伐前值；大树和林分聚集指数都不低于伐前值；竞争指数小于或等于伐前值。

1.3 模型建立

在确定目标与设置约束条件基础上，假设已取得样地每木调查因子如胸径、树高、树种和坐标等，建立林分择伐空间优化模型如下。

目标函数：

(4)

式中，M(g)林分混交度，, N样地林木株数(下同)，M_i第i株树木的点混交度，n最近邻木株数，n_i对象木i的n株最近邻木中不同树种个数，CI(g)林分竞争指数(Hegyi，1974)，CI(g)=，CI_i林木i的点竞争指数，L_ij对象木i与竞争木j之间的距离，d_i对象木i的胸径，d_j竞争木j的胸径，n竞争木株数；R(g)林分聚集指数(Donnelly，1978)，, r_i第i株树木到最近邻木的平均距离，F样地面积，P样地周长；Rb(g)大树聚集指数，计算公式同R(g)；σ_M林分混交度标准差；σ_CI林分竞争指数标准差；σ_R林分聚集指数标准差；σ_Rb大树聚集指数标准差；M₀伐前林分混交度；CI₀伐前林分竞争指数；R₀伐前林分聚集指数；Rb₀伐前大树聚集指数；d(g)伐后保留木径级个数；D₀伐前径级个数；q(g)伐后q值；q₁、q₂是q值上、下限；s(g)伐后树种个数；S₀伐前树种个数；t(g)树种多样性指数；p_i树种i的株数百分比；T₀伐前树种多样性指数；f(g)建群种的优势度；F₀建群种的伐前优势度；c(g)采伐量；Z林分生长量；v林木单株材积向量, v=(v₁, v₂, …, v_N)；1=(1, 1, …, 1); g保留木向量，g=(g₁, g₂, …, g_N)，。

目标函数(4)式取最大值使空间结构最优。约束(1)表示林分采伐后，径级个数不减少；约束(2)和(3)表示保留木保持倒J形分布；约束(4)表示树种个数不减少；约束(5)表示树种多样性指数不降低；约束(6)表示建群种的优势度不降低；约束(7)表示采伐量不超过生长量；约束(8)~(11)是初始空间结构约束；保留木向量g的每一个分量取值为0或1。显然，本模型是一个非线性多目标整数规划。

1.4 模型求解

由于模型中存在大量的整数变量，用穷举法难以求解，可采用Monte Carlo法求解此类问题。Monte Carlo法是根据随机抽样的原理，利用计算机语言所提供的随机数函数对约束优化问题的可行点进行快速随机抽样，经过对样本点的目标值过滤比较，找出全体样本点中目标值最优点，将该点视作原问题最优解的一个近似解，即次优解。Monte Carlo法求解的算法步骤如下：

(1) 样地初始保留木向量：g；

(2) 采伐控制参数：起始采伐直径D、采伐强度P、生长量Z；

(3) 计算非空间结构参数：D₀、S₀、F₀、T₀，以及空间结构指数M₀、CI₀、R₀、Rb₀。置初始目标值Q^*=-1000，目标无改善次数U=0，最大允许目标连续无改善次数U₀=500(调整产生多个解)；

(4) 在采伐量不超过生长量的控制下，随机选取采伐木，得到保留木向量g。置连续无可行解次数W=0；

(5) 计算非空间结构参数：d(g), q(g), s(g), t(t), f(g), c(g)；

(6) 如果至少有一个约束条件不成立，则转(4)；否则，转(7)；

(7) 计算空间结构指数：M、CI、R、Rb；

(8) 如果至少有一个空间结构质量降低，W=W+1，判断若W>10 000成立，认为无解，停止计算，若不成立转(4)；否则，找到一个可行解，W=0，转(9)；

(9) 计算可行解的目标值Q；

(10) 如果Q>Q^*，则转(11)；否则，转(12)；

(11) 保留此可行解为当前最优解，g^*=g，Q^*=Q，U=0，转(4)；

(12) 目标未改善次数：U=U+1；

(13) 如果U>U₀，转(14)；否则, 转(4)；

(14) 输出最优解：g^*、Q^*，结束。

根据以上算法可绘出流程图，并编制计算程序。

2 应用实例 2.1 样地调查

吉林省汪清林业局金沟岭林场的原始林20号样地建立于1986年，森林类型为天然云冷杉林(Form Picea jezoensis, Abies nephrolepis)，是该场实施检查法较早的固定样地。样地大小50 m×40 m。为制定下一经理期的采伐计划，根据1986和2002年两次调查结果，采用上述林分择伐空间优化模型进行择伐规划，经理期A=10 a(谢哲根等，1994)。样地用10 m×10 m的相邻网格进行调查，以每个网格为调查单元。每木调查因子包括胸径、树高、年龄、冠幅和坐标等。把外业数据录入计算机，建立数据文件，为规划作准备。

2.2 模型参数

模型共有11个参数，3个林木大小多样性参数，即径级数D₀、株数按径级倒J分布的q值上下限q₁和q₂；2个树种多样性参数树种数S₀和树种多样性指数T₀；建群种优势度参数F₀；生长量参数Z；4个伐前空间结构参数。这些参数均为约束条件的右端项。除计算生长量需要两次调查结果外，其它参数均根据本次调查结果确定。

2.2.1 林木大小多样性参数

本次调查得到株数按径级分布，径级数D₀=7。根据(1)式两边取自然对数线性化，拟合株数对数与胸径的关系曲线，求得参数a=0.0357。径阶宽度h=4 cm。把a和h代入(2)式计算出原始林20号样地的q=1.153 5，与de Liocourt(1899)所确定的q值范围1.2~1.5相比，q值略偏低，说明该林地立地质量较高，可支撑较宽的直径分布范围。据此结果，适当调整q值范围为1.15~1.5，即q₁=1.15，q₂=1.5。

2.2.2 树种多样性参数

树种数S₀=9，树种多样性指数T₀=1.917 5。

2.2.3 优势度参数

占优势的云杉(P. jezoensis)、冷杉(A. nephrolepis)和红松(Pinus koraiensis)3个针叶树种是建群种，优势度F₀=0.521 2。

2.2.4 生长量参数

总生长量根据调查间隔期内林木消长变化计算，平均生长量等于总生长量除以调查间隔期16 a。平均生长量乘以A就是下一经理期的预估生长量Z=19.395 m³，以此作为采伐限额。

2.2.5 伐前空间结构参数

伐前混交度M₀=0.661 9，伐前林分竞争指数CI₀=3.804 6，伐前林分聚集指数R₀=0.925 8，大树聚集指数Rb₀=1.065。

2.3 模型求解

把以上确定的参数代入模型(4)式，用编制的程序求解。起始采伐直径20 cm，控制择伐强度20%(小强度)。用Monte Carlo法求解整数规划，一般得到次优解，但可能求得最优解。当目标函数Q值连续无改善次数u从100~1 500，Q值出现波动。当u大于1 000次后，Q值基本稳定，且在u=1 000时，Q值最大，确定此时的可行解为次优解。

2.4 最优采伐方案

次优解所对应的采伐方案就是本次规划所确定的最优采伐方案(表 1、图 1)。

表 1记录了最优采伐方案中，采伐木编号、树种代码、坐标、胸径、树高、冠幅和材积等因子。共采伐16株林木：2株红松(代码1)、2株云杉(代码2)、5株冷杉(代码5)、6株椴木(Tilia amurensis，代码9)和1株枫桦(Betula costata，代码13)。总采伐材积19.197 0 m³。

图 1显示采伐木在样地中的位置。样地横边为X轴，样地纵边为Y轴。实心圆表示保留木，空心圆表示采伐木，圆的大小代表胸径，标注数字是采伐木号。可见，本规划可精确定位采伐木。

采伐前后林分结构变化见表 2。在非空间结构方面，径阶数和树种数均未减少；树种多样性指数和建群种优势度都增加；倒J形分布的q值在给定合理范围1.15~1.5之内；采伐量未超过生长量，采伐强度仅为18.247 1%。在空间结构方面，树种多样性混交度、大树聚集指数和全林聚集指数都增加。竞争指数减小，且减小幅度较大，达39. 72%，说明采伐能显著降低林分竞争水平。林分整体空间结构的稳定性普遍增强。林分空间结构总目标函数值比伐前提高近4倍。因此，本采伐方案最大限度地改善林分空间结构，又不破坏非空间结构，综合分析表明，此方案可作为制定采伐计划的依据。

3 结论

以空间结构为目标，非空间结构为主要约束，建立林分择伐空间优化模型的方法，突破传统的以功能优化为目标的建模思想局限性。模型集成现代森林经理学、近自然林业、森林生态系统可持续经营、生物多样性保护与计算机技术，并成功地与检查法相结合。

林分空间结构包括多个方面。所以，林分择伐空间优化问题必然是一个非线性多目标整数规划问题。Monte Carlo法是求解的可行方法之一。

用模型进行择伐规划，求解得到研究林分的最优采伐方案，包括16株采伐木。此方案采伐量没有超过生长量，采伐强度仅为18.2471%。伐后空间结构得到较大改善，非空间结构未被破坏。且方案提供了精确的决策信息，是科学经营的依据。

4 讨论

正如于政中(1993)所指出，在一次采伐中一般不能获得完全调整好的由q值所确定的最佳曲线。林分空间结构也必须通过多次调整才能趋于理想状态。调整次数取决于初始林分空间结构与理想空间结构之间的差距。在每次调整中，模型重视经营措施对空间结构的影响，避免产生对森林结构的不良影响。

模型中，大多数约束是针对林分结构多样性与稳定性设置的，更适用于生态公益林。模型没有直接涉及经济效益分析。对兼有重要生态意义的商品林，以经济效益为目标，以空间结构或非空间结构为约束，同样可以建立类似的空间优化经营模型。

对明显不具有培育前途的林木，如病虫害或火干扰后的残留木，以及确无保留价值的过熟木，最好把它们排除在优化之外。因为，我们事先就知道它们无需优化。

林分择伐空间优化模型的建立依赖于林木空间位置信息，如果能从数字化的航空影像获取林木空间信息，再与林分空间结构分析理论相结合，将成为低成本、高效率、大区域森林可持续经营的空间途径。