基于Landsat土地覆盖数据改进国家森林资源清查中乔木覆盖度的模型辅助估计研究

本文推荐一项利用Landsat卫星衍生的土地覆盖数据作为辅助信息，显著提升加拿大国家森林资源清查（NFI）中植被乔木（VT）覆盖度估计精度的研究。研究评估了基于Beta回归的模型辅助估计量（MAbeta）的统计特性，证明其在样本量充足时（n ≥ 100）能提供近乎无偏且覆盖率达95%的估计。与传统的设计基础比率估计量（DB）相比，MAbeta在时间匹配的辅助数据下效率提升显著（相对效率RE达3.25）。此外，研究通过模型引导的年度估计揭示了林分替代干扰是驱动生态区VT覆盖动态的关键因素。结果表明，整合遥感衍生的全覆盖数据产品能有效增强NFI的准确性与精确性，为多尺度森林资源监测提供有力支持。

　　

1. 引言

森林资源的可持续性管理需依赖多时空尺度下森林属性（如面积和蓄积量）的状态与变化监测。加拿大国家森林资源清查（NFI）项目旨在提供生态区及国家层面的森林生态系统状态与趋势信息。NFI依托永久性地面样点（地面样地）和遥感调查样地（照片样地）网络，其中照片样地源自航空摄影或空间分辨率优于1米×1米的卫星数据，作为国家报告的主要数据源。这些2公里×2公里的照片样地系统布设在20公里×20公里的网格上，由于加拿大森林广袤偏远，建立足够大且具代表性的地面样地网络成本过高。

NFI采用十年周期的复查制度，每个复查周期内每年测量的照片样地数量在时间和空间上变异显著。所有在复查周期内测量的照片样地被合并，采用时间无差异的设计基础（DB）估计量来获取给定复查周期内的森林资源状态信息。属性年度或周期性变化的估计则通过连续复查周期中测量的照片样地组合获得。时间无差异估计量假设样本单元间观测时间差异（由于源影像年代差异）及连续测量周期间样本单元的时间差异（或滞后）效应可忽略。然而，对于森林这类动态系统，这种忽略效应的假设可能不成立，时间差异会影响点估计和变化估计的偏差与精度。

与国家周期性抽样调查项目不同，由加拿大林务局运营的国家陆地生态系统监测系统（NTEMS）利用年度卫星数据收集，提供全国森林生态系统的空间全覆盖信息。现有多种公开的土地覆盖和森林结构相关数据产品。这种多时序、全覆盖信息可与样本数据（NFI照片样地数据）结合，采用更高效的统计方法以提高清查估计的精度。

模型辅助（MA）回归估计量利用辅助数据提高效率，同时保持近似设计无偏的特性，无论模型是否有效。这些估计量即使模型不正确也是渐近无偏的，但当模型未能准确表示调查变量与辅助变量间关系时会效率低下。影响估计量效率的因素包括辅助模型的选择、辅助数据与调查变量的相关程度、感兴趣区域内样本单元数量、数据源间位置匹配、调查与辅助数据间的时间差异以及辅助信息中的误差（如制图误差）等。此外，MA估计量易受小样本偏差影响。模型基础估计量则提供了另一种途径，其假设存在超总体模型，观测值是该模型的实现，但这需要准确的模型设定，否则可能引入模型偏差。

利用调查与辅助数据间存在时间差异的森林清查估计研究建议，仅使用与辅助数据同年的样本样地进行有效估计。当辅助数据与NFI调查间收集时间存在显著差异（最长可达10年）时，警告MA估计量可能因时间间隔内未检测到的森林损失而比简单随机抽样估计量效率更低。因此，辅助数据与调查数据间较长的时间间隔会引入人为噪声，从而增加辅助模型的残差方差。在森林清查中，噪声程度可能随时间间隔长短及区域内干扰动态导致的响应（或调查）变量变化程度而异。若辅助信息在复查周期内每年可用，则可评估调查与辅助数据时间不匹配对清查估计的影响。如此，借助遥感提供的年度、空间全覆盖、分类兼容的辅助变量，可开发更细致的方法将辅助数据（NTEMS数据）与NFI数据结合，以提高NFI估计的精度。

在此背景下，本研究的主要目标是开发并测试统计方法，将这些卫星衍生数据产品与NFI照片样地数据结合，以提高2007至2017年复查周期内大西洋海事生态区植被乔木（VT）覆盖度比例估计的精度。尽管Beta回归已被建议用于建模连续比例，但我们未发现在同行评审文献中评估以Beta回归为辅助模型的MA估计量统计特性的研究。逻辑广义回归估计量已被引入用于总体和子总体水平的类别频率、二元指标和计数型调查变量。使用逻辑广义回归估计量针对伯努利和二项响应变量的研究得出结论，传统标准方差估计量对小样本量会产生轻微低估。在林业中，多项研究使用MA估计量，以逻辑回归为辅助模型估计森林覆盖比例。

基于这些研究，我们预期以Beta回归为辅助模型的模型辅助回归估计量（MA_beta）与逻辑广义回归估计量表现相似。然而，确认这一预期至关重要。因此，我们为本研究制定了三个具体目标。第一个目标是利用人工模拟的总体数据评估以Beta回归为辅助模型的MA估计量的统计特性。第二个目标是调查调查与辅助数据间时间差异对MA估计量所用辅助模型性能的影响。当使用周期性清查进行点估计时，假设感兴趣变量在区域内无重要的空间和时间趋势。因此，第三个目标是检验这一假设的有效性。为此，我们旨在生成基于模型的生态区水平VT覆盖度比例年度估计，以评估趋势（增加、稳定或减少），这是监测森林资源可持续性的重要指标。此外，若周期内存在显著趋势，周期性估计可能存在偏差且可靠性较低。评估VT覆盖度比例趋势对干扰显著性的敏感性有助于识别该地区森林生态系统动态的关键驱动因素。

2. 材料与方法

2.1. 材料

2.1.1. 研究区

选择位于加拿大东南部的大西洋海事生态区作为试验点，其包含从大西洋海岸、阿巴拉契亚山脉到新不伦瑞克高地的多样森林生态系统。该生态区陆地面积为201,570.85平方公里。虽然新不伦瑞克和布雷顿角高地可见北方森林和泰加林，但温带阔叶林（主要为枫树和桦树）和混交林（落叶和针叶）是该生态区其他部分的主要森林类型。该生态区包括四个NFI分析单元：魁北克省东南部（单元10；QC）、新不伦瑞克省（单元6；NB）、新斯科舍省（单元4；NS）和爱德华王子岛（单元5；PE）。采伐是该生态区主要的林分替代干扰。

2.1.2. 国家森林资源清查设计与数据

加拿大NFI设计为单一系统位置（采样点）样本，位于20公里×20公里固定国家采样网格的节点上。每个节点中心是一个2公里×2公里的照片样地。本研究所用区域样本包括位于2016万公顷大西洋海事生态区内的NFI照片样地完整列表（n = 542，采样强度约1%）。1%的系统采样强度足以在国家和生态区层面生成关键森林属性的有效统计估计。按照标准森林清查实践，由经过认证的技术人员手动判读代表相对均质森林条件林分的多边形，并获取NFI所需属性。最初，NFI的相片判读基于常规无云立体航空照片，比例尺范围从1:10000到1:20000。照片被数字化至至少1米分辨率，用于软拷贝摄影测量系统。随着数码相机的出现，地面采样距离为30厘米的航空影像已普遍用于森林清查应用中的相片判读。在航空摄影不可用或采集逻辑困难的地区，如加拿大北部，使用高空间分辨率卫星影像（即≤1米空间分辨率）填补采样空白，并遵循相同的判读程序。NFI第一个复查周期中使用的卫星传感器示例包括WorldView和QuickBird。判读后，基于多边形的信息被聚合为表格结果，代表NFI采样网格上的给定节点。除97个从基线（2000-2006年）更新的照片样地外，所有照片样地的测量均在一个测量周期内完成，测量周期（2007至2017年）内每年测量的照片样地数量变异很大。在当前NFI估计过程中，位于NFI分析单元（省/地区边界与生态区边界的地理交集）内的所有照片样地数据用于估计，且忽略照片样地数据间的时间差异。对于NFI照片样地测量，在航空照片中 delineate 具有多边形的均匀区域，并按层次分类为植被或非植被，然后按土地覆盖类型（有林地或无林地）分类。若多边形至少10%的面积（按冠层覆盖计）由任何大小的树种组成，则视为有林地。本研究中，VT覆盖比例用作感兴趣变量（调查变量，即P_ij），每个照片样地的VT比例计算为VT类多边形总面积除以照片样地总面积。

