CN114821349A - 顾及谐波模型系数和物候参数的森林生物量估算方法 - Google Patents

顾及谐波模型系数和物候参数的森林生物量估算方法 Download PDF

Info

Publication number
CN114821349A
CN114821349A CN202210195246.3A CN202210195246A CN114821349A CN 114821349 A CN114821349 A CN 114821349A CN 202210195246 A CN202210195246 A CN 202210195246A CN 114821349 A CN114821349 A CN 114821349A
Authority
CN
China
Prior art keywords
model
forest
phenological
value
parameters
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.)
Granted
Application number
CN202210195246.3A
Other languages
English (en)
Other versions
CN114821349B (zh
Inventor
李明诗
张亚丽
蒋路凡
刘琴琴
王楠
杨博翔
叶鋆泓
彭钰雯
刘嘉薇
张银
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Forestry University
Original Assignee
Nanjing Forestry University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Nanjing Forestry University filed Critical Nanjing Forestry University
Priority to CN202210195246.3A priority Critical patent/CN114821349B/zh
Publication of CN114821349A publication Critical patent/CN114821349A/zh
Application granted granted Critical
Publication of CN114821349B publication Critical patent/CN114821349B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • 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
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G06N20/10—Machine learning using kernel methods, e.g. support vector machines [SVM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G06N20/20—Ensemble learning
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Mathematical Physics (AREA)
  • Medical Informatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开一种顾及谐波模型系数和物候参数的森林生物量估算方法，包括如下步骤：S1、影像获取及预处理；S2、建立基于多源遥感数据的MCCDC模型；S3、基于逻辑回归方程的森林年物候提取；S4、基于物候参数的森林生物量制图。本发明将基于遥感提取的森林物候参数和所构建的谐波模型MCCDC系数作为树种分类和生物量建模的输入变量，在提高森林树种组制图精度基础上，耦合光学遥感和雷达影像，分树种构建光学‑雷达‑物候模型，一定程度上弥补了传统AGB建模时不区分树种对建模精度的影响的缺陷，从而提高区域尺度森林生物量制图的精度。

Description

顾及谐波模型系数和物候参数的森林生物量估算方法
技术领域
本发明属于林业遥感技术领域，具体是涉及一种顾及多源遥感数据的不同地物长时 间序列反射率谐波模型模拟系数变化特征，以及据此提取的植被物候参数执行森林生物 量制图的新方法。
背景技术
森林生物量(AGB,Aboveground Biomass)是森林生态系统碳汇潜力评估的重要基础，其变化也反映了森林对生态演替、自然干扰、人类活动和气候变化等因素的响应， 是森林生态系统健康评价的重要指标。因此，快速准确地获取森林生物量信息对于森林 经营管理方案有效性以及森林的碳中和能力评价具有基础性支撑作用，也有助于加深对 森林生态系统碳循环和碳动态特征的认知。传统的森林生物量研究主要依赖于森林清查 方式进行。这种方式的主要缺陷是：人力物力投入大，清查结果的时效性较差，受地形 限制多以及难以大面积推广等。
针对传统清查方式的不足，当前主流研究方向是将地面观测数据与遥感卫星图像光 谱或散射、纹理及地形特征相结合进行统计建模来实现森林生物量估算或制图，目前已有一些关于森林生物量估算的研究，包括基于光学遥感数据(如AVHRR、MODIS、 Landsat、IKONOS、Quickbird)，基于雷达数据(如PALSAR数据)，采用非参数估算 模型(如RF、SVM、SGB模型)和参数模型如多元线性回归等多种建模方法的森林生 物量估算研究。其中基于光学遥感数据的参数与AGB具有较高的相关性，但因面临云 覆盖与信息饱和等问题，往往限制了其数据的可用性；雷达数据能够在一定程度上穿透 森林冠层，提供相对详细的三维植被结构，可以提高不同森林系统中植被结构参数甚至 生物量的估计精度，但在结构复杂的森林类型中，其后向散射系数或极化参数依然存在 饱和问题。基于非参数建模方法虽能很好地解决非线性和高维数等问题，但在实际应用 中存在过度拟合训练数据的风险且物理机理不明确等问题。
在这些现有的森林生物量估算研究中，综合考虑物候对AGB的影响研究较少。另外，当前大多数AGB建模往往忽略不同树种对建模精度的响应差异，将研究区的所有 树种整合进一个模型中，并且都是分单一时间点进行单独建模(未考虑时间融合策略)， 这会降低模型的估测能力和精度，并限制模型在时间尺度上外推的能力发挥。
发明内容
发明目的：本发明目的在于针对现有技术的不足，提供一种在输入特征方面考虑将 改进的连续变化检测与分类谐波模型(modified continuous change detection andclassification，MCCDC)系数，以及基于谐波模型所刻画的物候参数区域森林生物量的 区域森林生物量估算方法。
技术方案：本发明所述顾及谐波模型系数和物候参数的森林生物量估算方法，包括 如下步骤：
S1、影像获取及预处理；获取待分类区域一段时间内多个传感器的陆地卫星影像，并对获得的影像进行预处理；同时，采用最小二乘法拟合斜率和截距的方法，进行逐波 段的光谱归一化，最小化不同传感器间的光谱响应差异；
S2、建立基于多源遥感数据的MCCDC模型，具体步骤包括：
S21、建立SR时间序列模型，定义一个逐像元反射率时间序列模型，如式(1)所 示，以表达地表反射率变化的特征；
Figure BDA0003526955040000021
式(1)中包含一个傅立叶谐波模型和一个长期的趋势部分，x为儒略日期；T是周期天数，该值为365.25；
Figure BDA0003526955040000022
为第i波段在儒略日期x的地表反射率预测值；该时间序 列模型中有8个系数，系数a0,i为第i波段清晰观测的总体均值，系数a1,i、a2,i、b1,i、b2,i、 c1,i、c2,i用于模拟由物候引起的地表反射率的年内变化，它的最小正周期为一年，系数d1,i体现由气候变化和土壤退化等因素引起的地表反射率渐变的趋势；
S22、识别时间序列模型突变点，基于步骤S21建立的时间序列模型，看是否满足式(2)；
Figure BDA0003526955040000023
式(2)中，k为波段数量；n指需要判断的突变点后连续观测值的个数；ρ(i，x) 为第i波段在儒略日期x地表反射率的实际观测值；RMSEi是指Landsat第i波段的模型预 测值与观测值的均方根误差；
当满足式(2)时，证明区域内发生了土地覆盖变化，对变化的像元按照时间顺序排序，首先利用同一位置的前24个清晰像元，基于式(1)，建立初始化的时间序列模 型，然后判断每个波段接下来的4个像元的实际观测值与预测值的差异是否都超出模型 均方根误差的3倍，如果没有超出范围，则将这些观测值加入到之前参与建模的观测值 中去，重新拟合参数，再次判断后续的观测值，如果连续4个观测值都超出了范围，则 定义为突变点，即土地覆盖变化发生的时间；
S23、合成每日遥感图像，利用步骤S21～S22的方法生成每个像元不同波段的时间序列模型，对于每个像元的单一波段，通过给定变量儒略日期的值，依据时间序列模型， 生成该像元给定日期的SR预测值，将该方法应用到同一时间内所有像素点，从而生成 空间分辨率为30m的清晰合成图像；
S3、基于逻辑回归方程的森林年物候提取，具体方法为：
S31、基于获得的全部每天合成图像，进行逐年分析提取每年的物候参数，通过式(3)～式(5)依次计算三个植被指数EVI、NDVI以及LSWI获得不同的物候参数估 计值；
Figure BDA0003526955040000031
Figure BDA0003526955040000032
Figure BDA0003526955040000033
式中，ρBlue、ρRed、ρNIR、和ρswir1分别是蓝光，红光，近红外和短波红外1波段的地 表反射率值；
S32、基于三个植被指数，使用逻辑回归函数分别提取各自表达的年SOS，EOS和LOS，使用式(6)对每年卫星观测值导出的植被指数数据进行建模；
Figure BDA0003526955040000034
式中，t为观测时间，y(t)为植被指数值EVI、NDVI或者LSWI，a和b为拟合参 数，c+d为最大植被指数值，d为初始背景植被指数值；
通过式(7)计算拟合模型的曲率变化率；
Figure BDA0003526955040000041
式中，y(t)'和y(t)”分别是逻辑斯蒂克方程的一阶和二阶导数，ROC是拟合模型的曲率变化率；
当曲率的变化率首次达到极大值时对应的天数即为植被生长开始期SOS，当曲率的 变化率最后一次达到极小值时的天数为植被生长结束期EOS，取EOS与SOS的差值即 是植被生长季长度LOS；
S4、顾及物候参数和拟合MCCDC模型系数的森林生物量制图，具体包括如下步 骤：
S41、选择AGB预测模型变量来源，包括合成的光学遥感特征值、具有明显季相 节律特征的模型参数和物候期以及不受云雨干扰的雷达后向散射系数；
S42、在步骤S41的预测模型变量来源中筛选出重要变量，先使用随机森林分析包的特征重要性命令，通过均方差百分比增量和节点纯度增量两个测度对所有变量的重要性进行排序，剔除不重要的变量，均方差百分比增量和节点纯度增量值越高，对应的变 量越重要；
S43、基于光学和雷达数据的预测模型建立，对于所有来源的野外实测生物量数据， 随机将这些数据按照8:2的比例分为训练数据集和验证数据集，使用筛选后的预测变量 与三年的实地调查数据一一匹配，采用RF、SVM和SGB三种机器学习算法，基于森 林树种或树种组类型，依次建立AGB所有年份总预测模型。
本发明进一步优选地技术方案为，步骤S1中对获得的陆地卫星影像进行大气校正获得SR图像。
优选地，步骤S4中通过均方差百分比增量和节点纯度增量两个测度对所有变量的重要性进行排序，剔除不重要的变量后，最终选择16个变量，分别为Red_3，NIR， NIR_1，NIR_3，SWIR1，SWIR1_3，SWIR2_5，PCA1，PCA3，HV，HHHV，SOS，EOS，LOS，DEM和Wetness。
有益效果：本发明将MCCDC模型系数和物候期参数(SOS、EOS和LOS)用于 AGB估算，本发明中的方法并未使用中等空间分辨率遥感影像数据同更高分辨率的影 像数据相结合的分类方法，而是在仅使用中等空间分辨率遥感影像数据，利用MCCDC 模型合成每日遥感影像后，对研究区的森林树种组进行分类制图，然后分森林树种组建 模；同时本发明耦合光学遥感和雷达影像，并分树种构建光学-雷达-物候模型以提高 AGB估算精度；本发明将基于遥感提取的森林物候参数和所构建的谐波模型MCCDC 系数作为树种分类和生物量建模的输入变量，在提高森林树种组制图精度基础上，耦合 光学遥感和雷达影像，分树种构建光学-雷达-物候模型，一定程度上弥补了传统AGB 建模时不区分树种对建模精度的影响的缺陷，从而提高区域尺度森林生物量制图的精度。
本发明步骤S1中对获得的陆地卫星影像进行了逐波段的光谱归一化操作，最小化了不同传感器间光谱响应差异，从而提高了MCCDC模型拟合精度；本发明步骤S22 中利用后续4个连续观测值来定义真实的地物覆盖变化，显著降低了模型捕获虚假地物 变化的风险；本发明步骤S43，不同于现今众多的生物量建模方法，采用了分不同优势 树种(组)分别建模的策略；并且本发明还执行了时间融合(多年份数据集成在一起) 的建模策略，从而使得本发明具有时间内插和外推能力。
附图说明
图1为本发明森林生物量估算方法的流程图；
图2为本发明实施例中预测变量重要性排序图；
图3为本发明实施例中森林AGB建模性能散点图示例；
其中(a)为未分树种组建模的AGB散点图，(b)为分树种组建模后的AGB散点 图；
图4为本发明实施例中基于SVM的森林树种组时空分布图；
图5为本发明实施例中2000-2019年森林AGB空间分布图。
具体实施方式
下面通过附图对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所 述实施例。
实施例：一种顾及谐波模型系数和物候参数的森林生物量估算方法，包括如下步骤：
S1、影像获取及预处理；
S11、获取待分类区域一段时间内多个传感器的陆地卫星影像，本实施例的数据从美国地质调查局(USGS)官方网站免费下载1999年1月1日到2019年12月31日云 量小于80％的Landsat TM/ETM+/OLI的陆地卫星影像，下载的数据轨道号为131/042， 共计610景影像。TM影像为152景，ETM+影像为243景，OLI影像为119景，由于 天气原因夏季7月份的云量少于80％的影像相对较少。还从欧洲空间局下载Sentinel-2 数据，数据包括2017-2019年所有3月-11月所有云量低于80％的Sentinel-2的大气顶层 反射率(Top ofAtmosphere，TOA)辐射数据，轨道号为T47RPJ，包括2A数据58景， 2B数据38景。所选的Landsat和Sentinel-2图像基本信息如表1所示：
Figure BDA0003526955040000061
表1所用Landsat、Sentinel-2遥感影像基本信息
S12、对获得的影像进行预处理获得超分辨率图像：对于USGS提供的图像，本发 明直接申请了SR图像。其中TM和ETM+影像采用了LEDAPS算法进行大气校正从而 获得SR图像，而OLI影像采用了LaSRC算法进行大气校正从而获得SR图像。对于下 载的Sentinel-2TOA图像数据，使用Sen2Cor大气辐射传输校正代码(版本2.8，2019 年2月21日发布；http：//step.esa.int/main/third-party-plugins-2/sen2cor/)进行大气校正。
S13、为了最小化不同传感器间的光谱响应差异，本实施例中还采用最小二乘法拟合斜率和截距的方法，进行了逐波段的光谱归一化操作。
S2、建立基于多源遥感数据的MCCDC模型，具体步骤包括：
S21、建立SR时间序列模型：
基于傅立叶思想，任何复杂的周期函数均可表达成具有不同频率的正弦和余弦函数 加权和形式。本实施例首先定义一个逐像元反射率时间序列模型，如式(1)所示，以 表达地表反射率变化的特征；
Figure BDA0003526955040000071
式(1)中包含一个傅立叶谐波模型和一个长期的趋势部分，x为儒略日期；T是周期天数，该值为365.25；
Figure BDA0003526955040000072
为第i波段在儒略日期x的地表反射率预测值；该时间序 列模型中有8个系数，系数a0，i为第i波段清晰观测的总体均值，系数a1，i、a2，i、b1，i、b2，i、 c1，i、c2，i用于模拟由物候引起的地表反射率的年内变化，它的最小正周期为一年，系数d1，i体现由气候变化和土壤退化等因素引起的地表反射率渐变的趋势，这种趋势对于监测土 地的渐变非常重要；
S22、识别时间序列模型突变点，基于步骤S21建立的时间序列模型，看是否满足式(2)；
Figure BDA0003526955040000073
式(2)中，ρ(i，x)为第i波段在儒略日期x对应的地表反射率的实际观测值，k为 波段数量(本例中取值为6)，n为需要判断的突变点后连续观测值个数(本例中取值为4)，RMSEi是指第i波段的模型预测反射率值与实际遥感观测反射率值的均方根误差；
当满足式(2)时，则认为时间序列模型在整个观测期发生了土地覆盖变化。为了准确判断变化发生的具体日期，本实施例中采用模型CCDC的原理对变化的像元按照 时间顺序排序，首先利用同一位置的前24个清晰像元，基于式(1)，建立初始化的时 间序列模型，选择24个清晰观测值的原因是时间序列模型有8个系数，清晰观测值个 数是待拟合曲线系数个数的3倍能够使模型估计的系数更加准确。然后判断每个波段接 下来的4个像元的实际观测值与预测值的差异是否都超出模型均方根误差的3倍，如果 没有超出范围，则将这些观测值加入到之前参与建模的观测值中去，重新拟合参数，再 次判断后续的观测值，如果连续4个观测值都超出了范围，则定义为突变点，即土地覆 盖变化发生的时间；为了充分利用所有的光谱信息来准确定义断点，MCCDC模型的输 入是六个波段，所以需要平均所有波段的观测值和预测值的差异之和来定义“突变”。
S23、合成每日遥感图像，利用步骤S21～S23的方法生成每个像元不同波段的时间序列模型，对于每个像元的单一波段，通过给定变量儒略日期的值，依据时间序列模型， 生成该像元给定日期的SR预测值，将该方法应用到1999年至2019年每一天的所有像 素点，从而生成空间分辨率为30m的清晰合成图像。
S3、基于逻辑回归方程的森林年物候提取，具体方法为：
S31、基于获得的全部每天合成图像，进行逐年分析提取每年的物候参数，通过式(3)～式(5)依次计算三个植被指数EVI、NDVI以及LSWI获得不同的物候参数估 计值；
Figure BDA0003526955040000081
Figure BDA0003526955040000082
Figure BDA0003526955040000083
式中，ρBlue、ρRed、ρNIR和ρswir1分别是蓝光，红光，近红外和短波红外1波段的地 表反射率值；
S32、基于三个植被指数，使用逻辑回归函数分别提取各自表达的年SOS，EOS和LOS，使用式(6)对每年卫星观测值导出的植被指数数据进行建模；
Figure BDA0003526955040000091
式中，t为观测时间，y(t)为植被指数值EVI、NDVI或者LSWI)a和b为拟合参 数，c+d为最大植被指数值，d为初始背景植被指数值；
通过式(7)计算拟合模型的曲率变化率；
Figure BDA0003526955040000092
式中，y(t)′和y(t)″分别是逻辑斯蒂克方程的一阶和二阶导数，ROC是拟合模型的曲率变化率；
当曲率的变化率首次达到极大值时对应的天数即为植被生长开始期SOS，当曲率的 变化率最后一次达到极小值时的天数为植被生长结束期EOS，取EOS与SOS的差值即 是植被生长季长度LOS。
S4、顾及物候参数和拟合MCCDC模型系数的森林生物量制图：
野外实地调查方法是估算森林AGB最准确的方法，但是此方法并不适合在大面积和长时间段内计算生物量。作为替代方案，Landsat观测影像，可用于估算大面积的森 林AGB。雷达数据，如PALSAR数据，可以很好地反映森林AGB信息。因此将Landsat 光学传感器与PALSAR结合起来被认为能够更好地提高AGB的评估精度。
本发明顾及物候参数和拟合MCCDC模型系数的森林生物量制图，如下步骤：
S41、对于森林一类调查数据，可以直接提取样地点对应的像元值；对于森林二类调查数据，需要统计polygon覆盖所有像元的平均值，进而作为该样地点的预测变量。 模型预测变量主要有3个不同的来源，一是合成的光学遥感特征值；二是具有明显季相 节律特征的模型参数和物候期；三是不受云雨干扰的雷达后向散射系数，如表2。
Figure BDA0003526955040000093
Figure BDA0003526955040000101
表2 AGB建模的预测变量
S42、AGB建模模型的无关变量和冗余变量会降低模型的预测能力。因此，先使用随机森林分析包的特征重要性(feature importance)命令，通过均方差百分比增量(PercentIncMSE)和节点纯度增量(IncNodePurity)两个测度对所有变量的重要性进 行排序，从而剔除那些不重要的变量。PercentIncMSE定义为当从决策树中排除给定变 量后准确性的降低的幅度。IncNodePurity度量了由于给定变量的分裂而导致的节点杂 质的减少。PercentIncMSE和IncNodePurity值越高，对应的这些变量就越重要。基于这 两个指标的排序，本实施例中最终选择了16个变量，如图2所示，分别为Red_3，NIR， NIR_1，NIR_3，SWIR1，SWIR1_3，SWIR2_5，PCA1，PCA3，HV，HHHV，SOS， EOS，LOS，DEM，Wetness。NIR_1在这里指代的是NIR波段涉及的MCCDC的第一 个系数，其它的以此类推。
S43、基于光学和雷达数据的预测模型建立，对于所有来源的野外实测生物量数据， 随机将这些数据按照8:2的比例分为训练数据集和验证数据集，使用筛选后的预测变量 与三年的实地调查数据一一匹配，采用RF、SVM和SGB三种机器学习算法，基于森 林树种(组)类型，依次建立AGB所有年份总预测模型。RF是一种集成的非参数技术， 主要用于数据的回归和分类。SVM可以处理有限数量的训练样本。SGB与随机森林相 似，从伪残差(即前一棵树的损失函数的梯度)顺序构建多个小型分类或回归树。
本实施例森林AGB制图精度评价：
本发明需要评估AGB拟合模型的性能，以确定该模型是否可以反映实际情况。考虑不同的变量组合，不同的建模方法以及森林树种分类对AGB模型的影响，将20％的 实地AGB样本(样本来自森林一类/二类调查数据)带入创建的模型以获得AGB预测 值，然后将这些预测值与实际观测值进行比较。选择R2和RMSE，平均绝对误差(Mean Absolute Error，MAE)和均方根误差百分比(RMSE％)来评价生物量模型的性能，即 式8，式9，式10和式11。通常，R2越大，RMSE，MAE和RMSE％越小，模型预测 能力越好。
为了进一步比较不同处理方式对AGB模型预测的影响，本实施例使用非参数Wilcoxon符号秩检验方法比较森林树种分类建模前后，不同变量组合和机器学习算法 下的AGB预测值的差异。在Wilcoxon符号秩检验中，它把观测值和零假设的中心位置 之差的绝对值的秩分别按照不同的符号相加作为其检验统计量。它适用于T检验中的成 对比较，但并不要求成对数据之差服从正态分布。
Figure BDA0003526955040000111
Figure BDA0003526955040000112
Figure BDA0003526955040000121
Figure BDA0003526955040000122
本实施例的森林AGB验证结果：
在确定了参与AGB建模的16个变量后，为了进一步表明MCCDC系数，物候参 数和机器学习算法对AGB建模的影响，本实施例还对AGB建模变量进行了不同的组 合进行精度验证，如表3。仅使用传统的光谱变量+雷达变量，RF算法的R2最高，为 0.52，SVM算法的R2值最低，只有0.34，SGB的R2值略小于RF算法，值是0.42。新 加入MCCDC系数和物候参数建模后，发现三类算法的验证精度均有所提高，其中RF 得到的R2从0.52增加到0.56，RMSE从18.62t/ha下降到14.05t/ha，MAE略有下降， RMSE％从52.49％下降到32.91％。SGB的R2从0.42增加到0.49，RMSE从18.28t/ha 下降至15.95t/ha，MAE从15.39t/ha下降到14.85t/ha，RMSE％从58.55％下降到43.12％。 SVM的R2基本保持不变，RMSE从21.24t/ha下降到18.79t/ha，MAE略有下降，RMSE％ 从66.60％下降到49.96％。
Figure BDA0003526955040000123
Figure BDA0003526955040000131
表3不同变量组合的模型验证精度
总体上，对于相同的输入变量，RF算法的预测能力更好，这一方面是因为重要性变量是基于RF算法的重要性特征命令得到的，因此输入变量与RF算法的相关性会高 于其它算法，另一方面SGB算法的能力低于RF，是因为SGB涉及的参数设置要多于 RF，模型的性能可能因这些参数选择不同而有差异。对于同一机器学习算法，新加入 的变量对三类算法精度都有所改进，因此加入的MCCDC参数和物候参数可以提高AGB 模型的精度。
在确定光学-物候-雷达变量后，本实施例进一步探究树种分类对AGB建模的影响。结果如表4，树种组分类后，所有模型的精度都有所提高。对于RF算法建模，核桃的 R2最高，为0.78，它的RMSE，MAE和RMSE％相对比较小。松类，栎类和其它森林 树种的R2分别是0.61，0.62和0.65。虽然其它森林树种的R2值要略高于松类和栎类， 但是它的RMSE，MAE和RMSE％都是最高的，分别是24.45t/ha，18.62t/ha和56.67％。 对于SVM模型，核桃的验证精度也是最高的，其次是松类。其它森林树种的R2最低， 只有0.41，RMSE是26.43t/ha，MAE是22.32t/ha，最值得注意的是它的RMSE％达到 了最大值68.39％。对于SGB模型，其它森林树种的R2最高，为0.66，但是它的RMSE， MAE和RMSE％也是最高的，分别是22.43t/ha，16.63t/ha和59.97％。
Figure BDA0003526955040000132
Figure BDA0003526955040000141
表4森林树种组分类下基于光学-物候-雷达模型的AGB的验证精度
总体而言，这三类模型均在核桃的AGB验证模型当中精度最高，对于剩余森林树种组AGB展示了不同的预测能力。RF和SGB分类算法的R2均是在松类中值最低，但 是RMSE和RMSE％的值在其它森林树种中是最高的。SVM在其它森林树种的验证精 度最低。
图3展示了树种分类建模前后所有验证样本点实测AGB与预测AGB的散点图。 对于RF模型，在不分森林树种建模时，R2是0.56，RMSE是14.05t/ha。分树种建模后， R2增加至0.67，RMSE是12.13t/ha。对于SVM模型，不分树种建模的验证精度R2只 有0.35，RMSE是18.79t/ha。分树种建模后，SVM的建模验证精度R2是0.49，RMSE 下降至16.73t/ha。对于SGB建模而言，不分树种建模，它的AGB验证精度R2为0.49，RMSE是15.95t/ha，分树种建模后，R2是0.67，RMSE是14.52t/ha。从散点图可以看 出在树种分类前后，都存在一些样本点的低值高估和高值低估现象，尤其是在AGB<15 t/ha和AGB>60t/ha。但是进行树种组分类后，这种误差现象有所减缓。分类建模后， 三个模型的y＝x这条线与趋势线越来越接近。总体而言，所有机器学习算法在分树种建 模后，验证精度都有所提高，这些证明森林树种分类的必要性。
本实施例的森林树种组分布时空动态：
本实施例绘制2000-2019年森林树种组空间分布图，如图4。从图上可以看出松类主要分布在西边山体的右侧，其它山体的边缘以及右下角。栎类主要分布在西边山体的 左侧，右上角以及研究区的中上部。核桃主要分布在研究区的最左侧边缘且连成一片。 在研究区的右上角，栎类树种分布之间也夹杂着核桃，其它森林树种零星地分布在森林 边缘和内部。
本实施例的森林AGB时空动态：
采用分树种建模的光谱变量+MCCDC系数+物候参数+雷达的RF算法，本实施例 绘制了从2000-2019年的AGB空间分布图，如图5。从AGB时空分布图看，整个研究 区的AGB在空间和时间上均存在一定的差异。首先对于研究区的左侧偏下区域，它的 AGB小于30t/ha，这一部分主要是核桃分布区。对于研究区的东南角，AGB普遍也低 于30t/ha，结合图4，这一部分主要是松类。研究区靠左侧的中部，东北角，中部靠上 区域的AGB比较高，大于55t/ha，这一部分主要是栎类。从2000年到2019年大部分 区域呈现AGB小幅度增加，其中2000年到2007年之间AGB的颜色变化并不明显， 到了2010年AGB开始出现颜色替换区域。首先对于研究区的西北角区域，AGB主要 从30-55t/ha增加到大于55t/ha。研究区中部靠上的区域出现绿色被黄色替代的现象。 东南角随时间变化并不明显，AGB依旧主要集中在小于30t/ha。
如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为 对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。

Claims (3)

1.一种顾及谐波模型系数和物候参数的森林生物量估算方法，其特征在于，包括如下步骤：
S1、影像获取及预处理；获取待分类区域一段时间内多个传感器的陆地卫星影像，并对获得的影像进行预处理；同时，采用最小二乘法拟合斜率和截距的方法，进行逐波段的光谱归一化，最小化不同传感器间的光谱响应差异；
S2、建立基于多源遥感数据的MCCDC模型，具体步骤包括：
S21、建立SR时间序列模型，定义一个逐像元反射率时间序列模型，如式(1)所示，以表达地表反射率变化的特征；
Figure FDA0003526955030000011
式(1)中包含一个傅立叶谐波模型和一个长期的趋势部分，x为儒略日期；T是周期天数，该值为365.25；
Figure FDA0003526955030000012
为第i波段在儒略日期x的地表反射率预测值；该时间序列模型中有8个系数，系数a0,i为第i波段清晰观测的总体均值，系数a1,i、a2,i、b1,i、b2,i、c1,i、c2,i用于模拟由物候引起的地表反射率的年内变化，它的最小正周期为一年，系数d1,i体现由气候变化和土壤退化等因素引起的地表反射率渐变的趋势；
S22、识别时间序列模型突变点，基于步骤S21建立的时间序列模型，看是否满足式(2)；
Figure FDA0003526955030000013
式(2)中，k为波段数量；n指需要判断的突变点后连续观测值的个数；ρ(i，x)为第i波段在儒略日期x地表反射率的实际观测值；RMSEi是指Landsat第i波段的模型预测值与观测值的均方根误差；
当满足式(2)时，证明区域内发生了土地覆盖变化，对变化的像元按照时间顺序排序，首先利用同一位置的前24个清晰像元，基于式(1)，建立初始化的时间序列模型，然后判断每个波段接下来的4个像元的实际观测值与预测值的差异是否都超出模型均方根误差的3倍，如果没有超出范围，则将这些观测值加入到之前参与建模的观测值中去，重新拟合参数，再次判断后续的观测值，如果连续4个观测值都超出了范围，则定义为突变点，即土地覆盖变化发生的时间；
S23、合成每日遥感图像，利用步骤S21～S22的方法生成每个像元不同波段的时间序列模型，对于每个像元的单一波段，通过给定变量儒略日期的值，依据时间序列模型，生成该像元给定日期的SR预测值，将该方法应用到同一时间内所有像素点，从而生成空间分辨率为30m的清晰合成图像；
S3、基于逻辑回归方程的森林年物候提取，具体方法为：
S31、基于获得的全部每天合成图像，进行逐年分析提取每年的物候参数，通过式(3)～式(5)依次计算三个植被指数EVI、NDVI以及LSWI获得不同的物候参数估计值；
Figure FDA0003526955030000021
Figure FDA0003526955030000022
Figure FDA0003526955030000023
式中，ρBlue、ρRed、ρNIR和ρswir1分别是蓝光，红光，近红外和短波红外1波段的地表反射率值；
S32、基于三个植被指数，使用逻辑回归函数分别提取各自表达的年SOS，EOS和LOS，使用式(6)对每年卫星观测值导出的植被指数数据进行建模；
Figure FDA0003526955030000024
式中，t为观测时间，y(t)为植被指数值EVI、NDVI或者LSWI，a和b为拟合参数，c+d为最大植被指数值，d为初始背景植被指数值；
通过式(7)计算拟合模型的曲率变化率；
Figure FDA0003526955030000025
式中，y(t)'和y(t)″分别是逻辑斯蒂克方程的一阶和二阶导数，ROC是拟合模型的曲率变化率；
当曲率的变化率首次达到极大值时对应的天数即为植被生长开始期SOS，当曲率的变化率最后一次达到极小值时的天数为植被生长结束期EOS，取EOS与SOS的差值即是植被生长季长度LOS；
S4、顾及物候参数和拟合MCCDC模型系数的森林生物量制图，具体包括如下步骤：
S41、选择AGB预测模型变量来源，包括合成的光学遥感特征值、具有明显季相节律特征的模型参数和物候期以及不受云雨干扰的雷达后向散射系数；
S42、在步骤S41的预测模型变量来源中筛选出重要变量，先使用随机森林分析包的特征重要性命令，通过均方差百分比增量和节点纯度增量两个测度对所有变量的重要性进行排序，剔除不重要的变量，均方差百分比增量和节点纯度增量值越高，对应的变量越重要；
S43、基于光学和雷达数据的预测模型建立，对于所有来源的野外实测生物量数据，随机将这些数据按照8:2的比例分为训练数据集和验证数据集，使用筛选后的预测变量与三年的实地调查数据一一匹配，采用RF、SVM和SGB三种机器学习算法，基于森林树种或树种组类型，依次建立AGB所有年份总预测模型。
2.根据权利要求1所述的顾及谐波模型系数和物候参数的森林生物量估算方法，其特征在于，步骤S1中对获得的陆地卫星影像进行大气校正获得SR图像。
3.根据权利要求1所述的顾及谐波模型系数和物候参数的森林生物量估算方法，其特征在于，步骤S4中通过均方差百分比增量和节点纯度增量两个测度对所有变量的重要性进行排序，剔除不重要的变量后，最终选择16个变量，分别为Red_3，NIR，NIR_1，NIR_3，SWIR1，SWIR1_3，SWIR2_5，PCA1，PCA3，HV，HHHV，SOS，EOS，LOS，DEM和Wetness。
CN202210195246.3A 2022-03-01 2022-03-01 顾及谐波模型系数和物候参数的森林生物量估算方法 Active CN114821349B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210195246.3A CN114821349B (zh) 2022-03-01 2022-03-01 顾及谐波模型系数和物候参数的森林生物量估算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210195246.3A CN114821349B (zh) 2022-03-01 2022-03-01 顾及谐波模型系数和物候参数的森林生物量估算方法

Publications (2)

Publication Number Publication Date
CN114821349A true CN114821349A (zh) 2022-07-29
CN114821349B CN114821349B (zh) 2023-07-07

Family

ID=82529488

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210195246.3A Active CN114821349B (zh) 2022-03-01 2022-03-01 顾及谐波模型系数和物候参数的森林生物量估算方法

Country Status (1)

Country Link
CN (1) CN114821349B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115620151A (zh) * 2022-12-16 2023-01-17 中化现代农业有限公司 物候期识别方法、装置、电子设备和存储介质
CN116563716A (zh) * 2023-07-07 2023-08-08 吉林省林业科学研究院(吉林省林业生物防治中心站) 用于森林碳汇数据采集的gis数据处理系统
CN116579521A (zh) * 2023-05-12 2023-08-11 中山大学 产量预测时间窗口确定方法、装置、设备及可读存储介质
CN117292267A (zh) * 2023-11-27 2023-12-26 武汉大学 一种基于物候信息的水稻地上生物量分段估算方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200225075A1 (en) * 2019-01-14 2020-07-16 Wuhan University Method and system for optical and microwave synergistic retrieval of aboveground biomass
CN112395808A (zh) * 2020-11-18 2021-02-23 南京林业大学 一种结合随机森林和协同克里金的生物量遥感制图方法
CN113850139A (zh) * 2021-08-25 2021-12-28 南京林业大学 一种基于多源遥感的森林年际物候监测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200225075A1 (en) * 2019-01-14 2020-07-16 Wuhan University Method and system for optical and microwave synergistic retrieval of aboveground biomass
CN112395808A (zh) * 2020-11-18 2021-02-23 南京林业大学 一种结合随机森林和协同克里金的生物量遥感制图方法
CN113850139A (zh) * 2021-08-25 2021-12-28 南京林业大学 一种基于多源遥感的森林年际物候监测方法

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115620151A (zh) * 2022-12-16 2023-01-17 中化现代农业有限公司 物候期识别方法、装置、电子设备和存储介质
CN116579521A (zh) * 2023-05-12 2023-08-11 中山大学 产量预测时间窗口确定方法、装置、设备及可读存储介质
CN116579521B (zh) * 2023-05-12 2024-01-19 中山大学 产量预测时间窗口确定方法、装置、设备及可读存储介质
CN116563716A (zh) * 2023-07-07 2023-08-08 吉林省林业科学研究院(吉林省林业生物防治中心站) 用于森林碳汇数据采集的gis数据处理系统
CN116563716B (zh) * 2023-07-07 2023-09-08 吉林省林业科学研究院(吉林省林业生物防治中心站) 用于森林碳汇数据采集的gis数据处理系统
CN117292267A (zh) * 2023-11-27 2023-12-26 武汉大学 一种基于物候信息的水稻地上生物量分段估算方法及系统
CN117292267B (zh) * 2023-11-27 2024-02-02 武汉大学 一种基于物候信息的水稻地上生物量分段估算方法及系统

Also Published As

Publication number Publication date
CN114821349B (zh) 2023-07-07

Similar Documents

Publication Publication Date Title
Tian et al. Comparison of UAV and WorldView-2 imagery for mapping leaf area index of mangrove forest
Meng et al. Mapping canopy defoliation by herbivorous insects at the individual tree level using bi-temporal airborne imaging spectroscopy and LiDAR measurements
Yang et al. Daily Landsat-scale evapotranspiration estimation over a forested landscape in North Carolina, USA, using multi-satellite data fusion
Solano et al. A methodology based on GEOBIA and WorldView-3 imagery to derive vegetation indices at tree crown detail in olive orchards
Peña-Barragán et al. Object-based crop identification using multiple vegetation indices, textural features and crop phenology
Steininger Tropical secondary forest regrowth in the Amazon: age, area and change estimation with Thematic Mapper data
CN114821349B (zh) 顾及谐波模型系数和物候参数的森林生物量估算方法
Kuusinen et al. Effects of forest age on albedo in boreal forests estimated from MODIS and Landsat albedo retrievals
Olthof et al. Treeline vegetation composition and change in Canada's western Subarctic from AVHRR and canopy reflectance modeling
Yaney-Keller et al. Using Unmanned Aerial Systems (UAS) to assay mangrove estuaries on the Pacific coast of Costa Rica
CN112348812B (zh) 林分年龄信息测量方法及装置
CN113850139B (zh) 一种基于多源遥感的森林年际物候监测方法
US20220392215A1 (en) System and Method for Mapping Land Cover Types with Landsat, Sentinel-1, and Sentinel-2 Images
da Silveira et al. Use of MSI/Sentinel-2 and airborne LiDAR data for mapping vegetation and studying the relationships with soil attributes in the Brazilian semi-arid region
Launeau et al. Airborne hyperspectral mapping of trees in an urban area
Van Beek et al. Reducing background effects in orchards through spectral vegetation index correction
Geng et al. Vegetation coverage of desert ecosystems in the Qinghai-Tibet Plateau is underestimated
Badr et al. Estimating growing season length using vegetation indices based on remote sensing: A case study for vineyards in Washington State
López-García et al. Yield estimations in a vineyard based on high-resolution spatial imagery acquired by a UAV
Bindhu et al. Development of a spatio-temporal disaggregation method (DisNDVI) for generating a time series of fine resolution NDVI images
Pan et al. Using QuickBird imagery and a production efficiency model to improve crop yield estimation in the semi-arid hilly Loess Plateau, China
Fadaei Advanced land observing satellite data to identify ground vegetation in a juniper forest, northeast Iran
Danoedoro et al. Combining pan-sharpening and forest cover density transformation methods for vegetation mapping using Landsat-8 Satellite Imagery
Shao et al. Review of Selected Moderate-Resolution Imaging Spectroradiometer Algorithms, Data Products, and Applications
Mngadi et al. Quantifying carbon stock variability of species within a reforested urban landscape using texture measures derived from remotely sensed imagery

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
CB03 Change of inventor or designer information

Inventor after: Zhang Yali

Inventor after: Zhang Yin

Inventor after: Li Mingshi

Inventor after: Jiang Lufan

Inventor after: Liu Qinqin

Inventor after: Wang Nan

Inventor after: Yang Boxiang

Inventor after: Ye Junhong

Inventor after: Peng Yuwen

Inventor after: Liu Jiawei

Inventor before: Li Mingshi

Inventor before: Zhang Yin

Inventor before: Zhang Yali

Inventor before: Jiang Lufan

Inventor before: Liu Qinqin

Inventor before: Wang Nan

Inventor before: Yang Boxiang

Inventor before: Ye Junhong

Inventor before: Peng Yuwen

Inventor before: Liu Jiawei

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant