Vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics

Abstract

The application discloses a vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics, and relates to the technical field of remote sensing image recognition, and comprises the following steps: S1, acquiring multi-source remote sensing data of target area covering vegetation in a whole growth period, wherein the multi-source remote sensing data comprises optical remote sensing images and synthetic aperture radar images; S2, pre-processing the multi-source remote sensing data, and constructing a vegetation index time series curve based on the pre-processed multi-source remote sensing data. The vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics can effectively capture the internal differences of different vegetation types in growth rate, duration and curve shape, and realize fine differentiation of vegetation with similar spectra but different growth rhythms, by improving the vegetation classification from a static mode based on key phenological state differences to a fine description of the dynamic process in the whole growth period.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image recognition technology, specifically to a vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics. Background Technology

[0002] Remote sensing technology plays an increasingly crucial role in large-scale vegetation monitoring and fine classification. High temporal and spatial resolution remote sensing satellite data provides strong support for accurately acquiring spatiotemporal distribution information of vegetation. For example, Landsat series satellite data can effectively reflect the spectral differences of different vegetation types and has been widely used in the identification and mapping of forests, crops, etc. Sentinel-2 imagery, due to its higher spatial resolution than Landsat, performs better in characterizing vegetation boundaries and internal heterogeneity. Given that single data sources have their own advantages and disadvantages in terms of observation frequency and spatial detail, combining Landsat with multi-source optical data such as Sentinel-2 has become a reliable solution to improve the accuracy and temporal resolution of vegetation classification. However, the quality of optical imagery is highly susceptible to interference from clouds and shadows, especially in cloudy and rainy areas, and obtaining continuous, cloud-free, high-quality data throughout the entire vegetation growing season remains a significant challenge. Synthetic Aperture Radar (SAR) data, with its all-weather, cloud-penetrating capabilities, provides an effective supplement to overcome the limitations of optical data. Therefore, combining multi-source optical and radar remote sensing data is an important way to improve the accuracy and robustness of vegetation fine classification.

[0003] Phenological characteristics of vegetation are key to distinguishing different vegetation types. Traditional vegetation classification methods often rely on spectral information at a single time point or time-series curves constructed based on vegetation indices. In recent years, nuclear NDVI (kNDVI) has shown greater potential in capturing vegetation phenology due to its higher sensitivity to vegetation biophysical parameters.

[0004] However, most existing studies focus on the seasonal or stage-specific characteristics of vegetation growth. For example, the patent publication "CN119580106A" describes "A Classification Method and System for Double-Crop Rice Based on Multi-Source Remote Sensing Data and Improved Indices," which constructs an improved double-cropping rice identification index by calculating the difference or ratio between the kNDVI and LSWI indices at specific phenological stages, such as transplanting and maturity. While this method is effective, its core lies in capturing the "state differences" of vegetation at two key time points, namely transplanting and maturity. It fails to deeply explore and utilize the dynamic information of growth rate changes throughout the entire growth cycle, especially the "asymmetry" in rate and duration between the growth rise and decline phases. This asymmetry is an external manifestation of the intrinsic biological characteristics of different vegetation types, such as annual herbs and perennial woody plants, and different crop varieties, containing important distinguishing information. However, there is currently a lack of effective methods to utilize this feature for fine-grained vegetation classification. Summary of the Invention

[0005] The purpose of this invention is to provide a vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics to solve the problems mentioned in the background art.

[0006] To address the aforementioned technical problems, this invention provides the following technical solution: a vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics, comprising the following steps:

[0007] S1: Acquire multi-source remote sensing data of the target area covering the entire growth cycle of vegetation, including optical remote sensing images and synthetic aperture radar images.

[0008] S2: Preprocess the multi-source remote sensing data and construct a vegetation index time series curve based on the preprocessed multi-source remote sensing data;

[0009] S3: Based on the vegetation index time series curve, reconstruct the vegetation growth dynamic state trajectory of each pixel, wherein the vegetation growth dynamic state trajectory includes at least the vegetation greenness state, vegetation growth rate, vegetation growth acceleration and vegetation phase state.

[0010] S4: Based on the vegetation growth dynamics trajectory, extract vegetation growth asymmetry features, including rate asymmetry features, duration asymmetry features, and integral asymmetry features based on growth phase division.

[0011] S5: Input the vegetation growth asymmetry features and auxiliary features into the trained Siamese neural network to generate a growth trajectory embedding vector for each pixel; compare and verify the growth trajectory embedding vector with a preset category template library to obtain a preliminary classification result;

[0012] S6: Perform spatiotemporal context optimization on the preliminary classification results to obtain secondary classification results;

[0013] S7: Perform a temporal evolution trajectory consistency check on the secondary classification result to obtain the tertiary classification result;

[0014] S8: Perform spatial ecological logic verification on the three classification results to obtain the final classification result;

[0015] S9: Post-process the final classification results and output a spatial distribution map of vegetation in the target area.

[0016] Optical remote sensing images can use Landsat-8 and Sentinel-2 data, while synthetic aperture radar images use Sentinel-1 data. The vegetation index time series curve is constructed with the kernel normalized vegetation index as the core to ensure that the observation sequence covers the entire growth cycle of vegetation.

[0017] Preferably, the reconstruction step of the vegetation growth dynamics trajectory includes:

[0018] Based on the preprocessed multi-source remote sensing data, the kernel normalized vegetation index time series data for each pixel is calculated;

[0019] A weighted adaptive filtering algorithm is used to smoothly reconstruct the kernel normalized vegetation index time series data to eliminate the influence of noise.

[0020] Based on the reconstructed time-series data, the first derivative of each pixel is calculated as the vegetation growth rate, and the second derivative of each pixel is calculated as the vegetation growth acceleration.

[0021] The phenological period of each pixel is identified based on the dynamic threshold method, including the growth initiation period, peak period and growth end period, and the phenological period is encoded as vegetation phase state;

[0022] The vegetation greenness status, vegetation growth rate, vegetation growth acceleration, and vegetation phase status are combined in the time dimension to form a four-dimensional growth dynamics trajectory for each pixel.

[0023] The weighted adaptive filtering uses the Savitzky-Golay algorithm, and the dynamic threshold method uses the annual amplitude of the kernel normalized vegetation index of 20% as the standard to define the phenological period, ensuring the smoothness of the trajectory data and the accuracy of phenological period identification.

[0024] Preferably, the construction and training steps of the twin neural network include:

[0025] Two temporal convolutional neural networks with identical structures and shared weights are constructed as twin branches. The input of each branch is the growth dynamic state trajectory of the pixel, and the output is an embedding vector of fixed dimension.

[0026] Based on the growth dynamic state trajectories of a large number of unlabeled pixels, a contrastive learning strategy is used to pre-train the twin neural network.

[0027] During the pre-training process, the growth trajectories of the same pixel in different years, the trajectory of the same pixel after adding time-series enhancement, and the growth trajectories of neighboring pixels of the same plot are constructed as positive sample pairs, and the growth trajectories of randomly selected different pixels are constructed as negative sample pairs.

[0028] The network parameters of the Siamese neural network are optimized by maximizing the similarity between positive samples and their embedding vectors and minimizing the similarity between negative samples and their embedding vectors.

[0029] The kernel sizes of the temporal convolutional neural network are 5×1, 3×1, and 3×1, respectively. The pre-training batch size is set to 128, the number of iterations is 100 epochs, and the learning rate is dynamically adjusted through a cosine annealing strategy.

[0030] Preferably, the step of comparing and verifying the growth trajectory embedding vector with a preset category template library to obtain a preliminary classification result includes:

[0031] The growth dynamics state trajectory of each pixel is input into the pre-trained Siamese neural network to generate the corresponding growth trajectory embedding vector.

[0032] Calculate the similarity between the growth trajectory embedding vector and all template embedding vectors in a pre-defined template library for each category, which is constructed based on historical classification results or field survey samples;

[0033] If the average similarity between the current pixel and a certain category is higher than a preset first threshold, then the pixel is directly determined to belong to that category;

[0034] If the average similarity between the current pixel and all categories is lower than the first threshold, the growth trajectory embedding vector of the pixel is merged with the vegetation growth asymmetry feature and auxiliary features, and input into a deep neural network classifier to obtain a preliminary classification probability and a preliminary classification result.

[0035] The template vectors in the category template library are the average of the embedding vectors of typical pixels of the same category. The first threshold is set to 0.8. The deep neural network classifier contains three fully connected layers with 512, 256 nodes and the number of categories, respectively.

[0036] Preferably, the step of performing spatiotemporal context optimization on the preliminary classification result to obtain a secondary classification result includes:

[0037] A spatiotemporal conditional random field model is constructed based on the preliminary classification results. The unary term of the spatiotemporal conditional random field model is provided by the classification probability output by the deep neural network classifier.

[0038] The spatial smoothing term of the spatiotemporal conditional random field model is constructed based on the classification consistency of neighboring pixels, encouraging spatially adjacent pixels to have the same classification result;

[0039] The temporal consistency term of the spatiotemporal conditional random field model is constructed based on the classification results of the same pixel in different years, which encourages the classification results of the same pixel to remain consistent in different years.

[0040] The energy function of the spatiotemporal conditional random field model is optimized using the graph cut algorithm to obtain the optimized quadratic classification result and quadratic classification probability.

[0041] The spatial smoothing term of the spatiotemporal conditional random field adopts a four-neighbor system, the temporal consistency term is based on the classification results of the same pixel over many years, and the graph cut algorithm iterates no more than 50 times to ensure the efficiency and effectiveness of energy function optimization.

[0042] Preferably, the step of performing a temporal evolution trajectory consistency check on the secondary classification result to obtain the tertiary classification result includes:

[0043] For each pixel, extract its original growth dynamics trajectory;

[0044] The original growth dynamic state trajectory is compared with the typical growth trajectory of the corresponding category of the secondary classification result by dynamic time warping, and the temporal evolution distance between the two is calculated.

[0045] If the temporal evolution distance exceeds a preset second threshold, the pixel is determined to be an abnormal pixel;

[0046] For the anomalous pixels, their secondary classification results are used as prior information, and the typical growth trajectory of the corresponding category is used as a reference. The constrained curve fitting method is used to refit their kernel normalized vegetation index time series data.

[0047] Based on the refitted time-series data, vegetation growth asymmetry features were re-extracted and re-input into a deep neural network classifier to obtain updated classification results, which were used as the third classification results.

[0048] The curvature range of dynamic time regularization does not exceed 20% of the total trajectory time length. The second threshold is the mean distance within the class plus 1.5 times the standard deviation. The constraint curve fitting adopts a dual logic curve, and the phenological period deviation is controlled within 15 days.

[0049] Preferably, the step of performing spatial ecological logic verification on the three classification results to obtain the final classification result includes:

[0050] A vegetation niche constraint rule base is constructed. The rule base is based on topographic data, climate data and soil data, and includes the environmental suitability range of different vegetation types.

[0051] Pixels that violate the ecological niche constraint rule library in the three classification results are marked as logically abnormal pixels;

[0052] For the logically abnormal pixels, a weighted spatial interpolation method is used to reassign values. The weights of the weighted spatial interpolation are determined based on the spatial distance and terrain similarity of the non-logically abnormal pixels in the neighborhood.

[0053] The vegetation niche constraint rule library includes suitable ranges for altitude, slope, annual average temperature, annual precipitation, and soil type. Weighted spatial interpolation uses a 5×5 search window, and the weights are calculated by combining spatial distance and terrain similarity.

[0054] Preferably, the post-processing steps for the final classification result include:

[0055] The final classification results are then subjected to weighted voting fusion. The weights of the weighted voting fusion are determined by the confidence of each round of classification results in the local region. The confidence of the local region is calculated based on the average classification probability within the neighborhood.

[0056] The classification results after weighted voting fusion are subjected to mode filtering to eliminate isolated pixels;

[0057] The classification results after mode filtering are processed by morphological closing operation to fill small holes and smooth the boundaries of ground features.

[0058] The neighborhood window for weighted voting fusion is 3×3, the mode filter uses a 3×3 window, and the morphological closing operation uses a 3×3 structuring element. After one iteration, the spatial continuity of the classification results is effectively improved.

[0059] Vegetation classification systems based on multi-source remote sensing and growth asymmetry characteristics include:

[0060] The data acquisition module is used to acquire multi-source remote sensing data of the target area covering the entire growth cycle of vegetation. The multi-source remote sensing data includes optical remote sensing images and synthetic aperture radar images.

[0061] The data preprocessing module is used to preprocess the multi-source remote sensing data and construct a vegetation index time series curve based on the preprocessed multi-source remote sensing data.

[0062] The state reconstruction module is used to reconstruct the vegetation growth dynamic state trajectory of each pixel based on the vegetation index time series curve. The vegetation growth dynamic state trajectory includes at least the vegetation greenness state, vegetation growth rate, vegetation growth acceleration and vegetation phase state.

[0063] The feature extraction module is used to extract vegetation growth asymmetry features based on the vegetation growth dynamic state trajectory. The vegetation growth asymmetry features include rate asymmetry features, duration asymmetry features, and integral asymmetry features based on growth phase division.

[0064] The twin comparison module is used to input the vegetation growth asymmetry features and auxiliary features into a trained twin neural network to generate a growth trajectory embedding vector for each pixel, and to compare and verify the growth trajectory embedding vector with a preset category template library to obtain a preliminary classification result.

[0065] The spatiotemporal optimization module is used to optimize the spatiotemporal context of the preliminary classification results to obtain secondary classification results;

[0066] The timing verification module is used to verify the consistency of the timing evolution trajectory of the secondary classification result to obtain the tertiary classification result;

[0067] The logic verification module is used to perform spatial ecological logic verification on the three classification results to obtain the final classification result;

[0068] The post-processing module is used to post-process the final classification results and output a spatial distribution map of vegetation in the target area.

[0069] Each module achieves automated processing through data flow. The data acquisition module supports data download via FTP or API interface, and the post-processing module outputs vegetation spatial distribution maps in GeoTIFF or Shapefile format, adapting to common remote sensing application scenarios.

[0070] A computer-readable storage medium storing a computer program, which, when executed by a processor, implements the vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics as described above. The storage medium may be a non-volatile memory such as a solid-state drive or a hard disk drive. The program code is divided into functional modules, supporting execution on personal computers, servers, or cloud computing platforms, ensuring cross-device executability of the method.

[0071] This invention provides a vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics. It has the following beneficial effects:

[0072] This vegetation classification method and system, based on multi-source remote sensing and growth asymmetry features, elevates vegetation classification from a static model based on key phenological stage differences to a fine depiction of the dynamic process throughout the entire growth cycle by constructing vegetation growth dynamic state trajectories and extracting growth asymmetry features. It effectively captures the inherent differences in growth rate, duration, and curve morphology among different vegetation types, thus achieving precise differentiation between spectrally similar but differently growing rhythmic vegetation. Simultaneously, a Siamese neural network comparison verification mechanism is introduced, utilizing the growth trajectories of a large number of unlabeled pixels for self-supervised pre-training. This enables the model to automatically learn the essential characteristics of the vegetation growth process. Combined with multi-round cyclic feedback and spatiotemporal context optimization, abnormal pixels in the classification results are corrected step-by-step. The final generated vegetation spatial distribution map shows significant improvements in category boundary clarity, spatial continuity of land features, and conformity with natural geographical laws.

[0073] This vegetation classification method and system, based on multi-source remote sensing and growth asymmetry characteristics, combines optical imagery with synthetic aperture radar (SAR) data. It assesses the quality of optical data through the temporal stability of SAR data and constructs a confidence weight field, effectively reducing the interference of cloud and rain weather on temporal curve reconstruction. This allows the method to maintain stable classification performance even in cloudy and rainy areas. The spatiotemporal conditional random field model considers both spatial neighborhood consistency and multi-year temporal consistency. Spatial ecological logic verification incorporates multi-dimensional environmental constraints such as altitude, slope, climate, and soil, ensuring that the classification results match the geographical patterns of actual vegetation distribution. Weighted voting fusion and morphological post-processing further eliminate isolated pixels and holes. The output results can directly serve practical applications such as vegetation resource surveys, agricultural planting structure monitoring, and ecological red line assessment. The method as a whole exhibits good cross-regional transferability and reproducibility. Attached Figure Description

[0074] Figure 1 This is a flowchart illustrating the vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics of the present invention.

[0075] Figure 2 This is a flowchart illustrating the spatial ecological logic verification process of the vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics of this invention. Detailed Implementation

[0076] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0077] Please see Figure 1 and Figure 2 This invention provides a technical solution: a vegetation classification method and system based on multi-source remote sensing and growth asymmetry characteristics, comprising the following steps:

[0078] S1: Acquire multi-source remote sensing data of the target area covering the entire growth cycle of vegetation. The multi-source remote sensing data includes optical remote sensing images and synthetic aperture radar images.

[0079] S2: Preprocess the multi-source remote sensing data and construct a vegetation index time series curve based on the preprocessed multi-source remote sensing data;

[0080] S3: Based on the vegetation index time series curve, reconstruct the vegetation growth dynamic state trajectory of each pixel. The vegetation growth dynamic state trajectory includes at least the vegetation greenness state, vegetation growth rate, vegetation growth acceleration and vegetation phase state.

[0081] S4: Based on the dynamic state trajectory of vegetation growth, extract the asymmetric features of vegetation growth. The asymmetric features of vegetation growth include rate asymmetric features, duration asymmetric features and integral asymmetric features based on growth phase division.

[0082] S5: Input the vegetation growth asymmetry features and auxiliary features into the trained Siamese neural network to generate the growth trajectory embedding vector for each pixel; compare and verify the growth trajectory embedding vector with the preset category template library to obtain preliminary classification results;

[0083] S6: Perform spatiotemporal context optimization on the preliminary classification results to obtain secondary classification results;

[0084] S7: Perform temporal evolution trajectory consistency verification on the secondary classification results to obtain the tertiary classification results;

[0085] S8: Perform spatial ecological logic verification on the three classification results to obtain the final classification result;

[0086] S9: Post-process the final classification results and output a spatial distribution map of vegetation in the target area.

[0087] It should be further explained that, firstly, the complete growth cycle of the vegetation covering the target area should be obtained, such as multi-source remote sensing data for 12 months of the year or many years. These data include Landsat-8, Landsat-9, Sentinel-2 optical remote sensing images and Sentinel-1 synthetic aperture radar images.

[0088] The acquired remote sensing data is preprocessed as follows: optical images need to undergo cloud cover filtering and declouding, with Landsat data using the QA band for declouding and Sentinel-2 data using the MSK_CLDPRB probabilistic band for declouding, followed by radiometric calibration, atmospheric correction, and geometric fine correction; synthetic aperture radar data needs to undergo thermal noise removal, radiometric calibration, and terrain correction to generate backscattering coefficients in the VV and VH bands; all data are resampled to a uniform spatial resolution, such as 10 meters, and registered with Sentinel-2 imagery as a reference, and finally cropped according to the target area boundary.

[0089] Vegetation index time-series curves were constructed based on preprocessed multi-source remote sensing data: the kernel normalized vegetation index (kNDVI) was calculated scene by scene. The kNDVI uses a radial basis function (RBF) kernel as the kernel function, mapping the reflectance values ​​of the red and near-infrared bands to a high-dimensional space through the RBF kernel. During the mapping process, the RBF mapping rules in the field of vegetation remote sensing were used as a basis, and the mapping parameters were adjusted in combination with the vegetation type characteristics of the target area. After the mapping was completed, the index was calculated according to the conventional calculation logic of the normalized vegetation index, and the kernel normalized vegetation index time-series data for each pixel was obtained. The selection of the RBF kernel and the mapping operation adopted conventional parameter settings that can be directly implemented by those skilled in the art. The calculation formula uses the kernel method to map the reflectance of the red and near-infrared bands to a high-dimensional space and then calculates it to obtain the discrete observation sequence of each pixel in the time dimension. At the same time, the polarization ratio index VH / VV and its coefficient of variation were calculated in combination with synthetic aperture radar data to help evaluate the quality of optical time-series data.

