A method for extracting rubber forest leaf fall and compound leaf phenology based on kNDVI daily time series

Abstract

The present application relates to a kind of rubber forest leaf fall and compound leaf phenology extraction method based on kNDVI daily time series, comprising the following steps: S1, based on MODIS daily data constructs kNDVI time series and carries out cloud removal processing;S2, to the kNDVI time series after cloud removal, serious interpolation in time and space scale is carried out to the missing data, to ensure the integrity of time series;S3, to the kNDVI time series after interpolation, envelope line SG filtering is carried out, and the kNDVI time series after interpolation is smoothed;S4, for the kNDVI time series after processing of step S3, season amplitude method is used to extract rubber forest phenology information, obtains the four phenological indicators of rubber forest leaf fall beginning, leaf fall end, compound leaf beginning and compound leaf end.The method can accurately extract large-scale long time of rapidly changing rubber forest leaf fall and compound leaf phenology.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of phenology monitoring, and particularly relates to a rubber forest leaf fall and compound leaf phenology extraction method based on kNDVI daily time series. BACKGROUND

[0002] Rubber trees can produce latex after 5-7 years of growth, and rubber, as one of the four industrial raw materials, is an important resource-constrained strategic material, which is of great significance to national industrial development and strategic security. Rubber trees are evergreen broad-leaved species in the Amazon region, but after being introduced to Asia, they show a concentrated leaf fall phenology in winter from December to January, with a very low canopy density, a leafless period of about 2-4 weeks, and then enter the leafing period before the rainy season in February-April, with a rapid recovery of canopy. However, so far, no one has used remote sensing to retrieve the fall and spring phenology of rubber trees in large areas, and there is a lack of in-depth research.

[0003] The current means of studying phenology mainly include ground phenology observation network, near-ground observation, phenology model and remote sensing data. Ground phenology observation network is a first-hand data with reliable data quality, but its spatial distribution is uneven and there is no unified monitoring standard, which cannot guarantee the data quality. Near-ground observation such as unmanned aerial vehicles and phenology cameras can provide higher spatial and temporal resolution data compared with ground observation network, but their number and coverage are limited. In contrast, remote sensing data has global coverage and temporal continuity, which has a significant advantage in studying large-scale forest phenology and is most recognized. In the past few decades, satellite remote sensing data has been used to monitor vegetation phenology in large areas through various types of observations and methods. Most of the phenology monitoring research is through the construction of vegetation index time series.

[0004] Rubber forests are mostly planted in tropical and subtropical regions, which are cloudy and rainy, and remote sensing data is easily affected by weather, resulting in missing data and excessive noise in time series. At present, due to the limitation of remote sensing data quality, most of the literature is based on Landsat data, Sentinel data, MODIS data and other rough time resolution data for rubber classification research, such as "Phenology-Based Vegetation Index Differencing for Mapping of Rubber Plantations Using Landsat OLI Data". In the aspect of rubber forest phenology research, the existing literature is mainly based on 16-day resolution Landsat data and 8-day, 16-day resolution MODIS data, such as "Remote Sensing Monitoring of Spring Phenology of Natural Rubber Forests in Hainan Island" and based on site human observation phenology data for analysis of phenology and climate, such as "Responses of rubber leaf phenology to climatic variations in Southwest China". However, time series with such rough time resolution cannot accurately capture the rapid occurrence of leaf fall phenology start and end and compound leaf phenology start and end in rubber forests within about 2-4 weeks. Based on 16-day MODIS time series, the maximum error of 60 days per year and the average error of 13.5 days per year were found in the extraction of spring compound leaf phenology in Hainan Island, which verified the human observation phenology. It can be seen that the rough resolution time series has a great influence on the accuracy of phenology extraction. Therefore, how to construct high time resolution remote sensing time series such as daily scale to accurately extract the rapid change of rubber forest phenology is a technical difficulty that needs to be solved at present. SUMMARY

