Mountain area rubber forest expansion remote sensing detection method based on phenological characteristics

Abstract

The invention relates to the field of remote sensing image recognition, and discloses a mountainous area rubber forest expansion remote sensing detection method based on phenological characteristics. According to the method, the special geographical environment of the plateau mountain area is comprehensively considered, non-homologous data voting and random stratified sampling methods are combined, various remote sensing image features and the unique leaf falling phenomenon of the mountain area rubber forest are fully utilized, the phenological window is identified, the spatial distribution and expansion conditions of the rubber forest are monitored, time, manpower, material resources and financial resources are saved, and the method is suitable for popularization and application. And the method has mobility and universality in remote sensing recognition research of specific ground features in different regions, and has important significance in constructing regional ecological environment and protecting grain safety.

Description

Technical Field

[0001] The invention relates to the technical field of remote sensing image monitoring, and in particular to a remote sensing detection method for rubber forest expansion in mountainous areas based on phenological characteristics. Background Art

[0002] Traditional rubber plantation monitoring methods are mainly based on manual field surveys, which are costly, long-term, and cannot achieve continuous monitoring. Remote sensing is a science and technology based on the interaction between electromagnetic waves and the earth's surface materials to detect, analyze and study the earth's resources and environment, and reveal the spatial distribution characteristics and spatiotemporal changes of various surface elements. It has the advantages of large-scale synchronous observation, strong timeliness, and little human interference. It has been widely used in land use / cover classification, crop identification, and vegetation phenology monitoring.

[0003] Although the use of remote sensing technology to automatically identify rubber plantations has become a necessary and important means, most studies are still in the exploratory stage. There are still the following difficulties in scientifically and accurately obtaining its spatial distribution information: (1) Rubber is an evergreen broad-leaved vegetation. Its spectral characteristics are easily mixed with other vegetation types such as natural forests, tea gardens, orchards, and shrubs. The supervised or unsupervised classification methods that rely solely on spectral characteristics have great uncertainty. (2) With the development of remote sensing technology, the emergence of medium- and high-resolution remote sensing images has provided a new way to identify rubber plantations. However, most rubber plantations are distributed in tropical or subtropical rainforest areas with mixed vegetation. Frequent cloudy and rainy weather greatly reduces the availability of images. It is difficult to establish long-term time series data only through optical images of a single sensor. (3) Studies have confirmed that in the northern tropical marginal areas south of the Tropic of Cancer, rubber will experience a unique leaf fall phenomenon in the low temperature environment during the dry season. The emergence of this phenological feature provides a new idea for the effective identification of rubber plantations. How to obtain high-precision, long-term distribution information of rubber forests and then judge the expansion of rubber forests to provide a basis for balancing the economic benefits and environmental issues of rubber forests locally is an issue that technical personnel in this field urgently need to solve. Summary of the invention

[0004] In order to solve the above problems, the present invention provides a remote sensing detection method for rubber forest expansion in mountainous areas based on phenological characteristics;

[0005] To implement the above technical solution, the specific steps are as follows:

[0006] S1. Obtain multi-source satellite remote sensing data and ALOS12.5m DEM terrain data of the target area, and perform data preprocessing on the acquired remote sensing data;

[0007] Multi-source satellite remote sensing data are mainly obtained based on Google Earth Engine (GEE), including: multispectral images and radar data;

[0008] Multispectral images include: Landsat 7 data, Landsat 8 data, Sentinel-2 multispectral images, and Gaofen-2 data;

[0009] Radar data uses Sentinel-1 radar images, including vertical transmission / vertical reception VV, horizontal transmission / horizontal reception HH, vertical transmission / horizontal reception VH and horizontal transmission / vertical reception HV, four polarization combinations, all with a resolution of 10m;

[0010] ALOS12.5m DEM is collected by the ALOS satellite (Advanced Land Observing Satellite), with a horizontal / vertical accuracy of 12.5m, and can be used for all-day and all-weather land observation;

