Remote sensing-statistics-phenology coupled crop sample automatic generation method and device

Abstract

The application discloses a remote sensing-statistics-phenology coupling crop sample automatic generation method and equipment, and adopts harmonic regression to reconstruct a time curve when processing high-dimensional time characteristics of a pixel, effectively removes the influence of cloud and mist noise and data loss; and combines feature adaptive propagation clustering of similarity, compared with an existing K-Means algorithm, without presetting a cluster center number, can better capture crop subclass characteristics under a complex planting structure, and improves the robustness and discriminability of a feature space; in a label assignment link of a clustering cluster, the application introduces a phenology consistency constraint and county-level statistical data, and uses a greedy algorithm and local search optimization to obtain an optimal matching scheme of the clustering cluster and a crop category; finally, label optimization of the pixel is refined, a neighborhood weighted constraint and spectral iterative self-cleaning are adopted, and heterogeneous blocks and misclassified pixels are automatically removed, so that the purity of the pixel is significantly improved.

Description

Technical Field

[0001] This invention belongs to the field of agricultural remote sensing and information technology, specifically relating to an automatic method for generating crop samples by coupling remote sensing, statistics, and phenology, and also relates to computer equipment. Background Technology

[0002] Current crop sample acquisition mainly relies on manual field surveys, which are costly, time-consuming, and have limited sample coverage, making it difficult to meet the needs of large-scale and multi-type crop monitoring and classification. Existing remote sensing-based automatic sample generation methods mostly rely on single remote sensing features or simple clustering, without fully integrating phenological patterns and statistical prior knowledge, resulting in low accuracy and poor adaptability of sample labels, especially in agricultural areas with frequent replanting or fragmented plots.

[0003] Currently, manual field surveys (GPS point marking) are not only time-consuming, labor-intensive, and costly, but also difficult to cover large areas and have long update cycles. While some methods for automatically generating remote sensing samples exist, they suffer from the following significant problems: early studies often relied on a single data source (optical or radar only), resulting in weak identification capabilities when faced with phenomena such as different spectra of the same object or different objects having the same spectra; existing unsupervised clustering methods (such as K-Means, GMM) often generate pseudo-labeled samples containing mixed pixels or incorrect labels, lacking effective cleaning mechanisms; and most methods do not fully utilize government-released agricultural statistics (such as county-level planting area) as prior knowledge to calibrate remote sensing classification results, leading to significant deviations between the final extracted area and official data.

[0004] How to automatically generate high-purity and spatiotemporally representative training samples of various types of crops using multi-source remote sensing data in the absence of field survey samples, and ensure that the generated results are consistent with macroscopic statistical data, has become a common problem in existing solutions. Summary of the Invention

[0005] The purpose of this invention is to address the aforementioned problems in the existing technology by providing an automatic method for generating crop samples using a remote sensing-statistics-phenology coupling approach, and also to provide computer equipment.

[0006] The above-mentioned objectives of the present invention are achieved by the following technical means:

[0007] An automatic crop sample generation method based on remote sensing-statistics-phenology coupling includes the following steps:

[0008] Step 1: Calculate multiple vegetation indices based on Sentinel-2 time-series image data of the study area. Construct high-dimensional time-series feature sets for each pixel using the multispectral reflectance and multiple vegetation indices, and perform smoothing and quality screening processes to obtain pixel-level high-quality time-series feature sets for each pixel. Finally, downscale to the plot level and calculate plot-level high-quality time-series feature sets for each plot.

[0009] Step 2: Perform a weighted average of multiple vegetation indices in the high-quality temporal feature set of each plot to obtain the comprehensive phenological curve of each plot. Based on the comprehensive phenological curve and actual phenological conditions of each plot, divide the clustering time window of each plot.

[0010] Step 3: Using the high-quality temporal feature set of each land parcel located within the clustering time window of the land parcel as input, the feature-adaptive clustering algorithm for land parcels is used to iteratively optimize multiple cluster centers, and all land parcels are assigned to each cluster center to obtain unlabeled clusters centered at each cluster center.

[0011] Step 4: Obtain candidate cluster sets for each crop category based on phenological consistency constraints; then, based on the sown area of ​​each crop category in the county-level statistical data, use a greedy algorithm to assign initial crop category labels to each unlabeled cluster to obtain an initial matching scheme; finally, iteratively optimize the initial matching scheme through local search optimization to obtain the optimal matching scheme.

[0012] Step 5: Pass the crop category labels of the clusters in the optimal matching scheme to the effective pixels to construct the initial set of labeled pixels for each crop category;

[0013] The initial labeled cell set for each crop category is preliminarily screened using a neighborhood weighting constraint, and weights are assigned to the effective cells in the labeled cell set.

[0014] A spectral iterative self-cleaning method is used to construct class centers for each crop category based on the feature vectors and weights of effective pixels in the labeled pixel set of each crop category. Then, the labeled pixel set of each crop category is optimized and filtered through the class centers. Finally, the high-purity pixel set of each crop category is obtained by iteratively updating the labeled pixel set and class centers of each crop category.

[0015] The feature vectors of the effective pixels are high-quality temporal feature sets at the pixel level.

[0016] As mentioned above, the comprehensive phenological curve is obtained in the following way: multiple vegetation indices in the high-quality temporal characteristic set of the plot are weighted and fused into a fused vegetation index to obtain the comprehensive phenological curve of the plot. The comprehensive phenological curve includes the fused vegetation index at each observation time point.

[0017] Each vegetation index is weighted by its normalized Gini importance score.

[0018] The clustering time window for the land parcels is obtained in the following way:

[0019] The peak position of the comprehensive phenological curve of the plot is used as the peak time of the growing season of the corresponding plot, and the number of peaks is identified as the number of times the plot is replanted.

[0020] Step 2.2.1: Pre-divide the growing season based on the crop growth history of the study area;

[0021] Step 2.2.2: Calculate the target threshold based on the pre-divided growing season and comprehensive phenological curve. The target threshold is calculated based on the following formula:

[0022] ;

[0023] In the formula, The target threshold; This represents the minimum value of the integrated vegetation index of the comprehensive phenological curve within the pre-divided growing season period; This represents the maximum value of the vegetation index within the pre-divided growing season period of the composite phenological curve. The set scaling factor;

[0024] Step 2.2.3, when there is a continuous When the fused vegetation index at each observation time point is greater than or equal to the target threshold, then... The first day of each observation point marks the start of the growing season. When continuous When the fused vegetation index at each observation time point is less than or equal to the target threshold, then... The first day of each observation point marks the end of the growing season. Then the clustering time window of the land parcels .

[0025] As described above, obtaining the candidate cluster set for each crop category based on phenological consistency constraints specifically includes the following steps:

[0026] Step 4.1: Select a set of candidate clusters for each crop type, where the first... Candidate cluster set of crop-like organisms Specifically, it is obtained through the following steps:

[0027] Step 4.1.1: Compare the key phenological parameters of each cluster with the first... The key phenological parameters of similar crops are compared. When the comparison result is less than a set time tolerance, the cluster is determined to satisfy the phenological consistency constraint. By traversing all clusters, all clusters that satisfy the first phenological constraint are obtained. Clusters constrained by phenological consistency of similar crops;

[0028] The key phenological parameters include the start time of the growing season, the peak time of the growing season, and the end time of the growing season. The start time, peak time, and end time of the growing season for each cluster are the median of the corresponding parameters for all plots within the cluster.

[0029] Step 4.1.2, calculate the first... Initial prior feature centers of crop-like plants : In satisfying the first In the clusters constrained by phenological consistency of similar crops, calculate the key phenological parameters of each cluster and the first... The sum of the deviations between the standard key phenological parameters of similar crops, and the selection of the crop with the smallest sum of deviations. The first cluster is used as the initial reference cluster; then, the clusters in the initial reference cluster are... The eigenvectors of the i-th cluster are averaged to obtain the i-th cluster. Initial prior feature centers of crop-like plants ;

[0030] The feature vector of the cluster is obtained by taking the median of the high-quality features at each observation time point of all plots within the cluster to form the feature vector of the cluster.

