基于遥感-样地数据同化驱动的森林蓄积量估算及校正

Estimation and Calibration of Forest Stock Volume Based on Remote Sensing-Plot Data Assimilation

摘要

森林蓄积量作为衡量森林资源数量与质量的重要指标，其年度更新在森林资源管理中发挥着重要作用。然而，目前借助遥感数据完成森林蓄积量的年度更新存在较大的低值高估和高值低估问题。为此，以黑河市森林为研究对象，使用2010年森林资源一类清查固定样地数据与同期Landsat-7多光谱遥感影像数据，采用似乎不相关回归模型（seemingly unrelated regression，SUR）和最佳线性无偏预测器（best linear unbiased prediction，BLUP），构建“遥感-样地”协同优化的森林蓄积量估算体系。结果表明，“遥感-样地”协同优化后的森林蓄积量估算精度显著提高，同化校正后的模型调整决定系数（R²）从0.200提升至0.539，模型的拟合优度明显增加；均方根误差（RMSE）由23.207降低至17.541 m³/hm²，改进百分比达24.4%，表明蓄积量的预测值与真实值的绝对偏差程度显著降低；蓄积量估计值的标准差由17.367降低至12.146 m³/hm²，进一步证实优化后的模型参数估计的稳定性增强。研究表明所构建的数据同化模型显著提升了森林蓄积量遥感估算的精度和可靠性，为区域尺度森林碳汇评估和森林可持续经营决策提供可靠的技术支撑，同时支持了地级市尺度森林年度更新工作。

Abstract

Forest stock volume， a key indicator for measuring the quantity and quality of forest resources， plays a pivotal role in forest resource management through its annual updates. However， there remains a significant issue of overestimation for low values and underestimation for high values in the current annual update of forest stock volume using remote sensing data. Taking the forests in Heihe City as the research object， this study constructed a ‘remote sensing-plot’ synergistic optimization system for forest stock volume estimation by integrating seemingly unrelated regression （SUR） model and best linear unbiased prediction （BLUP）， based on 2010 national forest inventory permanent sample plot data and concurrent Landsat-7 multispectral remote sensing image data. The results indicated that the estimation accuracy of forest stock volume after ‘remote sensing-plot’ synergistic optimization was significantly improved， the coefficient of determination （R²） of the assimilated and calibrated model increased from 0.200 to 0.539， manifesting a marked enhancement in model goodness-of-fit； the root mean square error （RMSE） decreased from 23.207 to 17.541m³/hm²， with an improvement percentage as high as 24.4%， which signified a substantial reduction in the absolute deviation between predicted and actual values. The standard variance for volume estimates dropped from 17.367 to 12.146m³/hm²， further confirming the enhanced stability of parameter estimation in the optimized model. It is evident that the data assimilation model constructed in this study significantly improves the accuracy and reliability of remote sensing-based forest stock volume estimation， providing robust technical support for regional-scale forest carbon sink assessment and sustainable forest management decision-making. It also supports the annual forest update work at the city-level scale.

Graphical abstract

关键词

数据同化 / 似乎不相关回归模型 / 经验最佳线性无偏预测器 / 蓄积量 / 森林资源年度更新

Key words

Data assimilation / seemingly unrelated regression （SUR） model / best linear unbiased prediction （BLUP） / stock volume / annual updates for forest resources

0 引言

森林生态系统蓄积量作为评估森林结构、生物量以及碳汇潜力的重要参数，其年度变化的动态监测不仅是自然资源管理的技术难题，也逐渐成为国内外研究的前沿热点^［1］。作为陆地生态系统中最重要的碳库，森林在固碳释氧、减缓全球变暖和吸收温室气体等方面具有至关重要的作用^［2］。传统的森林蓄积量调查方法难以满足当今森林资源监测与管理高效与准确的要求^［3］。在衡量森林资源数量与质量时，森林蓄积量通常作为重要的指标，其年度更新在森林资源管理中更是发挥着重要的作用。目前我国森林资源采用5年一个周期、每年完成1/5的抽样调查方案，并通过移动平均法对森林资源总体情况进行联合估计^［4］，但该方法存在估计偏差大、方差大等问题。因此，结合森林资源清查数据与遥感技术，根据统计推断理论，实现森林资源高精度估计是目前研究的热点。随着3S技术和无人机摄影测量的快速发展，遥感影像与地面调查数据融合建模已成为蓄积量估测的重要手段^［5］。大尺度森林资源监测主要依靠优化传统抽样方法与融合遥感辅助数据2种方法^［6-8］来实现，如通过采用适当降低样地密度并结合滑动平均法，实现了森林资源的年度估算；同时，还有许多研究人员选择通过融合多源数据的方式，以达到提升效率的目的^［9-10］。在遥感辅助估测方面，已经形成了较为成熟的技术体系，包括设计推断法、模型辅助法、模型依赖法及混合估计法等^［11］，其中估测模型通过建立森林参数与遥感数据之间的多种数学关系来实现高精度预测^［12］。在设计推断法方面我国虽然进展良好，但在模型辅助法方面的研究与实际应用较薄弱；在模型依赖法上，虽然关于模型构建的研究丰富，但在对模型精度的评估和不确定性分析方面依旧存在着许多不足；混合估计法尚未引入数据同化等先进技术，仅对部分方法有初步探索，所以研究较为滞后；模型辅助法在国际上应用广泛，整合了地面调查数据与遥感数据，结合数据同化技术构建估计模型并有效地提升森林参数估测精度，对我国森林资源连续清查体系的优化具有重要参考价值^［11］。

数据同化技术在融合了实地观测数据与过程模型后，充分考虑了数据的时空异质性和误差特征，生成了具有物理一致性的数据集^［13］。这种方法不仅显著提高了估算精度，还可以反演出许多难以直接观测的生态参数，为森林资源管理提供技术路径。早期数据同化技术的研究主要借鉴了大气和海洋等科学算法的经验算法和统计推断方法。集合卡尔曼滤波（ensemble Kalman filter，EnKF）算法是当时最具突破性意义的方法^［14］，该算法通过引入蒙特卡罗模型的思想，克服了传统卡尔曼滤波在处理非线性问题时的局限性，但易受高斯分布假设约束，计算的复杂程度也在处理高维数据时明显上升；近年来发展的粒子滤波（particle filter，PF）技术则突破了这一限制，通过按照一定的顺序或者逻辑连贯性进行蒙特卡罗采样来近似后验概率分布，在处理高度非线性和非高斯系统时粒子滤波（PF）技术具有更大的优势^［15］。相较于集合卡尔曼滤波算法，粒子滤波算法避免了复杂的矩阵运算，计算效率更高，特别适用于森林生态系统这类具有明显时空异质性的研究对象^［13］。这些数据同化方法对于长时间序列参数的优化有明显的优势^［16］，但针对森林资源清查长时间序列数据缺失的问题，似乎不相关模型结合最佳线性无偏估计的方法更为适合。