[0011] Although radar data is not affected by cloudy and rainy weather, noise has a great impact on data quality. Sentinel-1 radar data is processed by thermal noise removal, radiation correction and speckle filtering, and terrain correction is performed using ALOS12.5mDEM terrain data.

[0012] Preprocessing includes: atmospheric correction, orthorectification, radiation correction, data clipping, mosaicking and cloud removal of multispectral data; the spatial resolution of Landsat 7 and Landsat 8 in the study is 30m, which is automatically sampled to 10m using the GEE platform to match the Sentinel-2 imagery, and the coordinate system is unified through an embedded algorithm so that each pixel can accurately represent the same range on the ground;

[0013] The Sentinel-1 radar data was processed by thermal noise removal, radiation correction, and speckle filtering, and terrain correction was performed using ALOS12.5mDEM terrain data;

[0014] Limited by the low image availability caused by cloudy and foggy weather in plateau mountainous areas, the present invention uses two data sources, multispectral and radar, to combine the advantages of different sensors in temporal and spatial resolution to obtain complete time series data;

[0015] In order to obtain images that can cover the entire study area, the images were screened. The specific steps are as follows:

[0016] 1. Spatial screening: Use the vector boundary of the study area to create a 10km buffer and filter out all images that intersect the boundary;

[0017] 2. Time screening: Taking into account the availability of Landsat 7, Landsat 8, and Sentinel-2 images, the time range was ultimately limited to 2014-2020;

[0018] 3. Attribute screening: Based on the cloud-containing pixels (LOUDY_PIXEL) and the QA60 cloud removal band, the cloud mask function is used to remove the pixels affected by clouds in each image in the study area, and the pixels with qualified quality are used in the subsequent calculation of surface parameters;

[0019] S2. Combine the GF-2 data and the existing three sets of high-resolution (10m) non-homologous land use / cover data to select samples as input samples for the random forest classification algorithm, as follows:

[0020] Specifically, three sets of high-resolution (10m) non-homologous land use / cover data products were used, combined with 1m resolution remote sensing images (obtained by the Gaofen-2 satellite), and sample data were selected by combining stratified sampling, non-homologous data voting and visual interpretation;

[0021] 3 sets of high-resolution (10m) non-homologous land use / cover data products include: ESA_2020_10m, Esri_Land_Cover_2020_10m and Google Dynamic World;

[0022] ESA_2020_10m is a data product produced by ESA in collaboration with a number of scientific research institutions around the world;

[0023] Esri_Land_Cover_2020_10m is a data product produced by Esri, an American company, using deep learning methods;

[0024] Google Dynamic World is a near-real-time 10m resolution global land use / cover dataset produced by Google using GEE and artificial intelligence based on Sentinel-2 images using a deep learning framework;

[0025] The input samples include six categories: natural forest, cultivated land, tea plantation, water body (river, lake, reservoir), impervious surface (construction land, road) and rubber plantation;

[0026] S3. Use the Harmonic Analysis of Time Series (HANTS) method to reconstruct the spectral index, and obtain the key phenological window period of rubber based on the reconstructed spectral index; at the same time, use the median synthetic image of the key phenological window period of rubber as the classification image of the random forest classification algorithm;

[0027] The expression of time series harmonic analysis method is as follows:

[0028]

[0029] In the formula, a0, a j , b jDifferent frequencies f j The sine and cosine waves, m is the number of harmonics, and t represents the corresponding time sequence;

[0030] Spectral indices include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI), Normalized Difference Backscatter Index (NDI_VV), Fractional Vegetation Cover (FVC); the corresponding phenological window period of rubber plantations can be calculated through each spectral index;

[0031] The specific steps of reconstruction are as follows:

[0032] Landsat 7 and Landsat 8 datasets were used to complete the missing Sentinel-2 data in the same time series after cloud removal, and the time series curves of NDVI, EVI and LSWI were reconstructed using the Harmonic Analysis of Time Series (HANTS) method to solve the problem of incomplete regional coverage of perennial time series multispectral data in plateau and mountainous areas.