[0005] The purpose of the present application is to provide a method for extracting rubber forest leaf fall and compound leaf phenology based on kNDVI daily time series, which can accurately extract the rapid change of rubber forest leaf fall and compound leaf phenology on a large scale and for a long time.

[0006] In order to achieve the above purpose, the technical scheme adopted by the present application is: a method for extracting rubber forest leaf fall and compound leaf phenology based on kNDVI daily time series, comprising the following steps:

[0007] S1, constructing kNDVI time series based on MODIS daily data and performing cloud removal processing;

[0008] S2, interpolating the kNDVI time series with serious missing data after cloud removal in time and space scales to ensure the integrity of the time series;

[0009] S3, envelope SG filtering is conducted on the interpolated kNDVI time series, and the interpolated kNDVI time series is smoothed;

[0010] S4, for the kNDVI time series processed in step S3, the seasonal amplitude method is used to extract the rubber forest phenology information, and four phenology indexes of rubber forest leaf fall start, leaf fall end, leaf regrowth start and leaf regrowth end are obtained.

[0011] Further, in step S1, the kNDVI time series is constructed based on the MOD09GQ and MOD09GA daily data and cloud removal processing is conducted; first, the sur_refl_b01_1 red (R) band and sur_refl_b02_1 near infrared (NIR) band two bands are extracted from the MOD09GQ data, and the daily kernel normalized vegetation index (kNDVI) of the study area is calculated, and the calculation formula is as follows;

[0012]

[0013] k(a,b)=exp(-(a-b) 2 / (2σ 2 )) (2)

[0014] Wherein, R, NIR are sur_refl_b01_1 band and sur_refl_b02_1 band of MOD09GQ image data; k(NIR, NIR), k(NIR, R) in formula (1) are calculated according to formula (2); σ is the length scale parameter, which is taken as the average of R and NIR two bands;

[0015] Then, the kNDVI daily time series calculated is subjected to cloud removal processing, and the MOD09GA quality control band is used to remove noise and retain good quality pixels; the main steps are: the 1 kilometer resolution daily quality evaluation flag is extracted from the MOD09GA product state_1km_1 band, which is resampled to 250m resolution, and the daily kNDVI data calculated is subjected to pixel quality screening according to the daily reliable pixels.

[0016] Further, in step S2, for the rubber forest target pixel affected by noise without surface reflection data, first, the kNDVI mean value of the adjacent date of the same pixel is used for replacement, which is specifically shown in formula (3);

[0017] Rubber(x,y,day)=mean(Rubber(x,y,day-1),Rubber(x,y,day+1)) (3)

[0018] Wherein, x, y are the row and column numbers of the target pixel, and day is the date;

[0019] If the target pixel of the rubber forest is affected by continuous noise for several days, and there is still no data after interpolation, further interpolation is performed using the mean value of the pixels in the same rubber square buffer zone, as shown in equation (4);

[0020]

[0021] Since the MOD09GQ data resolution used is 250 m, the pixel range of the square buffer zone is set to 10, 20, 30, and 40 pixels above and below the center pixel, i.e., 2.5 km, 5 km, 7.5 km, and 10 km, and the maximum buffer zone is set to 10 km. According to geographical knowledge, the smaller the buffer zone, the more similar the climate, geographical conditions, and vegetation characteristics. Therefore, if the 2.5 km buffer zone mean value has a value, it is replaced with the 2.5 km buffer zone mean value. If the 2.5 km buffer zone mean value has no value and the next level 5 km buffer zone mean value has a value, it is replaced with the 5 km buffer zone mean value, and so on, as shown in equation (5);

[0022]

[0023] Further, in step S3, first, the envelope line of each target rubber forest pixel is extracted year by year kNDVI time series, and the time series is reconstructed, as shown in equation (6);

[0024] Rubber_TS(x,y,day)=max(Rubber_TS(x,y,day-w:day+w)) (6)

[0025] where w is the envelope line sliding window value. If the kNDVI daily data is missing for several consecutive days, and after envelope line reconstruction, the time series still has missing kNDVI values for some dates, linear interpolation is used to interpolate the missing values in the time series to make the time series continuous and complete, as shown in equation (7);