本研究以黑龙江省黑河市森林为对象，基于2010年森林资源一类清查样地数据，采用多源数据协同分析框架，整合森林资源清查数据与Landsat-7多光谱遥感影像数据，运用似乎不相关回归（seemingly unrelated regression，SUR）模型进行联合估计，并采用最佳线性无偏预测（best linear unbiased prediction，BLUP）方法进行模型校准，探究数据同化（data assimilation，DA）技术在区域尺度森林蓄积量遥感估算中的优化作用，旨在提升传统森林蓄积量遥感估算的精度和可靠性，有效改善森林蓄积量遥感估算普遍存在的低值高估和高值低估问题，为区域尺度上森林年度更新估计问题提供了有效技术路径。

1 研究区概况与数据

1.1 研究区概况

黑河市位于黑龙江省西北部，地处小兴安岭北麓，地理坐标47°42′~51°03′N，124°45′~129°18′E；地势整体西北高、东南低，以中低山丘陵为主（64.3%），平原占35%。气候属寒温带大陆性季风气候，冬季寒冷漫长，最低气温达-48.1 ℃，夏季短促温暖，最高气温38.2℃，年降水量500~550 mm。土壤以暗棕壤（65.76%）和黑土（12.86%）为主，沼泽土、草甸土等非地带性土壤广泛分布。森林覆盖率达46.7%，主要树种为红松（Pinus koraiensis）、兴安落叶松（Larix gmelinii）和白桦（Betula platyphylla），另有16.08万hm²草原和85.77万hm²湿地，生态功能显著。

1.2 数据收集及预处理

1.2.1 森林资源清查数据

森林资源清查数据选用了黑河市2010年森林资源一类清查样地数据，样地大小为0.06 hm²，共设置493个样地，覆盖研究区域内的主要森林类型和地形条件，分别测量样地经度、纬度、地貌、海拔、坡度、坡向、植被类型、优势树种、平均树高、平均胸径和郁闭度等信息。在每个样地内，对所有胸径（diameter at breast height，DBH，式中记为D_BH）大于或等于5 cm的数据进行每木检尺。根据样地每木检尺数据，通过材积方程计算不同树种的材积，对于常见树种，材积计算公式如下。

V=a\times D_{\mathrm{B}\mathrm{H}}^b。（1）

式中：V为单株材积，m³，D_BH为胸径，cm；a，b为不同树种特定参数，具体参数值见表1^［17-18］。

将样地内所有树木的蓄积量相加，得到样地的总蓄积量，再除以样地面积，得到每块样地的单位面积蓄积量（m³/hm²），见表2。

1.2.2 遥感数据

本研究选用的Landsat-7多光谱数据是由美国地质调查局（united states geological survey，USGS）和美国国家航空航天局（national aeronautics and space administration，NASA）联合运营，于1999年4月15日发射，其搭载的增强型专题制图仪（enhanced thematic mapper，ETM+）是对早期Landsat-5 TM传感器的改进，具有优越的数据质量和良好的数据连续性。但需要注意的是由于Landsat-7卫星在2003年5月31日发生了故障，扫描线校正器（scan line corrector，SLC）永久失效，出现了一景影像中约有四分之一数据丢失的现象^［19］。所以在进行Landsat-7多光谱数据预处理时需要进一步进行条带修复处理以确保影像的一致性、完整性和数据质量。本研究下载了2010年7月份研究区域Landsat-7多光谱数据，并对其进行辐射定标、大气校正、几何校正、条带修复以及镶嵌裁剪等预处理。

1.2.3 数字高程模型（digital elevation model，DEM）数据

DEM数据来自中国科学院地理空间数据云平台（http：//www.gscloud.cn/）覆盖的黑龙江省黑河市17幅30 m×30 m分辨率DEM数据。通过提取坡度、坡向和地形起伏度等关键地形因子，为研究提供了重要的地形特征。

2 研究方法

数据同化校准程序需要3个技术步骤。步骤1为使用似不相关回归（SUR）模型构建森林蓄积量预测主模型和子模型，构成模型预测系统；步骤2为基于BLUP校准程序校准SUR模型预测结果；步骤3，从基于模型的推理角度，在步骤1和步骤2中分别获得的预测和校准预测结果的基础上估计总体均值估计器的方差。

2.1 特征因子筛选

研究采用Landsat-7的多光谱影像6个原始波段反射率（B1，B2，B3，B4，B5，B7）；6种常用植被指数，包括简单比值植被指数（ratio vegetation index，RVI）、归一化差值植被指数（normalized difference vegetation index，NDVI）、差值植被指数（difference vegetation index，DVI）、增强型植被指数（enhanced vegetation index，EVI）、大气抵抗植被指数（atmospherically resistant vegetation index，ARVI）、土壤调节植被指数（soil adjusted vegetation index，SAVI）^［5］；地形因子包括海拔（altitude）、坡度（slope）和坡向（aspect）；Landsat-7多光谱影像主成分分析后第一、第二主成分提取8种纹理特征^［20］，包括均值（mean，ME），方差（variance，VA），协同性（homogeneity，HO），熵（entropy，EN），对比度（contrast，CO），二阶矩（second moment，SM），相异性（dissimilarity，DI）和相关性（correlation，CC），滑动窗口大小设置为3×3，共计23个特征因子。筛选出与森林蓄积量存在显著相关关系，即P＜0.05的特征因子作为森林蓄积量遥感估算的候选变量，见图1。

为方便特征因子的后续使用，消除量纲和过大数值对结果造成的影响，需要将不同量纲和数量级的特征因子进行归一化处理，反射率、植被指数和海拔进行归一化处理见公式（2），对坡度进行归一化处理见公式（3），对坡向进行归一化处理见公式（4）。

