CN116091939A - 基于多尺度地理加权回归的森林地上生物量降尺度方法 - Google Patents

Abstract

本发明公开了基于多尺度地理加权回归的森林地上生物量降尺度方法，包括收集研究区资料，对预测变量进行提取；将全色波段与对应的多光谱波段进行融合，生成融合多光谱图像；确定与AGB值有更密切关系的特征变量，同时将多余的变量从建模过程中剔除；通过MGWR来捕捉各种预测变量的空间异质性水平的差异；利用粗分辨率数据集构建的MGWR统计回归模型直接应用于预测变量集，以完成降尺度的任务；并通过克里金法插值，将MGWRD得出的AGB残差的结构成分分离出来，并将分离出来的成分叠加到空间上对应的MGWRD预测AGB值上，形成AGB的最终分布模式。由较低成本获得的低空间分辨率遥感影像通过统计回归降尺度方法获得更高分辨率AGB分布图，并利用克里金法提高了预测精度。

Description

基于多尺度地理加权回归的森林地上生物量降尺度方法

技术领域

本发明涉及森林地上生物量估测技术领域，具体为基于多尺度地理加权回归的森林地上生物量降尺度方法。

背景技术

森林地上生物量(Above-ground Biomass，AGB)是评估森林生态系统生产力、衡量碳储量和碳汇潜力以及估计土地利用和气候变化引起的碳排放的基本指标，了解AGB的空间分布变化有助于了解全球气候变化对森林的影响，并为制定可持续森林管理战略提供基础数据支撑。尽管传统的森林AGB野外调查很准确，但这种方式收集AGB信息时存在费时、费力、成本高、破坏性大，而且在那些偏远和陡峭的山坡上很难获得等缺陷，因此很难在广泛的地区提供时空明确的AGB分布信息。遥感图像含有丰富的光谱、纹理信息，具有适当的时间和空间分辨率，覆盖面广，时效性好，在预测AGB时可以在一定程度上克服传统野外调查的缺陷。但由于数据限制、成本限制、时间和空间尺度等问题，用于AGB定量估算的遥感数据源各不相同。中等和高分辨率的遥感影像常被用于局域和区域尺度的AGB估算，对于全国甚至全球范围内的AGB估算，大多采用空间分辨率较低的遥感数据，但估算结果的精度相对较低。由于AGB明显的空间异质性以及遥感数据类型和估算方法的差异,不同的应用可能在AGB估算过程中出现不同的随机/系统误差。其中，野外实地采集样地大小和遥感图像像素大小之间的空间差异(或尺度不匹配)可能导致AGB估算误差。因此，如何从容易获得的低成本、低分辨率的遥感图像中导出更可靠、更高的空间分辨率的AGB分布图，空间降尺度是一个值得考虑的、有前景的方法。

降尺度将低分辨率的信息转化为高分辨率的信息，目前大量应用于气候变化、降水和方面的研究，但利用高分辨率多光谱数据对低分辨率AGB进行降尺度的研究很少，现有的研究主要是基于森林植被分布图给出每种植被类型的代表性生物量值，并将其降尺度为每种植被类型的平均生物量，以获得高分辨率的生物量数据。但是这种降尺度方法并不是统计上的空间降尺度方法。在地理学和空间分析领域，空间分布的变量之间关系的不稳定性被称为非平稳性。机器学习算法(如人工神经网络(ANN)、支持向量机(SVM)、随机森林(RF))可以模拟AGB与相关预测变量之间的非线性关系，从而实现从低分辨率到高分辨率的可靠降尺度。然而这些方法假设预测变量在不同空间位置对AGB的响应程度相同，由于森林生态系统具有自相关性和非平稳性，预测变量会随空间位置的变化而变化，现有大多数的AGB预测模型没有充分考虑到不同预测变量对AGB的空间分异效应的差异，得到的结果缺乏真实性。

发明内容

本发明提供基于多尺度地理加权回归的森林地上生物量降尺度方法，解决了现有大多数的AGB预测模型没有充分考虑到不同预测变量对AGB的空间分异效应的差异，得到的结果缺乏真实性的问题。

为实现上述目的，本发明提供如下技术方案：基于多尺度地理加权回归的森林地上生物量降尺度方法，包括以下步骤：

S1、收集研究区资料，对预测变量进行提取；

S2、将全色波段与对应的多光谱波段进行融合，生成融合多光谱图像；

S3、确定与AGB值有更密切关系的特征变量，同时将多余的变量从建模过程中剔除；

S4、通过MGWR模型来捕捉各种预测变量的空间异质性水平的差异；

S5、利用粗分辨率数据集构建的MGWR统计回归模型直接应用于预测变量集，以完成降尺度的任务；