[0090] Next, the vegetation growth dynamics trajectory of each pixel was reconstructed: Weighted adaptive Savitzky-Golay filtering was applied to the kNDVI time-series data. The filter weights were determined by a combination of cloud masking, sensor observation angle, and temporal consistency evaluated by synthetic aperture radar data. The weight fusion of the weighted adaptive Savitzky-Golay filtering used a linear weighting method. The cloud masking factor was quantized and assigned according to the probability of cloud interference in the pixel: pixels without cloud interference were assigned a cloud masking factor of 1; pixels with slight cloud interference were assigned a value of 0.6-0.8; pixels with moderate cloud interference were assigned a value of 0.3-0.5; and pixels with severe cloud interference were assigned a value of 0.5. The pixel value for cloud interference is assigned to 0-0.2; the quantization weight of the sensor observation angle is determined based on the observed zenith angle, with a value of 1 for zenith angles between 0-15°, 0.7-0.9 for 15-30°, 0.4-0.6 for 30-45°, and 0.1-0.3 for zenith angles greater than 45°; the SAR temporal stability is quantified based on the rate of change of the synthetic aperture radar backscattering coefficients adjacent to the corresponding time point of the pixel, with a value of 1 for rate of change between 0-5%, 0.7-0.9 for 5-10%, 0.4-0.6 for 10-20%, and 0.1-0.3 for zenith angles greater than 20%. The allocation coefficients of the three factors are fixed as follows: cloud mask factor accounts for 0.5, sensor observation angle factor accounts for 0.2, and SAR time series stability factor accounts for 0.3. The quantized weight values ​​of each factor are multiplied by the corresponding allocation coefficients and then summed to obtain the filter weight value of each pixel at each time point. This weight value is then applied to the smoothing calculation process of Savitzky-Golay filtering.

[0091] The construction of the confidence weight field combines three factors: cloud masking results from optical remote sensing data, sensor observation angle quality coefficient, and synthetic aperture radar temporal stability assessment value. First, each of the three factors is quantized to obtain the quantized value of each factor for each pixel at each time point. Then, the quantized values ​​of each factor are linearly weighted and fused according to the weight allocation coefficients of the weighted adaptive Savitzky-Golay filter. After fusion, the result is normalized within the range of 0 to 1 to obtain the confidence weight value for each pixel at each time point. The confidence weight values ​​of all pixels in the target area throughout the entire growth cycle are integrated according to the spatial and temporal dimensions to form a three-dimensional confidence weight field with two-dimensional space and one-dimensional time. The quantization index of this weight field is a value between 0 and 1; a higher value indicates better remote sensing data quality for the corresponding pixel at the corresponding time point. This value can be directly applied to filtering and curve fitting processes, providing a weight basis for corresponding operations.

[0092] This yields a smooth reconstructed curve. Based on this curve, the first derivative is calculated as the vegetation growth rate, and the second derivative as the vegetation growth acceleration. A dynamic threshold method or curve curvature analysis is used to automatically identify phenological periods: the start of growth (SOS), peak growth (POP), and end of growth (EOS). The relative positions of these phenological periods are encoded as vegetation phase states. For example, the phase values ​​between SOS and POP are linearly mapped to 0–1, between POP and EOS to 1–2, and the remaining time periods are fixed values. The linear mapping of vegetation phase states is based on the time nodes of vegetation growth: the start of growth is the mapping starting point, with a phase value of 0 at this time point; the peak growth is the mapping midpoint, with a phase value of 1 at this time point; and the end of growth is the mapping ending point, with a phase value of 0 at this time point. The phase value is assigned to 2. For each time point between the growth start period and the peak period, linear interpolation is performed according to the time interval between that time point and the growth start period and the peak period, and the result is in the range of 0 to 1. For each time point between the peak period and the growth end period, linear interpolation is performed according to the time interval between that time point and the peak period and the growth end period, and the result is in the range of 1 to 2. The phase value of all time points before the growth start period is fixed to 0, and the phase value of all time points after the growth end period is fixed to 2, thus completing the vegetation phase state encoding for all time points. The kNDVI value, growth rate, growth acceleration and phase state of each time point are combined into a four-dimensional vector to form the vegetation growth dynamic state trajectory of each pixel in the whole growth cycle.

[0093] Based on this state trajectory, vegetation growth asymmetry characteristics are extracted: First, the growth cycle is divided into five phases according to phenological stages and the zero point of growth acceleration: initiation acceleration phase, rapid growth phase, maturity and stabilization phase, senescence deceleration phase, and apoptosis phase. Then, the average growth rate, duration, and integral value under the curve for each phase are calculated, thereby constructing asymmetric characteristics, including the rate ratio between the rapid growth phase and the senescence deceleration phase (phase rate asymmetry index), the duration ratio between the rapid growth phase and the senescence deceleration phase (phase duration asymmetry index), the ratio of the integral of the rapid growth phase to the integral of the senescence deceleration phase (integral asymmetry index), and the curve skewness calculated through the third moment of the state trajectory.

[0094] The rate quantification index of growth asymmetry characteristics is the average vegetation growth rate within each growth phase, i.e., the arithmetic mean of the growth rate values ​​at all time points within the phase. The calculation range of the integral is the area under the kernel normalized vegetation index (NMR) time series curve corresponding to each growth phase. The area calculation accumulates and sums the NMR values ​​within the phase along the time dimension, and the time interval of the summation is consistent with the observation interval of the remote sensing data. The rate asymmetry characteristic is the ratio of the average growth rate during the rapid growth phase to the average growth rate during the senescence and deceleration phase. The duration asymmetry characteristic is the ratio of the duration of the rapid growth phase to the duration of the senescence and deceleration phase. The integral asymmetry characteristic is the ratio of the area under the curve during the rapid growth phase to the area under the curve during the senescence and deceleration phase. The values ​​obtained by calculating the above ratios are the specific quantification values ​​of each growth asymmetry characteristic.

[0095] Simultaneously, auxiliary features are extracted, including the red, near-infrared, and short-wave infrared reflectance of the median composite images of each phase, texture features extracted based on the first principal component of each phase image of Sentinel-2: entropy, second moment of angle, and inverse difference, terrain features of altitude, slope, and aspect extracted based on DEM, and the mean and coefficient of variation of the synthetic aperture radar backscattering coefficient.

[0096] After merging the above-mentioned vegetation growth asymmetry features with auxiliary features, the input is fed into a pre-trained Siamese neural network to generate a growth trajectory embedding vector for each pixel. The Siamese neural network consists of two temporal convolutional neural networks with identical structures and shared weights. Each branch takes the growth dynamic state trajectory of the pixel as input and outputs a 128-dimensional embedding vector.

[0097] All temporal convolutional layers in the Siamese neural network use the same padding method to ensure that the temporal length of the feature map after convolution is consistent with the input. The pre-training hyperparameters are set with a batch size of 128, and the learning rate is dynamically adjusted using a cosine annealing strategy. The initial learning rate is set to 0.001, and the number of iterations is preferably 100 epochs. Pre-training is stopped when the loss function value no longer decreases for 10 consecutive epochs. In the application of the NT-Xent loss function, the cosine similarity between the embedding vectors of positive and negative sample pairs is calculated first. Then, the similarity value is substituted into the conventional calculation logic of the loss function, and the loss function value is minimized through the backpropagation algorithm to complete the optimization of network parameters. The application of the loss function follows the conventional operating procedures in this field.

[0098] In the pre-training phase, a contrastive learning strategy is adopted. The trajectories of the same pixel in different years, the same pixel with added temporal enhancements such as time axis distortion and amplitude scaling, and the trajectories of neighboring pixels of the same plot are constructed as positive sample pairs. The trajectories of different randomly selected pixels are constructed as negative sample pairs. The network parameters are optimized by maximizing the cosine similarity between the embedding vectors of positive sample pairs and minimizing the similarity between negative sample pairs.

[0099] Preliminary classification results are obtained by comparing the growth trajectory embedding vector with a preset category template library: the category template library consists of the average embedding vectors of typical pixels of each category in historical classification results or field survey samples; for the current pixel, the cosine similarity between its embedding vector and each category template vector is calculated. If the average similarity with a certain category is higher than 0.8, it is directly determined to be of that category. Otherwise, the embedding vector is merged with all the features extracted in step four and input into a deep neural network classifier, which outputs the classification probability of each category. The one with the highest probability is taken as the preliminary classification result, while the probability value is retained for subsequent weighting.

[0100] The initial classification results are optimized using spatiotemporal context to obtain secondary classification results: a spatiotemporal conditional random field model is constructed, whose unary term is provided by the classification probability output by a deep neural network; the spatial smoothing term is constructed based on the classification consistency of four-neighbor or eight-neighbor pixels, and the Potts model is used to encourage neighborhood homonyms; the temporal consistency term is designed for multi-temporal data, such as multi-year imagery, to encourage consistent classification results for the same pixel in different years; the graph cut algorithm is used to minimize the energy function of the conditional random field to obtain the optimized secondary classification results and the corresponding secondary classification probabilities.

[0101] To obtain tertiary classification results, the temporal evolution trajectory consistency of the secondary classification results is verified: For each pixel, its original growth dynamics trajectory is compared with the typical growth trajectory of the corresponding category in the secondary classification results using Dynamic Time Warping (DTW). If the DTW distance exceeds 1.5 times the maximum DTW distance within the category, it is marked as an anomalous pixel. For anomalous pixels, the secondary classification result is used as a priori, and the typical growth trajectory of the corresponding category is used as a reference. Constrained dual-logic curve fitting is used to re-optimize the kNDVI temporal sequence. The constrained dual-logic curve fitting is based on the dual-logic growth law of vegetation growth to construct the basic curve. The curve construction adopts the conventional dual-logic curve construction approach in this field, and the morphological parameters of the curve are adjusted in combination with the full growth cycle characteristics of the vegetation in the target area. The loss function used for fitting is the mean squared error loss function. The mean squared error between the original kernel normalized vegetation index (KNOV) time-series observations and the calculated values ​​of the fitted curve was used as the criterion for evaluating the fitting effect. The phenological period constraints specifically stipulated that the deviation between the growth start, peak, and end periods of the fitted curve and the phenological period of the corresponding typical growth trajectory should not exceed 15 days. The amplitude constraints specifically stipulated that the deviation between the peak value of the KNOV of the fitted curve and the peak value of the typical growth trajectory should not exceed 20% of the peak value of the typical trajectory. During the fitting process, the curve parameters were iteratively adjusted to minimize the value of the mean squared error loss function while satisfying the above constraints, thus completing the refitting of the KNOV time-series data. Based on the optimized time series, asymmetric features were re-extracted and input into the deep neural network again to obtain updated classification results as the third classification result.