[0033] The remote sensing images in GEE of the key phenological window period of rubber were synthesized by mean synthesis to obtain the synthetic images of the key phenological window period of rubber, and the synthetic images of the key phenological window period of rubber were used as the input images of the random forest classification algorithm at the same time; in order to avoid the influence of feature redundancy on the classification results, the phenological period was subdivided with the help of GEE platform. Since the leaf fall of rubber in the dry and hot period of the dry season accounts for 73.87% of the annual leaf fall, the brightness (TCT-BRI), greenness (TCT-GRE), and humidity (TCT-WET) components obtained based on the Tasselled Cap Transformation (TCT) of remote sensing images were added to highlight the leaf fall characteristics during this period;

[0034] S4. Perform accuracy evaluation and classification processing on the random forest classification algorithm to verify the effectiveness of the method;

[0035] The accuracy evaluation indicators include: Overall classification accuracy (OA), producer accuracy (PA), user accuracy (UA) and Kappa coefficient; OA and Kappa coefficient are used to evaluate the overall classification accuracy, PA and UA are used to evaluate the misclassification and omission errors between different land cover types, and their expressions are as follows:

[0036]

[0037] (p0=OA, )

[0038] In the formula, i represents the pixel; n represents the total number of pixels; n ii represents the pixels on the diagonal of the confusion matrix, m represents the total number of pixels on the diagonal of the confusion matrix; q represents the number of categories in the confusion matrix; n i· Represents the sum of row pixels of a class in the confusion matrix; n ·i Represents the sum of column pixels of a certain class in the confusion matrix;

[0039] The classification processing is: combine the different feature combinations of normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), land surface water index (LSWI), normalized backscatter index (NDI_VV), vegetation coverage (FVC), brightness (TCT-BRI), greenness (TCT-GRE), and humidity (TCT-WET) to establish 5 classification schemes: Scheme 1: NDVI, EVI; Scheme 2: NDVI, EVI, LSWI; Scheme 3: NDVI, EVI, LSWI, TCT-BRI, TCT-GRE, TCT-WET; Scheme 4: NDVI, EVI, LSWI, TCT-BRI, TCT-GRE, TCT-WET, FVC; Scheme 5: NDVI, EVI, LSWI, TCT-BRI, TCT-GRE, TCT-WET, FVC, NDI_VV;

[0040] Beneficial effects of the present invention

[0041] 1. Aiming at the phenomenon that the terrain in mountainous areas is fragmented, the underlying surface is complex, and rubber forests are easily mixed with natural forests, the present invention proposes to use multi-source satellite remote sensing data to solve the problem of missing remote sensing image data in mountainous areas. By determining the phenological window of rubber forests and performing random forest classification, the spatial distribution of rubber forests in long time series can be quickly extracted to judge the expansion of rubber forests, saving time, manpower, material resources and financial resources.

[0042] 2. The present invention provides a remote sensing detection method for rubber forest expansion in mountainous areas based on phenological characteristics, which has transferability and versatility in remote sensing identification research of specific land types in different regions. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 It is a schematic diagram of the process of the present invention;

[0044] Figure 2 It is a schematic diagram of the framework of the present invention;

[0045] Figure 3It is the spectral curve reconstructed by the time series harmonic analysis method, among which part a is the reconstructed spectral curve of the normalized difference vegetation index (NDVI); part b is the reconstructed spectral curve of the enhanced vegetation index (EVI); part c is the reconstructed spectral curve of the land surface water index (LSWI); part d is the reconstructed spectral curve of the normalized backscatter index (NDI_VV);

[0046] Figure 4 This is the spatial distribution of Xishuangbanna rubber plantations in four phases: 2014, 2016, 2018 and 2020. DETAILED DESCRIPTION

[0047] See also Figure 1-2 In this embodiment, Xishuangbanna Dai Autonomous Prefecture is used as the research area to illustrate the method provided by the present invention;

