Rubber garden distribution and establishment time identification algorithm based on full-time-sequence image

Abstract

The invention discloses a rubber garden distribution and establishment time recognition algorithm based on a full-time-sequence image. The algorithm comprises the following steps: step 1, data acquisition and preprocessing; step 2, determining a plantation establishment year based on time sequence analysis and time correction of LandTrendr; 3, key features of the economic life cycle are extracted; 4, rubber is classified and recognized through the random forest; and step 5, precision verification and evaluation are carried out. According to the method, the key features in the whole economic life cycle of the rubber garden are analyzed, so that the spatial distribution of the rubber garden and the establishment year of the rubber garden are accurately identified; the method provided by the invention has relatively high adaptability and generalizability, and shows a remarkable identification advantage under a complex landscape condition.

Description

Technical Field

[0001] The invention relates to the technical field of remote sensing image processing, and in particular to a method for identifying the distribution of rubber plantations and their establishment years based on remote sensing images. Background Art

[0002] Using remote sensing to accurately monitor the spatiotemporal dynamics of rubber plantations is an urgent need to promote the sustainable development of the rubber industry. With the development of space-based observation technology, remote sensing methods are playing an increasingly important role in rubber plantation monitoring. Although remote sensing monitoring of rubber planting has undergone rapid development, the current mainstream algorithms are based on rubber forests in the northern tropical region, and the application scenarios are generally concentrated in provincial and municipal scale areas. There are significant inconsistencies in the distribution of large-scale rubber plantations and change detection applications. These problems limit the comprehensive and objective understanding of the dynamics of rubber planting in the current core planting areas and the impact of rubber planting on the regional ecological environment. Specifically, the existing technologies mainly have the following problems:

[0003] 1. Limitations of monitoring methods

[0004] (1) Existing methods mostly rely on single-phase or short-time series images. Currently, many remote sensing monitoring methods rely on single-phase images or only use short-time series image data. Early studies mainly relied on single or a small number of optical images. Although this approach can provide more intuitive information in some cases, it lacks continuity in the time dimension and cannot fully reveal the dynamic changes of the objects. Short-time series data may also not reflect the impact of long-term factors such as seasonal changes and climate change on the ecological environment and productivity of rubber plantations.

[0005] (2) It is difficult to fully utilize historical remote sensing data. With the development of remote sensing technology, the accumulation of historical remote sensing image data has become increasingly rich. Past image data contains a large amount of valuable information and can provide a basis for long-term trend analysis. However, existing monitoring methods often fail to fully tap the potential of these historical data, resulting in insufficient understanding of historical changes.

[0006] (3) Insufficient exploration of the characteristics of the economic life cycle of rubber plantations. Existing methods are still insufficient in exploring the characteristics of the economic life cycle of rubber plantations. The economic life cycle of rubber plantations includes multiple stages, such as planting, seedling stage, growth and maturity stage, harvest stage, etc., and the characteristics of these stages are significantly different. A single image or short time series data cannot fully reveal the land use changes and vegetation growth patterns in these different stages. In addition, the productivity of rubber plantations is affected by many factors, such as climate change, management mode (such as fertilization, irrigation), pest and disease control, etc., and existing methods are still insufficient in analyzing the interaction of these complex factors. Therefore, accurate monitoring and dynamic analysis of the economic life cycle of rubber plantations still requires more diverse and in-depth remote sensing technology support.

[0007] These limitations make it difficult for existing technologies to comprehensively and accurately monitor and assess the development status and environmental changes of rubber plantations. Therefore, more diversified monitoring methods are needed, combining long-term series, cross-source data, and high-resolution images to better explore the characteristics of various aspects of the rubber plantation life cycle.

[0008] 2. Incomplete feature extraction

[0009] (1) Failure to fully consider the characteristics of different growth stages of rubber plantations. During the growth process of rubber plantations, the characteristics of different growth stages show significant differences. Rubber trees have different growth characteristics and ground feature change patterns at various stages from planting to maturity (such as seedling stage, juvenile stage, and adult stage). Existing remote sensing monitoring methods may not fully consider the differences between these stages, resulting in the neglect of dynamic changes between growth stages when extracting features. For example, the coverage of rubber trees in the seedling stage is low and the vegetation signal is weak, while the mature stage has stronger canopy reflectance characteristics. Failure to consider the different characteristics of these growth stages may lead to deviations in monitoring results and fail to accurately reflect the status of rubber plantations at different growth stages.

[0010] (2) Insufficient recognition of surface disturbance features in the early stages of plantation construction. In the early stages of rubber plantation construction, large-scale land leveling and vegetation removal operations are usually carried out on the surface. These disturbances will leave obvious features in remote sensing images and affect the subsequent vegetation growth patterns. However, existing feature extraction methods may not be able to fully identify the surface disturbance features in the early stages of plantation construction. For example, exposed soil after land reclamation, traces of vegetation removal, and sparse vegetation in the early stages of reclamation are crucial for subsequent growth monitoring of rubber plantations. If the identification of these initial disturbance features is ignored, it may affect the subsequent growth trend analysis, which in turn leads to insufficient accuracy in judging the construction stage of the rubber plantation.

[0011] (3) Insufficient use of phenological characteristics. The phenological characteristics of plants (such as leaf area index, greenness index, seasonal changes, etc.) are important indicators that reflect the growth status of vegetation and can help analyze the health status of vegetation growth, productivity level, and response to environmental changes. However, existing remote sensing monitoring methods may not be fully utilized when extracting the phenological characteristics of rubber plantations. For example, the growth of rubber trees is affected by multiple factors such as seasonal changes, climatic conditions, fertilization and irrigation. Phenological characteristics can reflect the comprehensive effects of these factors on the growth of rubber plantations. If phenological characteristics are not fully utilized, in-depth analysis of vegetation growth patterns, yield prediction, and health diagnosis may be missed. In addition, the use of phenological characteristics usually relies on the support of long-term remote sensing data. By analyzing the seasonal changes of vegetation, key turning points in the growth cycle (such as leaf fall and leaf extraction) can be revealed. If the monitoring method lacks sufficient time series data or ignores the changes in different phenological stages, it may fail to accurately capture important information in the growth cycle of the rubber plantation, affecting the assessment of the growth health status.

[0012] 3. Uncertainty of application effect

[0013] (1) Poor accuracy under complex terrain conditions. Complex terrain (such as hills and mountains) will affect the acquisition and analysis of remote sensing images, especially when the terrain is undulating. The ground object signals in the remote sensing images are easily affected by shadows, occlusions, etc., thereby reducing the monitoring accuracy. For vegetation types such as rubber plantations, especially on mountainous or sloping land, the reflectivity of the tree canopy may have large deviations due to factors such as terrain slope and sun angle. This terrain effect makes land cover classification and feature extraction based on remote sensing more difficult, especially in the process of automated processing, which may result in misclassification and inaccurate analysis results. Therefore, under complex terrain conditions, the accuracy of remote sensing monitoring often fails to meet expectations, resulting in limitations in the accurate assessment of rubber plantations.

[0014] (2) Limited processing capabilities for mixed pixels. Pixels in remote sensing images usually represent a certain area of ​​the surface, but due to resolution limitations, many pixels may contain multiple types of land objects (for example, rubber plantations and surrounding farmland, woodlands, etc. are mixed in the same pixel). This phenomenon is called "mixed pixels". Currently, many remote sensing analysis methods have limited processing capabilities for mixed pixels, resulting in low accuracy in land object classification and feature extraction in this case. There may be other crops or natural vegetation around the rubber plantation, and these mixed signals are difficult to simply classify into a single category, especially in lower resolution images. If effective unmixing cannot be achieved (for example, through high-resolution images, band weighting methods, or machine learning techniques), it will lead to errors and uncertainties in the monitoring results.

[0015] (3) There are large errors in the estimation of the establishment year. When using remote sensing images to monitor rubber plantations, accurate estimation of the establishment year is very important because it helps analyze the growth stage, historical development and production potential of the rubber plantation. However, due to the limitations of the spatial and temporal resolution of remote sensing images, there are often large errors in estimating the establishment year. The establishment process of the rubber plantation may not be clearly visible in the image, especially in the early stage of plantation establishment, when there may be only a small amount of vegetation or land disturbance, resulting in the remote sensing data being unable to clearly capture the time point of plantation establishment. In addition, the initial growth of rubber trees is usually slow and may not be accurately identified by conventional features such as vegetation index in remote sensing images. Therefore, the estimation of the establishment year often relies on other indirect methods (such as ground surveys, historical records, etc.), which may lead to large errors due to lack of or inconsistent data. This error not only affects the accuracy of the plantation establishment time, but may also affect subsequent growth trends and yield forecasts. Summary of the invention