2.1.3. 辅助数据计算

使用公开可获取的、干扰信息化的年度土地覆盖图作为本研究辅助信息的主要来源，其土地覆盖类别（非植被、植被无林地、VT）的总体分类精度为90.2%（±1%），源自Landsat时间序列数据。整个研究区被划分为N个2公里×2公里的单元作为总体单元，建立系统网格。在N个总体单元中，NFI照片样地被视为n个样本单元。将系统网格和NFI照片样地数据的空间图层重投影至NAD1983加拿大兰伯特等角圆锥投影和北美1983基准，以匹配年度土地覆盖图栅格图层的投影。然后，通过将分类为VT的像素数（空间分辨率30米）除以网格单元或照片样地空间范围内的总像素数来计算聚合VT覆盖比例。为匹配调查变量与辅助变量的时间戳，使用2006年至2017年的年度土地覆盖图推导所有采样单元（即照片样地）的聚合VT覆盖比例。测试了两种不同的辅助信息聚合过程：（1）考虑范围内的所有像素；（2）应用排除边界像素的缓冲区，使用基于所有像素的聚合辅助信息的模型给出较低的均方根误差（结果未显示）。因此，在聚合辅助信息时考虑了照片样地或网格单元范围内的所有像素。基于Landsat时间序列数据衍生的年度土地覆盖图的聚合VT覆盖比例用作给定年份的x_ij（即第j个NFI单元中第i个照片样地的辅助变量值）。

如表1所示，在第一个复查周期的542个NFI照片样地数据中，97个照片样地（约占总照片样地的18%）从NFI基线（2000至2006年）更新。因此，不能保证这97个照片样地独立于本研究所用的辅助信息，因为2000至2006年的NFI照片样地数据是土地覆盖分类过程训练数据集的一部分。然而，保留这些照片样地在分析中以防止潜在的覆盖偏差。

使用Landsat衍生的森林采伐和森林野火干扰（均为年度）数据估计研究区的干扰显著性。给定年份的干扰显著性定义为该年估计的总干扰占生态区总陆地面积的百分比。总干扰结合了火灾烧毁区域13年累积和和采伐区域6年累积和，假设林分替代干扰后乔木覆盖（至少10%冠层闭合度）的恢复平均需要4.5年（采伐，标准差=3.43年）和8.5年（野火，标准差=5.92年）。为保守估计恢复到10%冠层闭合度（符合NFI对VT类的定义），我们使用光谱恢复的上四分位数，即大西洋海事生态区采伐干扰为6年，火灾干扰为13年。

2.2. 统计方法

2.2.1. 使用模拟数据评估MA估计量