[0031] Step 4.1.3, add all those that satisfy the first... The eigenvectors of clusters with phenological consistency constraints for crop-like crops and the first Initial prior feature centers of crop-like plants Calculate similarity, based on the number of... Initial prior feature centers of crop-like plants similarity From largest to smallest, for all conditions satisfying the first... The clusters of crops with phenological consistency constraints were sorted, and the top clusters were selected. The clusters constitute a candidate cluster set. ;

[0032] Step 4.1.4: Repeat step 4.1.3 to select candidate cluster sets for each type of crop, and then construct the candidate category set corresponding to each cluster by the crop category corresponding to the candidate cluster set to which each cluster belongs.

[0033] As described above, the initial crop category labels are assigned to each unlabeled cluster using a greedy algorithm, specifically including the following steps:

[0034] Step 4.2.1: Sort the clusters by area from largest to smallest, and initialize the remote sensing estimated area for each crop category. =0, where, the first The remote sensing estimated area of ​​crop-like crops is ,initialization , Pick , The total number of crop categories;

[0035] Step 4.2.2: Perform initial crop category label assignment for clusters: Each time initial crop category label assignment is performed for a cluster, the current cluster is sequentially assigned from the candidate category set, and the overall loss of the current assignment is calculated for each assignment. Finally, the total loss relative to the initial allocation is taken. The crop category with the largest decrease is the initial crop category label of the current cluster, if relative to the overall loss of the initial assignment. If none of them decrease, then the crop category initially classified is used as the initial crop category label for the current cluster.

[0036] Simultaneously update the remote sensing estimated area for the corresponding crop category. Then proceed with the initial crop category label assignment for the next cluster;

[0037] Total loss Calculated based on the following formula:

[0038] ;

[0039] In the formula, Crop category; The first in county-level statistical data Planting area of ​​similar crops; To take the absolute value;

[0040] Step 4.2.3: Traverse all clusters in descending order of area, repeating step 4.2.2 until all clusters have completed the initial crop category label assignment. All unlabeled clusters and their corresponding initial crop category labels constitute the initial matching scheme.

[0041] As described above, the initial screening of the labeled cell set for each crop category using neighborhood weighting constraints to obtain the labeled cell set for each crop category specifically includes the following steps:

[0042] Step 5.1.1: First, construct a spatial neighborhood window of a preset size centered on each valid cell in the initial labeled cell set. The window is located at the first valid cell in the initial labeled cell set. The spatial neighborhood window of each effective cell is ;

[0043] Then, count the number of valid pixels with the same crop category label as the corresponding valid pixel within the spatial neighborhood window of each valid pixel in the initial labeled pixel set; simultaneously, count the total number of valid pixels within the spatial neighborhood window of each valid pixel in the initial labeled pixel set. Then, the... The proportion of similar effective pixels Calculated based on the following formula:

[0044] ;

[0045] In the formula, For the first The number of valid pixels with the same crop category label as the corresponding valid pixel within the spatial neighborhood window of each valid pixel. For the first The total number of valid pixels within the spatial neighborhood window of each valid pixel;

[0046] Finally, based on the 8-neighborhood connectivity rule, connectivity analysis is performed on other valid pixels within the spatial neighborhood window of each valid pixel in the initial labeled pixel set that have the same crop category label as the valid pixel. This constructs the connectivity of each valid pixel, and the number of other valid pixels within the connectivity of each valid pixel is calculated. Among these, the 8-neighborhood connectivity rule... The number of other valid pixels within the connected component of each valid pixel is ;

[0047] Step 5.1.2: Set the first and second similarity ratio thresholds as follows: and The connectivity threshold is ;

[0048] Will and The effective pixels are determined to be high-confidence effective pixels, and and The valid pixels are determined to be valid pixels with medium confidence. and The valid pixels were determined to be noisy pixels with low confidence.

[0049] In the initial set of labeled pixels, retain the valid pixels with high and medium confidence, and remove the noise pixels with low confidence to obtain the set of labeled pixels.

[0050] The specific method for assigning weights to valid pixels in the labeled pixel set is as follows: assigning a weight of 1 to high-confidence valid pixels in the labeled pixel set, and assigning a weight equal to the proportion of pixels in the same category to medium-confidence valid pixels in the labeled pixel set; the first pixel in the labeled pixel set... The weight of each valid cell is ,but or ,in, For the first The proportion of similar effective pixels.

[0051] As mentioned above, the high-purity pixel sets for each crop category are obtained through the following steps:

[0052] Step 5.2.1: Construct class centers for each crop category based on the feature vectors and weights of effective pixels in the labeled pixel set for each crop category. Crop-like class center Calculated based on the following formula:

[0053] ;

[0054] in, For the first The weight of each effective cell, For the first Feature vectors of 1 valid pixel; Serial number , This represents the number of valid pixels in the tagged pixel set.

[0055] Step 5.2.2: Calculate the distance between each valid pixel in the labeled pixel set of each crop category and the class center of the corresponding crop category, where the distance is... The first tagged cell in the crop category set The effective pixel and the first Distance between the class centers of crop-like plants for: ,in, It is the Euclidean norm;

[0056] Filtering thresholds for labeled cell sets of each crop category for: ,in, and This represents the average and standard deviation of the distances from each valid cell in the labeled cell set to the class center.

[0057] Step 5.2.3: Iteratively update the labeled cell set: delete cells from the labeled cell set that meet the following conditions. The effective pixels are used to iteratively update the set of labeled pixels, thus obtaining the updated set of labeled pixels for the current iteration;

[0058] Step 5.2.4: Repeat steps 5.2.1 to 5.2.3 for the next iteration using the updated labeled cell set; until the updated labeled cell set no longer changes during iteration or the set maximum number of iterations is reached. At that time, the final high-purity pixel set is obtained.

[0059] As mentioned above, local search optimization includes swapping, moving, and fine-tuning;

[0060] The exchange involves selecting two clusters and swapping their initial crop category labels; if the overall loss... If the value decreases, the exchange will be accepted.

[0061] The movement involves replacing the current initial crop category label of a cluster with another crop category label from the corresponding candidate category set; if the overall loss... Decrease the value and allow movement;

[0062] The fine-tuning involves: under the premise of satisfying the consistency constraint of the phenological window, calculating the area error between the sown area and the remotely estimated area of ​​the crop category in the county-level statistical data under the matching scheme of the current round of local search, sorting the area errors from largest to smallest, and calculating the cumulative area error and total error in sequence according to the sorting order. Clusters of crop categories belonging to the category with a cumulative area error greater than a set percentage of the total error are reassigned; these reassigned clusters are then placed in the candidate set. The crop labels are reassigned to other crop category labels within the specified range, and the loss function is calculated. Change the label to select other crop category tags that would minimize losses;

[0063] Each local search operation executes the above three operations sequentially. When several consecutive rounds of local searches fail to further reduce the speed... Or, when the maximum number of iterations is reached, the optimization stops and the optimal matching solution is obtained.

[0064] As described above, step 3 specifically includes the following steps:

[0065] Step 3.1: Divide the growing season into different growth periods, and group the plots according to the number of times the plots are replanted to obtain different plot groups;

[0066] Step 3.2: Construct a similarity matrix between land parcels. Attraction Matrix and the degree of belonging matrix Similarity matrix Attraction matrix Belonging degree matrix ,in, For the plot of land With the plot of land Similarity between them For the plot of land land parcel The attractiveness For the plot of land land parcel degree of belonging Number of plots, serial number and All ;

[0067] Step 3.3: In each land parcel group, determine the cluster centers through iteration, and assign all land parcels to the respective cluster centers to form unlabeled clusters. The center point of each cluster is the corresponding cluster center. This specifically includes the following steps:

[0068] Step 3.3.1: Set the maximum number of iterations. After each iteration, set the bias parameter for the current iteration. The values ​​are directly assigned to the diagonal elements of the similarity matrix in the current iteration. This process yields the similarity matrix for the next iteration. The bias parameter for the next iteration is calculated based on the bias parameter from the current iteration and the median of the similarity matrix for the next iteration. The bias parameter from the first iteration uses the median of the similarity matrix before the first iteration. This process is then applied to the similarity matrix. Attraction Matrix and the degree of belonging matrix Perform iterative optimization; among which, For the plot of land Self-similarity;

