CN113221765B - 一种基于数字相机影像有效像元的植被物候期提取方法 - Google Patents

Abstract

本发明公开了一种基于数字相机影像有效像元的植被物候期提取方法，属于植被参数遥感反演方法技术领域，本发明针对物候相机或普通数字相机，利用逐日多时相观测影像计算像元尺度的绿度植被指数，基于绿度植被指数在年内变化的振幅大小，筛选振幅阈值，对影像感兴趣区进行地物识别，快速地标识植被叶片像元和土壤、枝干等非叶片像元，在求取逐日影像感兴趣区植被绿度指数均值时，剔除非叶片像元，提高植被绿度指数均值的信噪比，并利用降噪后的植被指数均值时序数据提取植被生长周期的关键物候期，属于植被定量遥感方法的研究领域。本发明可有效降低地表景观变化对感兴趣区植被物候提取的影响，对中低覆盖度落叶植被具有更好的适用性。

Description

一种基于数字相机影像有效像元的植被物候期提取方法

技术领域

本发明属于植被参数遥感反演方法技术领域，具体涉及一种基于数字相机影像有效像元的植被物候期提取方法。

背景技术

植物物候是指植物受气候等环境因素影响而出现的以年为周期的自然现象，是植物长期适应季节性变化的环境而形成的生长发育节律。通常而言，物候是指植物生长周期中某一事件发生的时间节点或不同状态之间的转折点，如落叶植物的发芽、常绿植物叶片从休眠状态进入光合作用状态的转折点等。Mark D. Schwartz等人2013年在《Phenology:An Integrative Environmental Science》一书中综述了物候研究进展，指出从植物个体到生态系统各个尺度，众多过程都由物候直接或间接调控，特别是与碳-氮-水耦合循环相关的生理生态过程，是地球科学研究的核心参数之一。物候对气候变化非常敏感，其变化可以反映出生物圈对气候、水文等环境因素的季节与年际响应与反馈，在全球变化研究中受到广泛关注。

早期人们通过野外目视观测来监测植物物候，主要针对某一个或多个物种进行的单株定点观测，来记录植被个体或群落的发芽、展叶、枯黄和落叶等物候日期。这种传统监测方法易于操作且时间精度高，但存在着观测范围有限、观测标准不一致、观测对象个体差异大等局限，无法实现大尺度的植物物候监测。随着遥感技术的发展，遥感数据能快速准确地反映地表植被信息，定量反映区域至全球尺度的植被生态特征，因此迅速应用于物候监测。与传统物候不同，遥感方法用于描述整个地表景观的物候变化，侧重表征植被结构或功能状态的转换，因此又称为遥感地表物候。物候遥感监测主要是利用植被指数(VegetationIndex)指示植被冠层绿度或覆盖度的周期性变化。研究者根据需求，在植被时序曲线上设定阈值或选取具有物理意义的节点，指示植被的关键物候事件或划分状态阶段。夏传福等人2013年在《遥感学报》发文，对近年来植物物候的遥感监测方法进行了综述与评价，并对遥感地表物候的验证方法、误差来源等进行了讨论，指出卫星遥感主要针对较大范围的植被覆盖区，但由于遥感器时空分辨率的限制，遥感地表物候监测在中小尺度的植被物候提取理论与方法还存在不足。因此，如何将遥感地表物候与真实植被物候进行关联，是物候遥感监测研究中的热点与难点。

近年来，遥感监测平台不断更新，以多光谱物候相机为代表的地基遥感平台迅速发展，在物候遥感监测领域日益受到关注。物候相机本质为数字相机，通常为红、绿、蓝(RGB)三波段，部分机型配有红外波段以及高光谱分辨率。数字相机影像分辨率高，可实现高频率连续自动观测，且布设灵活、成本较低，成为监测植被物候的一种新手段，已广泛应用于站点尺度植被物候动态的监测与研究。其中，Andrew D.Richardson等人2018年在《Scientific Data》上发表的《Tracking vegetation phenology across diverse NorthAmerican biomes using PhenoCam imagery》为最具代表性的工作之一。此外，物候相机的观测范围为单木到景观尺度，既可以和野外观测物候相对应，也可与遥感影像象元相匹配，因此被认为是地面观测与卫星遥感进行尺度转换的重要途径。因此，利用多光谱物候相机数据提取植被物候，是植被物候研究的重要发展趋势，将为区域尺度卫星遥感应用和生态系统过程模型改进提供理论基础和技术支持。

