Crop fine classification method based on random forest model

Abstract

The invention discloses a crop fine classification method based on a random forest model, and relates to the technical field of agricultural remote sensing. The method comprises the following steps of data acquisition and preprocessing, space mask processing, sample point selection, model training, classification and verification. Through smooth processing of MOD09Q1 NDVI data and time series data and a random forest classification model, the problems of classification efficiency and precision of a traditional method under a complex planting structure and a low-time-resolution scene are solved.

Description

Technical Field

[0001] The present invention relates to the field of agricultural remote sensing technology, and in particular to a crop fine classification method based on a random forest model. Background Art

[0002] High-precision classification of crop planting structures is a core technology in precision agriculture, crucial for agricultural production planning, optimal resource allocation, ecological and environmental monitoring, and food security assessment. With the advancement of remote sensing technology, crop classification methods based on multi-temporal and multi-spectral remote sensing data have become a research hotspot in agricultural remote sensing. Currently, mainstream classification methods include supervised classification based on spectral features, time series analysis based on phenological characteristics, and classification models based on machine learning.

[0003] However, existing methods still face many technical bottlenecks in practical applications, including: spatiotemporal resolution limitations: medium and low-resolution remote sensing data (such as MODIS) find it difficult to capture the detailed features of small-scale farmland, while high-resolution data (such as Sentinel-2 and Landsat) have insufficient temporal coverage, resulting in incomplete extraction of crop phenological information; noise interference: Affected by cloud cover, atmospheric conditions and sensor errors, vegetation index time series data such as NDVI often contain noise, which directly affects the input quality of the classification model; insufficient classification accuracy: in areas with complex planting structures (such as rotation areas and intercropping areas), traditional classification methods have limited ability to distinguish crop types and are prone to inter-class confusion; poor algorithm adaptability: existing methods mostly rely on fixed thresholds or single features, which are difficult to adapt to the classification needs of different regions and different crops, and have weak generalization capabilities. In particular, given the high degree of farmland fragmentation and diverse cropping patterns in my country, the effective integration of multi-source remote sensing data, particularly through data preprocessing to improve the quality of basic data, lay a good underlying data foundation for subsequent classification data models, and establish high-precision and robust classification models, has become a key scientific issue that needs to be addressed in the field of agricultural remote sensing monitoring. Therefore, developing a crop classification method that can accommodate spatial and temporal resolution, suppress data noise, and adapt to complex planting structures has important theoretical and applied value for promoting the development of precision agriculture and is a technical challenge that urgently needs to be addressed. Summary of the Invention

[0004] In order to solve the above problems, the present invention provides a crop fine classification method based on the random forest model. By combining high-resolution remote sensing data, time series data processing technology and machine learning algorithm (random forest), it aims to achieve fine classification of crop planting structure, thereby improving the efficiency and accuracy of crop classification and providing accurate data support for agricultural resource monitoring.

[0005] To achieve the above object, the technical solutions adopted by the present invention are as follows:

[0006] A crop fine classification method based on a random forest model includes the following steps:

[0007] S1 Data Acquisition and Preprocessing: The MOD09Q1 surface reflectance product of the MODIS data series was acquired through the Google Earth Engine platform, and the Normalized Difference Vegetation Index (NDVI) within a preset time range was extracted. An algorithm was used to smooth the acquired NDVI time series data. Specifically targeting the signals of long-term remote sensing data series, this algorithm decomposes periodic signals and filters noise to generate a smooth, continuous NDVI time series that reflects the true phenological characteristics of vegetation. This smoothed NDVI time series dataset serves as the data foundation for subsequent crop classification.

[0008] S2 spatial mask processing: Combined with China's 30-meter annual land cover data, the normalized vegetation index data is spatially masked to extract the spatial data of farmland and forest areas, while removing data interference from non-agricultural land areas;

[0009] S3 Sample Point Selection: Based on Google Earth Pro high-definition historical imagery, select sample points for major crops within the study area. These crops include, but are not limited to, winter wheat-summer corn rotation, spring corn monoculture, cotton, rice, fruit trees, and vegetable planting areas. Sample points must be selected evenly based on crop type and spatial distribution, ensuring a total annual sample size of no less than 700.

[0010] S4 Model Training: A crop classification model was trained using the random forest algorithm on the Google Earth Engine platform, using smoothed Normalized Difference Vegetation Index (NDVI) time series data and selected sample point data as model input features. The sample points were first matched with NDVI data, and the complete NDVI time series features were extracted at the sample point locations to construct a training sample set. A random forest classifier containing multiple decision trees was then used for training and modeling. The data was divided into training and test sets in an 8:2 ratio, and the number of decision trees in the random forest model was set to 500. This training effectively learned the relationship between NDVI variation patterns and crop types, generating the final crop classification model.

[0011] S5 Classification and Verification: Use the trained random forest model to classify the target area and generate a spatial distribution map of crops. Calculate the overall classification accuracy and test coefficient of the model by constructing a confusion matrix to verify the model performance.

[0012] As a preferred technical solution of the present invention, the selection of sample points in step S3 specifically includes:

[0013] When selecting sample points, based on Google Earth Pro high-definition historical images and combined with the dynamic change characteristics of NDVI time series data, major crop sample points are selected each year according to crop phenological periods (including the greening period of winter wheat and the boll-opening period of cotton) and planting structures (including mulching and greenhouses). These sample points cover various typical planting types, including winter wheat and summer corn rotation areas, spring corn, fruit and forestry, cotton, vegetables, and rice. Samples are selected according to crop type and spatial distribution to ensure uniform coverage of sample points, with a total annual sample size of ≥700. For missing years, sample points from neighboring years are used to supplement the image data, constructing a high-quality sample database that is continuous in time and space and balanced in categories to provide reliable input for model training.

[0014] As a preferred technical solution of the present invention, the specific implementation of the algorithm in step S1 includes:

[0015] For the proposed normalized vegetation index data within the preset time range, that is, the NDVI data of all MOD09Q1 images, each pixel corresponds to a time series, and the reference format of the time dimension parameter is set to year × n observation phases per year, and the reference interval is synthesized into days;

[0016] Then, the core is that, based on the significant annual periodicity of vegetation growth, namely the phenological period, the NDVI time series is regarded as a periodic signal with a period of one year, and short-term noise and long-term trends are superimposed. Then, the Fourier series can be introduced for data fitting: the NDVI sequence of each pixel is regarded as a function f(t), where t is time, in days, one scene for 8 days, and a total of 46 scenes in a year; it is decomposed into the sum of harmonic components of different frequencies through the Fourier series, and the following relationship is obtained: NDVI signal = basic value + multiple periodic fluctuations of different amplitudes + random noise; among them, low-order harmonics correspond to the main seasonal changes and correspond to the spring and summer growing seasons and autumn and winter decline periods of crops in the northern hemisphere, and high-order harmonics correspond to high-frequency noise and correspond to short-term weather fluctuations or other interference factors;

[0017] Then, on this basis, the least squares method is further introduced to calculate the amplitude and phase of each harmonic component to minimize the sum of square errors between the fitting curve and the original NDVI series; at the same time, signal reconstruction is performed through noise filtering. On the one hand, the 1st to 5th harmonics are retained by frequency selection, and higher-order harmonics are ignored to effectively filter out noise; then the retained low-frequency harmonic components are used to reversely synthesize a new NDVI time series, filtering out short-term abnormal fluctuations in the original data and retaining the interannual and seasonal scale variation trends; for the terminal data, a typical boundary processing algorithm is constructed to avoid discontinuities at the endpoints, and the edge data are supplemented by extrapolating the fitted harmonic model for the starting and ending years of the time series; in this way, a complete time series curve for each year within the preset years can be generated in the study area, which contains a smoothed NDVI time series dataset for each pixel, and the NDVI value of each phase reflects the actual growth status of vegetation after noise removal.

[0018] As a preferred technical solution of the present invention, in step S1:

[0019] The normalized vegetation index is calculated by subtracting the red light reflectance from the near infrared band reflectance, and then dividing the difference by the sum of the near infrared band reflectance and the red light band reflectance.

[0020] When the Fourier series is introduced for data fitting, a constant term is used to represent the average value of the normalized vegetation index time series data, and each order harmonic component is used for data reconstruction; the harmonic order is determined according to actual needs, the period is the total duration of the time series, and the time is a specific time point in the time series; on this basis, a typical data fitting paradigm is adopted, and the constant term of the Fourier series is divided by 2 and added to the sum of each order harmonic component. Each order harmonic component consists of a cosine term and a sine term. The cosine term is the harmonic coefficient multiplied by the cosine function value, and the sine term is the harmonic coefficient multiplied by the sine function value. The independent variables of the cosine function and the sine function are 2π multiplied by the harmonic order multiplied by time divided by the period.

[0021] As a preferred technical solution of the present invention, the training of the random forest model in step S4 specifically includes:

[0022] The normalized vegetation index data of 46 time series bands were used as input features. When dividing the training set and the test set, random number sequences were added to ensure the balance of the spatial distribution of samples. The number of decision trees was set to 500 to ensure the classification accuracy of the model.

[0023] As a preferred technical solution of the present invention, the model performance verification in step S5 specifically includes:

[0024] The confusion matrix is ​​used to compare the consistency of the classification results with the actual observations; the overall classification accuracy is calculated, that is, the proportion of all correctly classified samples to the total number of samples; the test coefficient is calculated by dividing the difference between the observed consistency and the expected consistency by 1 minus the expected consistency, which is used to evaluate the degree to which the classification result is better than random classification.

[0025] As a preferred technical solution of the present invention, the MOD09Q1 surface reflectance product used in the method has a temporal resolution of 8 days and a spatial resolution of 250 meters.

[0026] As a preferred technical solution of the present invention, the Chinese 30-meter annual land cover data used in the method includes land cover type information of farmland, forest land, water bodies, and construction land.

[0027] As a preferred technical solution of the present invention, when selecting sample points through Google Earth Pro high-definition historical images in the method, it is necessary to combine the phenological characteristics of the crop growing season to ensure the typicality and representativeness of the sample points.

[0028] As a preferred technical solution of the present invention, the training process of the random forest model in the method further includes evaluating the importance of input features and screening feature variables that have a significant impact on the classification results.

[0029] Compared with the prior art, the present invention has the following technical advances:

[0030] The present invention combines the random forest model with a high-precision sample library to solve the classification problem of complex planting areas; by combining the smoothing processing of MOD09Q1 NDVI data time series data with the random forest classification model, the fine classification of crop planting structure is achieved.

[0031] The algorithm eliminates NDVI time series noise and filters out abnormal short-term fluctuations in the original data. It is particularly effective in retaining and highlighting the changing trends on the interannual and seasonal scales (the NDVI curve of wheat after processing will show the typical form of "slow rise before winter → stable during the wintering period → rapid rise in spring → peak during the filling period → decline during maturity"; the curve of corn will show a single-peak feature of "rapid growth after spring sowing → peak in summer → decline before autumn harvest"), greatly improving data quality and laying the data foundation for subsequent improvement of classification accuracy through random forest model training.

[0032] The present invention utilizes complete temporal information to comprehensively reflect crop growth dynamics and improve classification accuracy. Traditional threshold methods may only use a few time points (such as the NDVI peak moment) and are susceptible to interference. Moreover, the model is based on automatic learning from data, reducing the subjectivity of manual threshold setting and having strong adaptability. Traditional methods require manual threshold setting, such as NDVI > 0.4 to identify a certain crop. These thresholds rely on expert experience and are difficult to generalize across regions or years. In addition, vegetation is greatly affected by climate, and the NDVI peak value of the same crop may vary in different years. Selecting sample points for training data sets each year can adapt to such changes. Finally, the present invention can handle multiple crop types and support automated classification processes, with good scalability and engineering application value.

[0033] Overall, we can see that the method we constructed can effectively identify the planting structure of different crops and provide accurate remote sensing data support for agricultural resource management and food security monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0035] In the attached figure:

[0036] Figure 1 The present invention proposes a flow chart of a method for fine classification of crops based on a random forest model.

[0037] Figure 2 This is the NDVI map of the North China Plain in 2013 after Hants smoothing by the present invention.

[0038] Figure 3 This is a graph of the NDVI of a certain pixel in the North China Plain in 2013 after Hants smoothing according to the present invention.

[0039] Figure 4 This is a graph obtained by spatially masking the NDVI time series data using land cover data (CLCD) in the present invention.

[0040] Figure 5 This is a spatial distribution map of crop sample points in 2013 marked by the present invention.

[0041] Figure 6 This is a graph showing the average NDVI corresponding to the crop sample points in 2013 marked by the present invention.

[0042] Figure 7 This is a crop classification spatial distribution map of the North China Plain in 2013, which was made by the present invention according to a crop fine classification method based on a random forest model. DETAILED DESCRIPTION

[0043] The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments. The embodiments of the present invention will be described below with reference to the accompanying drawings.

[0044] Example 1: Taking the classification of crops in the North China Plain in 2013 as an example, the specific steps are as follows:

[0045] like Figure 1 The present invention proposes a flow chart of a method for fine classification of crops based on a random forest model.

[0046] S1 Data acquisition and preprocessing: The MODIS (Moderate Resolution Imaging Spectroradiometer) data series MOD09Q1 surface reflectance product provided by NASA was obtained through the GEE (Google Earth Engine) platform, and the Normalized Difference Vegetation Index (NDVI) time series data for 2013 were obtained (see Figure 2 ), the constructed fitting algorithm is used to smooth the NDVI time series data, remove noise, and obtain a smooth time series curve to ensure data quality and stability (see Figure 3 ).

[0047] The mathematical expression of the normalized difference vegetation index (NDVI) in step S1 is:

[0048] (a)

[0049] Where: NIR is the near infrared band, and R is the red light band.

[0050] The mathematical expression of the fitting algorithm is:

[0051] (b)

[0052] Where: is a constant term that represents the average value or overall trend; and is the selected frequency component used for data reconstruction; is the period of the time series; is the number of harmonics, is the number of time series.

[0053] The core consideration in constructing the fitting algorithm is to use the algorithm to smooth the acquired normalized vegetation index time series data. Specifically targeting the signals of long-term remote sensing data series, through periodic signal decomposition and noise filtering, a smooth, continuous NDVI time series that can reflect the true phenological characteristics of vegetation is generated as a smoothed normalized vegetation index time series dataset, which serves as the data basis for subsequent crop classification.

[0054] For the proposed normalized vegetation index data within the preset time range, that is, the NDVI data of all MOD09Q1 images, each pixel corresponds to a time series, and the reference format of the time dimension parameter is set to year × n observation phases per year, and the reference interval synthesis days are specifically;

[0055] Basic remote sensing reflectance data were obtained using the MODIS MOD09Q1 product (250m spatial resolution, 8-day interval). This product is synthesized every 8 days and provides surface reflectance in the red and near-infrared bands. The data were processed using the following process:

[0056] (1) NDVI calculation: The normalized difference vegetation index (NDVI) was calculated for the data of each observation period of MOD09Q1 (see formula (a)), where NIR is the reflectance of the near-infrared band and RED is the reflectance of the red band.

[0057] (2) Constructing time series: Arrange all MOD09Q1 phases (about 46) in a complete year in chronological order to form the NDVI time series of each pixel.

[0058] (3) Time series smoothing and noise filtering: NDVI time series data are susceptible to short-term abnormal fluctuations due to factors such as clouds and shadows, and therefore require noise processing. A periodic signal decomposition method (fitting algorithm) is used to construct a smoother NDVI series. This smoothed time series data can more accurately reflect the phenological changes within the vegetation annual cycle.

[0059] (4) Time dimension structure: Each pixel corresponds to a time series with the structure of “year × approximately 46 8-day period observations per year”, which is used for subsequent classification model training.

[0060] Then, based on the significant annual periodicity of vegetation growth, namely the phenological period, the NDVI time series is considered as a periodic signal with a period of one year. By superimposing short-term noise and long-term trends, the Fourier series can be introduced for data fitting: the NDVI series of each pixel is regarded as a function f(t), where t is the time in days, with one scene every 8 days, and a total of 46 scenes in a year;

[0061] Through the Fourier series, it can be decomposed into the sum of harmonic components of different frequencies. Specifically, when the Fourier series is introduced for data fitting, the constant term is used to represent the average value of the normalized vegetation index time series data, and the harmonic components of each order are used for data reconstruction; the harmonic order is determined according to actual needs, the period is the total duration of the time series, and the time is the specific time point in the time series; on this basis, the typical data fitting paradigm is adopted, the constant term of the Fourier series is divided by 2 and added to the sum of the harmonic components of each order. The harmonic components of each order are composed of cosine terms and sine terms. The cosine term is the harmonic coefficient multiplied by the cosine function value, and the sine term is the harmonic coefficient multiplied by the sine function value. The independent variables of the cosine function and the sine function are 2π multiplied by the harmonic order multiplied by time divided by the period;

[0062] The relationship is as follows: NDVI signal = base value + multiple periodic fluctuations of varying amplitudes (corresponding to seasonal variations) + random noise. We associate low-order harmonics, including the 1st and 2nd harmonics, with major seasonal variations, corresponding to the spring and summer growing seasons and the autumn and winter decline periods of crops in the Northern Hemisphere. High-order harmonics correspond to high-frequency noise, corresponding to short-term weather fluctuations or other interference factors.

[0063] Then, on this basis, the least square method is further introduced to calculate the amplitude and phase of each harmonic component to minimize the sum of square errors between the fitting curve and the original NDVI sequence;

[0064] At the same time, signal reconstruction is performed through noise filtering. On the one hand, frequency selection is used to retain the 1st to 5th harmonics and ignore higher harmonics, thereby effectively filtering out noise. Then, the retained low-frequency harmonic components are used to reversely synthesize a new NDVI time series, filtering out short-term abnormal fluctuations in the original data and retaining the variation trends on the interannual and seasonal scales.

[0065] For the terminal data, a typical boundary processing algorithm is constructed to avoid discontinuity at the endpoints. For the starting and ending years of the time series, the edge data are supplemented by extrapolating the fitted harmonic model.

[0066] In this way, a complete time series curve for each year within the preset years in the study area can be generated, which contains a smoothed NDVI time series dataset for each pixel. The NDVI value of each phase reflects the true growth status of vegetation after noise removal.

[0067] In actual applications, it can be programmed as an executable program; during application, additional parameters can be set as needed according to actual conditions.

[0068] S2 spatial mask: Combined with the 2013 China 30-meter annual land cover data (CLCD) on the GEE platform, the NDVI time series data were spatially masked to extract farmland and forestland and remove non-agricultural land interference (see Figure 4 );

[0069] S3 sample point selection: Based on Google Earth Pro high-definition historical imagery and the dynamic characteristics of NDVI time series data, we selected sample points for major crops in 2013 in the North China Plain region according to crop phenology (e.g., winter wheat greening period, cotton boll opening period) and cropping structure (e.g., mulching, greenhouses). These sample points cover typical cropping types, including winter wheat-summer maize rotation, spring maize, fruit trees, cotton, vegetables, and rice. Samples were selected based on crop type and spatial distribution to ensure uniform coverage of sample points, with a total annual sample size of ≥700. For missing years, sample points from adjacent years were used to supplement the image data, thereby constructing a high-quality sample database that is continuous in time and space and balanced in categories, providing reliable input for model training (see [1]). Figure 5 ).

[0070] S4 Model Training and Classification: The classification model was trained using the random forest method on the GEE platform. The smoothed 46 NDVI time series data (representing a one-year time series) and sample data were used as input features for training. The sample points were matched with the NDVI data, and the complete NDVI time series features were extracted at the sample point locations. A training sample set was constructed. Subsequently, a random forest classifier containing multiple decision trees was used for training and modeling. 80% of the data was used as the training set and 20% of the data was used as the test set. Random number sequences were added to ensure a balanced spatial distribution of samples. The remote sensing images were matched and annotated with the sample point data. 500 decision trees were used for training to generate a crop classification model to ensure high classification accuracy (see [1]). Figure 6 );

[0071] S5 Classification Result Evaluation: The trained model was used for classification to generate a crop spatial distribution map. The overall classification accuracy and test coefficient were calculated using the confusion matrix. The overall classification accuracy was 92.26% and the test coefficient was 0.89, verifying that the model performed well. The 2013 crop classification map of the North China Plain was obtained, as shown in the figure below. Figure 7 shown.

[0072] Finally, it should be noted that the above descriptions are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art will be able to modify the technical solutions described in the aforementioned embodiments or substitute equivalents for some of the technical features. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the claims of the present invention.

Claims

1. A crop fine classification method based on random forest model, characterized in that: The following steps are involved: S1 Data Acquisition and Preprocessing: The MOD09Q1 surface reflectance product of the MODIS data series was acquired through the Google Earth Engine platform, and the Normalized Difference Vegetation Index (NDVI) within a preset time range was extracted. An algorithm was used to smooth the acquired NDVI time series data. Specifically targeting the signals of long-term remote sensing data series, this algorithm decomposes periodic signals and filters noise to generate a smooth, continuous NDVI time series that reflects the true phenological characteristics of vegetation. This smoothed NDVI time series dataset serves as the data foundation for subsequent crop classification. S2 spatial mask processing: Combined with China's 30-meter annual land cover data, the normalized vegetation index data is spatially masked to extract the spatial data of farmland and forest areas, while removing data interference from non-agricultural land areas; S3 Sample Point Selection: Based on Google Earth Pro high-definition historical imagery, select sample points for major crops within the study area. These crops include, but are not limited to, winter wheat-summer corn rotation, spring corn monoculture, cotton, rice, fruit trees, and vegetable planting areas. Sample points must be selected evenly based on crop type and spatial distribution, ensuring a total annual sample size of no less than 700. S4 model training: The crop classification model was trained using the random forest algorithm on the Google Earth Engine platform, with smoothed normalized vegetation index time series data and selected sample point data as model input features; First, the sample points were matched with NDVI data, and the complete NDVI time series features were extracted at the sample points to construct a training sample set. A random forest classifier containing multiple decision trees was then used for training and modeling. The entire data set was divided into a training set and a test set in an 8:2 ratio. The number of decision trees in the random forest model was set to 500. This model training effectively learned the relationship between NDVI variation patterns and crop types, and generated the final crop classification model. S5 Classification and Verification: Use the trained random forest model to classify the target area and generate a spatial distribution map of crops. Calculate the overall classification accuracy and test coefficient of the model by constructing a confusion matrix to verify the model performance.

2. The method according to claim 1, characterized in that The selection of sample points in step S3 specifically includes: When selecting sample points, based on Google Earth Pro high-definition historical images and combined with the dynamic change characteristics of NDVI time series data, major crop sample points are selected each year according to crop phenological periods (including the greening period of winter wheat and the boll-opening period of cotton) and planting structures (including mulching and greenhouses). These sample points cover various typical planting types, including winter wheat and summer corn rotation areas, spring corn, fruit and forestry, cotton, vegetables, and rice. Samples are selected according to crop type and spatial distribution to ensure uniform coverage of sample points, with a total annual sample size of ≥700. For missing years, sample points from neighboring years are used to supplement the image data, constructing a high-quality sample database that is continuous in time and space and balanced in categories to provide reliable input for model training.