S6、并通过克里金法插值，将MGWRD得出的AGB残差的结构成分分离出来，并将分离出来的成分叠加到空间上对应的MGWRD预测AGB值上，形成AGB的最终分布模式。

优选的，所述S2中，图像融合方法包括Brovey变换、GS变换、NNDiffuse方法和PCA融合，通过计算标准差、相关系数、平均梯度和信息熵等指标对融合效果进行定量评估，以确定最佳融合方法。

优选的，所述S3中，变量集的确定的方法包括多元逐步回归筛选法、随机森林重要性排序算法和皮尔森相关系数与方差膨胀因子相结合的方法。

优选的，

1)所述多元逐步回归筛选法，通过多元逐步筛选法分别用因变量建立一维回归模型，并计算每个变量对应的F值；

2)所述随机森林重要性排序算法，通过随机森林重要性排序算法进行特征重要性比较，以选择具有高重要性的预测因子；

3)所述皮尔森相关系数与方差膨胀因子相结合的方法，通过用Pearson相关系数筛选出与AGB高度相关的特征变量，然后对筛选出的特征变量进行是否违反多重共线性的检验，并用方差膨胀因子来判断：

其中，R_i是第i个变量对其余k-1个预测变量的复合相关系数，如果VIF在0到10之间，则不存在多重共线性；如果VIF≥10，说明变量之间存在较高的多重共线性，应将部分变量从模型中剔除。

优选的，所述随机森林算法包括第一指标和第二指标，用于衡量变量重要性，所述第一个指标为每棵树的预测的均方误差百分比的增加，所述第二个指标为所有树上的平均变量分裂的Inc Node Purity的总减少量。

优选的，所述S4中，MGWR模型表示如下：

其中y_i是响应变量的第i个观测值，x_ij是第j个解释变量在位置i的观测值，β_bwj代表不同变量j在不同带宽下的回归系数，(u_i,v_i)代表样本点的空间地理坐标是，k代表预测变量的数量，ε是模型回归残差。

优选的，其中AGB降尺度过程为：

AGB_high(OV)＝MGWRD(OV)

其中，AGB_high表示细尺度的AGB预测值；OV表示由MGWR分析确定的最佳预测变量，MGWRD表示由最佳预测变量集创建的AGB统计回归模型，分辨率较粗。

优选的，所述克里金插值方法包括普通克里金法和以坡向为协变量的协同克里金法。

优选的，所述S6中，AGB的最终分布模式公式为：

R(x_i)＝AGB(x_i)-AGB_MGWRD(x_i)

AGB_{MGWRD-OK/MGWRD-CK}(x_i)＝AGB_MGWRD(x_i)+R_OK/CK(x_i)

其中R(x_i)样本位置i的AGB残值，AGB(x_i)是样本位置i的AGB观测值，AGB_MGWRD(x_i)是样本位置i的MGWR预测值，AGB_{MGWRD-OK/MGWRD-CK}(x_i)是通过MGWRD-OK或MGWRD-CK模型得到的AGB预测值，R_OK/CK(x_i)是样本位置i处由OK或CK内插的残差结构成分。

与现有技术相比，本发明的有益效果：

1、本发明中，所使用的方法利用MGWR和克里金插值的整合模型，通过考虑森林生态系统复杂结构中存在的空间非平稳性和空间自相关性，由较低成本获得的低空间分辨率遥感影像通过统计回归降尺度方法获得更高分辨率AGB分布图，相较于现有许多机器学习模型，如支持向量机、随机森林、BP神经网络等，在进行AGB预测时它们往往没有考虑空间位置对AGB分布的影响，预测过程中产生了较大误差，本发明采用的方法可以改善这一缺陷，提高了预测精度。

2、本发明提出了一个考虑到空间非平稳性的统计降尺度模型框架，以解决基础数据有限和传统方法难以实现有效估算的问题，具体通过将30m AGB与预测变量之间的关系应用到15m预测变量集上，获得了较好的精度。相较于已有的AGB降尺度研究，本发明可以获得逐像元的精细AGB分布图，准确绘制该区域的AGB分布图将有助于林业和环保部门制定更具针对性的森林经营方案，更好地保护其境内的国家重点保护珍稀动植物，实现森林资源的可持续发展，主动应对全球气候变化。

3、本发明中进行了两阶段的变量选择策略，使得获得的AGB分布模式更具可靠性和真实性，通过在多元逐步回归、随机森林重要性排序和皮尔森相关系数三种变量筛选方法的基础上，进一步引入了作用尺度的概念，从空间非平稳性的角度筛选对特征变量，最终确定空间平稳性最高的变量集作为后续降尺度操作的输入预测变量集。