x_i=(x-x_{\mathrm{m}\mathrm{i}\mathrm{n}})/(x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}})。（2）

x_{\alpha }=\mathrm{s}\mathrm{i}\mathrm{n}\alpha。（3）

x_{\beta }=\frac{\mathrm{c}\mathrm{o}\mathrm{s}\beta +\mathrm{1}}{\mathrm{2}}。（4）

式中：x_i 为归一化后的特征因子值；x为原始特征因子值；x_min为该特征在数据集中的最小值；x_max为该特征在数据集中的最大值；x_{\alpha }为归一化后的坡度值；\alpha为原始坡度角度；x_{\beta }为归一化后的坡向值；\beta为原始坡向角度。

2.2 似乎不相关模型构建

2.2.1 多元线性回归主模型构建

本研究以Pearson相关性分析筛选出的特征因子为自变量，森林蓄积量数据为因变量，构建基于遥感数据的森林蓄积量多元线性回归模型，见公式（5）^［21］。

Y={\beta }_{\mathrm{0}}+{\beta }_{\mathrm{1}}{X}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{X}_{\mathrm{2}}+\cdots +{\beta }_k{X}_k+\epsilon。（5）

式中：Y为因变量，即森林蓄积量；X_{\mathrm{1}}，X_{\mathrm{2}}，…，X_k为自变量，即通过遥感数据获取的特征因子；{\beta }_{\mathrm{0}}为截距项；{\beta }_{\mathrm{1}}，{\beta }_{\mathrm{2}}，…，{\beta }_k为回归系数，量化了每个自变量对因变量的影响程度；\epsilon为误差项，反映了真实值和预测值之间的差异，也是筛选支持性子模型的重要依据。

2.2.2 支持性子模型构建

为预测与主模型残差高度相关的辅助变量，本研究从大量遥感或地形数据中（如NDVI、坡度和坡向等）提取潜在变量，计算每个潜在变量与主模型残差的相关系数，选择相关系数绝对值最大的变量作为支持性子模型的因变量，采用截距模型对支持性子模型进行拟合，见公式（6）。

z=\beta +\epsilon。（6）

式中：z为候选变量；\beta为变量的期望值；\epsilon为误差项。

支持性子模型的残差与主模型残差的协方差会被记录在矩阵的非对角元素中，为后续动态校准提供信息，利用观测数据校准主模型预测结果。

2.2.3 似乎不相关回归（SUR）模型联合估计

当目标是描述一个完整的变量系统时，多个模型是必要的，同时系统的多个因变量不一定是相互独立的^［22］。各个模型预测之间可以通过其随机残差相互作用，SUR可用于构建该模型系统，见公式（7）。该模型系统通过预测残差的联合方差协方差来捕获这些相关性^［23］。矩阵的非对角协方差元素可以提供关于过程的附加信息，这有助于使用BLUP进行校准程序^［22］。