[0102] Spatial ecological logic verification is performed on the three classification results to obtain the final classification result: a vegetation niche constraint rule library is pre-constructed, which is based on the digital elevation model and climate raster data of the study area, such as annual mean temperature, precipitation and soil type data, to establish the suitability range of different vegetation types for altitude, slope, aspect, temperature, precipitation and soil type; the three classification results are traversed, and if the classification category of a pixel does not match the niche rule of its location, it is judged as a logically abnormal pixel;

[0103] For logically abnormal pixels, a weighted spatial interpolation method is used to reassign values. The interpolation weights are determined by the inverse of the spatial Euclidean distance between the non-abnormal pixels in the neighborhood and the terrain similarity, based on the altitude difference.

[0104] Finally, the final classification results are post-processed: the local confidence scores of the initial, secondary, and tertiary classifications are defined as the average probability of each round of classification within the neighborhood of the pixel, and used as weights to perform weighted voting fusion of the three rounds of classification results to obtain the fusion result; the fusion result is then subjected to mode filtering with a 3×3 window to eliminate isolated pixels, and then morphological closing operation is performed, i.e., dilation followed by erosion, to fill small holes and smooth the boundaries, finally outputting the vegetation spatial distribution map of the target area.

[0105] The steps for reconstructing the dynamic trajectory of vegetation growth include:

[0106] Based on the preprocessed multi-source remote sensing data, the kernel normalized vegetation index time series data for each pixel is calculated;

[0107] A weighted adaptive filtering algorithm was used to smoothly reconstruct the kernel normalized vegetation index time series data to eliminate the influence of noise.

[0108] Based on the reconstructed time-series data, the first derivative of each pixel is calculated as the vegetation growth rate, and the second derivative of each pixel is calculated as the vegetation growth acceleration.

[0109] The phenological period of each pixel is identified based on the dynamic threshold method, including the growth initiation period, peak period and growth end period, and the phenological period is encoded as vegetation phase state;

[0110] By combining vegetation greenness status, vegetation growth rate, vegetation growth acceleration, and vegetation phase status in the time dimension, a four-dimensional growth dynamics trajectory for each pixel is formed.

[0111] It should be further explained that the reconstruction steps of the vegetation growth dynamics trajectory first involve calculating the kernel normalized vegetation index (kNDVI) of each pixel based on the preprocessed multi-source remote sensing data scene by scene, to obtain discrete time-series observations covering the entire growth cycle.

[0112] Then, a weighted adaptive Savitzky-Golay filtering algorithm is used to smooth and reconstruct the kNDVI time series data. The filter window width is set to 15 days or 30 days according to the data density. The filter weight is determined by three factors: the cloud mask band of the optical image, such as the Landsat QA band or Sentinel-2 MSK_CLDPRB band, which gives the probability of pixel interference by clouds, the quality coefficient corresponding to the zenith angle observed by the sensor, and the auxiliary weight based on the temporal stability assessment of synthetic aperture radar data. The auxiliary weight is obtained by calculating the rate of change of the backscattering coefficient of adjacent synthetic aperture radars before and after the time point. If the rate of change of synthetic aperture radar is lower than the preset threshold and the kNDVI decreases abnormally, the observation point is given a low weight, thereby effectively suppressing the influence of residual cloud noise on the curve morphology. After smooth reconstruction, a continuous kNDVI curve at the daily scale is obtained.

[0113] Next, the first derivative of the reconstructed curve is used to obtain the vegetation growth rate of each pixel, and the second derivative is used to obtain the vegetation growth acceleration.

[0114] Based on the reconstructed curve, a dynamic threshold method is used to identify phenological periods: the moment when the upward segment of the curve reaches 20% of the annual amplitude is taken as the growth start period, the moment when the downward segment of the curve decays to 20% of the annual amplitude is taken as the growth end period, and the moment corresponding to the maximum value of the curve is taken as the peak period; or the curvature analysis method is used to detect the extreme points of the curve curvature to obtain more accurate phenological transition moments.

[0115] The identified phenological periods are encoded as vegetation phase states: the phase value is set to 0 before the growth start period, the phase value is linearly interpolated in the range of 0 to 1 between the growth start period and the peak period, the phase value is linearly interpolated in the range of 1 to 2 between the peak period and the growth end period, and the phase value is fixed at 2 after the growth end period.

[0116] Finally, the kNDVI value, growth rate value, growth acceleration value, and phase state value at each time point are combined into a four-dimensional vector. The four-dimensional vectors at all time points are arranged in chronological order to form the vegetation growth dynamic state trajectory of that pixel.

[0117] The steps for constructing and training a Siamese neural network include:

[0118] Two temporal convolutional neural networks with identical structures and shared weights are constructed as twin branches. The input of each branch is the growth dynamic state trajectory of the pixel, and the output is an embedding vector of fixed dimension.

[0119] Based on the growth dynamic state trajectories of a large number of unlabeled pixels, a contrastive learning strategy is used to pre-train the Siamese neural network.

[0120] During the pre-training process, the growth trajectories of the same pixel in different years, the trajectory of the same pixel after adding time-series enhancement, and the growth trajectories of neighboring pixels of the same plot are constructed as positive sample pairs, and the growth trajectories of randomly selected different pixels are constructed as negative sample pairs.

[0121] The network parameters of the Siamese neural network are optimized by maximizing the similarity between positive samples and their embedding vectors and minimizing the similarity between negative samples and their embedding vectors.

[0122] Further explanation is needed. First, two identical temporal convolutional neural networks are constructed as twin branches. Each branch consists of three stacked temporal convolutional layers. Each convolutional layer is followed by a batch normalization layer and a ReLU activation function layer. Residual connections are introduced between each layer to alleviate the gradient vanishing problem. The first convolutional kernel size is set to 5×1 with a stride of 1 and 64 output channels. The second convolutional kernel size is set to 3×1 with a stride of 1 and 128 output channels. The third convolutional kernel size is set to 3×1 with a stride of 1 and 256 output channels. A global average pooling layer is connected after the third convolutional layer to compress the temporal dimension to 1, resulting in a 256-dimensional feature vector. Then, a fully connected layer maps the feature vector to 128 dimensions, which serves as the embedding vector output by the branch network.

[0123] The two branch networks share the same network structure and weight parameters. One branch receives the vegetation growth dynamics state trajectory of one cell as input, and the other branch receives the state trajectory of another cell as input.

[0124] A contrastive learning strategy was employed to pre-train the Siamese network: a large number of unlabeled pixels were randomly selected from the study area, and training sample pairs were constructed based on their growth dynamics trajectories. Positive sample pairs were constructed in three ways: growth trajectories of the same pixel from different years were used to form positive sample pairs; the original trajectory of the same pixel was combined with its trajectory after time-enhanced processing, including random scaling of the time axis, random scaling of amplitude, and random masking of local time periods; and the growth trajectories of neighboring pixels within the same plot with a spatial distance of less than 30 meters were used to form positive sample pairs. Negative sample pairs were constructed by randomly selecting growth trajectories from different pixels.

[0125] During pre-training, the NT-Xent loss function is used to calculate the contrastive loss. For a training batch containing N samples, a positive sample pair and several negative sample pairs are constructed for each sample. The loss function expression is: the negative logarithm of the similarity of the positive sample pair for each sample is divided by the sum of the similarities of all sample pairs, including both positive and negative samples. The similarity is calculated using cosine similarity, and the temperature coefficient is set to 0.5. The contrastive loss is minimized through backpropagation, and the weight parameters of the Siamese network are updated until the loss converges.

[0126] After pre-training, the weights of the two branches of the Siamese network are fixed. When used later, the state trajectory of a single pixel only needs to be input into any one branch to obtain the 128-dimensional growth trajectory embedding vector of that pixel. This embedding vector can reflect the similarity relationship of the growth trajectories of different pixels in the feature space. The embedding vectors of pixels with the same vegetation type are clustered together in the feature space, while the embedding vectors of pixels with different vegetation types are far apart.

[0127] The steps to obtain preliminary classification results by comparing and verifying the growth trajectory embedding vector with a pre-defined category template library include:

[0128] The growth dynamics state trajectory of each pixel is input into the pre-trained Siamese neural network to generate the corresponding growth trajectory embedding vector.

[0129] Calculate the similarity between the growth trajectory embedding vector and all template embedding vectors in the pre-defined category template library, which is constructed based on historical classification results or field survey samples;

[0130] If the average similarity between the current pixel and a certain category is higher than a preset first threshold, then the pixel is directly determined to belong to that category;

[0131] If the average similarity between the current pixel and all categories is lower than the first threshold, the growth trajectory embedding vector of the pixel is merged with the vegetation growth asymmetry feature and auxiliary features, and then input into the deep neural network classifier to obtain the preliminary classification probability and preliminary classification result.

[0132] It should be further explained that, firstly, a category template library is constructed. This template library is constructed based on pixels with a confidence level exceeding 0.9 in historical classification results or sample points collected from field surveys. For each vegetation type, the growth dynamics trajectory of all typical pixels under that type is selected and input into a pre-trained Siamese neural network to obtain a 128-dimensional embedding vector corresponding to each typical pixel. The average value of the embedding vectors of all typical pixels in the same type is calculated as the template vector for that type. At the same time, the standard deviation between the embedding vectors of all typical pixels in the same type and the template vector is calculated and used for subsequent threshold judgment. The template vectors of all categories and their standard deviations together constitute the category template library.

