Comparison and evaluation of several methods for estimating the average density of total forest volume in forest farm
-
摘要:目的 以林场或县森林资源总体为调查对象，及时、准确地调查监测总体平均每公顷蓄积量，对上级（如市、省）部门开展森林资源宏观管理、生态保护价值评价、森林碳储量计量、领导干部任期绩效考核等工作都有重要支撑作用。将卫星、无人机等多源遥感数据作为辅助数据，采用较少抽样调查样地数据，实现总体参数有效估测的新方法，已成为国内外重要的研究方向，但目前，国内尚无多种现有估计方法的比较评价研究。为了促进新一代遥感技术在森林资源调查业务中的应用，提高我国森林资源天空地一体化调查监测技术水平，亟需对现有林场或县总体参数估测方法进行比较评价研究。方法 以内蒙古旺业甸实验林场主要人工林树种为总体，基于2019年在林场获取的无人机激光雷达（LiDAR）抽样数据、Sentinel-2A多光谱数据（全覆盖）和少量样地数据，针对样地（p）、样地−卫星（ps）、样地−抽样无人机LiDAR（pl）以及样地−抽样无人机LiDAR-卫星（pls）4种模式，利用适合这4种模式的概率抽样法（DB）、模型辅助法（MA）、模型法（MD）和混合法（HY）4类共5种估测方法（DBp、MDps、MAps、HYpl以及MDpls）对总体森林蓄积量密度均值（MSVD）进行估计与对比分析。结果 （1）DBp、MDps、MAps、HYpl、MDpls 5种方法估测的MSVD分别为212.54、202.09、202.38、210.07以及208.96 m3/hm2，精度分别为90.44%、91.54%、91.69%、96.35%和96.45%，方差依次减小。（2）其他4种方法相对于MDpls方法的估计效率（RE）均大于1（REDBp,MDpls = 5.39，REMDps,MDpls = 3.82，REMAps,MDpls = 3.69，REHYpl,MDpls = 1.07）；HYpl相对于MDpls的RE略大于1，但几倍于其他3种方法（REDBp,HYpl = 5.02，REMAps,HYpl = 3.43，REMDps,HYpl = 3.56）。（3）包含LiDAR数据的HYpl与MDpls方法相对于包含Sentinel-2A数据的MDps与MAps方法均为正效率（REMAps,HYpl = 3.43，REMDps,HYpl = 3.56，REMDps,MDpls = 3.82，REMAps,MDpls = 3.69）；MDps与MAps方法之间的RE接近1，但MAps的效率微高于MDps（REMDps,MAps = 1.04）。结论 和只利用样地数据的估计方法相比，将遥感数据作为辅助变量建立估测模型，无论是采用对蓄积量不够敏感的林场全覆盖Sentinel-2A多光谱遥感数据，还是采用对蓄积量很敏感的抽样式获取的LiDAR数据，都可有效提高林场总体MSVD的估测精度。涉及遥感数据应用的4种方法，估计精度最高的为MDpls，其次为HYpl，这2种方法都包含了LiDAR遥感抽样观测数据的应用，都是适应于林场总体MSVD估计的年度监测方法。可实现天空地3种观测数据协同应用的MDpls估测精度和相对效率最高，可作为林场森林蓄积量年度监测的首选方法。Abstract:Objective Taking the overall forest resources of forest farms or counties as the object of investigation, timely and accurately investigating and monitoring of the mean stock volume density（MSVD）will play an important supporting role in the macro management of forest resources, the evaluation of ecological protection value, the measurement of forest carbon reserves, and the performance evaluation of the tenure of leading cadres by the superior departments (such as cities and provinces). It has become an important research direction at home and abroad using remote sensing data, such as satellites and unmanned aerial vehicles, combining with less sampling plot to effectively estimate the overall parameters. At present, there is no comparative validation of several estimation methods for mean stock volume based on multi-source data in domestic. In order to promote the application of remote sensing technology in the forest resource survey, there is an urgent need to compare and evaluate the methods for estimating overall parameters of forest farms or counties.Method The main plantation tree species of Wangyedian Forest Farm in Inner Mongolia of northern China were taken as the research object. Based on the sampling UAV LiDAR (herringbone system distribution), Sentinel-2A data (full coverage) and a small amount of sample plot data obtained in 2019, four patterns of sample plot (p), sample plot-satellite (ps), sample plot-sampling LiDAR (pl), and sample plot-sampling LiDAR-full coverage satellite (pls) were estimated and compared with five methods, including DBp, MDps, MDpls, HYpl and MAps, which were suitable for these four patterns and belong to design-based method (DB), model-assisted method (MA), model-dependent method (MD) and mixed or hybrid method (HY).Result (1)The mean stock volume densities of DBp, MDps, MAps, HYpl, and MDpls were 212.54, 202.09, 202.38, 210.07 and 208.96 m3/ha, respectively. The accuracy (P) was 90.44%, 91.54%, 91.69%, 96.35%, and 96.45%, respectively, with variances decreasing in turn. (2) The relative efficiency (RE) of other methods compared with MDpls was greater than 1 (REDBp,MDpls = 5.39, REMDps,MDpls = 3.82, REMAps,MDpls = 3.69, REHYpl,MDpls = 1.07), and the RE of HYpl compared with MDpls method was greater than 1, but close to 1. Compared with the other three methods, HYpl was more efficient (REDBp,HYpl = 5.02, REMAps,HYpl = 3.43, REMDps,HYpl = 3.56). (3) Both HYpl and MDpls methods containing LiDAR data had positive efficiency compared with MDps and MAps methods containing Sentinel-2A data (REMAps,HYpl = 3.43, REMDps,HYpl = 3.56, REMDps,MDpls = 3.82, REMAps,MDpls = 3.69). The RE of MDps and MAps was close to 1, but the efficiency of MAps was slightly higher than that of MDps (REMDps,MAps = 1.04).Conclusion Compared with the estimation method that only used the sample plot data, when the remote sensing data was used as an auxiliary variable to establish the estimation model, it can effectively improve the estimation accuracy of the MSVD of the forest farm, whether the Sentinel-2A multispectral remote sensing data, which is fully covered the forest farm but not sensitive enough to the stock volume, or the sampling LiDAR data which are sensitive to stock volume. Among four methods involving the application of remote sensing data, MDpls has the highest estimation accuracy, followed by HYpl. Both of these two methods including the application of LiDAR remote sensing sampling observation data are suitable for the MSVD estimation of forest farms. The MDpls method has the highest estimation accuracy and relative efficiency, which can realize the synergistic application of the Space-Air-Earth multi-source observation data, and can be used as the preferred method for annual monitoring of forest stock in forest farms.
-
以林场或县（后文统一称林场）森林资源总体为调查对象，及时、准确地调查监测总体平均每公顷蓄积量（mean stock volume density，MSVD），对上级（如市、省）部门开展森林资源宏观管理、生态保护价值评价、森林碳储量计量、领导干部任期绩效考核等工作都有重要支撑作用[1]。现有的森林资源一类、二类调查为实现林场总体MSVD估计提供了工作基础和数据基础[2-4]，但并不能满足当前对林场森林资源实现年度或按需监测（在二类调查的间隔期内任何一个年度完成调查）的国家需求。一类调查所依赖的系统抽样固定样地数据以省为总体设计，用于省、国家级总体出数可达到足够精度[5]，但若以林场为总体，只依赖一类调查固定样地无法达到精度要求。二类调查通过对每个小班进行角规或标准样地调查等可以得到每个小班的蓄积量，进而汇总计算得到总体MSVD，但这种调查方式工作量很大，通常间隔10年调查一次，不适合年度监测。根据二类调查技术规范[6]，为了控制小班蓄积量调查成果质量，需要通过系统、随机等概率抽样布设抽样样地，通过抽样样地蓄积量估计得到林场总体蓄积量。这是传统的概率抽样调查方法，但要做到年度出数或按需监测，就特别需要更加经济、高效的调查估计方法。而将辅助数据和样地数据综合用于总体参数的估计就是一种常用的方法，且各种遥感数据是最常用的辅助数据[7]。
Ståhl等[8]将森林资源调查方法归纳为4类：概率抽样法（design-based method，DB）、模型辅助法（model-assisted method，MA）、模型法（model-dependent method，MD）和混合法（mixed or hybrid method，HY）。DB依概率从总体单元中抽取样本，根据样本观测值估计总体参数，理论上具有设计无偏性。但由于DB只采用抽样样地数据进行推断，不采用辅助数据，因此在同样抽样精度下需要抽取更多的样本单元[8]。MA本质上属于DB方法，对其改进在于利用辅助数据建立估测模型，用残差校正估测模型的总体参数估计结果，若模型拟合得当，可以提高估计精度，但建立模型所用的样本必需采用概率抽样法获得[9]。MD将样本单元观测值作为因变量，将辅助数据（空间上对总体单元全覆盖）作为自变量，建立回归估测模型（如线性、非线性统计回归模型、各种机器学习模型等），通过辅助变量和估测模型预测每个总体单元的目标参数值，这些值的平均值就为总体参数估计值[10]。HY是指将DB和MD混合应用的方法[11-12]。HY的一种典型应用场景就是当MD方法所用的遥感辅助数据不是对总体全覆盖，而是采用概率抽样方式获取时。这时仍然需要建立一个估测模型，但对总体目标参数估计值的不确定性的估计则需要考虑遥感辅助数据是依概率抽样获取所引起的不确定性。国外已对以上4类方法开展了大量研究，包括对各种方法的比较评价[7-9,13-19]，所涉及的天（卫星）空（航空）地（样地）多源观测数据协同应用模式主要包括：样地−卫星遥感、样地−无/有人机遥感和样地−无/有人机遥感−卫星遥感等。
国内对这4类方法的研究也都有涉及，但很不均衡。DB是国家森林资源调查的基础方法[20]，对其研究较多，如曾伟生等[21]对群团抽样估计效率与群内样地间距的关系进行了研究，对两阶段群团抽样在森林调查中的综合效率进行了讨论；李春干等[22]通过布设抽样遥感样地，通过目视判读确定样地森林资源参数，进而估算总体参数；郑冬梅等[23]对比分析了群团样地判读和图班区划判读的精度，卫星数据采用了资源三号和高分一号数据，验证了两者对森林总体面积估计的差异较小。国内对MA开展研究很少，而且局限于地类面积的成数估计[24-25]，尚未见以遥感作为辅助变量采用MA估计总体蓄积量、平均高等参数的研究。国内森林参数遥感定量估测方向的研究很多[26-30]，所采用的方法绝大多数都属于MD类方法。在HY研究方面，国内主要采用了双重回归估计方法。如宋新民[31]介绍了双重回归抽样法，在方差估计上既考虑了回归估测模型误差，又考虑了抽样误差；针对二类调查业务的HY研究，遥感数据主要以航片为主[32-34]；在宏观监测应用方面，陈振雄等[35]基于遥感大样地双重抽样方法，以遥感样地为一重样本，用于区划判读，地面样地为二重样本，结合遥感样地判读结果，采用双重回归估计法对广东省主要地类面积进行了估计。国内相关研究所涉及的天空地多源数据协同应用模式以样地−卫星遥感、样地−无/有人机遥感为主，采用样地−无人机遥感−卫星遥感模式的较少。目前国内尚未见对现有4类估测方法的比较评价研究。
总之，国际上对区域总体森林参数的估计已形成一套方法体系，合理地利用各种方法将可最大限度发挥天空地多源观测数据协同应用优势，提高林场总体蓄积量的调查监测效率。为了促进新一代遥感技术在我国森林资源调查监测业务中的应用，提高森林资源天空地一体化调查监测技术水平，亟需对现有估测方法进行比较评价研究。
本研究以内蒙古自治区赤峰市喀喇沁旗旺业甸实验林场为研究区，以林场主要人工林的总体MSVD估计为目标，获取了标准样地数据、无人机LiDAR抽样数据，以及对林场全覆盖的Sentinel-2A卫星多光谱遥感数据，对适合以上数据应用的4类5种估测方法进行了比较评价，以期为国内林场总体蓄积量的年度监测提供技术参考。
1. 研究区和数据
1.1 研究区概况
研究区为内蒙古自治区赤峰市喀喇沁旗旺业甸实验林场所覆盖区域（图1），总面积约500 km2，地理坐标为118°09′ ~ 118°30′ E，41°21′ ~ 41°39′N，地势西南高东北低，以中山和低山山地为主，平均海拔800 ~ 1 890 m，年均降水量300 ~ 500 mm，年平均气温4.2 ℃，为大陆季风性气候。土壤类型以草甸土、棕壤、褐土和山地黑土为主。林场的人工林树种以落叶松（Larix gmelinii）、油松（Pinus tabuliformis）、樟子松（Pinus sylvestris）、云杉（Picea asperata）和红松（Pinus koraiensis）为主，天然次生林树种以白桦（Betula platyphylla）、山杨（Populus davidiana）和蒙古栎（Quercus mongolica）为主[36]。
1.2 数据获取与处理
研究所用数据包括卫星遥感数据、无人机LiDAR数据、标准样地数据、无人机CCD数码影像，以及林场的数字化林相图。图2给出了所采用的天空地数据获取方式。卫星遥感数据全覆盖研究区。无人机遥感观测区域和标准样地采用系统抽样方法布设，具体做法为：将研究区划分为22个4 km × 4 km大小的网格，左上角第一个网格的中心点的大地坐标采用随机方式确定，原则上以各网格的中心点作为无人机抽样中心点位置，按照“人”字形状获取3个长方形无人机LiDAR条带数据，3个条带在一端重合，组成一块LiDAR抽样观测数据；再在每个LiDAR抽样块覆盖的人工林内布设4块正方形标准样地，其中1块位于抽样块的中心位置，另外3块位于3个LiDAR条带上，每个条带1块（图2A）。同时还获取了高空间分辨率（0.05 m × 0.05 m）无人机CCD数码相机影像，每个LiDAR抽样块的覆盖范围如图2B所示。
![]()
图 2 无人机和样地数据抽样获取方案示意图Figure 2. Schematic diagram of UAV and sample plot data sampling acquisition scheme1.2.1 卫星遥感数据
卫星数据采用Sentinel-2A（S-2A）多光谱数据，需要2景才能对林场实现全覆盖。这2景数据的成像时间都为2019年10月7日。每景数据包括13个波段，不同波段范围空间分辨率不同，共含10、20和60 m 3种空间分辨率，本研究只采用了其中的10和20 m空间分辨率的多光谱数据，共9个波段（表1）。对每景数据分别进行了辐射定标、大气校正和正射校正等预处理，其中正射校正所需地面控制点的选取参考了正射校正后的无人机CCD数码影像，校正精度在一个像元以内。经以上处理后的S-2A数据的像元大小设为25 m × 25 m（和样地大小一致）；然后对这2景预处理后的数据进行镶嵌处理，进而利用林场矢量边界对镶嵌数据进行裁剪，得到覆盖整个林场的多光谱遥感数据。
表 1 S-2A卫星数据特征Table 1. S-2A satellite data characteristics特征名称
Feature name特征符号
Feature symbol计算公式
Calculation formula特征名称
Feature name特征符号
Feature symbol计算公式
Calculation formula光谱特征
Spectral featureB2、B3、B4、B5、B6、B7、B8a、B11、B12 角二阶矩
Angular second momentAN k∑i,j=1P2i,j 差值植被指数
Difference vegetation indexDVI B8−B4 熵
EntropyEN k∑i,j=1Pi,j(−lnPi,j) 增强植被指数
Enhanced vegetation indexEVI 2.5 × (B8 − B4)/(B8 + 6 × B4 − 7.5 × B2 + 1) 对比度
ContrastCON k∑i,j=1(i−j)2×Pi,j 归一化植被指数
Normalized difference vegetation indexNDVI (B8−B4)/(B8+B4) 均值
MeanME 1k2k∑i,j=1Pi,j 比值植被指数
Ratio vegetation indexSR B8/B4 方差
VarianceVAR k∑i,j=1Pi,j(i,j−ME) 转换归一化植被指数
Transformed normalized difference vegetation indexTNDVI √(B8−B4)/(B8+B4)+0.5 同质性
HomogeneityHOM k∑i,j=1Pi,j1+(i−j)2 叶绿素指数
Chlorophyll indexCIg B8/B3−1 相关性
CorrelationCOR k∑i,j=1P2i,j[(i−ui)−(j−uj)√σ2iσ2j] 反向红边叶绿素指数
Inverted red edge chlorophyll indexIRECI (B7−B4)/(B5/B6) 相异性
DissimilarityDIS k∑i,j=1Pi,j|i−j| 色素简单比值指数
Pigment specific simple ratioPSSRA B7/B4 修正叶绿素吸收反射指数
Modified chlorophyll absorption in reflectance indexMCARI [(B5 − B4) − 0.2 × (B5 − B3)] × (B5/B4) S-2A 红边位置指数
S-2A red edge positionS2REP 705 + 35 × [(B4 + B7)/2 − B5]/(B6 − B5) 修正窄边红边简单比值指数
Modified simple ratio red edge narrowMSRren [(B8a/B5)−1]/√(B8a/B5)+1 红边叶绿素指数
Red edge chlorophyll indexCIgre1 B5/B3−1 红边植被指数
Red edge vegetation indexNDVIre1 (B8−B5)/(B8+B5) CIgre2 B5/B3−1 NDVIre2 (B8−B6)/(B8+B6) CIgre3 B7/B3−1 NDVIre3 (B8−B7)/(B8+B7) 注：B2、B3、B4、B5、B6、B7、B8、B8a、B11、B12代表S-2A数据对应的波段，Pi,j=Di,j/∑ki,j=1Di,j，Di,j为i行j列对应的像元值，k代表计算纹理时窗口大小。Notes: B2, B3, B4, B5, B6, B7, B8, B8a, B11 and B12 represent the bands corresponding to S-2A data; Pi,j=Di,j/∑ki,j=1Di,j, Di,j is the pixel value corresponding to the i row and j column, k represents the window size when calculating the texture. 以S-2A卫星遥感数据为辅助数据建立估测模型时主要考虑了光谱反射率、光谱指数和纹理3类特征，每类包含的具体遥感特征名称、符号和计算公式见表1。纹理特征的计算都采用3 × 3窗口。这34个特征为估测模型建立的备选特征。
由于人力、物力限制，本研究只对优势树种为油松、落叶松和樟子松的人工林进行了样地调查，样地只对林场针叶人工林具有代表性，只适合用于林场针叶人工林总体MSVD的估测。因此林场总体单元只包括优势树种为针叶人工林的小班。采用林场小班矢量数据库中记录的优势树种属性信息，从卫星遥感辅助数据中去除优势树种非针叶人工林的像元，剩下的像元组成林场总体单元集U，包含的总体单元数N为159 535。
1.2.2 无人机数据
2019年9—10月完成了无人机LiDAR条带抽样获取，采用的无人机平台为8轴RC6-2，LiDAR系统为Riegl VUX-1激光雷达系统，飞行高度为58 ~ 340 m。由于空管政策限制，北部3个网格内的数据未能按计划获取，最终只获取了19块无人机LiDAR抽样观测点云数据，平均点云密度为107 ~ 152 pts/m2，垂直精度为−0.36 ~ 0.19 m，水平精度为0.25 m。原始LiDAR点云数据经姿态矫正、噪声点剔除、坐标转换、航带拼接、系统差改正等预处理后，利用TerraSolid软件进行点云分类、去噪以及归一化等处理得到归一化植被点云数据，并对强度信息进行了校正[37]。利用Fusion软件[38]提取常用的高度和强度相关特征参数，主要包括均值、方差、最大值、最小值、标准差、变异系数、四分距差、偏斜度、百分位数等[39]，如表2所示。这些特征的空间分辨率也设为25 m × 25 m。由于无人机飞行区域边界处有数据缺失值问题，将有效区域的边界向内收缩了一部分（图2-A）。无人机有效数据区域内的所有像元构成Sa集合，其中每个像元都和Ｕ的某个单元相对应。Sa包含像元（单元）总数M为11 902。
表 2 LiDAR数据特征Table 2. Feature of LiDAR data特征名称
Feature name特征符号
Feature symbol均值
MeanHmean, Imean 方差和标准差
Variance and standard deviationHvar, Ivar, Hsd, Isd 最大值和最小值
Maximum value and minimum valueHmax, Imax, Hmin, Imin 变异系数
Coefficient of variationHcv, Icv 四分距差
Interquartile distanceHiq, Iiq 偏斜度
SkewnessHsk, Isk 百分位数
PercentileHp05, Hp10, Hp20, Hp25, ···, Hp95, Hp99
Ip05, Ip10, Ip20, Ip25,···, Ip95, Ip99最小高度以上返回点
Count of return point above the minimum heightR1H, R2H,···, R8H, R9H 注：H代表高度，I代表强度，var代表方差，sd代表标准差，max与min分别代表最大值与最小值，cv、iq以及sk分别代表变异系数、四分距差以及偏斜度，p05—p99代表对应的百分位数。Notes: H represents height, I represents intensity, var represents variance, sd represents standard deviation, max and min represent maximum value and minimum value, respectively; cv, iq and sk represent coefficient of variation, interquartile distance and skewness, respectively; p05−p99 represent the corresponding percentiles. 1.2.3 标准样地数据
在获取无人机数据的同时，进行了标准样地调查。样地大小为25 m × 25 m，对样地内胸径大于5 cm的树木进行每木检尺，测量胸径、树高、东西冠幅、南北冠幅及枝下高等参数。为了减少后期遥感影像与样地位置偏差带来的影响，样地位置主要布设在范围较大且均匀的林分中心，样地4个角点位置及样地中心位置采用Zenith15R型GPS-RTK和ZT20型全站仪精确测量，测量精度优于0.15 m。共调查了88块样地，其中76块位于无人机抽样数据块内（图2）。样地优势树种类型包括落叶松、油松和樟子松，有些样地内也会出现少数白桦、黑桦（Betula davurica）等林木。根据表3所列材积计算公式[40]计算样地内所有单木材积，合计得到样地总蓄积量和每公顷蓄积量（蓄积量密度）。所有标准样地蓄积量密度构成样地数据集S，所含样本单元数n为76。
表 3 材积计算公式Table 3. Formula for volume calculation树种 Tree species 材积计算公式 Formula for volume calculation 白桦
Betula platyphyllaV=100.000397507075980248×d2×h+3.8121380807356×10−6×d3×h−0.000843446660194932×d2−0.000290088083727618×d2×h×lgd 黑桦
Betula davuricaV=0.0199215193890509+0.00038814516708027×d2−2.05059660776977×10−5×d2×h−
0.000870310875746131×h2+8.05362136949543×10−5×d×h2落叶松
Larix gmeliniiV=10−3.49890390728004+2.75504846502564×lgd−0.394050839410844×lgd2−1.37915356404294×lgh+1.14586681477751×lgh2 樟子松
Pinus sylvestrisV=0.000445504541007861×d1.66316523786429×exp(0.0989929357004981×h−1.66681646056217/h) 油松
Pinus tabuliformisV=2.53309290763168×10−5×d2×h−8.57378841160325×10−7×d3×h−6.00033429201054×10−5×d2+
2.93720677218884×10−5×d2×h×lgd注：引自文献[40]。V代表单木材积；d代表单木胸径；h代表单木树高。Notes: cited from reference [40]. V represents single tree stock volume, d represents single tree DBH, h represents single tree height. 2. 总体MSVD估计方法和评价指标
2.1 总体MSVD估计方法
以林场针叶人工林为总体，比较评价如下5种方法对总体MSVD的估计效果：
DBp：仅利用样地数据S，采用DB进行估计。这里下标p表示样地。
MDps：利用样地数据S和卫星遥感特征，建立估测模型，模型因变量为标准样地蓄积量密度，自变量为卫星遥感特征。这里，下标ps表示采用的数据包括标准样地蓄积量密度和卫星遥感特征，s表示卫星。
MAps：所用数据和估测模型与MDps相同。与MDps相比，其结果需用标准样地蓄积量密度与对应模型估测结果的残差校正估测模型对总体参数的估计结果。
HYpl：利用样地数据集S和无人机LiDAR遥感特征，用下标pl表示，其中l表示LiDAR；所采用的方法为HY，估测模型以标准样地蓄积量密度为因变量，以LiDAR遥感特征为自变量；和MDps相比，主要区别在总体平均值方差估计上，需要考虑LiDAR数据是按概率抽样获取的，而不是对总体全覆盖。
MDpls：利用了标准样地数据、无人机LiDAR遥感特征和卫星遥感特征；采用分层模型估测法（HMB），这是一种可综合应用天空地数据的估计方法，属于MD类方法。
图3给出了以上5种方法所利用的数据情况。该图最上面表示图例，最下面给出了5种方法的缩写名称，中间是对应每种方法所采用的数据类型，以及各数据类型之间的空间对应关系。
![]()
图 3 各估测方法所利用数据示意图Figure 3. Schematic diagram of the data used by each estimation method2.1.1 DBp
该方法直接采用标准样地数据进行总体均值的估计，根据系统抽样原理，其均值公式如（1）所示[17]。
¯uDB=1n∑i∈Syi (1) 式中：S为样地数据集，
¯uDB 代表总体均值，yi 为样本单元蓄积量密度（样地），n为样本单元数。对于有限总体，系统抽样的均值方差如（2）所示。
V(¯uDB)=1−fn2∑i∈S(yi−¯y)2 (2) 式中：
1−f=N−nN 为有限总体校正系数，¯y 为样本单元均值，N为总体单元数。2.1.2 MDps
该方法利用了全覆盖的卫星遥感特征与标准样地数据，总体均值的估计方程如（3）所示。
{\widehat{\mu }}_{{\mathrm{M}\mathrm{D}}_{\mathrm{p}\mathrm{s}}}={\boldsymbol{\iota }}_{U}^{\mathbf{T}}{\boldsymbol{Z}}_{U}{\widehat{\boldsymbol{\alpha }}}_{S} (3) 式中：
{\widehat{\boldsymbol{\alpha }}}_{{S}} 为模型参数，{\boldsymbol{\iota }}_{U} 为N维向量，每个元素维1/N，{\widehat{\boldsymbol{\alpha }}}_{S}=({{\boldsymbol{Z}}_{S}^{\mathbf{T}}{\boldsymbol{Z}}_{S})}^{-1}{\boldsymbol{Z}}_{S}^{\mathbf{T}}{\boldsymbol{y}}_{S} ，{\boldsymbol{Z}}_{S} 为样地对应的卫星遥感特征，{\boldsymbol{Z}}_{U} 为总体所有单元对应的卫星遥感特征，{\boldsymbol{y}}_{S} 为为样地数据集S的标准样地蓄积量密度。方差估计公式如（4）所示。
\begin{split}\\ V\left({\widehat{\mu }}_{{\mathrm{M}\mathrm{D}}_{\mathrm{p}\mathrm{s}}}\right)={\boldsymbol{\iota }}_{U}^{\mathbf{T}}{\boldsymbol{Z}}_{U}\widehat{\mathrm{C}\mathrm{o}\mathrm{v}}\left({\widehat{\boldsymbol{\alpha }}}_{S}\right){\boldsymbol{Z}}_{U}^{\mathbf{T}}{\boldsymbol{\iota }}_{U} \end{split} (4) 式中：
V\left({\widehat{\mu }}_{{\mathrm{M}\mathrm{D}}_{\mathrm{p}\mathrm{s}}}\right) 为均值方差，\widehat{\mathrm{C}\mathrm{o}\mathrm{v}}\left({\widehat{\boldsymbol{\alpha }}}_{S}\right)= \dfrac{{({\boldsymbol{y}}_{S}-{\boldsymbol{Z}}_{S}{\widehat{\boldsymbol{\alpha }}}_{\boldsymbol{S}})}^{\mathbf{T}}({\boldsymbol{y}}_{S}-{\boldsymbol{Z}}_{S}{\widehat{\boldsymbol{\alpha }}}_{S})}{n-p-1} ({{\boldsymbol{Z}}_{S}^{\mathbf{T}}{\boldsymbol{Z}}_{S})}^{-1} 代表模型参数协方差矩阵，p为遥感特征数[41]。2.1.3 MAps
该方法用MD估计的残差校正估测模型对总体参数的估测值，总体均值计算公式如（5）所示[42]。
{\widehat{\mu }}_{\mathrm{M}\mathrm{A}}=\frac{1}{n}\left({\sum} _{i\in S}{(y}_{i}-{\widehat{y}}_{i})\right) + \frac{1}{N}{\sum} _{i\in U}{\widehat{y}}_{i} (5) 式中：
{\widehat{\mu }}_{\mathrm{M}\mathrm{A}} 为总体均值，{\widehat{y}}_{i} 为i单元对应的估测值；总体均值方差的计算公式如（6）所示。V\left({\widehat{\mu }}_{\mathrm{M}\mathrm{A}}\right)=(1-\frac{n}{N})\frac{1}{n}\frac{1}{n-1}{\sum} _{i\in S}{{(y}_{i}-{\widehat{y}}_{i})}^{2} (6) 式中：
V\left({\widehat{\mu }}_{\mathrm{M}\mathrm{A}}\right) 为均值方差。2.1.4 HYpl
该方法采用地面样地调查数据和采用抽样方法获取的遥感数据估算总体参数[8]。总体均值计算公式如（7）所示[43]。
{\widehat{\mu }}_{\mathrm{H}\mathrm{Y}}={\boldsymbol{\iota }}_{{S}_{{\rm{a}}}}^{\mathbf{T}}{\boldsymbol{X}}_{{S}_{{\rm{a}}}}{\widehat{\boldsymbol{\beta }}}_{S} (7) 式中：Sa为无人机有效数据区域内的所有像元构成的集合。
{\boldsymbol{\iota }}_{{S}_{{\rm{a}}}} 为长度为M的列向量，元素为1/M，M为Sa包含的单元数。{\boldsymbol{X}}_{{S}_{{\rm{a}}}} 为M × （p + 1）的无人机遥感特征，{\widehat{\boldsymbol{\beta }}}_{S}={\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1}{\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{y}}_{S} ，为标准样地蓄积量密度{\boldsymbol{y}}_{S} 与其对应无人机LiDAR遥感特征{\boldsymbol{X}}_{S} 模拟的线性模型参数。计算方差的公式如（8）所示[15]。
\begin{split} V\left({\widehat{\mu }}_{\mathrm{H}\mathrm{Y}}\right)=&\;\left(1-\frac{M}{N}\right)\frac{1}{M}\frac{1}{M-1}{\sum }_{i=1}^{M}{{\Big(\widehat{y}}_{i}-\overline{\widehat{y}}\Big)}^{2} + \\&\;{\boldsymbol{\iota }}_{{S}_{{\rm{a}}}}^{\mathbf{T}}{\boldsymbol{X}}_{{S}_{{\rm{a}}}}\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\boldsymbol{\beta }}}_{S}\right){\boldsymbol{X}}_{{S}_{{\rm{a}}}}^{\mathbf{T}}{\boldsymbol{\iota }}_{{S}_{{\rm{a}}}} \end{split} (8) 式中：
\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\boldsymbol{\beta }}}_{S}\right)=\widehat{{\sigma }_{e}^{2}}{\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1} 为估计参数{\widehat{\boldsymbol{\beta }}}_{S} 的协方差矩阵，\widehat{{\sigma }_{e}^{2}}=\dfrac{{\hat{\boldsymbol{e}}}_{S}^{\mathrm{T}}{\hat{\boldsymbol{e}}}_{S}}{m-p-1} 为估计残差方差，{\hat{\boldsymbol{e}}}_{S}={\boldsymbol{y}}_{S}- {\boldsymbol{X}}_{S}{\hat{\boldsymbol{\beta }}}_{S} ，为样本集S的残差，长度为m的列向量，m为抽样单元数。在异方差的情况下，协方差矩阵\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\boldsymbol{\beta }}}_{S}\right)= {\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1}\left[\displaystyle {\sum }_{i=1}^{m}{\hat{e}}_{i}^{2}{{\boldsymbol{x}}}_{i}^{\mathbf{T}}{x}_{i}\right]{\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1} ，{{\boldsymbol{x}}}_{i} 为样本集中第i个样本对应的无人机遥感特征，同时，对残差平方{\hat{e}}_{i}^{2} 进行修正，修正参数为\dfrac{m}{m-p-1} [41]。2.1.5 MDpls
Saarela等[7]于2016年提出了分层建模估计方法，该方法综合利用了空天地3种数据源实现参数估算，总体均值的期望值估计公式如（12）所示。
{\widehat{\mu }}_{{\mathrm{M}\mathrm{D}}_{\mathrm{p}\mathrm{l}\mathrm{s}}}={\boldsymbol{\iota }}_{U}^{\mathbf{T}}{\boldsymbol{Z}}_{U}{\left({\boldsymbol{Z}}_{{S}_{\rm{a}}}^{\mathbf{T}}{\boldsymbol{Z}}_{{S}_{\rm{a}}}\right)}^{-1}{\boldsymbol{Z}}_{{S}_{\rm{a}}}^{\mathbf{T}}{\boldsymbol{X}}_{{S}_{\rm{a}}}{\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1}{\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{y}}_{S} (12) 式中：
{\widehat{\mu }}_{{\mathrm{M}\mathrm{D}}_{\mathrm{p}\mathrm{l}\mathrm{s}}} 表示分层建模估计的总体均值，{\boldsymbol{\iota }}_{U} 是一个N维向量，每个元素为1/N ，{\boldsymbol{Z}}_{U} 为所有单元对应的卫星数据特征，{\boldsymbol{Z}}_{{S}_{{\rm{a}}}} 为Sa样本集抽样样本对应的卫星遥感特征，{\boldsymbol{X}}_{{S}_{{\rm{a}}}} 为Sa样本集抽样样本对应的无人机遥感特征，{\boldsymbol{X}}_{S} 为样本集S对应的无人机遥感特征。均值方差的估计如（13）所示。V\left({\widehat{\mu }}_{{\mathrm{M}\mathrm{D}}_{\mathrm{p}\mathrm{l}\mathrm{s}}}\right)={\boldsymbol{\iota }}_{U}^{\mathbf{T}}{\boldsymbol{Z}}_{U}\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\boldsymbol{\alpha }}}_{{S}_{{\rm{a}}}}\right){\boldsymbol{Z}}_{U}^{\mathbf{T}}{\boldsymbol{\iota }}_{U} (13) 式中：
\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\boldsymbol{\alpha }}}_{{S}_{{\rm{a}}}}\right) 为模型参数{\widehat{\boldsymbol{\alpha }}}_{{S}_{{\rm{a}}}} 的协方差矩阵，计算公式如（14）所示。\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\boldsymbol{\alpha }}}_{{S}_{\rm{a}}}\right) = \frac{{({\boldsymbol{X}}_{{S}_{\rm{a}}}{\widehat{\boldsymbol{\beta }}}_{s} - {\boldsymbol{Z}}_{{S}_{\rm{a}}}{\widehat{\boldsymbol{\alpha }}}_{{S}_{\rm{a}}} )}^{\mathbf{T}}({\boldsymbol{X}}_{{S}_{\rm{a}}}{\widehat{\boldsymbol{\beta }}}_{S} - {\boldsymbol{Z}}_{{S}_{\rm{a}}}{\widehat{\boldsymbol{\alpha }}}_{{S}_{\rm{a}}} )}{M-( q + 1 )}{\left( {\boldsymbol{Z}}_{{S}_{\rm{a}}}^{\mathbf{T}}{\boldsymbol{Z}}_{{S}_{\rm{a}}} \right)}^{-1} + {\left({\boldsymbol{Z}}_{{S}_{ \rm{a}}}^{\mathbf{T}}{\boldsymbol{Z}}_{{S}_{ \rm{a}}}\right)}^{-1}{\boldsymbol{Z}}_{{S}_{ \rm{a}}}^{\mathbf{T}}\left[{\boldsymbol{X}}_{{S}_{ \rm{a}}}{\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1}{\boldsymbol{X}}_{{S}_{ \rm{a}}}^{\mathbf{T}}\right]{\boldsymbol{Z}}_{{S}_{ \rm{a}}}{\left({\boldsymbol{Z}}_{{S}_{ \rm{a}}}^{\mathbf{T}}{\boldsymbol{Z}}_{{S}_{ \rm{a}}}\right)}^{-1} (14) 式中：
{\widehat{\boldsymbol{\beta }}}_{S}={\left({\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{X}}_{S}\right)}^{-1}{\boldsymbol{X}}_{S}^{\mathbf{T}}{\boldsymbol{y}}_{S} 为模型参数。2.2 遥感特征选择方法
后4个估测方法都涉及以遥感特征为自变量的统计回归模型建立问题，而无论是卫星遥感数据，还是激光雷达数据都有比较多的备选特征。这些备选特征之间会存在共线性问题，需要进行特征的优选。使用R统计软件VSURF包分别选择LiDAR遥感特征、卫星遥感特征。该包提供了一种自动的方法来根据预测变量的重要性分数选择预测变量[44]，然后结合方差膨胀因子（VIF）以及全子集模型优化策略，利用交叉验证方法选择误差最小的特征组合。
2.3 估测效果评价指标
2.3.1 像元尺度估测效果评价指标
为评价后4个估测方法中所涉及的遥感估测模型的效果，采用留一交叉验证的方法对结果进行评价，选用平均误差（ME）、决定系数（R2）、均方根误差（RMSE）、相对均方根误差（RRMSE）作为像元尺度平均估测效果评价指标，计算公式见（15） ~ （18）。
\mathrm{M}\mathrm{E}=\frac{1}{n}{\sum} _{i=1}^{n}\left({y}_{i}-{\widehat{y}}_{i}\right) (15) {R}^{2}=1-\frac{\displaystyle {\sum} _{i=1}^{n}{\left({y}_{i}-{\widehat{y}}_{i}\right)}^{2}}{\displaystyle {\sum} _{i=1}^{n}{\left({y}_{i}-\overline{y}\right)}^{2}} (16) \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{\displaystyle {\sum} _{i=1}^{n}{\left({y}_{i}-{\widehat{y}}_{i}\right)}^{2}}{n-1}} (17) \mathrm{R}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\left(1-\frac{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}{\overline{y}}\right)\times 100 \% (18) 式中：
{y}_{i} 为样地蓄积量密度实测值；\overline{y} 为实测值对应的均值；{\widehat{y}}_{\mathrm{i}} 为估测值；n为样本点的个数。2.3.2 林场总体MSVD估计精度与相对效率
对5种方法总体结果的精度评价采用标准误（
{S}_{\overline{e}} ）和相对估测精度（P）指标，计算公式如（19）、（20）所示[14,45]。{S}_{\overline{\mathrm{e}}}=\sqrt{V\left(\overline{u}\right)} (19) P=\left(1-\frac{t \cdot {S}_{\overline{\mathrm{e}}}}{\overline{y}}\right)\times 100\% (20) 式中：
V\left(\overline{u}\right) 为均值方差，\overline{y} 为均值，t为可靠性指标，按95%置信区间，t = 1.96。将相对效率（relative efficiency，RE）作为评价不同总体结果估测方法相对效果的指标[14,45]，计算公式如（21）所示。
\mathrm{R}\mathrm{E}=\frac{V\left({\overline{u}}_{1}\right)}{V\left({\overline{u}}_{2}\right)} (21) 式中：
V\left({\overline{u}}_{1}\right) 和V\left({\overline{u}}_{2}\right) 代表两种方法对应的均值方差，当RE > 1，说明V\left({\overline{u}}_{1}\right) 对应的方法效率低于V\left({\overline{u}}_{2}\right) 对应的方法；当RE = 1，说明两种方法效率相同；当RE < 1，说明V\left({\overline{u}}_{1}\right) 对应的方法效率高于V\left({\overline{u}}_{2}\right) 对应的方法。3. 结果与分析
3.1 遥感特征选择结果
利用R语言VSURF包分别对S-2A卫星遥感特征、LiDAR遥感特征进行优选。对于无人机LiDAR数据，当特征数为3时，其RMSE最小（图4a），此时，对应特征为Hmean、Ip99以及R2H；对于S-2A数据，当特征数为4时，其RMSE最小（图4b），对应特征为B2、B2ME、MCARI指数以及B12COR。特征选择后保留的遥感特征见表4。
![]()
图 4 遥感特征选择折线代表值的变化趋势。The broken line represents the trend of value.Figure 4. Feature selection of remote sensing表 4 遥感特征优选后保留的特征Table 4. Remotes sensing features selected after optimum feature selection数据源
Data source特征选择结果
Result of feature selectionS-2A B2，MCARI，B2ME，B12COR LiDAR Hmean，Ip99，R2H 3.2 遥感估测模型的像元尺度估测精度
以S-2A多光谱遥感特征为自变量，以样地蓄积量密度为因变量建立的估测模型（S-V）的精度验证结果见图5b，决定系数R2为0.21，RMSE为81.23 m3/hm2，ME为−0.25 m3/hm2，RRMSE为59.81%。而以LiDAR遥感特征为自变量，以样地蓄积量密度为因变量建立的估测模型（L-V）的精度验证结果见图5a，决定系数R2为0.86，RMSE为33.12 m3/hm2，ME为−0.06 m3/hm2，RRMSE为84.42%。精度指标相比S-2A，R2提高了0.65，RMSE降低了48.11 m3/hm2，ME降低了0.19 m3/hm2，RRMSE提高了24.61%。以LiDAR估测的蓄积量密度为因变量，以S-2A多光谱遥感特征为自变量建立估测模型（S-L-V），估测精度验证结果见图5c，决定系数R2仅提高了0.07，估测精度远低于以LiDAR特征为自变量建立的模型。
![]()
图 5 遥感特征建模估测模型精度实线为1∶1验证线。图中RMSE和ME的单位均为m3/hm2。The solid line is 1∶1 verification line.The unit of RMSE and ME in the figure is m3/ha.Figure 5. Remote sensing feature modeling to estimate model accuracy总之，只要以S-2A多光谱遥感特征作为自变量建立模型估测森林蓄积量密度，无论因变量是样地蓄积量密度还是LiDAR估测得到的蓄积量密度，估测精度都很低。这主要是由于卫星多光谱遥感特征对森林蓄积量不太敏感，两者之间的相关性很低，而且光学遥感探测森林结构容易出现信号饱和问题。该结果与已有研究[15,46-47]的结论是一致的。
3.3 总体MSVD估计精度
如表5所示：5种方法总体结果有一定差别，DBp、MDps、MAps、HYpl、MDpls 5种方法估测总体蓄积均值分别为212.54、202.09、202.38、210.07以及208.96 m3/hm2，仅利用样地的DBp方法结果方差最大（107.51 m3/hm2），当加入全覆盖S-2A数据时，MDps与MAps方法的方差均减少，减少了31.33与33.96 m3/hm2，精度提高1.10%与1.25%。
相对于样地−全覆盖S-2A数据协同模式，采用样地-LiDAR抽样数据协同模式的方法精度更高，HYpl方法相比MDps与MAps方法，方差分别减少了54.76与52.13 m3/hm2，精度提高了4.81%与4.66%。
采用样地−抽样LiDAR数据−全覆盖S-2A数据协同模式的MDpls方法精度最高，均值方差为19.95 m3/hm2，精度比HYpl提高0.10%，比MDps与MAps提高4.91%与4.76%，比DBp提高6.01%，虽然该模式可提高总体均值估计精度，但与采用样地−抽样LiDAR数据协同的HYpl相比，提高幅度较小，这主要是由于S-2A遥感特征对蓄积量不够敏感，对应的估测模型精度低导致的。
表 5 5种方法的总体均值估计结果和精度Table 5. Results and accuracy of total mean estimation of the five methods方法
Method总体均值/（m3·hm−2）
Total mean/(m3·ha−1)均值方差/（m3·hm−2）
Mean variance/(m3·ha−1)标准误/（m3·hm−2）
Standard error/(m3·ha−1)估测精度
Estimation accuracyDBp 212.54 107.51 10.37 90.44% MDps 202.09 76.18 8.72 91.54% MAps 202.38 73.55 8.58 91.69% HYpl 210.07 21.42 4.63 96.35% MDpls 208.96 19.95 4.47 96.45% 3.4 总体MSVD估计相对效率
各个方法之间的相对效率结果如图6所示：DBp方法相对于其他4种方法均为负效率（RE均大于1），其中，MDpls方法估计效率最高（REDBp,MDpls = 5.39），HYpl方法次之（REDBp,HYpl = 5.02），DBp相对于MDps与MAps的RE分别为1.41与1.46，可见无论光学数据还是LiDAR数据的加入，估计效率均有所提高。5种方法中，其他4种方法相对于MDpls方法的RE均大于1（REDBp,MDpls = 5.39，REMDps,MDpls = 3.82，REMAps,MDpls = 3.69，REHYpl,MDpls = 1.07），可见MDpls方法效率最高。HYpl相对于MDpls的RE略大于1，说明估计效率低于MDpls，但比较接近；HYpl估计效率远远高于其他3种方法（REDBp,HYpl = 5.02，REMAps,HYpl = 3.43，REMDps,HYpl = 3.56）。包含LiDAR数据的HYpl与MDpls方法相对于包含S-2A数据的MDps与MAps方法均为正效率（REMAps,HYpl = 3.43，REMDps,HYpl = 3.56，REMDps,MDpls = 3.82，REMAps,MDpls = 3.69），显然LiDAR数据比光学数据更能有效提高总体估计效率。MDps与MAps方法之间的RE接近1，但MAps的效率微高于MDps（REMDps,MAps = 1.04）。
4. 讨 论
近年来，激光雷达数据被广泛用于森林资源监测，且精度较高[27,48-49]，但成本高问题限制了该技术的大范围应用，因此部分学者利用LiDAR抽样观测数据，通过天空地协同实现监测目标[11,14-15,50]，这为在有限经费支撑条件下，提高监测效率提供了可行思路。本研究对适用天空地多源数据应用的MDps、MAps、HYpl、MDpls 4种方法以及传统设计推断DBp方法共5种方法进行了林场总体MSVD估计效果比较评价。
比较评价结果表明5种方法的估计精度均在90%以上，其中HYpl与MDpls方法精度达96%以上，可满足森林资源监测的总体精度要求[6,51]。另外，5种方法中，利用了LiDAR数据的HYpl与MDpls 2种方法精度都高于仅基于样地、样地−卫星遥感数据结合的DBp、MDps以及MAps 3种方法，该结果与Saarela等和Puliti等的研究结果一致[7,15,17]。虽然样地蓄积量密度和S-2A卫星遥感特征之间的相关性不高，但在样地的基础上加入卫星遥感数据后，总体估计精度也有所提高（MDps与MAps精度分别提高1.10%与1.25%；REDBp,MDps = 1.41；REDBp,MAps = 1.46）。由于多光谱卫星遥感数据通常都可以免费获取，增加多光谱卫星数据并不会大幅度增加监测成本，值得在实际调查监测业务中推广应用。样地−抽样无人机LiDAR的协同模式对总体MSVD的估计精度达到96.35%，虽比样地−抽样无人机−卫星遥感协同模式的精度略低（REHYpl,MDpls = 1.07），但差别不大，在实际应用中，若只对总体蓄积量密度进行估计，采用样地−无人机LiDAR抽样的协同模式就足够。
需要特别指出的是，这5种方法中有3种方法（DBp、MAps、HYpl）需要以概率抽样获取观测数据。DBp、MAps都要求样地数据是以概率抽取的，而HYpl要求LiDAR数据是以概率抽取的（但并不要求样地数据按概率抽取）。本研究无人机和样地数据的获取就采用了系统抽样观测的思路，但在实际LiDAR和样地数据获取时做了一些方案调整，因此并不严格满足系统抽样条件，这会对这3种方法的应用产生一定负面影响，本研究假设这些不足不会对结论产生影响。但实际应用中应注意样地、无人机数据的获取应严格采用系统抽样方法抽取样本进行观测。
由于本研究获取的样地数量较少，在试验过程中并未考虑分别优势树种类型进行估计，而是将所有类型合并（不分层）进行总体参数估计。实际应用中可以按优势树种类型分层建立模型进行总体参数的估计，理论上可进一步提高总体参数估计精度，但需要对各优势树种类型分别抽取足够多的样本进行样地调查。本研究中DBp方法的方差为抽样误差，并未考虑样地蓄积量密度计算的不确定性；另外4种方法所采用的估测模型都属于线性模型，均值方差主要为模型带来的误差，同样未考虑样地蓄积量密度计算的不确定性，也未考虑由LiDAR估测的蓄积量密度之间的空间自相关带来的不确定性。国内外学者已研发了非线性估测模型和考虑了空间自相关性的分层估测模型和不确定性估计方法[18,30,52]，未来应有针对性地开展相关新技术新方法的应用试验和比较评价研究。
5. 结 论
本研究以旺业甸实验林场为研究区，以林场针叶人工林为总体，针对1种只应用样地数据的模式与3种将样地与遥感数据协同应用的模式（样地−卫星多光谱遥感、样地−无人机LiDAR抽样观测、样地−无人机LiDAR抽样观测−卫星多光谱遥感），对5种林场总体MSVD估计方法进行了比较评价。研究结果表明：（1）和只利用样地数据的估计方法相比，将遥感作为辅助变量建立估测模型，无论是采用对林场全覆盖但对蓄积量不够敏感的S-2A多光谱遥感数据，还是采用抽样式获取的对蓄积量很敏感的LiDAR数据，都可有效提高林场总体MSVD的估测精度。（2）涉及遥感数据应用的4种方法，估计精度最高的为MDpls，其次为HYpl，这2种方法都包含了LiDAR遥感抽样观测数据的应用，都是适合应用于林场总体MSVD估计的年度监测方法。（3）可实现天空地3种观测数据协同应用的方法（MDpls）估测精度和相对效率最高，可作为林场森林蓄积量年度监测的首选方法。
-
![]()
图 2 无人机和样地数据抽样获取方案示意图
Figure 2. Schematic diagram of UAV and sample plot data sampling acquisition scheme
![]()
图 3 各估测方法所利用数据示意图
Figure 3. Schematic diagram of the data used by each estimation method
![]()
图 4 遥感特征选择
折线代表值的变化趋势。The broken line represents the trend of value.
Figure 4. Feature selection of remote sensing
![]()
图 5 遥感特征建模估测模型精度
实线为1∶1验证线。图中RMSE和ME的单位均为m3/hm2。The solid line is 1∶1 verification line.The unit of RMSE and ME in the figure is m3/ha.
Figure 5. Remote sensing feature modeling to estimate model accuracy
表 1 S-2A卫星数据特征
Table 1 S-2A satellite data characteristics
特征名称
Feature name特征符号
Feature symbol计算公式
Calculation formula特征名称
Feature name特征符号
Feature symbol计算公式
Calculation formula光谱特征
Spectral featureB2、B3、B4、B5、B6、B7、B8a、B11、B12 角二阶矩
Angular second momentAN \displaystyle \sum \limits_{i,j = 1}^k {P_{i,j}^2} 差值植被指数
Difference vegetation indexDVI {B_8} - {B_4} 熵
EntropyEN \displaystyle \sum \limits_{i,j = 1}^k {P_{i,j}}\left( { - \ln {P_{i,j}}} \right) 增强植被指数
Enhanced vegetation indexEVI 2.5 × (B8 − B4)/(B8 + 6 × B4 − 7.5 × B2 + 1) 对比度
ContrastCON \displaystyle \sum \limits_{i,j = 1}^k {\left( {i - j} \right)^2} \times {P_{i,j}} 归一化植被指数
Normalized difference vegetation indexNDVI ({B_8} - {B_4})/({B_8} + {B_4}) 均值
MeanME \dfrac{1}{{{k^2}}}\displaystyle \sum \limits_{i,j = 1}^k {P_{i,j}} 比值植被指数
Ratio vegetation indexSR {B_8}/{B_4} 方差
VarianceVAR \displaystyle \sum \limits_{i,j = 1}^k {P_{i,j} }\left( {i,j - {\rm{ME}}} \right) 转换归一化植被指数
Transformed normalized difference vegetation indexTNDVI \sqrt {({B_8} - {B_4})/({B_8} + {B_4}) + 0.5} 同质性
HomogeneityHOM \displaystyle \sum \limits_{i,j = 1}^k \dfrac{ { {P_{i,j} } } }{ {1 + { {\left( {i - j} \right)}^2} } } 叶绿素指数
Chlorophyll indexCIg {B_8}/{B_3} - 1 相关性
CorrelationCOR \displaystyle\sum\limits_{i,j = 1}^k {P_{i,j}^2} \left[ {\frac{{\left( {i - {u_i}} \right) - \left( {j - {u_j}} \right)}}{{\sqrt {\sigma _i^2\sigma _j^2} }}} \right] 反向红边叶绿素指数
Inverted red edge chlorophyll indexIRECI ({B_7} - {B_4})/({B_5}/{B_6}) 相异性
DissimilarityDIS \displaystyle\sum\limits_{i,j = 1}^k {{P_{i,j}}} |i - j| 色素简单比值指数
Pigment specific simple ratioPSSRA {B_7}/{B_4} 修正叶绿素吸收反射指数
Modified chlorophyll absorption in reflectance indexMCARI [(B5 − B4) − 0.2 × (B5 − B3)] × (B5/B4) S-2A 红边位置指数
S-2A red edge positionS2REP 705 + 35 × [(B4 + B7)/2 − B5]/(B6 − B5) 修正窄边红边简单比值指数
Modified simple ratio red edge narrowMSRren [({B_{8{\text{a}}}}/{B_5}) - 1]/\sqrt {({B_{8{\text{a}}}}/{B_5}) + 1} 红边叶绿素指数
Red edge chlorophyll indexCIgre1 {{B}}_{{5}}{/}{{B}}_{{3}}{-1} 红边植被指数
Red edge vegetation indexNDVIre1 {(}{{B}}_{{8}}{-}{{B}}_{{5}}{)/(}{{B}}_{{8}}{ + }{{B}}_{{5}}{)} CIgre2 {B_5}/{B_3} - 1 NDVIre2 {(}{{B}}_{{8}}{-}{{B}}_{{6}}{)/(}{{B}}_{{8}}{ + }{{B}}_{{6}}{)} CIgre3 {{B}}_{{7}}{/}{{B}}_{{3}}{-1} NDVIre3 {(}{{B}}_{{8}}{-}{{B}}_{{7}}{)/(}{{B}}_{{8}}{ + }{{B}}_{{7}}{)} 注： {{B}}_{{2}} 、 {{B}}_{{3}} 、 {{B}}_{{4}} 、 {{B}}_{{5}} 、 {{B}}_{{6}} 、 {{B}}_{{7}} 、B8、{ {B} }_{ {8{\rm{a}}} }、 {{B}}_{{11}} 、 {{B}}_{{12}} 代表S-2A数据对应的波段，{ {P} }_{ {i,j} }{=}{ {D} }_{ {i,j} }{/}\displaystyle{ \sum} _{ {i,j}{=1} }^{ {k} }{ {D} }_{ {i,j} }， {{D}}_{{i,j}} 为i行j列对应的像元值，k代表计算纹理时窗口大小。Notes: {{B}}_{{2}} , {{B}}_{{3}} , {{B}}_{{4}} , {{B}}_{{5}} , {{B}}_{{6}} , {{B}}_{{7}} , B8, { {B} }_{ {8{\rm{a}}} }, {{B}}_{{11}} and {{B}}_{{12}} represent the bands corresponding to S-2A data; { {P} }_{ {i,j} }{=}{ {D} }_{ {i,j} }{/}\displaystyle {\sum} _{ {i,j}{=1} }^{ {k} }{ {D} }_{ {i,j} }, {{D}}_{{i,j}} is the pixel value corresponding to the i row and j column, k represents the window size when calculating the texture. 
下载: 导出CSV
表 2 LiDAR数据特征
Table 2 Feature of LiDAR data
特征名称
Feature name特征符号
Feature symbol均值
MeanHmean, Imean 方差和标准差
Variance and standard deviationHvar, Ivar, Hsd, Isd 最大值和最小值
Maximum value and minimum valueHmax, Imax, Hmin, Imin 变异系数
Coefficient of variationHcv, Icv 四分距差
Interquartile distanceHiq, Iiq 偏斜度
SkewnessHsk, Isk 百分位数
PercentileHp05, Hp10, Hp20, Hp25, ···, Hp95, Hp99
Ip05, Ip10, Ip20, Ip25,···, Ip95, Ip99最小高度以上返回点
Count of return point above the minimum heightR1H, R2H,···, R8H, R9H 注：H代表高度，I代表强度，var代表方差，sd代表标准差，max与min分别代表最大值与最小值，cv、iq以及sk分别代表变异系数、四分距差以及偏斜度，p05—p99代表对应的百分位数。Notes: H represents height, I represents intensity, var represents variance, sd represents standard deviation, max and min represent maximum value and minimum value, respectively; cv, iq and sk represent coefficient of variation, interquartile distance and skewness, respectively; p05−p99 represent the corresponding percentiles.
下载: 导出CSV
表 3 材积计算公式
Table 3 Formula for volume calculation
树种 Tree species 材积计算公式 Formula for volume calculation 白桦
Betula platyphylla{V}{=}{{10}}^{{0.000}\;{397}\;{507}\;{075}\;{980}\;{248 \;\times\; }{{d}}^{{2}}{ \;\times\; }{h}{ \;+\; 3.812}\;{138}\;{080}\;{735}\;{6 \;\times\; }{{10}}^{-{6}}{ \;\times\; }{{d}}^{{3}}{ \;\times\; }{h}\;-\;{0.000}\;{843}\;{446}\;{660}\;{194}\;{932 \;\times\; }{{d}}^{{2}}\;-\;{0.000}\;{290}\;{088}\;{083}\;{727}\;{618 \;\times\; }{{d}}^{{2}}{ \;\times\; }{h}{ \;\times\; }{\lg}\;{d}} 黑桦
Betula davurica{V}{=0.019}\;{921}\;{519}\;{389}\;{050}\;{9 \;+\; 0.000}\;{388}\;{145}\;{167}\;{080}\;{27 \;\times\; }{{d}}^{{2}}\;-\;{2.050}\;{596}\;{607}\;{769}\;{77 \;\times\; }{{10}}^{-{5}}{ \;\times\; }{{d}}^{{2}}{ \;\times\; }{h} -
{0.000}\;{870}\;{310}\;{875}\;{746}\;{131 \;\times\; }{{h}}^{{2}}{ \;+\; 8.053}\;{621}\;{369}\;{495}\;{43 \;\times\; }{{10}}^{-{5}}{ \;\times\; }{d}{ \;\times\; }{{h}}^{{2}}落叶松
Larix gmelinii{V}{=}{{10}}^{-{3.498}\;{903}\;{907}\;{280}\;{04 \;+\; 2.755}\;{048}\;{465}\;{025}\;{64 \;\times\; }{\lg}\;{d}\;-\;{0.394}\;{050}\;{839}\;{410}\;{844 \;\times\; }{\lg}\;{{d}}^{{2}}\;-\;{1.379}\;{153}\;{564}\;{042}\;{94 \;\times\; }{\lg}\;{h}{ \;+\; 1.145}\;{866}\;{814}\;{777}\;{51 \;\times\; }{\lg}\;{{h}}^{{2}}} 樟子松
Pinus sylvestris{V}{=0.000}\;{445}\;{504}\;{541}\;{007}\;{861 \;\times\; }{{d}}^{{1.663}\;{165}\;{237}\;{864}\;{29 \;\times\; \exp(0.098}\;{992}\;{935}\;{700}\;{498}\;{1 \;\times\; }{h}\;-\;{1.666}\;{816}\;{460}\;{562}\;{17/}{h}{)}} 油松
Pinus tabuliformis{V}{=2.533}\;{092}\;{907}\;{631}\;{68 \;\times\; }{ {10} }^{ {-5} }{ \;\times\; }{ {d} }^{ {2} }{ \;\times\; }{h}\;-\;{8.573}\;{788}\;{411}\;{603}\;{25 \;\times\; }{ {10} }^{ {-7} }{ \;\times\; }{ {d} }^{ {3} }{ \;\times\; }{h}\;-\;{6.000}\;{334}\;{292}\;{010}\;{54 \;\times\; }{ {10} }^{ {-5} }{ \;\times\; }{ {d} }^{ {2} } \;+\;
{ 2.937}\;{206}\;{772}\;{188}\;{84 \;\times\; }{ {10} }^{ {-5} }{ \;\times\; }{ {d} }^{ {2} }{ \;\times\; }{h}{ \;\times\; }{ {\lg} }\; {d}注：引自文献[40]。V代表单木材积；d代表单木胸径；h代表单木树高。Notes: cited from reference [40]. V represents single tree stock volume, d represents single tree DBH, h represents single tree height.
下载: 导出CSV
表 4 遥感特征优选后保留的特征
Table 4 Remotes sensing features selected after optimum feature selection
数据源
Data source特征选择结果
Result of feature selectionS-2A B2，MCARI，B2ME，B12COR LiDAR Hmean，Ip99，R2H
下载: 导出CSV
表 5 5种方法的总体均值估计结果和精度
Table 5 Results and accuracy of total mean estimation of the five methods
方法
Method总体均值/（m3·hm−2）
Total mean/(m3·ha−1)均值方差/（m3·hm−2）
Mean variance/(m3·ha−1)标准误/（m3·hm−2）
Standard error/(m3·ha−1)估测精度
Estimation accuracyDBp 212.54 107.51 10.37 90.44% MDps 202.09 76.18 8.72 91.54% MAps 202.38 73.55 8.58 91.69% HYpl 210.07 21.42 4.63 96.35% MDpls 208.96 19.95 4.47 96.45%
下载: 导出CSV
-
[1] 宋进春. 省级森林资源面积年度出数研究:以广西壮族自治区为例[J]. 中南林业调查规划, 2021, 40(2): 9−14, 26. Song J C. Research on annual results for provincial forest resource area: a case study in Guangxi Zhuang Autonomous Region[J]. Central South Forest Inventory and Planning, 2021, 40(2): 9−14, 26.
[2] 张晓丽, 游先祥. 应用“3S”技术进行北京市森林立地分类和立地质量评价的研究[J]. 遥感学报, 1998, 2(4): 292−295. doi: 10.11834/jrs.19980411 Zhang X L, You X X. Application of 3S technology to forest site type classification and site quality evaluation in Beijing[J]. Journal of Remote Sensing, 1998, 2(4): 292−295. doi: 10.11834/jrs.19980411
[3] 黄国胜, 马炜, 王雪军, 等. 基于一类清查数据的福建省立地质量评价技术[J]. 北京林业大学学报, 2014, 36(3): 1−8. Huang G S, Ma W, Wang X J, et al. Forestland site quality evaluation of Fujian Province based on continuous forest inventory data[J]. Journal of Beijing Forestry University, 2014, 36(3): 1−8.
[4] 王雪军, 张煜星, 黄国胜, 等. 全国森林面积和森林蓄积年度出数方法探讨[J]. 江西农业大学学报, 2016, 38(1): 9−18. Wang X J, Zhang Y X, Huang G S, et al. Discussion on methods for annual national producing estimates of forest area and forest stock in China[J]. Acta Agriculturae Universitatis Jiangxiensis, 2016, 38(1): 9−18.
[5] 聂祥永, 姚顺彬, 张敏, 等. 国家森林资源连续清查数据处理统计规范(LY/T 1957—2011)[S]. 北京: 中国标准出版社, 2011. Nie X Y , Yao S B, Zhang M, et al. Specification for data processing and statistic in national forest inventory(LY/T 1957−2011)[S]. Beijing: China Standard Press, 2011.
[6] 唐小平, 翁国庆, 陈雪峰, 等. 森林资源规划设计调查技术规程(GB/T26424—2010)[S]. 北京: 中国标准出版社, 2010. Tang X P, Weng G Q, Chen X F, et al. Technical regulations for inventory for forest management planning and design (GB/T26424−2010)[S]. Beijing: China Standard Press, 2010.
[7] Saarela S, Holm S, Grafström A, et al. Hierarchical model-based inference for forest inventory using three sources of information[J]. Annals of Forest Science, 2016, 73: 895−910. doi: 10.1007/s13595-016-0590-1
[8] Ståhl G, Saarela S, Schnell S, et al. Use of models in large-area forest surveys: comparing model-assisted, model-based and hybrid estimation[J]. Forest Ecosystems, 2016, 3: 1−11. doi: 10.1186/s40663-016-0060-0
[9] Ene L T, Næsset E, Gobakken T, et al. Assessing the accuracy of regional LiDAR-based biomass estimation using a simulation approach[J]. Remote Sensing of Environment, 2012, 123: 579−592. doi: 10.1016/j.rse.2012.04.017
[10] Graubard B I, Korn E L. Inference for superpopulation parameters using sample surveys[J]. Statistical Science, 2002, 17(1): 73−96.
[11] Ståhl G, Holm S, Gregoire T G, et al. Model-based inference for biomass estimation in a LiDAR sample survey in Hedmark County, Norway[J]. Canadian Journal of Forest Research, 2011, 41(1): 96−107. doi: 10.1139/X10-161
[12] Corona P, Fattorini L, Franceschi S, et al. Estimation of standing wood volume in forest compartments by exploiting airborne laser scanning information: model-based, design-based, and hybrid perspectives[J]. Canadian Journal of Forest Research, 2014, 44(11): 1303−1311. doi: 10.1139/cjfr-2014-0203
[13] McRoberts R E, Westfall J A. The effects of uncertainty in model predictions of individual tree volume on large area volume estimates[J]. Forest Science, 2014, 60: 34−43. doi: 10.5849/forsci.12-141
[14] Puliti S, Ene L T, Gobakken T, et al. Use of partial coverage UAV data in sampling for large scale forest inventories[J]. Remote Sensing of Environment, 2017, 194: 115−126. doi: 10.1016/j.rse.2017.03.019
[15] Puliti S, Saarela S, Gobakken T, et al. Combining UAV and Sentinel-2 auxiliary data for forest growing stock volume estimation through hierarchical model-based inference[J]. Remote Sensing of Environment, 2018, 204: 485−497. doi: 10.1016/j.rse.2017.10.007
[16] Condés S, McRoberts R E. Updating national forest inventory estimates of growing stock volume using hybrid inference[J]. Forest Ecology and Management, 2017, 400: 48−57. doi: 10.1016/j.foreco.2017.04.046
[17] Saarela S, Holm S, Healey S P, et al. Generalized hierarchical model-based estimation for aboveground biomass assessment using GEDI and landsat data[J]. Remote Sensing, 2018, 10: 1832. doi: 10.3390/rs10111832
[18] Saarela S, Wästlund A, Holmström E, et al. Mapping aboveground biomass and its prediction uncertainty using LiDAR and field data, accounting for tree-level allometric and LiDAR model errors[J]. Forest Ecosystems, 2020, 7(1): 1−17. doi: 10.1186/s40663-019-0212-0
[19] Richard W Guldin. A systematic review of small domain estimation research in forestry during the twenty-first century from outside the United States[J/OL]. Frontiers in Forests and Global Change, 2021, 4: 695929[2022−01−05]. https://doi.org/10.3389/ffgc.2021.695929.
[20] 肖兴威. 中国森林资源清查[M]. 北京: 中国林业出版社, 2005. Xiao X W. National forest inventory of China[M]. Beijing: China Forestry Publishing House, 2005.
[21] 曾伟生, 骆期邦, 彭长清. 两阶群团抽样在森林调查中的估计效率研究[J]. 林业科学研究, 1995, 8(5): 483−488. doi: 10.13275/j.cnki.lykxyj.1995.05.003 Zeng W S, Luo Q B, Peng C Q. Study on efficiency of two-stage cluster sampling in forest inventory[J]. Forest Research, 1995, 8(5): 483−488. doi: 10.13275/j.cnki.lykxyj.1995.05.003
[22] 李春干, 陈琦, 谭必增. 基于卫星遥感数据空中抽样的大尺度森林资源动态监测[J]. 林业资源管理, 2009, 2: 106−110,127. doi: 10.3969/j.issn.1002-6622.2009.02.020 Li C G, Chen Q, Tan B Z. Large-scale forest resources monitoring by means of spatial sampling of Hi-resolution remote sensing images[J]. Forest Resources Management, 2009, 2: 106−110,127. doi: 10.3969/j.issn.1002-6622.2009.02.020
[23] 郑冬梅, 智长贵, 黄国胜, 等. 遥感大样地点面判读方法在森林资源宏观监测中的应用分析[J]. 林业资源管理, 2017, 6(6): 137−142. Zheng D M, Zhi C G, Huang G S, et al. Analysis of forest resources monitoring by zoning interpretation and point interpretation of remote sensing large plot[J]. Forest Resources Management, 2017, 6(6): 137−142.
[24] 唐守正. 关于两相抽样面积蓄积统计的原则[J]. 林业资源管理, 1996, 4: 18−22. Tang S Z. The principles of two-phase sampling area accumulation statistics[J]. Forest Resources Management, 1996, 4: 18−22.
[25] 张宗秀, 高天雷, 张文. 双重二阶抽样提高森林资源抽样精度的研究[J]. 四川林业科技, 2013, 34(5): 8−12. Zhang Z X, Gao T L, Zhang W. Study on improving sampling accuracy of forest resources by double two-stage sampling[J]. Journal of Sichuan Forestry Science and Technology, 2013, 34(5): 8−12.
[26] 李增元, 赵磊, 李堃, 等. 合成孔径雷达森林资源监测技术研究综述[J]. 南京信息工程大学学报(自然科学版), 2020, 12(2): 150−158. Li Z Y, Zhao L, Li K, et al. A survey of developments on forest resources monitoring technology of synthetic aperture radar[J]. Journal of Nanjing University of Information Science and Technology (Natural Science Edition), 2020, 12(2): 150−158.
[27] 庞勇, 李增元, 陈博伟, 等. 星载激光雷达森林探测进展及趋势[J]. 上海航天, 2019, 36(3): 20−28. Pang Y, Li Z Y, Chen B W, et al. Status and development of spaceborne lidar applications in forestry[J]. Aerospace Shanghai, 2019, 36(3): 20−28.
[28] 黄华国. 林业定量遥感研究进展和展望[J]. 北京林业大学学报, 2019, 41(12): 1−14. Huang H G. Progress and perspective of quantitative remote sensing of forestry[J]. Journal of Beijing Forestry University, 2019, 41(12): 1−14.
[29] 张王菲, 陈尔学, 李增元, 等. 干涉、极化干涉SAR技术森林高度估测算法研究进展[J]. 遥感技术与应用, 2017, 32(6): 983−997. Zhang W F, Chen E X, Li Z Y, et al. Development of forest height estimation using InSAR/PolInSAR technology[J]. Remote Sensing Technology and Aplication, 2017, 32(6): 983−997.
[30] Zhao J P, Zhao L, Chen E X, et al. An improved generalized hierarchical estimation framework with geo statistics for mapping forest parameters and its uncertainty: a case study of forest canopy height[J/OL]. Remote Sensing, 2022, 14: 568[2022−10−22]. https://doi.org/10.3390/rs14030568.
[31] 宋新民. 抽样调查技术(第2版)[M]. 北京: 中国林业出版社, 2007. Song X M. Sampling survey techniques[M]. Beijing: China Forestry Publishing House, 2007.
[32] 赵凯峰. 双重回归估计法测算森林生长量的探讨[J]. 林业资源管理, 1996, 5: 48−52. Zhao K F. Study on estimating forest growth by double regression estimation method[J]. Forest Resources Management, 1996, 5: 48−52.
[33] 鲁赛尼·阿特马维扎扎, 劳可遒. 用航空照片双重抽样进行森林资源清查[J]. 中南林业调查规划, 1982, 1: 53−55, 43. Lusaini A, Lao K Q. Double sampling using aerial photo for forest resource inventory[J]. Central South Forest Inventory and Planning, 1982, 1: 53−55, 43.
[34] 张秋江, 袁亚东, 丛修宝. 用遥感及实测回归估计森林蓄积量[J]. 东北林业大学学报, 1997, 25(2): 18−23. Zhang Q J, Yuan Y D, Cong X B. The method for estimating forest volume by remote sensing and regression of practical measurement[J]. Journal of Northeast Forestry University, 1997, 25(2): 18−23.
[35] 陈振雄, 熊智平, 曾伟生, 等. 基于大样地双重抽样方法的广东省森林资源监测研究[J]. 中南林业调查规划, 2014, 33(3): 28−33. doi: 10.3969/j.issn.1003-6075.2014.03.009 Chen Z X, Xiong Z P, Zeng W S, et al. Forest resources monitoring based on double sampling with large plot in Guangdong[J]. Central South Forest Inventory and Planning, 2014, 33(3): 28−33. doi: 10.3969/j.issn.1003-6075.2014.03.009
[36] 陈中超, 刘清旺, 李春干, 等. 基于无人机激光雷达的人工林碳储量线性与非线性估测模型比较[J]. 北京林业大学学报, 2021, 43(12): 9−16. doi: 10.12171/j.1000-1522.20200417 Chen Z C, Liu Q W, Li C G, et al. Comparison in linear and nonlinear estimation models of carbon storage of plantations based on UAV LiDAR[J]. Journal of Beijing Forestry University, 2021, 43(12): 9−16. doi: 10.12171/j.1000-1522.20200417
[37] Donoghue D N M, Watt P J, Cox N J, et al. Remote sensing of species mixtures in conifer plantations using LiDAR height and intensity data[J]. Remote Sensing of Environment, 2007, 110(4): 509−522. doi: 10.1016/j.rse.2007.02.032
[38] McGaughey R J. FUSION/LDV: Software for LiDAR data analysis and visualization, January 2021-FUSION Version 4.20[Z/OL]. Washington D C: United Stated of Department of Agriculture, 2021. http://forsys.cfr.washington.edu/software/fusion/FUSION_manual.pdf.
[39] 曹林, 徐婷, 申鑫, 等. 集成Landsat OLI和机载LiDAR条带数据的亚热带森林生物量制图[J]. 遥感学报, 2016, 20(4): 665−678. doi: 10.11834/jrs.20165091 Cao L, Xu T, Shen X, et al. Mapping biomass by integrating Landsat OLI and airborne LiDAR transect data in subtropical forests[J]. Journal of Remote Sensing, 2016, 20(4): 665−678. doi: 10.11834/jrs.20165091
[40] 何诚. 森林精准计测关键技术研究[D]. 北京: 北京林业大学, 2013. He C. The key technology fore precision measurement in forest surveying[D]. Beijing: Beijing Forestry University, 2013.
[41] Davidson R, MacKinnon J G. Estimation and Inference in Econometrics[M]. New York: Oxford University Press, 1993.
[42] McConville K S, Moisen G G, Frescino T S. A tutorial on model-assisted estimation with application to forest inventory[J]. Forests, 2020, 11: 244. doi: 10.3390/f11020244
[43] Ståhl G, Heikkinen J, Petersson H, et al. Sample-based estimation of greenhouse gas emissions from forests-a new approach to account for both sampling and model errors[J]. Forest Science, 2014, 60: 3−13. doi: 10.5849/forsci.13-005
[44] Genuer R, Poggi J M, Tuleau-Malot C. VSURF: an R package for variable selection using random forests[J]. R Journal, 2015, 7(2): 19−33. doi: 10.32614/RJ-2015-018
[45] Næsset E, Ørka H O, Solberg S, et al. Mapping and estimating forest area and aboveground biomass in Miombo woodlands in Tanzania using data from airborne laser scanning, tandem-X, Rapideye, and global forest maps: a comparison of estimated precision[J]. Remote Sensing of Environment, 2016, 175: 282−300. doi: 10.1016/j.rse.2016.01.006
[46] Lu D, Chen Q, Wang G, et al. Aboveground forest biomass estimation with Landsat and LiDAR data and uncertainty analysis of the estimates[J]. International Journal of Forestry Research, 2012, 2012: 1−16.
[47] Chen L, Wang Y Q, Ren C Y, et al. Assessment of multiwavelength SAR and multispectral instrument data for forest aboveground biomass mapping using random forest kriging[J]. Forest Ecology and Management, 2019, 447: 12−25. doi: 10.1016/j.foreco.2019.05.057
[48] 赵峰, 李增元, 王韵晟, 等. 机载激光雷达(LiDAR)数据在森林资源调查中的应用综述[J]. 遥感信息, 2008, 1: 106−110, 53. doi: 10.3969/j.issn.1000-3177.2008.01.021 Zhao F, Li Z Y, Wang Y S, et al. The application of LiDAR data in forest[J]. Remote Sensing Information, 2008, 1: 106−110, 53. doi: 10.3969/j.issn.1000-3177.2008.01.021
[49] 韩婷婷. 激光雷达数据在森林垂直结构参数反演中的应用综述[J]. 北京测绘, 2020, 34(8): 1061−1065. Han T T. Summarization of the application of LiDAR data in forest vertical structure parameters inversion[J]. Beijing Surveying and Maping, 2020, 34(8): 1061−1065.
[50] 赵俊鹏, 赵磊, 陈尔学, 等. ZY3立体像对和机载LiDAR抽样数据协同估测森林平均高[J]. 林业科学, 2020, 57(9): 66−75. Zhao J P, Zhao L, Chen E X, et al. Synergy application of ZY3 stereo imaging pairs and airborne LiDAR sampling data for estimating mean height of forest[J]. Scientia Silvae Sinicae, 2020, 57(9): 66−75.
[51] 张煜星, 闫宏伟, 黄国胜, 等. 森林资源连续清查技术规程(GB/T 38590—2020)[S]. 北京: 中国标准出版社, 2020. Zhang Y X, Yan H W, Huang G S, et al. Technical regulations for continuous forest inventory(GB/T 38590−2020)[S]. Beijing: China Standard Press, 2020.
[52] Chen Q, Mcroberts R E, Wang C W, et al. Forest aboveground biomass mapping and estimation across multiple spatial scales using model-based inference[J]. Remote Sensing of Environment, 2016, 184: 350−360. doi: 10.1016/j.rse.2016.07.023
-
期刊类型引用(4)
1. 高郯，张铎，卢杰，王超，李江荣. 色季拉山高山松林降雨再分配及重金属元素的时空特征研究. 西南林业大学学报(自然科学). 2022(01): 115-123 . 
百度学术
2. 陈新宇，孟景祥，周先清，袁虎威，钮世辉，李悦. 油松地理种群针叶形态解剖与生理指标遗传变异分析. 北京林业大学学报. 2019(07): 19-30 . 本站查看
3. 金微微，张会慧，滕志远，孙广玉，许楠. 乡土风箱果和紫叶风箱果及其杂交种F_1叶片的光合功能研究. 中南林业科技大学学报. 2018(04): 33-39 . 百度学术
4. 周丽，陈诗，朱存福，许玉兰，李悦，李伟，蔡年辉. 高山松不同种源苗木在云南松生境下的生物量与生长分析. 广东农业科学. 2017(02): 68-75 . 百度学术
其他类型引用(4)
计量
- 文章访问数: 719
- HTML全文浏览量: 208
- PDF下载量: 89
- 被引次数: 8