[0026]

[0027] where Rubber_TS(x,y,day) is the kNDVI time series of the target rubber forest pixel on day day after moving window reconstruction, day1 is the first date with kNDVI value after day day of the target pixel, and day0 is the first date with kNDVI value before day day of the target pixel;

[0028] Then, the SG filter is used to smooth the time series of kNDVI for the whole year. The rubber forest phenology mainly occurs from December of the previous year to April of the current year. In order to smooth and eliminate the influence of non-phenological period, the mean value smoothing is performed on the kNDVI data of the target rubber forest pixel in the non-phenological period, so as to facilitate the fitting of the subsequent phenological curve and the extraction of the phenological indicators. The specific implementation is shown in formula (8).

[0029] Rubber_TS(x,y,day Oct.～Nov.,May )=mean(Rubber_TS(x,y,day Dec.～Apr. )) (8)

[0030] Rubber_TS(x,y,day Oct.～Nov.,May ) is the daily kNDVI value of the target pixel in the non-rubber forest phenological period (i.e. October and November of the previous year and May of the current year) of a certain year, and Rubber_TS(x,y,day Dec.～Apr. ) is all daily kNDVI values of the target pixel in the rubber forest phenological period (i.e. from December of the previous year to April of the current year) of a certain year. The mean value of kNDVI of all dates in the rubber forest phenological period is used to replace the daily kNDVI value in the non-phenological period, so as to achieve the purpose of smoothing the phenological curve.

[0031] Further, in step S4, the four phenological indicators are defined as follows:

[0032] 1) Beginning of defoliation: it refers to the turning point when the rubber forest enters the rapid defoliation from the stable growth period. The vegetation index value decreases, indicating that the greenness decreases, which means that the rubber forest enters the defoliation period.

[0033] 2) End of defoliation: it refers to the point when the vegetation index value of the rubber forest decreases to the minimum and no longer decreases, which means that the defoliation period has ended.

[0034] 3) Beginning of germination: it refers to the point when the vegetation index of the rubber forest starts to increase after a short period of stability, which means that the rubber forest enters the leaf regrowth period.

[0035] 4) End of leaf expansion / peak leaf period: it refers to the point when the vegetation index value of the rubber forest increases to a certain value and no longer increases, which means that the rubber forest enters a relatively stable stage, and the leaf is in the peak period, and the vegetation index value is high.

[0036] For the kNDVI daily time series of the target rubber forest pixel after step S3, the asymmetric Gaussian function fitting (AG) algorithm is used for curve fitting, and then the seasonal amplitude method is used to extract the four phenological indicators of the rubber forest in each year.

[0037] Compared with the prior art, the present application has the following beneficial effects: the present application provides a rubber forest leaf fall and compound leaf phenology extraction method based on kNDVI daily time series, which is based on remote sensing and geographic knowledge to construct a reasonable and complete MODIS kNDVI daily time series, and can accurately extract large-scale long-time rubber forest leaf fall and compound leaf phenology which changes rapidly, significantly improves the accuracy of remote sensing data in extracting rubber forest phenology, solves the problems of how to reasonably interpolate the missing MODIS daily time series of rubber forest in the multi-cloud and rainy tropical region, the insufficient accuracy of rubber forest phenology extracted by the existing remote sensing time series, and the incomplete spatial and temporal continuous extraction of rubber forest leaf fall and compound leaf phenology, and the like. Therefore, the present application has strong practicability and broad application prospect. BRIEF DESCRIPTION OF DRAWINGS

[0038] Figure 1 is a method implementation flowchart of an embodiment of the present application;

[0039] Figure 2 is a rubber forest kNDVI time series phenology extraction schematic diagram in an embodiment of the present application. DETAILED DESCRIPTION

[0040] The present application will be further described below in combination with the drawings and embodiments.

[0041] It should be noted that the following detailed description is exemplary and is intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the present application belongs.

