机载激光雷达（LiDAR）量化冰暴干扰与抢救性采伐后森林结构动态变化研究

本刊推荐：本研究创新性地运用多时相高密度机载激光扫描（ALS）数据，结合野外实测验证，建立了冰暴干扰后森林结构变化的量化评估体系。研究表明，基于体素（voxel）的生物量损失估算与叶面积密度（LAD）剖面能有效反映冰暴损害格局，而单木分割方法虽低估林木密度减少，但胸高断面积（BA）与蓄积量损失估算与野外测量高度吻合（即使低点云密度下亦表现稳健）。该技术为区域尺度森林干扰监测与管理决策提供了可扩展的解决方案。

　　

摘要

冰暴作为温带森林常见的干扰因子，常导致复杂的树冠部分损毁，形成异质性破坏格局，为生产林管理带来严峻挑战。传统的灾后评估依赖耗时费力的野外调查，难以捕捉损害的空间异质性。遥感技术为此提供了高效替代方案。2014年，一场强冰暴袭击了迪纳拉山脉北部的混交林。本研究利用多时相高密度机载激光扫描（ALS）数据，以野外数据为基准，验证了量化欧洲云杉、银冷杉和欧洲水青冈主导林分冰暴损害的流程。LiDAR衍生的叶面积密度（LAD）剖面与基于体素（voxel）的生物量损失估算有效反映了野外观测模式。基于单木分割的方法虽低估了扰动后林木密度的减少，但胸高断面积（basal area, BA）与蓄积量（volume）损失估算与野外测量高度吻合，即使在低点云密度下亦然。这些方法为损害评估提供了可扩展的途径，增进了对冰暴影响空间变异性的理解，并为拥有区域双时相LiDAR数据集的土地管理者带来广阔应用前景。

引言

中度严重干扰（如冰暴）在多种空间尺度上对森林结构、动态产生持久影响，涉及碳储存、初级生产、结构复杂性、群落组成及生物多样性栖息地等多方面。冰暴因其能导致树木折断、倾倒及树冠部分损毁（从小枝脱落至近全冠丧失）而形成显著复杂的损害格局。其气象条件常覆盖数百至数千平方公里，致使森林损害范围极广。

在木材生产林中，冰暴的大范围与复杂损害模式带来了独特挑战。灾后管理决策（如卫生伐与抢救性采伐）通常基于仓促收集的野外数据，成本高、耗时长且难以准确量化异质性损害。遥感技术，尤其是机载激光扫描（ALS），以其高分辨率、大范围覆盖能力，为量化冰暴后森林结构变化提供了理想解决方案。然而，现有研究很少探讨ALS在量化冰暴损害中的应用，更缺乏对冰暴损害即刻后及抢救性采伐后精细尺度结构变化的评估。

2014年，斯洛文尼亚发生大规模强冰暴，对迪纳拉山脉的温带森林造成中度异质损害。该区域长期实施连续覆盖森林经营，注重可持续木材生产。自然灾害常破坏常规经营计划中的大量木材，因此准确评估扰动后损害格局对于调整长期规划与抢救性采伐操作至关重要。本研究主要目标是评估强冰暴事件对经营混交山林分结构变化的影响，并比较基于野外测量与ALS数据的结构变化推导。研究具体目的包括：评估主要树种林木密度、胸高断面积和立木蓄积量的灾后减少；评估林分水平结构变化；检验灾前灾后叶面积密度（LAD）用于检测冠层结构变化的可行性；评估体素估算方法在量化冰暴及抢救性采伐后地上生物量（aboveground biomass, AGB）损失中的可靠性；评估不同点云密度下LiDAR衍生的林木密度、胸高断面积和蓄积量变化检测的准确性（以野外测量为参考）。

材料与方法

研究区概况

