本发明属于遥感反演技术领域，特别涉及森林地上生物量估算方法。
背景技术：
森林地上生物量是表征森林碳储量的重要参数，准确估算森林地上生物量对深入研究森林演替、人类活动、自然干扰和气候变化等方面具有重要意义。传统地上生物量的监测手段主要通过实地测量实现，虽然该方法精度高，但耗时耗力，成本高，采样通常离散、稀疏、难以大范围覆盖。遥感技术作为一种新型获取数据的手段，因其覆盖范围广、获取数据周期短、成本低、不受时空限制等优点使得大尺度变化监测变成可能。
目前发展了诸多森林生物量遥感估算法，如lidar、sar、光学遥感等。其中，星载lidar目前难以实现大范围时空连续覆盖，星载sar虽具有较高时间和空间分辨率，但目前sar卫星寿命较短，很少有同一系列卫星具有超过10年的观测周期，且部分sar卫星并未免费公开，其历史存档数据很少，极大地限制了大范围森林地上生物量(agb)的连续监测。相比之下，在开展大范围森林agb长时间序列监测方面，光学数据landsat因其高时空分辨率以及免费公开的长时间历史存档数据而具有得天独厚的优势。如何基于单一光学数据开展森林agb准确估算，是实现近几十年以来大范围森林生物量及相应碳收支深入研究的关键。
然而，现有生物量光学估算方法多基于单时相、或者部分多时相数据开展，这导致森林生物量估算还存在两大问题。一是数据缺乏时空连续性问题，光学数据自身容易受云雾干扰，导致光学数据时空连续性受到严重限制。二是森林生物量估算精度低的问题，单时相光学数据忽略了森林季节变化特征，单一基于地表瞬时观测光谱反射率开展agb估算，光谱反射率主要敏感于森林冠层的叶片特征，然而森林叶片只占据森林agb的2～11％，加之单时相数据中不可避免的噪声问题，导致森林agb估算精度始终受到限制。为此，本发明拟利用所有可获取的landsat数据开展时序建模，提取森林时序变化特征，过滤信号中的噪声影响，并创新性地将森林时序特征参数引入到传统森林地上生物量估算中，达到提高森林agb估算精度的目的，为大范围agb制图提供技术支撑。
技术实现要素：
针对现有基于单时相光学数据的agb估算模型存在数据时空连续性差、易饱和、易受噪声干扰、估算精度差等问题，本发明提供了一种基于光学数据时间序列建模的森林agb估算方法，并以landsat数据为例展开方法说明。其中创新性地将森林时序变化参数引入到agb估算中，并构建一套完整的agb估算指标体系，借此有效降低噪声影响，缓解饱和问题，提高基于光学数据的森林agb估算精度与普适性，为大范围森林agb估算提供技术支持。
本发明提供的技术方案为：一种基于landsat时间序列建模的森林地上生物量估算方法，该方法包括：
步骤1：收集不同地块的经度、纬度、高程、东西朝向坡度、南北朝向坡度，这些地块不同时间点的生物量和对应时间点如下波段的反射率：蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2)，热红外波段(tir)和质量评估波段(qa)；
步骤2：建立时间序列模型；
步骤2.1：根据qa波段筛除步骤1数据中的云、云阴影、雪的观测，并基于清晰观测值和如下公式分别建立蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2)和热红外波段(tir)的时间序列模型；
其中：
步骤2.2：变化监测；
对比绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2)的时间序列模型预测的反射率与真实反射率，如果两者差值连续在6个观测时相上都大于2倍均方根误差，则视其发生突变，并完成突变前的时间序列建模。基于突变时相后的所有清晰像元，重复步骤2.1和上述步骤，直至最后一个时相处的清晰观测完成迭代。至此，完成该像元整个时间序列上的变化监测，获取所有时间序列建模。
步骤3：构建森林生物量估算模型；
步骤3.1：计算七个波段的总体反射率
采用步骤2得到的时间序列模型，分别计算蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2)和热红外波段(tir)由时序模型衍生的总体反射率(overallreflectance)。计算公式如下。
tyear是朱利安日期，设定为所选取年份的7月2日(年内中间天数)。
步骤3.2：基于机器学习方法训练模型
该模型的输入包括三类指标地理指标，经度、维度、数字地表模型(dsm)、东西朝向坡度、南北朝向坡度；时序变化参数指标，时序变化参数中的a1，b1，c1，rmse四种时序变化参数；七个波段的总体反射率，蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2)和热红外波段(tir)。输出为森林生物量。本发明此处选取、但不限于随机森林(rf)这一机器学习方法开展agb估算模型构建。
步骤3.3：基于步骤3.2训练的模型，应用到目标日期目标地块生物量估算。
本发明在随机森林模型训练时，创新性地将landsat时序建模特征引入到森林agb估算中，构建一套完整的森林agb估算指标体系，可从三大方面提高现有传统单时相森林agb估算方法的精度。一是landsat时序建模可有效去除单时相中噪声的影响，以此更加稳定且准确地表征森林光谱反射率，去除外界干扰对agb估算精度的影响。第二，将时序特征引入到agb估算中，弥补现有传统agb估算方法忽略森林时间维变化特征的不足，提高agb估算精度；其中特别忽略a2，b2，a3，b3年间多周期变化系数以达到提高本方法计算精度与效率的目的。第三，构建完整的agb估算指标体系。引入了地理位置、高程、坡度坡向等地理信息，间接引入了气候信息，提高了生物量在不同地理位置的计算精度。综合空间维光谱反射率指标、时间维的时序变化参数指标、以及地理指标三大指标类型进一步提高森林agb估算精度。总体而言，本发明可有效去除光学数据缺乏时空连续性问题，缓解森林agb饱和问题，有效提高agb估算精度与普适性，适用于大范围森林agb长时序制图。
附图说明
图1为代表样点的时间序列观测数据及时序拟合曲线，右下角展示了该像元在不同年份的googleearth图像。
图2为针对所有样点在其采样年份四个季节的模拟反射率与观测反射率散点图对比。
图3为基于rf模型不同季节下综合单时相所有波段反射率的agb估算精度；(a)观测反射率估算结果，(b)模拟反射率估算结果，(c)将(b)中春季估算模型运用于其他季节的估算结果。
图4为基于rf模型不同指标体系下agb估算精度；(a)传统基于六个波段(蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2))的agb估算结果，(b)基于时序模型衍生的七个波段overallreflectance的agb估算结果，(c)引入时序变化特征的agb估算结果，(d)引入地理指标后的agb估算结果，即本发明所提出的最终agb估算模型。柱状图展示了不同输入变量的重要性排序，此处只列出了排名前15的参数。
具体实施方式
下面结合具体实施例和说明书附图对本发明提供的基于时间序列landsat估算森林地上生物量的方法作进一步说明：
(1)数据集介绍
实测生物量数据集为thebiomassplot142library数据，由tern2汇编和出版的现场生物量清单数据库。该数据库主要分布于澳大利亚全国范围，收集了联邦、州和地方政府部门、大学、私营公司和其他机构的stem库存数据。主要的空间分辨率0.05-1ha，空间覆盖为110.00to155.001329e,-10.00to-45.000512n，空间参考wgs84。此次共收集采样时间在2000年以后的数据，一共7771个点，此外额外加入了5759个扩充点(其agb为0tons/ha)。长时间序列landsat数据集通过gee平台获取，输入采样点中心坐标提取采样点从2000年至2019年的usgs官方发布的地表反射率二级产品，包括蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)、短波红外2(swir2)，热红外波段(tir)和质量评估波段(qa)。
(2)反射率时序模拟效果展示
图1展示了某一代表样点的反射率时序拟合结果。该图展示了时序模型可准确捕捉森林反射率的时序变化状况。googleearth高分辨图像展示了该像元中的森林于2013年发生了扰动，ccdc模型有效监测到了该扰动。
图2展示了采样点所在年份四个季节下(夏季：1月15，秋季：4月15，冬季：7月15，春季：10月15日)观测反射率与模拟反射率之间的精度对比。结果显示ccdc模拟反射率与观测反射率高度相关，且能有效去除观测反射率中的部分异常值(未去除的云雾干扰、噪声等影响)。
(3)单时相森林agb估算结果
表1展示了不同季节下，单时相模拟反射率与观测反射率与agb之间的相关性。可看出不同季节之间的相关度差异极大，其中秋季效果最差，其次冬季，其次夏季，春季相关度最高。其中，所有波段的模拟反射率与agb相关度都高于观测反射率，这表明时间序列建模可去除噪声的影响，提高光学数据的agb表征能力
图3展示了基于rf模型不同季节下，综合所有波段反射率的agb估算精度，其中rf基于matlab中的treebagger函数实现，并采用4-folder进行交叉验证。可看出基于模拟反射率的agb精度都高于观测反射率的估算精度。图3c展示了将春季构建的模型运用到夏、秋、冬季中，结果显示基于某一季节或时相构建的模型难以用于其他时相，否则agb估算精度会急剧下降。而在实际应用中几乎不可能在某一固定时相获取大范围空间连续的观测数据，换言之，难以基于观测反射率构建出一个普适通用的agb估算模型，这极大限制了agb的大范围制图。而时序模型可模拟任一时相的反射率，可有效解决该问题。
(4)基于时间序列模型的森林agb估算结果
图4展示了不同agb估算模型的估算精度。图4a为目前主流agb估算法所采用的波段(蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(nir)、短波红外1(swir1)和短波红外2(swir2)，共六个波段)，其估算精度为r2＝0.43，rmse＝87.4tons/ha。图4b为本发明引入热红外(tir)波段后的结果，即采用本发明所构建指标体系中的第一类—光谱反射率指标，可看出tir波段可有效提高agb估算精度，r2＝0.49，rmse＝82.7tons/ha。图4a和b都为基于单一时相的估算精度，尽管已利用时序模型去除噪声等影响，但其估算精度已达到上限。图4c为引入时序变化参数后的估算结果，即综合了本发明所构建指标体系中的第一类和第二类指标。可看出时序变化参数的引入显著提高了agb估算精度，r2＝0.56，rmse＝76.5tons/ha，其中，明显降低了森林agb估算饱和问题这一重要难题。指标重要性结果显示，b1、a1、rmse等时序参数在agb模型中占据较高重要性，表明时序参数具有单时相数据没有的丰富信息，能弥补传统单时相agb估算的不足。图4d为在图4c基础上引入地理指标后的结果，即综合了本发明所构建指标体系中所有的第一、二和三类指标。agb估算精度得到了进一步提高，r2＝0.59，rmse＝73.4tons/ha。且lat、lon和dsm在参数重要性排序靠前，这表明地理参数能提供独立于光学数据的特征，能进一步提高agb光学估算法的精度。因此，相比于传统单时相观测数据的agb估算法(图3a，最优精度为春季，r2＝0.41，rmse＝88.8tons/ha)，本发明提出的基于时序建模的agb估算方法(图4d)将r2提高了0.18，rmse降低了15.4tons/ha。
表1针对所有样点在其采样年份四个季节的模拟反射率与观测反射率散点图对比。
1.一种基于landsat时间序列建模的森林地上生物量估算方法，该方法包括：
步骤1：收集不同地块的经度、纬度、高程、东西朝向坡度、南北朝向坡度，这些地块不同时间点的生物量和对应时间点如下波段的反射率：蓝波段、绿波段、红波段、近红外波段、短波红外1、短波红外2、热红外波段和质量评估波段；
步骤2：建立时间序列模型；
步骤2.1：根据qa波段筛除步骤1数据中的云、云阴影、雪的观测，并基于清晰观测值和如下公式分别建立蓝波段、绿波段、红波段、近红外波段、短波红外1、短波红外2、热红外波段的时间序列模型；
其中：
步骤2.2：变化监测；
对比绿波段、红波段、近红外波段、短波红外1、短波红外2的时间序列模型预测的反射率与真实反射率，如果两者差值连续在6个观测时相上都大于2倍均方根误差，则视其发生突变，并完成突变前的时间序列建模。基于突变时相后的所有清晰像元，重复步骤2.1和上述步骤，直至最后一个时相处的清晰观测完成迭代；至此，完成该像元整个时间序列上的变化监测，获取所有时间序列建模；
步骤3：构建森林生物量估算模型；
步骤3.1：计算七个波段的总体反射率
采用步骤2得到的时间序列模型，分别计算蓝波段、绿波段、红波段、近红外波段、短波红外1、短波红外2和热红外波段由时序模型衍生的总体反射率；计算公式如下；
tyear是朱利安日期，设定为年内的一天；
步骤3.2：基于机器学习方法训练模型
该模型的输入包括三类指标地理指标，经度、维度、数字地表模型、东西朝向坡度、南北朝向坡度；时序变化参数指标，时序变化参数中的a1，b1，c1，rmse四种时序变化参数；七个波段的总体反射率，蓝波段、绿波段、红波段、近红外波段、短波红外1、短波红外2和热红外波段；输出为森林生物量。本发明此处选取、但不限于随机森林这一机器学习方法开展agb估算模型构建；
步骤3.3：基于步骤3.2训练的模型，应用到目标日期目标地块生物量估算。
技术总结