[0048] Xishuangbanna Dai Autonomous Prefecture is an internationally important tropical rainforest distribution area. The rubber plantation area has increased from 1.3% of the prefecture's area in 1976 to 22.14% in 2014. Rubber plantations have become one of the typical and dominant land use types in the area. The continuous expansion of its planting area has also caused a series of problems such as the intensified fragmentation of natural forest habitats, the weakening of ecosystem carbon sequestration capacity, the loss of biodiversity, drought and water shortage. Therefore, accurate acquisition of spatial distribution information and expansion of rubber plantations is crucial to regional ecological security and sustainable economic and social development.

[0049] A remote sensing detection method for rubber forest expansion in mountainous areas based on phenological characteristics, the specific steps are as follows:

[0050] S1. Obtain multi-source satellite remote sensing data and ALOS12.5m DEM terrain data of the target area, and perform data preprocessing on the acquired terrain data;

[0051] Among them, multi-source satellite remote sensing data are obtained based on Google Earth Engine (GEE);

[0052] As shown in Table 1, the multispectral images include: Landsat 7 data, Landsat 8 data, and Sentinel-2 multispectral images;

[0053] Landsat 7 and Landsat 8 are surface reflectance data products obtained after atmospheric correction, covering visible light (VIS), near infrared (NIR), and shortwave infrared (SWIR) bands, with a resolution of 30m, a revisit period of 16 days, and a swath width of 185km;

[0054] Sentinel-2 is a surface reflectance data product obtained after atmospheric and orthorectification correction. It provides 13 spectral bands, with a resolution of 10m for blue (B2), green (B3), red (B4) and near-infrared (B8) bands; 20m for red end (B5), near-infrared NIR (B6, B7 and B8A) and short-wave infrared SWIR (B11 and B12); 60m for coastal atmospheric aerosol (B1), water vapor band (B9) and cirrus band (B10); the revisit period of one satellite is 10 days, and the revisit period of two satellites is 5 days.

[0055] The radar data uses the Sentinel-1 data product, which is an image dataset that has undergone Doppler centroid estimation, single-view composite focusing, and post-processing (calling the data set that has been pre-processed by the GEE data processing platform). Each image has 4 bands, corresponding to 4 polarization combinations: vertical transmission / vertical reception VV; horizontal transmission / horizontal reception HH; vertical transmission / horizontal reception VH; horizontal transmission / vertical reception HV; the resolution is 10m; although radar data is not affected by cloudy and rainy weather, noise has a greater impact on data quality. Sentinel-1 radar data is processed by thermal noise elimination, radiation correction and speckle filtering, and terrain correction is performed using ALOS12.5m DEM terrain data; ALOS12.5m DEM is collected by the ALOS satellite (Advanced Land Observing Satellite) with a horizontal / vertical accuracy of 12.5m, which can be used for all-day and all-weather land observations;

[0056] Preprocessing includes: atmospheric correction, orthorectification, radiation correction, data clipping, mosaicking and cloud removal for multispectral data; among them, cloud removal includes: thermal noise removal, radiation correction, speckle filtering for Sentinel-1 radar data, and terrain correction using ALOS12.5m DEM terrain data.

[0057] The Sentinel-1 radar data was processed by thermal noise removal, radiation correction, and speckle filtering, and terrain correction was performed using ALOS12.5mDEM terrain data;

[0058] The spatial resolution of Landsat 7 and Landsat 8 in this study is 30m, which is automatically sampled to 10m using the GEE platform to match the Sentinel-2 imagery, while the coordinate system is unified through an embedded algorithm so that each pixel can accurately represent the same range on the ground;

[0059] Table 1 Multi-source satellite remote sensing data

[0060]

[0061] Limited by the low image availability caused by cloudy and foggy weather in plateau mountainous areas, the present invention uses two data sources, multispectral and radar, to combine the advantages of different sensors in temporal and spatial resolution to obtain complete time series data;

[0062] In order to obtain images that can cover the entire study area, the images were screened. The specific steps are as follows:

[0063] 1. Spatial screening: Use the vector boundary of the study area to create a 10km buffer and filter out all images that intersect the boundary;

[0064] 2. Time screening: Taking into account the availability of Landsat 7, Landsat 8, and Sentinel-2 images, the time range was ultimately limited to 2014-2020;

[0065] 3. Attribute screening: Based on the cloud-containing pixels (LOUDY_PIXEL) and the QA60 cloud removal band, the cloud mask function is used to remove the pixels affected by clouds in each image in the study area, and the pixels with qualified quality are used in the subsequent calculation of surface parameters;

[0066] S2. Combine the GF-2 data and the existing three sets of high-resolution (10m) non-homologous land use / cover data to select samples as input samples for the random forest classification algorithm, as follows:

[0067] like Figure 2 As shown, specifically, three sets of high-resolution (10m) non-homologous land use / cover data products were used, combined with 1m resolution remote sensing images (GF-2), and sample data were selected by combining stratified sampling, non-homologous data voting and visual interpretation;

[0068] 3 sets of high-resolution (10m) non-homologous land use / cover data products include: ESA_2020_10m, Esri_Land_Cover_2020_10m and Google Dynamic World;

[0069] ESA_2020_10m is a data product produced by ESA in collaboration with a number of scientific research institutions around the world;

[0070] Esri_Land_Cover_2020_10m is a data product produced by Esri, an American company, using deep learning methods;

[0071] Google Dynamic World is a near-real-time 10m resolution global land use / cover dataset produced by Google using GEE and artificial intelligence based on Sentinel-2 images using a deep learning framework;

[0072] The specific data information of this embodiment is shown in Table 2;

[0073] Table 2 Three sets of high-resolution (10 m) land use / cover data from different sources

[0074]

[0075] Furthermore, the sample data were selected by combining 1m resolution remote sensing imagery (GF-2) with three sets of high-resolution (10m) non-homologous land use / cover data, using a combination of stratified sampling, non-homologous data voting and visual interpretation;

[0076] The input samples are pixels, including six categories: natural forest, cultivated land, tea plantation, water body (river, lake, reservoir), impervious surface (construction land, road) and rubber plantation;

[0077] In the selection process of non-rubber plantation sample points (except tea plantations), pure pixels with consistent classification of natural forest, cultivated land, water body and impervious surface in three sets of high-resolution (10m) non-homologous land use / cover data were selected as the selection range, and sample points were randomly generated within the range;

[0078] Since none of the above three sets of high-resolution (10m) non-homologous land use / cover data subdivided tea gardens, the sample points of tea gardens and rubber gardens were selected point by point using manual visual interpretation in combination with 1m resolution remote sensing images (GF-2) as an aid in the selection process;

[0079] Finally, this example obtained 6569 sample points, including 1051 natural forests, 843 cultivated lands, 1000 tea gardens, 216 water bodies (rivers, lakes, reservoirs), 459 impervious surfaces (construction land, roads), and 3000 rubber plantations ( Figure 3 The sample data of 2014, 2016, 2018 and 2020 are shown); the sample data are randomly divided into 70% and 30% as training samples and validation sample data of the random forest classification algorithm respectively;

[0080] S3. Use the Harmonic Analysis of Time Series (HANTS) method to reconstruct the spectral index to solve the problem of incomplete regional coverage of perennial time-series multispectral data in plateau mountainous areas, and obtain the key phenological window period of rubber based on the reconstructed spectral index; at the same time, use the median synthetic image of the key phenological window period of rubber as the classification image of the random forest classification algorithm;

[0081] like Figure 2As shown in the figure, the time series harmonic analysis is an algorithm based on improved Fourier analysis, which can unify the temporal and spatial variation laws of remote sensing images according to their characteristics in the time and space dimensions. The expression is as follows:

[0082]

[0083] In the formula, a0, a j , b j Different frequencies f j The sine and cosine waves, m is the number of harmonics, and t represents the corresponding time sequence;