[0042] It should be noted that the terms used herein are only for the purpose of describing specific embodiments, and are not intended to limit the exemplary embodiments according to the present application. As used herein, unless the context clearly indicates otherwise, the singular form is intended to include the plural form, and in addition, it should be understood that when the terms "comprise" and / or "include" are used in the specification, they indicate the presence of a feature, step, operation, device, component and / or combination thereof.

[0043] As shown in Figure 1 The present embodiment provides a rubber forest leaf fall and compound leaf phenology extraction method based on kNDVI daily time series, which comprises the following steps:

[0044] S1, constructing kNDVI time series based on MODIS daily data and performing cloud removal processing.

[0045] Because the rubber forest is mostly planted in tropical areas, the remote sensing surface radiation data is easily affected by clouds and rainy seasons, and the remote sensing data noise is more, which further affects the quality and reliability of the daily kNDVI data calculated, so it is necessary to remove the cloud according to the MODIS product QA quality band, and to remove the unreliable pixel value.

[0046] S2, the kNDVI time series after cloud removal with serious missing data is interpolated in time and space scale to ensure the integrity of the time series.

[0047] Because the kNDVI daily time series of the rubber forest after cloud removal has a serious missing rate, most of the pixels have a missing rate of more than 50% for about 240 days from October 1 of the previous year to May 31 of the current year, that is, nearly 120 days of remote sensing data are affected by clouds and rainy seasons, and in order to build a more complete time series and more accurately extract the leaf fall and leaf recovery phenology of the rubber forest, it is extremely necessary to interpolate the rubber forest pixels with serious missing data. First, for a rubber forest target pixel with missing data on a certain day, search for whether there is kNDVI data on the previous day and the next day, if there is, interpolate by adjacent dates; because the missing is too serious, only time interpolation cannot extract the rubber forest phenology, therefore, based on the annual distribution of rubber, for the target pixel with missing data, the mean value of all rubber forest pixels with kNDVI value within the 2.5km square buffer area centered on the target pixel is searched for interpolation, if there is still no data, continue to search for the mean value of the 5km, 7.5km and 10km square buffer area. Through testing, the missing rate of the time series can be reduced to about 20% through superposition of time and space interpolation, which is enough to build a more complete kNDVI time series of the rubber forest and extract the phenology.

[0048] S3, the envelope SG filtering of the interpolated kNDVI time series is carried out to smooth the interpolated kNDVI time series.

[0049] The envelope SG filtering of the kNDVI time series can better smooth the rubber forest vegetation index time series and better reflect the rubber forest phenology, therefore, the envelope SG filtering smoothing curve is adopted for the time series in the method. The non-phenological period mean value is used for smoothing, so as to eliminate the interference of non-phenological period noise and be beneficial to the fitting of the phenology curve.

[0050] S4, for the kNDVI time series after step S3, the seasonal amplitude method is used to extract the rubber forest phenology information and obtain four phenological indexes of the rubber forest, including the beginning of leaf fall, the end of leaf fall, the beginning of leaf recovery and the end of leaf recovery.

[0051] The method is based on programming batch reading of the kNDVI time series of the rubber forest in each year, and fitting of the phenology curve for each pixel, and then extracting the rubber forest phenological indexes in each year based on the seasonal amplitude method.

[0052] The related technical content involved in the implementation process of the method will be further described in combination with specific embodiments.

[0053] The rubber forest leaf fall and compound leaf phenology extraction method based on kNDVI daily time series provided in the embodiment is specifically executed as follows in actual application: steps S1 to S4.

[0054] S1, obtain the MODIS kNDVI vegetation index daily remote sensing image time series corresponding to each year of the target area in the main occurrence period of rubber phenology (i.e. from October of the previous year to May of the current year), and perform corresponding cloud removal processing according to the remote sensing quality image to eliminate noise pixels. Then go to step S2.

[0055] The data used by the method is the MOD09GQ (Terra) land surface reflection data set downloaded free of charge by NASA USGS data center, with a time resolution of daily and a spatial resolution of 250 meters. Since the rubber forest phenology in the study area mainly occurs from November to the following April, the time range of the remote sensing data of the method is from October of the previous year to May of the current year each year, and there are nearly 240 images (8 months) per year, which are continuous time images. The data set is an atmospheric correction and radiation correction data product that can be directly cropped for use.