研究在斯洛文尼亚南部迪纳拉山脉的Sne?nik森林经营单元（FMU Sne?nik）进行。该景观以典型高喀斯特地貌为特征，主要土壤类型为Leptosols、Cambisols和Luvisols。气候属湿润温带，年均降雨量1642毫米，年均温6.6°C。冬季（11月至3月）降水占全年42%（690毫米）。森林群落以欧洲水青冈（Fagus sylvatica）、银冷杉（Abies alba）和欧洲云杉（Picea abies）为主，零星分布槭树（Acer pseudoplatanus）、山榆（Ulmus glabra）、欧洲白蜡（Fraxinus excelsior）和椴树（Tilia cordata）。研究森林平均立木蓄积量为442 m³ ha^?1，平均年生长量8.3 m³ ha^?1 yr^?1。过去40年，因银冷杉慢性衰退及青睐异龄林分的经营实践，森林结构与组成发生变化，导致高结构多样性。研究位点（45.672° N, 14.460° E）海拔750–870米，面积约70公顷。

2014年冰暴及其对FMU Sne?nik的影响

2014年1月30日至2月5日，一场大范围强冰暴对斯洛文尼亚超过60万公顷森林造成损害，在克罗地亚也有发生。最严重区域降水超350毫米，平坦表面冰厚达50-90毫米。森林损害在海拔300-1200米最为严重。抢救性采伐于2014年2月冰暴后立即开始，并因树皮甲虫爆发在2015年持续。

研究样地设立与冰暴前后树木测量

2011年，作为Life+ManFor C.BD项目一部分，在Sne?nik试验地设立了六个研究样地。三个样地为稠密（郁闭）冠层，三个为稀疏（开阔）冠层。样地分布于不同树种组成的林分中：两个位于云杉主导林分，两个位于冷杉主导林分，两个位于水青冈主导林分。样地为圆形，面积4000平方米（半径35.7米）。每个样地内，用Leica TS12全站仪测量了所有活树（DBH ≥ 10厘米）的位置（x, y）和胸径（DBH，1.3米高处），并记录树种（共681株）。在南北走向穿过样地中心的10米宽带状区域内，对树木子样本测量了树高（H），用于建立种特异性树高曲线。所有样地在冰暴前四个月（2013年）进行了重测。2014年冰暴后，抢救性采伐完成后，再次重测样地，并记录了采伐树木的 stump 位置。主要对冰暴期间倾倒或折断的树木进行了抢救性采伐，仅树冠受损的树木则被保留。

野外测量树木属性估算

基于每个样地树木子样本的树高数据，建立了特定位点、特定树种的树高曲线，以估算野外未测的树高值。每棵树的胸高断面积（BA）用标准方程计算，而材积（V）和地上生物量（AGB）则利用DBH、H和树种数据，通过R软件中的“rBDAT”库函数 getVolume（从地表至7厘米去皮直径的粗木材带皮材积）和 getBiomass（总地上生物量）分别计算。

机载激光扫描数据获取

ALS数据获取于2013年11月、2014年4月和2014年10月。三次获取均在落叶树种欧洲水青冈无叶期（leaf-off conditions）下进行。研究区用搭载Riegl LMS5600激光扫描仪的Eurocopter EC 120B直升机勘测。该扫描仪相对水平精度10厘米，相对垂直精度3厘米，工作波长1550纳米，最大有效测量距离达1100米（取决于大气条件与目标反射率），可记录多次回波。降尺度前的原始点云密度分别为130.9点 m^?2（2013年11月）、319.4点 m^?2（2014年4月）和164.7点 m^?2（2014年10月）。

点云二次采样与密度变异

后续处理与分析基于预处理点云（即分类并地面归一化）。为比较野外测量与ALS数据评估冰暴及抢救性采伐后林分结构变化，应用了相同分析流程，但使用三种不同点云密度：高、中、低点密度。降尺度使用CloudCompare软件中的Subsample函数，采用空间二次采样方法。该算法通过强制保留点间最小距离来细化原始点云。处理顺序进行，保留第一个点，然后评估每个后续点；仅当无先前保留点位于指定距离内时，才保留该点。输出为所需密度的减少点云。最小点距离设置为：高密度0.1米（平均点密度97.6点 m^?2）、中密度0.8米（平均点密度5.1点 m^?2）、低密度1.2米（平均点密度2.4点 m^?2）。

LiDAR-based森林结构变化分析

本研究采用多种指标比较LiDAR与野外测量的冰暴事件及后续管理后森林变化。这些适用于大区域森林结构表征的指标，涵盖从较粗尺度的冠层结构变化估算（如叶面积密度LAD变化）到更精细尺度的属性（通过体素指标和单木树冠分割捕获生物量、直径、胸高断面积和材积的个体水平变化）。

叶面积密度估算

叶面积密度（LAD）从最高密度点云计算，采用Bouvier等人描述的方法，仅基于首次和单次LiDAR回波。LAD按1米高度分箱计算，以生成三个不同时间点的详细冠层结构垂直剖面。由于LiDAR数据在落叶季收集，欧洲水青冈的LAD值应解释为反映木质元素（枝、干）的结构替代指标，而非实际叶面积密度。未对叶季LAD进行估算校正，该局限在讨论中 acknowledged。为确保准确性，仅使用扫描角±23°内的点（Liu等人推荐）。超出天底角的大角度与树冠 lateral 穿透相关，不适用于代表垂直结构的指标。LAD计算在R软件“lidR”库中实现，基于 gap fraction 剖面，估算到达给定层 z + dz 的激光脉冲数与穿过层 [z, z + dz] 的脉冲数。两者比值的对数除以消光系数 k，得到每个垂直段的LAD估算值。

体素评估地上生物量变化

为评估地上生物量（AGB）变化，使用R软件“lidR”包的 voxel_metrics 函数，以1×1×1米体素分辨率，计算了归一化ALS点云的体素表示。体素指标分别基于2013年11月（冰暴前）和2014年10月（抢救性采伐后）获取的ALS数据集（仅最高密度点云）计算。通过比较2013年11月和2014年10月的体素层，识别出三类体素：a) 两年均存在的体素（vx_common）；b) 仅2013年冰暴前存在、2014年抢救性采伐后缺失的体素（vx_lost）；c) 2014年新出现、2013年缺失的体素（vx_new）。

树木分割与ALS数据树木属性估算

点云数据的单木分割采用Dalponte和Comes提出的方法，分别应用于2013年11月和2014年10月的ALS数据集。分割基于每次ALS数据获取的冠层高度模型（CHM），像素大小0.5米，采用Khosravipour等人开发的 pit-free 算法。单木分割使用“lidR”包的 dalponte2016 函数在R软件中完成。对于ALS数据中每个分割出的树木，树高（H）直接从归一化LiDAR点云导出。冠幅从分割树冠的凸包估算（同样从ALS数据获得）。DBH则使用Jucker等人的异速生长模型估算，该模型将树高和冠幅与DBH关联。基于估算的DBH，计算每棵树的胸高断面积（BA）。材积估算采用2.3.1节描述流程。每个ALS分割树木的树种身份，根据野外数据相应样地的优势树种指定。为研究LiDAR点云密度对基本林分结构特征（林木密度、胸高断面积和立木蓄积量）变化评估的影响，使用了三种不同点密度（见2.4.1节），并分别对每种密度进行了树木分割和属性提取。

统计建模与数据分析

线性混合模型（LMMs）用于评估LiDAR衍生的体素损失对解释野外计算地上生物量变化（ΔAGB）变异的程度。标准化丢失体素体积（vx_lost）作为固定效应，样地（plot）作为随机截距，以考虑六个样地间的变异性。尽管每个样地分为四个子样地，但子样地水平无重复，应用了简化层次模型：ΔAGB_ij = β₀ + β₁vx_lost + b_j + ε_ij，其中ΔAGB_ij为样地j子样地i的预测地上生物量变化，β₀为模型截距，β₁为体素体积损失vx_lost的固定效应，b_j为样地的随机效应，ε_ij为残差。