\left[\begin{array}{c}{\mathit{y}}_{\mathrm{1}}\\ {\mathit{y}}_{\mathrm{2}}\\ ⋮\\ {\mathit{y}}_K\end{array}\right]=\left[\begin{array}{cccc}{\mathit{X}}_{\mathrm{1}}& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& {\mathit{X}}_{\mathrm{2}}& \cdots & \mathrm{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & {\mathit{X}}_K\end{array}\right]\left[\begin{array}{c}{\mathit{\beta }}_{\mathrm{1}}\\ {\mathit{\beta }}_{\mathrm{2}}\\ ⋮\\ {\mathit{\beta }}_K\end{array}\right]+\left[\begin{array}{c}{\mathit{\epsilon }}_{\mathrm{1}}\\ {\mathit{\epsilon }}_{\mathrm{2}}\\ ⋮\\ {\mathit{\epsilon }}_K\end{array}\right]。（7）

式中：K为子模型的数量；{\mathit{y}}_k为第k个因变量的向量，k=\left\{\begin{array}{c}\mathrm{1},\begin{array}{c}\mathrm{2},\begin{array}{c}\cdots ,\begin{array}{c}K\end{array}\end{array}\end{array}\end{array}\right\}，{\mathit{y}}_k^{\mathrm{\text{'}}}=\left\{\begin{array}{c}{\mathit{y}}_{k\mathrm{1}},\begin{array}{c}{\mathit{y}}_{k\mathrm{2}},\begin{array}{c}\cdots ,\begin{array}{c}{\mathit{y}}_{kn}\end{array}\end{array}\end{array}\end{array}\right\}n为观测值的数量；{\mathit{X}}_k为块对角矩阵形式的设计矩阵；{\mathit{\beta }}_k是模型参数的向量；{\mathit{\epsilon }}_k是误差的向量，均对应于第k个因变量。系统的误差项服从多元正态分布，且具有以下性质

\mathrm{E}\left(\mathbf{\epsilon }\right)=\mathrm{0}。（8）

\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathit{\epsilon }\right)=\mathrm{E}\left(\mathit{\epsilon }\mathit{\epsilon }\mathrm{\text{'}}\left|\mathit{X}\right.\right)=

\left[\begin{array}{cccc}{\sigma }_{\mathrm{11}}{\mathit{I}}_n& {\sigma }_{\mathrm{12}}{\mathit{I}}_n& \dots & {\sigma }_{\mathrm{1}K}{\mathit{I}}_n\\ {\sigma }_{\mathrm{21}}{\mathit{I}}_n& {\sigma }_{\mathrm{22}}{\mathit{I}}_n& \dots & {\sigma }_{\mathrm{2}K}{\mathit{I}}_n\\ ⋮& ⋮& \ddots & ⋮\\ {\sigma }_{K\mathrm{1}}{\mathit{I}}_n& {\sigma }_{K\mathrm{2}}{\mathit{I}}_n& \dots & {\sigma }_{KK}{\mathit{I}}_n\end{array}\right]={\mathit{\Sigma }}_{K\times K}\otimes {\mathit{I}}_n。（9）

式中：{\mathit{\Sigma }}_{K\times K}是由子模型特定误差方差和交叉子模型误差协方差组成的正定对称矩阵； I_n 为n阶单位矩阵；n为样本容量；{\sigma }_{ij}为方程间误差协方差矩阵的元素。本研究采用广义最小二乘法（feasible generalized least squares，FGLS）法估计参数向量\widehat{\mathit{\beta }}，首先对每个子模型使用普通最小二乘法（ordinary least squares，OLS）拟合一个初始模型，以获取残差信息{\widehat{\mathit{\epsilon }}}_i，并根据残差信息，估计残差项的方差-协方差矩阵{\mathit{\Sigma }}_{K\times K}。然后使用此矩阵对原始模型进行加权最小二乘估计。

{\widehat{\mathit{\sigma }}}_{ij}=\frac{\mathrm{1}}{n}{\widehat{\mathit{\epsilon }}}_i\mathrm{\text{'}}{\widehat{\mathit{\epsilon }}}_j。（10）

\widehat{\mathit{\beta }}={\left(\mathit{X}\mathrm{\text{'}}\left({\widehat{\mathit{\Sigma }}}_{K\times K}^{-\mathrm{1}}\otimes {\mathit{I}}_n\right)\mathit{X}\right)}^{-\mathrm{1}}\mathit{X}\mathrm{\text{'}}\left({\widehat{\mathit{\Sigma }}}_{K\times K}^{-\mathrm{1}}\otimes {\mathit{I}}_n\right)\mathit{y}。（11）

式中：\widehat{\mathit{\beta }}为模型参数估计向量； X 为块对角形式的设计矩阵；{\mathit{X}}^{\text{'}}为 X 的转置矩阵；{\widehat{\mathit{\Sigma }}}_{K\times K}^{-\mathrm{1}}为{\widehat{\mathit{\Sigma }}}_{K\times K}^{}的逆矩阵； y 为堆叠的因变量观测值向量； I_n 为n阶单位矩阵；n为样本容量；{\widehat{\mathit{\sigma }}}_{ij}为方程间误差协方差矩阵元素{\sigma }_{ij}的估计值；{\widehat{\mathit{\epsilon }}}_i\mathrm{\text{'}}{\widehat{\mathit{\epsilon }}}_j为两方程残差向量的内积。

本研究借助R语言环境中的systemfit包联合估计每个子模型，构建黑河市森林蓄积量的模型预测系统，获取的误差项之间的相关性和模型之间的协方差矩阵是下一步最佳线性无偏预测器BLUP校准的关键输入，对整个预测系统的精度和可靠性具有至关重要的影响。

2.3 最佳线性无偏估计同化校正

本研究基于Henderson（1975）建立的线性混合效应模型理论体系，采用频率学派框架下的最佳线性无偏预测器（best linear unbiased prediction，BLUP）这一具有严格理论保证的统计方法，构建了面向森林蓄积量估计的多源数据同化系统。从概率论的角度看，BLUP方法本质上是多元正态分布条件下条件期望的显式解，对于一个L维的随机向量 y 被分为两部分，其数学基础可表述为

\mathit{y}=\left[\begin{array}{c}{\mathit{y}}_{\mathrm{1}}\\ {\mathit{y}}_{\mathrm{2}}\end{array}\right]。（12）

式中： y₁为待预测的变量； y₂为观测的变量或已知的辅助变量，在联合分布\mathit{y}~N\left(\mathit{\mu },\mathit{\Sigma }\right)的前提下，通过求解条件{\mathit{y}}_{\mathrm{1}}\left|{\mathit{y}}_{\mathrm{2}}\right.实现最优预测。具体而言，设随机向量 y 具有如下分布特性。

\mathrm{E}\left({\mathit{y}}_{\mathrm{1}}\right)={\mathit{\mu }}_{\mathrm{1}}。（13）

\mathrm{E}\left({\mathit{y}}_{\mathrm{2}}\right)={\mathit{\mu }}_{\mathrm{2}}。（14）

\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathit{y}}_{\mathrm{1}}\right)={\mathit{\Omega }}_{\mathrm{1}}。（15）

\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathit{y}}_{\mathrm{2}}\right)={\mathit{\Omega }}_{\mathrm{2}}。（16）

\mathrm{C}\mathrm{o}\mathrm{v}\left({\mathit{y}}_{\mathrm{1}},{\mathit{y}}_{\mathrm{2}}\right)={\mathit{\Omega }}_{\mathrm{12}}。（17）

基于上述理论框架，随机向量 y 的总体协方差结构可严格表述为以下分块矩阵形式。

\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathit{y}\right)=\left[\begin{array}{cc}{\mathit{\Omega }}_{\mathrm{1}}& {\mathit{\Omega }}_{\mathrm{12}}\\ {\mathit{\Omega }}_{\mathrm{12}}^{\mathrm{T}}& {\mathit{\Omega }}_{\mathrm{2}}\end{array}\right]。（18）

式中：{\mathit{\mu }}_{\mathrm{1}}为目标变量的均值；{\mathit{\mu }}_{\mathrm{2}}为辅助变量的均值；{\mathit{\Omega }}_{\mathrm{1}}为目标变量的协方差矩阵；{\mathit{\Omega }}_{\mathrm{12}}是目标变量与辅助变量的协方差；{\mathit{\Omega }}_{\mathrm{2}}是辅助变量的协方差矩阵。当 y₂被观测后，可以使用最佳线性无偏预测器（式中记为B_LUP）来对 y₁进行预测。最优预测器为

B_{\mathrm{L}\mathrm{U}\mathrm{P}}\left({\mathit{y}}_{\mathrm{1}}\right)={\mathit{\mu }}_{\mathrm{1}}+{\mathit{\Omega }}_{\mathrm{12}}{\mathit{\Omega }}_{\mathrm{2}}^{-\mathrm{1}}\left({\mathit{y}}_{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}\right)。（19）

BLUP预测器的预测误差方差矩阵可严格表述为以下形式。

\mathrm{V}\mathrm{a}\mathrm{r}(B_{\mathrm{L}\mathrm{U}\mathrm{P}}({\mathit{y}}_{\mathrm{1}})-y_{\mathrm{1}})={\mathit{\Omega }}_{\mathrm{1}}-{\mathit{\Omega }}_{\mathrm{12}}{\mathit{\Omega }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{\Omega }}_{\mathrm{12}}^{\mathbf{T}}。（20）

基于观测变量 y₂的信息，BLUP能够对目标变量 y₁的预测进行最优调整，精准修正预测偏差，并显著降低预测误差方差，从而实现预测精度的显著提升。本研究采用R语言lme4包中的lmer函数构建了包含随机效应的SUR模型系统，并系统性地实施了最佳线性无偏预测（BLUP）校准算法，实现了森林蓄积量的精准校准。

2.4 结果验证

模型结果采用拟合优度R²和均方根误差RMSE（式中记为R_MSE）来评估校准前和校准后的预测精度。

R^{\mathrm{2}}=\mathrm{1}-\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(y_i-{\widehat{y}}_{\mathrm{i}}\right)}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^n}\left(y_i-\overline{y}\right)}。（21）

R_{\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt[]{\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}{\left({\widehat{y}}_{\mathrm{i}}-y_i\right)}^{\mathrm{2}}}{n}}。（22）

R_{\mathrm{M}\mathrm{S}\mathrm{E}\mathrm{\%}}=\frac{R_{\mathrm{M}\mathrm{S}\mathrm{E}}}{\overline{y}}\times \mathrm{100}。（23）

式中：n为样本量；{\widehat{y}}_i为第i个蓄积量的预测值；y_i为第i个蓄积量的实测值；\overline{y}为蓄积量实测值的平均值；R_MSE%为相对均方根误差。

2.5 基于模型的总体均值推断

在森林调查中，模型及模型预测是以总体参数的置信区间形式进行推断的中间产物。因此，有必要证明所提出的DA校准程序对总体参数置信区间的影响。本研究主要关注模型估计量\widehat{\mathit{\mu }}，基于总体全部单元的总体均值真值\mathit{\mu }可以表示为

\mathit{\mu }=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{y}}_i。（24）

