刘 峰 李建军 谭 畅 杨志高

(1.中南林业科技大学理学院， 长沙 410004； 2.数字洞庭湖南省重点实验室， 长沙 410004；3.中南林业科技大学经济学院， 长沙 410004)

0 引言

林窗(Canopy gap)指森林群落中主林层受人为或自然干扰在林地上形成的不连续林中空隙地，是促进森林植被更新的重要空间[1]。林窗是森林中最普遍的干扰形式，是森林群落演替的主要驱动力，对森林生态系统的恢复和保护具有不可或缺的作用。迄今，国内外学者对林窗开展了广泛的应用研究，主要在林窗测量、林窗更新、林窗环境异质性、林窗动态以及森林经营等方面取得了显著成果[2-5]。相比之下，对于分布在长江中下游流域的亚热带阔叶林的林窗干扰实证研究很少，一些问题尚不明确，例如该区域林窗的结构特征、空间分布模式以及影响林窗发生的主要空间变量等，解决好这些问题对以林窗干扰为主要模式的近自然森林经营管理具有重要意义。

林窗特征主要包括林窗大小、形状、垂直结构以及形成方式[6]。尽管野外人工测量能在小尺度上准确获取相关信息，但将调查结果外推至大尺度甚至景观水平时往往不可靠。已有学者利用遥感技术对林窗特征进行了多尺度监测，但对林窗空间分布模式缺乏进一步研究。在成熟林的森林结构动态变化中，林窗的形成和闭合过程大致平衡，可将林窗假设为与周围冠层有明显区别的离散个体[7-8]，进而可利用空间自相关(例如Ripley’sK、Moran’sI)对林窗中心位置进行多尺度空间分布模式分析。林窗的形成受地形、植被、气候和人为活动等因素(在地理学范畴都能以空间变量的形式存在)的综合影响，探究这些空间变量对林窗发生的影响规律，不仅可以明确林窗干扰的机制，还可为制定适应气候变化的森林经营决策提供科学依据。

本文以亚热带阔叶林为研究对象，利用机载激光雷达和野外调查相结合的方式，收集并处理1块次生林和1块人工林样地的林窗信息，采用描述性统计方法量化不同样地和干扰类型(自然干扰和人为干扰)的林窗特征并检验其差异；利用Ripley’sK函数解析多尺度下不同林窗的空间分布模式；利用空间点模式分析方法解析林窗发生的主要影响因素，以期为亚热带阔叶林林窗特征提取和空间分布模式研究提供方法借鉴和技术支撑。

1 研究区与研究方法

1.1 研究区概况

研究区位于湖南省云山国家森林公园(26°25′～27°00′N、110°22′～113°3′E)，地处祁邵丘陵区，植被类型以常绿阔叶林为主，海拔在550～850 m之间，地形起伏较大，三面环山，南高北低。次生林样地呈现异龄、林层多、间距不规则等成熟林特点，优势树种有青冈栎(Cyclobalanopsisglauca(Thunb.) Oerst. )、檫树(Sassafrastzumu(Hemsl.)Hemsl)、锥栗(Castaneahenryi)和麻栎(Quercusacutissima)等。人工林样地呈现近熟林特点，优势树种有杉木(Cunninghamialanceolata)、马尾松(Pinusmassoniana)、桤木(AlnuscremastogyneBurk.)和木荷(Schimasuperba)等。研究区内林窗类型以自然干扰(如大风、冰雪、病虫害等)为主，主要表现为伏倒、折枝折干以及枯木等形式。人为干扰林窗则主要表现为采伐、放牧和林火等形式。

1.2 试验数据

2016年7月采用机载激光雷达系统(LiteMapper5600型)获取试验区植被点云信息，航带间重叠率约20%，点云密度为3～4个/m2，点云数据采用LAS格式，投影方式为UTM，参考椭球为WGS84。对点云进行预处理并提取林冠点云，生成2 m分辨率的数字冠层模型(Canopy height models，CHMs)，具体步骤见文献[9-10]。利用多阈值组合方法在CHMs上识别林窗并生成专题图[11](图1)，根据试验数据设定CHMs高度阈值大于等于5 m、水平面积大于等于25 m2、CHMs坡度阈值小于等于60°。