利用数字相机影像提取地表物候的技术流程主要包括区域选取与像素值计算、植被指数计算、时序数据过滤、数据拟合、物候参数提取方法、不确定性评估与对比研究等步骤。Michael Toomey等人在2015年《Ecological Applications》发表题为《Greennessindices from digital cameras predict the timing and seasonal dynamics ofcanopy-scale photosynthesis》的论文，对主流的数字相机数据及预处理、分析方法等进行了详尽说明，为物候监测应用提供有益参考。当前研究中，通常采用自定义感兴趣区(Area of Interest，AOI)作为植被物候提取的研究区域：首先选取一张质量较好的照片作为该站点的参考，加载参考影像后选取多个点围成封闭的AOI，然后将选取的AOI存储为多边形，继而用该多边形自动裁剪时间序列照片获取相应的AOI，对整个AOI中的像素值进行平均，提取AOI的数值组成连续的时间序列，并应用于后续物候参数提取过程。上述方法效率高、易于操作，且精度较好。

然而，对于非均匀下垫面，特别是植被稀疏区域，AOI中可能会存在大量的土壤、植被枝干等无效像元。如将这些无效像元与植被冠层(叶片)有效像元共同参与AOI像素值计算，必然会降低影像信噪比，在植被指数时序数据中引入噪声，进而影响地表物候提取精度。如引入人工目视解译或影像分类技术识别 AOI有效像元，则可能会增加人工或运算成本；同时，易受到待解译图像拍摄时间的影响，因为不同季节、时相的相机影像，可能得到的有效像元解译结果不一致，导致新的不确定性。

发明内容

发明目的：本发明的目的在于提供一种基于数字相机影像有效像元的植被物候期提取方法，用于数字相机影像中植被叶片有效像元识别，剔除影像感兴趣区域土壤、枝干等非叶片像元，进而降低目标植被绿度指数均值噪声，提取植被关键物候期，作为生态系统物候过程模拟与预测的基础。

技术方案：为了实现上述目的，本发明提供一种基于数字相机影像有效像元的植被物候期提取方法，包括如下步骤：

1)选定影像感兴趣区；

根据数字相机观测影像，选取感兴趣区AOI；

2)影像数据预处理；

设相机拍摄频率为n景/小时，每天拍摄时长为m小时，则针对每日多时相影像n×m景，利用红、绿、蓝三波段数据，逐景逐像元计算绿度指数Gcc；

3)植被指数振幅计算；

针对AOI的Gcc三维矩阵，进行逐像元处理；

4)有效像元识别；

针对AOI的AMP二维矩阵，将所有像元的AMP进行排序，得到AMP升序曲线，其横坐标为像元的排列序号，纵坐标为对应像元的AMP；基于分段回归模型，采用逐点测试的方式，识别AMP升序曲线的突变点；

5)植被指数均值时序曲线计算；

为降低误分风险，对有效像元AMP再次进行排序，选取AMP大于AMP序列下四分位数的有效像元，为最终有效像元；对Gcc三维矩阵中每一天的所有最终有效像元的Gcc求平均，得到AOI有效像元逐日Gcc均值的时序曲线；

6)遥感地表物候提取；

利用五参数Logistic方程，对筛选后的逐日Gcc均值时序曲线进行拟合，并根据拟合方程计算曲率变化率RCC，对AOI有效像元逐日Gcc均值时序数据提取开始展叶期与成熟叶片期的物候日期，实现关键物候期识别。