LMMs也用于评估LiDAR衍生的林木密度、胸高断面积和立木蓄积量与野外测量相比的性能。这使用3种不同点云密度衍生的指标进行，以确定点云密度的重要性。清查方法（野外基础估算和使用高、中、低点密度的ALS衍生估算）作为分类固定效应。随机效应考虑了数据的嵌套结构：Y_ijk = β₀ + β₁method_i + b_k + b_kj + ε_ijk，其中Y_ijk为清查方法i、子样地j、样地k的预测结构指标（林木密度、胸高断面积、立木蓄积量），β₀为模型截距，β₁为清查方法（低、中、高）的固定效应，b_k, b_kj为样地和子样地（嵌套）的随机效应，ε_ijk为残差。

模型在R中使用lme4包实现，使用 lmer 函数配合限制性最大似然估计和Satterthwaite近似求p值。广义混合效应模型的伪R平方使用MuMIn包的 r.squaredGLMM 函数计算。为检查点云密度因子水平间的差异，使用emmeans包的 emmeans 函数计算了估计边际均值（EMMs）。

结果

冰暴及抢救性采伐后树木水平结构变化

总体上，活立木数量从2013年到2014年下降， due to the combined effects of the ice storm and salvage logging。欧洲云杉相对下降最大（50.3%），而银冷杉和欧洲水青冈下降较小（分别为26.6%和20.9%）。2013年，银冷杉的平均胸径（DBH=39.3厘米）和平均材积（V=2.0 m3）最高，欧洲水青冈最低（DBH=30.6厘米，V=1.3 m3）。欧洲云杉的平均树高最高（H=26.0米）。2014年，由于冰暴损害树木的死亡和移除，平均值发生变化。所有三个树种的平均树木胸径和材积均增加。平均树高最大增加发生在欧洲云杉（H=28.5米），而银冷杉和欧洲水青冈种群仅检测到平均树高的微小变化，但两者均显示平均树木材积增加。因冰暴损害移除的树木，其平均胸径和树高低于存活树木。抢救性采伐的树木中，欧洲云杉平均胸径31.2厘米，平均材积1.1 m3，显著低于2014年存活树木的数值。银冷杉和欧洲水青冈记录的值甚至更低。

冰暴及抢救性采伐后林分水平结构变化

所有林分类型中，林木密度、胸高断面积和立木蓄积量从2013年到2014年均下降， due to mortality and salvage logging following ice storm damage。下降最大发生在云杉主导林分，林木密度从2013年371株 ha^?1降至2014年208株 ha^?1（下降44%）。云杉主导林分的胸高断面积和立木蓄积量也分别下降35%和33%。银冷杉和水青冈主导林分的林分参数也有减少，但程度较轻。银冷杉主导林分的林木密度下降24%，胸高断面积和立木蓄积量下降19%。水青冈主导林分显示最小总体下降：林木密度和胸高断面积下降14%，立木蓄积量下降13%。抢救性采伐移除树木数量最低记录于水青冈主导林分（39株 ha^?1），其次是冷杉主导林分（65株 ha^?1）和云杉主导林分（164株 ha^?1）。

冰暴及抢救性采伐前后叶面积密度（LAD）剖面

LAD剖面显示，冰暴前、冰暴后及抢救性采伐后，林分间存在 considerable variation。在一些林分中，LAD明显集中在垂直剖面上部（如水青冈主导稀疏林分），而银冷杉在稠密和稀疏林分中均显示更均匀分布。冷杉主导林分受冰暴损害最小，抢救性采伐砍伐的树木很少。水青冈主导林分遭受 considerable 树冠损害，但大多数情况下，受损树木未被移除。对于云杉主导林分，冰暴直接影响相对较小，但后续抢救性采伐强度大，导致LAD剖面显著下降。

冰暴及抢救性采伐后体素估算地上生物量损失