收集研究区内森林资源连续清查固定样地数据，根据遥感范围进行同期野外调查，采用鱼眼镜头照片分析法[12]实测林窗大小，DGPS测量单木位置，伸缩式测高器测量树高。选取71株单木和56个林窗的野外测量数据验证CHMs精度，冠层位置、高度和林窗面积的拟合精度R2分别为0.72(p=0.045)、0.71(p=0.05)和0.67(p=0.05)，结果表明利用CHMs数据能有效提取林窗水平和垂直方向的信息。为适应本文点空间模式分析方法的要求，提取林窗多边形的中心点作为林窗空间位置，并剔除面积超大(≥1 000 m2)和形状特殊(线条状)的林窗。

1.3 研究方法

选取林窗面积(AREA)、形状复杂指数(Gap shape complexity index，GSCI)、植被高度指数(Shannon’s height diversity index，HDI)和林窗径高比(Gap diameter to perimeter tree height ratio，DTH)等4个公认指数[11,13]表征林窗特征，4个指数均从CHMs上提取。

AREA为CHMs上提取的林窗多边形面积。形状复杂指数计算公式为

(1)

式中P——林窗多边形周长

A——林窗面积

植被高度指数计算公式为

(2)

其中，Hi是林窗内以2 m×2 m为栅格单元计算林下植被不同高度等级(0.3～1.5 m、1.5～3 m、>3 m)的比例，HDI(0～1)取值越大，表示林窗植被高度异质性越明显，反之亦然。

林窗径高比计算公式为

(3)

式中R——4个方向林窗半径的平均值

T——林窗边缘木高度

DTH取值越大，表示林窗光照增强越明显，反之亦然。

利用t检验法(t-test)分别判断AREA、GSCI、HDI和DTH等指数在不同样地林窗上的差异。利用单因素方差分析(One-way ANOVA)判断自然干扰和人为干扰林窗在不同特征上的差异。采用广义Ripley’sK-function方法研究林窗分布中的非齐次泊松过程(Inhomogeneous Poisson process，IPP)，定义广义Ripley’sK-function[14]为

(4)

在研究区内每一点对象以距离d为半径划定圆，在“完全空间随机”(Completely spatial random, CSR)的零假设条件下，计算落在圆内的点对象期望值与该区域的点对象密度的商，即为K(d)，采用蒙特卡洛模拟并构建95%置信区间包络线，当L(d)值超过包络线的上、下限时，分别表示显著聚集和规则分布模式，当L(d)值在包络线的范围内，表示随机分布模式。

采用空间点模式分析方法[15]解析空间变量(Spatial covariates)对林窗发生密度的影响程度。空间变量分为离散型和连续型，其中离散型空间变量包括林分密度、坡位和坡向；连续型空间变量包括坡度、海拔、与最近道路的距离。根据实测数据和历史记录，以胸径(DBH)大于5 cm为标准划分林分密度等级，其中次生林林分密度划分为高(>1 500株/hm2)、中(500～1 500株/hm2)、低(0～500株/hm2)等级，人工林林分密度划分为高(>3 000株/hm2)、中(1 000～3 000株/hm2)、低(0～1 000株/hm2)等级。以激光点云生成的数字高程模型(Digital elevation model，DEM)为基础[16]，将坡位划分为平坡、坡中、坡谷和山脊；坡向划分为平地、阳坡和阴坡。连续型空间变量直接从DEM计算得出。以上空间变量均以分辨率5 m的栅格数据类型存储，并与林窗专题图进行配准、叠置。

观测值是实际发生的林窗密度，期望值则是在CSR零假设条件下，模拟林窗发生的密度。对于不同离散型空间变量，若林窗发生观测值大于期望值，说明该空间变量或相应区域更容易促使林窗发生。通过卡方检验判断两者之间差异的显著性。

采用潜变量作图法[17-18]揭示连续型空间变量对林窗发生密度的影响，该图通过皮尔逊残差进行分析，若残差接近于零，表明CSR能模拟林窗发生；若残差大于2倍标准差，表明在该连续型空间变量的影响下，林窗发生概率显著高于空间随机条件下的林窗发生概率。

林窗发生在空间分布上具有泊松分布的特点，因此，本文采用IPP模型拟合林窗发生和空间变量的关系，即

λ(u)=exp(θ0+θ1v1+…+θnvn)

(5)

式中λ(u)——林窗分布密度，即单位面积林窗发生的概率

v1,v2,…,vn是影响林窗发生的空间变量，θ0,θ1,…,θn是其对应的参数。利用R软件中的Spatstat程序包进行模型拟合，根据AIC(Akaike information criterion)[18]标准对空间变量的重要性进行筛选，若筛选某一变量前后AIC的变化程度越大，则表示该变量对林窗发生的影响越大。