进一步地，所述的步骤1)具体为：AOI为矩形或闭合多边形边框，须完整覆盖目标地物，并剔除非目标地物；记录AOI的影像坐标，用于时序影像提取；观测期间数字相机影像视场发生变化，须先对影像进行配准，以保证时序影像 AOI空间范围一致。

进一步地，所述的步骤2)绿度指数Gcc，如公式(I)所示：

其中，i为观测日期，取值为1～365；j为每天中的第j景影像，取值为1～n×m； p和q分别为目标像元在影像中的行号和列号，用于确定像元位置；R、G、B分别为红、绿、蓝波段的影像像素值；之后，再对每天中每一个像元的n×m个Gcc 求取中值，作为当日该像元的Gcc值，进而得到一年中AOI逐日Gcc三维矩阵。

进一步地，所述的步骤2)中，考虑到太阳角度对植被冠层光谱反射的影响，选取每日上午9:00以后至下午3:00之前作为拍摄时间，且以正午12:00为中心时间，前后的拍摄时长相等，以体现每日Gcc中值的代表性。

进一步地，所述的步骤3)中，所述的进行逐像元处理具体为：提取每一像元Gcc时序数据，即该像元按时间排列的逐日Gcc；利用奇异谱分析SSA方法对Gcc时序数据进行噪声去除、缺值插补；计算时序曲线生长季内的Gcc振幅 AMP，采用Gcc夏季高值与春季低值之差表征振幅；经AOI逐像元计算，得到 AOI中每一个像元的AMP，即AMP二维矩阵。

进一步地，所述的步骤4)中，具体方法为：选取AMP升序曲线上的任一点为测试点bp，即假设的突变点，以该测试点构建分段回归模型，如公式(II) 所示：

其中，x为AMP序列中像元的排列序号(即横坐标值)，x_bp为待测试点的排列序号(即横坐标值)，AMP_x为x所对应的AMP，k₁、b₁、k₂、b₂分别为待测试点前后两段线性回归方程的参数；根据序号x与对应的AMP，利用最小二乘法求得k₁、b₁、k₂、b₂四个参数；得到分段回归模型参数之后，利用均方根误差 (Root Mean Squared Error，RMSE)评价所建立的模型对AMP升序曲线的拟合精度。

进一步地，所述的利用均方根误差(Root Mean Squared Error，RMSE)评价所建立的模型对AMP升序曲线的拟合精度

其中，RMSE_bp表示以bp为测试点(假设的突变点)时，所构建的分段回归模型的均方根误差；AMP_x为AMP升序曲线中第x个像元对应的AMP，AMP_x’为在构建的分段回归模型中x对应的AMP预测值，S为像元总数；基于上述方法，以AMP升序序列中每一个点作为待测试点，构建分段回归模型并计算RMSE；最后，我们选取所有测试点中，RMSE最小的测试点，作为实际突变点，而突变点所对应的AMP即为有效像元分类阈值(TD_AMP)；之后，在AOI的AMP数据二维矩阵中，逐像元判断AMP与TD_AMP之间关系，将AOI像元分为有效(叶片)和无效像元，即大于TD_AMP则为有效像元，反之亦然。

进一步地，所述的对筛选后的逐日Gcc均值时序曲线进行拟合，根据公式 (IV)，所述的根据拟合方程计算曲率变化率RCC，根据公式(V)：

其中，t是年积日，y(t)是第t天Gcc值，a和b是拟合参数，c+d是年内Gcc 最大值，d是Gcc初始值，z＝e^a+b·t；在1-9月，取RCC达到第一个局部最大值的时间为开始展叶期、第二个局部最大值为成熟叶片期。

进一步地，所述的RCC达到第一个局部最大值，具体判断方法如下：

对RCC曲线求取一节导度RCC’，RCC’第一次由正变负时所对应的RCC，即为第一个局部最大值；第二次由正变负时所对应的RCC，即为第二个局部最大值。

