Wind-sand land planting bearing capacity evaluation method and system based on multispectral remote sensing

Abstract

The application provides a wind-sand land planting carrying capacity evaluation method and system based on multispectral remote sensing, and particularly relates to the fields of remote sensing image processing and agricultural information technology. The method acquires a time sequence image stack composed of multi-temporal multispectral remote sensing images of a target wind-sand land evaluation area in a preset time period, the adjacent time phases of which are less than a preset threshold. A dynamic interference recognition model is used to extract an instantaneous dynamic sand particle mask and calculate a pixel-level ground surface dynamic interference index. An evaluation time phase is determined, a background signal feature map is generated by separating the background signal from the signal using a signal decoupling model, and a preliminary evaluation result is obtained by inputting the anti-dynamic interference evaluation model. A clean time phase is determined according to the instantaneous dynamic sand particle mask, and a historical carrying capacity evaluation result is obtained. The preliminary result is corrected in the time domain according to the historical result, and a final planting carrying capacity distribution map is output, so that the essential planting carrying capacity of the land can be stably and accurately evaluated, and the reliability and actual value of the evaluation result are improved.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing image processing and agricultural information technology, and more specifically, to a method and system for assessing the planting carrying capacity of sandy land based on multispectral remote sensing. Background Technology

[0002] With the deep integration of multispectral remote sensing and agricultural information technology, large-scale and rapid assessment of land planting carrying capacity using remote sensing methods has become an important direction for modern agricultural management. Especially in ecologically fragile arid and semi-arid regions, accurate assessment of the planting potential of wind-blown sandy land is crucial for combating desertification and promoting sustainable agricultural development. However, due to the characteristics of wind and sand, the surface sand cover of wind-blown sandy land often changes rapidly and irregularly within days or even hours. This high dynamism causes strong transient interference to key inversion results based on remote sensing spectra, such as vegetation indices and soil parameters, severely affecting the stability and accuracy of subsequent carrying capacity assessment models.

[0003] To improve the reliability of remote sensing assessments of wind-blown sandy land, traditional methods mainly focus on two aspects: first, pursuing higher spatiotemporal resolution to try to "see" dynamic details; and second, constructing more complex statistical or intelligent models to "fit" stable patterns from mixed signals. For example, high-frequency imagery is used for time-series analysis to extract trends, or deep learning models are used to directly learn the complex mapping between spectral features and carrying capacity. However, these methods fail to effectively distinguish between "stable components representing the inherent properties of land" and "interference components caused by transient sand cover" in the spectral signal at the underlying model level, instead mixing all variations. This makes the assessment results susceptible to being dominated by short-term surface conditions, misjudging "fertile land temporarily covered by sand" as having low carrying capacity and "barren land temporarily free of sand cover" as having high carrying capacity. Essentially, this misinterprets dynamic noise signals as essential land features, which is the root cause of model distortion.

[0004] In summary, how to effectively overcome the spectral signal aliasing caused by rapid dynamic sand cover in the assessment of planting carrying capacity of aeolian sandy land based on multispectral remote sensing, and thus avoid misjudging instantaneous interference as an inherent property of the land, is an urgent problem to be solved. Summary of the Invention

[0005] The main objective of this invention is to provide a method and system for assessing the planting carrying capacity of aeolian sandy land based on multispectral remote sensing. This method aims to at least solve the technical problem of spectral signal aliasing caused by rapid dynamic sand cover in the assessment of planting carrying capacity of aeolian sandy land based on multispectral remote sensing, thereby avoiding the misjudgment of transient interference as an inherent property of the land. This enables a stable and accurate assessment of the land's intrinsic planting carrying capacity, effectively eliminates assessment noise caused by transient aeolian activity, and significantly improves the reliability and practical guiding value of the assessment results in the time domain.

[0006] To achieve the above objectives, the present invention provides a method and system for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing.

[0007] In a first aspect, the present invention provides a method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing, the method comprising:

[0008] Acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a preset time period, with the time interval between adjacent time phases being less than a preset threshold. The multi-temporal multispectral remote sensing images constitute a time-series image stack, which is composed of multiple pixels.

[0009] Based on the spectral temporal variation characteristics of each pixel in the temporal image stack, a dynamic interference identification model for identifying instantaneous wind and sand cover based on the spectral temporal variation characteristics is used to extract the instantaneous dynamic sand mask of each time phase of the target wind and sand assessment area, and the pixel-level surface dynamic interference index is calculated.

[0010] An evaluation phase is determined from the time-series image stack; based on the image of the evaluation phase and the instantaneous dynamic sand mask corresponding to the image of the evaluation phase, the background signal characterizing the land stability attribute is separated from the original spectral signal of each pixel through a signal decoupling model, and a background signal feature map is generated.

[0011] The background signal feature map and the surface dynamic interference index are input into the pre-trained anti-dynamic interference assessment model to obtain the preliminary planting carrying capacity assessment results of the target sandy land assessment area at the assessment time phase.

[0012] Based on the instantaneous dynamic sand mask, at least one clean phase is determined from the time-series image stack; based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference evaluation model, the corresponding historical carrying capacity evaluation results are obtained;

[0013] The historical carrying capacity assessment results are used to perform time-domain consistency correction on the preliminary planting carrying capacity assessment results, and the final planting carrying capacity distribution map of the target sandy land assessment area is output.

[0014] Specifically, the acquisition of multi-temporal multispectral remote sensing images of the target sandstorm assessment area within a preset time period, with adjacent time intervals less than a preset threshold, and the formation of a time-series image stack from these multi-temporal multispectral remote sensing images, includes:

[0015] Acquire multi-temporal, multispectral remote sensing images of the target aeolian sandy land assessment area within a complete windy season cycle, with the time interval between adjacent temporal phases being less than or equal to 7 days;

[0016] The multi-temporal multispectral remote sensing images are subjected to radiometric calibration and atmospheric correction to obtain preprocessed multi-temporal multispectral remote sensing images.

[0017] The preprocessed multi-temporal multispectral remote sensing images are stacked in chronological order to form the time-series image stack.

[0018] Specifically, based on the spectral temporal variation characteristics of each pixel in the temporal image stack, a dynamic interference identification model for identifying instantaneous wind and sand cover based on the spectral temporal variation characteristics is used to extract the instantaneous dynamic sand mask of each temporal phase of the target wind and sand assessment area, and to calculate the pixel-level surface dynamic interference index, including:

[0019] For each pixel in the time-series image stack, the normalized vegetation index and normalized sand index are calculated for each time phase to form the spectral temporal variation characteristics of the pixel.

[0020] Based on the spectral temporal variation characteristics, the abrupt change amplitude of the normalized sand grain index of each pixel between adjacent time phases is detected, and pixels with abrupt change amplitude exceeding a first threshold are marked as instantaneous dynamic sand grain pixels, thereby generating the instantaneous dynamic sand grain mask for each time phase.

[0021] For each pixel, the total frequency of the pixel marked as the instantaneous dynamic sand grain pixel within the corresponding time period of the time-series image stack is counted, and the total frequency is normalized and used as the surface dynamic interference index of the pixel.

[0022] Specifically, the step involves determining an evaluation phase from the temporal image stack; based on the image of the evaluation phase and the corresponding instantaneous dynamic sand mask, using a signal decoupling model, separating the background signal characterizing land stability attributes from the original spectral signal of each pixel, and generating a background signal feature map, including:

[0023] The latest image is selected from the time-series image stack as the evaluation time phase;

[0024] The original multispectral reflectance values ​​of each pixel in the image at the evaluation time phase constitute the original spectral signal;

[0025] Using a signal decoupling model based on sparse representation, and based on the instantaneous dynamic sand grain mask corresponding to the image of the evaluation time phase, an interference spectral dictionary is extracted from the region marked as instantaneous dynamic sand grain pixels, and a background spectral dictionary is extracted from the unmarked region.

[0026] Based on the interference spectral dictionary and the background spectral dictionary, the original spectral signal of each pixel is sparsely decomposed to obtain the decomposed background signal;

[0027] The background signals of all pixels are arranged according to their spatial positions to generate the background signal feature map.

[0028] Specifically, the step of inputting the background signal feature map and the surface dynamic disturbance index into a pre-trained anti-dynamic disturbance assessment model to obtain the preliminary planting carrying capacity assessment results of the target aeolian sandy land assessment area at the assessment time phase includes:

[0029] The background signal feature map is input into the first convolutional neural network branch of the anti-dynamic interference evaluation model to extract spatial spectral features;

[0030] Spatially expand the surface dynamic disturbance index to generate a dynamic attention weight map with the same size as the spatial spectral feature;

[0031] The spatial spectral features are multiplied element-wise with the dynamic attention weight map to obtain the weighted spatial spectral features;

[0032] The weighted spatial spectral features are input into the fully connected layer of the anti-dynamic interference assessment model, and the preliminary planting carrying capacity assessment result corresponding to each pixel is output. The preliminary planting carrying capacity assessment result is a preset carrying capacity level label.

[0033] Specifically, the step of determining at least one clean phase from the temporal image stack based on the instantaneous dynamic sand mask; and obtaining the corresponding historical carrying capacity assessment results based on the images of the clean phase and through the signal decoupling model and the anti-dynamic interference assessment model, includes:

[0034] For each phase in the time-series image stack, the proportion of pixels marked as instantaneous dynamic sand particles in the instantaneous dynamic sand particle mask is calculated, and the phases with the proportion lower than the second threshold are determined as the clean phases.

[0035] For each of the clean phases, based on the image of the clean phase and its corresponding instantaneous dynamic sand mask, the background signal feature map corresponding to the clean phase is obtained through the signal decoupling model.

[0036] The background signal feature map and the surface dynamic interference index corresponding to each clean time are input into the anti-dynamic interference assessment model to obtain the historical carrying capacity assessment result corresponding to each clean time.

[0037] Specifically, the step of using the historical carrying capacity assessment results to perform temporal consistency correction on the preliminary planting carrying capacity assessment results, and outputting the final planting carrying capacity distribution map of the target sandy land assessment area, includes:

[0038] Calculate the mean and standard deviation of the historical bearing capacity assessment results corresponding to all clean times at the same spatial location pixel;

[0039] For each pixel value in the preliminary planting bearing capacity assessment results at the assessment time phase, calculate the absolute value of the deviation between it and the mean value of the results at the corresponding pixel location;

[0040] When the absolute value of the deviation is greater than a preset multiple of the standard deviation of the result, the pixel value in the preliminary planting bearing capacity assessment result is replaced with the mean value of the result; otherwise, the original value is retained, and a corrected bearing capacity assessment result matrix is ​​generated.

[0041] The corrected bearing capacity assessment result matrix is ​​output as the final planting bearing capacity distribution map.

[0042] Secondly, the present invention provides a system for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing. The assessment system applies the assessment method described in the first aspect, and the assessment system includes:

[0043] The data acquisition and construction module is used to acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a preset time period, with the time interval between adjacent temporal phases being less than a preset threshold. The multi-temporal multispectral remote sensing images form a time-series image stack, which is composed of multiple pixels.

[0044] The dynamic interference identification module is connected to the data acquisition and construction module. The dynamic interference identification module is used to extract the instantaneous dynamic sand mask of each phase of the target sandy land assessment area based on the spectral temporal change characteristics of each pixel in the temporal image stack, and to calculate the pixel-level surface dynamic interference index.

[0045] The background signal decoupling module is connected to the data acquisition and construction module and the dynamic interference identification module. The background signal decoupling module is used to determine an evaluation phase from the time-series image stack, and based on the image of the evaluation phase and its corresponding instantaneous dynamic sand mask, the background signal characterizing the land stability attribute is separated from the original spectral signal of each pixel through the signal decoupling model to generate a background signal feature map.

[0046] The preliminary bearing capacity assessment module is connected to the background signal decoupling module and the dynamic interference identification module. The preliminary bearing capacity assessment module is used to input the background signal feature map and the surface dynamic interference index into the pre-trained anti-dynamic interference assessment model to obtain the preliminary planting bearing capacity assessment results of the target sandy land assessment area in the assessment phase.

[0047] The historical results acquisition module is connected to the data acquisition and construction module, the background signal decoupling module, the dynamic interference identification module, and the preliminary bearing capacity assessment module. The historical results acquisition module is used to determine at least one clean phase from the time-series image stack based on the instantaneous dynamic sand mask, and to obtain the corresponding historical bearing capacity assessment results based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference assessment model.

[0048] The time-domain correction and output module is connected to the preliminary bearing capacity assessment module and the historical result acquisition module. The time-domain correction and output module is used to perform time-domain consistency correction on the preliminary planting bearing capacity assessment results using the historical bearing capacity assessment results, and output the final planting bearing capacity distribution map of the target sandy land assessment area.

[0049] Specifically, the data acquisition and construction module includes:

[0050] The image acquisition submodule is used to acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a complete windy season cycle, with the time interval between adjacent temporal phases being less than or equal to 7 days.

[0051] The image preprocessing submodule is connected to the image acquisition submodule. The image preprocessing submodule is used to perform radiometric calibration and atmospheric correction on the multi-temporal multispectral remote sensing image to obtain the preprocessed multi-temporal multispectral remote sensing image.

[0052] The stacking construction submodule is connected to the image preprocessing submodule. The stacking construction submodule is used to stack the preprocessed multi-temporal multispectral remote sensing images in chronological order to form the temporal image stack.

[0053] Specifically, the dynamic interference identification module includes:

[0054] The feature calculation submodule is used to calculate the normalized vegetation index and normalized sand grain index of each pixel in the temporal image stack at each time phase, so as to form the spectral temporal variation feature of the pixel.

[0055] A mask generation submodule is connected to the feature calculation submodule. The mask generation submodule is used to detect the abrupt change amplitude of the normalized sand grain index of each pixel between adjacent time phases based on the spectral temporal change characteristics, mark the pixels with abrupt change amplitude exceeding a first threshold as instantaneous dynamic sand grain pixels, and generate the instantaneous dynamic sand grain mask for each time phase.

[0056] An index calculation submodule is connected to the mask generation submodule. The index calculation submodule is used to calculate the total frequency of each pixel that is marked as the instantaneous dynamic sand grain pixel in the corresponding time period of the time-series image stack, and normalize the total frequency as the surface dynamic interference index of the pixel.

[0057] This application provides a method and system for assessing the planting carrying capacity of aeolian sandy land based on multispectral remote sensing. The method first acquires multi-temporal multispectral remote sensing images of the target aeolian sandy land assessment area within a preset time period, with adjacent temporal intervals less than a preset threshold, forming a temporal image stack composed of multiple pixels. Using a dynamic interference identification model, the instantaneous dynamic sand grain mask for each temporal phase is extracted based on the temporal variation characteristics of each pixel's spectral data, and the pixel-level surface dynamic interference index is calculated. After determining the assessment temporal phase, a background signal feature map is generated by separating the background signal from the original spectral signal using a signal decoupling model, and input into the anti-dynamic interference assessment model to obtain preliminary assessment results. Then, based on the instantaneous dynamic sand grain mask, a clean temporal phase is determined, historical carrying capacity assessment results are obtained, and temporal consistency correction is performed on the preliminary results. Finally, a planting carrying capacity distribution map of the target aeolian sandy land assessment area is output, effectively solving the problem of spectral signal aliasing and improving the reliability and value of the assessment. Attached Figure Description

[0058] The accompanying drawings, which form part of this application, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0059] Figure 1 A flowchart illustrating the method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing provided in this application;

[0060] Figure 2 A schematic diagram of the connection of the wind-blown sandy land planting carrying capacity assessment system based on multispectral remote sensing provided in this application.

[0061] The accompanying drawings illustrate specific embodiments of this application, which will be described in more detail below. These drawings and descriptions are not intended to limit the scope of the concept in any way, but rather to illustrate the concept of this application to those skilled in the art through reference to particular embodiments. Detailed Implementation

[0062] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0063] The terms "first," "second," "third," "fourth," etc. (if present) in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein.

[0064] In this invention, the terms "exemplary" or "for example" are used to indicate examples, illustrations, or descriptions. Any embodiment or design described as "exemplary" or "for example" in this application should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of terms such as "exemplary" or "for example" is intended to present the relevant concepts in a concrete manner.

[0065] This application provides a method and system for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing. The method first acquires multi-temporal multispectral remote sensing images of the target area within a preset time period, with short intervals between adjacent phases, forming a time-series image stack. A dynamic interference identification model is used to extract instantaneous dynamic sand grain masks and calculate the interference index. A signal decoupling model separates the background signal from the assessment phase images to generate a feature map, which, combined with the interference index, yields preliminary assessment results. Then, clean phases are identified to obtain historical results, the preliminary results are corrected, and a final carrying capacity distribution map is output, thus solving the problem of spectral signal aliasing and accurately assessing the intrinsic carrying capacity of the land.

[0066] The technical solution of this application and how the technical solution of this application solves the above-mentioned technical problems are described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments. The embodiments of this application will now be described with reference to the accompanying drawings.

[0067] Figure 1The flowchart of the method for assessing the planting carrying capacity of aeolian sandy land based on multispectral remote sensing provided in this application is shown below. Figure 1 As shown, this embodiment provides a method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing. The method includes:

[0068] S101: Acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a preset time period, with adjacent time intervals less than a preset threshold. The multi-temporal multispectral remote sensing images form a time-series image stack, which consists of multiple pixels.

[0069] Specifically, the acquisition of multi-temporal multispectral remote sensing images of the target sandstorm assessment area within a preset time period, with the time interval between adjacent time phases being less than a preset threshold, and the formation of a time-series image stack from the multi-temporal multispectral remote sensing images, includes: acquiring multi-temporal multispectral remote sensing images of the target sandstorm assessment area within a complete windy season cycle, with the time interval between adjacent time phases being less than or equal to 7 days; performing radiometric calibration and atmospheric correction processing on the multi-temporal multispectral remote sensing images to obtain preprocessed multi-temporal multispectral remote sensing images; and stacking the preprocessed multi-temporal multispectral remote sensing images in chronological order to form the time-series image stack.

[0070] The specific steps of implementation S101 include:

[0071] 1. Determine the data source and time window: Select Landsat 8 satellite OLI sensor data covering the target aeolian sandy land assessment area as the data source for multispectral remote sensing imagery. Set the preset time period to a complete local windy season cycle, for example, from March 1, 2023 to May 31, 2023. Within this time period, filter all usable images with cloud cover below 10%, ensuring that the time interval between adjacent phases of the filtered images is less than or equal to 7 days. The final result is a list of images arranged in chronological order, for example, a date sequence of: 2023-03-05, 2023-03-12, 2023-03-20, …, 2023-05-28.

[0072] 2. Perform radiometric calibration: Perform radiometric calibration on each original multispectral remote sensing image in the above list. Specifically, read the digital quantization (DN) values ​​for each band of the image, and use the radiometric calibration coefficients (including gain coefficients and offsets) provided in the Landsat 8 data product, along with the formula... = ML*Qcal + AL converts the DN value to the radiance value of the top layer of the atmosphere at the sensor's entrance pupil. Where ML is the band-specific gain coefficient, AL is the band-specific offset, and Qcal is the DN value of the pixel.

[0073] 3. Perform Atmospheric Correction: Perform atmospheric correction on each image after radiometric calibration to eliminate the effects of atmospheric scattering and absorption, converting the radiance values ​​into true surface reflectance. This step uses the FLAASH atmospheric correction module in ENVI software. In the module, input the sensor type as "Landsat-8 OLI", input the radiance image obtained in the previous step, set the latitude and longitude of the image center point, imaging time, average altitude, and other geographical parameters, select the mid-latitude summer atmospheric model and the rural aerosol model, and run the correction. After processing, the surface reflectance image of each image is obtained. The value of each pixel in each band is the reflectance value, ranging from 0 to 1.

[0074] 4. Constructing a Temporal Image Stack: Spatially register all surface reflectance images (i.e., pre-processed multi-temporal multispectral remote sensing images) obtained after radiometric calibration and atmospheric correction to ensure complete spatial alignment of images from different time phases. Then, stack all these images into a single multidimensional data volume according to the order of image acquisition, thus forming the temporal image stack. This temporal image stack can be viewed as a four-dimensional array with dimensions (T, H, W, B), where T represents the number of time phases (i.e., the number of image scenes), H represents the number of rows in the image, W represents the number of columns in the image, and B represents the number of multispectral bands. The temporal image stack consists of H×W pixels, and each pixel contains a spectral temporal sequence of length T×B.

[0075] This step involves acquiring high-frequency, temporally continuous multispectral remote sensing images and performing precise radiometric calibration and atmospheric correction to construct a time-series image stack with high radiometric consistency and a clear temporal dimension. This data foundation ensures that subsequent steps can effectively capture and analyze the dynamic changes in surface spectral characteristics over time, especially rapidly occurring wind and sand cover changes. It provides reliable, high-quality input data for accurately identifying transient dynamic disturbances and separating stable land background attributes, and is a prerequisite for the implementation of the entire assessment method.

[0076] S102: Based on the spectral temporal variation characteristics of each pixel in the temporal image stack, a dynamic interference identification model for identifying instantaneous sandstorm cover based on the spectral temporal variation characteristics is used to extract the instantaneous dynamic sand mask of each phase of the target sandstorm assessment area and calculate the pixel-level surface dynamic interference index.

[0077] Specifically, based on the spectral temporal variation characteristics of each pixel in the temporal image stack, a dynamic interference identification model for identifying instantaneous wind and sand cover based on the spectral temporal variation characteristics is used to extract the instantaneous dynamic sand grain mask of each time phase of the target wind and sand assessment area and calculate the pixel-level surface dynamic interference index. This includes: for each pixel in the temporal image stack, calculating its normalized vegetation index and normalized sand grain index in each time phase to form the spectral temporal variation characteristics of the pixel; based on the spectral temporal variation characteristics, detecting the abrupt change amplitude of the normalized sand grain index of each pixel between adjacent time phases, marking pixels with abrupt change amplitude exceeding a first threshold as instantaneous dynamic sand grain pixels, and generating the instantaneous dynamic sand grain mask for each time phase; for each pixel, counting the total frequency of being marked as instantaneous dynamic sand grain pixels in the corresponding time period of the temporal image stack, and normalizing the total frequency as the surface dynamic interference index of that pixel.

[0078] The specific steps of implementation S102 include:

[0079] 1. Calculate pixel-level spectral temporal variation characteristics: For each pixel in the temporal image stack, the dimensions of the temporal image stack are (T, H, W, B), where T is the number of time phases, H is the image height (number of rows), W is the image width (number of columns), and B is the number of bands. Traverse each time phase t (t=1, 2, …, T), and use the pixel reflectance value of the corresponding image for that time phase to calculate two spectral indices:

[0080] Normalized Difference Vegetation Index (NDVI): Calculated using the formula NDVI(t) = (ρ_nir(t) - ρ_red(t)) / (ρ_nir(t) + ρ_red(t)). Where ρ_nir(t) represents the reflectance of the pixel in the near-infrared band (band 5 for Landsat 8 OLI) at time t, and ρ_red(t) represents the reflectance of the pixel in the red band (band 4 for Landsat 8 OLI) at time t.

[0081] Normalized Difference Sand Index (NDSI): Calculated using the formula NDSI(t) = (ρ_swir1(t) - ρ_green(t)) / (ρ_swir1(t) + ρ_green(t)). Where ρ_swir1(t) represents the reflectance of the pixel in the shortwave infrared band 1 (band 6 for Landsat 8 OLI) at time t, and ρ_green(t) represents the reflectance of the pixel in the green band (band 3 for Landsat 8 OLI) at time t.

[0082] For a pixel at spatial location (i, j), its NDVI and NDSI values ​​over the entire time series together constitute the spectral temporal variation characteristics of the pixel, which can be represented as two vectors of length T: NDVI_vector(i,j)=[NDVI(1), NDVI(2), …, NDVI(T)] and NDSI_vector(i, j) = [NDSI(1), NDSI(2), …, NDSI(T)].

[0083] 2. Detecting abrupt changes and generating a transient dynamic sand grain mask: Based on the NDSI time-series vector NDSI_vector(i, j) of each pixel, detect spectral abrupt changes of the pixel between adjacent time phases. Specifically, for time phase t (t from 2 to T), calculate the NDSI difference (abrupt change magnitude) of the pixel between time phase t and time phase t-1: ΔNDSI(i, j, t) = |NDSI(t) -NDSI(t-1)|. Set a first threshold, for example, Threshold1 = 0.2. Iterate through all pixels (i, j) and all time phases t (t>=2) and make a judgment: if ΔNDSI(i, j, t) > Threshold1, then mark the pixel at time phase t as a "transient dynamic sand grain pixel", otherwise mark it as a "non-transient dynamic sand grain pixel". For each time phase t, a binary image is generated based on the labeling results of all pixels in that time phase, where the position labeled "instantaneous dynamic sand grain pixel" has a value of 1, otherwise it has a value of 0. This binary image is the instantaneous dynamic sand grain mask Mask(t) corresponding to that time phase t, and its size is H × W.

[0084] 3. Calculate the pixel-level surface dynamic disturbance index: For each pixel (i, j), count the total number of times it is marked as a "transient dynamic sand grain pixel" across all T time phases. Specifically, iterate through all Mask(t) images from time phase 2 to time phase T, checking if the value of pixel (i, j) in each Mask(t) is 1, and sum the results to obtain the total count_sand(i, j). Since mutation detection starts from time phase 2, the maximum possible count is (T-1). Then, normalize Count_sand(i, j) to map it to the [0, 1] interval, which serves as the surface dynamic disturbance index for that pixel. The normalization formula is: DI(i, j) = Count_sand(i, j) / (T - 1). Here, DI(i, j) is the pixel-level surface dynamic disturbance index; the closer its value is to 1, the higher the frequency of transient wind and sand cover disturbance at that pixel location during the observation period. Finally, a surface dynamic disturbance index map DI_Map with a size of H×W is generated.

[0085] This step is crucial for identifying and quantifying dynamic surface disturbances. By calculating the NDSI and detecting its abrupt changes between adjacent time phases, it creatively distinguishes between spectral changes caused by slow variations such as vegetation growth and soil moisture, and rapid spectral abrupt changes caused by transient windblown sand cover / removal, thus accurately generating a transient dynamic sand mask for each moment. Furthermore, by statistically analyzing and normalizing the frequency of pixels marked as interfering pixels throughout the entire time series, a quantified surface dynamic disturbance index is generated. This index objectively reflects the frequency of windblown sand activity at various points on the surface, providing direct and critical quantitative evidence for subsequent steps to distinguish between "stable background" and "dynamic noise" in the signal. It is the primary and key step in solving the problem of misjudging transient disturbances.

[0086] S103: Determine an evaluation phase from the time-series image stack; based on the image of the evaluation phase and the instantaneous dynamic sand mask corresponding to the image of the evaluation phase, separate the background signal characterizing the land stability attribute from the original spectral signal of each pixel through a signal decoupling model, and generate a background signal feature map.

[0087] Specifically, the step of determining an evaluation time phase from the temporal image stack; based on the image of the evaluation time phase and the instantaneous dynamic sand mask corresponding to the image of the evaluation time phase, using a signal decoupling model, separating the background signal characterizing land stability attributes from the original spectral signals of each pixel, and generating a background signal feature map, includes: selecting the latest image from the temporal image stack as the evaluation time phase; constructing the original spectral signal from the original multispectral band reflectance values ​​of each pixel in the image of the evaluation time phase; using a signal decoupling model based on sparse representation, extracting an interference spectral dictionary from the region marked as instantaneous dynamic sand grain pixels and extracting a background spectral dictionary from the unmarked region based on the instantaneous dynamic sand grain mask corresponding to the image of the evaluation time phase; performing sparse decomposition on the original spectral signal of each pixel based on the interference spectral dictionary and the background spectral dictionary to obtain the decomposed background signal; and arranging the background signals of all pixels according to their spatial positions to generate the background signal feature map.

[0088] The specific steps in step S103 during implementation include:

[0089] 1. Determine the evaluation time phase: From the time-series image stack constructed in step S101, select the image with the latest acquisition time (i.e., the newest) as the evaluation time phase image. Assume the time-series image stack has T time phases, then the evaluation time phase corresponds to the time phase index T. The evaluation time phase image Image_T is a data set with dimensions (H, W, B), where H is the number of rows, W is the number of columns, B is the number of bands, and each element is the surface reflectance.

[0090] 2. Extract the original spectral signal: Extract the reflectance values ​​of all B bands for each pixel at spatial location (i, j) in the evaluation temporal image Image_T, forming a B-dimensional column vector. The symbol “^T” is a standard symbol in linear algebra, representing transpose. This vector x(i,j) is the original spectral signal of the pixel at the evaluation phase, which contains mixed information of inherent land properties and possible transient wind and sand cover properties.

[0091] 3. Construct a sample set and learn a spectral dictionary: Using the instantaneous dynamic sand grain mask Mask_T (a binary image of size H×W, where a value of 1 represents "instantaneous dynamic sand grain pixel" and a value of 0 represents "non-instantaneous dynamic sand grain pixel") generated in step S102 and corresponding to the evaluation time phase T, two types of training samples are separated from the evaluation time phase image:

[0092] Interference spectral sample set: Collect the original spectral signal vectors x(i,j) corresponding to all pixel positions marked as 1 in Mask_T, and form a matrix X_noise.

[0093] Background spectral sample set: Collect the original spectral signal vectors x(i,j) corresponding to all pixel positions marked as 0 in Mask_T, and form a matrix X_base.

[0094] A dictionary learning algorithm based on sparse representation is adopted, specifically using the K-SVD algorithm, to train X_noise and X_base respectively:

[0095] Using X_noise as input, an interference spectrum dictionary D_n is trained, with a size of B x K_n, where K_n is the number of atoms in the dictionary (for example, K_n=50). Each atom is a B-dimensional vector, representing a typical instantaneous wind and sand cover spectrum pattern.