2 结果与分析

2.1 林窗特征

林窗主要特征如表1所示，林窗平均密度为12.8个/hm2，其中次生林的密度为13.8个/hm2，人工林的密度为11.7个/hm2。林窗平均面积为(78.9±34.2) m2(平均值±标准方差)，其中次生林林窗平均面积(AREA)指数(70.54 m2)小于人工林林窗(102.11 m2,t=2.408,p<0.05)。林窗平均形状复杂指数(GSCI)为2.82±1.27，其中次生林林窗平均GSCI(3.21)大于人工林林窗(2.29,t=1.121,p<0.05)。林窗平均HDI和DTH指数分别为0.58±0.36和0.42±0.05。根据t检验，HDI、DTH指数在次生林和人工林林窗特征上无明显差异。

表1 云山国家森林公园林窗特征统计结果Tab.1 Descriptive statistics for canopy gap characteristics within subtropical broad-leaved forest

根据Jenks自然断裂法，将各个林窗特征指数分成大小两类。就整体林窗而言，面积小于250 m2的林窗个数接近80%，大面积(≥350 m2)的林窗个数不足10%(图2a)。近70%林窗具有较大的GSCI指数(2.75～3.68)(图2b)，具有较大HDI指数(0.56～0.91)的林窗个数与具有较小HDI指数的林窗个数基本相当(图2c)，DTH指数也有类似情况(图2d)。

图2 次生林和人工林样地林窗特征数量统计Fig.2 Percentages of gaps in secondary forest and managed sites within subtropical broad-leaved forest based on four characteristics

根据单因素方差分析，自然干扰林窗和人为干扰林窗AREA、GSCI和HDI指数有显著差异(表2)。具体地，次生林中由自然干扰形成的GSCI、HDI指数分别显著大于人为干扰形成的GSCI指数(p≤0.028)、HDI指数(p≤0.05)，人工林中由人为干扰形成的林窗面积大于自然干扰林窗面积(p≤0.032)，DTH指数在不同样地和干扰类型林窗上的差异不显著。

表2 次生林和人工林样地的林窗平均AREA、GSCI、HDI、DTH指数Tab.2 Mean values for AREA, GSCI, HDI and DTH for gap damage types at three sites

注：*表示单因素方差分析差异显著(p<0.05)。

2.2 林窗分布模式

尽管林窗分布模式有差异，但整体而言，在大尺度上林窗呈现聚集分布模式。其中，次生林林窗在15 m以内呈规则分布模式，在70～90 m尺度呈聚集分布模式，其他尺度呈随机分布模式(图3a)；人工林样地中的林窗呈聚集分布模式，大于60 m尺度时呈显著聚集模式，一些林窗之间的最短距离为10 m左右(图3b)；自然干扰林窗在2块样地大尺度距离时呈聚集分布模式，但均未达到显著程度(图3c、3d)；人为干扰林窗在次生林样地呈随机分布模式(图3e)，但在人工林样地的大于38 m尺度时呈显著聚集分布模式(图3f)。

图3 不同类型林窗在次生林和人工林块样地的分布模式Fig.3 Distributions of canopy gaps of different damage types at three sites tested for deviations from complete spatial randomness

2.3 空间变量对林窗发生密度的影响

在离散空间变量方面，林分密度(高)显著影响林窗在人工林(p<0.05，图 4b)发生的概率，在林分密度大的林地，林窗发生的观测值大于随机分布的零假设期望值，说明该林分更容易受到林窗干扰。坡位(山脊)显著影响林窗在次生林和人工林发生的概率(p<0.05，图4c、4d)，林窗在山脊的发生概率要大于零假设期望值，说明山脊处的单木更容易受到大风损害。虽然阴坡可能引发更多的林窗发生，但总体上坡向对林窗发生密度的影响不显著(图4e、4f)。

在连续空间变量方面，坡度对林窗发生密度有显著影响。在次生林样地中，累积皮尔逊残差曲线在30°坡度附近出现最小值，并超过了两倍标准误差边界，随后有显著拉升(图5a，图中虚线表示2倍标准误差边界，下同)，该结果表明，在缓坡时零假设模型高估了林窗发生的概率，即陡坡(例如大于30°)容易造成林窗，在人工林样地中也有类似情况(图5b)。海拔虽然也对林窗发生密度有影响，但影响不显著(未超过2倍标准误差边界，图5c、5d)，主要原因可能是研究区海拔高度差不大，局部气候变化不明显。与最近道路的距离显著影响人工林样地中的林窗发生密度(图5f)，当距离在30～70 m之间时，累积皮尔逊残差曲线超过边界线，表明零假设模型低估了林窗发生的概率，即林窗更容易发生在离道路较近的区域。主要原因是择伐需要道路运输，其他人为活动也受可达度的影响，所以人为干扰林窗大都集在道路附近。与最近道路的距离对次生林的林窗发生也有影响，但不显著(图5e)。