线性混合效应模型显示，标准化vx_lost与AGB损失间存在强且统计显著的正相关关系。具体而言，vx_lost的固定效应估计值为2645.1（p=0.0022），表明标准化体素损失每增加一单位，地上生物量减少约2645.1千克或每个体素损失1.52千克。模型截距估计为4954.5（p=0.0030）。随机效应显示，森林样地间存在 considerable variability。边际R2（R2m=0.470）表明固定效应单独解释了生物量变化约47%的方差，而条件R2（R2c=0.635）表明完整模型（包括随机效应）解释了总方差约64%。

冰暴及抢救性采伐后林分结构变化：野外基础与LiDAR基础估算比较

林木密度减少估算的比较分析揭示了野外测量与LiDAR基础方法间的显著差异。截距值-9.0（p=0.0023）代表野外测量确定的平均林木密度减少。与野外测量相比，所有点密度方法 systematically 低估了林木密度减少，而LiDAR衍生估算间的差异无统计显著性（p>0.05）。结果还显示不同样地间损害严重性的变异性（plot_group，方差=16.72），表明林木密度减少程度在研究区 vary。此外，子样地水平也观察到变异性（plot_group:subplot，方差=2.99），表明损害即使在单个样地内也非空间均匀。

胸高断面积减少的分析显示，LiDAR基础方法与野外参考测量间无统计显著差异。截距值-0.78（p=0.002）反映野外测量记录的平均胸高断面积减少。无一LiDAR基础方法与野外估算有显著差异，而且不同点云密度的LiDAR衍生估算间也未检测到统计显著差异（p>0.05）。结果还显示不同样地间胸高断面积减少严重性的变异性（plot_group，方差=0.10），以及样地内子样地水平的变异性（plot_group:subplot，方差=0.10）。

立木蓄积量减少的分析显示，LiDAR基础方法与野外参考测量间统计不显著差异。截距值-10.27（p=0.007）反映野外测量记录的平均材积减少。结果还显示样地间材积减少严重性的变异性（plot_group，方差=22.22），以及样地内变异性（plot_group:subplot，方差=21.92）。此外，事后成对比较显示，三种点云密度间的立木蓄积量估算无统计显著差异（所有p>0.05）。这证实了即使从较低密度点云也能可靠估算蓄积量变化。

讨论

机载激光扫描在捕捉冰暴及灾后抢救性采伐引起的森林结构变化方面表现良好，但性能因特定LiDAR衍生指标而异，凸显了LiDAR方法的效用与局限。此外，对于某些指标（如胸高断面积和蓄积量变化），低分辨率LiDAR数据与高分辨率数据表现相当，这是一个鼓舞人心的结果，因为低分辨率数据更可能在大范围森林区域可用，且需要更少的处理和数据存储计算资源。

主要树种的灾后林木密度、胸高断面积和立木蓄积量减少，与2014年冰暴事件损害格局的其他独立评估一致。 notably 水青冈及阔叶树种整体最易受严重冰损害，而冷杉和云杉最不易感。Nagel等人也发现冠层损害在水青冈中常见，尤其较小尺寸树木和中低风暴强度下。他们还报道较大树木更易倾倒，这与研究位点野外观测一致，一些较大云杉和冷杉树木在冰暴期间倾倒或基部折断。这些损害格局反映在2014年4月叶面积密度（LAD）垂直剖面中，水青冈林分倾向于显示冠层高度最高水平冠层损失，破碎冠层材料在冰暴后移至近地面较低高度。然而，重要的是注意LiDAR数据在落叶季获取，这意味着报告的LAD值主要反映木质结构元素（如枝、干），不能直接与叶季数据绝对值比较。尽管叶季与落叶季植物面积密度剖面间绝对差异可能很大，但相对变化（如扰动前后）仍然有效且信息丰富。