发明原理：利用数字相机植被指数和有效像元识别算法进行相机影像中植被叶片的识别，优化感兴趣区域(Area of Interest，AOI)植被指数均值的信噪比，进而提高地基遥感地表物候的提取精度。该方法主要有一下算法组成：利用数字相机获取的完整年每日多时相RGB波段影像，计算AOI逐像元逐日绿度植被指数(Green Chromatic Coordinate，Gcc)，获得AOI的Gcc三维矩阵；逐像元计算 Gcc时序曲线的振幅(AMP)，并对AOI所有像元的AMP矩阵进行排序，利用分段线性回归算法，识别AMP序列的突变点，作为区分植被叶片与非叶片像元的阈值；针对叶片像元，进行AOI有效像元Gcc逐日均值计算，再依据时序数据突变点识别关键物候期。

有益效果：与现有技术相比，本发明的一种基于数字相机影像有效像元的植被物候期提取方法，利用数字相机植被指数时序数据进行叶片有效像元自动识别，可根据数字相机影像生成的植被指数时序数据，对植被叶片像元和背景像元进行快速区分，实现影像感兴趣区有效像元的自动识别，相比于常规监督/非监督分类方法，人工干预更少、效率更高，且不受待分类影像时相差异对分类精度的影响。本发明通过有效像元的识别，可有效降低土壤、枝干等非叶片像元对感兴趣区植被指数均值的干扰，提高影像时序数据的信噪比，进而提高植被关键物候期提取精度。本发明利用数字相机时序影像，快速、准确的区分植被叶片有效像元和土壤、枝干等无效像元，提高数字相机影像信噪比。利用有效像元，计算影像感兴趣区植被指数均值时序数据，研究对象更加具体、真实、可靠，可进一步提高遥感地表物候识别精度，与观测的植被物候信息具有更佳的一致性。本方法可逐年识别影像有效像元，有效降低因自然扰动(如火灾)或人工管理(如修枝、砍伐)引起的地表景观变化对感兴趣区植被物候提取的影响，自动化程度高，对中低覆盖度落叶植被具有更好的适用性。利用MATLAB编写的程序可以有效处理海量的数据。

附图说明

图1为一种基于数字相机影像有效像元的植被物候期提取方法的流程图；

图2为观测站点相机影像视场与感兴趣区(AOI)示意图；

图3是Gcc振幅升序曲线与突变点(即有效像元阈值)示意图；

图4为AOI有效像元标识示意图；

图5为AOI有效像元检验样本点分布及子区域划分图。

具体实施方式

下面结合附图和具体实施方式来详细说明本发明：为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

本发明的一种利用数字相机植被指数和有效像元识别算法优化地基遥感地表物候参数提取的方法，用于数字相机影像中植被叶片有效像元识别，剔除影像感兴趣区域土壤、枝干等非叶片像元，进而降低目标植被绿度指数均值噪声，提取植被关键物候期，作为生态系统物候过程模拟与预测的基础。

本发明的技术方案主要包括以下步骤：

(1)选定影像感兴趣区。

根据数字相机观测影像，选取感兴趣区(Area of Interest，AOI)，AOI可为矩形或闭合多边形边框，须完整覆盖目标地物，并剔除天空、建筑物等明显非目标地物；记录AOI的影像坐标，用于时序影像提取。值得注意的是，如观测期间数字相机影像视场发生变化，须先对影像进行配准，以保证时序影像AOI空间范围一致。

(2)影像数据预处理。

设相机拍摄频率为n景/小时，每天拍摄时长为m小时，则针对每日多时相影像(n×m景)，利用红、绿、蓝三波段数据，逐景逐像元计算绿度指数Gcc(Green ChromaticCoordinate)，如公式(1)所示：

其中，i为观测日期，取值为1～365；j为每天中的第j景影像，取值为1～n×m； p和q分别为目标像元在影像中的行号和列号，用于确定像元位置；R、G、B分别为红、绿、蓝波段的影像像素值。之后，再对每天中每一个像元的n×m个Gcc 求取中值，作为当日该像元的Gcc值，进而得到一年中AOI逐日Gcc三维矩阵。在这一步骤中需要注意的是，考虑到太阳角度对植被冠层光谱反射的影响，建议选取每日上午9:00以后至下午3:00之前作为拍摄时间，且以正午12:00为中心时间，前后的拍摄时长相等，以体现每日Gcc中值的代表性。