[0133] For the target pixel to be classified, its growth dynamics state trajectory is first input into the Siamese neural network to generate its 128-dimensional growth trajectory embedding vector.

[0134] Next, the cosine similarity between the embedded vector and each category template vector in the template library is calculated. The cosine similarity calculation formula is the dot product of the two vectors divided by the product of their respective magnitudes. The similarity value ranges from negative one to one. The closer the value is to one, the more consistent the directions of the two vectors are.

[0135] For the similarity results between the target pixel and each category, if the similarity between the target pixel and a certain category is greater than 0.8 and the similarity value exceeds the threshold corresponding to the standard deviation within that category, such as exceeding the average similarity of that category minus twice the standard deviation, then the target pixel is directly determined to belong to that category, and there is no need to proceed to the subsequent deep neural network classification process.

[0136] If the similarity between the target pixel and all categories is below 0.8, or if there is a similarity above 0.8 but the intra-class standard deviation is too large, indicating insufficient template representativeness, then the growth trajectory embedding vector of the target pixel is concatenated with the previously extracted vegetation growth asymmetry features and auxiliary features. The merged feature vector is then fed into a deep neural network classifier as input data. This deep neural network classifier consists of three fully connected layers: the first layer has 512 nodes, the second layer has 256 nodes, and the third layer has nodes equal to the number of vegetation categories to be classified. Each layer is followed by a batch normalization layer and a ReLU activation function, and the last layer is followed by a Softmax function to output the classification probability of each category. The category corresponding to the maximum classification probability is taken as the preliminary classification result of the pixel, and this maximum probability value is saved as the preliminary classification confidence for subsequent weighted fusion steps.

[0137] The steps for optimizing the spatiotemporal context of the preliminary classification results to obtain secondary classification results include:

[0138] A spatiotemporal conditional random field model is constructed based on the preliminary classification results. The unary term of the spatiotemporal conditional random field model is provided by the classification probability output by the deep neural network classifier.

[0139] The spatial smoothing term of the spatiotemporal conditional random field model is constructed based on the classification consistency of neighboring pixels, encouraging spatially adjacent pixels to have the same classification result;

[0140] The temporal consistency term of the spatiotemporal conditional random field model is constructed based on the classification results of the same pixel in different years, which encourages the classification results of the same pixel to remain consistent in different years;

[0141] The energy function of the spatiotemporal conditional random field model is optimized using the graph cut algorithm to obtain the optimized quadratic classification result and the quadratic classification probability.

[0142] It should be further explained that, firstly, a spatiotemporal conditional random field model is constructed. This model is defined on the cell grid of the target area, with each cell corresponding to a random variable whose value is the vegetation category label.

[0143] The unary term of the model is provided by the initial classification probabilities output by the deep neural network classifier. Specifically, for pixel i, the unary potential function is defined as the negative logarithm of the probability that the pixel belongs to each category, i.e. ,in y represents the classification probability output by the deep neural network. i x represents the category label of pixel i. i This represents the feature vector of pixel i.

[0144] The spatial smoothing term of the model is constructed based on the classification consistency of neighboring pixels, adopting the Potts model form. For neighboring pixels i and j in a four-neighbor or eight-neighbor system, the spatial smoothing potential function is defined as:

[0145] , where β s Set the spatial smoothing coefficient to 1.0, [y i ≠y j ] is an indicator function when y i With y j Take 1 if they are not equal, otherwise take 0, Δ ij The spectral difference between pixels i and j is calculated based on the Euclidean distance between their reflectance in the red and near-infrared bands and normalized to the interval of 0 to 1, so that neighboring pixels with similar spectra are subject to stronger penalties when their categories are inconsistent.

[0146] The time consistency term of the model is constructed for regions with multi-year remote sensing data. For pixels i and k in the same geographic location but from different years, the time consistency potential function is defined as:

[0147] , where β t Set the time consistency coefficient to 0.8, [y i ≠y k ] is an indicator function, T ik To measure the temporal stability of pixels at the same location in different years, T is calculated and normalized based on the dynamic temporal regularization distance between the annual kNDVI curves of the two graphs. When the vegetation growth curves are similar in shape between years, T... ik If the value is relatively small, then a stronger penalty will be imposed if the classification results are inconsistent.

[0148] The complete conditional random field energy function is obtained by weighted summation of the univariate term, the spatial smoothing term, and the temporal consistency term:

[0149] The energy function of a spacetime conditional random field (CRF) is constructed based on a linear weighted sum of a univariate term, a spatial smoothing term, and a temporal consistency term. The Δ in the spatial smoothing term... ij Normalize the Euclidean distance values ​​based on the reflectance range of the red and near-infrared bands, converting the Euclidean distance values ​​to values ​​between 0 and 1. The larger the distance, the greater the Δ value. ij The larger the value, the more significant the T value in the time consistency term. ik Normalize the dynamic time warping distance range of the kNDVI curve throughout the year, converting the dynamic time warping distance values ​​into values ​​between 0 and 1. The greater the difference in curve shape, the greater the T... ik The larger the value.

[0150] Weighting coefficient λ s , λ tβ s β t The determination is based on the complexity of vegetation types in the target area and the time series length of remote sensing data. When the vegetation types in the area are complex and the time series length is long, the correlation coefficient λ of the time consistency term should be appropriately increased. t and β t When the spatial heterogeneity of regional vegetation is high, the correlation coefficient λ of the spatial smoothing term should be appropriately increased. s and β s The value of λ is determined by analyzing the vegetation characteristics and remote sensing data characteristics of the target area before cross-regional adaptation. Then, based on the above criteria, the coefficient value is adjusted within a range of 0.2 above and below the example value. The adapted coefficient value needs to be verified through a small-scale sample to determine the final value. s and λ t To balance the contributions of each factor, the weighting coefficients were set to 0.5 and 0.3, respectively. s N represents the set of spatial neighborhood cell pairs. t This represents the set of time neighborhood pixel pairs.

[0151] The α-expansion method in the graph cut algorithm is used to minimize the energy function. The algorithm iterates through all class labels, selects a class α each time, binarizes all non-α labels in the current solution and updates them by minimum cut calculation, until the energy no longer decreases or the maximum number of iterations of 50 is reached.

[0152] After the graph cut optimization is completed, each pixel obtains an updated category label. At the same time, the unary probability value corresponding to the pixel during the optimization process is retained as the secondary classification probability. The optimized category label set constitutes the secondary classification result.

[0153] The steps to obtain the tertiary classification result by performing temporal evolution trajectory consistency verification on the secondary classification result include:

[0154] For each pixel, extract its original growth dynamics trajectory;

[0155] The original growth dynamic state trajectory is compared with the typical growth trajectory of the corresponding category of the secondary classification result by dynamic time warping, and the temporal evolution distance between the two is calculated.

[0156] If the temporal evolution distance exceeds the preset second threshold, the pixel is determined to be an abnormal pixel;

[0157] For anomalous pixels, their secondary classification results are used as prior information, and the typical growth trajectory of the corresponding category is used as a reference. The constraint curve fitting method is used to refit their kernel normalized vegetation index time series data.

[0158] Based on the refitted time-series data, vegetation growth asymmetry features were re-extracted and re-input into a deep neural network classifier to obtain updated classification results, which were used as the third classification results.

[0159] It should be further explained that, firstly, a typical growth trajectory library for each vegetation category is constructed based on the secondary classification results. For each category, pixels with a classification confidence score higher than 0.9 and not located in the category boundary area are selected from the secondary classification results as representative pixels of the category. The original growth dynamic state trajectory of these representative pixels is extracted, and the mean of the trajectory of all representative pixels of the same category at each time point is calculated to form the typical growth trajectory of the category. At the same time, the standard deviation at each time point is calculated as the trajectory variation range.

[0160] For each pixel within the target area, its original vegetation growth dynamics trajectory is extracted, and this trajectory is dynamically time-warped (DTW) aligned with the typical growth trajectory of the corresponding category in the secondary classification result. The DTW algorithm constructs a cumulative distance matrix and searches for the optimal curvature path, calculating the minimum cumulative distance between the two trajectories as the temporal evolution distance. The DTW cumulative distance matrix is ​​constructed based on the four-dimensional feature values ​​of the vegetation growth dynamics trajectory. The four-dimensional feature values ​​of the two comparison trajectories are extracted sequentially in chronological order, and the Manhattan distance of the feature values ​​at corresponding time points is calculated. The Manhattan distances at each time point are then accumulated sequentially to form the initial cumulative distance matrix. The rows and columns of the matrix correspond to the time dimensions of the two comparison trajectories, respectively. The dynamic programming algorithm for the optimal path searches according to the principle of minimum cumulative distance from the starting point to the ending point of the trajectory. During the search process, local bending matching of the time dimension of the trajectory is allowed, and the bending range does not exceed 20% of the total time length of the trajectory. Starting from the starting position of the cumulative distance matrix, the cumulative distance values ​​of adjacent positions are compared in turn, and the position with the smallest distance value is selected as the next node of the path until the ending position of the matrix is ​​reached. The cumulative distance value corresponding to the final optimal path is the temporal evolution distance between the two trajectories. This distance value reflects the degree of deviation between the actual growth trajectory of the pixel under test and the typical growth pattern of the category.

[0161] The dynamic temporal warping distances of all representative pixels within a class to their typical trajectories are statistically analyzed. The mean and standard deviation of the distances within the class are calculated, and the mean plus 1.5 times the standard deviation is used as the anomaly detection threshold. If the temporal evolution distance of the current pixel exceeds this threshold, the pixel is marked as a temporal anomalous pixel.

[0162] For a temporally anomalous pixels, their secondary classification results are used as prior information. The typical growth trajectory of the corresponding category is used as a reference benchmark. The kNDVI time series data of the pixel is refitted using a constrained double-logic curve fitting method. The constraints are set as follows: the variation range of the phenological parameters (growth start period, peak period, growth end period) of the fitted curve relative to the corresponding parameters of the typical trajectory does not exceed ±15 days, and the variation range of the amplitude of the fitted curve does not exceed ±20% of the amplitude of the typical trajectory. Under the constraint conditions, the weighted sum of squared residuals between the fitted curve and the original observations is minimized. The weights are provided by the confidence weight field.