[0056] First, extract the sur_refl_b01_1 red (R) band and sur_refl_b02_1 near-infrared (NIR) band from the MOD09GQ data, and batch calculate the kernel normalized vegetation index (kNDVI) of the study area each day, the calculation formula is as follows:

[0057]

[0058] k(a,b)=exp(-(a-b) 2 / (2σ 2 )) (2)

[0059] Where R and NIR are the sur_refl_b01_1 band (620-670 nm) and sur_refl_b02_1 band (841-876 nm) of the MOD09GQ image data. k(NIR, NIR) and k(NIR, R) in formula (1) are calculated according to formula (2); σ is the length scale parameter, and the reasonable value is the average of the two bands a and b.

[0060] Then, the results are de-clouded by the MOD09GA quality control band of the corresponding date to ensure the reliability of the calculated kNDVI values. According to the product guide, MOD09GQ should be used in combination with MOD09GA product. MOD09GA stores important quality information corresponding to the MOD09GQ surface reflection band. Therefore, the 1-kilometer resolution daily quality assessment (QA) flag is extracted from the MOD09GA product state_1km_1 band, which is resampled to 250m resolution. This band stores information about the atmospheric state (clouds, cloud shadows, cirrus clouds) of the corresponding data product of the surface reflection band, based on which the pixel quality of the MOD09GQ surface reflection data is screened. Using the MOD09GA quality control band, data affected by cloud, snow, aerosol and other noise in the surface reflection band are removed, and good quality pixels are retained.

[0061] S2, the time series with missing data is spatiotemporal scale interpolated to ensure the integrity of the time series. Then step S3 is entered.

[0062] For missing pixels affected by noise such as clouds, first, time interpolation is performed using adjacent dates, and then the mean value of all rubber forest pixels within the 50m elevation range and within the spatial buffer range (2.5km, 5km, 7km, 10km, the search range is continuously expanded until there is a replacement value) is searched according to the annual rubber forest classification data, and the interpolation is controlled within a reasonable range by twice the standard deviation.

[0063] Assuming that the target rubber forest pixel has no surface reflection data due to weather effects on a certain day, the mean value of the surface reflection data with good quality of the same pixel on adjacent dates is used to replace it, as shown in formula (3);

[0064] Rubber(x,y,day)=mean(Rubber(x,y,day-1),Rubber(x,y,day+1)) (3)

[0065] Where x, y are the row and column numbers of the target pixel, and day is the date.

[0066] If the target pixel of the rubber forest is affected by continuous noise such as weather for several days, and there is still no data after adjacent interpolation, the mean value of the pixels within the same rubber forest in the square buffer is considered for interpolation, as shown in formula (4);

[0067]

[0068] In the method, the MOD09GQ data used has a resolution of 250 m, and the square buffer pixel range is set to 10, 20, 30 and 40 pixels, i.e. 2.5 km, 5 km, 7.5 km and 10 km, to the left and right of the target pixel, and the maximum buffer is set to 10 km, and it is considered that the rubber forest physiological and geographical characteristics within the range of 10 km are similar. According to geographical knowledge, the closer to the buffer, the more similar, if the 2.5 km buffer mean value has a value, the 2.5 km buffer mean value is replaced, if the 2.5 km buffer mean value has no value and the 5 km buffer mean value has a value, the 5 km buffer mean value is replaced, and so on, as shown in formula (5);

[0069]

[0070] S3, the envelope line SG filtering method (ED-SG) is used to smooth the time series after interpolation, and reliable phenological information can be extracted. Then step S4 is entered.

[0071] The kNDVI time series envelope line SG filtering method (Envelope Detection and the Savitzky-Golay filter, ED-SG) is used to smooth the time series curve, and then the data of October and November of the previous year and the data of May of the current year are smoothed, and the mean value is taken as the mean value of all dates of kNDVI value in the period from December 1 of the previous year to April 30 of the current year, so as to ensure that the time series has no noise interference outside the phenological period, and the phenological date can be better extracted.