3. The method according to claim 1, characterized in that The specific implementation of the algorithm in step S1 includes: For the proposed normalized difference vegetation index data within the preset time range, that is, the NDVI data of all MOD09Q1 images, each pixel corresponds to a time series, and the reference format of the time dimension parameter is set to year × n observation phases per year, and the reference interval is synthesized into days; Then, the core is that, based on the significant annual periodicity of vegetation growth, namely the phenological period, the NDVI time series is regarded as a periodic signal with a period of one year. By superimposing short-term noise and long-term trends, it is possible to consider introducing the Fourier series for data fitting: the NDVI series of each pixel is regarded as a function f(t), where t is time, in days, one scene for 8 days, and a total of 46 scenes in a year; it is decomposed into the sum of harmonic components of different frequencies through the Fourier series, and the following relationship is obtained: NDVI signal = basic value + multiple periodic fluctuations of different amplitudes + random noise; among them, low-order harmonics correspond to the main seasonal changes and correspond to the spring and summer growing seasons and autumn and winter decline periods of crops in the northern hemisphere, and high-order harmonics correspond to high-frequency noise and correspond to short-term weather fluctuations or other interference factors; Then, on this basis, the least squares method is further introduced to calculate the amplitude and phase of each harmonic component to minimize the sum of square errors between the fitting curve and the original NDVI series; at the same time, signal reconstruction is performed through noise filtering. On the one hand, the 1st to 5th harmonics are retained by frequency selection, and higher-order harmonics are ignored to effectively filter out noise; then the retained low-frequency harmonic components are used to reversely synthesize a new NDVI time series, filtering out short-term abnormal fluctuations in the original data and retaining the interannual and seasonal scale variation trends; for the terminal data, a typical boundary processing algorithm is constructed to avoid discontinuities at the endpoints, and the edge data are supplemented by extrapolating the fitted harmonic model for the starting and ending years of the time series; in this way, a complete time series curve for each year within the preset years can be generated in the study area, which contains a smoothed NDVI time series dataset for each pixel, and the NDVI value of each phase reflects the actual growth status of vegetation after noise removal.

4. The method according to claim 1, wherein In step S1: The normalized vegetation index is calculated by subtracting the red light reflectance from the near infrared band reflectance, and then dividing the difference by the sum of the near infrared band reflectance and the red light band reflectance. When the Fourier series is introduced for data fitting, a constant term is used to represent the average value of the normalized vegetation index time series data, and each order harmonic component is used for data reconstruction; the harmonic order is determined according to actual needs, the period is the total duration of the time series, and the time is a specific time point in the time series; on this basis, a typical data fitting paradigm is adopted, and the constant term of the Fourier series is divided by 2 and added to the sum of each order harmonic component. Each order harmonic component consists of a cosine term and a sine term. The cosine term is the harmonic coefficient multiplied by the cosine function value, and the sine term is the harmonic coefficient multiplied by the sine function value. The independent variables of the cosine function and the sine function are 2π multiplied by the harmonic order multiplied by time divided by the period.

5. The method according to claim 1, wherein The training of the random forest model in step S4 specifically includes: The normalized vegetation index data of 46 time series bands were used as input features. When dividing the training set and the test set, random number sequences were added to ensure the balance of the spatial distribution of samples. The number of decision trees was set to 500 to ensure the classification accuracy of the model.

6. The method according to claim 1, characterized in that The model performance verification in step S5 specifically includes: The confusion matrix is ​​used to compare the consistency of the classification results with the actual observations; the overall classification accuracy is calculated, that is, the proportion of all correctly classified samples to the total number of samples; the test coefficient is calculated by dividing the difference between the observed consistency and the expected consistency by 1 minus the expected consistency, which is used to evaluate the degree to which the classification result is better than random classification.

7. The method according to any one of claims 1 to 6, characterized in that The MOD09Q1 surface reflectance product used in the method has a temporal resolution of 8 days and a spatial resolution of 250 meters.

8. The method according to any one of claims 1 to 6, characterized in that The Chinese 30-meter annual land cover data used in the method contains land cover type information of farmland, forest land, water bodies, and construction land.

9. The method according to any one of claims 1 to 6, characterized in that When selecting sample points using Google Earth Pro high-definition historical images in the method, the selection must be made in combination with the phenological characteristics of the crop growing season to ensure the typicality and representativeness of the sample points.

10. The method according to any one of claims 1 to 6, characterized in that The training process of the random forest model in the method also includes evaluating the importance of input features and screening feature variables that have a significant impact on the classification results.

Images