[0069] Step 3.3.2: After each iteration, based on the attraction matrix of the current iteration... With the degree of belonging matrix Constructing a decision matrix According to the decision matrix diagonal elements Determine candidate cluster centers; All plots with a value greater than 0 were recorded as candidate cluster centers, among which, For the plot of land Degree of belonging to oneself For the plot of land Attraction to oneself;

[0070] When in succession In each iteration, if the candidate cluster centers remain unchanged, then the continuous... Candidate cluster centers that remain unchanged throughout the iterations are determined as cluster centers, and the iteration stops. If no cluster centers are found by the maximum number of iterations, the attraction matrix at the maximum number of iterations is used. With the degree of belonging matrix Calculate all Plots with a value greater than 0 are cluster centers;

[0071] Step 3.3.3: Assign all land parcels to various cluster centers, forming unlabeled clusters centered at each cluster center. Specifically, this is based on the following formula: Assigned to cluster center:

[0072] ;

[0073] In the formula, In order to make To obtain the maximum value The value of , The index of the cluster center. Select from the numbers of the plots identified as cluster centers. For the plot of land Cluster centers degree of belonging Cluster center land parcel Its attractiveness.

[0074] The smoothing process described above is based on the following procedure:

[0075] First, the number of globally optimal harmonic combinations for the harmonic regression model is obtained by minimizing the root mean square error;

[0076] Secondly, a harmonic regression model with the globally optimal number of harmonic combinations is used to smooth the high-dimensional time series feature set of all pixels, so as to obtain the harmonic regression-based time series feature set of each pixel. The harmonic regression time series feature set of each pixel includes multiple harmonic regression features at each observation time point, and the number of harmonic regression features is equal to the number of features in the high-dimensional time series feature set.

[0077] The quality screening process is based on the following steps:

[0078] Step 1.4.1: Normalize the harmonic regression time series feature set of each pixel;

[0079] Step 1.4.2: Based on Mahalanobis distance, remove the outlier observation time points and corresponding harmonic regression features from the normalized harmonic regression time series feature set of each pixel to obtain the preprocessed harmonic regression time series feature set of each pixel.

[0080] Step 1.4.3: Use the Pearson correlation coefficient to remove harmonic regression features with a correlation higher than the set correlation threshold from the time series features of the preprocessed harmonic regression time series feature set;

[0081] Then, using univariate statistical tests, all significant harmonic regression features with P values ​​less than the set P value threshold are retained to obtain the significant harmonic regression time series feature set for each pixel. The significant harmonic regression time series feature set includes multiple significant harmonic regression features at each observation time point.

[0082] Step 1.4.4: Finally, input the saliency harmonic regression time series feature set of each pixel into the random forest model to calculate the Gini importance, sort them from high to low, and select the top... Each significant harmonic regression feature is a high-quality feature, resulting in a pixel-level high-quality temporal feature set for each pixel. The pixel-level high-quality temporal feature set for each pixel includes features from each observation time point. One high-quality feature;

[0083] The high-quality temporal feature sets at the plot level for each plot are obtained through the following steps:

[0084] Step 1.5.1: Count the number of observations that pass quality control for each pixel. If the proportion of the number of observations that pass quality control to the total number of observations at all time points of the corresponding pixel is not less than the set proportion threshold, then the corresponding pixel is determined to be a valid pixel, and all valid pixels in each plot are obtained.

[0085] Step 1.5.2: For the pixel-level high-quality temporal feature set of all valid pixels within each plot, the median of the high-quality features at each observation time point is used for fusion to obtain the plot-level high-quality temporal feature set for each plot. The plot-level high-quality temporal feature set includes the features at each observation time point. High-quality characteristics at the plot level.

[0086] A computer device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of the remote sensing-statistics-phenology coupled crop sample automatic generation method as described above.

[0087] Compared with the prior art, the present invention has the following advantages:

[0088] (1) This invention introduces agricultural statistical data constraints into the label allocation process of unsupervised clustering; by minimizing the difference between the cluster area and the statistical area, it achieves automatic and accurate matching between unlabeled clusters and real crop categories.

[0089] (2) The present invention uses harmonic regression to reconstruct time series curves, effectively removing the influence of cloud noise and missing data; and combined with similarity adaptive propagation clustering (AP clustering), compared with the existing K-Means algorithm, it does not require a preset number of cluster centers, can better capture the characteristics of crop subclasses under complex planting structures, and improves the robustness and discriminative power of the feature space.

[0090] (3) By introducing a dynamic target threshold phenology extraction window related to phenology and combining it with neighborhood weighted constraint (NWC) and spectral iterative self-cleaning (SIS) strategies, this invention can automatically remove heterogeneous blocks and misclassified pixels, significantly improving the purity of the pixel set and providing high-quality training data for subsequent deep learning models. Attached Figure Description

[0091] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0092] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to embodiments. The embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0093] Example 1:

[0094] like Figure 1 As shown, the automatic generation method for crop samples using remote sensing-statistics-phenology coupling includes the following steps:

[0095] Step 1: Acquire Sentinel-2 time-series image data of the study area and calculate multiple vegetation indices. Construct high-dimensional time-series feature sets for each pixel in the Sentinel-2 time-series image data based on the multispectral reflectance and multiple vegetation indices. Perform smoothing and quality filtering on each pixel to obtain a high-quality pixel-level time-series feature set. Finally, downscale to the plot level and calculate the plot-level high-quality time-series feature set for each plot. This includes the following steps:

[0096] Step 1.1: Obtain high-quality Sentinel-2 time-series image data of the study area. The Sentinel-2 time-series image data includes Sentinel-2 image data of multiple observation time points arranged in chronological order. Using the multispectral reflectance of the Sentinel-2 image data of each observation time point as input features, extract multiple vegetation indices of the Sentinel-2 image data of each observation time point.

[0097] In this embodiment, the high-quality Sentinel-2 image data specifically uses Sentinel-2 Level-2A (surface reflectance) image data with less than 30% cloud cover from the Google Earth Engine (GEE) cloud platform, which has undergone radiometric calibration and atmospheric correction. The Sentinel-2 Level-2A image data at multiple consecutive observation time points are acquired and arranged sequentially in chronological order to obtain Sentinel-2 time-series image data.

[0098] The multispectral reflectance of Sentinel-2 Level-2A image data represents the surface reflectance across multiple bands. In this embodiment, surface reflectance from the red band (Band 4), near-infrared band (Band 8), and shortwave infrared bands (Band 11 and Band 12) of the Sentinel-2 Level-2A image data is selected as input features. The following indicators are extracted: Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI), Green Chlorophyll Vegetation Index (GCVI), Normalized Difference Senescent Vegetation Index (NDSVI), Normalized Multi-band Drought Index (NMDI), and Normalized Difference Yellowness Index (NDY). Multiple vegetation indices, including NDYI.

[0099] Step 1.2: Arrange all observation time points and corresponding high-dimensional feature sets of each pixel in the Sentinel-2 time series image data in chronological order to form a high-dimensional time series feature set for a long time series. This yields the high-dimensional time series feature set of each pixel in the Sentinel-2 time series image data. The high-dimensional feature set of each observation time point includes the multispectral reflectance and various vegetation indices of the corresponding observation time point.

[0100] Step 1.3: To characterize the periodicity of crop growth, a harmonic regression model is used to decompose the high-dimensional time-series feature set of each pixel. First, some pixels within the study area are randomly selected to fit different harmonic combinations. Second, based on minimizing the root mean square error (RMSE), the number of globally optimal harmonic combinations is determined by comparing the accuracy of different fitting results, ensuring the adaptability and reliability of the decomposition model. Finally, based on the globally optimal harmonic combination number obtained, a unified harmonic decomposition process is performed on the high-dimensional time-series feature set of all pixels within the study area. By separating the trend term, periodic term, and noise term in the data, the harmonic regression time-series feature set of each pixel is obtained. The harmonic regression time-series feature set of each pixel includes multiple harmonic regression features at each observation time point. The number of harmonic regression features is equal to the number of features in the high-dimensional time-series feature set, thus preserving the core information of the crop growth periodicity.

[0101] Step 1.4: Perform outlier removal preprocessing on the harmonic regression time series feature set of each pixel, as follows:

[0102] Step 1.4.1: Normalize the harmonic regression time series feature set of each pixel to accelerate the convergence of subsequent calculations.

[0103] Step 1.4.2: Based on Mahalanobis distance (such as the 3-standard deviation criterion), remove outlier observation time points and corresponding harmonic regression features from the normalized harmonic regression time series feature set of each pixel to obtain the preprocessed harmonic regression time series feature set of each pixel, so as to enhance feature robustness.

[0104] Step 1.4.3: Then, Pearson correlation coefficient is used to remove redundant harmonic regression features with correlation higher than 0.9 from the preprocessed harmonic regression time series feature set; then, univariate statistical test (t-test) is used to retain all significant harmonic regression features with P value less than 0.05, to obtain the significant harmonic regression time series feature set of each pixel. The significant harmonic regression time series feature set includes multiple significant harmonic regression features at each observation time point.

[0105] Step 1.4.4: Finally, input the saliency harmonic regression time series feature set of each pixel into the random forest model to calculate the Gini Importance, sort them from high to low, and select the top scores. Each significant harmonic regression feature is used as a high-quality feature to obtain a pixel-level high-quality temporal feature set for each pixel. The pixel-level high-quality temporal feature set for each pixel includes features from each observation time point. A high-quality feature.

[0106] Step 1.5: Downscale to farmland plots, and fuse the pixel-level high-quality temporal feature sets using the statistics of effective pixels within each plot to generate plot-level high-quality temporal feature sets for each plot. This includes the following steps:

[0107] Step 1.5.1: Determine the validity of pixels to obtain all valid pixels within each plot; whereby valid pixels are determined in the following way: count the number of observations that have passed quality control for each pixel. If the proportion of the number of observations that have passed quality control to the total number of observation points of the corresponding pixel is not less than 70%, then the corresponding pixel is determined to be a valid pixel.

[0108] Quality control involves identifying and eliminating observations that are affected by factors such as clouds, cloud shadows, aerosols, sensor noise, and geometric errors, which may prevent them from accurately reflecting the surface conditions.

[0109] Step 1.5.2: For the pixel-level high-quality temporal feature set of all valid pixels within each plot, the median of the high-quality features corresponding to all valid pixels at each observation time point is used for fusion to obtain the plot-level high-quality temporal feature set for each plot. The plot-level high-quality temporal feature set includes the features from each observation time point. High-quality features at the plot level are used to enhance the stability of time-series curves and reduce the impact of outliers and mixed pixels.

[0110] Step 2: Calculate the weighted average of multiple vegetation indices from the high-quality temporal feature sets at the plot level for each plot to obtain the comprehensive phenological curve for each plot. Then, based on the comprehensive phenological curve and actual phenological conditions, define the clustering time window for each plot. This specifically includes the following steps:

[0111] Step 2.1: Obtain the integrated phenological curve: For multiple vegetation indices in the high-quality temporal feature set at the plot level for each plot, weights are assigned according to the Gini importance score corresponding to each vegetation index to construct an integrated phenological curve; where each vegetation index is weighted by its corresponding normalized Gini importance score; based on the weights, the values ​​of each vegetation index at each observation time point are weighted and fused to obtain the fused vegetation index at each observation time point, and the fused vegetation index at all observation time points constitutes the integrated phenological curve; the peak position of the integrated phenological curve is taken as the peak time of the growing season for the plot. The number of peaks is used to determine the number of times a plot of land can be replanted.

[0112] The number of vegetation indices used to construct the comprehensive phenological curve is denoted as . , No. The Gini importance score for each vegetation index is: Serial Number Then the first Weight of each vegetation index Defined as the normalized Gini importance score, i.e., the... Weight of each vegetation index Based on the following formula:

[0113] (1);

[0114] in, Indicates the first The weight of each vegetation index.

[0115] Step 2.2: The present invention divides the growing season in the following way:

[0116] Step 2.2.1 First, the growing season is pre-divided based on the existing crop growth history of the study area: First, the crop growing season within a year is pre-divided according to the local crop planting system, sowing and harvesting time statistics (growing history or agricultural data) to limit the time search range for phenological extraction and clustering.

[0117] Step 2.2.2: Next, the target threshold is calculated using the dynamic threshold formula based on the pre-divided growing season and comprehensive phenological curve, according to the following formula:

[0118] (2);

[0119] In the formula, The target threshold; This represents the minimum value of the integrated vegetation index of the comprehensive phenological curve within the pre-divided growing season period; This represents the maximum value of the integrated phenological curve, which incorporates vegetation indices within the pre-defined growing season period. The scaling factor is set to 0.3 in this embodiment.

[0120] Step 2.2.3: Next, the pre-divided growing season boundaries are corrected and refined using the target threshold to obtain the growing season, in order to reduce cross-seasonal confusion and improve the temporal consistency of the plots. Specifically, this is based on the following judgment rule: when there is continuous... The fused vegetation index at each observation time point is greater than or equal to the target threshold. At that time, then The first day of each observation point marks the start of the growing season. When continuous The fused vegetation index at each observation time point is less than or equal to the target threshold. At that time, then The first day of each observation point marks the end of the growing season. In this embodiment, .

[0121] To generate clustering time windows that match the planting system for each plot, the clustering time windows are defined directly from the start and end of the growing season calculated in step 2.2.3. The clustering time window limits the time range of subsequent clustering analysis. By processing only the data within the growing season, it can effectively eliminate the interference of data from non-growing seasons (such as bare soil and snow cover periods), focusing on the period when crop growth is most vigorous and the characteristics are most obvious, thereby improving the accuracy of clustering. The clustering time window will serve as a constraint on subsequent time-series slicing and clustering.

[0122] Step 3: Construct a feature-adaptive clustering algorithm oriented towards land parcels. Using the feature vectors of each land parcel as input, iteratively optimize to generate multiple cluster centers, and assign all land parcels to these centers, resulting in unlabeled clusters centered at each center. The feature vectors of the land parcels are high-quality temporal feature sets at the land parcel level within the clustering time window. This specifically includes the following steps:

[0123] Step 3.1, Spatiotemporal Traversal and Initialization: The growing season is divided into three growth periods: emergence, vigorous growth, and maturity. The plots are then grouped based on the number of cropping cycles, resulting in different plot groups. Plots with different cropping systems participate in clustering within their respective groups and time structures. The emergence period is... Wang is long-term Maturity period is ,in, The time buffer window is 15 days in this embodiment (i.e., 15 observation time points).

[0124] Step 3.2, Feature Adaptive Propagation Clustering: Using the feature similarity between land parcels as a metric, an improved similarity adaptive propagation clustering method is performed: the bias parameter is initialized with the median of the similarity matrix, and the attraction matrix is ​​iteratively updated by analyzing the distribution relationship of feature sets within the clustering time window. and membership matrix The automatic determination of cluster centers includes the following steps:

[0125] Step 3.2.1, Construction Cluster similarity matrix The feature vector of each land parcel is used as input. The feature vector of each land parcel is one or more high-quality features among the high-quality time-series features of the land parcel level located within the clustering time window.

[0126] For any two plots and Let the eigenvectors be respectively and The similarity is then defined as the negative of the distance, to ensure that the closer the distance, the greater the similarity, based on the following formula:

[0127] (3);

[0128] In the formula, For the plot of land With the plot of land similarity, It is a Euclidean norm; it can be used to normalize or weightedly fuse features from different sources to eliminate dimensional differences.

[0129] Construct a similarity matrix based on the similarity between all plots. ,but , This refers to the number of land parcels.

[0130] Step 3.2.2: Construct the attraction matrix and membership matrix Attraction matrix and membership matrix They are obtained based on the following formulas:

[0131] (4);

[0132] (5);

[0133] In the formula, For the plot of land land parcel Responsibility, attractiveness Indicates land parcel Suitable as a plot of land The degree of cluster centering; For the plot of land attractiveness to oneself For the plot of land land parcel The attractiveness; For the plot of land land parcel Availability Indicates land parcel Select land parcel The suitability of a cluster center; For the plot of land land parcel degree of belonging; For the plot of land With the plot of land Similarity; and All are plot numbers. and Take 1~ , and It also serves as the plot number. Remove current Value and The serial numbers of all other plots besides the value. Remove current The serial numbers of all other plots besides the value; To find the maximum value function, This is a function that takes the minimum value.