[0072] Firstly, the envelope line of the kNDVI time series of each target rubber forest pixel is extracted year by year, and the time series is reconstructed, as shown in formula (6);

[0073] Rubber_TS(x, y, day) = max(Rubber_TS(x, y, day-w:day+w)) (6)

[0074] In the method, the MOD09GQ data used has a resolution of 250 m, and the square buffer pixel range is set to 10, 20, 30 and 40 pixels, i.e. 2.5 km, 5 km, 7.5 km and 10 km, to the left and right of the target pixel, and the maximum buffer is set to 10 km, and it is considered that the rubber forest physiological and geographical characteristics within the range of 10 km are similar. According to geographical knowledge, the closer to the buffer, the more similar, if the 2.5 km buffer mean value has a value, the 2.5 km buffer mean value is replaced, if the 2.5 km buffer mean value has no value and the 5 km buffer mean value has a value, the 5 km buffer mean value is replaced, and so on, as shown in formula (5);

[0075]

[0076] Rubber_TS(x,y,day) is the kNDVI time series of the target rubber forest pixel on the day after the moving window reconstruction, day1 is the date of the first kNDVI value after the day of the target pixel, and day0 is the date of the first kNDVI value before the day of the target pixel.

[0077] Then, the SG filter is applied to the kNDVI time series which has been completely continuous, and the time series is smoothed. The rubber forest phenology mainly occurs from December of the previous year to April of the current year. In order to smooth and eliminate the influence of the non-phenological period, the kNDVI data of the target rubber forest pixel in the non-phenological period is smoothed by mean value, so as to facilitate the fitting of the subsequent phenological curve and the extraction of the phenological index. The specific implementation is shown in formula (8).

[0078] Rubber_TS(x,y,day Oct.～Nov.,May )=mean(Rubber_TS(x,y,day Dec.～Apr. )) (8)

[0079] Rubber_TS(x,y,day Oct.～Nov.,May ) is the daily kNDVI value of the target pixel in the non-rubber forest phenological period (i.e. October and November of the previous year and May of the current year) of a year, and Rubber_TS(x,y,day Dec.～Apr. ) is all daily kNDVI values of the target pixel in the rubber forest phenological period (i.e. from December of the previous year to April of the current year) of a year. The mean value of kNDVI of all dates in the rubber forest phenological period is used to replace the daily kNDVI value in the non-phenological period, so as to achieve the purpose of smoothing the phenological curve.

[0080] S4, the rubber forest phenological information is extracted based on the seasonal amplitude method, and four phenological indexes of the rubber forest, i.e. the beginning of leaf fall, the end of leaf fall, the beginning of leaf recovery and the end of leaf recovery, are obtained. In order to accurately describe the phenological change of the rubber forest, the present application adopts the key phenological node as the phenological characteristic parameter to analyze the rubber forest phenology. The rubber forest mainly has two phenological periods (leaf fall period and leaf recovery period). In a complete growth cycle of the rubber forest, the present method defines four phenological indexes, as shown in formula (9), including: Figure 2

[0081] 1) The beginning of leaf fall is the turning point of the rubber forest from the stable growth period to the rapid leaf fall period, the vegetation index value decreases, indicating that the greenness decreases, and it is indicated that the rubber forest enters the leaf fall period.

[0082] 2) The end of leaf fall is that the vegetation index value of the rubber forest decreases to the minimum and no longer decreases, indicating that the leaf fall period has ended.

[0083] ​3) Germination start: refers to the vegetation index of rubber forest after a short time of stabilization, starts to rise, that is, the rubber forest enters the compound leaf period.

[0084] 4) End of leaf expansion / peak leaf period: refers to the vegetation index value of rubber forest rising to a certain value and no longer rising, entering a relatively stable stage, indicating that the rubber forest compound leaf ends and the leaf is in the peak period, and the vegetation index value is high.