(3)植被指数振幅计算。

针对AOI的Gcc三维矩阵，进行逐像元处理：提取每一像元Gcc时序数据，即该像元按时间排列的逐日Gcc；利用奇异谱分析(Singular Spectrum Analysis， SSA)方法对Gcc时序数据进行噪声去除、缺值插补；计算时序曲线生长季内的Gcc振幅(AMP)，采用Gcc夏季高值与春季低值之差表征振幅；经AOI逐像元计算，得到AOI中每一个像元的AMP，即AMP二维矩阵(或称AMP空间分布矩阵)。

(4)有效像元识别。

针对AOI的AMP二维矩阵，将所有像元的AMP进行排序，得到AMP升序曲线，其横坐标为像元的排列序号，纵坐标为对应像元的AMP。本发明中，基于分段回归模型，采用逐点测试的方式，识别AMP升序曲线的突变点。具体方法为：选取AMP升序曲线上的任一点为测试点bp，即假设的突变点，以该测试点构建分段回归模型，如公式(2)所示：

其中，x为AMP序列中像元的排列序号(即横坐标值)，x_bp为待测试点的排列序号(即横坐标值)，AMP_x为x所对应的AMP，k₁、b₁、k₂、b₂分别为待测试点前后两段线性回归方程的参数。根据序号x与对应的AMP，利用最小二乘法求得k₁、b₁、k₂、b₂四个参数。得到分段回归模型参数之后，利用均方根误差 (Root Mean Squared Error，RMSE)评价所建立的模型对AMP升序曲线的拟合精度。

其中，RMSE_bp表示以bp为测试点(假设的突变点)时，所构建的分段回归模型的均方根误差；AMP_x为AMP升序曲线中第x个像元对应的AMP，AMP_x’为在构建的分段回归模型中x对应的AMP预测值，S为像元总数。基于上述方法，以AMP升序序列中每一个点作为待测试点，构建分段回归模型并计算RMSE。最后，我们选取所有测试点中，RMSE最小的测试点，作为实际突变点，而突变点所对应的AMP即为有效像元分类阈值(TD_AMP)。之后，在AOI的AMP数据二维矩阵中，逐像元判断AMP与TD_AMP之间关系，将AOI像元分为有效(叶片)和无效像元，即大于TD_AMP则为有效像元，反之亦然。

(5)植被指数均值时序曲线计算。

为降低误分风险，对有效像元AMP再次进行排序，选取AMP大于AMP序列下四分位数的有效像元，为最终有效像元；对Gcc三维矩阵中每一天的所有最终有效像元的Gcc求平均，得到AOI有效像元逐日Gcc均值的时序曲线。

(6)遥感地表物候提取。

利用五参数Logistic方程，公式(4)对筛选后的逐日Gcc均值时序曲线进行拟合，并根据拟合方程计算曲率变化率RCC，公式(5)，对AOI有效像元逐日Gcc均值时序数据提取开始展叶期与成熟叶片期的物候日期，实现关键物候期识别。

其中，t是年积日，y(t)是第t天Gcc值，a和b是拟合参数，c+d是年内Gcc 最大值，d是Gcc初始值，z＝e^a+b·t。在1-9月，取RCC达到第一个局部最大值的时间为开始展叶期、第二个局部最大值为成熟叶片期。其中，为判断RCC局部最大值，可对RCC曲线求取一节导度(RCC’)，RCC’第一次由正变负时所对应的RCC，即为第一个局部最大值；第二次由正变负时所对应的RCC，即为第二个局部最大值。

实施例