[0016] In view of this, the present invention provides a rubber plantation distribution and establishment time recognition algorithm based on full temporal images.

[0017] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0018] The rubber plantation distribution and establishment time recognition algorithm based on full temporal images includes the following steps:

[0019] Step 1: Data acquisition and preprocessing:

[0020] Data acquisition includes: long-term optical imaging data, high-resolution optical data, SAR data, and auxiliary data;

[0021] Data preprocessing includes: optical data preprocessing, SAR data preprocessing, and vegetation index calculation;

[0022] Step 2: Determine the establishment time of the plantation based on LandTrendr’s time series analysis and time correction: First, use the LandTrendr algorithm to analyze forest disturbance, then perform time correction based on rubber growth characteristics, establish year determination criteria, and determine the year of plantation establishment;

[0023] Step 3: Extraction of key features of the economic life cycle: including juvenile characteristics, mature characteristics, spectrum, phenology and structural characteristics of the mapping year;

[0024] Step 4: Random forest classification and recognition: Create a classification data feature set and perform classification;

[0025] Step 5: Conduct accuracy verification and evaluation.

[0026] Preferably, in step 1, the long-time optical image data includes: all Landsat 4 / 5 / 7 / 8 / 9 Level 2 Collection 2 Tier 1 data provided by the United States Geological Survey;

[0027] High-resolution optical data include: monthly and semi-annual mosaics of Planet images with a spatial resolution of 4.77 meters and four spectral bands: blue, red, green, and near-infrared;

[0028] SAR data includes horizontal transmission horizontal reception (HH) and horizontal transmission vertical reception (HV) polarization bands, stored as 16-bit digital values;

[0029] Auxiliary data include DEM elevation data, land cover data, and administrative boundary data.

[0030] Preferably, in step 1, the optical data preprocessing includes radiation correction and geometric correction;

[0031] Radiometric correction is: both Landsat and Sentinel-2 convert the digital values ​​of the sensor into radiance through calibration coefficients, and then perform atmospheric correction to obtain ground reflectance;

[0032] Geometric correction is: image registration and projection correction are performed on both Landsat and Sentinel-2 to ensure that the geographic coordinates of the remote sensing images are consistent with the real ground.

[0033] Preferably, in step 1, the SAR data preprocessing includes orthorectification and slope correction;

[0034] The orthorectification is to correct the geometric distortion in the SAR image by using a digital elevation model or a digital surface model to ensure that each pixel is consistent with the actual ground position; convert the SAR image from its original sensor coordinate system to the required map coordinate system to adjust the spatial position;

[0035] The slope correction is performed using a theoretical model according to the slope, azimuth and angle of incidence of the object; the influence of the slope is removed by adjusting the reflection intensity of the image.

[0036] Preferably, in step 1, the vegetation index calculation refers to the calculation of the normalized difference vegetation index, surface water index and normalized burn index in land use change monitoring, and the calculation formulas are as shown in (1), (2) and (3):

[0037]

[0038] In the formula, ρ red , ρ NIR , ρ SWIR1, and ρ SWIR2 They respectively represent the reflectivity of the red band, near infrared band, shortwave infrared band 1 and shortwave infrared band 2 of the remote sensing image.

[0039] Preferably, in step 2, the LandTrendr algorithm parameter optimization value setting range is: maxSegments: 8-15, spikeThreshold: 0.3-0.7, vertexCountOvershoot: 1-2, recoveryThreshold: 0.3-0.7, minObservationsNeeded: 3-5; preferably, the parameter optimization values ​​are set to: maxSegments=11, spikeThreshold=0.5, vertexCountOvershoot=1, recoveryThreshold=0.5, minObservationsNeeded=4.

[0040] Preferably, in step 2, the establishment year of the rubber plantation is determined as the year that satisfies the following four criteria in reverse order:

[0041] Below the minimum NBR threshold is Ts NBR_MIN ;

[0042] NBR START_IRP is the lowest NBR within the three-year buffer, with the starting vertex of the gain fitted by the LandTrendr algorithm as reference;

[0043] Gain amplitude NBR MAG_IRP Above the minimum amplitude threshold Ts ΔNBR ;