Beta回归模型在响应变量为连续比例时表现良好。使用人工模拟的总体数据评估以Beta回归为辅助模型的MA回归估计量（称为MA_beta）的统计特性。创建了两个不同的有限总体，各包含N = 100,000个单元。第一个基于随机变量（y_i）的对称分布（累积分布函数遵循逻辑分布），而第二个使用变量的非对称分布（累积分布函数遵循Gumbel分布）。从每个总体中随机无放回抽取七种不同样本量（n = 15, 30, 50, 100, 500, 2000, 和5000）。该过程重复5000次（Γ = 1,…,5000），每次迭代中从每个样本生成MA估计。MA估计量可表达为特定公式。

对于β参数估计，假设y ～ Beta（μ, φ），0 < μ < 1，且φ > 0，其中E（y）= μ，且方差有特定形式，φ是精度参数。然后，响应变量y_i被建模，其中β是未知参数列向量，x_i是自变量行向量，η_i是线性预测因子。

回归估计量的方差可估计。其中g_i是g-权重，估计方式确保样本在校准已知总体辅助变量均值时，g-权重的样本均值等于1，且辅助变量x的样本均值乘积等于辅助变量的总体均值。

从模拟研究的实现中评估了均值相对偏差比例、标准误差估计相对偏差比例和经验覆盖度。经验覆盖度指的是估计量95%置信区间的实际覆盖度。计算置信区间时，假定样本量n的样本均值分布服从中心极限定理下的正态分布。传统标准方差估计量（即取g_i=1 ? i）已知对小样本量会产生轻微低估。因此，分别评估了在方差估计量中使用g-权重（来自公式）和不使用g-权重（g_i=1 ? i）时的覆盖率（CR），以检验g-权重是否减轻小样本量对95%置信区间CR的影响。

2.2.2. 国家森林资源清查数据分析

2.2.2.1. 直接设计基础估计

遵循当前NFI估计程序获取直接DB估计。在估计过程中，NFI分析单元被视为分层。该程序使用比率估计量估计每个NFI单元中VT覆盖的平均面积比例及相关方差。在系统抽样中无设计无偏方差估计量，因此最常使用简单随机抽样方差估计量作为“保守估计”。NFI单元“j”中VT面积比例平均值的近似方差估计量由特定公式给出。

通过将每个NFI单元视为独立分层，使用分层随机抽样公式获得生态区水平平均值。对于大西洋海事生态区，均值的方差估计为特定形式。

2.2.2.2. 使用经验数据的VT覆盖比例模型辅助估计

提出的MA估计量用于估计给定NFI单元j的VT覆盖面积比例，构建方式考虑了部分照片样地跨越NFI单元或生态区边界的影响，通过在估计量中纳入基于照片样地面积（M_ij）的权重。为评估样本依赖的g-权重是否会改善结果，给出了包含g-权重的近似方差估计量。g_ij是第j个NFI单元中第i个样本单元的g-权重，并使用公式估计。为获得测量周期的生态区水平估计，使用公式组合NFI单元水平估计值用于均值和方差估计。为评估估计量的可靠性，估计了相对标准误（RSE）。通过比较替代估计量与直接DB估计量的效率，估计了相对效率（RE），其中特定方差分别为使用DB直接估计量和MA估计量估计的方差。

Beta回归用于估计公式中的参数向量（β）。估计过程及与Beta回归相关的假设详见第2.2.1节。在初步试验中，我们比较了三种不同的链接函数（logit、probit和互补双对数），发现互补双对数链接函数更适合我们的数据。使用R包“betareg”拟合模型。β参数的估计值见附录B。

在存在分层的情况下，回归估计量可以通过为每个分层单独拟合回归或合并所有分层拟合回归来形成。我们将所有NFI单元在第一个复查周期（2007年至2017年）内测量的照片样地数据合并拟合回归，以从NFI单元外的照片样地借用强度。更具体地说，我们通过组合来自相关区域和时间段的样本单元来借用强度，以使用大且空间代表性的样本估计回归系数（β参数），使回归模型提供更准确可靠的预测。此外，当我们将拟合到生态区所有单元的模型应用于单个NFI分析单元时，辅助模型变为准外部，且MA回归估计量具有更稳健的统计特性。航空照片基础的VT覆盖比例用作响应变量，Landsat衍生的VT覆盖比例以及分类变量NFI单元ID（Prov）用作解释变量。协变量Prov（分类变量）用于模型中，为每个NFI单元提供不同的截距。分别使用2017年（期末）辅助数据和使用时间戳匹配的辅助数据（即辅助数据日期与航空照片获取日期相同）拟合模型，以评估调查与辅助数据间时间差异对MA回归估计量所用辅助模型性能的影响。