选取江苏省宿迁市泗洪县陈圩林场的杨树人工林生态观测站，技术流程如图 1所示。泗洪站(118°36'E，33°32'N)植被为加拿大杨(Populus canadensis Moench)，2007年造林，土壤为淤积土、黏壤土。加拿大杨在我国分布广泛，是具有代表性的人工林树种。泗洪站于2016年4月建成，建有35米高通量塔一座，装备有多光谱物候相机(PhenoCam)、涡动相关通量观测系统及其他观测仪器，可以进行连续自动观测，设备运行稳定。

根据技术方案步骤(1)，对观测数据进行预处理。相机系统可拍摄每日多时相影像，影像大小为1296×960像素。本例中选取2019年逐日9:00-15:00时间范围内影像为数据源，每半小时1景，每日12景。选择多景不同日期且质量较好的影像作为相机系统视场的参考影像，确定感兴趣研究区(AOI)。为测试本发明所提出方法的效果，以2019年5月20日11:30拍摄的影像为例选取AOI(如图2所示)，并记录AOI行列号。AOI大小为400×550像素矩形区域，其中包括杨树叶片、树干和裸土像元。

根据技术方案步骤(2)，针对单日的12景影像的AOI，逐景逐像元计算Gcc，并计算每日Gcc中值，得到AOI逐日Gcc的三维矩阵。

根据技术方案步骤(3)，利用奇异谱分析(SSA)方法对Gcc三维矩阵进行逐像元去噪、插补处理，本例中，针对日尺度Gcc数据，SSA滤波窗口大小选为10天。针对预处理后的Gcc三维矩阵，本例中生长季振幅定义为夏季(6-8 月)Gcc中值与春季(3-4月)Gcc下四分位数的差，逐像元计算振幅，得到AOI 的AMP二维矩阵。

根据技术方案步骤(4)，进行AOI有效像元识别。将矩阵中220000个像元 AMP升序排列(图3)。利用逐点测试的分段回归模型方法，识别AMP升序曲线的突变点，并确定有效像元阈值。在本例中，被识别为突变点的像元排序编号为21653，AOI有效像元AMP阈值为0.0537(图3)。之后，逐像元判断AMP 与阈值关系，并在AOI中对有效像元与无效像元进行标记，完成有效像元的识别，结果如图4中 的 (a)所示。图4中 的 (a)中，标记为红色的像元即为被识别为冠层叶片的有效像元。

为验证有效像元识别精度，基于2019年5月20日11:30影像，利用目视判读方式在AOI内选择200个样本点，如图5中 的 (a)所示。图5中 的 (a)中红色的点即为选择的样本点，包括140个叶片像元、60个非叶片(裸土、枝干)像元。之后，利用Kappa系数(公式6-7)判断识别效。

其中，Coe_kappa为Kappa系数，p_o是每一类正确分类的样本数量之和除以总样本数；u₁、u₂分别为叶片像元和非叶片像元的本个数，v₁、v₂分别为算法识别出来的叶片像元和非叶片像元的样本个数，w为总样本个数。结果表明，表征叶片的有效像元识别率为133/140，非有效像元识别率为53/60，Kappa系数为0.8333，识别精度较高。为与本发明的方法对比，基于该日期影像的RGB段数据，利用K-Means聚类算法对AOI进行自动聚类，类别数量设置为2，识别结果如图4中 的 (b) 所示，其中标记为红色的像元即为被识别为冠层叶片的有效像元。结果表明，基于K-Means聚类算法的有效像元识别率为80/140，非有效像元识别率为40/60， Kappa系数为0.2，识别效果一般。此外，考虑到K-Means聚类结果与选择的影像时间与质量有关，在林冠未达到郁闭时，分类结果不具有参考价值。因此，本发明建议的方法，具有更高精度，也更广泛的适用性。

在本例中，除整个AOI外，我们分别选取2个具有代表性的子感兴趣区 (sub-AOI，如图5中 的 (b)所示，其中红框范围即为子感兴趣区域)，用于验证在近单木尺度上，本发明方法提取植被物候的优势。根据技术方案步骤(5)，分别求取AOI和sub-AOI有效像元日均Gcc时序曲线，以及AOI和sub-AOI全部像元的日均Gcc时序曲线。