[0084] The spectral indexes include: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI), Normalized Backscatter Index (NDI_VV), Fractional Vegetation Cover (FVC);

[0085] Normalized Difference Vegetation Index (NDVI), which can effectively reflect the density and intensity of vegetation growth through calculations between near-infrared bands and red light bands;

[0086] The Enhanced Vegetation Index (EVI) is an optimized vegetation index that increases sensitivity to areas of high vegetation cover and enhances the ability to monitor vegetation canopy changes;

[0087] The Land Surface Water Index (LSWI) can effectively reflect the changes in vegetation canopy, soil moisture and soil surface water content;

[0088] The normalized difference backscatter index (NDI_VV) is very sensitive to the orientation and structural changes of vegetation canopy and the water content of leaves, and can effectively distinguish vegetation, bare land and built-up land;

[0089] Fractional vegetation cover (FVC) refers to the proportion of the ground surface covered by vegetation within a certain range;

[0090] Since the study area is located in the tropical rainforest, the vegetation is mixed, and it is affected by cloud, rain and fog weather, resulting in a low number of available images in the area. In addition, the quality of remote sensing images is difficult to guarantee during the rainy season from May to October each year (the peak season for vegetation growth). Therefore, the Landsat 7, Landsat 8 and Sentinel 2 datasets were reconstructed, and the missing Sentinel 2 data of the same time series after cloud removal were mosaicked using the Landsat 7 and Landsat 8 datasets. The time series curves of the three indexes, NDVI, EVI and LSWI, were reconstructed using the Harmonic Analysis of Time Series (HANTS) method to solve the problem of incomplete regional coverage of perennial time series multispectral data in plateau mountainous areas.

[0091] The NDVI, EVI, LSWI and NDI_VV indices were calculated using four datasets from Landsat 7, Landsat 8, Sentinel 2 and Sentinel 1 radars, and then sorted in time series to form a new dataset. The HANTS algorithm was then used to reconstruct the spectral index curves, as shown in the following figure. Figure 3 As shown;

[0092] The window period is the leaf-falling period and the leaf-full period. The calculation formulas of the reconstructed spectral indices and the corresponding phenological window periods of rubber plantations are shown in Table 3.

[0093] Table 3 Spectral index calculation formula and rubber plantation phenological window period

[0094]

[0095] In the table, ρ nir is the near-infrared band, ρ red is the red light band, ρ blue is the blue light band, ρ swir2 The image is from the Sentinel-2 shortwave infrared band 2. VV is the polarization mode of vertical transmission and reception. VH is the polarization mode of vertical transmission and horizontal reception. NDVI max is the maximum value of NDVI, NDVI max is the minimum NDVI value;

[0096] Natural forest is the most easily confused vegetation type with rubber, but the vegetation index of rubber in the deciduous period and rapid growth period is obviously different from that of natural forest. The multi-source remote sensing images of rubber in the key phenological window period are synthesized by mean in GEE to obtain the synthetic image of rubber in the key phenological window period, and it is also used as the input feature of the random forest classification algorithm to avoid the influence of feature redundancy on the classification result. With the help of GEE platform, the phenological period is subdivided by median synthesis. Since the leaf fall of rubber in the dry and hot period of the dry season accounts for 73.87% of the annual leaf fall, the brightness (TCT-BRI), greenness (TCT-GRE), and humidity (TCT-WET) components obtained based on the Tasselled Cap Transformation (TCT) of remote sensing images are added to highlight the leaf fall characteristics during this period; the transformation coefficients used are shown in Table 4.

[0097] Table 4 Transformation coefficients

[0098]

[0099]

[0100] S4, performing accuracy evaluation and classification processing on the random forest classification algorithm;

[0101] The accuracy evaluation indicators include: Overall classification accuracy (OA), producer accuracy (PA), user accuracy (UA) and Kappa coefficient; OA and Kappa coefficient are used to evaluate the overall classification accuracy, PA and UA are used to evaluate the misclassification and omission errors between different land cover types, and their expressions are as follows:

[0102]

[0103] (p0=OA, )

[0104] In the formula, i represents the pixel; n represents the total number of pixels; n ii represents the pixels on the diagonal of the confusion matrix, m represents the total number of pixels on the diagonal of the confusion matrix; q represents the number of categories in the confusion matrix; n i· Represents the sum of row pixels of a class in the confusion matrix; n ·i Represents the sum of column pixels of a certain class in the confusion matrix;

[0105] The classification processing is: combine the different feature combinations of normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), land surface water index (LSWI), normalized backscatter index (NDI_VV), vegetation coverage (FVC), brightness (TCT-BRI), greenness (TCT-GRE), and humidity (TCT-WET) to establish 5 classification schemes: Scheme 1: NDVI, EVI; Scheme 2: NDVI, EVI, LSWI; Scheme 3: NDVI, EVI, LSWI, TCT-BRI, TCT-GRE, TCT-WET; Scheme 4: NDVI, EVI, LSWI, TCT-BRI, TCT-GRE, TCT-WET, FVC; Scheme 5: NDVI, EVI, LSWI, TCT-BRI, TCT-GRE, TCT-WET, FVC, NDI_VV;

[0106] Different feature combinations were input into the random forest classifier to obtain the classification results, as shown in Table 5;

[0107] Table 5 Classification results

[0108]

[0109] Scheme 5 has the most input features involved in classification, and the final classification accuracy is also the highest. After analyzing the optimal combination results, the user accuracy (UA) of rubber forest is 88.8%, the producer accuracy (PA) is 95.2%, the overall classification accuracy (OA) is 90%, and the Kappa coefficient is 0.86, which meets the production needs.

[0110] The following conclusions can be drawn from the analysis:

[0111] 1. There is a mixed classification phenomenon between rubber plantations and tea plantations, cultivated land and natural forests. The classification accuracy of tea plantations is relatively low, mainly because the spectral and texture characteristics of tea plantations and rubber plantations are very similar. Although the phenological period obtained based on the time series curve can effectively distinguish most tea plantations, they have been in an environment of fragmented terrain and mixed vegetation for a long time, which increases the difficulty and uncertainty of identification;

[0112] 2. Except for tea gardens, the PA and UA values ​​of other land use / land cover types are all above 88%, among which the PA values ​​of water bodies, artificial surfaces, cultivated land and rubber plantations exceed 94%, which can achieve better classification results;

[0113] 3. The overall classification accuracy is 90.0%, the Kappa coefficient is 0.86, and the PA and UA of the rubber plantation can reach more than 85%, meeting the needs of spatial analysis and practical applications;

[0114] Using the feature combination of scheme 5 and inputting it into the random forest classifier, we can obtain the spatial distribution of rubber plantations in Xishuangbanna in 2014, 2016, 2018 and 2020, as shown in the following figure: Figure 4 Rubber plantations are mainly concentrated in Jinghong City and Mengla County in the central and eastern parts of the country, and are relatively less distributed in Menghai County in the west. Three typical rubber plantation planting areas were randomly selected according to different geomorphic features for comparison with Google Earth high-resolution remote sensing images. The comparison results showed that the rubber plantations in these three areas could be accurately identified, and the identified boundary information had a high degree of consistency with the image. This also confirmed the accuracy and effectiveness of the method, and it can determine the expansion of the four-phase rubber forest.

[0115] In summary, the example of the present invention provides a remote sensing detection method for the expansion of rubber forests in mountainous areas based on phenological characteristics. The method combines non-homologous data voting and random stratified sampling, uses the spectral characteristic time series curve of rubber forests, judges the key phenological window, combines the phenological window period with other input characteristics, and monitors the spatial distribution and expansion of rubber forests. Compared with the traditional visual interpretation method, the invention comprehensively considers the special geographical environment characteristics of plateau mountainous areas, makes full use of multiple input characteristics, and combines the unique leaf-falling phenomenon of rubber forests in mountainous areas. It can provide strong support for identifying the expansion of rubber forests and protecting the local environment. The invention can quickly and accurately identify the spatial distribution of rubber forests, saving time, manpower, material resources and financial resources, and the remote sensing classification and identification research of specific objects in different regions has mobility and versatility, which is of great significance for building regional ecological environment and protecting food security.