[0085] For the target rubber forest pixel kNDVI daily time series after step S3 processing, first, an asymmetric Gaussian function fitting (AG) algorithm is used for curve fitting, and then a seasonal amplitude method is used to extract the phenological index of the rubber forest. After testing, the threshold values of 10% for the start of leaf fall and the start of compound leaf, and 90% for the end of leaf fall and the end of compound leaf are taken. According to the above process, the annual rubber forest leaf fall and compound leaf phenology can be extracted.

[0086] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0087] The present application is described with reference to flowcharts and / or block diagrams according to the methods, devices (systems), and computer program products of the embodiments of the present application. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, and the combination of flows and / or blocks in the flowcharts and / or block diagrams can be implemented by computer program instructions. These computer program instructions can be provided to a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing apparatus to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing apparatus produce a device that implements the functions specified in the flowcharts and / or block diagrams. Figure 1 The functions specified in one or more flows and / or blocks Figure 1 The means for performing the functions specified in one or more blocks or flows.

[0088] These computer program instructions can also be stored in a computer-readable memory that can cause the computer or other programmable data processing apparatus to work in a specific manner, so that the instructions stored in the computer-readable memory produce a manufactured product including instruction means, which implements the functions specified in the flowcharts and / or block diagrams. Figure 1 The functions specified in one or more flows and / or blocks Figure 1 The means for performing the functions specified in one or more blocks or flows.

[0089] These computer program instructions can also be loaded into a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 1 one or more flowcharts and / or blocks Figure 1 one or more flowcharts and / or blocks

[0090] The above descriptions are only the preferred embodiments of the present application, not intended to limit the present application to other forms described. Any person skilled in the art may make changes or modifications to the above-described technical contents as equivalent embodiments without departing from the technical solutions of the present application. However, any simple modification, equivalent change and modification of the above embodiments according to the technical essence of the present application without departing from the technical solutions of the present application still belongs to the protection scope of the technical solutions of the present application.

Claims

1. A method for extracting rubber forest deciduous leaf and compound leaf phenology based on kNDVI daily time series, characterized by, The method comprises the following steps: S1, constructing kNDVI time series based on MODIS daily data and performing cloud removal processing; S2, interpolating the kNDVI time series after cloud removal in the time and space scales to ensure the integrity of the time series; S3, performing envelope line SG filtering on the interpolated kNDVI time series to smooth the interpolated kNDVI time series; S4, for the kNDVI time series processed in step S3, extracting rubber forest phenology information by using the seasonal amplitude method to obtain four phenology indexes of rubber forest, i.e., leaf fall start, leaf fall end, leaf regrowth start and leaf regrowth end; In step S2, for a target pixel of the rubber forest affected by noise without ground reflection data on a certain day, the mean value of the kNDVI of the adjacent dates of the same pixel is used for replacement, as shown in formula (3); Rubber(x,y,day)=mean(Rubber(x,y,day-1),Rubber(x,y,day+1)) (3) Wherein, x and y are the row and column numbers of the target pixel, and day is the date; If the target pixel of the rubber forest is affected by noise for consecutive days and still has no data after adjacent interpolation, the mean value of the pixels of the rubber forest in the square buffer is used for interpolation, as shown in formula (4); Wherein, since the resolution of the MOD09GQ data used is 250m, the range of the pixels in the square buffer is set to 10, 20, 30 and 40 pixels above, below, left and right of the center pixel, i.e., 2.5km, 5km, 7.5km and 10km, and the maximum buffer is 10km; according to geographical knowledge, the smaller the buffer, the more similar the climate, geographical conditions and vegetation characteristics; therefore, if the mean value of the 2.5km buffer has a value, the mean value of the 2.5km buffer is used for replacement, if the mean value of the 5km buffer has a value, the mean value of the 5km buffer is used for replacement, and so on, as shown in formula (5); In step S3, first, the envelope line of the kNDVI time series of each target rubber forest pixel is extracted year by year, and the time series is reconstructed, as shown in formula (6); Rubber_TS(x,y,day)=max(Rubber_TS(x,y,day-w:day+w)) (6) Wherein, w is the value of the envelope line sliding window; if the kNDVI daily data is missing for consecutive days, and the time series still has missing kNDVI values after envelope line reconstruction, linear interpolation is used to interpolate the missing values of the time series to make the time series continuous and complete, as shown in formula (7); IF Rubber_TS(x,y,day)=novalue, Wherein, Rubber_TS(x,y,day) is the kNDVI time series of the target rubber forest pixel on day day after moving window reconstruction, day1 is the first date with kNDVI value after day day of the target pixel, and day0 is the first date with kNDVI value before day day of the target pixel; Then, the SG filter is used to smooth the time series of kNDVI. The rubber forest phenology mainly occurs from December of the previous year to April of the current year. In order to smooth and eliminate the influence of non-phenological period, the mean value smoothing is performed on the kNDVI data of the target rubber forest pixel in the non-phenological period, so as to facilitate the fitting of the subsequent phenological curve and the extraction of the phenological index. The specific implementation is shown in formula (8). Rubber_TS(x, y, day Oct.～Nov.,May ) = mean(Rubber_TS(x, y, day Dec.～Apr. )) (8) Rubber_TS(x, y, day) = Rubber_TS(x, y, day) - Rubber_TS(x, y, day Oct.～Nov.,May ) for each day in the non-rubber forest phenological period of the target pixel in a certain year, Rubber_TS(x, y, day Dec.～Apr. ) is the kNDVI value of the target pixel in a certain year Rubber forest phenological period of all daily kNDVI values, using the average of all dates of kNDVI of rubber forest phenological period instead of daily kNDVI value in non-phenological period, so as to smooth the phenological curve.