[0134] Step 3.3: In each land parcel group, cluster centers are determined iteratively. All land parcels are assigned to these cluster centers to form unlabeled clusters. The center point of each unlabeled cluster is the corresponding cluster center. This step includes the following steps:

[0135] Step 3.3.1: Set the maximum number of iterations for the similarity matrix. Attraction Matrix and the degree of belonging matrix Iterative optimization is performed: To avoid the sensitivity of manually setting bias parameters, an adaptive adjustment mechanism is introduced. During the iteration process, the bias parameters are dynamically adjusted according to the following rules to control the number of clusters:

[0136] (6);

[0137] In the formula, For the first The bias parameter after the next iteration. For the first The bias parameter after the next iteration. Similarity matrix The median similarity. This represents the current iteration number; the above update rule causes the bias parameter to gradually approach the target value with a step size that decreases with each iteration. This reduces the sensitivity to bias parameters and stabilizes the number of clusters.

[0138] After each iteration, the bias parameter of the current iteration is directly assigned to the diagonal elements of the similarity matrix of the current iteration. This yields the similarity matrix for the next iteration, and thus the similarity matrix for the next iteration step. the median of The bias parameter for the next iteration is calculated, and this process is repeated; the bias parameter for the first iteration is the median of the similarity matrix before the iteration.

[0139] exist In clustering algorithms, diagonal elements Represents the plot of land The prior probability (preference) of becoming a cluster center is used to update the bias parameters according to the formula. Changes have occurred; Changes will directly alter attractiveness. The calculation, and thus the change in affiliation. The calculation.

[0140] Step 3.3.2: After each iteration, based on the attraction matrix of the current iteration... With the degree of belonging matrix Constructing a decision matrix According to the decision matrix diagonal elements Determine candidate cluster centers; All plots with a value greater than 0 were recorded as candidate cluster centers, among which, For the plot of land Degree of belonging to oneself For the plot of land Attraction to oneself;

[0141] When in succession In each iteration, if the candidate cluster centers remain unchanged, then the continuous... Candidate cluster centers that remain unchanged throughout the iterations are determined as cluster centers, and the iteration stops. If no cluster centers are found by the maximum number of iterations, the attraction matrix at the maximum number of iterations is used. With the degree of belonging matrix Calculate all Plots with a value greater than 0 are cluster centers;

[0142] The above methods can achieve automatic determination of cluster centers and automatic generation of cluster affiliation relationships without the need for manual designation of centers.

[0143] Step 3.3.3: Finally, all plots are assigned to cluster centers, forming unlabeled clusters centered at each cluster center. Specifically, the plots are assigned using the following formula. Assigned to cluster center:

[0144] (7);

[0145] In the formula, To maximize the function The value of , In order to make To obtain the maximum value The possible values ​​of , where, The index of the cluster center. Select from the numbers of the plots identified as cluster centers. For the plot of land Cluster centers degree of belonging Cluster center land parcel Its attractiveness.

[0146] Step 4: Obtain candidate cluster sets for each crop category based on phenological consistency constraints; then, introduce the sown area of ​​each crop category from county-level statistical data, and use a greedy algorithm to assign initial crop category labels to each unlabeled cluster to obtain an initial matching scheme; finally, use local search optimization to iteratively optimize the initial matching scheme to obtain the optimal matching scheme; specifically, the following steps are included:

[0147] By introducing statistical data on crop sowing area in counties, comparing the difference between the remote sensing area of ​​clusters under phenological windows and the statistical data, and based on the principle of loss function, setting a judgment rule for crop type matching clusters, the crop type and area information of agricultural statistical data are transmitted to unlabeled clusters.

[0148] To avoid combinatorial explosion caused by fully traversing all cluster-crop type matching schemes, this invention adopts an implementable solution strategy of "greedy initialization + local search optimization".

[0149] Step 4.1: Select a set of candidate clusters for each crop type, where the first... Candidate cluster set of crop-like organisms Specifically, it includes the following steps:

[0150] Step 4.1.1: For any cluster The key phenological parameters of all plots within a cluster are aggregated to obtain the key phenological parameters of the cluster. These key phenological parameters include the start time of the growing season, the peak time of the growing season, and the end time of the growing season. The start time, peak time, and end time of the growing season are defined as the median of the corresponding parameters for all plots within the cluster. The start time, peak time, and end time of the growing season are denoted as: , as well as .

[0151] Clustering Key phenological parameters and the first Key phenological parameters in the standard growth calendar of similar crops , as well as A comparison is made, and if the result of the comparison is less than the set time tolerance, a cluster is determined. Satisfy the first Phenological consistency constraints for crop-like organisms, namely:

[0152] (8);

[0153] (9);

[0154] (10);

[0155] in, , as well as The first The start time of the growing season, the peak time of the growing season, and the end time of the growing season in the standard growing season calendar for similar crops; In this embodiment, time tolerance is indicated. A 15-day period is used; the standard growth period of crops is determined by summarizing the sowing period, vigorous growth period, and maturity and harvest period recorded in agricultural data of the study area, local crop planting system data, or statistical yearbooks.

[0156] Finally, all satisfactions are obtained. Clusters constrained by phenological consistency of crops will not proceed to the next screening step if they do not meet all of the above conditions.

[0157] Step 4.1.2: For clusters that satisfy the phenological consistency constraint, the feature vector of the cluster is further constructed: the feature vector of the cluster is obtained by aggregating the high-quality temporal feature sets of all plots within the cluster; specifically, the median of each plot-level high-quality feature at each observation time point is taken to form the cluster. The cluster-level feature vector, where the cluster is a subset of the cluster. The eigenvectors are denoted as .

[0158] To construct the first The prior feature centers of crop-like organisms first satisfy the first... In the clusters constrained by phenological consistency of similar crops, calculate the key phenological parameters of each cluster and the first... The sum of deviations among key phenological parameters in the standard growth calendar of similar crops is calculated, and the crop with the smallest sum of deviations is selected. Using one cluster as the initial reference cluster, this embodiment... Take 3; then, for the initial reference cluster... The eigenvectors of the i-th cluster are averaged to obtain the i-th cluster. Initial prior feature centers of crop-like plants .

[0159] Step 4.1.3, after obtaining the first Initial prior feature centers of crop-like plants Then, calculate the condition that satisfies the first... The various clusters of crops with phenological consistency constraints and the first Initial prior feature centers of crop-like plants The similarity.

[0160] Among them, clusters With the Similarity of similar crops Defined as:

[0161] (11);

[0162] in, This represents the Euclidean norm.

[0163] According to the Initial prior feature centers of crop-like plants similarity From largest to smallest, for all conditions satisfying the first... The clusters of crops with phenological consistency constraints were sorted, and the top clusters were selected. The clusters constitute a candidate cluster set. ;in, As a preset constant, this embodiment Take 5.

[0164] Step 4.1.4: Repeat step 4.1.3 to select candidate cluster sets for each type of crop, and then construct the candidate category set corresponding to each cluster by the crop category corresponding to the candidate cluster set to which each cluster belongs.

[0165] Step 4.2: Obtain the initial matching scheme using a greedy algorithm, which includes the following steps:

[0166] Step 4.2.1: First, sort the clusters by area from largest to smallest; initialize the remote sensing estimated area of ​​each crop category to 0, where the... The remote sensing estimated area of ​​crop-like plants is denoted as This is called initialization. , Pick , This represents the total number of crop categories.

[0167] Step 4.2.2: Initial category assignment for clusters: Each time an initial category assignment is performed on a cluster, the current cluster is assigned a category from the candidate category set (which it belongs to). The allocation is performed sequentially within the categories, and the total loss for each allocation is calculated. Finally, the total loss relative to the initial allocation is taken. The crop category with the largest decrease is used as the initial crop category label for the current cluster; if relative to the overall loss of the initial assignment... If none of the crop categories decrease, the crop category initially classified is used as the initial crop category label for the current cluster; simultaneously, the remote sensing estimated area of ​​the corresponding crop category is updated. Then proceed with the initial category assignment for the next cluster.

[0168] Among them, the overall loss Calculated based on the following formula:

[0169] (12);

[0170] In the formula, The first in county-level statistical data Planting area of ​​similar crops To take the absolute value.