最后，根据技术方案步骤(6)，利用五参数Logistic方程(公式4)，对Gcc 曲线进行拟合，并根据拟合方程计算曲率变化率RCC(公式5)，对考虑和不考虑有效像元的AOI和sub-AOI日均Gcc时序曲线提取开始展叶期与成熟叶片期的物候日期。本例中，参考中国科学院地理科学与资源研究所2010年12月发布的《中国物候观测规范质量要求与观测报表》，对2019年研究区关键物候期进行目视观测，其中开始展叶期约为4月6日、成熟叶片期约为4月22日，并与遥感提取的关键物候期进行比较(表1)。结果表明，顾及影像有效像元时所提取的开始展叶期与成熟叶片期平均误差分别为4.3±0.9天和3.3±0.9，相比于传统的 AOI平均方法，本发明方法的物候结果更接近真是物候期，特别是在土壤背景占比较大的近单木尺度(图5中 的 (b)-②)。存在土壤像元会导致物候期提前，主要是因为土壤像元区域在春季有草本植物生长：草本植物返青早于乔木，因此初春时土壤像元区Gcc值高于有效像元区，因此整个AOI的开始展叶期早于有效像元区的开始展叶期；而随着乔木快速生长，林冠逐渐郁闭，对草地生长造成水分和辐射限制，造成草地Gcc峰值早于林冠Gcc峰值，因此整个AOI的成熟叶片期早于有效像元区的成熟叶片期。

根据本发明提出的方法，可获取高精度的地基遥感地表物候，用于植被生理生态过程监测，也可用于卫星遥感地表物候尺度转化，或其他植被遥感方面的研究，该方法可实现跨尺度迁移，通用性强。适用于季节性节律植被，特别是中低郁闭度生态系统，通过对植被有效像元的识别，精确监测植被冠层生长与衰落过程，提取植被关键物候期，结合卫星遥感地表物候，可直接用于改进地球系统模型，提高陆地生态系统碳、水通量估算的精度。

表1本发明方法与不考虑有效像元物候提取方法的结果比较

物候期研究区	开始展叶期	成熟叶片期
			AOI	3月30日	4月17日
有效AOI	4月1日	4月18日
			sub-AOI①	3月29日	4月15日
有效sub-AOI①	4月1日	4月18日
			sub-AOI②	3月25日	4月12日
有效sub-AOI②	4月3日	4月20日

Claims (5)

1.一种基于数字相机影像有效像元的植被物候期提取方法，其特征在于，包括如下步骤：

1)选定影像感兴趣区

根据数字相机观测影像，选取感兴趣区AOI；

2)影像数据预处理

设相机拍摄频率为n景/小时，每天拍摄时长为m小时，则针对每日多时相影像n×m景，利用红、绿、蓝三波段数据，逐景逐像元计算绿度指数Gcc，绿度指数Gcc，如公式(I)所示：

其中，i为观测日期，取值为1～365；j为每天中的第j景影像，取值为1～n×m；p和q分别为目标像元在影像中的行号和列号，用于确定像元位置；R、G、B分别为红、绿、蓝波段的影像像素值；再对每天中每一个像元的n×m个Gcc求取中值，作为当日该像元的Gcc值，进而得到一年中AOI逐日Gcc三维矩阵；

3)植被指数振幅计算

针对AOI的Gcc三维矩阵，进行逐像元处理；

4)有效像元识别

针对AOI的AMP二维矩阵，将所有像元的AMP进行排序，得到AMP升序曲线，其横坐标为像元的排列序号，纵坐标为对应像元的AMP；基于分段回归模型，采用逐点测试的方式，识别AMP升序曲线的突变点，具体方法为：选取AMP升序曲线上的任一点为测试点bp，即假设的突变点，以该测试点构建分段回归模型，如公式(II)所示：