灾后抢救性采伐也与斯洛文尼亚及中欧自然灾害后管理建议一致。倾倒或折断树木，无论树种，通常因经济原因被移除，而云杉受损活立木则常规移除以限制后续树皮甲虫爆发。抢救性采伐后结构被LAD剖面方法相对较好地捕捉。云杉显示冠层较低值，反映大量云杉从林分中移除，而水青冈主导林分显示冰暴后与抢救性采伐后LAD剖面间微小差异。冷杉主导稠密林分灾后干预 minimal，但稀疏银冷杉林分的LAD剖面变化更类似于云杉主导林分，可归因于较多大型冷杉树木在冰暴期间倾倒、基部折断或失去顶梢。

体素方法在捕捉生物量损失方面也表现良好，证实体素指标是地上生物量（AGB）变化的稳健预测因子。然而，所选体素大小（即1 m3）用于记录体素与AGB损失间关系的文档，可能过于粗糙无法捕捉冰暴造成的部分损害，这结合点云密度，可能是效应估算高变异性的原因之一。体素大小是点云体素化的关键参数，因为它强烈影响形状识别、表面识别和计算效率。过细体素大小可能导致数据冗余和高处理开销，而过粗大小可能导致信息丢失并将不同结构合并 within a single voxel。因此，体素大小选择取决于研究目标、对象复杂性和LiDAR平台。为本研究，对三种体素大小（即0.5米、1米和2米）进行了敏感性分析，显示跨三种分辨率的一致结果。为清晰和与其他研究可比，所有分析使用了1 m³体素大小。总体，这些结果 promising 并支持测试方法的适用性，但更细分辨率在某些情况下可能 yield improved outcomes。

比较基于单木分割的LiDAR衍生林木密度、胸高断面积和立木蓄积量变化与野外测量的结果好坏参半。LiDAR在所有点云密度下低估了林木密度减少。估算林木密度可能受许多误差源影响， particularly 分割过程。自上而下算法，如基于分水岭的，已知在针叶林内表现更好 due to their conical architecture。与Li等人提出的方法相比，分水岭算法遭受严格2D方法 related to the use of a rasterized canopy model to segment the point-cloud below，且缺乏识别共优和下层木的能力。在受损林分中识别更具挑战性，其中树木不再有明确定义冠层， whereby 受损树木可能失去种子点（即树梢），这可能导致不同分割输出。

与密度估算相反，LiDAR衍生的胸高断面积和材积损失估算与野外测量高度匹配 across all point densities。林木密度估算通常相比其他标准森林清查参数具有最低精度。胸高断面积和材积等属性对与单木识别相关的问题较不敏感， likely because the small trees in the understory that were difficult to identify do not have a large influence on these estimates。鉴于胸高断面积和材积损失估算也适用于较低分辨率，它们为土地管理者带来希望，他们可能拥有国家或区域尺度双时相LiDAR数据集，这些正变得更易获取但通常具有较低点密度。

尽管允许直接比较我们结果的研究很少，机载激光扫描已显示通过降低冠层高度、垂直结构变化和冠隙形成等指标有效检测雪诱导森林扰动。这些变化可使用双时相或单时相ALS数据集量化，并通过冠层高度模型（CHMs）或百分位基础高度分布衍生的指标表达。鉴于冰暴引起与重湿雪类似的机械性树冠损害，这些ALS衍生指标可能可转移用于检测冰相关森林损害。

尽管LiDAR显示检测冰暴损害及后续抢救性采伐的前景，本研究有几个值得关注的局限。最 notably，我们的案例研究基于相对较小采样区（即70公顷研究区内六个0.4公顷样地）和异质森林结构，因此推广结果 beyond our study area 应谨慎。还有一些方法学挑战需要进一步考虑。用于估算叶面积密度（LAD）的点云密度可能不足以完全解析冠层的垂直结构， especially in complex, multi-layered, dense forest stands and in leaf-off conditions。类似地，冠层变化的体素分析依赖于结构细节水平，难以用中密度点云实现，这可能损害地上生物量损失衍生估算的准确性。由于本研究重点不是测试不同单木分割算法，我们未评估单木分割准确性，因此无法判断ALS衍生树木水平结构指标反映实际条件的程度， especially in damaged