式中：N为研究区内总体单元个数； y_i为第i个单元的蓄积量真值。

\widehat{\mathit{\mu }}=\frac{\mathrm{1}}{n}{\displaystyle \sum _i^n}{\widehat{\mathit{y}}}_i。（25）

式中：{\widehat{\mathit{y}}}_i为模型对第i个总体单元的蓄积量预测值。

其置信区间通过自助重抽样法（Bootstrap）过程获得的\mathrm{V}\widehat{\mathrm{a}}\mathrm{r}\left(\widehat{\mathit{\mu }}\right)来构造。Bootstrap方法是基于重采样的概念，应用Bootstrap残差程序来构造多元线性回归模型的重采样，重点是对残差\mathit{\epsilon }进行重采样，并从N\left(\mathrm{0},{\widehat{\mathit{\Sigma }}}_{K\times K}\right)中抽取随机样本，然后将其添加回多元线性回归模型中，以形成重采样或Bootstrap样本。对于每个重采样，重新估计主模型参数，然后估计第k个因变量的重采样后的\widehat{\mathit{\mu }}。Bootstrap方差估计为

\mathrm{V}\widehat{\mathrm{a}}\mathrm{r}\left(\widehat{\mathit{\mu }}\right)=\frac{\mathrm{1}}{B-\mathrm{1}}{\displaystyle \sum _{b=\mathrm{1}}^B}{\left({\widehat{\mathit{\mu }}}^{\mathrm{*}}-{\widehat{\mathit{\mu }}}_b^{\mathrm{*}}\right)}^{\mathrm{2}}。（26）

{\widehat{\mathit{\mu }}}^{\mathrm{*}}=\frac{\mathrm{1}}{B}{\displaystyle \sum _{b=\mathrm{1}}^B}\left({\widehat{\mathit{\mu }}}_b^{\mathrm{*}}\right)。（27）

式中：{\widehat{\mathit{\mu }}}_b^{\mathrm{*}}是第b个重采样的Bootstrap估计，b=1，2，…，B，且B=1 000。校准后BLUP的Bootstrap方差可能向下偏移，这将产生过窄且偏离中心的置信区间；对于小样本量这种影响更加明显。为了克服这个问题，研究使用了拓展参数Bootstrap程序^［22］。具体而言，首先将SUR模型拟合到初始样本数据 y 上，以获得协方差矩阵{\widehat{\mathit{\Sigma }}}_{K\times K}和回归系数\widehat{\mathit{\beta }}的估计值。随后，生成一个与{\mathit{y}}_{\mathrm{2}}维度相同的误差向量{\mathit{\epsilon }}_{O\mathrm{b}\mathrm{s}}^{\mathrm{*}}，{\mathit{\epsilon }}_{O\mathrm{b}\mathrm{s}}^{\mathrm{*}}的元素独立地从N\left(\mathrm{0},{\widehat{\mathit{\Omega }}}_{\mathrm{2}}\right)分布中抽取。基于{\mathit{y}}_{\mathrm{2}}^{\mathrm{*}}={\mathit{X}}_{y_{\mathrm{2}}}{\widehat{\mathit{\beta }}}_{y_{\mathrm{2}}}+{\mathit{\epsilon }}_{Obs}^{\mathrm{*}}构建Bootstrap数据{\mathit{y}}_{\mathrm{2}}^{\mathrm{*}}，并进一步构建完整的Bootstrap样本{\mathit{y}}^{\mathrm{*}}。在Bootstrap样本{\mathit{y}}^{\mathrm{*}}上，再次拟合SUR模型，以获得新的估计值{\widehat{\mathit{\beta }}}^{\mathrm{*}}\left({\widehat{\mathit{\Sigma }}}_{K\times K}\right)和{\widehat{\mathit{\Sigma }}}_{K\times K}^{\mathrm{*}}。基于这些估计值，计算自助法最佳线性无偏预测（Bootstrap BLUP），得到{\mathit{\theta }}^{\mathrm{*}}\left(\widehat{\mathit{\Omega }}\right)和{\widehat{\mathit{\Sigma }}}_{K\times K}^{\mathrm{*}}，{\mathit{\theta }}^{\mathrm{*}}(\widehat{\mathit{\Omega }})为第b次Bootstrap的固定效应参数估计，其生成依赖于原始协方差\widehat{\mathit{\Omega }}。这一过程重复B次，每次重复都计算{\widehat{\mathit{\mu }}}^{\mathrm{*}}={\mathit{I}}^{\mathrm{\text{'}}}{\mathit{\theta }}^{\mathrm{*}}\left(\widehat{\mathit{\Omega }}\right){\left({\mathit{I}}^{\mathrm{\text{'}}}\mathit{I}\right)}^{-\mathrm{1}}和{\widehat{\mathit{\mu }}}_b^{\mathrm{*}}={\mathit{I}}^{\mathrm{\text{'}}}{\mathit{\theta }}^{\mathrm{*}}\left({\widehat{\mathit{\Omega }}}^{\mathrm{*}}\right){\left({\mathit{I}}^{\mathrm{\text{'}}}\mathit{I}\right)}^{-\mathrm{1}}，其中，b表示第b次Bootstrap重复， I 是n×1的单位矢量。对模型参数的不确定性进行量化，基于Bootstrap的方差估计\mathrm{V}\widehat{\mathrm{a}}\mathrm{r}(\widehat{\mathit{\mu }})及其标准误\mathrm{S}\widehat{\mathrm{D}}(\widehat{\mathit{\mu }})表示如下。