4、通过考虑AGB的空间自相关性，对比了对比普通克里金法和协同克里金方法两种插值方法的不同模型，包括指数模型、球状模型、高斯模型，对于AGB残差的拟合效果，通过计算各预测变量在空间上平稳程度，选取了坡向作为协变量，并将其叠加到之前的降尺度结果上，进一步提高了本发明的预测精度，这对未来预测复杂山区的AGB分布提供了参考。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。

在附图中：

图1是本发明预测变量带宽影响下的降尺度像元变化示意图；

图2是本发明前15个变量的重要性排序图；

图3是本发明所选特征变量与AGB之间的相关性热力图；

图4是本发明通过不同的变量筛选方法图：(a)随机森林重要性排序；(b)多元逐步回归；(c)Pears-Vif.GWR和MGWR确定的不同预测变量的带宽；

图5是本发明MGWRD模型的训练和验证性能图：(a)训练样本(n＝118)和(b)验证样本(n＝50)；

图6是本发明基于MGWR模型的AGB降尺度结果图；

图7是本发明MGWRD预测的AGB残差的分布直方图；

图8是本发明基于MGWRD-CK模型的改进的AGB预测结果图；

图9是本发明的方法流程图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

实施例：如图9所示，基于多尺度地理加权回归的森林地上生物量降尺度方法，包括以下步骤：

S1、收集研究区资料，对预测变量进行提取；

S2、将全色波段与对应的多光谱波段进行融合，生成融合多光谱图像；图像融合方法包括Brovey变换、GS变换、NNDiffuse方法和PCA融合，通过计算标准差、相关系数、平均梯度和信息熵等指标对融合效果进行定量评估，以确定最佳融合方法；

S3、确定与AGB值有更密切关系的特征变量，同时将多余的变量从建模过程中剔除；变量集的确定的方法包括多元逐步回归筛选法、随机森林重要性排序算法和皮尔森相关系数与方差膨胀因子相结合的方法；

1)所述多元逐步回归筛选法，通过多元逐步筛选法分别用因变量建立一维回归模型，并计算每个变量对应的F值；

2)所述随机森林重要性排序算法，通过随机森林重要性排序算法进行特征重要性比较，以选择具有高重要性的预测因子；所述随机森林算法包括第一指标和第二指标，用于衡量变量重要性，所述第一个指标为每棵树的预测的均方误差百分比的增加，所述第二个指标为所有树上的平均变量分裂的Inc Node Purity的总减少量；

3)所述皮尔森相关系数与方差膨胀因子相结合的方法，通过用Pearson相关系数筛选出与AGB高度相关的特征变量，然后对筛选出的特征变量进行是否违反多重共线性的检验，并用方差膨胀因子来判断：

其中，R_i是第i个变量对其余k-1个预测变量的复合相关系数，如果VIF在0到10之间，则不存在多重共线性；如果VIF≥10，说明变量之间存在较高的多重共线性，应将部分变量从模型中剔除；

S4、通过MGWR模型来捕捉各种预测变量的空间异质性水平的差异；MGWR模型表示如下：

其中y_i是响应变量的第i个观测值，x_ij是第j个解释变量在位置i的观测值，β_bwj代表不同变量j在不同带宽下的回归系数，(u_i,v_i)代表样本点的空间地理坐标是，k代表预测变量的数量，ε是模型回归残差；

S5、利用粗分辨率数据集构建的MGWR统计回归模型直接应用于预测变量集，以完成降尺度的任务；AGB降尺度过程为：

AGB_high(OV)＝MGWRD(OV)

其中，AGB_high表示细尺度的AGB预测值；OV表示由MGWR分析确定的最佳预测变量，MGWRD表示由最佳预测变量集创建的AGB统计回归模型，分辨率较粗；

S6、并通过克里金法插值，所述克里金插值方法包括普通克里金法和以坡向为协变量的协同克里金法，将MGWRD得出的AGB残差的结构成分分离出来，并将分离出来的成分叠加到空间上对应的MGWRD预测AGB值上，形成AGB的最终分布模式，AGB的最终分布模式公式为：

R(x_i)＝AGB(x_i)-AGB_MGWRD(x_i)

AGB_{MGWRD-OK/MGWRD-CK}(x_i)＝AGB_MGWRD(x_i)+R_OK/CK(x_i)

其中R(x_i)样本位置i的AGB残值，AGB(x_i)是样本位置i的AGB观测值，AGB_MGWRD(x_i)是样本位置i的MGWR预测值，AGB_{MGWRD-OK/MGWRD-CK}(x_i)是通过MGWRD-OK或MGWRD-CK模型得到的AGB预测值，R_OK/CK(x_i)是样本位置i处由OK或CK内插的残差结构成分。

1、实现功能或结果的具体方法步骤：

1.1数据源和预处理