2.2.2.3. 大西洋海事生态区乔木覆盖动态

使用2007年至2017年测量的所有照片样地，可以为整个复查周期生成VT覆盖比例的MA估计。然而，对于VT覆盖比例的年度估计，由于样本量较小且年度面板空间覆盖度差，无法维持基于MA方法的推断有效性，因此我们依赖基于模型的推断。年度测量面板在样本量和空间覆盖上变异很大。然而，在基于模型的方法中，总体估计（给定年份VT覆盖比例的平均估计）被视为随机变量，并通过使用给定年份每个总体单元的模型预测生成。在这种情况下，推断的有效性基于正确的模型设定而非概率样本。因此，使用辅助模型生成基于模型的VT覆盖比例年度估计，使用卫星衍生的年度土地覆盖数据作为大西洋海事生态区的输入。对于非线性回归模型，基于均方误差估计量的模型预测不确定性估计可能存在显著偏差。另一方面，研究表明自助法（Bootstrap）针对模型设定错误提供了保障。因此，我们使用“配对自助法”计算自助均值（MB_BAE）及相关标准误（SE）。我们通过重复重采样原始样本（有放回抽样）选择了1000个自助样本（重复）。对于每个自助样本，我们拟合了一个Beta回归模型，估计了模型参数，使用该模型预测总体参数（即平均VT比例），然后使用公式分别估计了每年的总体参数（均值和SE）。

3. 结果

3.1. 基于模拟研究的MA_beta统计特性

分析表明，MA_beta估计量近乎无偏。在对称和非对称总体中，所有样本量的相对均值偏差均可忽略不计。然而，估计相对误差的分布对于小样本量较大，并随着样本量增加逐渐减小。发现非对称分布总体中估计相对误差的分布更大。

SE估计量显示，从对称和非对称分布总体中抽取的中等和大样本（n ≥ 100）具有较小的相对偏差（-0.05 < SE < 0.05）。这证实了使用MA_beta方差估计量推导的SE估计量是渐近近似无偏的。发现小样本（n < 100）存在负的相对偏差，表明MA_beta估计量低估了小样本的方差。负的相对偏差随着样本量减小而迅速增加。在方差估计量中纳入g-权重有助于略微减轻小样本效应。

对于中等至大样本，经验CR发现接近0.95名义覆盖度，适用于两个总体。然而，小样本量未达到0.95名义CR。使用提出的g-权重仅能略微缓解此问题。带有g-权重的方差估计量始终为小至大样本量的两个总体产生更好的CR（CR > 90%）。

3.2. MA_beta改进2007–2017年期间VT覆盖比例估计

不同NFI单元中调查变量与辅助变量间的关系表明，总体上调查变量与辅助变量正相关。然而，相关程度在NFI单元间各异，在QC较强，在NS和PE较弱。

发现当调查变量与辅助变量间时间戳匹配时，辅助模型的预测能力（即较小的均方根预测误差）对所有NFI单元均一致改善。

在生态区水平，发现MA_beta估计量对VT覆盖比例估计的效率比直接DB估计量高2.75–3.25倍。当使用时间戳匹配数据估计辅助模型参数（β′s）时，估计量的效率从3.18提高到3.25。此外，在方差估计量中纳入g-权重并未提高估计量的效率。然而，仍使用g-权重以校正与小样本量方差估计相关的潜在偏差。使用g-权重后，当使用时间不匹配数据进行估计时，估计量的效率从3.18降至2.75。类似地，当使用时间匹配数据进行估计时，估计量的效率从3.25降至3.04。