[0116] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A remote sensing monitoring method for rubber forest expansion in mountainous areas based on phenological characteristics, characterized in that: The following steps are involved: S1. Obtain multi-source satellite remote sensing data and ALOS12.5mDEM terrain data of the target area, and perform data preprocessing on the acquired remote sensing data; The multi-source satellite remote sensing data includes: multispectral images and radar data; wherein the multispectral images include: Landsat 7 data, Landsat 8 data, Sentinel 2 multispectral images, and Gaofen 2 data; the radar data uses Sentinel 1 radar images; S2, combining GF-2 data and existing high-resolution non-homologous land use / cover data to select samples as input samples for the random forest classification algorithm; S3. Use the time series harmonic analysis method to reconstruct the spectral index, and obtain the key phenological window period of rubber based on the reconstructed spectral index; at the same time, use the median synthetic image of the key phenological window period of rubber as the classification image of the random forest classification algorithm; S4. Perform accuracy evaluation and classification processing on the random forest classification algorithm to verify the effectiveness of the method.

2. A remote sensing monitoring method for rubber forest expansion in mountainous areas based on phenological characteristics according to claim 1, characterized in that: The preprocessing includes: performing atmospheric correction, orthorectification, radiation correction, data clipping, mosaicking and cloud removal on multispectral data; performing thermal noise elimination, radiation correction, speckle filtering on Sentinel-1 radar data, and performing terrain correction using ALOS12.5mDEM terrain data.

3. The remote sensing monitoring method for rubber forest expansion in mountainous areas based on phenological characteristics according to claim 1 is characterized in that: The sample selection was carried out in combination with the existing high-resolution non-homologous land use / cover data. Specifically, the sample data were screened by combining stratified random sampling, non-homologous data voting and visual interpretation; The input samples are pixels, including natural forests, cultivated land, tea plantations, water bodies, impervious surfaces and rubber plantations.

4. The remote sensing monitoring method for rubber forest expansion in mountainous areas based on phenological characteristics according to claim 1 is characterized in that: The reconstructed spectral index includes: normalized vegetation index, enhanced vegetation index, surface water index, normalized backscatter index, and vegetation coverage; the expressions are as follows: In the formula, ρ nir is the near-infrared band, ρ red is the red light band, ρ blue is the blue light band, ρ swir2 The image is from the Sentinel-2 shortwave infrared band 2. VV is the polarization mode of vertical transmission and reception. VH is the polarization mode of vertical transmission and horizontal reception. NDVI max is the maximum value of NDVI, NDVI max is the minimum NDVI value.

5. The remote sensing monitoring method for rubber forest expansion in mountainous areas based on phenological characteristics according to claim 1 is characterized in that: The classification process combines the normalized vegetation index, enhanced vegetation index, surface water index, normalized backscatter index, vegetation coverage, brightness, greenness, and humidity; the brightness, greenness, and humidity are obtained through the tassel cap transformation of the synthetic image of the rubber forest phenological window period.

6. The remote sensing monitoring method for rubber forest expansion in mountainous areas based on phenological characteristics according to claim 1 is characterized in that: The accuracy evaluation includes: overall classification accuracy, producer accuracy, user accuracy and Kappa coefficient, and the expressions are as follows: In the formula, OA represents the overall classification accuracy; PA represents the producer accuracy; UA represents the user accuracy; i represents the pixel; n represents the total number of pixels; n ii represents the pixels on the diagonal of the confusion matrix, m represents the total number of pixels on the diagonal of the confusion matrix; q represents the number of categories in the confusion matrix; n i· Represents the sum of row pixels of a class in the confusion matrix; n ·i Represents the sum of the column pixels for a class in the confusion matrix.

Images