图5 连续空间变量与零假设模型比较的累积皮尔逊残差曲线Fig.5 Cumulative Pearson residuals plots against three continuous spatial covariates for null model of gap occurrence at three sites

对林窗发生观测数据进行非齐次IPP模型拟合，根据标准判断空间变量(包括离散和连续变量)对林窗发生密度的影响程度(表3)。在次生林中，坡度和坡位是主要影响因素。在人工林中坡度、与最近道路距离和林分密度是主要影响因素。

表3 AIC统计值的非齐次IPP模型拟合参数Tab.3 Coefficients of parameters in inhomogeneous Poisson process (IPP) models fitted according to AIC statistics

3 讨论

本研究发现湖南省亚热带阔叶林自然干扰林窗面积的平均值和变化范围都小于一些温带阔叶林和热带雨林[19]，这可能与成熟单木的间距和自然干扰方式有关。湖南省很少有暴风雪或飓风等强干扰气候，但冬天常出现的冻雨或中等强度的大风可能会导致树冠折枝或折干，造成的林窗面积通常不大。这种林窗在闭合过程中由于边缘木偏冠生长通常快于林下植被的纵向生成，林窗形状呈不规则状态。人工林林窗通常面积较大，理论上应该更能促进林下植被的生长，HDI指数也应较高，但本文结果并不能支持这一假设。主要原因可能是此次机载激光雷达获取的点云密度不高，加之亚热带阔叶林的多林层特性，导致遥感技术对林下幼苗或幼树高度生长状况监测不够，这也是今后进行此类研究应注意的问题。

亚热带阔叶林林窗在大尺度下呈聚集分布模式，这与其他类型森林的林窗分布模式相似[20]。在次生林的中小尺度上，本研究也发现一些小林窗聚集在大林窗附近，这可能是林窗形成木在伏倒过程伤害了周围单木的枝干[21]，造成了新的冠层空隙。伏倒木在人工林中也较常见，特别是在平均胸径较小的林分中呈聚集模式。这与单木的锥状树形有关，例如当冠幅与胸径比例超过一定阈值时，杉木则更容易受外力伏倒。

利用林窗观测数据与非齐次IPP模型拟合，表明影响亚热带阔叶林林窗发生的主要因素包括坡度、坡位、与最近道路的距离和林分密度。在此基础上加入气候因素分析，模拟灾害天气下林窗干扰的趋势，完善亚热带森林经营适应性措施。例如在持续冰冻或暴雪的自然灾害天气下，在坡度大的林分内伏倒的单木要多于坡度小的，因为坡度越大，林分根系发育越不均衡，在外力下越容易倒向一边，这点在杉木人工林中尤为明显[22]。在山脊处林分受到大风干扰容易产生较大面积的林窗，而低的坡位有利于林分抵御大风和冰雪等自然灾害干扰。大的林分密度意味着单木平均胸径较小，由于缺乏大径阶单木的庇护，加之自身健康状况不佳，小径阶林木往往容易受到灾害影响。因此，在实际森林经营过程中宜遵循林窗干扰机制，选择立体条件的林地，避开陡坡、山脊和风口等地形，降低林分密度，促使林木生长健康粗壮，优先营造混交林，促进森林生态系统的恢复和保护。

4 结论

(1)利用机载激光雷达和野外调查相结合的方式能有效提取大范围、多尺度的林窗结构特征。结果表明，林窗平均面积为78.9 m2，林窗平均密度为12.8个/hm2，自然干扰和人为干扰的林窗在面积、形状和林下植被高度等特征上存在差异。

(2)利用Ripley’sK对林窗进行多尺度分布模式分析，自然干扰和人为干扰林窗在中、小尺度上的分布模式有差异，但在大尺度上都呈现出聚集模式。

(3)次生林林窗分布的主要影响因素是坡度和坡位，而人工林林窗分布的主要影响因素是坡度、林分密度和可达度(与最近道路距离)。