[0044] The duration N between the start and end points DUR_IRP More than two years.

[0045] Preferably, in step 3, the juvenile stage characteristics refer to the period from the establishment year of the rubber plantation to the mature year, in the immature stage or immature rubber plantation, the selection characteristics NBR START_IRP 、NBR MAG_IRP and N DUR_IRP To classify;

[0046] The starting year of the maturity characteristic is calculated by adding the plantation establishment year to the number of years normally required for maturity, N. DUR_IRP The LSWI for each maturity year is generated using imagery acquired between day 1 and day 100 of the year. MIN LSWI STD The frequency F at which LSWI is less than 0.15 LSWI＜0.15, the seasons corresponding to rubber leaf senescence, leaf fall and leaf emergence are called SDF seasons; the average values ​​of these characteristics in the mature stage are calculated and recorded as LSWI MIN_AVG_SDF_MRP LSWI STD_AVG_SDF_MRP and F LSWI＜0.15_AVG_SDF_MRP ;

[0047] The mapping year characteristics include optical characteristics and SAR characteristics.

[0048] Preferably, in step 3, the optical characteristics include annual statistical characteristics, growth period characteristics, leaf senescence, leaf fall and leaf extraction period characteristics, and leaf bud development period characteristics;

[0049] Among them, the statistical characteristics of the annual period are the annual median values ​​of the SWIR1, NIR and red bands of Sentinel-2 and NDVI, which are marked as S2_swir1 and S2_swir2 respectively. MED_ANN 、S2_nir MED_ANN 、S2_red MED_ANN and S2_NDVI MED_ANN ; For Planet images, calculate the median of the blue, green, red, NIR bands and NDVI, denoted as Pl_blue MED_ANN Pl_green MED_ANN Pl_red MED_ANN Pl_nir MED_ANN and Pl_NDVI MED_ANN ;

[0050] The growing season is characterized by calculating the average NDVI, denoted as NDVI AVG_GRN ; The average values ​​of NDVI and LSWI between 0 and 10 quantiles are recorded as NDVI IM010_GRN and LSWI IM010_GRN ; and the frequency of NDVI values ​​less than 0.6 in LS2, denoted as F NDVI＜0 . 6_GRN ;

[0051] The characteristics of leaf senescence, leaf fall and leaf emergence are calculated by using the images of the survey year and the previous year to calculate the mean, minimum and standard deviation of LSWI, which are recorded as LSWI AVG_SDF LSWI MIN_SDF , and LSWI MIN_SDF ;

[0052] The characteristics of leaf bud development period were calculated by calculating the median values ​​of blue, green, red, re1, re2, NIR, and SWIR1 bands of Sentinel-2;

[0053] The SAR characteristics are calculated as the annual median values ​​of the HH and HV polarizations of PALSAR-2 during the survey year.

[0054] Preferably, in step 4, the specific method is: summarizing all the characteristics of the immature rubber plantation and mature rubber plantation stages, as well as the characteristics of the annual period, leaf green period, leaf senescence, leaf fall and leaf emergence period, and leaf budding and development period of the mapping year, and the altitude and slope of the DEM data to create a comprehensive feature set.

[0055] Compared with the prior art, the present invention has achieved the following technical effects:

[0056] The present invention realizes accurate identification of the spatial distribution of rubber plantations and their establishment years; the method proposed by the present invention has strong adaptability and scalability, and shows significant identification advantages under complex landscape conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 It is a schematic diagram of the rubber plantation distribution and establishment time recognition algorithm based on full time series images of the present invention;

[0058] Figure 2 It is a LandTrendr fitting sequence diagram in the rubber plantation distribution and establishment time recognition algorithm based on full time series images of the present invention;

[0059] a) Ideal LandTrendr fit; b) Corrected LandTrendr fit; c) Searching for the starting point in the backward direction; d) Searching for the starting point in the forward direction; e) Extra long gain that requires moving the starting point forward; f) Extreme fit with only one gain;

[0060] Figure 3 This is a scatter plot of the Vietnamese rubber distribution and establishment time accuracy assessment based on the rubber plantation distribution and establishment time recognition algorithm of the present invention based on the full time series image. DETAILED DESCRIPTION

[0061] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0062] The present invention discloses a rubber plantation distribution and establishment time recognition algorithm based on full temporal images, comprising the following steps:

[0063] Step 1: Data acquisition and preprocessing:

[0064] Data acquisition includes: long-term optical imaging data, high-resolution optical data, SAR data, and auxiliary data;

[0065] Data preprocessing includes: optical data preprocessing, SAR data preprocessing, and vegetation index calculation;

[0066] Long-term optical image data include: all Landsat 4 / 5 / 7 / 8 / 9 Level 2 Collection 2 Tier 1 data provided by the United States Geological Survey;

[0067] High-resolution optical data include: monthly and semi-annual mosaic images with a spatial resolution of 4.77 meters and four spectral bands: blue, red, green, and near-infrared;

[0068] SAR data includes horizontal transmission horizontal reception (HH) and horizontal transmission vertical reception (HV) polarization bands, stored as 16-bit digital values;

[0069] Auxiliary data include DEM elevation data, land cover data, and administrative boundary data;

[0070] Optical data preprocessing includes radiation correction and geometric correction;

[0071] Radiometric correction is: both Landsat and Sentinel-2 convert the digital values ​​of the sensor into radiance through calibration coefficients, and then perform atmospheric correction to obtain ground reflectance;

[0072] Geometric correction is: image registration and projection correction are performed on both Landsat and Sentinel-2 to ensure that the geographic coordinates of the remote sensing images are consistent with the real ground;

[0073] SAR data preprocessing includes orthorectification and slope correction;

[0074] The orthorectification is to correct the geometric distortion in the SAR image by using a digital elevation model or a digital surface model to ensure that each pixel is consistent with the actual ground position; convert the SAR image from its original sensor coordinate system to the required map coordinate system to adjust the spatial position;

[0075] The slope correction is performed using a theoretical model according to the slope, azimuth and angle of incidence of the object; the influence of the slope is removed by adjusting the reflection intensity of the image;

[0076] Vegetation index calculation refers to the calculation of normalized difference vegetation index, surface water index and normalized burn index in land use change monitoring. The calculation formulas are shown in (1), (2) and (3):

[0077]

[0078] In the formula, ρ red , ρ NIR , ρSWIR1 , and ρ SWIR2 Respectively represent the reflectance of the red band, near infrared band, shortwave infrared band 1 and shortwave infrared band 2 of the remote sensing image;

[0079] Step 2: Time series analysis based on LandTrendr: First optimize the parameters of the LandTrendr algorithm, and then establish the year determination criteria;

[0080] Among them, the parameter optimization value setting range of the LandTrendr algorithm is: maxSegments: 8-15, spikeThreshold: 0.3-0.7, vertexCountOvershoot: 1-2, recoveryThreshold: 0.3-0.7, minObservationsNeeded: 3-5; preferably, the parameter optimization value is set to: maxSegments=11, spikeThreshold=0.5, vertexCountOvershoot=1, recoveryThreshold=0.5, minObservationsNeeded=4.

[0081] The establishment year of a rubber plantation was determined as the year that met the following four criteria in reverse order:

[0082] Below the minimum NBR threshold is Ts NBR_MIN ;

[0083] NBR START_IRP is the lowest NBR within the three-year buffer, with the starting vertex of the gain fitted by the LandTrendr algorithm as reference;

[0084] Gain amplitude NBR MAG_IRP Above the minimum amplitude threshold Ts ΔNBR ;

[0085] The duration N between the start and end points DUR_IRP More than two years;

[0086] Step 3: Extraction of economic life cycle characteristics: including juvenile characteristics, mature characteristics, spectrum, phenology and structure characteristics of the mapping year;

[0087] Among them, the juvenile characteristics refer to the period from the establishment year of the rubber plantation to the mature year. In the immature stage or immature rubber plantation, the selection characteristics NBR START_IRP 、NBR MAG_IRP and N DUR_IRP To classify;

[0088] The starting year of the maturity characteristic is calculated by adding the plantation establishment year to the number of years normally required for maturity, N. DUR_IRP The LSWI for each maturity year is generated using imagery acquired between day 1 and day 100 of the year. MIN LSWI STD The frequency F at which LSWI is less than 0.15 LSWI＜0.15 , the seasons corresponding to rubber leaf senescence, leaf fall and leaf shedding are called SDF seasons; the average values ​​of these characteristics in the mature stage are calculated and recorded as LSWI MIN_AVG_SDF_MRP LSWI STD_AVG_SDF_MRP and F LSWI＜0.15_AVG_SDF_MRP ;