\mathrm{V}\widehat{\mathrm{a}}\mathrm{r}\left(\widehat{\mathit{\mu }}\right)=\frac{\mathrm{1}}{B-\mathrm{1}}{\displaystyle \sum _{b=\mathrm{1}}^B}({\widehat{\mathit{\mu }}}^{*}-{\widehat{\mathit{\mu }}}_b^{*}{)}^{\mathrm{2}}。（28）

\mathrm{S}\widehat{\mathrm{D}}\left(\widehat{\mathit{\mu }}\right)=\mathrm{V}\widehat{\mathrm{a}}\mathrm{r}{\left(\widehat{\mathit{\mu }}\right)}^{\mathrm{1}/\mathrm{2}}=\sqrt[]{\frac{\mathrm{1}}{B-\mathrm{1}}{\displaystyle \sum _{b=\mathrm{1}}^B}({\widehat{\mathit{\mu }}}^{\mathrm{*}}-{\widehat{\mathit{\mu }}}_b^{\mathrm{*}}{)}^{\mathrm{2}}}。（29）

3 结果与分析

3.1 似乎不相关回归（SUR）建模

3.1.1 多元线性回归模型拟合结果

采用逐步回归的方法来构建多元线性回归模型，得到基于Landsat-7遥感影像的多元线性回归方程，见公式（30）和表3，残差分布图显示无明显异方差现象，见图2。

V_{\mathrm{2010}}=\mathrm{95.702}+\mathrm{30.204}\times D_{\mathrm{E}\mathrm{M}}-\mathrm{130.510}\times V_{\mathrm{A}}-\mathrm{102.388}\times B_{\mathrm{7}}。（30）

式中：D_EM为海拔；V_A为方差纹理特征；B₇为Landsat-7第7波段反射率。

模型P<0.001，说明回归模型是成立的；调整决定系数R²=0.200，模型的拟合效果不太理想；RMSE为23.207 m³/hm²，均方根误差较大，模型预测精度较低。

3.1.2 支持性子模型拟合结果

研究以Pearson相关性分析筛选出与主模型残差具有高度相关性的潜在变量作为因变量，在R语言环境中应用截距式构建支持性子模型，得到的支持性子模型截距（S_lope）方程如下。

S_{\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}=\mathrm{0.102}\mathrm{1}。（31）

2010年的支持性子模型截距方程中参与建模的为Slope因子，其与多元线性回归模型的残差相关性最高，模型P<0.001，说明支持性子模型是成立的。

3.1.3 似乎不相关回归（SUR）联合估计结果

研究基于广义最小二乘理论框架，在R语言中的systemfit包的似乎不相关回归算法，实现了森林蓄积量估测双子系统模型的联合优化估计。不仅显著提升了参数估计效率，更重要的是精确推导出模型的残差项的方差-协方差矩阵，其数学表达形式为

\sum =\left[\begin{array}{cc}{\sigma }_{\mathrm{11}}& {\sigma }_{\mathrm{12}}\\ {\sigma }_{\mathrm{21}}& {\sigma }_{\mathrm{22}}\end{array}\right]。（32）

式中：对角线元素{\sigma }_{\mathrm{11}}和{\sigma }_{\mathrm{22}}分别表征2个子系统模型的独立变异程度，而非对角元素{\sigma }_{\mathrm{12}}和{\sigma }_{\mathrm{21}}则量化了模型间的协同变异性特征。该矩阵作为关键输入参数被引入后续的最佳线性无偏预测校准过程，其具体结果如下。

{\sum }_{\mathrm{2010}}=\left[\begin{array}{cc}\mathrm{538.525}& -\mathrm{0.147}\\ -\mathrm{0.147}& \mathrm{0.004}\end{array}\right]。（33）

通过上述估计的方差-协方差矩阵捕捉模型间的残差相关性，有效地解决了传统单方程估计中的信息损失问题；为BLUP算法提供了严格的数学基础，确保随机效应的最优线性预测满足E\left(\mathit{\epsilon }\right)=\mathrm{0}的无偏性条件；基于最大似然准则的模型拟合优度检验，该联合估计方法较传统OLS估计显著提高了参数估计效率，为多源遥感数据协同反演森林参数提供了可靠的统计框架。

3.2 最佳线性无偏预测器BLUP数据同化结果

研究借助R语言环境中的lme4包的lmer函数，对似乎不相关回归（SUR）模型中的每个子模型的因变量进行最佳线性无偏预测（BLUP）校准，构建了校正后的森林蓄积量估算模型，见公式（34）。

V_{\mathrm{2010}}=\mathrm{99.700}+\mathrm{30.204}\times D_{\mathrm{E}\mathrm{M}}-\mathrm{130.510}\times V_{\mathrm{A}}-\mathrm{102.388}\times B_{\mathrm{7}}-\mathrm{39.157}\times S_{\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}。（34）

BLUP数据同化模型的统计检验见表4，该模型的显著性检验P均为0.000（P<0.001），表明BLUP数据同化模型具有显著的统计学意义。在模型拟合效果方面，调整决定系数R²达到0.539，模型解释力和拟合效果良好。值得注意的是，该模型的预测精度显著优于多元线性回归模型，其均方根误差（RMSE）降低至17.541 m³/hm²，显示出更好的预测性能。

3.3 基于模型的总体总量推断结果