NFI单元水平周期性VT覆盖比例估计及相关抽样误差（方差）如表4所示。使用带有g-权重的方差估计量获取NFI单元水平方差估计，因为模拟研究表明使用g-权重的方差估计量提高了CR。估计量的RE从1.9到8.8不等，在NFI单元5（PE）最高，在NFI单元4（NS）最低。出乎意料的是，NFI单元5的效率增益相对较高（即8.8）。这需要进一步审查以确认这是真实的效率增益，还是与MA_beta估计量相关的人为现象，因为我们的模拟研究表明MA_beta估计量通常低估小样本量的方差。基于模拟研究结果，可以假设其余NFI单元的方差估计近似无偏或具有可忽略的偏差，因为它们有n > 100个样本单元。

3.3. 大西洋海事生态区VT覆盖动态

依赖基于模型的自助估计量，我们生成了每年的VT覆盖比例估计。基于模型的生态区VT覆盖比例年度估计与DB周期性估计比较良好。尽管基于模型的年度估计与DB周期性（时间无差异）估计在平均VT覆盖比例上无显著差异，但基于模型的年度估计更精确且随时间变化不同。总体上观察到上升趋势，表明生态区VT覆盖比例增加。此外，VT覆盖比例轨迹中的非线性时间模式显示，2007年至2015年间VT覆盖比例逐渐增加，随后随着干扰显著性开始增加而下降至测量周期结束。发现估计的VT覆盖动态轨迹对大西洋海事生态区干扰显著性高度敏感，存在强负相关（Pearson r = -0.94，p值 < 0.001）。

4. 讨论

我们提出了一个基于Beta回归模型（MA_beta）的MA估计量，用于计算连续数据的比例调查变量（例如，VT覆盖面积除以VT覆盖面积与非VT覆盖面积之和），以提高周期性估计的精度。通过模拟评估了MA_beta估计量的统计特性，并将此估计量用于经验数据（NFI照片样地数据）。我们发现，与当前使用的DB估计量相比，MA_beta估计量显著提高了2007年至2017年期间VT覆盖比例估计的精度。当辅助模型使用时间匹配的调查和辅助数据拟合时，估计量的效率进一步提高，证明了多时序辅助数据的有用性，特别是在干扰强度和模式非常多样的地区。回归模型还使我们能够使用可作为解释变量的每个总体单元的辅助信息（多时序Landsat衍生土地覆盖数据）生成基于模型的VT覆盖比例年度估计。估计的年度VT比例与生态区期间干扰显著性估计呈强负相关，表明年度VT覆盖比例变化主要由采伐驱动，因为采伐占生态区总林分替代干扰的99%以上。后续章节详细讨论了提出的MA_beta估计量在使用加拿大NFI数据及其他可用辅助信息生成改进的森林属性估计方面的有用性。

4.1. MA_beta估计量的统计特性

本研究的第一个目标是使用人工模拟的总体数据评估以Beta回归为辅助模型的MA估计量的统计特性。文献中有一些使用非线性模型的MA估计量示例。值得注意的是，McRoberts等人建议在MA估计量中使用不同非线性模型形式时谨慎行事，因为非线性回归的MA估计量的统计特性相对较少文献记载。此外，据报道，使用逻辑回归的MA估计量低估小样本量的方差。因此，我们在将MA_beta用于经验数据之前，使用模拟人工总体数据评估了其统计特性。结果表明，MA_beta均值估计量是无偏的，而方差估计量是渐近近似无偏的，对于样本量≥100的样本，真实均值的置信区间覆盖度约为0.95，这与Kennel和Valliant以及McRoberts等人的研究一致。如Myrskyl?所述，分析表明，对于较小样本量（n < 100），方差估计量低估了真实方差。然而，对于中等（n ≥