[0089] The mapping year characteristics include optical characteristics and SAR characteristics;

[0090] Optical characteristics include annual statistical characteristics, growth period characteristics, leaf senescence, leaf fall and leaf emergence characteristics, and leaf bud development characteristics;

[0091] Among them, the statistical characteristics of the annual period are the annual median values ​​of the SWIR1, NIR and red bands of Sentinel-2 and NDVI, which are marked as S2_swir1 and S2_swir2 respectively. MED_ANN 、S2_nir MED_ANN 、S2_red MED_ANN and S2_NDVI MED_ANN ; For Planet images, calculate the median of the blue, green, red, NIR bands and NDVI, denoted as Pl_blue MED_ANN Pl_green MED_ANN Pl_red MED_ANN Pl_nir MED_ANN and Pl_NDVI MED_ANN ;

[0092] The growing season is characterized by calculating the average NDVI, denoted as NDVI AVG_GRN ; The average values ​​of NDVI and LSWI between 0 and 10 quantiles are recorded as NDVI IM010_GRN and LSWI IM010_GRN ; and the frequency of NDVI values ​​less than 0.6 in LS2, denoted as F NDVI＜0.6_GRN ;

[0093] The characteristics of leaf senescence, leaf fall and leaf emergence are calculated by using the images of the survey year and the previous year to calculate the mean, minimum and standard deviation of LSWI, which are recorded as LSWI AVG_SDF LSWI MIN_SDF , and LSWI MIN_SDF ;

[0094] The characteristics of leaf bud development period were calculated by calculating the median values ​​of blue, green, red, re1, re2, NIR, and SWIR1 bands of Sentinel-2;

[0095] The SAR characteristics are calculated as the annual median values ​​of the HH and HV polarizations of PALSAR-2 in the mapping year;

[0096] Step 4: Random forest classification and recognition:

[0097] All characteristics of the immature and mature rubber plantation stages, as well as the characteristics of the annual period, leaf green period, leaf senescence, leaf fall and leaf emergence period, and leaf budding and development period in the mapping year, and the elevation and slope of the DEM data were summarized to create a comprehensive feature set;

[0098] Step 5: Accuracy verification and evaluation.

[0099] Embodiment 1:

[0100] like Figure 2 As shown in Figure 1, there are six typical LandTrendr fitting sequences during the identification process. Ideally, the starting point should closely correspond to the starting vertex of the last gain fitted by the LandTrendr algorithm ( Figure 2 a) In the second case, the LandTrendr algorithm identifies several consecutive short gains, each of which lasts less than six years. These consecutive short gains should be combined to determine the short-term fluctuations that meet Ts. ΔNBR The appropriate endpoint of the standard ( Figure 2 b). The starting year is determined by the minimum NBR value found around the most recent gain start vertex within the three-year buffer ( Figure 2 b). You can go backwards ( Figure 2 c) and forward ( Figure 2 d) Search in both directions for the minimum NBR value near the gain start vertex. Another case involves an extremely long gain that meets all criteria, but the establishment year is determined too early because the immature stage of the rubber plantation is too long ( Figure 2 e) The immature stage usually lasts about 7 years, but can be as long as 10 years in marginal areas. Considering that the LT fitting period may include more points, gains exceeding 12 years are defined as super long gains. In this case, the first one of the latest long gains that meets the NBR is determined. START_IRP <Ts NBR_MIN The reverse point is then searched forward with a buffer of 3 years for the minimum NBR value. If this minimum NBR value is used as the starting point and all criteria are met, the establishment year is updated to the year of the searched minimum NBR ( Figure 2 e). An extreme case occurs when there is only one gain in the entire series ( Figure 2 f). Take the first one that is less than Ts NBR_MINThe minimum NBR value within the 3-year buffer is determined using the reverse point of the 20-year buffer as a reference. This point will be used as the starting point, and a forward search will be performed with a 12-year buffer (threshold for excessively long gains) to find the highest NBR as the end point. If the two points meet all five criteria, the starting point year is set to the year of plantation establishment. If no gain meets all four criteria (special case handling), the plantation establishment year is set to 1986 - the year when most areas began to have multiple LandsatTM images. Ts NBR_MIN and Ts ΔNBR The values ​​of are determined by analyzing the distribution of ground samples and are set to rounded 520 and 150, respectively.