本实施例使用了两景2014年7月22日Landsat 8OLI影像，条带号分别为119/040和119/041，云量分别为6.89％和4.02％。Landsat 8OLI地表反射率数据是由美国地质调查局EROS数据中心从Landsat地表反射率代码内部程序(LaSRC)直接生成的,该程序使用MODIS辅助气候数据和独特的辐射传输模型,利用沿海气溶胶带进行气溶胶反演测试。此外，通过使用C校正模型对Landsat表面反射率图像进行地形校正，以尽量减少地形阴影对地表反射率的影响。本发明还使用了(http://earthexplorer.usgs.gov)分辨率为30米的免费数字高程模型(DEM)，并在ArcGIS10.3软件包中导出了坡度和坡向信息。基于DEM数据，通过使用ENVI 5.3软件包中的RPC模块对Landsat 8OLI图像进行了正射校正，以进一步细化像元的定位精度，将研究区域内地形起伏的影响降至最低。

本实施例使用的地面森林清查数据是由浙江省2014年全国森林资源清查(NFI)结果编译而成的。浙江省的NFI始于1979年，采用系统抽样的方法建立了一系列样本面积为0.08公顷(28.28米×28.28米)的正方形样地。每个样地的AGB是由样地内每棵被测树木的AGB加总得出的，每棵被测树木的AGB是由该树种对应的异速生长方程计算出来的。对于某些没有生物量异速生长方程的树种，使用相同种属的方程作为替代来进行近似计算。为了获得更准确的反演结果，在分析中排除了Landsat 8OLI图像中被云层像素覆盖的样地。此外,虽然OLI图像的像素大小(30m×30m)与每个样地的大小(28.28m×28.28m)非常接近,但在空间尺度上仍有小的差异,为了便于建模和绘图过程,有必要将单位面积的AGB转换为像素级别的AGB,单位为t/900m²。完成后，用SPSS软件(21.0版)分析过滤后样本的AGB值的统计特征，那些大于或小于样地AGB平均值加/减3倍标准差的AGB值被视为异常值，并在分析中进一步排除。最终，共保留了168块样地，其中70％被随机确定用于模型训练，其余30％用于模型验证。下表筛选后样地AGB值的统计描述，总结了最终样本地块的AGB值的统计特征：

1.2自变量提取

本实施例中提取了四类建模变量，包括原始波段及其转换、植被指数、纹理特征和地形特征。光谱特征转换包括主成分分析(PCA)、缨帽变换(TC)和最小噪声变换(MNF)的结果，实施这些转换是为了去除原始波段之间的冗余信息，并产生那些可能与AGB高度相关的综合特征。如下表中使用的AGB建模变量，显示了生成特征的技术细节，总共有132个特征变量参与建模分析。

1.3图像融合

为了获得更高分辨率的预测变量集以支持随后的降尺度建模，采用了包括Brovey变换、GS变换、NNDiffuse方法和PCA融合在内的四种图像融合方法，将全色波段与对应的多光谱波段融合，生成空间分辨率为15米的融合多光谱图像。通过计算标准差、相关系数、平均梯度和信息熵等指标对融合效果进行定量评估，以确定最佳融合方法。在此基础上，从融合后的图像和重采样的DEM中相应地产生了那些可能包含在降尺度模型中的空间分辨率为15米的特征。

1.4变量集的确定

变量选择是指通过一定的标准来确定那些与AGB值有更密切关系的特征变量，同时将多余的变量从建模过程中剔除。为了比较不同特征或变量组合对AGB的总体影响，实施并比较了三种变量选择方法，包括多元逐步回归、随机森林重要性排序和皮尔森相关系数与方差膨胀因子相结合的方法。

1.4.1多元逐步回归(Multiple Stepwise Regression,MSR)

多元逐步筛选法结合了变量正向选择法和反向选择法的优点，分别用因变量建立一维回归模型，并计算每个变量对应的F值。在建立逐步回归模型时，在目前未纳入模型的预测变量中选择F值最大的自变量纳入模型，然后对每个被选的预测变量进行t检验，如果有被选的预测变量不再显著，则将其剔除。重复上述步骤，直到既没有选择显著的解释变量纳入回归方程，也没有从回归方程中删除不显著的解释变量，保留的最佳子集被用作建立模型的特征变量。在本分析中，逐步回归是在SPSS软件(Version 25,Armonk,NY,USA)中进行的，F检验的显著性水平被设定为0.05和0.10，用于确定变量进入和删除。

1.4.2随机森林重要度排序(Random Forest Importance Ranking,RFR)

随机森林重要性排序算法被用来进行特征重要性比较，以选择具有高重要性的预测因子。随机森林算法有两个衡量变量重要性的指标。第一个指标(从袋外数据的排名中计算)是每棵树的预测的均方误差百分比(％IncMSE)的增加，第二个指标是所有树上的平均变量分裂的Inc Node Purity的总减少量。更高的％IncMSE和Inc Node Purity值表示更重要的预测变量。

1.4.3皮尔逊相关系数(Pearson-Vif)

用Pearson相关系数筛选出与AGB高度相关的特征变量，然后对筛选出的特征变量进行是否违反多重共线性的检验，本发明用方差膨胀因子(VIF)来判断。

其中，R_i是第i个变量对其余k-1个预测变量的复合相关系数，如果VIF在0到10之间，则不存在多重共线性。如果VIF≥10，说明变量之间存在较高的多重共线性，应将部分变量从模型中剔除。

1.5多尺度地理加权回归模型(Multiscale Geographic Weighted Regression,MGWR)

与经典的GWR相比，MGWR的优势在于消除了GWR对所有变量的单一带宽的假设，从而产生了更接近实际状态的空间过程。考虑到空间异质性和作用尺度，MGWR可以应用于研究各种影响因素及其相应的过程。理论上，MGWR允许不同的预测因子有不同的空间尺度，表明每个预测因子有不同的空间平滑范围。在MGWR模型中，这种空间平滑范围由每个预测因子的带宽反映。带宽小意味着关系在相对局部范围内变化，而带宽大则表明关系在较大范围内稳定。带宽越接近全域，代表空间平滑度越高，表明预测变量在研究区域内产生的结果变化越小。在这项工作中，MGWR被用来捕捉各种预测变量的空间异质性水平的差异。MGWR模型表示如下：

这里y_i是响应变量的第i个观测值，x_ij是第j个解释变量在位置i的观测值，β_bwj代表不同变量j在不同带宽下的回归系数，(u_i,v_i)代表样本点的空间地理坐标是，k代表预测变量的数量，ε是模型回归残差。

1.6AGB降尺度

1.6.1多尺度地理加权回归降尺度模型(MGWRD模型)

MGWR模型可以通过考虑预测变量的空间平稳性确定预测变量的最佳组合，分析确定每个预测变量的空间作用尺度(带宽)。预测变量在空间位置上的变化较小，或具有较大的作用尺度(带宽)，说明这种影响在空间上是相似的，存在较小的异质性，这将导致更平滑的降尺度输出结果，如图1(b)所示。反之，当预测变量的作用尺度(带宽)较小时，将产生更多的异质性降尺度结果，如图1(a)所示。一旦确定了预测变量的最佳组合，就将利用粗分辨率数据集(30米分辨率)构建的MGWR统计回归模型直接应用于15米分辨率的预测变量集，以完成降尺度的任务。AGB降尺度过程如下：

AGB_high(OV)＝MGWRD(OV)

这里AGB_high表示细尺度的AGB预测值；OV表示由MGWR分析确定的最佳预测变量，MGWRD表示由最佳预测变量集创建的AGB统计回归模型，分辨率较粗；

1.6.2多尺度地理加权回归降尺度-协同克里金模型(MGWRD-CK模型)

由于MGWRD模型没有考虑到AGB样本的空间自相关，我们采用了MGWRD和克里金插值相结合的方法来更准确地确定AGB的空间分布。本文使用的克里金插值方法包括普通克里金法(OK)和带有协变量的协同克里金法(CK)。具体步骤如下：(1)通过普通克里金法和协同克里金法插值，将MGWRD得出的15米分辨率AGB残差的结构成分分离出来；(2)将分离出来的成分叠加到空间上对应的MGWRD预测AGB值上，形成15米AGB的最终分布模式，并通过如下两式说明了这两个步骤的具体原理：

R(x_i)＝AGB(x_i)-AGB_MGWRD(x_i)

AGB_{MGWRD-OK/MGWRD-CK}(x_i)＝AGB_MGWRD(x_i)+R_OK/CK(x_i)

其中R(x_i)样本位置i的AGB残值，AGB(x_i)是样本位置i的AGB观测值，AGB_MGWRD(x_i)是样本位置i的MGWR预测值。AGB_{MGWRD-OK/MGWRD-CK}(x_i)是通过MGWRD-OK或MGWRD-CK模型得到的AGB预测值，R_OK/CK(x_i)是样本位置i处由OK或CK内插的残差结构成分。在此，我们重点比较了普通克里金(OK)插值和以坡向为协变量的协同克里金(CK)插值的残差预测精度，我们选择了指数模型、球形模型和高斯模型的半方差函数来拟合残差。为了测试MGWR算法(公式2)的性能，我们将其分析结果与随机森林模型(RF)、地理加权回归(GWR)和传统的普通最小二乘回归(OLS)方法进行了比较。

1.7交叉验证

这里，采用十倍交叉验证的方法来评估四个模型基于不同预测变量的平均预测性能。在此基础上，包括平均绝对误差(MAE)、均方根误差(RMSE)、决定系数平均值(R2)和偏差在内的统计数据从十倍交叉验证中得出，以显示模型的表现。其公式分别如下：

其中n是验证观测值的总数，

是模型预测的AGB，y_i是地面观测的AGB，

是所有观测的AGB值的算术平均值。

2、获得产品或结果的性能的实验数据：

2.1预测变量集的确定

2.1.1多元逐步回归(MSR)

如下表多元逐步回归结果，总结了通过MSR方法选择的预测变量，以及相应的拟合回归方程。MSR最终保留了12个预测变量，纹理的数量占比最大，一些光谱特征、地形因子(包括坡度和坡度)也被选中；

注:Y代表AGB,X₁,X₂,X₃,…,X₁₁分别为slope,cor_B1_3,B15,EVI,sec_B4_3,aspect,var_B4_3,dis_B4_3,homo_B6_5,SAVI,B4和mean_B5_3

2.1.2随机森林重要性排序(RFR)

图2显示了RF模型运行100次后前15个变量的重要性得分。通过比较IncMSE和IncNodePurity指标，将con_B4_3、B2和con_B2_5从前15名中删除，最后选择elevation、aspect、sec_B3_5、B357、B14、B7、var_B5_3、mean_B7_5、亮度、B16、B5和cor_B5_5作为输入变量进行后续分析。

2.1.3皮尔森相关系数(Pearson-Vif)

图3显示了通过Pearson-Vif方法得出的前12个预测变量与AGB之间的相关性热力图，包括两个地形因素、五个纹理特征和五个光谱特征。其中，两个地形因素，即坡向和海拔，与AGB的相关性最强，而其余的预测变量与AGB的相关性没有明显差异。

2.2作用尺度

在构建MGWR模型之前，对样本之间进行了空间自相关检验，以确定是否存在显著关系。结果，计算出Moran'I指数为0.23，Z-score为3.6，P值为0.0003。Z-score和P-value的结果共同表明，研究区域内的AGB呈现出统计学上显著的聚集模式，而不是随机分布模式。如下表三种变量筛选方法确定的不同自变量集下OLS、RF、GWR和MGWR模型的预测性能指标，显示了三组预测变量下，四种AGB建模算法的不同预测性能。在实施OLS、GWR、MGWR和RF回归分析时，RFR选择的变量集获得的RSS、AIC、AICc和CV值最低，RFR选择的变量集的拟合度(R平方值)也相应最高。因此，我们认为RFR选择的变量集是最佳的变量组合，对AGB变化具有最高的解释能力。在此基础上，比较OLS、GWR、MGWR和RF的回归结果，MGWR在RSS、AIC、AICc、CV和R平方统计方面优于OLS、GWR和RF。因此，进一步得出结论，MGWR回归是目前分析中的最佳回归；

注：RSS为残差平方和；AIC为赤池信息量准则；AICc为更正的赤池信息量准则；CV为变异系数；R square为决定系数("*"表示在5％水平上有统计学意义)。

如图4显示了GWR和MGWR对不同预测变量所确定的作用尺度或带宽。在不同的自变量组中，经典GWR的识别带宽在121到130公里之间。这种单一的带宽假定所有变量在同一区域范围内影响AGB，这具有有很大的局限性。相反，MGWR可以捕捉到预测变量的空间异质性水平的差异。例如，在图5(a)中，海拔和坡向的带宽较小，表明这两个变量分别在相对局部的尺度上影响AGB。sec_B3_5、NDVI、var_B5_3、B7、mean_B4_5、B16和AGB这六个变量之间的关系表现出空间平滑性，但在广泛的区域范围内，其过程有所不同。其他变量在全球范围内影响AGB，因为它们的最佳带宽接近最大可能的邻居数168(总样本量)。RFR选择的预测变量集表现出更大的作用规模，平均为130公里，中位数为139公里，更接近全球范围或总样本量(图5(a))。这些结果进一步证实了RFR选择的预测变量集是最好的，可以产生较少的空间异质性降尺度结果。

2.3多尺度地理加权回归降尺度模型(MGWRD模型)

在当前方法中，一直采用十倍交叉验证的策略。图5(a)、(b)显示了基于RFR选择的变量建立的MGWRD模型的拟合精度和验证精度。如图5所示，MGWRD模型的拟合精度和验证精度(R²)分别为0.57和0.58；偏差值分别为-2.78t/ha和-2.74t/ha。这说明本发明构建的MGWRD模型在这里没有出现过拟合现象，取得了较高的AGB估计精度。

如图6显示了MGWRD预测的AGB空间分布模式。预测的AGB值在12.53-306.81吨/公顷之间，预测的AGB低值主要分布在研究区的北部、西南部、东南部和中部，预测的高值主要集中在研究区的南部、中东部和中西部地区。

2.3多尺度地理加权回归降尺度-协同克里金模型(MGWRD-CK模型)

如图7显示了MGWRD预测的AGB残差的柱状图。残差的标准差为63.36吨/公顷，绝对峰度值为-0.36，绝对偏度接近1，表明残差近似于正常分布。因此，满足了进行克里金插值的前提假设，可以对AGB预测残差进行后续的克里金插值分析。从图4中可以看出，三组特征变量的作用尺度显示，坡向具有很强的空间异质性，因此，选择它作为协变量进行克里金差值。

如下表OK和CK方法下不同模型的半方差函数统计参数，显示了基于指数模型、球形模型和高斯模型的半方差函数模型残差拟合结果和相关统计。根据GS+建模的结果，总体而言，指数模型在OK和CK两种方法下都优于其他两个模型，对于CK方法，指数模型取得了更高的精度，在下表中R square为0.81，高于0.78(OK)。因此，CK指数函数模型被选择用以进一步提高降尺度AGB的精度；

如图8显示了基于MGWRD-CK模型的改进的AGB预测结果。显然，MGWRD-CK的预测(图8)与MGWRD的预测(图6)之间有明显的空间相似性，而在MGWRD-CK的预测中，其最小AGB下降到5.82吨/公顷，最大AGB增加到328.18吨/公顷。与MGWR模型的预测相比，MGWRD-CK模型的AGB预测的范围扩大了，表明该模型对AGB预测的适应性更强。利用30％的独立样本对模型进行了验证，相关统计数据见下表基于验证数据的MGWRD-OK和MGWRD-CK模型的验证精度统计，显示了MGWRD-OK和MGWRD-CK模型在MAE、RMSE、R square和Bias等具体统计数据方面比MGWRD模型的精度提高；

目前常见的传感器提供的数据往往存在空间分辨率高而时间分辨率低，或时间分辨率高而空间分辨率低的矛盾，这种时空分辨率间的权衡限制了同时具有高频率和高空间分辨率遥感数据的获取。例如，在估算较小区域尺度森林地上生物量时，学者们大多采用高分系列卫星等数据，这些数据往往存在获取成本高、工作量大的特点。而对于国家尺度甚至是全球尺度范围内的生物量估算，多采用时间分辨率高但空间分辨率低的遥感数据，其估算结果的精度比较低，不能细致地呈现森林地上生物量的空间分布格局，不能满足森林调查及动态变化监测的需要。基于此，本发明将统计回归降尺度的思想应用于森林地上生物量数据来以较低成本生成高空间分辨率的AGB空间模式图。近年来，地表温度、降水等领域的降尺度方法研究成果丰富，但针对森林地上生物量的降尺度研究很少，已有的研究也多是基于森林植被分布图，给出各植被类型的生物量代表值，将其作为各植被类型在大尺度上的平均生物量降尺度到格网上，从而得到高分辨率生物量数据，这并非是传统意义上的统计降尺度，而是基于赋值策略的降尺度，不能有效表达AGB局部细节特征。本发明提出了一个考虑到空间非平稳性的统计降尺度模型框架，以解决基础数据有限和传统方法难以实现有效估算的问题。具体来说，将30m AGB与预测变量之间的关系应用到15m预测变量集上，获得了较好的精度。相较于已有的AGB降尺度研究，本发明可以获得逐像元的精细AGB分布图，准确绘制该区域的AGB分布图将有助于林业和环保部门制定更具针对性的森林经营方案，更好地保护其境内的国家重点保护珍稀动植物，实现森林资源的可持续发展，主动应对全球气候变化。

由于森林生态系统的复杂结构以及遥感数据类型和估算方法的差异，不同的应用可能在AGB估算过程中出现不同的随机/系统误差。其中，原地样本地块大小和遥感图像像素大小之间的空间差异(或尺度不匹配)可能导致AGB估算误差。例如，在使用标准尺寸估计较大区域的AGB时，局部空间变化会导致较大的抽样误差。本发明中，遥感影像数据的分辨率为30m×30m，将28.28m×28.28m的地面调查数据按单位面积转换为30m×30m的生物量实测数据，通过该转换，尽可能的减少尺度不匹配带来的误差。

本发明中进行了两阶段的变量选择策略，使得获得的AGB分布模式更具可靠性和真实性。森林生态系统的复杂结构产生的空间非平稳性是不能忽视的问题，为了获得更加平滑的降尺度结果，在多元逐步回归、随机森林重要性排序和皮尔森相关系数三种变量筛选方法的基础上，进一步引入了作用尺度的概念，从空间非平稳性的角度筛选对特征变量，最终确定空间平稳性最高的变量集作为后续降尺度操作的输入预测变量集。现有的AGB分布预测模型，大都选择一种变量筛选方法，可能会存在不同预测变量之间出现共线性、冗余性的问题。同时，现有的AGB预测模型大都没有将特征变量的空间属性纳入考虑之中，而本方法对比了三种特征变量筛选方法，同时将空间位置属性考虑在内，会有效提高AGB预测的精度。除此之外，通过考虑特征变量的作用尺度，该方法可以统计各特征变量在空间上对AGB的响应程度，以此来分析哪些为局部变量，哪些为全局变量，可以为今后关于AGB预测变量的筛选提供参考。

本发明通过考虑AGB的空间自相关性，对比了对比普通克里金法和协同克里金方法两种插值方法的不同模型(指数模型、球状模型、高斯模型)对于AGB残差的拟合效果，通过计算各预测变量在空间上平稳程度，选取了坡向作为协变量，并将其叠加到之前的降尺度结果上，进一步提高了本发明的预测精度。这对未来预测复杂山区的AGB分布提供了参考。

最后应说明的是：以上所述仅为本发明的优选实例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.基于多尺度地理加权回归的森林地上生物量降尺度方法，包括以下步骤：

S1、收集研究区资料，对预测变量进行提取；

S2、将全色波段与对应的多光谱波段进行融合，生成融合多光谱图像；

S3、确定与AGB值有更密切关系的特征变量，同时将多余的变量从建模过程中剔除；

S4、通过MGWR模型来捕捉各种预测变量的空间异质性水平的差异；

S5、利用粗分辨率数据集构建的MGWR统计回归模型直接应用于预测变量集，以完成降尺度的任务；

S6、并通过克里金法插值，将MGWRD得出的AGB残差的结构成分分离出来，并将分离出来的成分叠加到空间上对应的MGWRD预测AGB值上，形成AGB的最终分布模式。

2.根据权利要求1所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述S2中，图像融合方法包括Brovey变换、GS变换、NNDiffuse方法和PCA融合，通过计算标准差、相关系数、平均梯度和信息熵等指标对融合效果进行定量评估，以确定最佳融合方法。

3.根据权利要求1所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述S3中，变量集的确定的方法包括多元逐步回归筛选法、随机森林重要性排序算法和皮尔森相关系数与方差膨胀因子相结合的方法。

4.根据权利要求3所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：

1)所述多元逐步回归筛选法，通过多元逐步筛选法分别用因变量建立一维回归模型，并计算每个变量对应的F值；