[0096] Using X_base as input, a background spectral dictionary D_b is trained, with a size of B x K_b, where K_b is the number of dictionaries (e.g., K_b=100), and each atom represents a typical stable surface spectral pattern (such as vegetation, bare soil, water bodies, etc.).

[0097] 4. Solving the background signal using joint dictionary sparse decomposition: For each pixel in the evaluation temporal image, its original spectral signal x is modeled as a joint linear representation of the background dictionary D_b and the interference dictionary D_n, plus a small residual term e, i.e., x ≈ D_b * α_b + D_n * α_n + e. Here, α_b is a sparse coefficient vector of dimension K_b (most elements are 0), and α_n is a sparse coefficient vector of dimension K_n. The signal is separated by solving the following sparse coding problem:

[0098] .

[0099] in, represents the L0 norm (the number of non-zero elements), where L is the sparsity constraint (e.g., L=10). argmin represents the parameter minimization (Argument of the Minimum). st indicates the subject to, which introduces the condition that the optimization problem must satisfy, i.e., the sparsity constraint. and Given constraints, this is the set of sparse coefficients that minimizes the objective function. This optimization problem can be approximated using the orthogonal matching pursuit algorithm. The sparse coefficients are then obtained. and Then, the background signal s_base of that pixel is transmitted through... * The calculation yields s_base, a B-dimensional vector representing the spectral reflectance of the land at the pixel location after removing transient wind and sand interference.

[0100] 5. Generating the Baseline Signal Feature Map: Traverse each pixel (i, j) in the evaluation temporal image Image_T, repeating step 4 to calculate the baseline signal vector s_base(i,j) for each pixel. Then, rearrange the s_base(i,j) of all pixels according to their original spatial positions (i, j) to form a new data cube with dimensions (H, W, B). This data cube is the Baseline Signal Feature Map Base_Feature_Map. This feature map is spatially aligned with the original image, and each channel in the spectral dimension represents the "cleaned" reflectance information of a band.

[0101] This step is the core creative element in solving the problem of spectral signal aliasing. It utilizes the identified transient dynamic sand grain mask as prior knowledge to guide the separation of pure land background attributes from the mixed spectral signals. By learning physically meaningful interference and background spectral dictionaries from both the interfered and undisturbed areas, and using a sparse decomposition model to resolve the spectrum of each pixel into background and interference components, the dynamic noise of "transient wind and sand cover" is quantitatively removed. The generated "background signal feature map" eliminates the spectral disturbances caused by rapidly changing wind and sand cover, providing stable, reliable input features that only reflect the inherent attributes of the land for subsequent carrying capacity assessment. This is crucial to ensuring that the final assessment results are not misled by transient phenomena.

[0102] S104: Input the background signal feature map and the surface dynamic interference index into the pre-trained anti-dynamic interference assessment model to obtain the preliminary planting carrying capacity assessment results of the target sandy land assessment area in the assessment phase.

[0103] Specifically, the step of inputting the background signal feature map and the surface dynamic interference index into a pre-trained anti-dynamic interference assessment model to obtain the preliminary planting carrying capacity assessment result of the target sandy land assessment area at the assessment time phase includes: inputting the background signal feature map into the first convolutional neural network branch of the anti-dynamic interference assessment model to extract spatial spectral features; spatially expanding the surface dynamic interference index to generate a dynamic attention weight map with the same size as the spatial spectral features; performing element-wise multiplication of the spatial spectral features and the dynamic attention weight map to obtain weighted spatial spectral features; inputting the weighted spatial spectral features into the fully connected layer of the anti-dynamic interference assessment model to output the preliminary planting carrying capacity assessment result corresponding to each pixel, wherein the preliminary planting carrying capacity assessment result is a preset carrying capacity level label.

[0104] The specific steps in step S104 during implementation include:

[0105] 1. Constructing a dynamic interference robustness evaluation model: The dynamic interference robustness evaluation model is a deep convolutional neural network, whose structure includes an input layer, a feature extraction branch, an attention fusion module, and an output layer. The specific model construction is as follows:

[0106] Input Layer: Two input ports are designed. The first input port receives the background signal feature map Base_Feature_Map with dimensions (H, W, B), where H is the height, W is the width, and B is the number of bands (feature dimension). The second input port receives the surface dynamic interference index map DI_Map with dimensions (H, W).

[0107] Feature extraction branch (first convolutional neural network branch): This branch processes the Base_Feature_Map. The branch structure is as follows:

[0108] The first convolutional layer uses 64 3×3 convolutional kernels with a stride of 1, employing "same" padding, followed by a ReLU activation function. After the input Base_Feature_Map passes through this layer, the output is a primary feature map with dimensions (H, W, 64).

[0109] First max pooling layer: Uses a 2×2 pooling window with a stride of 2. Input the primary feature map, output a feature map with dimensions (H / 2, W / 2, 64).

[0110] The second convolutional layer uses 128 3×3 convolutional kernels with a stride of 1, employing "same" padding, followed by a ReLU activation function. It takes the pooled feature map as input and outputs a deep feature map with dimensions (H / 2, W / 2, 128).

[0111] The second max-pooling layer uses a 2×2 pooling window with a stride of 2. It takes a deep feature map as input and outputs a feature map with dimensions (H / 4, W / 4, 128). This feature map is denoted as the spatial spectral feature F_spatial.

[0112] Attention Weight Map Generation Module: This module processes the DI_Map. First, the DI_Map is copied and expanded along the channel dimension to generate an intermediate map with dimensions (H, W, 128). Then, this intermediate map is subjected to two 2×2 max pooling operations (with a stride of 2) to reduce its spatial size to (H / 4, W / 4). Next, a linear transformation is performed through a 1×1 convolutional layer (using 128 convolutional kernels, no activation function) to finally generate a dynamic attention weight map W_attention with the exact same dimensions as F_spatial. The value of each position and each channel in W_attention is obtained from the surface dynamic disturbance index of the corresponding position in the original DI_Map after the above transformation. Its value range is maintained between 0 and 1. The larger the value, the higher the potential risk of dynamic disturbance at that position, and the more attention the model should give.

[0113] Feature weighted fusion: The spatial spectral feature F_spatial is multiplied element-wise (Hadamard product) with the dynamic attention weight map W_attention. Specifically, for each spatial location (i,j) and each channel c, F_weighted(i,j,c) = F_spatial(i,j,c) * W_attention(i,j,c). This yields the weighted spatial spectral feature F_weighted, which still has dimensions (H / 4, W / 4, 128). This operation enables the model to adaptively enhance or suppress the responses of different regions and feature channels in the spatial spectral feature based on the surface dynamic disturbance index.

[0114] Fully Connected Layer and Output Layer: The F_weighted tensor is flattened into a one-dimensional vector in spatial dimensions. This vector is input to a fully connected layer with a preset number of bearing capacity levels C (e.g., C=5, corresponding to bearing capacity levels: low, lower, medium, higher, high). A Softmax activation function follows the fully connected layer, outputting a C-dimensional probability vector for each pixel (corresponding to a spatial location in the original F_weighted tensor). This probability vector represents the confidence level of the pixel belonging to each bearing capacity level. The level with the highest probability is taken as the preliminary planting bearing capacity assessment result for that pixel, generating a preliminary assessment result image Prelim_Result with dimensions (H / 4, W / 4). Since the model is downsampled by 4 times, bilinear interpolation is used to upsample Prelim_Result back to the original image size (H, W) to obtain a preliminary planting bearing capacity assessment result Result_T that is spatially aligned with the image at the assessment time.

[0115] 2. Pre-trained dynamic interference resistance evaluation model: Before using the above model for evaluation, it needs to be pre-trained. The training process is as follows:

[0116] Prepare training data: Select several sample areas with known stable surface properties (i.e., unaffected by transient wind and sand disturbances) across multiple historical periods and different wind-blown sand scenarios. For each sample area, perform steps S101 to S103 as described above to generate its background signal feature map and surface dynamic disturbance index map. Simultaneously, through field surveys or high-precision image interpretation, label these sample areas with true-value maps of planting carrying capacity levels (same size as the assessment results, with one level label per pixel).

[0117] Model Training: Using the background signal feature map and the surface dynamic disturbance index map as input, and the bearing capacity level ground value map as the target, the difference between the model's prediction and the ground value is calculated using the classification cross-entropy loss function. The Adam optimizer is used to minimize the loss function, and the model parameters are iteratively updated. During training, the dataset is divided into training, validation, and test sets, and performance is monitored on the validation set to prevent overfitting. When the model's accuracy on the validation set no longer improves, training is stopped, and the optimal model parameters are saved, thus obtaining the pre-trained dynamic disturbance resistance evaluation model.

[0118] 3. Perform preliminary assessment: For the target assessment area, load the pre-trained anti-dynamic interference assessment model described above. Input the background signal feature map of the assessment time phase generated in step S103 and the surface dynamic interference index map generated in step S102 into the model according to the model input requirements. After forward propagation, the model directly outputs the preliminary planting carrying capacity assessment result Result_T of the target sandy land assessment area at the assessment time phase. This result is a two-dimensional matrix, where each element is an integer representing the carrying capacity level label of the corresponding pixel location.

[0119] This step is the core reasoning stage of the evaluation model. It employs a meticulously designed, attention-integrated, dynamic disturbance-resistant evaluation model that organically combines background signal features representing land stability with a surface dynamic disturbance index that quantifies wind and sand disturbance risk. The model uses the dynamic disturbance index as spatial attention, guiding the network to focus more on areas with high disturbance risk or those identified as dynamic zones during feature learning, thus making the evaluation process "resistant" to potential instantaneous wind and sand cover noise. This step ultimately outputs preliminary planting carrying capacity assessment results. These results are based on the "purified" background signal and a comprehensive judgment considering spatial heterogeneity disturbance risks, laying the foundation for obtaining a stable carrying capacity distribution map that reflects the land's intrinsic potential.

[0120] S105: Based on the instantaneous dynamic sand mask, determine at least one clean phase from the time-series image stack; based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference evaluation model, obtain the corresponding historical carrying capacity evaluation result.