[0101] Threshold determination during the recognition process: NBRmin can be adjusted within the range of 400-600, and Magmin can be adjusted within the range of 100-200.

[0102] Example 2: Spatial distribution accuracy verification and evaluation

[0103] The rubber map generated by this method was evaluated using a confusion matrix consisting of 40% randomly selected ground samples separated from the original dataset to ensure representative and unbiased verification; the classification achieved high accuracy in all indicators, with a producer accuracy (PA) of 96.92%, a user accuracy (UA) of 89.87%, and an F1 score of 0.93; the overall accuracy (OA) reached 93.75%, with a kappa coefficient of 0.87, indicating a high degree of consistency with the ground truth data.

[0104] The non-rubber classification performed better, with 97.35% for UA and 89.87% for rubber. The errors for rubber and non-rubber were 10.13% and 2.65% respectively, while the omission errors were 3.08% for rubber and 8.80% for non-rubber.

[0105] Significantly outperforms existing methods (e.g., 75.00% accuracy in WangR_2021).

[0106] Example 3: Establishing Year Accuracy Verification and Evaluation

[0107] The algorithm was evaluated on the establishment year data using 100 random samples generated on the classified rubber distribution map using the “ee.Image.sample” method in the GEE platform. The observations were taken from historical images of Google Earth to visually interpret the establishment year of the rubber plantations. Figure 3 As shown in the figure, the observed age of rubber plantations was compared with the estimated age of establishment, and a linear fit was calculated to obtain the coefficient of determination (R 2 ): 0.99, root mean square error: 0.25 years.

[0108] Obviously, it is highly consistent with the actual statistical data.

[0109] The above description is only a preferred embodiment of the present invention and does not limit the technical scope of the present invention. Therefore, any slight modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. The rubber plantation distribution and establishment time recognition algorithm based on full-time image is characterized by: The following steps are involved: Step 1: Data acquisition and preprocessing: Data acquisition includes: long-term optical imaging data, high-resolution optical data, SAR data, and auxiliary data; Data preprocessing includes: optical data preprocessing, SAR data preprocessing, and vegetation index calculation; Step 2: Determine the year of plantation establishment based on LandTrendr’s time series analysis and time correction: First, use the LandTrendr algorithm to analyze forest disturbance, then perform time correction based on rubber growth characteristics, establish year determination criteria, and determine the year of plantation establishment; Step 3: Extraction of key characteristics of the economic life cycle: including juvenile characteristics, mature characteristics, phenological, spectral and structural characteristics of the mapping year; Step 4: Random forest classification to identify rubber plantations: Create a classification data feature set and perform classification; Step 5: Conduct accuracy verification and evaluation.

2. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1 is characterized in that: In the step 1, Long-term optical image data include: all Landsat 4 / 5 / 7 / 8 / 9 Level 2 Collection 2 Tier 1 data provided by the United States Geological Survey; High-resolution optical data include: monthly and semi-annual composite Planet images with a spatial resolution of 4.77 meters and four spectral bands: blue, red, green, and near-infrared; SAR data includes PALSAR-2 horizontal transmit horizontal receive (HH) and horizontal transmit vertical receive (HV) polarization bands, stored as 16-bit digital values; Auxiliary data include DEM elevation data, land cover data, and administrative boundary data.

3. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In the step 1, the optical data preprocessing includes radiation correction and geometric correction; Radiometric correction is: both Landsat and Sentinel-2 convert the digital values ​​of the sensor into radiance through calibration coefficients, and then perform atmospheric correction to obtain ground reflectance; Geometric correction is: image registration and projection correction are performed on both Landsat and Sentinel-2 to ensure that the geographic coordinates of the remote sensing images are consistent with the real ground.

4. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In the step 1, the SAR data preprocessing includes orthorectification and slope correction; The orthorectification is to correct the geometric distortion in the SAR image by using a digital elevation model or a digital surface model to ensure that each pixel is consistent with the actual ground position; convert the SAR image from its original sensor coordinate system to the required map coordinate system to adjust the spatial position; The slope correction is performed using a theoretical model according to the slope, azimuth and angle of incidence of the object; the influence of the slope is removed by adjusting the reflection intensity of the image.

5. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In step 1, the vegetation index calculation refers to the calculation of the normalized difference vegetation index, surface water index and normalized burn index in land use change monitoring, and the calculation formulas are shown in (1), (2) and (3): In the formula, ρ red , ρ NIR , ρ SWIR1 , and ρ SWIR2 They respectively represent the reflectivity of the red band, near infrared band, shortwave infrared band 1 and shortwave infrared band 2 of the remote sensing image.

6. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In step 2, the parameter optimization value setting range of the LandTrendr algorithm is: maxSegments: 8-15, spikeThreshold: 0.3-0.7, vertexCountOvershoot: 1-2, recoveryThreshold: 0.3-0.7, minObservationsNeeded: 3-5; preferably, the parameter optimization value is set to: maxSegments=11, spikeThreshold=0.5, vertexCountOvershoot=1, recoveryThreshold=0.5, minObservationsNeeded=4.

7. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In step 2, the establishment year of the rubber plantation is determined as the year that meets the following four criteria in reverse order: Below the minimum NBR threshold is Ts NBR_MIN ; NBR START_IRP is the lowest NBR within the three-year buffer, with the starting vertex of the gain fitted by the LandTrendr algorithm as reference; Gain amplitude NBR MAG_IRP Above the minimum amplitude threshold Ts ΔNBR ; The duration N between the start and end points DUR_IRP More than two years.

8. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1 is characterized in that: In step 3, the juvenile stage characteristics refer to the period from the establishment year of the rubber plantation to the mature year, in the immature stage or immature rubber plantation, the selection characteristics NBR START_IRP 、NBR MAG_IRP and N DUR_IRP To classify; The starting year of the maturity characteristic is calculated by adding the plantation establishment year to the number of years normally required for maturity, N. DUR_IRP The LSWI for each maturity year is generated using imagery acquired between day 1 and day 100 of the year. MIN LSWI STD The frequency F at which LSWI is less than 0.15 LSWI＜0.15 , the seasons corresponding to rubber leaf senescence, leaf fall and leaf emergence are called SDF seasons; the average values ​​of these characteristics in the mature stage are calculated and recorded as LSWI MIN_AVG_SDF_MRP LSWI STD_AVG_SDF_MRP and F LSWI＜0.15_AVG_SDF_MRP ; The mapping year characteristics include optical characteristics and SAR characteristics.

9. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In step 3, the optical characteristics include annual statistical characteristics, growth period characteristics, leaf senescence, leaf fall and leaf extraction period characteristics, and leaf bud development period characteristics; Among them, the statistical characteristics of the annual period are the annual median values ​​of the SWIR1, NIR and red bands of Sentinel-2 and NDVI, which are marked as S2_swir1 and S2_swir2 respectively. MED_ANN 、S2_nir MED_ANN 、S2_red MED_ANN and S2_NDVI MED_ANN ; For Planet images, calculate the median of the blue, green, red, NIR bands and NDVI, denoted as Pl_blue MED_ANN Pl_green MED_ANN Pl_red MED_ANN Pl_nir MED_ANN and Pl_NDVI MED_ANN ; The growing season is characterized by calculating the average NDVI, denoted as NDVI AVG_GRN ; The average values ​​of NDVI and LSWI between 0 and 10 quantiles are recorded as NDVI IM010_GRN and LSWI IM010_GRN ; and the frequency of NDVI values ​​less than 0.6 in LS2, denoted as F NDVI＜0.6_GRN ; The characteristics of leaf senescence, leaf fall and leaf emergence are calculated by using the images of the survey year and the previous year to calculate the mean, minimum and standard deviation of LSWI, which are recorded as LSWI AVG_SDF LSWI MIN_SDF , and LSWI MIN_SDF ; The characteristics of leaf bud development period were calculated by calculating the median values ​​of blue, green, red, re1, re2, NIR, and SWIR1 bands of Sentinel-2; The SAR characteristics are calculated as the annual median values ​​of the HH and HV polarizations of PALSAR-2 during the survey year.

10. The rubber plantation distribution and establishment time recognition algorithm based on full temporal imagery according to claim 1, characterized in that: In step 4, the specific method is as follows: all the characteristics of the immature rubber plantation and mature rubber plantation stages, as well as the characteristics of the annual period, leaf green period, leaf senescence, leaf fall and leaf emergence period, and leaf budding and development period of the mapping year, and the altitude and slope of the DEM data are summarized to create a comprehensive feature set.

Images