Method for near real-time monitoring of phenology and type of mid-season crops based on dense time series

Abstract

This invention discloses a near-real-time monitoring method for mid-season crop phenology and types based on dense temporal series data, belonging to the interdisciplinary field of satellite remote sensing Earth observation and agricultural information technology. The method includes: acquiring HLS imagery to construct a dense temporal series dataset; recursively updating the dataset using the Stochastic Continuous Change Detection (S-CCD) algorithm; fitting the vegetation index temporal characteristic curve for each pixel based on the S-CCD results to predict the phenological stage of the current year; classifying pixels into stable, uncertain, and variable classes, with the stable and uncertain classes inheriting historical planting types; training a random forest model using automatically selected stable class samples to classify the variable classes, and performing Bayesian temporal fusion with previous probability maps to determine the final crop type; fusing the results with vector plots and automatically publishing them through a network service. This invention achieves high-precision and high-efficiency dynamic monitoring of phenology and types within the crop growth cycle, supporting near-real-time weekly updates.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of satellite remote sensing Earth observation and agricultural information technology, specifically involving a near real-time identification and monitoring method for mid-season crop phenological stages and crop types based on dense time-series remote sensing images and machine learning algorithms. Background Technology

[0002] Crop phenology and type identification refers to the process of automatically identifying and analyzing the types and growth stages of crops planted in farmland using remote sensing technology. Crop type identification can clearly identify the specific crop varieties planted in the field, such as wheat, rapeseed, and rice; phenology identification monitors key stages of crop growth and development (such as emergence and maturity). These two technologies are core components of precision agriculture and have significant application value for food security early warning, agricultural subsidy distribution, and farmland disaster assessment. To implement scientific agricultural management decisions, timely and accurate information on crop planting structure and growth status is necessary. Accurate identification of crop phenology and types is not only a crucial link in agricultural remote sensing monitoring but also an important foundation for yield forecasting, pest and disease control, and water and fertilizer regulation.

[0003] Compared to traditional crop monitoring methods based on manual field surveys, satellite remote sensing technology offers advantages such as wide coverage, strong periodicity, and low cost, making it a primary means of modern agricultural monitoring. Remote sensing images contain rich spectral features and temporal variation information, which can be used to effectively identify the growth characteristics of different crops. Early crop identification relied mainly on expert visual interpretation, analyzing image features such as color and texture and combining this with field experience. However, this method is inefficient and requires highly specialized expertise from interpreters, making it difficult to meet the needs of large-scale, high-frequency crop monitoring. With the development of remote sensing technology, digital methods such as time-series vegetation index analysis, phenological feature extraction, and statistical classification methods have been gradually applied to crop identification. While these methods have improved identification efficiency, they suffer from poor adaptability when dealing with complex planting structures and insufficient utilization of temporal features, making it difficult to fully extract crop growth information from high-resolution satellite imagery.

[0004] Despite significant advancements in existing technologies, challenges remain in achieving near real-time, high-precision monitoring of crop phenology and crop type dynamics. First, efficiently processing massive amounts of time-series remote sensing data and accurately capturing subtle changes in crop growth is a key challenge. Second, crop type changes and natural phenological changes exhibit similarities in time-series curves, making effective differentiation difficult. Furthermore, existing methods typically require reclassification across the entire region, resulting in low computational efficiency and failing to meet the demands of operational real-time monitoring. Therefore, developing an intelligent monitoring method capable of automatically, accurately, and efficiently distinguishing changing areas and fully utilizing historical prior knowledge is crucial for promoting the development of smart agriculture. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a near real-time monitoring method for crop phenology and type based on dense temporal remote sensing imagery and an improved Stochastic Continuous Change Detection (S-CCD) algorithm. This invention can accurately monitor and determine the phenological stages and planting types of crops, updating weekly with the latest imagery to help understand real-time crop growth and crop rotation. Compared to traditional interference monitoring algorithms, the recursive computation characteristic of S-CCD based on a Kalman filter eliminates the need to refit the entire time series each time new data is acquired. This significantly improves the computational efficiency of the algorithm, meeting the needs of large-scale real-time monitoring. Simultaneously, it automatically extracts random forest samples for determining change types from stable classes, enabling the entire process to be automated. Furthermore, this invention innovatively introduces an "inheritance-update" strategy, which can effectively distinguish between crop type changes and natural phenological changes, reclassifying only areas where real changes have occurred, thereby significantly improving monitoring efficiency while maintaining high accuracy.

[0006] The specific technical solution adopted in this invention is as follows: This invention provides a near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences, as detailed below: S1: Acquire multispectral remote sensing image data of the area to be monitored from the past to the current time phase, divide the image data into blocks and stack them in chronological order to construct a time series dataset, process the time series dataset using a random continuous change detection algorithm, and recursively update the detection result for newly added image data; S2: Based on the harmonic coefficients in the random continuous change detection results described in S1, construct and predict the characteristic curve for each pixel for the current year, analyze the phenological stage based on the threshold method and the enhanced vegetation index estimate, and generate a raster phenological product that includes multiple cropping index, crop phenological stage and change confidence. S3: Optimize the parameters of the random continuous change detection described in S1, and then classify the change class pixels, uncertain class pixels, and stable class pixels based on the obtained breakpoints and outliers; combine the uncertain class pixels and stable class pixels with the multiple cropping index of S2 to inherit the corresponding historical crop type; S4: Automatically select samples from the stable class pixels described in S3, extract the current time period features and train a random forest classification model to classify crop types for the changing class pixels, and perform Bayesian fusion with the historical crop type posterior probability map to obtain the current crop type map. S5: Spatially fuse and assign attribute values ​​to the raster phenology product described in S2 and the current crop type map described in S4 with the vector plot data of the area to be monitored, construct a large-scale, efficient, block-based parallel computing module, and realize the automated release and weekly update of monitoring results based on the block-based parallel computing and network map service interface.

[0007] Preferably, in S1, the multispectral bands include visible light, near-infrared band, short-wave infrared band, and quality assessment band.

[0008] Preferably, in step S1, to enhance the expressive power of temporal features, the enhanced vegetation index (EVI) of the multispectral remote sensing image data is first calculated; then, a single remote sensing image data is divided into multiple spatial blocks and stacked in chronological order to construct a dense temporal dataset for temporal analysis and phenological identification.

[0009] Preferably, in step S1, the random continuous change detection algorithm is performed after cloud masking. Based on the state-space model and the Kalman filter recursive update mechanism, harmonic modeling and change detection are performed on the temporal vegetation index of each pixel. The harmonic model obtained through harmonic modeling is used to simulate the annual and semi-annual cycle signals of crop growth. The specific formula for the harmonic coefficients of the harmonic model is as follows: ; ; in, The periodic frequency; For time indexing; The main periodic state variables at time t; These are the principal periodic state variables at time t+1; Main periodic noise; For the main periodic state variables The process noise variance; The auxiliary state variable at time t Let be the auxiliary state variable at time t+1; As auxiliary periodic noise; Auxiliary periodic state variables The process noise variance.

[0010] Preferably, in step S3, the parameter optimization method is as follows: A loose threshold is used to perform preliminary breakpoint detection on historical time-series data, and positive samples caused by changes in actual crop type and negative samples caused by noise or phenological changes are selected to form a parameter optimization sample set. Using the F1 score as the evaluation index, the optimal parameter combination is obtained by searching in the preset random continuous change detection parameter space through the grid search method for the final pixel category classification.

[0011] Preferably, in step S4, when automatically selecting samples from stable class pixels, a grid sampling method based on spatial heterogeneity is used, as follows: The area to be monitored is divided into R×C uniform grids. The coefficient of variation (CV) of the EVI values ​​of all pixels in each grid is calculated and classified according to heterogeneity. More sample points are allocated in highly heterogeneous grids to ensure the spatial representativeness of the samples and the balance of class distribution.

[0012] Preferably, in step S4, the method for training the random forest classification model to classify crop types for change-type pixels is as follows: First, based on the stable class pixels divided by S3, the latest observable values ​​of the six bands of green light, red light, near infrared, SWIR1, SWIR2 and enhanced vegetation index (EVI) and their corresponding crop category labels are extracted to construct a training sample set; then, the random forest classification model is trained using the training sample set, and the optimal split point is determined according to the principle of maximizing the reduction of Gini index when splitting nodes. For the change class pixels to be classified, extract their six common band feature values, input them into the trained random forest classification model, and determine their final crop category through the majority voting mechanism of multiple decision trees; in view of the large number of change class pixels in large-scale remote sensing data, a batch prediction mechanism is adopted to avoid memory overflow.

[0013] Preferably, in S4, the Bayesian fusion method is as follows: Using the category probability distribution of each pixel in the historical crop type map as the prior probability and the current observation likelihood probability map as the likelihood probability, the posterior probability distribution of all pixels in the region is calculated using the Bayesian update formula; the Bayesian update formula is used to calculate the posterior probability distribution of all pixels in the region for each pixel belonging to category C. i Updated posterior probability P posterior (C i The calculation formula is: ; Among them, P prior (C i Let P be the prior probability. likelihood (C i) represents the likelihood probability. j Iterate through all crop categories.

[0014] Preferably, in step S5, when spatially fusing and assigning attributes to raster phenological products and vector plot data, the dominant crop type within each vector polygon is determined using the mode voting method, and the average confidence level, most recent observation time, and crop rotation type of the pixels within that vector polygon are simultaneously calculated as attribute information; the dominant crop type C dominant Determined by the following formula: ; Where, Count Ci For the vector polygon, the polygon belongs to category C. i The number of raster cells, and the average confidence level is the arithmetic mean of the confidence levels of all cells within the vector polygon.

[0015] Preferably, in S5, automated publishing is achieved by calling the REST API of the geographic information server, including the entire process of automatically deleting old version data storage, uploading new version raster result files, configuring and publishing new layers.

[0016] Compared with the prior art, the present invention has the following advantages: This invention possesses five significant advantages: near real-time performance, high efficiency and high precision, strong anti-interference and discrimination capabilities, full-process automation, and strong practical applicability of the results. First, through the incremental update mechanism of the Random Continuous Change Detection (S-CCD) algorithm, the latest remote sensing image data can be quickly integrated to achieve near real-time dynamic monitoring of crop growth status and planting structure on a weekly basis, effectively meeting the high standards of timeliness required by operational applications.

[0017] Second, it innovatively proposes an "inheritance-update" strategy. By accurately dividing the pixels, most unchanged areas directly inherit reliable historical results, while only a small number of changed areas are reclassified. This significantly reduces the amount of computation and improves processing efficiency, while ensuring the accuracy and reliability of the overall monitoring results.

[0018] Third, relying on a specialized parameter optimization process, the improved S-CCD algorithm can effectively distinguish between "crop type changes" and "natural phenological changes / noise," significantly reducing the false alarm rate and improving the accuracy of change detection.

[0019] Fourth, it achieves fully automated processing of the entire chain from data acquisition, preprocessing, time series analysis, change detection, machine learning classification to result publication and visualization, which greatly reduces manual intervention, lowers operating costs, enhances the repeatability and stability of the method, and is easy to integrate into business operation systems.

[0020] Fifth, the results are deeply integrated with vector plots and can be directly published through online map services, greatly facilitating intuitive querying and use by agricultural management departments and related users, resulting in outstanding decision support effects. Attached Figure Description

[0021] Figure 1 This is the overall flowchart of the method of the present invention; Figure 2 This is a time series diagram stage classification (taking the wheat + rice rotation pattern as an example); Figure 3 It is an S-CCD fitting image of pixels in the changing, unstable and stable categories; Figure 4 This is a webpage display of the phenological results of the method of this invention; Figure 5 This is a webpage display of the classification results of the method of this invention. Detailed Implementation

[0022] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in various embodiments of the present invention can be combined accordingly without mutual conflict.

[0023] Crop phenology and type monitoring is fundamental to precision agricultural management. However, due to insufficient utilization of temporal characteristics and difficulty in distinguishing between crop type changes and phenological changes, real-time monitoring of crop growth dynamics based on remote sensing is currently difficult. Therefore, this invention constructs a near-real-time monitoring method for mid-season crop phenology and type based on dense temporal sequences. Its core lies in achieving accurate temporal change detection and pixel partitioning through an improved S-CCD algorithm, and adopting differentiated processing strategies for different regions.

[0024] This invention innovatively proposes a classification strategy based on inheriting invariant points. Using an improved S-CCD algorithm, pixels are precisely divided into three categories: stable, uncertain, and change. For stable and uncertain categories, historical planting types are directly inherited; for change categories, a high-quality sample selected from the stable category drives a random forest model for fine-grained classification, achieving crop type identification. This innovative method fully utilizes the continuous characteristics of farmland planting in the spatiotemporal dimensions. By inheriting reliable classification results from the past, it effectively reduces the number of samples requiring reclassification, significantly improving identification efficiency while maintaining classification accuracy. Experimental results show that this inheritance strategy can effectively utilize historical prior knowledge, especially in areas with relatively stable crop planting structures, significantly improving the stability and reliability of classification.

[0025] like Figure 1 As shown, the specific steps of the method of the present invention are as follows: S1: Acquire multispectral remote sensing image data of the area to be monitored from historical to current time phases, divide the obtained image data into blocks and stack them in chronological order to construct a time series dataset, process the obtained time series dataset using the Random Continuous Change Detection (S-CCD) algorithm, and recursively update the detection result for newly added image data.

[0026] S-CCD is a novel forest disturbance monitoring algorithm that incorporates two mathematical tools: a state-space model and a Kalman filter. Unlike traditional methods, the state-space model can better capture the complex nonlinear dynamic changes in vegetation time series, rather than simply linear trends and periodic variations. This makes S-CCD more sensitive and accurate in monitoring forest disturbance events that cause subtle spectral changes. Simultaneously, the recursive computational nature of the Kalman filter allows S-CCD to quickly and efficiently update model parameters each time new observation data arrives. This significantly improves the algorithm's computational efficiency, making it more suitable for real-time monitoring needs.

[0027] In a preferred embodiment of the present invention, the area to be monitored is Haiyan County, Jiaxing City, Zhejiang Province, China.

[0028] In a preferred embodiment of the present invention, the multispectral bands include visible light, near-infrared, shortwave infrared (SWIR), and quality assessment (QA) bands, with a spatial resolution of 30 m. To enhance the expressive power of temporal features, an Enhanced Vegetation Index (EVI) is further calculated based on the near-infrared, red, and blue light bands to reflect crop growth status. A single remote sensing image is divided into multiple spatial blocks and stacked chronologically to construct a dense temporal dataset for temporal analysis and phenological identification. The EVI calculation formula is as follows: ; Where, ρ NIR ρ RED ρ BLUE These are the surface reflectances in the near-infrared, red, and blue light bands after atmospheric correction; 2.5 is the gain coefficient, 6 and 7.5 are the atmospheric impedance coefficients, and 10000 is the soil adjustment parameter (all band values ​​have been magnified ten thousand times).

[0029] In practical use, the original image should first be divided into blocks. For example, if a single image has a resolution of 3660 * 3660, it should be cropped into 30 * 30 blocks of 122 * 122 resolution images and stacked in chronological order to build a dense temporal data cube. This provides a unified data structure for subsequent temporal analysis and modeling, and new image data should be added for recursive update processing.

[0030] In a preferred embodiment of the present invention, to accurately capture crop growth time sequence characteristics, the parameters in the S-CCD processing of S1 are as follows: tmask_b1 is set to the green band, tmask_b2 is set to the SWIR1 band, lam is set to 20 to balance model complexity and goodness of fit, p_cg is set to 0.9999 to control the false detection rate, cons is set to 8 to ensure the reliability of change detection, anomaly_pcg is set to 0.90 to control the false detection rate of anomalies, and anomaly_conse is set to 3 to control the reliability of anomaly detection.

[0031] As a preferred embodiment of the present invention, after the cloud masking process is completed, the S-CCD algorithm based on the state-space model and the Kalman filter recursive update mechanism is used to perform harmonic modeling and change detection on the temporal vegetation index of each pixel.

[0032] Harmonic modeling is used to simulate periodic signals in crop growth (such as annual and semi-annual cycles). The state equation of the harmonic model is defined by a pair of coupled state variables, and its update formula is as follows: ; ; in, The periodic frequency; For time indexing; The main periodic state variables at time t; These are the principal periodic state variables at time t+1; Main periodic noise; For the main periodic state variables The process noise variance; The auxiliary state variable at time t Let be the auxiliary state variable at time t+1; As auxiliary periodic noise; Auxiliary periodic state variables The process noise variance.

[0033] In this harmonic model, the angular frequency Used to define the period length (e.g., annual period T=1, semi-annual period T=0.5), while independent noise terms and (variance This endows the model with the ability to dynamically adjust, through and The linear combination simulates the annual / seasonal periodic signals of crop growth, while the noise term allows the periodic amplitude to vary over time to accommodate non-stationary time-series data; auxiliary variables It is used only to maintain the mathematical closure of recursive calculations and does not directly participate in observation prediction (its coefficients in the corresponding observation matrix Z are 0). In practical applications, the state vector... It can simultaneously model trends, annual cycles, and semi-annual cycles, and update them recursively. To accurately estimate and track key phenological stages of crop growth (such as the greening-up stage and the maturity stage).

[0034] This differs significantly from traditional methods (such as COLD), which require rebuilding the model each time. In traditional methods, the entire time series model needs to be reconstructed using the LASSO regression algorithm every time new observation data is added. S-CCD, however, only needs to update the state variables and covariance matrix, without reloading all historical image data. This recursive update method not only greatly improves computational efficiency but also enables S-CCD to achieve near real-time monitoring, as it can detect changes without waiting for all historical data to be loaded.

[0035] S2: Based on the harmonic coefficients in the random continuous change detection results of S1, construct and predict the characteristic curve of the year for each pixel, analyze the phenological stage based on the threshold method and the enhanced vegetation index estimate, and generate a raster phenological product that includes multiple cropping index, crop phenological stage and change confidence.

[0036] In a preferred embodiment of the present invention, an EVI feature curve is fitted to each pixel based on the S-CCD results and predicted up to the end of the year. Based on a threshold method combined with the predicted seasonal saturation EVI value and the latest observed value, the crop rotation stage and phenological stage of the year are determined, including four stages: dormancy, greening-up stage, maturity stage, and senescence stage. Figure 2 It calculates the confidence level of the current phenological stage and finally outputs a raster phenological product (i.e., raster result) that includes the multiple cropping index, crop phenological stage and change confidence level.

[0037] S3: Optimize the parameters of the random continuous change detection in S1, and then classify the change class pixels, uncertain class pixels, and stable class pixels based on the obtained breakpoints and outliers; combine the uncertain class pixels and stable class pixels with the multiple cropping index of S2 to inherit the corresponding historical crop type.

[0038] As a preferred embodiment of the present invention, the parameter optimization method is specifically as follows: To optimize the S-CCD algorithm's sensitivity to detecting "actual crop type changes" and reduce false alarms caused by "phenological changes" and "noise," parameter sensitivity analysis is required, specifically as follows: A relaxed threshold is used to perform preliminary breakpoint detection on historical time-series data. Positive samples caused by actual crop type changes and negative samples caused by noise or phenological changes are manually selected to form a parameter optimization sample set. Using the F1 score as the evaluation index, a grid search method is used to optimize within a pre-defined random continuous change detection parameter space, obtaining the optimal parameter combination for final pixel category classification. The formula for calculating the F1 score is: ; In the formula, True Positive (TP): The number of samples that are actually "Type Change" and are also judged as "Type Change" by the algorithm; False Positive (FP): The number of samples that are actually "Other Factors" but are misjudged as "Type Change" by the algorithm; False Negative (FN): The number of samples that are actually "Type Change" but are missed by the algorithm (i.e., judged as "Other Factors").

[0039] The specific method is as follows: First, a set of low detection thresholds (such as low const and p_cg) is set to initially screen out a large number of pixels with breakpoints. Then, through manual visual interpretation or time-series curve morphology analysis, 100-200 samples with breakpoints caused by changes in real crop types, as well as samples with breakpoints caused by other factors (cloud noise, phenological fluctuations), are selected to form a parameter optimization dataset (e.g., 100 positive samples "caused by crop changes" and 50 negative samples "caused by noise"). A grid search method is used to traverse within a preset parameter space (e.g., const: 5-8, p_cg: 0.9-0.9999), using the F1 score as the evaluation index, to find the optimal parameter combination that can most effectively distinguish real type changes, thereby improving the sensitivity and specificity of S-CCD to changes in real planting systems. Then, S-CCD is run with these parameters, and pixels of the change class are divided according to whether breakpoints occur, and pixels of the uncertain class are divided according to whether there are outliers. The remaining pixels are the stable class pixels. Figure 3 ).

[0040] S4: Automatically select samples from the stable class pixels divided by S3, extract the features of the current time period, and train a random forest classification model to classify crop types for the changing class pixels. Then, perform Bayesian fusion with the historical crop type posterior probability map to obtain the current crop type map.

[0041] Specifically, a random forest model is trained by automatically generating sample features for the current time period using historical stable points. This model then classifies and predicts the changing class pixels, resulting in a current observation likelihood probability map covering the entire map. Subsequently, a unified Bayesian update framework is constructed, using historical crop probability maps as prior knowledge, integrating current observations, calculating the updated posterior probability, and determining the crop type with the highest probability to complete the planting type map update.

[0042] In a preferred embodiment of the present invention, when automatically selecting samples from stable class pixels, a grid sampling method based on spatial heterogeneity is used to automatically generate training samples from historical stable class pixels. Specifically: The area to be monitored is divided into R×C uniform grids. The coefficient of variation (CV) of all pixel EVI values ​​within each grid is calculated and classified according to heterogeneity. More sample points are allocated in highly heterogeneous grids to ensure spatial representativeness and balanced class distribution. The formula for calculating the coefficient of variation (CV) is: ; Where σ is the standard deviation of the EVI values ​​of all pixels in the grid, and μ is its average value.

[0043] For example, a 3660×3660 monitored area can be divided into a 30×30 grid (each grid contains 122×122 pixels), and the coefficient of variation (EVI) value within each grid can be calculated. Grids with an EVI greater than 0.3 are marked as high heterogeneity, 0.1-0.3 as medium heterogeneity, and less than 0.1 as low heterogeneity, thus classifying the grids into three heterogeneity levels: high, medium, and low. Double samples are drawn from high-heterogeneous grids, baseline samples are drawn from medium-heterogeneous grids, and half samples are drawn from low-heterogeneous grids.

[0044] Sample sizes are allocated based on heterogeneity, with more samples assigned to grids exhibiting high heterogeneity. Considering that other categories typically account for a larger proportion (over 80%) than major crops like rice, wheat, and rapeseed in actual monitoring, a moderately weighted strategy is employed during sample extraction: the sample proportion is appropriately increased for a few categories (rice, wheat, rapeseed) to ensure the model can fully learn their features, while sufficient samples are extracted for most categories (others) at a proportion close to reality to cover their internal diversity. The impact of class imbalance on the classification model is mitigated through class weight adjustment or oversampling techniques. For example, the sample size is set at 2000 for rice, 1000 for wheat, 300 for rapeseed, and 5000 for other categories, totaling 8300 sample points, ensuring representativeness in both spatial and class distribution. Only stable class pixels are used as training samples to ensure the reliability of sample labels. These are then used to train a random forest classifier, generating probability maps for rice, wheat, rapeseed, and other categories to classify areas of cropping system change.

[0045] In a preferred embodiment of the present invention, the machine learning classification method used for the variable-class pixels is the random forest algorithm. Specifically, based on the stable-class pixels divided by S3, the latest observable values ​​of their six bands (green, red, near-infrared, SWIR1, SWIR2, and Enhanced Vegetation Index (EVI)) and their corresponding crop category labels are extracted to construct a training sample set. Subsequently, the random forest classification model is trained using the training sample set, and the optimal split point is determined based on the principle of maximizing the reduction of the Gini index when splitting nodes. The formula for calculating the reduction of the Gini index ΔGini is: ; ; Where S is the node dataset, and K is the total number of categories. It represents the proportion of category k in the current node, A is the candidate feature, v is the split point, and S is the split point. left and S right The left and right subsets are generated after splitting. For pixels of the change class to be classified, their six common band feature values ​​are extracted and input into a pre-trained random forest classification model. The final crop class is determined through a majority voting mechanism of multiple decision trees. To address the large number of change class pixels in large-scale remote sensing data, a batch prediction mechanism is adopted to avoid memory overflow. Furthermore, the model can evaluate the importance of each band feature based on out-of-bag data.

[0046] For example, the input feature matrix X∈R N×6, where N is the number of training samples, and the six features are the latest observable values ​​of six bands: green, red, near-infrared, SWIR1, SWIR2, and EVI. The total number of decision trees is set to T=150. Feature subsampling is performed during single-tree modeling, randomly selecting m=⌊√6⌋=2 candidate features from the six features. Each decision tree constructs a training set through Bootstrap sampling, where approximately 63.2% of the original samples are selected for training, and the remaining 36.8% are used as out-of-bag data for model evaluation. At each node of the decision tree, the optimal split point is determined based on the principle of minimizing the Gini index from the two randomly selected candidate features. The largest feature and the split point are used as the optimal split. Recursively split until the stopping condition is met: the number of node samples is less than min_samples_split=10 or the maximum depth max_depth=20 is reached.

[0047] In practical applications, the hyperparameters of the random forest classifier can be set as follows: number of decision trees n_estimators=150, maximum depth max_depth=20, minimum number of split samples per node min_samples_split=10, minimum number of samples per leaf node min_samples_leaf=5, and the number of features randomly selected per tree split max_features='sqrt' (i.e., ⌊√6⌋=2 features). Parallel computation is enabled during model training n_jobs=-1 to fully utilize multi-core CPU resources.

[0048] As a preferred embodiment of the present invention, to address the problem that direct prediction may lead to memory overflow due to the large number of change-type pixels, a batch prediction mechanism is adopted. Specifically, the change-type pixels to be classified are divided into multiple batches, and the latest observable values ​​of their six bands are extracted to form a batch feature matrix X. batch The input is a pre-trained random forest model for prediction. Predictions for batch samples by a single decision tree are achieved through recursive traversal of the decision tree, and the ensemble prediction results are determined by majority voting. ; in, For sample x in batch b b The final predicted category, h t Let be the prediction function for the t-th decision tree. After all batches of predictions are completed, the results are stitched together to form a complete crop type map.

[0049] For example, by setting the batch size to 50,000, the pixels of the change type are divided into multiple batches, and the observation values ​​of their six bands are extracted batch by batch to form the feature matrix X. batch ∈R 50000×6Input the pre-trained random forest model. Internally, the model recursively traverses 150 decision trees, each independently predicting the category. Finally, the crop type for each pixel is determined by majority vote. ; Where, x i h represents the six-band feature vectors of the i-th change class pixel. t Let I(·) be the prediction function for the t-th decision tree, and let I(·) be the indicator function.

[0050] In a preferred embodiment of the present invention, after the model training is completed, the importance of the six band features is evaluated based on out-of-bag data. For each decision tree, the classification accuracy Acc is calculated using its out-of-bag data. t Subsequently, feature X j The values ​​are randomly arranged, and the accuracy is calculated again. Feature X j Importance is measured by the degree of decrease in accuracy: ; Importance scores can identify key bands that contribute most to crop classification, providing a basis for feature selection and model interpretation. For example, in this embodiment, the importance analysis results show that the near-infrared band and SWIR1 value contribute the most to crop classification, with importance scores of 0.21 and 0.18, respectively. This analysis result provides a basis for subsequent feature optimization.

[0051] In a preferred embodiment of the present invention, the category probability distribution of each pixel in the historical crop type map is used as the prior probability, and the current observed likelihood probability map is used as the likelihood probability. The posterior probability distribution of all pixels in the region is calculated using the Bayesian update formula. For each pixel, the Bayesian update formula determines its category C. i Updated posterior probability P posterior (C i The calculation formula is: ; Among them, P prior (C i ) is category C i The prior probability, P likelihood (C i ) is category C i The likelihood probability, j Iterate through all crop categories.

[0052] The specific fusion method is as follows: The posterior probability distribution map of crop categories calculated and generated in the previous monitoring period is used as the historical prior probability benchmark P. priorThis image records the probability value of each pixel belonging to each crop category. Simultaneously, the pixel stability labels from the previous period are retrieved. For pixels labeled "stable," their historical classification results have high confidence; therefore, their prior probabilities are reset to a deterministic distribution, i.e., the probability of their historical category is set to 1, and the probability of other categories is 0. For pixels labeled "uncertain" or "changing," their prior probabilities remain as their original probability distribution vectors. Using the remote sensing image features of the current time period, the probability of each pixel belonging to each crop category is obtained through a random forest model, forming the likelihood probability map P of the current observation. likelihood For each pixel in the monitored area, its specific prior probability vector and the current likelihood probability vector are substituted into Bayes' theorem to calculate the updated posterior probability vector P. posterior This calculation is performed in parallel across the entire map. The effect of this step is that stable pixels largely inherit historical results due to their deterministic priors, while the probability distribution of unstable pixels is dynamically adjusted based on new and old evidence. The updated posterior probability map undergoes spatial consistency processing (e.g., 3×3 Gaussian filtering) to smooth noise, and then the maximum probability method is used, selecting the category with the highest posterior probability value in each pixel as its final crop type, generating the final crop type classification map. Simultaneously, the posterior probability distribution of crop categories for each pixel calculated in this cycle (stored in multi-band raster form, with each band corresponding to the probability value of one crop) will serve as the prior probability input for the next monitoring cycle, thereby achieving iterative accumulation and transmission of classification knowledge at the probabilistic level.

[0053] This step constructs a unified Bayesian probability update framework, the core of which is to rigorously integrate historical knowledge with current observations at the probabilistic level. Specifically: (1) Setting the prior probability: For stable class pixels in the historical classification results, their crop type has been determined with high confidence. Therefore, their prior probability is set to a deterministic distribution - that is, the probability of their class is 1 and the probability of other classes is 0. This setting ensures that after the Bayesian update, unless the current observation provides strong contrary evidence (likelihood probability is 0), the posterior probability of the pixel will remain stable and inherit the historical classification results.

[0054] (2) Obtaining the likelihood probability: Using the remote sensing image features of the current time period, a random forest classifier is used to predict all pixels to obtain the probability that each pixel belongs to a certain type of crop, which is used as the likelihood probability of the current observation.

[0055] (3) Bayesian update calculation: For each pixel in the entire region, its specific prior probability distribution and the current likelihood probability distribution are substituted into Bayes' formula to calculate the updated posterior probability distribution. This process is mathematically represented as: ; in Category C i The posterior probability, P likelihood (C i ) is category C i The likelihood probability, P prior (C i ) is category C i The prior probability.

[0056] (4) Update Results: The ingenuity of this framework lies in the fact that it naturally achieves differentiated update strategies through different prior probability settings. Stable pixels remain stable due to their prior certainty; pixels with uncertain or changing histories, because their priors are probability distributions, can fully absorb current observation evidence and achieve substantial probability updates. This ensures the temporal consistency of classification results while allowing real changes to be detected and corrected in a timely manner.

[0057] S5: Spatially fuse and assign attributes to the raster phenology products generated in S2 and the current crop type map generated in S4 with the vector plot data of the area to be monitored, construct a large-scale, efficient block-based parallel computing module, and realize the automated release and weekly update of monitoring results based on the obtained block-based parallel computing and network map service interface.

[0058] Specifically, this step involves adding geographic coordinate information to the obtained raster results of phenological stages and crop types, and then fusing them with existing vector plot information of the area to be monitored using the Majority Vote method to generate a phenological stage and crop type map with fine spatial units (such as plots and fields). Subsequently, the requests library of Python is used to call the GeoServer REST API to automatically upload the fused result in GeoTIFF format and replace the old layer, realizing real-time updates and visualization of the monitoring results, and presenting the current phenological stage and crop type map on the web page. Figure 4 , Figure 5 ).

[0059] In a preferred embodiment of the present invention, when spatially fusing and assigning attributes to raster phenological products and vector plot data, a mode voting method is used to determine the dominant crop type within each vector polygon, and auxiliary information such as the average confidence level, most recent observation time, and crop rotation type of the pixels within that vector polygon are simultaneously calculated as attribute information; the dominant crop type C dominant Determined by the following formula: ; Where, Count Ci For the vector polygon, the polygon belongs to category C. iThe number of raster cells, and the average confidence level is the arithmetic mean of the confidence levels of all cells within the vector polygon.

[0060] As a preferred embodiment of the present invention, the Majority Vote method used to count the values ​​of all raster cells covered within each vector polygon is as follows: First, the mode of the main category raster is quickly calculated using zonal_stats; then, for the determined main category region, auxiliary indicators such as confidence (mean), lasttime (mode), and rotation (mode) are calculated.

[0061] As a preferred embodiment of the present invention, automated publishing is achieved by calling the REST API of the geographic information server, including the entire process of automatically deleting old version data storage, uploading new version raster result files, configuring and publishing new layers.

[0062] As a preferred embodiment of the present invention, the automatic update of vector layers first obtains server administrator access through HTTP Basic Auth authentication; then, it uploads a new GeoTIFF file to a specified storage path via requests.put(); after the new layer is successfully published and verified to be usable, it deletes the old version's data storage and associated layers via an HTTP DELETE request; and the above steps are encapsulated into an automated function, integrating a mechanism for capturing and handling network anomalies, authentication failures, file operation failures, etc., to ensure that the operation can be rolled back or an alarm can be issued if any stage of the update process fails, thus ensuring the continuous availability of the map service.

[0063] Example This embodiment selects Haiyan County, Jiaxing City, Zhejiang Province, China as the study area, and uses a near real-time monitoring method for mid-season crop phenology and types based on dense time series provided by this invention for analysis. The specific steps are as follows: Step 1) Data Acquisition and Preprocessing: HLS (Harmonized Landsat Sentinel) satellite remote sensing images of the study area were obtained from the NASA EarthData platform, spanning from 2016 to 2024. The data includes blue, green, red, near-infrared, SWIR1, SWIR2, and QA bands, with a spatial resolution of 30m. Each 3660x3660 image covering the study area was regularly cropped into 30x30 grids, with each grid (or block) measuring 122x122 pixels. All blocks were arranged chronologically and stacked into a dense time-series data cube. Simultaneously, the EVI index for each pixel was calculated and added to the time series. The S-CCD algorithm, run in Python, was used to process each 122x122 pixel block. The initial run used historical data (2016-2024), with the following parameters set: tmask_b1='GREEN', tmask_b2='SWIR1' (using green light and the SWIR1 band for cloud detection), lam=20 (Lasso regularization intensity), p_cg=0.9999 (change confidence threshold), and conse=8 (continuous observation number threshold).

[0064] Step 2) Phenological Identification: When new image data arrives, S-CCD update processing is performed. Based on the updated time-series data, the harmonic model parameters are refitted to obtain a coefficient matrix reflecting the latest vegetation growth status. Using the updated harmonic model, the system fits the EVI time-series curve and predicts up to the end of the year. By identifying extreme points within the growth cycle and calculating thresholds for each phenological stage (15%, 90%, etc.), the current stage's position within the growth cycle is determined. Phenological stage determination is based on the comparison between the latest observed values ​​and the predicted curve, dividing the stage into four phases: pre-emergence stage (0%-15%), greening stage (15%-90%), maturity stage (near the peak), senescence stage (EVI value drops from the peak to the 90%-15% range), and pre-emergence stage (EVI value drops to 15%-0%). The stage confidence level is calculated based on the stability of continuous observations. The final output is raster data containing the current phenological stage, confidence level, and multiple cropping index.

[0065] Step 3) Parameter Optimization: To optimize the S-CCD's ability to detect changes in cropping systems, we performed parameter tuning. First, we ran the S-CCD with a set of relaxed parameters (e.g., conse=3, p_cg=0.90) to quickly identify a large number of candidate change points. Then, we manually interpreted the time-series curves of these points and selected 150 sample points: 100 of these were "positive samples" caused by abrupt changes in curve morphology (e.g., amplitude and phase changes) due to actual crop type changes, and 50 were "negative samples" caused by false alarms due to cloud pollution or normal fluctuations within the phenological period. We used a grid search method to search the parameter space (conse: [6, 7, 8], p_cg: [0.99, 0.999, 0.9999]), using the F1 score as the evaluation index, and finally determined the optimal parameter combination to be conse=6, p_cg=0.99. Rerunning the S-CCD with these optimized parameters allowed for more accurate marking of actual cropping type change areas. By improving the breakpoint and outlier output of the S-CCD algorithm, the pixels in the study area are divided into three categories: stable, uncertain, and changing.

[0066] Step 4) Training Sample Construction and Classification Prediction: From the 2024 classification results of stable class pixels, 8300 training samples were obtained using a grid sampling method based on spatial heterogeneity. The sample sizes were set as follows: 2000 for rice, 1000 for wheat, 300 for rapeseed, and 5000 for other categories. The latest observable values ​​(green, red, near-infrared, SWIR1, SWIR2, and EVI) of six bands were extracted for each sample point to form a feature matrix X∈R. 8300×6 A random forest classifier was trained using these features, with the following parameters: n_estimators=150, max_depth=20, min_samples_split=10, min_samples_leaf=5, max_features='sqrt', n_jobs=-1. A batch processing method (batch_size=50,000) was used to input the six band features of all variable class pixels into the model for prediction, obtaining the likelihood probability map of each variable pixel belonging to three crop categories and other categories. Feature importance calculations based on out-of-bag data showed that the near-infrared band and EVI value were the two features that contributed the most to classification.

[0067] Step 5: Probability Update and Classification Based on a Bayesian Framework. This step employs a Bayesian probability update framework, integrating historical knowledge with current observations. First, the posterior probability map P of the crop category output from the previous period (e.g., week T-1) is loaded. priorThis image contains the specific probability values ​​for each pixel belonging to one of four categories: rice, wheat, rapeseed, and other. It also retrieves the pixel stability labels (stable / uncertain / changing) from the previous period. For pixels labeled "stable," their prior probability vectors are modified to a deterministic distribution. For example, if a pixel was previously identified as rice, its probability vector is set to [rice: 1.0, wheat: 0.0, rapeseed: 0.0, other categories: 0.0]. For pixels labeled "uncertain" or "changing," their probability vectors remain unchanged. Then, using the HLS image features from the current period (week T), the pre-trained random forest classification model is input to obtain the probability of each pixel in the entire study area belonging to one of the four categories, forming the current likelihood probability map P. likelihood Next, Bayesian update calculations are performed in parallel on all pixels of the entire image, and the calculated posterior probability map P is then processed. posterior A 3×3 window Gaussian filter is applied to improve spatial smoothness. Finally, each pixel is traversed, and its posterior probability values ​​for the four crop categories are compared. The category with the highest probability is determined as the crop type for that pixel in this period, generating the crop type distribution map for week T. The key outputs of this period include the final classification label map and a crop category probability raster file (each band stores the probability value of one crop). This probability raster file fully records the posterior probability distribution of each pixel belonging to each crop category, which will serve as the prior probability input for Bayesian updates in the next monitoring period (week T+1), thus forming a monitoring closed loop of continuous iterative optimization at the probability level.

[0068] Step 6) Result Fusion and Visualization: The Majority Vote method is used to fuse raster results with vector plot data. Zonal_stats is used to calculate the mode of raster cells within polygons, obtaining the dominant crop type within each vector plot. Average confidence scores are calculated for the main category areas, and auxiliary indicators such as the most recent observation date and crop rotation type are calculated using the mode. An automated script is written to call GeoServer's REST API using the requests library: first, authentication and login are performed; then, a PUT request is used to upload the newly generated fused vector data file; next, a DELETE request is sent to delete the old layer and datastore under the workspace; finally, the new layer is configured and published. Once published successfully, the front-end visualization webpage can retrieve and display the latest monitoring results in real time by calling the WMS service.

[0069] This invention utilizes S-CCD recursive updates and random forest classification to achieve near real-time dynamic updates on a weekly basis. Simultaneously, the monitoring process is highly automated, significantly reducing the workload of manual interpretation and mapping. The related results can provide agricultural management departments with scientific, accurate, and real-time information on crop planting structure and phenological dynamics, significantly improving the efficiency and reliability of agricultural production monitoring and decision support, demonstrating broad application prospects and practical value.

[0070] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences, characterized in that, Specifically as follows: S1: Acquire multispectral remote sensing image data of the area to be monitored from historical to current time, divide the image data into blocks and stack them in chronological order to construct a time series dataset; process the time series dataset using a random continuous change detection algorithm, and recursively update the detection result for newly added image data; S2: Based on the harmonic coefficients in the random continuous change detection results described in S1, construct and predict the characteristic curve for each pixel for the current year, analyze the phenological stage based on the threshold method and the enhanced vegetation index estimate, and generate a raster phenological product that includes multiple cropping index, crop phenological stage and change confidence. S3: Optimize the parameters of the random continuous change detection described in S1, and then classify the change class pixels, uncertain class pixels, and stable class pixels based on the obtained breakpoints and outliers; combine the uncertain class pixels and stable class pixels with the multiple cropping index of S2 to inherit the corresponding historical crop type; S4: Automatically select samples from the stable class pixels described in S3, extract the current time period features and train a random forest classification model to classify crop types for the changing class pixels, and perform Bayesian fusion with the historical crop type posterior probability map to obtain the current crop type map. S5: Spatially fuse and assign attribute values ​​to the raster phenology product described in S2 and the current crop type map described in S4 with the vector plot data of the area to be monitored, construct a large-scale, efficient, block-based parallel computing module, and realize the automated release and weekly update of monitoring results based on the block-based parallel computing and network map service interface.

2. The near-real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In S1, the multispectral bands include visible light, near-infrared band, short-wave infrared band, and quality assessment band.

3. The near-real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In step S1, to enhance the expressive power of temporal features, the enhanced vegetation index (EVI) of the multispectral remote sensing image data is first calculated; then, a single remote sensing image data is divided into multiple spatial blocks and stacked in chronological order to construct a dense temporal dataset for temporal analysis and phenological identification.

4. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In step S1, the random continuous change detection algorithm is performed after cloud masking. Based on the state-space model and Kalman filter recursive update mechanism, harmonic modeling and change detection are performed on the temporal vegetation index of each pixel. The harmonic model obtained through harmonic modeling is used to simulate the annual and semi-annual cycle signals of crop growth. The specific formula for the harmonic coefficients of the harmonic model is as follows: ； ； in, The periodic frequency; For time indexing; The main periodic state variables at time t; These are the principal periodic state variables at time t+1; Main periodic noise; For the main periodic state variables The process noise variance; The auxiliary state variable at time t Let be the auxiliary state variable at time t+1; As auxiliary periodic noise; Auxiliary periodic state variables The process noise variance.

5. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In S3, the parameter optimization method is specifically as follows: A loose threshold is used to perform preliminary breakpoint detection on historical time series data, and positive samples caused by changes in actual crop type and negative samples caused by noise or phenological changes are selected to form a parameter optimization sample set. Using the F1 score as the evaluation index, the optimal parameter combination is obtained by searching in the preset space of random and continuously changing detection parameters through the grid search method, and then used for the final pixel category classification.

6. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In step S4, when automatically selecting samples from stable class pixels, a grid sampling method based on spatial heterogeneity is used, as follows: The area to be monitored is divided into R×C uniform grids. The coefficient of variation (CV) of the EVI values ​​of all pixels in each grid is calculated and classified according to heterogeneity. More sample points are allocated in highly heterogeneous grids to ensure the spatial representativeness of the samples and the balance of class distribution.

7. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In S4, the method for training the random forest classification model to classify crop types for changing class pixels is as follows: First, based on the stable class pixels divided by S3, the latest observable values ​​of the six bands of green light, red light, near infrared, SWIR1, SWIR2 and enhanced vegetation index (EVI) and their corresponding crop category labels are extracted to construct a training sample set; then, the random forest classification model is trained using the training sample set, and the optimal split point is determined according to the principle of maximizing the reduction of Gini index when splitting nodes. For the change class pixels to be classified, extract their six common band feature values, input them into the trained random forest classification model, and determine their final crop category through the majority voting mechanism of multiple decision trees; in view of the large number of change class pixels in large-scale remote sensing data, a batch prediction mechanism is adopted to avoid memory overflow.

8. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In S4, the Bayesian fusion method is as follows: Using the category probability distribution of each pixel in the historical crop type map as the prior probability and the current observation likelihood probability map as the likelihood probability, the posterior probability distribution of all pixels in the region is calculated using the Bayesian update formula; the Bayesian update formula is used to calculate the posterior probability distribution of all pixels in the region for each pixel belonging to category C. i Updated posterior probability P posterior (C i The calculation formula is: ； Among them, P prior (C i Let P be the prior probability. likelihood (C i ) represents the likelihood probability. j Iterate through all crop categories.

9. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequences according to claim 1, characterized in that, In step S5, when spatially fusing and assigning attributes to raster phenological products and vector plot data, the dominant crop type within each vector polygon is determined using the mode voting method, and the average confidence level, most recent observation time, and crop rotation type of the pixels within that vector polygon are simultaneously calculated as attribute information; the dominant crop type C dominant Determined by the following formula: ； Where, Count Ci For the vector polygon, the polygon belongs to category C. i The number of raster cells, and the average confidence level is the arithmetic mean of the confidence levels of all cells within the vector polygon.

10. The near real-time monitoring method for mid-season crop phenology and types based on dense temporal sequence according to claim 1, characterized in that, In S5, automated publishing is achieved by calling the REST API of the geographic information server, including the entire process of automatically deleting old version data storage, uploading new version raster result files, configuring and publishing new layers.

Images