[0171] Calculate the total loss When calculating the loss function, only the remote sensing estimated area corresponding to the clusters that have been assigned is considered; unassigned clusters are not included in the current loss calculation.

[0172] Step 4.2.3: Traverse all clusters in descending order of area, repeating step 4.2.2 until all clusters have completed the initial category assignment and the initial matching scheme is obtained.

[0173] Step 4.3, Local Search Optimization and Optimal Match Determination, specifically:

[0174] Based on the initial matching scheme, several rounds of local search are performed to further reduce the area difference loss. The local search optimization includes swapping, moving, and fine-tuning:

[0175] Swap: Select two clusters and swap the initial crop category labels; if the overall loss is... If the value decreases, an exchange will be accepted.

[0176] Move: Replace a cluster with another crop category label from the corresponding candidate category set; if the overall loss is... If the size is reduced, the movement will be accepted.

[0177] Fine-tuning: Under the premise of satisfying the consistency constraint of the phenological window, calculate the area error between the sown area and the remotely estimated area of ​​the crop category in the county-level statistical data under the matching scheme of the current round of local search, sort the area errors from largest to smallest, and calculate the cumulative area error and total error in the order of sorting. Reassign the clusters of crop categories whose cumulative area error is greater than 70% of the total error; then, add the reassigned clusters to the candidate set. The crop labels are reassigned to other crop category labels within the specified range, and the loss function is calculated. The change involves selecting labels for other crop categories that would result in the greatest reduction in losses, in order to quickly reduce the area variation within that category.

[0178] Each local search operation executes the above three operations sequentially. When several consecutive rounds of local searches fail to further reduce the speed... Or, when the maximum number of iterations is reached, the optimization stops and the optimal matching solution is obtained.

[0179] Step 5: Transfer the crop labels of the clusters in the optimal matching scheme to the effective pixels to construct the initial labeled pixel set for each crop category; use neighborhood weighting constraints to perform preliminary screening of the initial labeled pixel set for each crop category to obtain the labeled pixel set for each crop category, and assign weights to the effective pixels in the labeled pixel set for each crop category; use spectral iterative self-cleaning to construct the class center for each crop category based on the feature vectors and weights of the effective pixels in the labeled pixel set for each crop category, and then optimize and screen the labeled pixel set for each crop category through the class center; finally, obtain the high-purity pixel set for each crop category by iteratively updating the labeled pixel set and class center for each crop category; specifically including the following steps:

[0180] Step 5.1: After completing the cluster-crop category matching in Step 4, firstly, based on the crop category labels of each cluster in the optimal matching scheme, the crop category labels corresponding to each cluster are passed to each plot within the corresponding cluster, and further passed to the effective pixels within each plot, constructing an initial labeled pixel set by constructing all effective pixels and their corresponding crop category labels.

[0181] Neighborhood weighted constraint (NWC) is used to characterize the consistency constraint between a pixel and its spatial neighborhood. Subsequently, the initial labeled pixel set for each crop category is preliminarily screened using the neighborhood weighted constraint to obtain the labeled pixel set for each crop category. Weights are then assigned to the effective pixels in the labeled pixel set for each crop category. The specific steps include:

[0182] Step 5.1.1, Neighborhood Construction: Construct a fixed-size spatial neighborhood window centered on each effective pixel. In this embodiment, the neighborhood window size is 3×3. Let the size of the neighborhood window be denoted as... The spatial neighborhood window of each effective cell is .

[0183] Similarity ratio calculation: Count the number of valid pixels within the neighborhood window that have the same crop category label as valid pixels. At the same time, the total number of valid cells within the neighborhood window is counted. Then the ratio of similar items is defined as:

[0184] (13);

[0185] In the formula, For the first The proportion of similar effective pixels, For the first The number of valid pixels with the same crop category label as the corresponding valid pixel within the spatial neighborhood window of each valid pixel. For the first The total number of valid pixels within the spatial neighborhood window of each valid pixel.

[0186] Spatial connectivity calculation: Within the neighborhood window, based on the 8-neighborhood connectivity rule, connectivity analysis is performed on other valid pixels within the neighborhood window that match the crop category label of the valid pixel. This constructs the connected components of the valid pixel, and the number of other valid pixels within the connected components of the valid pixel is calculated. The number of other valid pixels within the connected components of the valid pixel is then calculated. The number of other valid pixels within the connected component of each valid pixel is .

[0187] Step 5.1.2: Filter the valid pixels in the initial labeled pixel set through consistency judgment: Let the first and second similarity ratio thresholds be respectively... and The connectivity threshold is .

[0188] Will and The effective pixels are determined to be high-confidence effective pixels, and and The effective pixels are determined as the medium confidence effective pixels; and The valid pixels were identified as noisy pixels with low confidence.

[0189] In the initial set of labeled pixels, retain the valid pixels with high and medium confidence, and remove the noise pixels with low confidence to obtain the set of labeled pixels.

[0190] In this embodiment and Take 0.6 and 0.4, Take 3.

[0191] Step 5.1.3: Assign a weight of 1 to the high-confidence valid pixels in the labeled pixel set, and assign a weight equal to the proportion of pixels in the same category to the valid pixels in the labeled pixel set; denote the first pixel in the labeled pixel set as the first pixel in the labeled pixel set. The weight of each valid cell is ,but or ,in, For the first The proportion of similar effective pixels.

[0192] Neighborhood weighting constraints can suppress interference caused by mixed pixels at plot boundaries, scattered patches, and local misclassification.

[0193] Step 5.2: Spectral Iterative Self-Cleaning (SIS) is used to perform spectral consistency iterative filtering on the retained effective pixels. SIS constructs class centers for each crop category based on the feature vectors and weight coefficients of the effective pixels in the labeled pixel set of each crop category. The pixel sets of each crop category are then optimized and filtered using these class centers. By iteratively updating the labeled pixel set and class centers of each crop category, a high-purity pixel set for each crop category is obtained. This specifically includes the following steps:

[0194] Step 5.2.1: Construct the class center for each crop category using the feature vectors and weights of the effective pixels in the pixel set of each crop category, and define the first... Crop-like class center For the first The weighted average of the feature vectors and weights of all valid pixels in the labeled pixel set of crop-like classes is then calculated. Crop-like class center Calculated based on the following formula:

[0195] (14);

[0196] in, For the first The weight of each effective cell, For the first The feature vectors of each effective pixel are derived from a high-quality temporal feature set at the pixel level; (Sequence number) , This represents the number of valid pixels in the tagged pixel set.

[0197] Step 5.2.2: Calculate the distance between each valid cell in the labeled cell set of each crop category and the class center of the corresponding crop category, and construct a filtering threshold to filter the valid cells in the labeled cell set.

[0198] Among them, the The first tagged cell in the crop category set The effective pixel and the first Distance between the class centers of crop-like plants for:

[0199] (15);

[0200] In the formula, It is the Euclidean norm.

[0201] Filtering threshold Constructed based on the following formula:

[0202] (16);

[0203] in, and This represents the average and standard deviation of the distances from each valid pixel in the labeled pixel set to the class center.

[0204] Step 5.2.3: Iteratively update the labeled cell set: delete cells from the labeled cell set that meet the following conditions. The effective pixels are used to iteratively update the set of labeled pixels, thus obtaining the updated set of labeled pixels for the current iteration.

[0205] Step 5.2.4: Repeat steps 5.2.2 and 5.2.3 for the next iteration with the updated labeled cell set; until the updated labeled cell set no longer changes during iteration or the set maximum number of iterations is reached. In this embodiment Taking 5, we obtain the final high-purity pixel set.

[0206] By employing county-level statistical area constraints, greedy initialization and local search optimization of cluster-category mapping, and joint purification of NWC and SIS, this invention can automatically generate high-purity and time-consistent multi-type crop sample sets without human intervention, providing reliable sample support for subsequent crop identification model training and mapping applications.

[0207] In one embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above method embodiments.

[0208] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon that, when executed by a processor, implements the steps in the above method embodiments.