2. The method according to claim 1, wherein, In step S1, the kNDVI time series is constructed based on the daily data of MOD09GQ and MOD09GA, and cloud removal is performed. First, the sur_refl_b01_1 red band and sur_refl_b02_1 near-infrared band are extracted from the MOD09GQ data, and the daily kernel normalized vegetation index kNDVI of the study area is calculated. The calculation formula is as follows: k(a,b) = exp(-(a-b) 2 / (2σ 2 )) (2) Wherein, R, NIR are sur_refl_b01_1 band and sur_refl_b02_1 band of MOD09GQ image data; k(NIR, NIR), k(NIR, R) in formula (1) are calculated according to formula (2); σ is the length scale parameter, which is taken as the mean value of R and NIR; Then, the kNDVI daily time series obtained by calculation is subjected to cloud removal, and the MOD09GA quality control band is used to remove noise and retain good quality pixels. The main steps are as follows: the 1-kilometer resolution daily quality evaluation flag is extracted from the state_1km_1 band of MOD09GA product, which is resampled to 250m resolution. According to this, the daily reliable pixels of the calculated daily kNDVI data are screened, and the reliable pixels are retained.

3. The method according to claim 1, wherein, In step S4, the four phenological indexes are defined as follows: 1) The beginning of defoliation: it is the turning point of the rubber forest from the stable growth period to the rapid defoliation period, the vegetation index value decreases, indicating that the greenness decreases, and it means that the rubber forest enters the defoliation period; 2) The end of defoliation: it is the point at which the vegetation index value of the rubber forest decreases to the minimum and no longer decreases, indicating that the defoliation period has ended; 3) The beginning of germination: it is the point at which the vegetation index of the rubber forest starts to rise after a short period of stability, that is, the rubber forest enters the leaf regrowth period; 4) The end of leaf expansion / leafing period: it is the point at which the vegetation index value of the rubber forest rises to a certain value and no longer rises, entering a relatively stable stage, indicating that the rubber forest has completed leaf regrowth, and the leaves are in the peak period, and the vegetation index value is high; For the kNDVI daily time series of the target rubber forest pixel after step S3, the asymmetric Gaussian function fitting algorithm is used for curve fitting, and then the seasonal amplitude method is used to extract the four phenological indexes of the rubber forest in each year.

Images