[0163] Based on the refitted kNDVI time-series data, vegetation growth rate, growth acceleration, and phase state were recalculated according to the vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics. The updated growth dynamics trajectory was reconstructed, and vegetation growth asymmetry characteristics were re-extracted according to the vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics.

[0164] The newly extracted feature vector is merged with the original features of the pixel and then input into the deep neural network classifier again to obtain the updated classification probability distribution. The category corresponding to the maximum probability is taken as the three-class classification result of the pixel, and the probability value is saved as the three-class classification confidence.

[0165] The steps to perform spatial ecological logic verification on the three classification results and obtain the final classification result include:

[0166] A vegetation niche constraint rule base was constructed. The rule base was built based on topographic data, climate data and soil data, and included the environmental suitability range of different vegetation types.

[0167] Pixels that violate the niche constraint rule library in the three classification results are marked as logically abnormal pixels.

[0168] For logically abnormal pixels, a weighted spatial interpolation method is used to reassign values. The weights of the weighted spatial interpolation are determined based on the spatial distance and terrain similarity of non-logically abnormal pixels in the neighborhood.

[0169] It should be further explained that, firstly, a vegetation niche constraint rule base was constructed. The establishment of this rule base is based on multi-source geographic environmental data of the study area, including digital elevation model and its derived slope and aspect data, multi-year average temperature and precipitation climate raster data, and soil type distribution map.

[0170] For each vegetation type, the suitability ranges for altitude, slope, aspect, average annual temperature, annual precipitation, and soil type are determined through literature review and field survey sample statistics. For example, the suitable altitude range for a certain vegetation type is 200 meters to 800 meters, the suitable slope range is 0 degrees to 15 degrees, the suitable average annual temperature range is 15 degrees Celsius to 22 degrees Celsius, the suitable annual precipitation range is 1000 mm to 1600 mm, and the soil type is limited to paddy soil or alluvial soil. These suitability ranges are stored in the form of rules to form a rule library containing multi-dimensional constraints.

[0171] Each pixel in the three classification result images is traversed, and the values ​​of altitude, slope, aspect, annual average temperature, annual precipitation, and soil type of the pixel's location are read. These environmental factor values ​​are then compared item by item with the suitability range corresponding to the current classification category of the pixel in the rule base. If all environmental factor values ​​fall within the corresponding suitability range, the pixel is determined to have passed the ecological logic check, and its classification result is retained. If any environmental factor value exceeds the suitability range, the pixel is marked as a logically abnormal pixel. For example, a pixel classified as double-cropping rice located in a high-altitude area at 1200 meters above sea level is outside the suitable altitude range for double-cropping rice and is therefore judged as logically abnormal.

[0172] For all pixels marked as logically abnormal, a weighted spatial interpolation method is used to reassign values. Specifically, a 5×5 pixel search window is set centered on each logically abnormal pixel. Within this window, all non-logically abnormal pixels that pass the ecological logic verification are selected as reference pixels. The spatial Euclidean distance and terrain similarity between each reference pixel and the central abnormal pixel are calculated. Terrain similarity is calculated based on the absolute value of the elevation difference between the two pixels; the smaller the elevation difference, the higher the terrain similarity. The overall weight is obtained by multiplying the spatial distance weight and the terrain similarity weight. The spatial distance weight is calculated using a Gaussian kernel function; the closer the distance, the greater the weight. The terrain similarity weight... Defined as:

[0173] , where Δh is the absolute value of the altitude difference.

[0174] The width of the Gaussian kernel function for weighted spatial interpolation is determined based on the spatial resolution of the remote sensing data of the target area. When the spatial resolution is 10 meters, the kernel width is preferably 3 pixels, and other spatial resolutions can be adjusted proportionally. The spatial distance normalization method uses the interpolation center pixel as the origin and normalizes the spatial Euclidean distance between the neighboring reference pixels and the center pixel according to the spatial resolution of the target area pixels, converting the distance value into a value between 0 and 1. The comprehensive weight is calculated by first obtaining the spatial distance weight and the terrain similarity weight separately, and then multiplying the two weight values ​​directly to obtain the final comprehensive weight of each reference pixel. The spatial distance weight is calculated by combining the Gaussian kernel function with the normalized spatial distance, and the terrain similarity weight is calculated based on the altitude difference. The smaller the altitude difference, the higher the terrain similarity weight.

[0175] The categories of all reference pixels within the window are weighted and voted on. The weight is the comprehensive weight mentioned above. The category with the most votes is the final category of the logically abnormal pixel. If there are fewer than 3 reference pixels in the search window, the window is gradually expanded to 7×7 or 9×9 until the reference pixel count requirement is met.

[0176] After all logically abnormal pixels are reassigned through weighted spatial interpolation, they are combined with the pixels that passed the original verification to form the classification result after spatial ecological logical verification. This result is then output as the final classification result.

[0177] The post-processing steps for the final classification results include:

[0178] The final classification results are weighted and fused by voting. The weight of the weighted voting is determined by the confidence of each round of classification results in the local region. The confidence of the local region is calculated based on the average classification probability in the neighborhood.

[0179] The classification results after weighted voting fusion are subjected to mode filtering to eliminate isolated pixels;

[0180] The classification results after mode filtering are processed by morphological closing operation to fill small holes and smooth the boundaries of ground features.

[0181] It should be further explained that, firstly, a weighted voting fusion operation is performed on the final classification results obtained after spatial ecological logic verification. This operation integrates the preliminary classification results, the secondary classification results, and the tertiary classification results.

[0182] For each pixel, the classification confidence scores for the pixel in the initial classification, secondary classification, and tertiary classification are read from the previous steps. The classification confidence scores are the maximum classification probability output by the deep neural network during the initial classification, the unary probability value retained after optimization by the spatiotemporal conditional random field, and the maximum classification probability output by the deep neural network during reclassification after temporal trajectory consistency verification.

[0183] A 3×3 pixel neighborhood window is constructed centered on the given pixel. The mean initial classification confidence, mean secondary classification confidence, and mean tertiary classification confidence of all pixels within this window are calculated. The neighborhood window size for calculating the confidence of local areas is preferably 3×3 pixels. A weighted average is used for averaging, with the classification probability weight of the central pixel within the window being 2, and the weights of the other 8 neighboring pixels being 1. The calculation involves multiplying the classification probability of each pixel by its corresponding weight to obtain the weighted classification probability of each pixel. Then, the sum of the weighted classification probabilities of all pixels within the window is calculated. Finally, this sum is divided by the sum of the weights within the window. The result is the confidence value of the local area corresponding to the central pixel, ranging from 0 to 1. A higher value indicates a higher confidence level in the classification result of that area. These three means are used as the local confidence weights for the classification results in the corresponding round.

[0184] The optimal neighborhood window size for weighted voting fusion is 3×3 pixels. The center pixel and neighboring pixels within the window are assigned non-equal weights, with the center pixel having a weight of 2 and all other neighboring pixels having a weight of 1. The weighted votes are calculated by first calculating the classification frequency within the neighborhood window for the initial, secondary, and tertiary classifications. The frequency of each classification result is multiplied by the corresponding pixel's weight value to obtain the weighted frequency of each classification result within the window. This weighted frequency is then multiplied by the local region confidence score of the classification result in that round to obtain the final weighted vote count. The local region confidence score is used as the weight coefficient. After summing the weighted vote counts of the three rounds of classification, the category with the highest total weighted vote count is taken as the fused classification result.

[0185] According to the weighted voting rules, for the central pixel, the frequency of all pixels in the neighborhood window in the preliminary classification result is counted, multiplied by the local confidence weight of the preliminary classification, to obtain the weighted votes for the preliminary classification. Similarly, the weighted votes for the secondary and tertiary classifications are calculated. The weighted votes of the three categories are added together, and the category with the highest total weighted votes is taken as the category after weighted voting fusion for that pixel. After traversing all pixels to complete the weighted voting fusion, the fused classification result image is obtained.

[0186] Next, the fused classification result image is subjected to mode filtering. A 3×3 pixel square filtering window is set, and the category of the center pixel of the window is replaced with the category that appears most frequently among all pixels in the window. If multiple categories with the same highest frequency appear, one of them is randomly selected. The entire image is traversed by sliding the window pixel by pixel to eliminate isolated salt-and-pepper noise pixels in the classification result that are inconsistent with the surrounding pixels.

[0187] Finally, morphological closing operations are performed on the classification results image after mode filtering. The closing operation is a dilation followed by erosion operation. The dilation operation uses a 3×3 structuring element to replace the category of each pixel with the category that appears most frequently in the neighborhood, so that the boundaries of similar features expand outward and fill the small holes left by the previous processing inside the features. The erosion operation also uses a 3×3 structuring element to replace the category of each pixel with the category that appears most frequently in the neighborhood, so that the expanded boundaries shrink inward and restore the original general outline of the features, while making the feature boundaries smoother. The number of iterations for the closing operation is set to one to prevent excessive smoothing from causing the loss of feature details.

[0188] After weighted voting fusion, mode filtering, and morphological closing operation, the final post-processed result image is output, which is the final output spatial distribution map of vegetation in the target area.

[0189] Vegetation classification systems based on multi-source remote sensing and growth asymmetry characteristics include:

[0190] The data acquisition module is used to acquire multi-source remote sensing data of the target area covering the entire growth cycle of vegetation. The multi-source remote sensing data includes optical remote sensing images and synthetic aperture radar images.

[0191] The data preprocessing module is used to preprocess multi-source remote sensing data and construct a vegetation index time series curve based on the preprocessed multi-source remote sensing data.

[0192] The state reconstruction module is used to reconstruct the vegetation growth dynamic state trajectory of each pixel based on the vegetation index time series curve. The vegetation growth dynamic state trajectory includes at least the vegetation greenness state, vegetation growth rate, vegetation growth acceleration and vegetation phase state.

[0193] The feature extraction module is used to extract vegetation growth asymmetry features based on the vegetation growth dynamic state trajectory. The vegetation growth asymmetry features include rate asymmetry features, duration asymmetry features, and integral asymmetry features based on growth phase division.