[0121] Specifically, the step of determining at least one clean phase from the time-series image stack based on the instantaneous dynamic sand mask, and obtaining the corresponding historical carrying capacity assessment result based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference assessment model, includes: for each phase in the time-series image stack, calculating the proportion of pixels marked as instantaneous dynamic sand particles in the corresponding instantaneous dynamic sand mask, and determining the phase with the proportion lower than a second threshold as the clean phase; for each clean phase, processing the image of the clean phase and its corresponding instantaneous dynamic sand mask through the signal decoupling model to obtain the background signal feature map corresponding to the clean phase; inputting the background signal feature map corresponding to each clean phase and the surface dynamic interference index into the anti-dynamic interference assessment model to obtain the historical carrying capacity assessment result corresponding to each clean phase.

[0122] The specific steps in step S105 during implementation include:

[0123] 1. Calculate the instantaneous interference ratio for each time phase: For each time phase t (t=1, 2, …, T) in the time-series image stack constructed in step S101, obtain the instantaneous dynamic sand grain mask Mask(t) generated in step S102 for the corresponding time phase. Mask(t) is a binary image of size H×W, where a pixel with a value of 1 represents a "instantaneous dynamic sand grain pixel", and a pixel with a value of 0 represents a "non-instantaneous dynamic sand grain pixel". For each time phase t, calculate its instantaneous interference ratio P(t), using the formula:

[0124] P(t) = (sum(Mask(t) == 1)) / (H * W).

[0125] Where sum(Mask(t) == 1) represents the total number of pixels in Mask(t) with a value of 1, and (H * W) is the total number of pixels in the image. The value of P(t) is between 0 and 1, representing the proportion of the area of ​​pixels in the image affected by instantaneous wind and sand cover.

[0126] 2. Determine the clean time phase list: Set a second threshold, denoted as Threshold_clean, to define the degree of "cleanliness". For example, set Threshold_clean = 0.05. Iterate through the instantaneous disturbance ratio P(t) of all time phases, and filter out the time phases t that satisfy the condition P(t) < Threshold_clean to form a clean time phase list Clean_List = [t_c1, t_c2, …, t_cM], where M is the total number of filtered clean time phases, and t_cm represents the index of the m-th clean time phase. The images of these time phases are considered to have minimal impact from transient wind and sand cover on the land surface at the time of acquisition, and can better reflect the relatively stable background state of the land at that time.

[0127] 3. Decouple signals for each clean phase: For each clean phase t_cm in the Clean_List, perform signal decoupling processing similar to but independent of step S103:

[0128] Image and Mask Acquisition: Extract the image Image(t_cm) corresponding to phase t_cm from the temporal image stack, and obtain the corresponding instantaneous dynamic sand grain mask Mask(t_cm).

[0129] Constructing background and interference sample sets: Using Mask(t_cm), extract the spectra of pixels marked as 1 from Image(t_cm) as the interference sample set, and the spectra of pixels marked as 0 as the background sample set.

[0130] Training the spectral dictionary: Using the K-SVD algorithm, based on the interference sample set and the background sample set, the interference spectral dictionary D_n(t_cm) and the background spectral dictionary D_b(t_cm) corresponding to the clean time are trained. This step is the same in principle and parameter settings as the dictionary training for the evaluation phase in S103, but the data of the clean time phase t_cm are used.

[0131] Sparse decomposition to generate feature maps: For each pixel in Image(t_cm), sparse decomposition is performed using its original spectral signal, D_n(t_cm), and D_b(t_cm) with the orthogonal matching pursuit algorithm to solve for the background sparsity coefficients and calculate the background signal of that pixel. After traversing all pixels, a background signal feature map Base_Feature_Map(t_cm) corresponding to the clean phase t_cm is generated, with dimensions (H, W, B).

[0132] 4. Carrying capacity assessment for each clean phase: For each clean phase t_cm, the Base_Feature_Map(t_cm) generated in step 3 and the time-invariant surface dynamic disturbance index map DI_Map generated in step S102 are input into the anti-dynamic disturbance assessment model constructed and pre-trained in step S104. The forward propagation process of the model is exactly the same as the processing of the assessment phase described in S104: the Base_Feature_Map(t_cm) obtains spatial spectral features through the feature extraction branch, the DI_Map generates a dynamic attention weight map, the two are weighted and fused and then passed through a fully connected layer to finally output the carrying capacity assessment result of the clean phase t_cm. This result is denoted as the historical carrying capacity assessment result Hist_Result(t_cm), with dimensions (H, W), and each pixel value is a carrying capacity level label.

[0133] 5. Summarize historical assessment results: Save the Hist_Result(t_cm) obtained from all clean time phases t_cm (m=1, 2, …, M) to form a historical results set {Hist_Result(t_c1), Hist_Result(t_c2), …, Hist_Result(t_cM)}. This set contains the spatial distribution of planting carrying capacity of the target area at multiple "relatively clean" historical moments, providing a multi-temporal reference benchmark for subsequent temporal consistency correction.

[0134] This step is crucial for constructing a time-domain correction benchmark. It uses a "clean time phase" screening mechanism to automatically identify moments in the time series that are minimally affected by transient wind and sand disturbances. The images of these moments undergo independent background signal extraction and carrying capacity assessment processes that are completely consistent with the current assessment time phase. This yields a set of historical carrying capacity "truths" or "baselines" reflecting the land's carrying capacity under "relatively undisturbed" conditions. These historical results are not simply a stack of raw images, but have undergone the same "purification" process of dynamic disturbance identification, signal decoupling, and anti-interference assessment. Therefore, they are more representative of the land's stability potential than transient appearances. This high-quality set of historical benchmarks provides a reliable and consistent time-domain reference framework for the next step of identifying and correcting outliers in the current assessment results that may be caused by residual transient disturbances or assessment errors.

[0135] S106: Use the historical carrying capacity assessment results to perform time-domain consistency correction on the preliminary planting carrying capacity assessment results, and output the final planting carrying capacity distribution map of the target sandy land assessment area.

[0136] Specifically, the step of using the historical carrying capacity assessment results to perform temporal consistency correction on the preliminary planting carrying capacity assessment results and outputting the final planting carrying capacity distribution map of the target sandy land assessment area includes: calculating the mean and standard deviation of the historical carrying capacity assessment results corresponding to all clean times at the same spatial location pixel; for each pixel value in the preliminary planting carrying capacity assessment results of the assessment time phase, calculating the absolute value of the deviation between it and the mean of the result at the corresponding pixel location; when the absolute value of the deviation is greater than a preset multiple of the standard deviation of the result, replacing the pixel value in the preliminary planting carrying capacity assessment results with the mean of the result, otherwise retaining the original value, generating a corrected carrying capacity assessment result matrix; and outputting the corrected carrying capacity assessment result matrix as the final planting carrying capacity distribution map.

[0137] The specific steps in step S106 during implementation include:

[0138] 1. Calculate the statistics of historical results: Suppose that step S105 yields M cleanliness phases, corresponding to M historical carrying capacity assessment results, denoted as the set {Hist_Result_1, Hist_Result_2, …, Hist_Result_M}, where each Hist_Result_m is a matrix of dimension (H, W), and each element is an integer carrying capacity level label. For each pixel location (i, j) in space, based on the values ​​of these M historical results at that location, calculate two statistics:

[0139] Historical Result Mean μ(i, j): Calculates the arithmetic mean of the rank values ​​of M historical results at cell (i, j). The formula is: μ(i, j) = (1 / M) * Σ_{m=1}^{M} Hist_Result_m(i, j). Here, Σ_{m=1}^{M} is a mathematical summation symbol, representing the traversal and summation of historical results from the 1st to the Mth clean phase. Traversing all cells generates a historical result mean matrix Mean_Map of size (H, W).