[0209] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0210] It should be noted that the embodiments described in this invention are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains can make various modifications or additions to the described embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. An automatic method for generating crop samples using a remote sensing-statistics-phenology coupling approach, characterized in that, Includes the following steps: Step 1: Calculate multiple vegetation indices based on Sentinel-2 time-series image data of the study area. Construct high-dimensional time-series feature sets for each pixel using the multispectral reflectance and multiple vegetation indices, and perform smoothing and quality screening processes to obtain pixel-level high-quality time-series feature sets for each pixel. Finally, downscale to the plot level and calculate plot-level high-quality time-series feature sets for each plot. Step 2: Perform a weighted average of multiple vegetation indices in the high-quality temporal feature set of each plot to obtain the comprehensive phenological curve of each plot. Based on the comprehensive phenological curve and actual phenological conditions of each plot, divide the clustering time window of each plot. Step 3: Using the high-quality temporal feature set of each land parcel located within the clustering time window of the land parcel as input, the feature-adaptive clustering algorithm for land parcels is used to iteratively optimize multiple cluster centers, and all land parcels are assigned to each cluster center to obtain unlabeled clusters centered at each cluster center. Step 4: Obtain candidate cluster sets for each crop category based on phenological consistency constraints; then, based on the sown area of ​​each crop category in the county-level statistical data, use a greedy algorithm to assign initial crop category labels to each unlabeled cluster to obtain an initial matching scheme; finally, iteratively optimize the initial matching scheme through local search optimization to obtain the optimal matching scheme. Step 5: Pass the crop category labels of the clusters in the optimal matching scheme to the effective pixels to construct the initial set of labeled pixels for each crop category; The initial labeled cell set for each crop category is preliminarily screened using a neighborhood weighting constraint, and weights are assigned to the effective cells in the labeled cell set. A spectral iterative self-cleaning method is used to construct class centers for each crop category based on the feature vectors and weights of effective pixels in the labeled pixel set of each crop category. Then, the labeled pixel set of each crop category is optimized and filtered through the class centers. Finally, the high-purity pixel set of each crop category is obtained by iteratively updating the labeled pixel set and class centers of each crop category. The feature vectors of the effective pixels are high-quality temporal feature sets at the pixel level.

2. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 1, characterized in that, The comprehensive phenological curve is obtained by weighting and fusing multiple vegetation indices from the high-quality temporal characteristic set of the plot to form a fused vegetation index, thereby obtaining the comprehensive phenological curve of the plot. The comprehensive phenological curve includes the fused vegetation index at each observation time point. Each vegetation index is weighted by its normalized Gini importance score. The clustering time window for the land parcels is obtained in the following way: The peak position of the comprehensive phenological curve of the plot is used as the peak time of the growing season of the corresponding plot, and the number of peaks is identified as the number of times the plot is replanted. Step 2.2.1: Pre-divide the growing season based on the crop growth history of the study area; Step 2.2.2: Calculate the target threshold based on the pre-divided growing season and comprehensive phenological curve. The target threshold is calculated based on the following formula: ； In the formula, The target threshold; This represents the minimum value of the integrated vegetation index of the comprehensive phenological curve within the pre-divided growing season period; This represents the maximum value of the vegetation index within the pre-divided growing season period of the composite phenological curve. The set scaling factor; Step 2.2.3, when there is a continuous When the fused vegetation index at each observation time point is greater than or equal to the target threshold, then... The first day of each observation point marks the start of the growing season. When continuous When the fused vegetation index at each observation time point is less than or equal to the target threshold, then... The first day of each observation point marks the end of the growing season. Then the clustering time window of the land parcels .

3. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 2, characterized in that, The process of obtaining the candidate cluster set for each crop category based on phenological consistency constraints specifically includes the following steps: Step 4.1: Select a set of candidate clusters for each crop type, where the first... Candidate cluster set of crop-like organisms Specifically, it is obtained through the following steps: Step 4.1.1: Compare the key phenological parameters of each cluster with the first... The key phenological parameters of similar crops are compared. When the comparison result is less than a set time tolerance, the cluster is determined to satisfy the phenological consistency constraint. By traversing all clusters, all clusters that satisfy the first phenological constraint are obtained. Clusters constrained by phenological consistency of similar crops; The key phenological parameters include the start time of the growing season, the peak time of the growing season, and the end time of the growing season. The start time, peak time, and end time of the growing season for each cluster are the median of the corresponding parameters for all plots within the cluster. Step 4.1.2, calculate the first... Initial prior feature centers of crop-like plants : In satisfying the first In the clusters constrained by phenological consistency of similar crops, calculate the key phenological parameters of each cluster and the first... The sum of the deviations between the standard key phenological parameters of similar crops, and the selection of the crop with the smallest sum of deviations. The first cluster is used as the initial reference cluster; then, the clusters in the initial reference cluster are... The eigenvectors of the i-th cluster are averaged to obtain the i-th cluster. Initial prior feature centers of crop-like plants ; The feature vector of the cluster is obtained by taking the median of the high-quality features at each observation time point of all plots within the cluster to form the feature vector of the cluster. Step 4.1.3, add all that satisfy the first... The eigenvectors of clusters with phenological consistency constraints for crop-like crops and the first Initial prior feature centers of crop-like plants Calculate similarity, based on the number of... Initial prior feature centers of crop-like plants similarity From largest to smallest, for all conditions satisfying the first... The clusters of crops with phenological consistency constraints were sorted, and the top clusters were selected. The clusters constitute a candidate cluster set. ; Step 4.1.4: Repeat step 4.1.3 to select candidate cluster sets for each type of crop, and then construct the candidate category set corresponding to each cluster by the crop category corresponding to the candidate cluster set to which each cluster belongs.

4. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 3, characterized in that, The process of assigning initial crop category labels to each unlabeled cluster using a greedy algorithm specifically includes the following steps: Step 4.2.1: Sort the clusters by area from largest to smallest, and initialize the remote sensing estimated area for each crop category. =0, where, the first The remote sensing estimated area of ​​crop-like crops is ,initialization , Pick , The total number of crop categories; Step 4.2.2: Perform initial crop category label assignment for clusters: Each time initial crop category label assignment is performed for a cluster, the current cluster is sequentially assigned from the candidate category set, and the overall loss of the current assignment is calculated for each assignment. Finally, the total loss relative to the initial allocation is taken. The crop category with the largest decrease is the initial crop category label of the current cluster, if relative to the overall loss of the initial assignment. If none of them decrease, then the crop category initially classified is used as the initial crop category label for the current cluster. Simultaneously update the remote sensing estimated area for the corresponding crop category. Then proceed with the initial crop category label assignment for the next cluster; Total loss Calculated based on the following formula: ； In the formula, Crop category; The first in county-level statistical data Planting area of ​​similar crops; To take the absolute value; Step 4.2.3: Traverse all clusters in descending order of area, repeating step 4.2.2 until all clusters have completed the initial crop category label assignment. All unlabeled clusters and their corresponding initial crop category labels constitute the initial matching scheme.

5. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 4, characterized in that, The process of using neighborhood weighting constraints to initially filter the initial labeled pixel set for each crop category to obtain the labeled pixel set for each crop category specifically includes the following steps: Step 5.1.1: First, construct a spatial neighborhood window of a preset size centered on each valid cell in the initial labeled cell set. The window is located at the first valid cell in the initial labeled cell set. The spatial neighborhood window of each effective cell is ; Then, count the number of valid pixels with the same crop category label as the corresponding valid pixel within the spatial neighborhood window of each valid pixel in the initial labeled pixel set; simultaneously, count the total number of valid pixels within the spatial neighborhood window of each valid pixel in the initial labeled pixel set. Then, the... The proportion of similar effective pixels Calculated based on the following formula: ； In the formula, For the first The number of valid pixels with the same crop category label as the corresponding valid pixel within the spatial neighborhood window of each valid pixel. For the first The total number of valid pixels within the spatial neighborhood window of each valid pixel; Finally, based on the 8-neighborhood connectivity rule, connectivity analysis is performed on other valid pixels within the spatial neighborhood window of each valid pixel in the initial labeled pixel set that have the same crop category label as the valid pixel. This constructs the connectivity of each valid pixel, and the number of other valid pixels within the connectivity of each valid pixel is calculated. Among these, the 8-neighborhood connectivity rule... The number of other valid pixels within the connected component of each valid pixel is ; Step 5.1.2: Set the first and second similarity ratio thresholds as follows: and The connectivity threshold is ; Will and The effective pixels are determined to be high-confidence effective pixels, and and The valid pixels are determined to be valid pixels with medium confidence. and The valid pixels were determined to be noisy pixels with low confidence. In the initial set of labeled pixels, retain the valid pixels with high and medium confidence, and remove the noise pixels with low confidence to obtain the set of labeled pixels. The specific method for assigning weights to valid pixels in the labeled pixel set is as follows: assigning a weight of 1 to high-confidence valid pixels in the labeled pixel set, and assigning a weight equal to the proportion of pixels in the same category to medium-confidence valid pixels in the labeled pixel set; the first pixel in the labeled pixel set... The weight of each valid cell is ,but or ,in, For the first The proportion of similar effective pixels.

6. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 5, characterized in that, The high-purity pixel sets for each crop category are obtained through the following steps: Step 5.2.1: Construct class centers for each crop category based on the feature vectors and weights of effective pixels in the labeled pixel set for each crop category. Crop-like class center Calculated based on the following formula: ； in, For the first The weight of each effective cell, For the first Feature vectors of 1 valid pixel; Serial number , This represents the number of valid pixels in the tagged pixel set. Step 5.2.2: Calculate the distance between each valid pixel in the labeled pixel set of each crop category and the class center of the corresponding crop category, where the distance is... The first tagged cell in the crop category set The effective pixel and the first Distance between the class centers of crop-like plants for: ,in, It is the Euclidean norm; Filtering thresholds for labeled cell sets of each crop category for: ,in, and This represents the average and standard deviation of the distances from each valid cell in the labeled cell set to the class center. Step 5.2.3: Iteratively update the labeled cell set: delete cells from the labeled cell set that meet the following conditions. The effective pixels are used to iteratively update the set of labeled pixels, thus obtaining the updated set of labeled pixels for the current iteration; Step 5.2.4: Repeat steps 5.2.1 to 5.2.3 for the next iteration using the updated labeled cell set; until the updated labeled cell set no longer changes during iteration or the set maximum number of iterations is reached. At that time, the final high-purity pixel set is obtained.

7. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 6, characterized in that, The local search optimization includes swapping, moving, and fine-tuning; The exchange involves selecting two clusters and swapping their initial crop category labels; if the overall loss... If the value decreases, the exchange will be accepted. The movement involves replacing the current initial crop category label of a cluster with another crop category label from the corresponding candidate category set; if the overall loss... Decrease the value and allow movement; The fine-tuning involves: under the premise of satisfying the consistency constraint of the phenological window, calculating the area error between the sown area and the remotely estimated area of ​​the crop category in the county-level statistical data under the matching scheme of the current round of local search, sorting the area errors from largest to smallest, and calculating the cumulative area error and total error in sequence according to the sorting order. Clusters of crop categories belonging to the category with a cumulative area error greater than the total error by a set percentage are reassigned; the reassigned clusters are then placed in the candidate cluster set. The crop labels are reassigned to other crop category labels within the specified range, and the loss function is calculated. Change the label to select other crop category tags that will minimize losses; Each local search operation executes the above three operations sequentially. When several consecutive rounds of local searches fail to further reduce the speed... Or, when the maximum number of iterations is reached, the optimization stops and the optimal matching solution is obtained.

8. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 7, characterized in that, Step 3 specifically includes the following steps: Step 3.1: Divide the growing season into different growth periods, and group the plots according to the number of times the plots are replanted to obtain different plot groups; Step 3.2: Construct a similarity matrix between land parcels. Attraction Matrix and the degree of belonging matrix Similarity matrix Attraction matrix Belonging degree matrix ,in, For the plot of land With the plot of land Similarity between them For the plot of land land parcel The attractiveness For the plot of land land parcel degree of belonging Number of plots, serial number and All ; Step 3.3: In each land parcel group, determine the cluster centers through iteration, and assign all land parcels to the respective cluster centers to form unlabeled clusters. The center point of each cluster is the corresponding cluster center. This specifically includes the following steps: Step 3.3.1: Set the maximum number of iterations. After each iteration, set the bias parameter for the current iteration. The values ​​are directly assigned to the diagonal elements of the similarity matrix in the current iteration. This process yields the similarity matrix for the next iteration. The bias parameter for the next iteration is calculated based on the bias parameter from the current iteration and the median of the similarity matrix for the next iteration. The bias parameter from the first iteration uses the median of the similarity matrix before the first iteration. This process is then applied to the similarity matrix. Attraction Matrix and the degree of belonging matrix Perform iterative optimization; among which, For the plot of land Self-similarity; Step 3.3.2: After each iteration, based on the attraction matrix of the current iteration... With the degree of belonging matrix Constructing a decision matrix According to the decision matrix diagonal elements Determine candidate cluster centers; All plots with a value greater than 0 were recorded as candidate cluster centers, among which, For the plot of land Degree of belonging to oneself For the plot of land Attraction to oneself; When in succession In each iteration, if the candidate cluster centers remain unchanged, then the continuous... Candidate cluster centers that remain unchanged throughout the iterations are determined as cluster centers, and the iteration stops. If no cluster centers are found by the maximum number of iterations, the attraction matrix at the maximum number of iterations is used. With the degree of belonging matrix Calculate all Plots with a value greater than 0 are cluster centers; Step 3.3.3: Assign all land parcels to various cluster centers, forming unlabeled clusters centered at each cluster center. Specifically, this is based on the following formula: Assigned to cluster center: ； In the formula, In order to make To obtain the maximum value The value of , The index of the cluster center. Select from the numbers of the plots identified as cluster centers. For the plot of land Cluster centers degree of belonging Cluster center land parcel Its attractiveness.

9. The method for automatically generating crop samples using remote sensing-statistics-phenology coupling according to claim 1, characterized in that, The smoothing process is based on the following procedure: First, the number of globally optimal harmonic combinations for the harmonic regression model is obtained by minimizing the root mean square error; Secondly, a harmonic regression model with the globally optimal number of harmonic combinations is used to smooth the high-dimensional time series feature set of all pixels, so as to obtain the harmonic regression-based time series feature set of each pixel. The harmonic regression time series feature set of each pixel includes multiple harmonic regression features at each observation time point, and the number of harmonic regression features is equal to the number of features in the high-dimensional time series feature set. The quality screening process is based on the following steps: Step 1.4.1: Normalize the harmonic regression time series feature set of each pixel; Step 1.4.2: Based on Mahalanobis distance, remove the outlier observation time points and corresponding harmonic regression features from the normalized harmonic regression time series feature set of each pixel to obtain the preprocessed harmonic regression time series feature set of each pixel. Step 1.4.3: Use the Pearson correlation coefficient to remove harmonic regression features with a correlation higher than the set correlation threshold from the time series features of the preprocessed harmonic regression time series feature set; Then, using univariate statistical tests, all significant harmonic regression features with P values ​​less than the set P value threshold are retained to obtain the significant harmonic regression time series feature set for each pixel. The significant harmonic regression time series feature set includes multiple significant harmonic regression features at each observation time point. Step 1.4.4: Finally, input the saliency harmonic regression time series feature set of each pixel into the random forest model to calculate the Gini importance, sort them from high to low, and select the top... Each significant harmonic regression feature is a high-quality feature, resulting in a pixel-level high-quality temporal feature set for each pixel. The pixel-level high-quality temporal feature set for each pixel includes features from each observation time point. One high-quality feature; The high-quality temporal feature sets at the plot level for each plot are obtained through the following steps: Step 1.5.1: Count the number of observations that pass quality control for each pixel. If the proportion of the number of observations that pass quality control to the total number of observations at all time points of the corresponding pixel is not less than the set proportion threshold, then the corresponding pixel is determined to be a valid pixel, and all valid pixels in each plot are obtained. Step 1.5.2: For the pixel-level high-quality temporal feature set of all valid pixels within each plot, the median of the high-quality features at each observation time point is used for fusion to obtain the plot-level high-quality temporal feature set for each plot. The plot-level high-quality temporal feature set includes the features at each observation time point. High-quality characteristics at the plot level.

10. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the automatic generation method for crop samples by remote sensing-statistics-phenology coupling as described in any one of claims 1 to 9.

Images