CN112395808A - 一种结合随机森林和协同克里金的生物量遥感制图方法 - Google Patents
一种结合随机森林和协同克里金的生物量遥感制图方法 Download PDFInfo
- Publication number
- CN112395808A CN112395808A CN202011299207.5A CN202011299207A CN112395808A CN 112395808 A CN112395808 A CN 112395808A CN 202011299207 A CN202011299207 A CN 202011299207A CN 112395808 A CN112395808 A CN 112395808A
- Authority
- CN
- China
- Prior art keywords
- agb
- model
- variables
- value
- variable
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000007637 random forest analysis Methods 0.000 title claims abstract description 81
- 239000002028 Biomass Substances 0.000 title claims abstract description 42
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000013507 mapping Methods 0.000 title claims abstract description 25
- 238000005070 sampling Methods 0.000 claims abstract description 14
- 230000003595 spectral effect Effects 0.000 claims abstract description 13
- 230000009466 transformation Effects 0.000 claims description 12
- 238000012549 training Methods 0.000 claims description 11
- 238000012216 screening Methods 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 8
- 238000002310 reflectometry Methods 0.000 claims description 8
- 238000012795 verification Methods 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 238000000605 extraction Methods 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 6
- 230000010287 polarization Effects 0.000 claims description 6
- 238000005259 measurement Methods 0.000 claims description 5
- 238000012163 sequencing technique Methods 0.000 claims description 5
- 229910052799 carbon Inorganic materials 0.000 claims description 4
- 238000001228 spectrum Methods 0.000 claims description 4
- 230000000873 masking effect Effects 0.000 claims description 3
- 238000000513 principal component analysis Methods 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 238000012544 monitoring process Methods 0.000 abstract description 4
- 238000011835 investigation Methods 0.000 abstract description 3
- 230000007613 environmental effect Effects 0.000 abstract description 2
- 238000002360 preparation method Methods 0.000 abstract 1
- 238000012360 testing method Methods 0.000 abstract 1
- 230000006870 function Effects 0.000 description 8
- 238000004422 calculation algorithm Methods 0.000 description 7
- 230000003287 optical effect Effects 0.000 description 6
- 238000011160 research Methods 0.000 description 6
- 238000012937 correction Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000007781 pre-processing Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000012417 linear regression Methods 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 3
- 238000010200 validation analysis Methods 0.000 description 3
- 241000196324 Embryophyta Species 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 230000008520 organization Effects 0.000 description 2
- 238000012827 research and development Methods 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 101150021453 ARI1 gene Proteins 0.000 description 1
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 102100031383 Fibulin-7 Human genes 0.000 description 1
- 101000846874 Homo sapiens Fibulin-7 Proteins 0.000 description 1
- 240000002044 Rhizophora apiculata Species 0.000 description 1
- 208000027066 STING-associated vasculopathy with onset in infancy Diseases 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000003066 decision tree Methods 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 101150009632 prx2 gene Proteins 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/243—Classification techniques relating to the number of classes
- G06F18/24323—Tree-organised classifiers
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/40—Analysis of texture
- G06T7/41—Analysis of texture based on statistical description of texture
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/90—Determination of colour characteristics
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Mathematical Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computational Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Mathematical Optimization (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Pure & Applied Mathematics (AREA)
- Artificial Intelligence (AREA)
- Mathematical Physics (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Remote Sensing (AREA)
- Computer Graphics (AREA)
- Medical Informatics (AREA)
- Computer Hardware Design (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种结合随机森林和协同克里金的生物量遥感制图方法，属于森林资源监测和环境因子调查领域。本发明主要包括：利用国家森林资源连续清查数据，并结合提取样地点在多源遥感影像中对应的各个光谱信息、纹理信息、地形因子作为生物量建模的备选变量；设计了一种随机森林环境下重要变量的选择视图，同时将两种重要度类型综合反映出变量对生物量的预测能力，便于本方案的建模准备；对每个样地点随机森林反演得到的结果与实测数据的残差值，将高程作为协变量，残差作为主变量进行协同克里金插值，模拟得到整个区域随机森林预测结果的残差；采用划分好的测试样本集验证生物量反演精度；区域生物量空间分布制图。
Description
技术领域
本发明属于森林资源监测和环境因子调查领域，更具体地说，涉及一种结合随机森林和协同克里金的生物量遥感制图方法。
背景技术
光学遥感数据是最广泛获得的数据类型，Landsat(美国陆地卫星)属于其中的中分辨率卫星，自1972年发射以来，已经积累了大量的地球表面遥感影像数据，Landsat凭借连续观测周期、适中的空间分辨率、高质量、科学的数据处理及存档等优势，在土地覆盖变化分析，植被生长监测，农作物估产等方面发挥着重要甚至是不可替代的作用。其影像数据自2008年以来可免费获得，有16天的覆盖周期，覆盖范围广泛，30m的空间分辨率与中国国家森林清查的25m×25m样地相近，其历史档案影像数据可以从网站上(https：//glovis.usgs.gov/)免费下载。在结构复杂的森林中单一使用光学遥感图像时，容易遇到光谱饱和，合成孔径雷达(SAR)的高穿透能力允许对植物的结构参数进行广泛的信息提取。2007年，先进陆地观测卫星(ALOS)相控阵L波段合成孔径雷达(PALSAR)的成功发射使得对雷达数据的使用进一步增加。25m空间分辨率的ALOS/PALSAR数据由日本宇宙航空研究开发机构(JAXA)进行预处理，公开在网站(https：//www.eorc.jaxa.jp/ALOS/en/)供免费下载。
多源遥感图像结合适当的统计建模算法被认为是一种有效和可靠的手段，以模拟得到大范围森林地上生物量(以下简称生物量或AGB)。利用遥感数据和实地森林调查数据来反演AGB的方法有很多，主要有线性回归模型、非参数估计模型或者地统计模型。Foody等在2001年发表文章，利用Landsat TM数据反演了婆罗洲地区生物量，并对比了人工神经网络法与多元逐步回归法在森林生物量估测中的效果；Lefsky等在2005年发表文章，基于机载雷达数据，使用多元线性归回归模型估测了喀斯喀拉山脉西坡的森林生物量；曹庆先等2011年发表的文章基于Landsat TM影像，使用KNN算法估算了海南和广西地区红树林的生物量；宋茜、王晓宇等利用ALOS PALSAR双极化雷达数据分别反演了大兴安岭和云南山区的森林生物量。Chen等2010年发表的文章结合雷达数据和Quickbird影像，使用支持向量回归的方法对加拿大温哥华岛的地上生物量进行了反演；Powell等在2010年发表文章结合样地数据与Landsat时间序列数据，比较了随机森林与另两种模型对于预测森林地上生物量的效果。沈文娟等以1988-2017年的广东省连清数据为基础，结合Landsat数据与ALOS数据，使用随机森林算法估测了粤北地区的森林生物量；王云飞等基于随机森林算法对景洪市的橡胶林地上生物量进行了反演估测，分析了随机森林算法在光学遥感估测生物量上的表现；张峰等基于地统计分析和浙江省连清数据对浙江省森林碳空间分布进行了研究；贺鹏等基于二类调查的局级固定样地数据，比较了地统计学中普通克里金法和协同克里金法对吉林省金仓林场森林地上生物量估计的效果。
随机森林(random forest，RF)模型是由Leo Breiman和Adele Cutler在2001年提出的一种基于决策树的分类和回归的机器学习算法。克里金插值(Kriging)也称空间局部估计或空间局部插值，它是建立在变异函数理论及结构分析基础上的，在有限区域内对区域化变量的取值进行无偏最优估计的一种方法。协同克里金(Co-Kring，CK)是克里金差值的扩展，在克里金差值的基础上，协同克里金还结合了相关变量，将不同变量间的交叉相关性融入其中。
现有技术的缺点：
(1)在结构复杂的森林中，在使用单一光学遥感图像时，会遇到光谱饱和，对AGB较大的区域造成估计误差；
(2)多元线性回归对原始采样数据的统计分布有正态分布的要求，且由于森林AGB与遥感特征之间存在复杂的非线性关系，这使得在实际的AGB估算中应用此方法有限；
(3)大面积的森林中AGB取值范围较大，用有限的地面观测数据训练的随机森林模型来预测AGB容易产生过度拟合现象；
(4)单一的随机森林模型进行AGB制图时忽略了样本之间的空间自相关关系，不能充分描述AGB空间分布局部特征。
发明内容
针对现有技术存在的上述问题，本发明所要解决的技术问题在于提供一种结合随机森林和协同克里金的生物量遥感制图方法。
为了解决上述技术问题，本发明所采用的技术方案如下：
一种结合随机森林和协同克里金的生物量遥感制图方法，包括如下步骤：
(1)获取样地点在多源遥感影像中对应的各个光谱信息、纹理信息、地形因子，作为随机森林预测AGB的备选变量；
(2)获取样地点的AGB实测值，并划分为训练集和验证集；使用随机森林对备选变量进行重要性分析并筛选出重要变量，以重要变量和训练集进行随机森林初步反演AGB，得到RF模型预测的AGB值；
(3)将每个样地点RF模型预测的AGB值与AGB实测值的残差值作为主变量，将高程作为协变量，进行协同克里金插值，模拟得到整个区域随机森林预测结果的CK残差，将RF模型预测的AGB值与CK残差进行相加运算，即得RFCK模型预测的AGB；
(4)采用验证集验证RFCK模型预测AGB的精度；
(5)制作RFCK模型预测AGB的空间分布图。
进一步地，步骤(1)具体为：
将Landsat的地表反射率数据进行主成分分析、缨帽变换、光谱特征变换、植被指数计算和纹理参数提取，提取出5组特征变量：原始地表反射率波段变量、光谱组合变量、植被指数变量、纹理变量和缨帽变换变量；
用DEM数据提取三个地形因子特征变量：高程、坡度和坡向；
对PALSAR的HH、HV极化数据进行组合计算，提取2组特征变量：向后散射变量和纹理变量；
将以上特征变量作为随机森林模型预测AGB的备选变量；将以上所有特征变量对应的特征层图像堆栈后，将样地点的坐标与堆栈的图像的像素点坐标进行匹配，使每个样地点的AGB对应于一组特征变量的数值，作为随机森林预测AGB的备选变量。
进一步地，步骤(2)中所述使用随机森林对备选变量进行重要性分析并筛选出重要变量的方法具体为：利用随机森林进行将备选变量随机替换，然后评价替换前后预测器的预测正确率变化，正确率大的备选变量在AGB预测中起到的作用大，属于重要变量。
进一步地，步骤(2)中所述使用随机森林对备选变量进行重要性分析并筛选出重要变量的方法具体为：导出每个备选变量的％IncMSE和IncNodePurity两个重要性因子，选择％IncMSE排前25位的备选变量，再结合其IncNodePurity进行综合制图排序展示，目视筛选，选择同时满足排序较高的变量作为RF建模变量。
进一步地，步骤(3)具体包括如下步骤：
1)根据样地点的AGB实测值及RF模型预测的AGB值，计算得到样点的残差值：
其中Z(xi)代表i点的AGB残差值，C(xi)是i点实测得到的AGB值，
是i点RF模型预测的AGB值；
2)将残差作为主变量，高程为协变量，对残差值进行半方差函数分析，选择最优函数模型参数，并基于最优函数模型参数对AGB的残差值进行协同克里金插值，得到全区域的残差预测结果：
其中ZCK *(x0)是任意一位置待估计的AGB残差；Z1(x1i)是i点的AGB残差值；λ1i是i点的权重系数，Z2(x2j)是j点的高程，λ2j是协变量高程的权重系数，N1是主变量训练样本的数量，N2是协变量高程的采样点数，要求N1≥N2；
3)将RF模型预测的AGB预测值与基于协同克里金法的残差估计值进行相加运算，即得到RFCK模型AGB预测结果：
其中，
是用RFCK模型得到的该像元位置AGB值，
是该点通过CK插值得到的AGB估计残差。
进一步地，步骤(4)具体为：将验证集中AGB实测值与RFCK模型预测的AGB值对比，用R2、均方根误差、平均绝对误差和偏差评价RFCK模型预测AGB的精度。
进一步地，步骤(5)具体为：将RF模型预测的AGB值的空间分布预测图与CK的残差插值的结果图进行图像加运算，再对非森林区域进行掩膜，即得RFCK模型预测AGB的空间分布图。
相比于现有技术，本发明的有益效果为：
(1)本发明融合了光学、雷达等多源遥感数据来提取备选变量，将多源遥感数据中衍生出的诸多信息因子进行重要性综合排序作图，重要的变量用于建模。减少了单一光学遥感数据在估算AGB时出现的光谱饱和现象，将多源遥感数据的优势有效发挥；
(2)本发明将残差进行协同克里金插值用于修正初步预测的AGB，与只考虑遥感数据变量的RF模型相比，RFCK模型的AGB预测精度更高；
(3)本发明将高程作为协同克里金插值的协变量，充分考虑了残差与AGB在地形复杂的大区域中空间分布的关系；
(4)本发明提出了结合地统计学方法和机器学习方法的优势进行AGB空间制图的思想，考虑了数据之间的空间自相关，将协同克里金对训练样本残差的插值结果纳入AGB估算及制图策略中，减少了使用单一随机森林模型时出现的过拟合现象，特别是AGB较高的区域。
附图说明
图1为本发明结合了随机森林和协同克里金的AGB估算方法及制图的流程图；
图2为建模变量选择时的重要性综合排序图；
图3为协同克里金对训练样本点残差插值的结果图；
图4为RF与RFCK模型在验证样本中实测AGB与预测AGB的散点图；
图5为粤北地区(p122/r43)森林AGB空间分布图；
图6为一样地点小块区域的真实地物影像和预测放大展示对比图；(a)粤北一块小区域山地位置；(b)该区域的真彩色影像；(c)该区域DEM数据；(d)单一RF模型的AGB制图；(e)该区域残差的CK插值图；(f)通过(d)+(e)得到的RFCK模型的AGB制图。
具体实施方式
下面结合具体实施例对本发明进一步进行描述。
实施例1：
本发明在现有AGB估算模型的基础上，发展了一种如图1所示的结合随机森林和协同克里金的生物量遥感制图的方法，即RFCK(Random forest-CoKriging)模型。该模型充分考虑到AGB在地理空间位置上的相关性和光学遥感传感器中的光谱饱和的影响，采用了协同克里金法对随机森林预测中的残差进行空间插值，在一定程度上修正随机森林预测下的过拟合，提高区域AGB估算和制图的精度和效率。具体包括如下步骤：
1、数据获取和预处理
本实施例采用从美国地质调查局(USGS)网站下载的2010年9月26日云量为11％的Landsat TM的陆地卫星影像数据，下载的数据轨道号为122/43，本研究利用NASA开发的一种图像预处理算法LEDAPS通过几何校正、辐射定标、大气校正、云检测，获得地表反射率数据，并对地表反射率数据进行地形校正；同时下载NASA研制的航天飞机雷达地形测绘任务(SRTM)在分辨率为30米数字高程模型(DEM)数据；下载得到的ALOS/PALSAR数据由日本宇宙航空研究开发机构(JAXA)进行预处理，再根据公式(1)利用ENVI中的波段计算器将其代表SAR信号振幅的DN值转换为HH和HV的后向散射系数。
σ0(dB)＝10log10(DN2)+CF#(1)
式1中，DN为HH和HV后向散射振幅的数值(表示为无符号短整型)，CF是绝对校准因子，设置为-83。并对通过该式计算得到的HH和HV后向散射图像使用7×7Lee滤波器来减少散斑噪声。所有预处理后遥感影像参照Landsat数据的轨道号的边界进行剪裁。
采用广东省森林资源监测中心提供的2007年和2012年国家森林资源连续清查数据在Landsat影像122/43轨道号内的样方，由广东省森林资源监测中心发布的树种特异性异速生长方程计算出各个样地的AGB。将2007年和2012年对应的样地AGB进行线性插值得到2010年样地实测AGB。再对样地进行筛选，采用的方法为将样地叠加在遥感影像上目视解译，手动删除被云层覆盖的样地，再将初次筛选的样地AGB值进行统计分析，采用三倍标准差法剔除作为离群值的样地，得到最终的样地数据。
将Landsat的地表反射率进行主成分分析、缨帽变换、植被指数计算、纹理指数提取；用DEM数据提取坡度和坡向；对PALSAR的HH、HV极化数据进行一系列组合计算。
(1)主成分变换
本实施例对经过预处理的Landsat影像的原始波段在ENVI软件中进行主成分变换，生成PC1、PC2、PC3三个分量，各期的三个主成分均携带了原图像90％以上的光谱信息，其中PC1包含了80％以上的光谱信息。
(2)缨帽变换
本实施例在ENVI软件中对Landsat影像进行缨帽变换，将生成的三个波段命名为TCB(Brightness)，TCG(Greenness)和TCW(Wetness)。同时对缨帽变换结果的三个波段进行了组合计算得到TCA(Tasseled Cap distance)和TCD(Tasseled Cap disturbanceindex)。
(3)光谱特征变换
除了Landsat影像的6个原始波段外(去除了热红外波段)，我们还提取波段反射率的倒数，和以多种组合计算得到的多个单波段作为森林生物量反演的备选变量，共计24个变量：1/TM1，1/TM2，1/TM3，1/TM4，1/TM5，1/TM7，TM12，TM13，TM14，TM15，TM17，TM23，TM24，TM25，TM27，TM34，TM35，TM37，TM45，TM47，TM57，TM2357，TM357，TM3574。
因为HH极化数据对森林冠层下的水敏感，而HV极化数据对森林体积散射和生物量敏感，对ALOS/PALSAR预处理后用公式转化得到的HH、HV极化下的后向散射波段计算出一些组合波段，其中有两者相加(HH+HV)和两者比值(HH/HV)，以及雷达森林退化指数(RFDI)。
(4)植被指数
植被指数是利用卫星不同波段探测数据组合而成的，能反映植物生长状况的指数。参照已有研究并结合每种植被指数的特性，本实施例计算了NDVI、DVI、SAVI、SLAVI、TVI、EVI、ARI1、ARVI、RSI、VRI、RDVI、NRI、SIPI、OSAVI共14种植被指数作为备选变量。
(5)纹理参数
纹理是遥感影像的重要特征，它不仅反映地表的粗糙程度，也揭示了影像中地物的结构信息及其与周围环境的关系，折射出地表覆盖类型空间变化的重要信息。本文对第一主成分采(PC1)用8个纹理测度，滞后一个距离，进行纹理提取，窗口大小依次设置为3×3，5×5，7×7，9×9。
(6)地形因子
林木的生长往往会受到海拔、坡度、坡向等的影响，因此对实施例所在研究区的DEM数据进行地理因子的提取。由于Landsat影像为30米分辨率，为保持空间分辨率的一致性，首先对90米分辨率的SRTM DEM数据进行重采样，再在ArcGIS软件中提取坡度、坡向。
表1展示了上述所有从遥感数据中得到的特征变量。将以上遥感数据诸多可用的特征变量所对应的特征图像在ENVI中进行堆栈，再将样方数据中样地点的经纬度坐标导入，叠加在堆栈后的图像上，提取出每个点对应的各个遥感特征光谱特征值并导出，得到每个样本点的AGB值对应于一组特征变量的数值。将所有样地采用分层抽样的不放回抽样分为两组分别作为建模数据和验证数据，其中建模数据占总数的80％，验证数据占总数的20％。
表1多源遥感数据中提取的AGB建模备选变量合集
2、变量选择和随机森林建模
在R软件中使用randomForest包进行变量的重要性分析和AGB初步预测。利用随机森林进行特征选择的基本原理是将待考察属性随机替换，然后评价替换前后预测器的预测正确率变化，从而对改变量进行评分，依据评分，可以对变量进行排序从而判断在预测中起到的作用的大小，并筛选出相对重要的变量。其中评价变量重要性的两个因子分别为％IncMSE和IncNodePurity，导出每个变量的两个重要性因子先进行粗略排序，选择％IncMSE排前25位的变量，再结合其IncNodePurity进行综合制图展示，如图2所示将两个重要性因子表示在一张图上以便目视筛选，选择同时满足排序较高且图标圆圈大偏浅色的变量作为RF建模变量。
筛选出10个综合重要性排序最高的变量作为RF建模的自变量，使用划分出的训练数据建立新的随机森林模型，调整合适的模型参数：mtry＝3、ntree＝500和Nodesize＝5，预测得到每个样地初步的AGB预测值。同时将建模自变量对应的10个图层重新堆栈，用同样的RF模型得到初步的区域AGB空间分布制图。
3、残差的协同克里金插值
复杂地形区域的森林AGB会受到高程变化的影响，即高山地区森林茂密生物量偏高、平原地区城市聚集而森林稀薄生物量偏低的现象。结合了CK和RF从而得到最终AGB预测值的步骤如下：
(1)根据样地数据中AGB实测值及RF模型预测值，按式2计算得到样点的残差值：
其中Z(xi)代表i点的AGB残差值，C(xi)是i点实测得到的AGB值，
是i点RF模型预测得到的AGB值。
(2)将残差作为主变量，高程为协变量，利用GS+9.0软件对残差值进行半方差函数分析，选择最优函数模型参数，并基于最优函数模型参数在ArcGIS软件的地统计工具中对AGB的残差值进行协同克里金插值(式3)，得到如图3全区域的残差预测结果：
其中ZCK *(x0)是任意一位置待估计的AGB残差；Z1(x1i)是i点的AGB残差值；λ1i是i点的权重系数，Z2(x2j)是j点的高程，λ2j是协变量高程的权重系数，N1是主变量训练样本的数量，N2是协变量高程的采样点数，要求N1≥N2。
(3)将基于RF模型的AGB预测值与基于协同克里金法的残差估计值进行相加运算(式4)，即得到基于RFCK模型的AGB预测结果：
其中，
是用RFCK模型得到的该像元位置AGB值，
是该点通过CK插值得到的AGB估计残差。
4、AGB预测的精度验证及空间制图
将总样本中分出的20％训练样本数据实测值，分别与经过残差相加前后的RF初步预测值对比，用R2、均方根误差(RMSE)、平均绝对误差(MAE)和偏差(Bias)来评价模型精度，如表2所示。从验证精度上，结合了协同克里金插值修正残差后的RFCK模型各项指标表现都优于单一的RF模型，其中R2从0.49提升到0.55，RMSE相比提高了3.8％，偏差绝对值更接近0，MAE也更低。将RF模型和RFCK模型的预测值对应实测值生成散点图如图4，RFCK的蓝色散点相比RF的橙色散点更靠近1∶1的标准趋势线，即预测值更接近实测值，AGB＞100t/ha的几个样本点更加明显，由此可以证实加入CK插值对残差进行一部分修正，可以一定程度上减少RF模型在森林茂密的AGB值较高区域进行预测时发生的过拟合现象。
表2验证样本中随机森林模型在加入残差插值前后的精度对比
| 模型 | R<sup>2</sup> | MAE(t/ha) | RMSE(t/ha) | Bias(t/ha) |
| RF | 0.49 | 24.56 | 31.21 | 2.09 |
| RFCK | 0.55 | 22.14 | 29.97 | 1.07 |
运用ArcGIS软件中的栅格计算功能(raster calculator)将步骤“1、数据获取和预处理”中基于RF模型的AGB初步空间分布预测图与步骤“3、残差的协同克里金插值”中CK的残差插值结果图(图3)进行图像加运算，再通过对Landsat影像在ENVI软件的监督分类结果对非森林区域进行掩膜，即得到本实例区域基于RFCK模型预测AGB的空间分布图，如图5所示。
图6为其中一小块区域的放大展示，通过目视解读可以发现，RFCK的AGB制图结果可以很好得与真实地物影像(图6a)中山谷与植被覆盖茂密地区的纹理对应，同时对应了(图6c)图中山脚和平原人类活动多的地带AGB低，山顶人类活动较少植被茂密的地区AGB高。
本发明在现有机器学习模型基础上，发展了一种适合于复杂地形区域反演生物量的RFCK模型。采用地统计分析中的克里金插值法后，可以将样地点生物量真实测量值与随机森林预测值的残差外推到整个研究区，用于修正使用单一随机森林模型预测产生的过拟合情况。普通克里金法对残差进行插值时没有考虑到研究区的地形因素，即高山地区森林茂密生物量偏高、平原地区城市聚集森林稀薄生物量偏低的情况。本发明将研究区的高程作为插值的协变量，采用了协同克里金对随机森林产生的残差进行地统计插值，从而修正一部分预测生物量时的高值低估、低值高估情况，提高了生物量反演的精度和生物量空间制图的可靠性。
Claims (7)
1.一种结合随机森林和协同克里金的生物量遥感制图方法，其特征在于，包括如下步骤：
(1)获取样地点在多源遥感影像中对应的各个光谱信息、纹理信息、地形因子，作为随机森林预测AGB的备选变量；
(2)获取样地点的AGB实测值，并划分为训练集和验证集；使用随机森林对备选变量进行重要性分析并筛选出重要变量，以重要变量和训练集进行随机森林初步反演AGB，得到RF模型预测的AGB值；
(3)将每个样地点RF模型预测的AGB值与AGB实测值的残差值作为主变量，将高程作为协变量，进行协同克里金插值，模拟得到整个区域随机森林预测结果的CK残差，将RF模型预测的AGB值与CK残差进行相加运算，即得RFCK模型预测的AGB；
(4)采用验证集验证RFCK模型预测AGB的精度；
(5)制作RFCK模型预测AGB的空间分布图。
2.根据权利要求1所述结合随机森林和协同克里金的基于多源遥感影像生物量制图新方法，其特征在于，步骤(1)具体为：
将Landsat的地表反射率数据进行主成分分析、缨帽变换、光谱特征变换、植被指数计算和纹理参数提取，提取出5组特征变量：原始地表反射率波段变量、光谱组合变量、植被指数变量、纹理变量和缨帽变换变量；
用DEM数据提取三个地形因子特征变量：高程、坡度和坡向；
对PALSAR的HH、HV极化数据进行组合计算，提取2组特征变量：向后散射变量和纹理变量；
将以上特征变量作为随机森林模型预测AGB的备选变量；将以上所有特征变量对应的特征层图像堆栈后，将样地点的坐标与堆栈的图像的像素点坐标进行匹配，使每个样地点的AGB对应于一组特征变量的数值，作为随机森林预测AGB的备选变量。
3.根据权利要求1所述结合随机森林和协同克里金的基于多源遥感影像生物量制图新方法，其特征在于，步骤(2)中所述使用随机森林对备选变量进行重要性分析并筛选出重要变量的方法具体为：利用随机森林进行将备选变量随机替换，然后评价替换前后预测器的预测正确率变化，正确率大的备选变量在AGB预测中起到的作用大，属于重要变量。
4.根据权利要求3所述结合随机森林和协同克里金的基于多源遥感影像生物量制图新方法，其特征在于，步骤(2)中所述使用随机森林对备选变量进行重要性分析并筛选出重要变量的方法具体为：导出每个备选变量的％IncMSE和IncNodePurity两个重要性因子，选择％IncMSE排前25位的备选变量，再结合其IncNodePurity进行综合制图排序展示，目视筛选，选择同时满足两个重要性排序都较高的变量作为RF建模变量。
5.根据权利要求1所述结合随机森林和协同克里金的基于多源遥感影像生物量制图新方法，其特征在于，步骤(3)具体包括如下步骤：
1)根据样地点的AGB实测值及RF模型预测的AGB值，计算得到样点的残差值：
其中Z(xi)代表i点的AGB残差值，C(xi)是i点实测得到的AGB值，
是i点RF模型预测的AGB值；
2)将残差作为主变量，高程为协变量，对残差值进行半方差函数分析，选择最优函数模型参数，并基于最优函数模型参数对AGB的残差值进行协同克里金插值，得到全区域的残差预测结果：
其中ZCK *(x0)是任意一位置待估计的AGB残差；Z1(x1i)是i点的AGB残差值；λ1i是i点的权重系数，Z2(x2j)是j点的高程，λ2j是协变量高程的权重系数，N1是主变量训练样本的数量，N2是协变量高程的采样点数，要求N1≥N2；
3)将RF模型预测的AGB预测值与基于协同克里金法的残差估计值进行相加运算，即得到RFCK模型AGB预测结果：
其中，
是用RFCK模型得到的该像元位置AGB值，
是该点通过CK插值得到的AGB估计残差。
6.根据权利要求1所述结合随机森林和协同克里金的基于多源遥感影像生物量制图新方法，其特征在于，步骤(4)具体为：将验证集中AGB实测值与RFCK模型预测的AGB值对比，用R2、均方根误差、平均绝对误差和偏差评价RFCK模型预测AGB的精度。
7.根据权利要求1所述结合随机森林和协同克里金的基于多源遥感影像生物量制图新方法，其特征在于，步骤(5)具体为：将RF模型预测的AGB值的空间分布预测图与CK的残差插值的结果图进行图像加运算，再对非森林区域进行掩膜，即得RFCK模型预测AGB的空间分布图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011299207.5A CN112395808A (zh) | 2020-11-18 | 2020-11-18 | 一种结合随机森林和协同克里金的生物量遥感制图方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011299207.5A CN112395808A (zh) | 2020-11-18 | 2020-11-18 | 一种结合随机森林和协同克里金的生物量遥感制图方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112395808A true CN112395808A (zh) | 2021-02-23 |
Family
ID=74607516
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011299207.5A Pending CN112395808A (zh) | 2020-11-18 | 2020-11-18 | 一种结合随机森林和协同克里金的生物量遥感制图方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112395808A (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113205565A (zh) * | 2021-04-13 | 2021-08-03 | 武汉大学 | 基于多源数据的森林生物量遥感制图方法 |
| CN113569488A (zh) * | 2021-08-04 | 2021-10-29 | 中国科学院地理科学与资源研究所 | 一种基于随机森林回归的体感温度预测方法及系统 |
| CN113609940A (zh) * | 2021-07-26 | 2021-11-05 | 南通智能感知研究院 | 一种广覆盖的土壤重金属遥感混合制图方法 |
| CN113902580A (zh) * | 2021-10-18 | 2022-01-07 | 四川农业大学 | 一种基于随机森林模型的历史耕地分布重建方法 |
| CN114694036A (zh) * | 2022-03-18 | 2022-07-01 | 南京农业大学 | 一种基于高分影像和机器学习的高海拔地区作物分类识别的方法 |
| CN114722646A (zh) * | 2022-06-10 | 2022-07-08 | 西安交通大学 | 基于克里金模型的自给能探测器三维测点布置优化方法 |
| CN114780599A (zh) * | 2022-04-06 | 2022-07-22 | 四川农业大学 | 基于小麦品比试验数据的综合分析系统 |
| CN114821349A (zh) * | 2022-03-01 | 2022-07-29 | 南京林业大学 | 顾及谐波模型系数和物候参数的森林生物量估算方法 |
| CN114937029A (zh) * | 2022-06-21 | 2022-08-23 | 西南林业大学 | 森林碳储量抽样估测方法、装置、设备及存储介质 |
| CN115859032A (zh) * | 2022-12-16 | 2023-03-28 | 西北农林科技大学 | 一种逐日气象要素自适应插值方法 |
| CN117290910A (zh) * | 2023-09-13 | 2023-12-26 | 爬山虎科技股份有限公司 | 一种基于空间插值算法的土壤图件自动化编制方法 |
Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107247809A (zh) * | 2017-07-19 | 2017-10-13 | 南京林业大学 | 一种人工林森林年龄空间制图的新方法 |
| US20180096146A1 (en) * | 2015-11-18 | 2018-04-05 | Tencent Technology (Shenzhen) Company Limited | Method and apparatus for identifying malicious software |
| WO2018098190A1 (en) * | 2016-11-28 | 2018-05-31 | The Climate Corporation | Determining intra-field yield variation data based on soil characteristics data and satellite images |
| WO2018156778A1 (en) * | 2017-02-22 | 2018-08-30 | The United States Of America, As Represented By The Secretary, Department Of Health And Human Services | Detection of prostate cancer in multi-parametric mri using random forest with instance weighting & mr prostate segmentation by deep learning with holistically-nested networks |
| CN108876917A (zh) * | 2018-06-25 | 2018-11-23 | 西南林业大学 | 一种森林地上生物量遥感估测通用模型构建方法 |
| CN108921885A (zh) * | 2018-08-03 | 2018-11-30 | 南京林业大学 | 一种综合三类数据源联合反演森林地上生物量的方法 |
| CN109063657A (zh) * | 2018-08-08 | 2018-12-21 | 武汉大学 | 面向均值地域光谱单元的地上生物量估算和尺度转换方法 |
| CN109142679A (zh) * | 2018-08-13 | 2019-01-04 | 中国科学院东北地理与农业生态研究所 | 基于人工神经网络克里金插值的森林土壤养分的空间预测方法 |
| CN109213964A (zh) * | 2018-07-13 | 2019-01-15 | 中南大学 | 一种融合多源特征地理参数的卫星aod产品校正方法 |
| CN109583516A (zh) * | 2018-12-24 | 2019-04-05 | 天津珞雍空间信息研究院有限公司 | 一种基于地基和卫星观测的时空连续pm2.5反演方法 |
| CN109657363A (zh) * | 2018-12-24 | 2019-04-19 | 天津珞雍空间信息研究院有限公司 | 一种时空连续的pm2.5反演方法 |
| CN110610258A (zh) * | 2019-08-20 | 2019-12-24 | 中国地质大学（武汉） | 融合多源时空数据的城市空气质量精细化估测方法及装置 |
| WO2020139344A1 (en) * | 2018-12-27 | 2020-07-02 | Halliburton Energy Services, Inc. | Hydraulic fracturing operation planning using data-driven multi-variate statistical machine learning modeling |
-
2020
- 2020-11-18 CN CN202011299207.5A patent/CN112395808A/zh active Pending
Patent Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20180096146A1 (en) * | 2015-11-18 | 2018-04-05 | Tencent Technology (Shenzhen) Company Limited | Method and apparatus for identifying malicious software |
| WO2018098190A1 (en) * | 2016-11-28 | 2018-05-31 | The Climate Corporation | Determining intra-field yield variation data based on soil characteristics data and satellite images |
| WO2018156778A1 (en) * | 2017-02-22 | 2018-08-30 | The United States Of America, As Represented By The Secretary, Department Of Health And Human Services | Detection of prostate cancer in multi-parametric mri using random forest with instance weighting & mr prostate segmentation by deep learning with holistically-nested networks |
| CN107247809A (zh) * | 2017-07-19 | 2017-10-13 | 南京林业大学 | 一种人工林森林年龄空间制图的新方法 |
| CN108876917A (zh) * | 2018-06-25 | 2018-11-23 | 西南林业大学 | 一种森林地上生物量遥感估测通用模型构建方法 |
| CN109213964A (zh) * | 2018-07-13 | 2019-01-15 | 中南大学 | 一种融合多源特征地理参数的卫星aod产品校正方法 |
| CN108921885A (zh) * | 2018-08-03 | 2018-11-30 | 南京林业大学 | 一种综合三类数据源联合反演森林地上生物量的方法 |
| CN109063657A (zh) * | 2018-08-08 | 2018-12-21 | 武汉大学 | 面向均值地域光谱单元的地上生物量估算和尺度转换方法 |
| CN109142679A (zh) * | 2018-08-13 | 2019-01-04 | 中国科学院东北地理与农业生态研究所 | 基于人工神经网络克里金插值的森林土壤养分的空间预测方法 |
| CN109583516A (zh) * | 2018-12-24 | 2019-04-05 | 天津珞雍空间信息研究院有限公司 | 一种基于地基和卫星观测的时空连续pm2.5反演方法 |
| CN109657363A (zh) * | 2018-12-24 | 2019-04-19 | 天津珞雍空间信息研究院有限公司 | 一种时空连续的pm2.5反演方法 |
| WO2020139344A1 (en) * | 2018-12-27 | 2020-07-02 | Halliburton Energy Services, Inc. | Hydraulic fracturing operation planning using data-driven multi-variate statistical machine learning modeling |
| CN110610258A (zh) * | 2019-08-20 | 2019-12-24 | 中国地质大学（武汉） | 融合多源时空数据的城市空气质量精细化估测方法及装置 |
Non-Patent Citations (3)
| Title |
|---|
| HUIYI SU .ETAL: "Machine learning and geostatistical approaches for estimating aboveground biomass in Chinese subtropical forests", 《WWW.REAEARCHAQUARE.COM/ARTICLE/RS-25148/V3》 * |
| RUILI .ETAL: "Satellite-based prediction of daily SO2 exposure across China using a high-quality random forest-spatiotemporal Kriging (RF-STK) model for health risk assessment", 《ATMOSPHERIC ENVIRONMENT》 * |
| 陆秋琴 等: "基于网格划分的空间关联区域VOCs浓度预测研究", 《西安理工大学学报》 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113205565A (zh) * | 2021-04-13 | 2021-08-03 | 武汉大学 | 基于多源数据的森林生物量遥感制图方法 |
| CN113609940A (zh) * | 2021-07-26 | 2021-11-05 | 南通智能感知研究院 | 一种广覆盖的土壤重金属遥感混合制图方法 |
| CN113569488A (zh) * | 2021-08-04 | 2021-10-29 | 中国科学院地理科学与资源研究所 | 一种基于随机森林回归的体感温度预测方法及系统 |
| CN113902580B (zh) * | 2021-10-18 | 2023-04-07 | 四川农业大学 | 一种基于随机森林模型的历史耕地分布重建方法 |
| CN113902580A (zh) * | 2021-10-18 | 2022-01-07 | 四川农业大学 | 一种基于随机森林模型的历史耕地分布重建方法 |
| CN114821349B (zh) * | 2022-03-01 | 2023-07-07 | 南京林业大学 | 顾及谐波模型系数和物候参数的森林生物量估算方法 |
| CN114821349A (zh) * | 2022-03-01 | 2022-07-29 | 南京林业大学 | 顾及谐波模型系数和物候参数的森林生物量估算方法 |
| CN114694036A (zh) * | 2022-03-18 | 2022-07-01 | 南京农业大学 | 一种基于高分影像和机器学习的高海拔地区作物分类识别的方法 |
| CN114780599A (zh) * | 2022-04-06 | 2022-07-22 | 四川农业大学 | 基于小麦品比试验数据的综合分析系统 |
| CN114722646B (zh) * | 2022-06-10 | 2022-08-26 | 西安交通大学 | 基于克里金模型的自给能探测器三维测点布置优化方法 |
| CN114722646A (zh) * | 2022-06-10 | 2022-07-08 | 西安交通大学 | 基于克里金模型的自给能探测器三维测点布置优化方法 |
| CN114937029A (zh) * | 2022-06-21 | 2022-08-23 | 西南林业大学 | 森林碳储量抽样估测方法、装置、设备及存储介质 |
| CN115859032A (zh) * | 2022-12-16 | 2023-03-28 | 西北农林科技大学 | 一种逐日气象要素自适应插值方法 |
| CN115859032B (zh) * | 2022-12-16 | 2023-09-22 | 西北农林科技大学 | 一种逐日气象要素自适应插值方法 |
| CN117290910A (zh) * | 2023-09-13 | 2023-12-26 | 爬山虎科技股份有限公司 | 一种基于空间插值算法的土壤图件自动化编制方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112395808A (zh) | 一种结合随机森林和协同克里金的生物量遥感制图方法 | |
| Gkikas et al. | ModIs Dust AeroSol (MIDAS): a global fine-resolution dust optical depth data set | |
| Lu et al. | A survey of remote sensing-based aboveground biomass estimation methods in forest ecosystems | |
| Villoslada et al. | Fine scale plant community assessment in coastal meadows using UAV based multispectral data | |
| Vaudour et al. | Temporal mosaicking approaches of Sentinel-2 images for extending topsoil organic carbon content mapping in croplands | |
| Cheng et al. | Detecting diurnal and seasonal variation in canopy water content of nut tree orchards from airborne imaging spectroscopy data using continuous wavelet analysis | |
| Nelson et al. | Estimating Siberian timber volume using MODIS and ICESat/GLAS | |
| Roberts et al. | Relationships between dominant plant species, fractional cover and land surface temperature in a Mediterranean ecosystem | |
| Cutler et al. | Estimating tropical forest biomass with a combination of SAR image texture and Landsat TM data: An assessment of predictions between regions | |
| Luo et al. | Retrieving aboveground biomass of wetland Phragmites australis (common reed) using a combination of airborne discrete-return LiDAR and hyperspectral data | |
| Lin et al. | Monitoring sugarcane growth using ENVISAT ASAR data | |
| Olmedo et al. | Nine years of SMOS sea surface salinity global maps at the Barcelona Expert Center | |
| CN114460013B (zh) | 滨海湿地植被地上生物量gan模型自学习遥感反演方法 | |
| CN110687053B (zh) | 一种基于高光谱影像的区域有机质含量估算方法和装置 | |
| Skowronek et al. | Transferability of species distribution models for the detection of an invasive alien bryophyte using imaging spectroscopy data | |
| Hojatimalekshah et al. | Tree canopy and snow depth relationships at fine scales with terrestrial laser scanning | |
| Bayat et al. | Retrieval of land surface properties from an annual time series of Landsat TOA radiances during a drought episode using coupled radiative transfer models | |
| CN113205014B (zh) | 一种基于图像锐化的时序数据耕地提取方法 | |
| CN113466143B (zh) | 土壤养分反演方法、装置、设备及介质 | |
| Ni et al. | Retrieval of forest biomass from ALOS PALSAR data using a lookup table method | |
| GB2620469A (en) | Spatial prediction and evaluation method of soil organic matter content based on partition algorithm | |
| Quiñones et al. | Exploration of factors limiting biomass estimation by polarimetric radar in tropical forests | |
| Bolun et al. | Estimating rice paddy areas in China using multi-temporal cloud-free normalized difference vegetation index (NDVI) imagery based on change detection | |
| Yang et al. | Forest canopy height mapping over China using GLAS and MODIS data | |
| Khudhur et al. | Comparison of the accuracies of different spectral indices for mapping the vegetation covers in Al-Hawija district, Iraq |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20210223 |