[0140] Historical Result Standard Deviation σ(i, j): Calculates the sample standard deviation of the rank values ​​of M historical results at pixel (i, j), used to measure the dispersion of historical evaluation results. The formula is: σ(i, j) = sqrt( (1 / (M-1)) * Σ_{m=1}^{M} (Hist_Result_m(i, j) - μ(i, j)) 2 Iterate through all pixels and generate a historical results standard deviation matrix Std_Map of size (H, W).

[0141] 2. Calculate the deviation of the preliminary results: Obtain the preliminary planting capacity assessment result matrix Prelim_Result generated in step S104, which also has dimensions (H, W). For each pixel value P(i, j) in Prelim_Result, calculate the absolute value of the deviation D(i, j) from the historical result mean μ(i, j) at the corresponding position in Mean_Map. The formula is: D(i, j) = |P(i, j) - μ(i, j)|. Traverse all pixels to generate a deviation absolute value matrix Dev_Map of size (H, W).

[0142] 3. Perform temporal consistency correction judgment and replacement: Set a constant K greater than 1 as a preset multiple, for example, set K=2.5. For each pixel position (i, j), perform the following judgment and operation:

[0143] If D(i, j) > K * σ(i, j) holds true, meaning the deviation between the current preliminary assessment result and the historical average value at this location exceeds the threshold of the historical result fluctuation range (K times the standard deviation), then the preliminary result is considered to be affected by residual transient interference or random error of the assessment model at this point, and is an "outlier" that needs correction. In this case, the value of this pixel in the final result is corrected to the historical result mean μ(i, j).

[0144] If D(i, j) <= K * σ(i, j) holds, then the current preliminary assessment result is considered to be consistent with the historical trend and within a reasonable range of fluctuation. Therefore, the preliminary result value P(i, j) is retained as the final value.

[0145] The above judgment and assignment operations are applied to all cells to generate a corrected bearing capacity assessment result matrix Corrected_Result of size (H, W). This operation can be expressed as: for each (i, j), Corrected_Result(i, j) = μ(i, j) if D(i, j) > K*σ(i, j) else P(i, j).

[0146] 4. Output the final distribution map: Output the Corrected_Result matrix generated in step 3 as the final planting capacity distribution map. This map is a two-dimensional spatial data, where each cell contains an integer label representing the planting capacity level after temporal consistency correction. The data format can be a common geospatial raster data format such as GeoTIFF, and it includes geographic coordinate reference information consistent with the original input image, facilitating viewing, analysis, and application in GIS software.

[0147] This step is the final crucial process to ensure the temporal stability and reliability of the assessment results. It utilizes the bearing capacity "baseline" (mean) derived from multiple clean historical moments and its natural fluctuation range (standard deviation) to statistically validate the current assessment results. By replacing outliers in the current results that significantly deviate from historical trends with more representative historical means, "noise points" or "singular values" in the assessment results caused by incompletely eliminated transient wind and sand disturbances, local model misjudgments, or other random factors are effectively filtered out. This correction process ensures that the final output bearing capacity distribution map is not only based on the "cleaning" information of the current time phase but also consistent with the historical performance of regional land attributes, thus exhibiting stronger consistency and robustness in the time dimension and significantly improving the credibility and guiding value of the assessment results in practical applications.

[0148] This embodiment provides a method for assessing the planting carrying capacity of aeolian sandy land based on multispectral remote sensing. The method first acquires multi-temporal multispectral remote sensing images of the target aeolian sandy land assessment area within a preset time period, with adjacent temporal intervals less than a preset threshold, constructing a temporal image stack composed of multiple pixels. Using a dynamic interference identification model, based on the temporal spectral variation characteristics of each pixel in the temporal image stack, the instantaneous dynamic sand mask of each temporal phase is extracted, and the pixel-level surface dynamic interference index is calculated. After determining the assessment temporal phase, a signal decoupling model is used to separate the background signal from the original spectral signals of the assessment temporal phase image and the corresponding mask, generating a background signal feature map. This feature map, along with the surface dynamic interference index, is input into an anti-dynamic interference assessment model to obtain preliminary assessment results. Based on the instantaneous dynamic sand mask, a clean temporal phase is determined, historical carrying capacity assessment results are obtained, and temporal consistency correction is performed on the preliminary results. Finally, a planting carrying capacity distribution map of the target aeolian sandy land assessment area is output, resolving the spectral signal aliasing problem and improving the reliability and value of the assessment.

[0149] Figure 2 A connection diagram of the multispectral remote sensing-based planting carrying capacity assessment system for wind-blown sandy land provided in this application is shown below. Figure 2 As shown, this embodiment provides a multispectral remote sensing-based system for assessing the planting carrying capacity of wind-blown sandy land. This system applies... Figure 1 The method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing described in the embodiment includes an assessment system comprising:

[0150] The data acquisition and construction module is used to acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a preset time period, with the time interval between adjacent temporal phases being less than a preset threshold. The multi-temporal multispectral remote sensing images form a time-series image stack, which is composed of multiple pixels.

[0151] The dynamic interference identification module is connected to the data acquisition and construction module. The dynamic interference identification module is used to extract the instantaneous dynamic sand mask of each phase of the target sandy land assessment area based on the spectral temporal change characteristics of each pixel in the temporal image stack, and to calculate the pixel-level surface dynamic interference index.

[0152] The background signal decoupling module is connected to the data acquisition and construction module and the dynamic interference identification module. The background signal decoupling module is used to determine an evaluation phase from the time-series image stack, and based on the image of the evaluation phase and its corresponding instantaneous dynamic sand mask, the background signal characterizing the land stability attribute is separated from the original spectral signal of each pixel through the signal decoupling model to generate a background signal feature map.

[0153] The preliminary bearing capacity assessment module is connected to the background signal decoupling module and the dynamic interference identification module. The preliminary bearing capacity assessment module is used to input the background signal feature map and the surface dynamic interference index into the pre-trained anti-dynamic interference assessment model to obtain the preliminary planting bearing capacity assessment results of the target sandy land assessment area in the assessment phase.

[0154] The historical results acquisition module is connected to the data acquisition and construction module, the background signal decoupling module, the dynamic interference identification module, and the preliminary bearing capacity assessment module. The historical results acquisition module is used to determine at least one clean phase from the time-series image stack based on the instantaneous dynamic sand mask, and to obtain the corresponding historical bearing capacity assessment results based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference assessment model.

[0155] The time-domain correction and output module is connected to the preliminary bearing capacity assessment module and the historical result acquisition module. The time-domain correction and output module is used to perform time-domain consistency correction on the preliminary planting bearing capacity assessment results using the historical bearing capacity assessment results, and output the final planting bearing capacity distribution map of the target sandy land assessment area.

[0156] Specifically, the data acquisition and construction module includes:

[0157] The image acquisition submodule is used to acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a complete windy season cycle, with the time interval between adjacent temporal phases being less than or equal to 7 days.

[0158] The image preprocessing submodule is connected to the image acquisition submodule. The image preprocessing submodule is used to perform radiometric calibration and atmospheric correction on the multi-temporal multispectral remote sensing image to obtain the preprocessed multi-temporal multispectral remote sensing image.

[0159] The stacking construction submodule is connected to the image preprocessing submodule. The stacking construction submodule is used to stack the preprocessed multi-temporal multispectral remote sensing images in chronological order to form the temporal image stack.

[0160] Specifically, the dynamic interference identification module includes:

[0161] The feature calculation submodule is used to calculate the normalized vegetation index and normalized sand grain index of each pixel in the temporal image stack at each time phase, so as to form the spectral temporal variation feature of the pixel.

[0162] A mask generation submodule is connected to the feature calculation submodule. The mask generation submodule is used to detect the abrupt change amplitude of the normalized sand grain index of each pixel between adjacent time phases based on the spectral temporal change characteristics, mark the pixels with abrupt change amplitude exceeding a first threshold as instantaneous dynamic sand grain pixels, and generate the instantaneous dynamic sand grain mask for each time phase.

[0163] An index calculation submodule is connected to the mask generation submodule. The index calculation submodule is used to calculate the total frequency of each pixel that is marked as the instantaneous dynamic sand grain pixel in the corresponding time period of the time-series image stack, and normalize the total frequency as the surface dynamic interference index of the pixel.

[0164] This embodiment provides a multispectral remote sensing-based system for assessing the planting carrying capacity of wind-blown sandy land, used to perform all the steps described in the aforementioned method embodiments. The assessment system can be physically deployed on one or more interconnected servers, workstations, or high-performance computing terminals, and logically includes multiple functional modules. These modules interact and control processes through a system bus, internal communication interfaces, or preset data exchange protocols. Figure 1 A schematic diagram of the logical architecture of the evaluation system in this embodiment is shown.

[0165] The evaluation system includes: a data acquisition and construction module, a dynamic interference identification module, a background signal decoupling module, a preliminary load-bearing capacity assessment module, a historical result acquisition module, and a time-domain correction and output module.

[0166] 1. Data Acquisition and Construction Module

[0167] The data acquisition and construction module is used to acquire and prepare the raw remote sensing data required for evaluation. This module includes an image acquisition submodule, an image preprocessing submodule, and a stacking construction submodule, which are connected sequentially.

[0168] Image Acquisition Submodule: This submodule connects to a remote sensing data server or local database via a network interface or local storage interface. The image acquisition submodule is configured to initiate query and download requests to the data source based on the geographic boundary coordinates of the target aeolian sand assessment area, either user-input or preset. Request parameters include: spatial extent, temporal range (set to a complete aeolian season cycle, e.g., March to May), sensor type (e.g., Landsat 8 OLI), and cloud cover filtering conditions (e.g., below 10%). The image acquisition submodule ensures that in the final acquired multi-temporal multispectral remote sensing image sequence, the acquisition time interval between any two adjacent images is less than or equal to 7 days. This high-frequency data acquisition strategy is to ensure effective capture of rapid spectral changes on the land surface caused by aeolian sand activity.

[0169] Image Preprocessing Submodule: This submodule connects to the output of the image acquisition submodule and receives raw multi-temporal multispectral remote sensing images. The image preprocessing submodule is configured to perform radiometric calibration and atmospheric correction processing sequentially on each raw image. Radiometric calibration uses calibration coefficients provided in the sensor metadata to convert the image's digital quantization values ​​into atmospheric top-layer radiance values. Atmospheric correction employs the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) atmospheric correction algorithm, taking into account parameters such as radiance imagery, imaging time, center point latitude and longitude, mean altitude, atmospheric model (e.g., mid-latitude summer model), and aerosol model (e.g., rural model), and finally outputs a true surface reflectance image. This processing eliminates the influence of differences in atmospheric and illumination conditions, ensuring radiometric consistency across images from different temporal phases, laying the foundation for temporal series analysis.

[0170] The stacking construction submodule connects to the output of the image preprocessing submodule and receives all preprocessed surface reflectance images. First, the stacking construction submodule calls an image registration algorithm (such as an automatic registration algorithm based on feature points) to precisely align all images to the same coordinate system and pixel grid in space. Then, the stacking construction submodule stacks the registered multi-temporal multispectral remote sensing images in the time dimension according to the order of image acquisition, forming a four-dimensional time-series image stack. The data structure dimensions of the time-series image stack are (T, H, W, B), where T represents the number of temporal phases, H represents the image height (number of rows), W represents the image width (number of columns), and B represents the number of bands. This time-series image stack is the data foundation for all subsequent analysis and processing.

[0171] 2. Dynamic Interference Identification Module

[0172] The dynamic interference identification module is connected to the data acquisition and construction module and receives the time-series image stack. This module is used to identify and quantify instantaneous wind and sand cover interference from temporal changes, and includes a feature calculation submodule, a mask generation submodule, and an index calculation submodule, which are connected sequentially.

[0173] Feature Calculation Submodule: This submodule iterates through each pixel in the temporal image stack. For each pixel, at each time phase t, it calculates two spectral indices using the reflectance value of the image at that time phase: the Normalized Difference Vegetation Index (NDVI) and the Normalized Difference Sand Index (NDSI). The NDVI is calculated using the formula (ρ_nir - ρ_red) / (ρ_nir + ρ_red). The NDSI is calculated using the formula (ρ_swir1 - ρ_green) / (ρ_swir1 + ρ_green). For a pixel, its NDVI and NDSI values ​​across all T time phases together constitute the temporal spectral variation characteristics of that pixel, i.e., two temporal sequences of length T. Calculating these indices highlights the spectral response characteristics of vegetation and sand grains, which is a prerequisite for detecting dynamic changes.

[0174] Mask Generation Submodule: This submodule connects to the feature calculation submodule to obtain the NDSI temporal sequence for each pixel. For each pixel, the mask generation submodule calculates the absolute value of the NDSI difference between adjacent time phases (t and t-1), i.e., the abrupt change amplitude. A first threshold (e.g., 0.2) is set. It iterates through all time phases t (t≥2) and all pixels. If the NDSI abrupt change amplitude of a pixel at time phase t exceeds the first threshold, then in the mask image of time phase t, this pixel location is marked as a "transient dynamic sand grain pixel" (assigned a value of 1); otherwise, it is marked as a "non-transient dynamic sand grain pixel" (assigned a value of 0). For each time phase t, a corresponding binary image is generated, called the transient dynamic sand grain mask for that time phase. This mask accurately indicates the spatial location of the land surface covered by transient sand grains at each moment.

[0175] The index calculation submodule, connected to the mask generation submodule, acquires instantaneous dynamic sand grain masks for all time phases. For each pixel, the submodule counts the total number of times it is marked as a "instantaneous dynamic sand grain pixel" across all masks from time phase 2 to time phase T. This total number is then divided by the maximum possible number (T-1) and normalized to obtain the surface dynamic disturbance index for that pixel, with a value between 0 and 1. A higher index value indicates a higher frequency of instantaneous wind and sand cover disturbance at that pixel location throughout the observation period. After traversing all pixels, a surface dynamic disturbance index map with the same size as the single-time phase image is generated. This index map quantifies the long-term risk of dynamic disturbance at various points on the surface, providing important prior information for subsequent assessment.

[0176] 3. Background signal decoupling module

[0177] The background signal decoupling module is connected to both the data acquisition and construction module and the dynamic interference identification module, receiving the time-series image stack and the corresponding instantaneous dynamic sand grain mask during evaluation. This module is used to separate pure land background information from the interfered mixed spectral signals.

[0178] The background signal decoupling module first selects the latest (newest) image from the time-series image stack as the evaluation time-phase image.

[0179] Next, the background signal decoupling module uses a signal decoupling model based on sparse representation and dictionary learning for processing. Specifically, this module utilizes the instantaneous dynamic sand grain mask of the evaluation phase to separate two types of samples from the evaluation phase image: a set of pixel samples labeled as interference by the mask, and a set of unlabeled pixel samples. Then, the background signal decoupling module calls the K-SVD dictionary learning algorithm, using these two sets of samples as training data respectively, to learn two spectral dictionaries: an interference spectral dictionary (whose atoms represent typical instantaneous sand grain cover spectra) and a background spectral dictionary (whose atoms represent typical stable surface spectra, such as soil and vegetation).

[0180] For each pixel in the assessment time-phase image, the background signal decoupling module models its original spectral signal as a joint linear representation of the background dictionary and the interference dictionary. The module finds a set of sparse coefficients by solving an optimization problem with sparsity constraints (e.g., using an orthogonal matching pursuit algorithm) such that the original signal can be accurately represented with as few dictionary atoms as possible. The signal reconstructed using the background dictionary and the corresponding sparse coefficients is the background signal, representing the land stability attribute, separated from the original signal of that pixel.

[0181] The background signal decoupling module traverses all pixels of the evaluation time-phase image, repeating the above sparse decomposition and reconstruction process. Finally, it rearranges the "purified" background signals of all pixels according to their original spatial positions, generating a background signal feature map. This feature map removes the spectral contribution of instantaneous wind and sand cover and serves as a reliable input for subsequent carrying capacity assessment.

[0182] 4. Preliminary Load-Bearing Capacity Assessment Module

[0183] The preliminary bearing capacity assessment module is connected to both the background signal decoupling module and the dynamic interference identification module, receiving background signal feature maps and surface dynamic interference index maps. This module integrates a pre-trained dynamic interference resistance assessment model for preliminary determination of planting bearing capacity level.

[0184] The dynamic interference resistance assessment model deployed in the preliminary carrying capacity assessment module is a deep convolutional neural network. This model includes two parallel input processing paths. The first path is a feature extraction branch, consisting of several convolutional layers, activation function layers, and pooling layers, used to extract deep spatial spectral features from the input background signal feature map. The second path is an attention weight generation path, used to process the input surface dynamic interference index map: first, it is copied and expanded along the channel dimension; then, its spatial size is adjusted through pooling layers; and finally, it is transformed through a 1x1 convolutional layer, ultimately generating a dynamic attention weight map with the exact same spatial spectral feature size as the output of the first path.

[0185] The model performs element-wise multiplication (Hadamard product) of the extracted spatial spectral features with the dynamic attention weight map. This operation enables the model to adaptively adjust the level of attention given to different regions in the feature map based on the risk of dynamic interference at each location. The contribution of high-interference-risk regions is appropriately suppressed, thereby enhancing the model's robustness against interference.

[0186] The weighted features are fed into a fully connected layer for classification. The number of neurons in the fully connected layer is set to a preset number of planting capacity levels (e.g., 5 levels). The network then outputs the probability of each pixel belonging to each capacity level using the Softmax function, and the level with the highest probability is taken as the preliminary evaluation result for that pixel. The preliminary capacity evaluation module outputs a preliminary planting capacity evaluation result image aligned with the input image space.

[0187] 5. Historical Results Acquisition Module

[0188] The historical results acquisition module is connected to the data acquisition and construction module, the background signal decoupling module, the dynamic interference identification module, and the preliminary bearing capacity assessment module. This module is used to acquire the assessment results of historical "clean" moments as a benchmark for time-domain correction.

[0189] The historical results acquisition module first calculates the proportion of pixels marked as interference in each phase mask based on the instantaneous dynamic sand grain masks generated by the dynamic interference identification module for all phases. A second threshold (e.g., 0.05) is set, and phases with a proportion of interference pixels below this threshold are identified as "clean phases".

[0190] For each defined clean phase, the historical results acquisition module calls the function of the background signal decoupling module (using the image of the clean phase and its corresponding mask) to generate the background signal feature map of that clean phase.

[0191] Then, the historical results acquisition module calls the pre-trained anti-dynamic interference assessment model (the same model used by this module for assessing time phases) in the preliminary carrying capacity assessment module. It inputs the background signal characteristic map of each clean time phase and the unified surface dynamic interference index map into the model, performs forward inference, and obtains the historical carrying capacity assessment results corresponding to each clean time phase. All these historical results are saved, forming a historical assessment result set.

[0192] 6. Time Domain Correction and Output Module

[0193] The time-domain correction and output module is connected to both the preliminary bearing capacity assessment module and the historical results acquisition module, receiving the preliminary planting bearing capacity assessment results and all historical bearing capacity assessment results. This module performs time-domain consistency correction and outputs the final result.

[0194] For each pixel location, the temporal correction and output module calculates the historical result mean (μ) and historical result standard deviation (σ) based on M historical evaluation results provided by the historical result acquisition module. The mean represents the average level of historical carrying capacity at that point, and the standard deviation represents its historical fluctuation range.

[0195] Next, the module calculates the absolute value of the deviation between the value of each pixel in the preliminary assessment results and the historical mean of the corresponding location. A constant K (e.g., 2.5) greater than 1 is set as the multiplier. For each pixel, it is determined whether its absolute deviation is greater than K times the historical standard deviation.

[0196] If the deviation is greater than K times the standard deviation, the preliminary assessment result is considered an outlier, potentially deviating from the historical normal range due to residual interference or random errors. The time-domain correction and output module replaces the value of this pixel with its historical mean. If the deviation is less than or equal to K times the standard deviation, the preliminary result is considered reasonable and retained.

[0197] After completing the traversal, judgment, and correction of all pixels, the temporal correction and output module generates the final corrected planting carrying capacity distribution map. This module has a result output interface, which can output the final distribution map in standard geospatial raster data formats such as GeoTIFF to a file system, database, or send it to a client via the network, thereby completing the entire evaluation process and providing users with a stable, reliable, and temporally consistent spatial distribution result of planting carrying capacity in aeolian sandy areas.

[0198] This embodiment discloses a multispectral remote sensing-based system for assessing the planting carrying capacity of aeolian sandy land, addressing the technical problem of spectral signal aliasing and inaccurate assessment caused by rapid dynamic sand cover in aeolian sandy land. The system includes a sequential collaborative data acquisition and construction module, a dynamic interference identification module, a background signal decoupling module, a preliminary carrying capacity assessment module, a historical results acquisition module, and a temporal correction and output module. The system acquires high-frequency time-series remote sensing images, extracts instantaneous dynamic sand grain masks using the dynamic interference identification module, and calculates the surface dynamic interference index. Through the background signal decoupling module, a signal decoupling model based on K-SVD dictionary learning and sparse representation is used to separate the stable background signal of the land from the mixed spectrum. A preliminary carrying capacity assessment is then performed using an anti-dynamic interference assessment model integrating a dynamic attention mechanism. Finally, the system uses assessment results obtained from historical clean phases to perform temporal consistency statistical correction on the preliminary results. This system effectively isolates instantaneous aeolian sand interference, outputting a stable and reliable spatial distribution map of the land's intrinsic planting carrying capacity, significantly improving the accuracy and practicality of remote sensing assessment in dynamic desertification environments.

[0199] Other embodiments of this application will readily occur to those skilled in the art upon consideration of the specification and practice of the invention disclosed herein. This application is intended to cover any variations, uses, or adaptations of this application that follow the general principles of this application and include common knowledge or customary techniques in the art not disclosed herein. The specification and examples are to be considered exemplary only.

[0200] It should be understood that this application is not limited to the precise structure described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope.

Claims

1. A method for evaluating the planting bearing capacity of aeolian sandy land based on multispectral remote sensing, characterized in that, The method includes: Acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a preset time period, with the time interval between adjacent time phases being less than a preset threshold. The multi-temporal multispectral remote sensing images constitute a time-series image stack, which is composed of multiple pixels. Based on the spectral temporal variation characteristics of each pixel in the temporal image stack, a dynamic interference identification model for identifying instantaneous wind and sand cover based on the spectral temporal variation characteristics is used to extract the instantaneous dynamic sand mask of each time phase of the target wind and sand assessment area, and the pixel-level surface dynamic interference index is calculated. An evaluation phase is determined from the time-series image stack; based on the image of the evaluation phase and the instantaneous dynamic sand mask corresponding to the image of the evaluation phase, the background signal characterizing the land stability attribute is separated from the original spectral signal of each pixel through a signal decoupling model, and a background signal feature map is generated. The background signal feature map and the surface dynamic interference index are input into the pre-trained anti-dynamic interference assessment model to obtain the preliminary planting carrying capacity assessment results of the target sandy land assessment area at the assessment time phase. Based on the instantaneous dynamic sand mask, at least one clean phase is determined from the time-series image stack; based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference evaluation model, the corresponding historical carrying capacity evaluation results are obtained; The historical carrying capacity assessment results are used to perform time-domain consistency correction on the preliminary planting carrying capacity assessment results, and the final planting carrying capacity distribution map of the target sandy land assessment area is output.

2. The method for evaluating the wind-blown sand area planting carrying capacity based on multispectral remote sensing according to claim 1, characterized in that, The acquisition of multi-temporal multispectral remote sensing images of the target sandstorm assessment area within a preset time period, with adjacent time intervals less than a preset threshold, constitutes a time-series image stack, including: Acquire multi-temporal, multispectral remote sensing images of the target aeolian sandy land assessment area within a complete windy season cycle, with the time interval between adjacent temporal phases being less than or equal to 7 days; The multi-temporal multispectral remote sensing images are subjected to radiometric calibration and atmospheric correction to obtain preprocessed multi-temporal multispectral remote sensing images. The preprocessed multi-temporal multispectral remote sensing images are stacked in chronological order to form the time-series image stack.

3. The method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing according to claim 1, characterized in that, Based on the spectral temporal variation characteristics of each pixel in the temporal image stack, a dynamic interference identification model is used to identify instantaneous wind and sand cover according to the spectral temporal variation characteristics. This model extracts the instantaneous dynamic sand mask of each temporal phase in the target wind and sand assessment area and calculates the pixel-level surface dynamic interference index, including: For each pixel in the time-series image stack, the normalized vegetation index and normalized sand index are calculated for each time phase to form the spectral temporal variation characteristics of the pixel. Based on the spectral temporal variation characteristics, the abrupt change amplitude of the normalized sand grain index of each pixel between adjacent time phases is detected, and pixels with abrupt change amplitude exceeding a first threshold are marked as instantaneous dynamic sand grain pixels, thereby generating the instantaneous dynamic sand grain mask for each time phase. For each pixel, the total frequency of the pixel marked as the instantaneous dynamic sand grain pixel within the corresponding time period of the time-series image stack is counted, and the total frequency is normalized and used as the surface dynamic interference index of the pixel.

4. The method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing according to claim 1, characterized in that, The process involves determining an evaluation time phase from the time-series image stack; based on the image of the evaluation time phase and the corresponding instantaneous dynamic sand mask, using a signal decoupling model, separating the background signal characterizing land stability attributes from the original spectral signal of each pixel, and generating a background signal feature map, including: The latest image is selected from the time-series image stack as the evaluation time phase; The original multispectral reflectance values ​​of each pixel in the image at the evaluation time phase constitute the original spectral signal; Using a signal decoupling model based on sparse representation, and based on the instantaneous dynamic sand grain mask corresponding to the image of the evaluation time phase, an interference spectral dictionary is extracted from the region marked as instantaneous dynamic sand grain pixels, and a background spectral dictionary is extracted from the unmarked region. Based on the interference spectral dictionary and the background spectral dictionary, the original spectral signal of each pixel is sparsely decomposed to obtain the decomposed background signal; The background signals of all pixels are arranged according to their spatial positions to generate the background signal feature map.

5. The method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing according to claim 1, characterized in that, The step of inputting the background signal feature map and the surface dynamic disturbance index into a pre-trained anti-dynamic disturbance assessment model to obtain the preliminary planting carrying capacity assessment results of the target sandy land assessment area at the assessment time phase includes: The background signal feature map is input into the first convolutional neural network branch of the anti-dynamic interference evaluation model to extract spatial spectral features; Spatially expand the surface dynamic disturbance index to generate a dynamic attention weight map with the same size as the spatial spectral feature; The spatial spectral features are multiplied element-wise with the dynamic attention weight map to obtain the weighted spatial spectral features; The weighted spatial spectral features are input into the fully connected layer of the anti-dynamic interference assessment model, and the preliminary planting carrying capacity assessment result corresponding to each pixel is output. The preliminary planting carrying capacity assessment result is a preset carrying capacity level label.

6. The method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing according to claim 1, characterized in that, The step involves determining at least one clean phase from the temporal image stack based on the instantaneous dynamic sand mask; and obtaining the corresponding historical carrying capacity assessment results based on the images of the clean phase and through the signal decoupling model and the anti-dynamic interference assessment model, including: For each phase in the time-series image stack, the proportion of pixels marked as instantaneous dynamic sand particles in the instantaneous dynamic sand particle mask is calculated, and the phases with the proportion lower than the second threshold are determined as the clean phases. For each of the clean phases, based on the image of the clean phase and its corresponding instantaneous dynamic sand mask, the background signal feature map corresponding to the clean phase is obtained through the signal decoupling model. The background signal feature map and the surface dynamic interference index corresponding to each clean time are input into the anti-dynamic interference assessment model to obtain the historical carrying capacity assessment result corresponding to each clean time.

7. The method for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing according to claim 1, characterized in that, The step of using the historical carrying capacity assessment results to perform time-domain consistency correction on the preliminary planting carrying capacity assessment results, and outputting the final planting carrying capacity distribution map of the target sandy land assessment area, includes: Calculate the mean and standard deviation of the historical bearing capacity assessment results corresponding to all clean times at the same spatial location pixel; For each pixel value in the preliminary planting bearing capacity assessment results at the assessment time phase, calculate the absolute value of the deviation between it and the mean value of the results at the corresponding pixel location; When the absolute value of the deviation is greater than a preset multiple of the standard deviation of the result, the pixel value in the preliminary planting bearing capacity assessment result is replaced with the mean value of the result; otherwise, the original value is retained, and a corrected bearing capacity assessment result matrix is ​​generated. The corrected bearing capacity assessment result matrix is ​​output as the final planting bearing capacity distribution map.

8. A system for assessing the planting carrying capacity of wind-blown sandy land based on multispectral remote sensing, characterized in that, The evaluation system applies the evaluation method according to any one of claims 1-7, and the evaluation system comprises: The data acquisition and construction module is used to acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a preset time period, with the time interval between adjacent temporal phases being less than a preset threshold. The multi-temporal multispectral remote sensing images form a time-series image stack, which is composed of multiple pixels. The dynamic interference identification module is connected to the data acquisition and construction module. The dynamic interference identification module is used to extract the instantaneous dynamic sand mask of each phase of the target sandy land assessment area based on the spectral temporal change characteristics of each pixel in the temporal image stack, and to calculate the pixel-level surface dynamic interference index. The background signal decoupling module is connected to the data acquisition and construction module and the dynamic interference identification module. The background signal decoupling module is used to determine an evaluation phase from the time-series image stack, and based on the image of the evaluation phase and its corresponding instantaneous dynamic sand mask, the background signal characterizing the land stability attribute is separated from the original spectral signal of each pixel through the signal decoupling model to generate a background signal feature map. The preliminary bearing capacity assessment module is connected to the background signal decoupling module and the dynamic interference identification module. The preliminary bearing capacity assessment module is used to input the background signal feature map and the surface dynamic interference index into the pre-trained anti-dynamic interference assessment model to obtain the preliminary planting bearing capacity assessment results of the target sandy land assessment area in the assessment phase. The historical results acquisition module is connected to the data acquisition and construction module, the background signal decoupling module, the dynamic interference identification module, and the preliminary bearing capacity assessment module. The historical results acquisition module is used to determine at least one clean phase from the time-series image stack based on the instantaneous dynamic sand mask, and to obtain the corresponding historical bearing capacity assessment results based on the image of the clean phase and through the signal decoupling model and the anti-dynamic interference assessment model. The time-domain correction and output module is connected to the preliminary bearing capacity assessment module and the historical result acquisition module. The time-domain correction and output module is used to perform time-domain consistency correction on the preliminary planting bearing capacity assessment results using the historical bearing capacity assessment results, and output the final planting bearing capacity distribution map of the target sandy land assessment area.

9. The wind-blown sand land planting carrying capacity assessment system based on multispectral remote sensing according to claim 8, characterized in that, The data acquisition and construction module includes: The image acquisition submodule is used to acquire multi-temporal multispectral remote sensing images of the target sandy land assessment area within a complete windy season cycle, with the time interval between adjacent temporal phases being less than or equal to 7 days. The image preprocessing submodule is connected to the image acquisition submodule. The image preprocessing submodule is used to perform radiometric calibration and atmospheric correction on the multi-temporal multispectral remote sensing image to obtain the preprocessed multi-temporal multispectral remote sensing image. The stacking construction submodule is connected to the image preprocessing submodule. The stacking construction submodule is used to stack the preprocessed multi-temporal multispectral remote sensing images in chronological order to form the temporal image stack.

10. The wind-blown sand land planting carrying capacity assessment system based on multispectral remote sensing according to claim 8, characterized in that, The dynamic interference identification module includes: The feature calculation submodule is used to calculate the normalized vegetation index and normalized sand grain index of each pixel in the temporal image stack at each time phase, so as to form the spectral temporal variation feature of the pixel. A mask generation submodule is connected to the feature calculation submodule. The mask generation submodule is used to detect the abrupt change amplitude of the normalized sand grain index of each pixel between adjacent time phases based on the spectral temporal change characteristics, mark the pixels with abrupt change amplitude exceeding a first threshold as instantaneous dynamic sand grain pixels, and generate the instantaneous dynamic sand grain mask for each time phase. An index calculation submodule is connected to the mask generation submodule. The index calculation submodule is used to calculate the total frequency of each pixel that is marked as the instantaneous dynamic sand grain pixel in the corresponding time period of the time-series image stack, and normalize the total frequency as the surface dynamic interference index of the pixel.

Images