其中，x为AMP序列中像元的排列序号，x_bp为待测试点的排列序号，AMP_x为x所对应的AMP，k₁、b₁、k₂、b₂分别为待测试点前后两段线性回归方程的参数；根据序号x与对应的AMP，利用最小二乘法求得k₁、b₁、k₂、b₂四个参数；得到分段回归模型参数之后，利用均方根误差RMSE评价所建立的模型对AMP升序曲线的拟合精度；

所述的利用均方根误差(Root Mean Squared Error，RMSE)评价所建立的模型对AMP升序曲线的拟合精度

其中，RMSE_bp表示以bp为测试点时，所构建的分段回归模型的均方根误差；AMP_x为AMP升序曲线中第x个像元对应的AMP，AMP_x’为在构建的分段回归模型中x对应的AMP预测值，S为像元总数；以AMP升序序列中每一个点作为待测试点，构建分段回归模型并计算RMSE；最后，选取所有测试点中，RMSE最小的测试点，作为实际突变点，而突变点所对应的AMP即为有效像元分类阈值TD_AMP；在AOI的AMP数据二维矩阵中，逐像元判断AMP与TD_AMP之间关系，将AOI像元分为有效和无效像元，即大于TD_AMP则为有效像元，反之亦然；

所述的对筛选后的逐日Gcc均值时序曲线进行拟合，根据公式(IV)；所述的根据拟合方程计算曲率变化率RCC，根据公式(V)：

其中，t是年积日，y(t)是第t天Gcc值，a和b是拟合参数，c+d是年内Gcc最大值，d是Gcc初始值，z＝e^a+b·t；在1-9月，取RCC达到第一个局部最大值的时间为开始展叶期、第二个局部最大值为成熟叶片期；

5)植被指数均值时序曲线计算

为降低误分风险，对有效像元AMP再次进行排序，选取AMP大于AMP序列下四分位数的有效像元，为最终有效像元；对Gcc三维矩阵中每一天的所有最终有效像元的Gcc求平均，得到AOI有效像元逐日Gcc均值的时序曲线；

6)遥感地表物候提取

利用五参数Logistic方程，对筛选后的逐日Gcc均值时序曲线进行拟合，并根据拟合方程计算曲率变化率RCC，对AOI有效像元逐日Gcc均值时序数据提取开始展叶期与成熟叶片期的物候日期，实现关键物候期识别。

2.根据权利要求1所述的一种基于数字相机影像有效像元的植被物候期提取方法，其特征在于，所述的步骤1)具体为：AOI为矩形或闭合多边形边框，须完整覆盖目标地物，并剔除非目标地物；记录AOI的影像坐标，用于时序影像提取；观测期间数字相机影像视场发生变化，先对影像进行配准，保证时序影像AOI空间范围一致。

3.根据权利要求1所述的一种基于数字相机影像有效像元的植被物候期提取方法，其特征在于，所述的步骤2)中，选取每日上午9:00以后至下午3:00之前作为拍摄时间，且以正午12:00为中心时间，前后的拍摄时长相等，以体现每日Gcc中值的代表性。

4.根据权利要求1所述的一种基于数字相机影像有效像元的植被物候期提取方法，其特征在于，所述的步骤3)中，所述的进行逐像元处理具体为：提取每一像元Gcc时序数据，该像元按时间排列的逐日Gcc；利用奇异谱分析SSA方法对Gcc时序数据进行噪声去除、缺值插补；计算时序曲线生长季内的Gcc振幅AMP，采用Gcc夏季高值与春季低值之差表征振幅；经AOI逐像元计算，得到AOI中每一个像元的AMP，即AMP二维矩阵。

5.根据权利要求1所述的一种基于数字相机影像有效像元的植被物候期提取方法，其特征在于，所述的RCC达到第一个局部最大值，具体判断方法如下：

对RCC曲线求取一节导度RCC’，RCC’第一次由正变负时所对应的RCC，即为第一个局部最大值；第二次由正变负时所对应的RCC，即为第二个局部最大值。

Images