[0194] The twin comparison module is used to input the asymmetric features of vegetation growth and auxiliary features into the trained twin neural network, generate the growth trajectory embedding vector of each pixel, and compare and verify the growth trajectory embedding vector with the preset category template library to obtain the preliminary classification result.

[0195] The spatiotemporal optimization module is used to optimize the spatiotemporal context of the preliminary classification results to obtain secondary classification results;

[0196] The temporal verification module is used to verify the consistency of the temporal evolution trajectory of the secondary classification results and obtain the tertiary classification results.

[0197] The logic verification module is used to perform spatial ecological logic verification on the three classification results to obtain the final classification result.

[0198] The post-processing module is used to post-process the final classification results and output a spatial distribution map of vegetation in the target area.

[0199] It should be further explained that the system consists of nine functional modules connected in sequence: data acquisition module, data preprocessing module, state reconstruction module, feature extraction module, twin comparison module, spatiotemporal optimization module, time series verification module, logic verification module, and post-processing module.

[0200] The data acquisition module is used to access multi-source remote sensing data covering the entire growth cycle of vegetation in the target area, including Landsat-8 optical imagery, Landsat-9 optical imagery, Sentinel-2 optical imagery, and Sentinel-1 synthetic aperture radar imagery. This module automatically downloads raw remote sensing data of a specified time and spatial range from data distribution platforms such as USGS and ESA via FTP or API interfaces, and transmits the raw data packets to the data preprocessing module.

[0201] After receiving the raw data, the data preprocessing module performs cloud filtering and declouding. For Landsat data, cloud pixels are identified and masked using the QA band. For Sentinel-2 data, strict declouding is performed using the MSK_CLDPRB probability band with a threshold of 0. For Sentinel-1 synthetic aperture radar data, thermal noise removal and boundary removal are performed. Then, all data are cropped according to the target area vector boundary and image registration is performed based on Sentinel-2 imagery. Finally, all registered data are resampled to a spatial resolution of 10 meters to generate a preprocessed standardized data stack.

[0202] The state reconstruction module receives the preprocessed data stack and first calculates the kernel normalized vegetation index (kNDVI) pixel by pixel to form a time-series curve. A weighted adaptive Savitzky-Golay filter is used to smooth and reconstruct the time-series curve. The filter weights are determined by a combination of cloud mask results, sensor observation angles, and auxiliary weights calculated based on the time-series stability of synthetic aperture radar data. The first derivative is calculated as the growth rate and the second derivative as the growth acceleration based on the reconstructed curve. The dynamic threshold method is used to identify the growth initiation period, peak period, and growth termination period and encode them as phase states. The kNDVI values, growth rate values, growth acceleration values, and phase state values ​​at each time point are combined into a four-dimensional vector and arranged in chronological order to form the vegetation growth dynamic state trajectory of each pixel.

[0203] The feature extraction module receives the state trajectory generated by the state reconstruction module. Based on the phenological period and the zero acceleration point, the growth cycle is divided into five phases: the initiation acceleration phase, the rapid growth phase, the maturity and stabilization phase, the aging and deceleration phase, and the apoptosis phase. The average growth rate, duration, and integral value under the curve of each phase are calculated. Then, the phase rate asymmetry index, the phase duration asymmetry index, and the integral asymmetry index are constructed. At the same time, the red, near-infrared, and short-wave infrared reflectance of the median image of each phase is extracted as spectral features. The entropy, second moment of angle, and inverse difference of the first principal component of each phase image of Sentinel-2 are extracted as texture features. The altitude, slope, and aspect of the DEM data are extracted as terrain features. The mean and coefficient of variation of the synthetic aperture radar backscattering coefficient are extracted as polarization features.

[0204] The Siamese contrast module incorporates a pre-trained Siamese neural network, which consists of two identical temporal convolutional neural networks with shared weights. During the pre-training phase, a contrastive learning strategy is employed, using the growth trajectories of unlabeled pixels to construct positive and negative sample pairs for training. The module receives all features output from the feature extraction module and the state trajectory output from the state reconstruction module. First, the state trajectory is input into the Siamese neural network to generate a 128-dimensional growth trajectory embedding vector. Then, this embedding vector is compared with a pre-defined category template library using cosine similarity. If the similarity to a certain category is higher than 0.8 and exceeds the intra-class standard deviation threshold, the category is directly output as the preliminary classification result. Otherwise, the embedding vector is concatenated with all features output from the feature extraction module and fed into a deep neural network classifier. The classification probability is output, and the category corresponding to the maximum value is taken as the preliminary classification result. The classification probability value is also saved.

[0205] The spatiotemporal optimization module receives the preliminary classification results and classification probabilities output by the twin comparison module, constructs a spatiotemporal conditional random field model, the unary term is composed of the negative logarithm of the classification probability, the spatial smoothing term is based on the classification consistency of four neighboring pixels and weighted by spectral difference, the temporal consistency term is constructed for multi-year data and weighted by temporal stability, and the α-expansion method of the graph cut algorithm is used to minimize the energy function to obtain the optimized secondary classification results and the corresponding secondary classification probabilities.

[0206] The temporal verification module receives the secondary classification results output by the spatiotemporal optimization module. First, it constructs a typical growth trajectory library for each category based on the high-confidence pixels in the secondary classification results. Then, it performs dynamic time warping comparison between the original state trajectory of each pixel and the typical trajectory of the corresponding category, calculates the temporal evolution distance and compares it with the intra-class distance threshold. If the distance exceeds the threshold, it is marked as an abnormal pixel. For abnormal pixels, the kNDVI temporal sequence is reconstructed using constrained double logistic curve fitting with the secondary classification results as priors. Based on the reconstructed temporal sequence, features are re-extracted and input into the deep neural network again to obtain the updated tertiary classification results and tertiary classification probabilities.

[0207] The logic verification module receives the three classification results output by the time-series verification module. Internally, it stores a vegetation niche constraint rule library constructed based on the digital elevation model of the study area, climate raster data, and soil type data. It reads the environmental factor values ​​of each pixel and compares them with the suitability range of the current classification category. If any factor is out of range, it is marked as a logically abnormal pixel. For logically abnormal pixels, a 5×5 window is used to search for reference pixels. The module calculates a comprehensive weight based on spatial distance and terrain similarity, performs weighted interpolation, and reassigns the values ​​to obtain the final classification result that passes the ecological logic verification.

[0208] The post-processing module receives the final classification result output by the logic verification module. First, it reads the confidence level of the pixel in the preliminary classification, secondary classification, and tertiary classification. It calculates the average confidence level of each round in the neighborhood as the weight, performs weighted voting fusion on the three rounds of classification results, then performs 3×3 window mode filtering on the fusion result to eliminate isolated pixels, and finally performs morphological closing operation to fill small holes and smooth the boundaries, outputting the processed target area vegetation spatial distribution map.

[0209] The above nine modules are connected in sequence, with the output data of the previous module serving as the input data for the next module. Together, they complete the fully automated processing from the input of raw remote sensing data to the output of the final vegetation distribution map.

[0210] A computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements a vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics, as described above. It should be further noted that the storage medium is in the form of a non-volatile memory, including but not limited to any one of a solid-state drive, hard disk drive, optical disk, read-only memory, or universal serial bus flash memory drive, and stores computer-executable program code thereon.

[0211] The program code is written into multiple functional modules, including a data acquisition module, a data preprocessing module, a state reconstruction module, a feature extraction module, a twin comparison module, a spatiotemporal optimization module, a timing verification module, a logic verification module, and a post-processing module. The software implementation of each module is completely consistent with the logical function of the corresponding module in the system.

[0212] When the processor loads and executes the program code, it first calls the data acquisition module to read multi-source remote sensing data of the target area covering the entire growth cycle of vegetation from external storage or network interface, including Landsat-8 optical imagery, Landsat-9 optical imagery, Sentinel-2 optical imagery, and Sentinel-1 synthetic aperture radar imagery.

[0213] Next, the data preprocessing module is invoked to perform cloud cover filtering, cloud removal, radiometric calibration, atmospheric correction, geometric fine correction, boundary clipping, image registration, and resampling operations on the read raw remote sensing data, generating a standardized data stack and storing it in memory or a temporary file.

[0214] Then, the state reconstruction module is called to calculate the kernel normalized vegetation index time series curve pixel by pixel based on the preprocessed data. Weighted adaptive filtering is used for smooth reconstruction. The first and second derivatives are calculated as growth rate and growth acceleration, respectively. Phenological periods are identified and encoded as phase states, and combined to form a four-dimensional vegetation growth dynamic state trajectory.

[0215] Then, the feature extraction module is called to divide the growth phase based on the state trajectory and extract the phase rate asymmetry index, phase duration asymmetry index, integral asymmetry index, as well as spectral features, texture features, terrain features and polarization features.

[0216] Then, the twin comparison module is invoked. This module contains pre-trained twin neural network weight parameters. The module inputs the state trajectory into the twin network to generate a 128-dimensional embedding vector, which is then compared and verified with a preset category template library before outputting preliminary classification results and classification confidence.

[0217] Next, the spatiotemporal optimization module is called to construct a spatiotemporal conditional random field model and optimize it using the graph cut algorithm, outputting the secondary classification result and the secondary classification probability.

[0218] Then, the time-series verification module is called to perform dynamic time warping comparison between the original state trajectory and the typical trajectory of the corresponding category in the secondary classification. Constraint curve refitting and reclassification are performed on the time-series abnormal pixels, and the tertiary classification results and tertiary classification probabilities are output.

[0219] Then, the logic verification module is called to load the vegetation niche constraint rule library stored in the medium, compare the suitability range of environmental factors and classification categories pixel by pixel, perform weighted spatial interpolation on logically abnormal pixels, and output the final classification result.

[0220] Finally, the post-processing module is called to perform weighted voting fusion, mode filtering, and morphological closing operations on the final classification results, and outputs a spatial distribution map of vegetation in the target area.

[0221] All intermediate calculation results involved in the execution of the above program code are temporarily stored in memory or cache, and the final generated vegetation spatial distribution map is persistently stored in GeoTIFF or Shapefile format on the computer-readable storage medium or the specified output path.