研究比较分析了多元线性回归模型和最佳线性无偏预测（BLUP）校正模型的预测性能。通过决定系数（R²）来量化模型对变异程度的解释能力、均方根误差（RMSE）和相对均方根误差（RMSE_%）来分别量化预测值与真实值之间的绝对偏差和相对偏差程度。同时，还采用了Bootstrap重抽样方法（b=1 000）分别估计了2个模型的参数方差并量化了森林蓄积量均值估计的不确定性，有效地验证了模型的稳健性和可靠性，见表5，残差分布图显示无明显异方差现象，见图3。

模型评估指标显示，BLUP校正后的模型调整决定系数（R²）从0.200提升至0.539，模型的拟合优度明显增加，模型对变异程度的解释能力显著增加；均方根误差（RMSE）由23.207降低至17.541 m³/hm²，精度改进幅度达24.4%，这说明预测值与真实值的绝对偏差程度显著降低；相对均方根误差（RMSE_%）从29.43%降至22.24%，相对偏差程度显著改善。蓄积量估计值标准差由17.367减小至12.146 m³/hm²，进一步证实了优化后的模型参数估计的稳定性增强。

经最佳线性无偏预测（BLUP）校正后的模型能够更准确地预测黑河市森林蓄积量的均值，有效校正了多元线性回归模型存在的系统性高估、低估偏差问题。这一结果表明，结合似乎不相关回归（SUR）模型与最佳线性无偏预测器（BLUP）方法构建的校正模型在森林蓄积量的预测中具有良好的可靠性和实用性。

3.4 黑龙江省黑河市森林蓄积量空间分布可视化结果

在构建BLUP校正模型的基础上，对黑龙江省黑河市2010年森林蓄积量进行可视化分析，对非林地像元进行掩膜剔除处理，绘制了相应的森林蓄积量空间分布图，见图4。黑河市森林蓄积量空间分布呈现出明显的规律性和异质性，高蓄积量区域与低蓄积量区域空间分布异质性显著，蓄积量高值地区集中分布在黑河行政区内东部和北部的小兴安岭主脉及其余脉范围，蓄积量低值地区则分散分布于中西部平原产粮地区、黑龙江等主要河流沿岸与城镇周边地区，这直观体现了人为农业生产、城镇建设活动对森林植被生长分布的干扰。空间分布图直观地、清晰地、有效地展示了研究区域内森林蓄积量的空间分布格局，为实现我国森林资源的动态监测、管理和政策制定提供了可靠的依据。

4 结论与讨论

本研究以黑龙江省黑河市各市县林区为研究对象，通过整合2010年固定样地实测数据与Landsat-7遥感影像数据，构建基于多源数据同化的森林蓄积量估算校正模型。研究采用似乎不相关回归（SUR）模型框架结合最佳线性无偏预测器（BLUP）校准方法，实现了地面实测数据与遥感数据的优势互补，建立适配黑河市的森林蓄积量年度动态估算与更新体系。

数据同化模型的调整决定系数（R²）为0.539，均方根误差（RMSE）为17.541 m³/hm²，相对均方根误差（RMSE_%）为22.24%。这些指标表明模型能够较为准确地估计森林蓄积量。经过BLUP校准后的模型在预测精度上显著优于未经校准的模型，验证了SUR-BLUP同化建模体系的科学性与可靠性。研究采用的似乎不相关回归（SUR）获取多个模型之间的潜在相关性，提高了模型的估计效率；结合最佳线性无偏预测器（BLUP）校准程序，解决了遥感估算模型空间异质性强、样地-遥感数据匹配偏差大的技术难题。研究获得黑河市2010年森林蓄积量总体均值78.859 m³/hm²，接近森林资源清查数据的年度蓄积量均值真值，实现了森林蓄积量年度动态估算和森林资源的动态监测评估，同时也为黑河市森林资源的精细化管理和可持续发展决策提供了重要的科学依据。

尽管研究在模型预测精度方面取得了较好的效果，但仍存在着不可忽视的局限性：在遥感数据选择方面，研究仅选用了Landsat-7多光谱遥感影像数据，数据维度单一、时空分辨率存在局限；在算法假设方面，BLUP校正过程所依赖的多元正态分布假设，在复杂异质的森林生态系统中无法完全满足，模型假设存在固有约束，未来可以尝试引入非参数贝叶斯方法或集成学习框架等相关内容进行补充。未来，如果能够获取黑河市“5年为周期、每年调查1/5”的样地数据，可以通过使用研究中的结合似乎不相关回归（SUR）模型与最佳线性无偏预测器（BLUP）方法进行真正意义上的年度更新。另外，SUR-BLUP模型可以进行其他区域的外推，只要研究区域有少量随机分布调查的样地数据即可。

综上，相较于依赖移动平均法估算的传统轮查体系，研究构建的SUR-BLUP联合模型有效地解决了传统模型残差系统性偏差累积、估算精度偏低的问题，通过协方差结构建模（SUR）识别系统误差源，进而采用BLUP进行最佳线性无偏估计，形成具有修正自身误差能力的动态监测体系。该模型不仅完善了我国森林资源连续清查（NFI）的技术框架，也实现了森林蓄积量监测从5年周期静态估算向年度动态精准估算的技术升级。

参考文献

[1] 吴恒.基于时空演变分析的森林资源调查抽样设计优化［D］.昆明：西南林业大学，2022.

[2] WU H.Optimization for sampling based on spatial and temporal analysis of forest resources［D］.Kunming：Southwest Forestry University，2022.

[3] 吕浩，彭春良，吴惠俐，等.森林生态系统不同碳库碳储量估算方法综述［J］.湖南林业科技，2022，49（6）：90-98.

[4] LYU H， PENG C L， WU H L，et al.A review of methods to estimate carbon storage of forest ecosystems［J］.Hunan Forestry Science & Technology，2022，49（6）：90-98.

[5] 赵升.基于森林资源遥感动态监测研究对生态环境和生态资源的影响［J］.中国林业产业，2023（11）：89-91.

[6] ZHAO S.Effects of remote sensing-based dynamic monitoring of forest resources on the ecological environment and resources［J］.China Forestry Industry，2023（11）：89-91.

[7] 曾伟生，夏锐.全国森林资源调查年度出数统计方法探讨［J］.林业资源管理，2021（2）：29-35.