2)所述随机森林重要性排序算法，通过随机森林重要性排序算法进行特征重要性比较，以选择具有高重要性的预测因子；

3)所述皮尔森相关系数与方差膨胀因子相结合的方法，通过用Pearson相关系数筛选出与AGB高度相关的特征变量，然后对筛选出的特征变量进行是否违反多重共线性的检验，并用方差膨胀因子来判断：

其中，R_i是第i个变量对其余k-1个预测变量的复合相关系数，如果VIF在0到10之间，则不存在多重共线性；如果VIF≥10，说明变量之间存在较高的多重共线性，应将部分变量从模型中剔除。

5.根据权利要求4所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述随机森林算法包括第一指标和第二指标，用于衡量变量重要性，所述第一个指标为每棵树的预测的均方误差百分比的增加，所述第二个指标为所有树上的平均变量分裂的Inc Node Purity的总减少量。

6.根据权利要求1所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述S4中，MGWR模型表示如下：

其中y_i是响应变量的第i个观测值，x_ij是第j个解释变量在位置i的观测值，β_bwj代表不同变量j在不同带宽下的回归系数，(u_i，v_i)代表样本点的空间地理坐标是，k代表预测变量的数量，ε是模型回归残差。

7.根据权利要求1所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述S5中，AGB降尺度过程为：

AGB_high(OV)＝MGWRD(OV)

其中，AGB_high表示细尺度的AGB预测值；OV表示由MGWR分析确定的最佳预测变量，MGWRD表示由最佳预测变量集创建的AGB统计回归模型，分辨率较粗。

8.根据权利要求1所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述克里金插值方法包括普通克里金法和以坡向为协变量的协同克里金法。

9.根据权利要求1所述的基于多尺度地理加权回归的森林地上生物量降尺度方法，其特征在于：所述S6中，AGB的最终分布模式公式为：

R(x_i)＝AGB(x_i)-AGB_MGWRD(x_i)

AGB_{MGWRD-OK/MGWRD-CK}(x_i)＝AGB_MGWRD(x_i)+R_OK/CK(x_i)

其中R(x_i)样本位置i的AGB残值，AGB(x_i)是样本位置i的AGB观测值，AGB_MGWRD(x_i)是样本位置i的MGWR预测值，AGB_{MGWRD-OK/MGWRD-CK}(x_i)是通过MGWRD-OK或MGWRD-CK模型得到的AGB预测值，R_OK/CK(x_i)是样本位置i处由OK或CK内插的残差结构成分。

Images