[0222] This storage medium can be installed in general-purpose personal computers, servers, embedded systems, or cloud computing platforms, and implements a vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics through processor calls.

[0223] This system elevates vegetation classification from a static model based on key phenological stage differences to a fine depiction of the dynamic process throughout the entire growth cycle by constructing vegetation growth dynamic state trajectories and extracting growth asymmetry features. It effectively captures the inherent differences in growth rate, duration, and curve morphology among different vegetation types, thus achieving precise differentiation between spectrally similar but differently growing rhythmic vegetation. Simultaneously, a Siamese neural network comparison and verification mechanism is introduced, utilizing the growth trajectories of a large number of unlabeled pixels for self-supervised pre-training. This allows the model to automatically learn the essential characteristics of the vegetation growth process. Combined with multi-round cyclic feedback and spatiotemporal context optimization, abnormal pixels in the classification results are corrected step-by-step. The final generated vegetation spatial distribution map shows significant improvements in category boundary clarity, spatial continuity of land features, and conformity with natural geographical laws.

[0224] This system combines optical imagery with synthetic aperture radar (SAR) data. It assesses the quality of optical data through the temporal stability of SAR data and constructs a confidence weight field, effectively reducing the interference of cloud and rain weather on temporal curve reconstruction. This allows the method to maintain stable classification performance even in cloudy and rainy areas. The spatiotemporal conditional random field (SFR) model considers both spatial neighborhood consistency and multi-year temporal consistency. Spatial ecological logic verification incorporates multi-dimensional environmental constraints such as altitude, slope, climate, and soil, ensuring that the classification results match the geographical patterns of actual vegetation distribution. Weighted voting fusion and morphological post-processing further eliminate isolated pixels and holes. The output results can directly serve practical applications such as vegetation resource surveys, agricultural planting structure monitoring, and ecological red line assessment. The method as a whole possesses good cross-regional transferability and reproducibility.

[0225] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0226] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics, characterized in that, Includes the following steps: S1: Acquire multi-source remote sensing data of the target area covering the entire growth cycle of vegetation, including optical remote sensing images and synthetic aperture radar images. S2: Preprocess the multi-source remote sensing data and construct a vegetation index time series curve based on the preprocessed multi-source remote sensing data; S3: Based on the vegetation index time series curve, reconstruct the vegetation growth dynamic state trajectory of each pixel, wherein the vegetation growth dynamic state trajectory includes at least the vegetation greenness state, vegetation growth rate, vegetation growth acceleration and vegetation phase state. S4: Based on the vegetation growth dynamics trajectory, extract vegetation growth asymmetry features, including rate asymmetry features, duration asymmetry features, and integral asymmetry features based on growth phase division. S5: Input the vegetation growth asymmetry features and auxiliary features into the trained Siamese neural network to generate a growth trajectory embedding vector for each pixel; compare and verify the growth trajectory embedding vector with a preset category template library to obtain a preliminary classification result; S6: Perform spatiotemporal context optimization on the preliminary classification results to obtain secondary classification results; S7: Perform a temporal evolution trajectory consistency check on the secondary classification result to obtain the tertiary classification result; S8: Perform spatial ecological logic verification on the three classification results to obtain the final classification result; S9: Post-process the final classification results and output a spatial distribution map of vegetation in the target area.

2. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 1, characterized in that: The steps for reconstructing the vegetation growth dynamics trajectory include: Based on the preprocessed multi-source remote sensing data, the kernel normalized vegetation index time series data for each pixel is calculated; A weighted adaptive filtering algorithm is used to smoothly reconstruct the kernel normalized vegetation index time series data to eliminate the influence of noise. Based on the reconstructed time-series data, the first derivative of each pixel is calculated as the vegetation growth rate, and the second derivative of each pixel is calculated as the vegetation growth acceleration. The phenological period of each pixel is identified based on the dynamic threshold method, including the growth initiation period, peak period and growth end period, and the phenological period is encoded as vegetation phase state; The vegetation greenness status, vegetation growth rate, vegetation growth acceleration, and vegetation phase status are combined in the time dimension to form a four-dimensional growth dynamics trajectory for each pixel.

3. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 2, characterized in that: The construction and training steps of the twin neural network include: Two temporal convolutional neural networks with identical structures and shared weights are constructed as twin branches. The input of each branch is the growth dynamic state trajectory of the pixel, and the output is an embedding vector of fixed dimension. Based on the growth dynamic state trajectories of a large number of unlabeled pixels, a contrastive learning strategy is used to pre-train the twin neural network. During the pre-training process, the growth trajectories of the same pixel in different years, the trajectory of the same pixel after adding time-series enhancement, and the growth trajectories of neighboring pixels of the same plot are constructed as positive sample pairs, and the growth trajectories of randomly selected different pixels are constructed as negative sample pairs. The network parameters of the Siamese neural network are optimized by maximizing the similarity between positive samples and their embedding vectors and minimizing the similarity between negative samples and their embedding vectors.

4. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 3, characterized in that: The steps for comparing and verifying the growth trajectory embedding vector with a pre-defined category template library to obtain preliminary classification results include: The growth dynamics state trajectory of each pixel is input into the pre-trained Siamese neural network to generate the corresponding growth trajectory embedding vector. Calculate the similarity between the growth trajectory embedding vector and all template embedding vectors in a pre-defined template library for each category, which is constructed based on historical classification results or field survey samples; If the average similarity between the current pixel and a certain category is higher than a preset first threshold, then the pixel is directly determined to belong to that category; If the average similarity between the current pixel and all categories is lower than the first threshold, the growth trajectory embedding vector of the pixel is merged with the vegetation growth asymmetry feature and auxiliary features, and input into a deep neural network classifier to obtain a preliminary classification probability and a preliminary classification result.

5. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 4, characterized in that: The steps for optimizing the spatiotemporal context of the preliminary classification results to obtain secondary classification results include: A spatiotemporal conditional random field model is constructed based on the preliminary classification results. The unary term of the spatiotemporal conditional random field model is provided by the classification probability output by the deep neural network classifier. The spatial smoothing term of the spatiotemporal conditional random field model is constructed based on the classification consistency of neighboring pixels, encouraging spatially adjacent pixels to have the same classification result; The temporal consistency term of the spatiotemporal conditional random field model is constructed based on the classification results of the same pixel in different years, which encourages the classification results of the same pixel to remain consistent in different years. The energy function of the spatiotemporal conditional random field model is optimized using the graph cut algorithm to obtain the optimized quadratic classification result and quadratic classification probability.

6. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 5, characterized in that: The steps for performing temporal evolution trajectory consistency verification on the secondary classification results to obtain the tertiary classification results include: For each pixel, extract its original growth dynamics trajectory; The original growth dynamic state trajectory is compared with the typical growth trajectory of the corresponding category of the secondary classification result by dynamic time warping, and the temporal evolution distance between the two is calculated. If the temporal evolution distance exceeds a preset second threshold, the pixel is determined to be an abnormal pixel; For the anomalous pixels, their secondary classification results are used as prior information, and the typical growth trajectory of the corresponding category is used as a reference. The constrained curve fitting method is used to refit their kernel normalized vegetation index time series data. Based on the refitted time-series data, vegetation growth asymmetry features were re-extracted and re-input into a deep neural network classifier to obtain updated classification results, which were used as the third classification results.

7. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 6, characterized in that: The steps for performing spatial ecological logic verification on the three classification results to obtain the final classification result include: A vegetation niche constraint rule base is constructed. The rule base is based on topographic data, climate data and soil data, and includes the environmental suitability range of different vegetation types. Pixels that violate the ecological niche constraint rule library in the three classification results are marked as logically abnormal pixels; For the logically abnormal pixels, a weighted spatial interpolation method is used to reassign values. The weights of the weighted spatial interpolation are determined based on the spatial distance and terrain similarity of the non-logically abnormal pixels in the neighborhood.

8. The vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics according to claim 7, characterized in that: The post-processing steps for the final classification results include: The final classification results are then subjected to weighted voting fusion. The weights of the weighted voting fusion are determined by the confidence of each round of classification results in the local region. The confidence of the local region is calculated based on the average classification probability within the neighborhood. The classification results after weighted voting fusion are subjected to mode filtering to eliminate isolated pixels; The classification results after mode filtering are processed by morphological closing operation to fill small holes and smooth the boundaries of ground features.

9. A vegetation classification system based on multi-source remote sensing and growth asymmetry characteristics, characterized in that, include: The data acquisition module is used to acquire multi-source remote sensing data of the target area covering the entire growth cycle of vegetation. The multi-source remote sensing data includes optical remote sensing images and synthetic aperture radar images. The data preprocessing module is used to preprocess the multi-source remote sensing data and construct a vegetation index time series curve based on the preprocessed multi-source remote sensing data. The state reconstruction module is used to reconstruct the vegetation growth dynamic state trajectory of each pixel based on the vegetation index time series curve. The vegetation growth dynamic state trajectory includes at least the vegetation greenness state, vegetation growth rate, vegetation growth acceleration and vegetation phase state. The feature extraction module is used to extract vegetation growth asymmetry features based on the vegetation growth dynamic state trajectory. The vegetation growth asymmetry features include rate asymmetry features, duration asymmetry features, and integral asymmetry features based on growth phase division. The twin comparison module is used to input the vegetation growth asymmetry features and auxiliary features into a trained twin neural network to generate a growth trajectory embedding vector for each pixel, and to compare and verify the growth trajectory embedding vector with a preset category template library to obtain a preliminary classification result. The spatiotemporal optimization module is used to optimize the spatiotemporal context of the preliminary classification results to obtain secondary classification results; The timing verification module is used to verify the consistency of the timing evolution trajectory of the secondary classification result to obtain the tertiary classification result; The logic verification module is used to perform spatial ecological logic verification on the three classification results to obtain the final classification result; The post-processing module is used to post-process the final classification results and output a spatial distribution map of vegetation in the target area.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the vegetation classification method based on multi-source remote sensing and growth asymmetry characteristics as described in any one of claims 1 to 8.

Images