[8] ZENG W S， XIA R.Discussion on methodology for generating annual estimates in national forest inventory［J］.Forest Resources Management，2021（2）：29-35.

[9] 李海月.基于遥感纹理特征的北京市森林蓄积量估测模型研究［D］.北京：北京林业大学，2020.

[10] LI H Y.Research on Beijing forest volume estimation model based on remote sensing texture features［D］.Beijing：Beijing Forestry University，2020.

[11] ROESCH F A， SCHROEDER T A， VOGT J T.Effects of cycle length and plot density on estimators for a national-scale forest monitoring sample design［J］.Forests，2017，8（9）：325.

[12] 吴恒，吕大伟，汪圣云，等.“双碳”战略背景下的多尺度森林碳储量年度监测研究进展与思考［J］.世界林业研究，2024，37（1）：59-64.

[13] WU H， LYU D W， WANG S Y，et al.Research progress in multi-scale annual monitoring of forest carbon storage under the background of “dual carbon” goal［J］.World Forestry Research，2024，37（1）：59-64.

[14] FRANKLIN S E， AHMED O S， WULDER M A，et al.Large area mapping of annual land cover dynamics using multitemporal change detection and classification of Landsat time series data［J］.Canadian Journal of Remote Sensing，2015，41（4）：293-314.

[15] SMITH W B.Forest inventory and analysis：A national inventory and monitoring program［J］.Environmental Pollution，2002，116：S233-S242.

[16] NI W J， YU T Y， PANG Y，et al.Seasonal effects on aboveground biomass estimation in mountainous deciduous forests using ZY-3 stereoscopic imagery［J］.Remote Sensing of Environment，2023，289：113520.

[17] 丁相元，陈尔学，李增元，等.国家森林资源清查遥感应用主要技术进展［J］.南京林业大学学报（自然科学版），2023，47（1）：1-12.

[18] DING X Y， CHEN E X， LI Z Y，et al.A review of remote sensing application in national forest inventory［J］.Journal of Nanjing Forestry University （Natural Science Edition），2023，47（1）：1-12.

[19] FERNÁNDEZ-DELGADO M， SIRSAT M S， CERNADAS E，et al.An extensive experimental survey of regression methods［J］.Neural Networks，2019，111：11-34.

[20] 秦思娴.陆面数据同化中的智能算法研发与实验［D］.北京：中国科学院大学，2013.

[21] QIN S X.Development and experiment of intelligent algorithm in land surface data assimilation［D］.Beijing：University of Chinese Academy of Sciences，2013.

[22] 谢元富，朱舒君，和玉君，等.数据同化百年史：演化、里程碑和未来方向［J］.气象学报，2025，83（3）：536-551.

[23] XIE Y F， ZHU S J， HE Y J，et al.A century of data assimilation：Evolution，milestones，and future directions［J］.Acta Meteorologica Sinica，2025，83（3）：536-551.

[24] 李科，苑福利，刘厂.基于海气耦合模式的均权重粒子滤波与集合卡尔曼滤波比较研究［J］.海洋通报，2021，40（4）：417-424.

[25] LI K， YUAN F L， LIU C.Comparison of particle filter and EAKF based on the ocean-atmosphere coupled model［J］.Marine Science Bulletin，2021，40（4）：417-424.

[26] 王仲睿，孙浩昊，雷荔傈，等.数据同化组成部分对初始条件和误差增长的影响［J］.中国科学：地球科学，2024，54（1）：110-121.

[27] WANG Z R， SUN H H， LEI L L，et al.The importance of data assimilation components for initial conditions and subsequent error growth［J］.Science China Earth Sciences，2024，54（1）：110-121.

[28] 曾伟生，孙乡楠，王六如，等.东北林区10种主要森林类型的蓄积量、生物量和碳储量模型研建［J］.北京林业大学学报，2021，43（3）：1-8.

[29] ZENG W S， SUN X N， WANG L R，et al.Developing stand volume，biomass and carbon stock models for ten major forest types in forest region of Northeastern China［J］.Journal of Beijing Forestry University，2021，43（3）：1-8.

[30] 辽阳市林业和草原局. 辽阳地区主要树种一元、二元立木材积表：DB2110/T 0004—2020 ［S］.辽阳：辽阳市市场监督管理局，2020.

[31] Liaoyang Municipal Forestry and Grassland Bureau. One-way and two-way tree volume tables of main tree species in Liaoyang region：DB2110/T 0004—2020 ［S］.Liaoyang：Liaoyang Municipal Bureau of Market Supervision and Administration，2020.

[32] 马友平.基于Landsat-7 ETM SLC-OFF影像的岩溶区石漠化研究［J］.湖北民族大学学报（自然科学版），2020，38（1）：13-18.

[33] MA Y P.Study on Karst rocky desertification in Karst area based on Landsat-7 SLC-OFF ETM image［J］.Journal of Hubei Minzu University （Natural Sciences Edition），2020，38（1）：13-18.

[34] 吴艳双，张晓丽.结合多尺度纹理特征的高光谱影像面向对象树种分类［J］.北京林业大学学报，2020，42（6）：91-101.

[35] WU Y S， ZHANG X L.Object-oriented tree species classification with multi-scale texture features based on airborne hyperspectral images［J］.Journal of Beijing Forestry University，2020，42（6）：91-101.

[36] 李靖.黑龙江穆棱市林业局森林蓄积量估测模型研究［D］.北京：北京林业大学，2012.

[37] LI J.Study on quantitative estimation of forest volume in Muling Forestry Bureau［D］.Beijing：Beijing Forestry University，2012.

[38] HOU Z Y， MEHTÄTALO L， MCROBERTS R E，et al.Remote sensing-assisted data assimilation and simultaneous inference for forest inventory［J］.Remote Sensing of Environment，2019，234：111431.

[39] 王鹏飞，董利虎，谢龙飞，等.基于层次贝叶斯似乎不相关回归的人工长白落叶松生物量模型构建［J］.应用生态学报，2025，36（5）：1298-1308.

[40] WANG P F， DONG L H， XIE L F，et al.Construction of biomass models for Larix olgensis plantation using hierarchical Bayesian seemingly unrelated regression［J］.Chinese Journal of Applied Ecology，2025，36（5）：1298-1308.

基金资助

“十四五”国家重点研发计划项目(2023YFD22017040302)

2025年高教强省专项项目