A soil salinity inversion method and system based on thresholded NDVI curve area

Abstract

This invention relates to the field of remote sensing inversion technology, and in particular to a method and system for soil salinity inversion based on the area under the thresholded NDVI curve. The method includes acquiring multi-source optical satellite surface reflectance data of the target growing season in the study area; calculating the Normalized Difference Vegetation Index (NDVI) pixel by pixel after quality control; performing temporal resampling and missing data interpolation on non-equal-interval NDVI observation sequences to generate equal-interval NDVI sequences with a preset time step; using local least squares smoothing to denoise and reconstruct the equal-interval NDVI sequences; traversing and calculating the area under the threshold (AUC) of the smoothed NDVI sequences at each threshold; determining the optimal threshold by combining ground samples and generating a pixel-level optimal AUC raster for the study area; and constructing an SSC-AUC function model based on the measured soil salinity content (SSC) of ground samples and the corresponding optimal AUC value. This invention uses the NDVI time series as the core observation object throughout the process, effectively improving the spatial consistency and accuracy of the inversion results.

Description

Technical Field

[0001] This invention relates to the field of remote sensing inversion technology, and in particular to a method and system for soil salinity inversion based on the area under the thresholded NDVI curve. Background Technology

[0002] Soil salinization is a core obstacle to the efficient use of land resources and sustainable agricultural production. The accumulation of salt in the soil surface and root activity layer reduces soil water potential, causes ion toxicity, inhibits plant water absorption and photosynthesis, directly leading to crop yield reduction, vegetation community degradation, and consequently a decline in ecosystem service functions. Due to the significant spatial heterogeneity and seasonal fluctuations in salinization, traditional monitoring methods combining ground sampling and laboratory measurements are insufficient to meet the needs of dynamic assessment and zoned management of salinization at the regional scale.

[0003] In recent years, optical remote sensing technology has become an important means of monitoring salinization. The Normalized Difference Vegetation Index (NDVI), as a core indicator characterizing vegetation growth status, combined with soil salinity index (SSC) to form a remote sensing inversion technique, has become a research focus in the field of digital soil mapping. Existing soil salinity remote sensing inversion techniques are mainly divided into two categories: one is based on spectral anomalies of bare soil and salt crust to construct a salinity index and establish a direct inversion model; the other uses vegetation signals as an indirect proxy for salt stress, inferring soil salinity levels through vegetation indices. However, both techniques have significant limitations in complex surface environments. Land cover mosaics and mixed pixels are common in areas such as coastal deltas and farmland wetlands. The rapid changes in vegetation cover can obscure the direct spectral signal of soil salinity as vegetation cover increases, significantly reducing the stability of the inversion. Furthermore, factors such as soil moisture content and surface roughness, along with salinity, jointly affect spectral observations, making it difficult for local empirical models to maintain parameter consistency and transferability across different years, land types, and sensors. Existing technologies introduce machine learning to integrate multiple covariates such as topography and meteorology to improve accuracy, but due to their reliance on large amounts of external data and complex model structures, they suffer from insufficient model interpretability and difficulty in ensuring the comparability of features and parameters across sensors. At present, there is a need for a soil salinity inversion method and system based on the area under the thresholded NDVI curve. Summary of the Invention

[0004] To address the technical problems of existing soil salinity inversion methods in complex terrains, such as the presence of time-series observation gaps and empirical threshold settings, which lead to characteristic contamination and difficulty in achieving stable regional-scale inversion, this invention provides a soil salinity inversion method and system based on the area under the thresholded NDVI curve.

[0005] In a first aspect, the present invention provides a soil salinity inversion method based on the area under the thresholded NDVI curve, which adopts the following technical solution: A method for soil salinity inversion based on the area under the thresholded NDVI curve includes: Multi-source optical satellite surface reflectance data of the target growing season in the study area were acquired. After quality control, the Normalized Difference Vegetation Index (NDVI) was calculated pixel by pixel to form a non-equal interval NDVI observation sequence at the pixel level. Time resampling and missing data interpolation are performed on non-equal interval NDVI observation sequences to generate equal interval NDVI sequences with a preset time step; Local least squares smoothing is used to denoise and reconstruct equally spaced NDVI sequences to obtain smoothed NDVI time series; Construct a candidate threshold set for NDVI, iterate through and calculate the area under the curve (AUC) of the smoothed NDVI sequence under each threshold, combine ground samples to determine the optimal threshold and generate a pixel-level optimal AUC raster for the study area. Based on the measured soil salinity content (SSC) at ground sampling points and the corresponding optimal AUC value, an SSC-AUC function model is constructed to obtain the spatial distribution results of SSC and output pixel-level uncertainty labeling information.

[0006] Further, the calculation of the Normalized Difference Vegetation Index (NDVI) pixel by pixel after quality control includes: Acquire at least one optical satellite surface reflectance data covering the target growing season of the study area, excluding months affected by snow cover in the study area; the surface reflectance data is a standardized product that has been radiatively and atmospherically corrected. Pixel-level quality control was performed on the surface reflectance data to remove low-quality observation pixels such as clouds, cirrus clouds, cloud shadows, snow cover, and water bodies, resulting in quality-controlled non-equal interval effective observation data. Normalized Difference Vegetation Index (NDVI) was calculated pixel-by-pixel from the valid observation data. The NDVI values ​​of each pixel were then sorted chronologically by the satellite imaging date to form a non-equidistant NDVI observation sequence at the pixel level for the study area. ,in, For imaging time, The expression for the NDVI value corresponding to the imaging time is: ; In the formula, For near-infrared band reflectivity, This refers to the reflectivity in the red light band.

[0007] Further, generating the equally spaced NDVI sequence with a preset time step includes: A continuous time axis is constructed using the target growing season of the study area as the time range. A preset time step is set and equally spaced time nodes are generated. The preset time step is selected as a daily scale. The start and end times of the target growing season are calculated based on the historical NDVI time series of the study area to obtain the average NDVI sequence. The start date of the target growing season is when the NDVI value in the multi-year average NDVI sequence first continuously exceeds a preset threshold. The end date is the date when the NDVI value last falls below a preset threshold for the longest consecutive period. The date; For each pixel's non-equal interval NDVI observation sequence, time interpolation is used to estimate the NDVI value of the missing time nodes, realizing the mapping of non-equal interval observations to an equal interval time axis; For consecutive missing periods beyond the effective observation range, no extrapolation interpolation is performed. Simultaneously, the missing percentage and longest consecutive missing length for each pixel are statistically analyzed and used as quality markers. Finally, pixel-level NDVI sequences with pre-defined time steps and equal intervals for the study area are generated. ; Among them, threshold and The background NDVI values ​​were set based on the bare soil and typical vegetation in the study area.

[0008] Furthermore, the step of using local least squares smoothing to denoise and reconstruct equally spaced NDVI sequences includes setting the sliding window length m to an odd number, and the polynomial order p being less than the sliding window length m. Based on the equally spaced NDVI sequences, a p-order polynomial is fitted within a sliding window of length m according to the least squares criterion. The polynomial fitting value at the center point of the sliding window is used as the smoothed output value. Full-sequence sliding fitting processing is then performed on the equally spaced NDVI sequences to obtain a smoothed NDVI time series. The expression for the smoothed NDVI time series is: ; in, The sliding window is half the window length. The convolution coefficients are determined by the sliding window length m and the polynomial order p. This represents the smoothed NDVI value at the k-th time node on an equally spaced time axis. This represents the original, equally spaced NDVI values ​​at the (k+j)th time node.

[0009] Furthermore, the step of using local least squares smoothing to denoise and reconstruct equally spaced NDVI sequences also includes setting multiple sets of conditions where m is an odd number and The (m,p) parameter combinations are used to perform local least squares smoothing on the equally spaced NDVI sequences. The curve shapes of the smoothed NDVI time series under different parameter combinations are compared, and the optimal (m,p) parameter combination is selected. The final smoothed NDVI time series is determined based on the smoothing result of the optimal (m,p) parameter combination.

[0010] Furthermore, the step of traversing and calculating the area under the curve (AUC) of the smoothed NDVI sequence at each threshold includes setting the range and step size of the NDVI threshold based on the ecological characteristics of vegetation growth in the study area, and constructing a set of candidate NDVI thresholds. For the smoothed NDVI time series, with candidate thresholds As constraints, calculate the conditions that must be met during the growing season. The area under the curve above the threshold for each time period is calculated using the trapezoidal rule to achieve numerical integration, yielding the pixel-level AUC value for each candidate threshold. The area under the threshold time integral is defined as: ; In the formula, , The start and end times of the target growing season in the study area. It is a function for maximizing the value.

[0011] Furthermore, the step of determining the optimal threshold by combining ground sampling points and generating the optimal AUC raster at the pixel level for the study area includes fitting the AUC values ​​of ground sampling point pixels under each candidate threshold with the measured soil salinity content (SSC) and calculating the coefficient of determination. The optimal threshold is determined by considering the root mean square error (RMSE). Based on the optimal threshold The area under the curve above the threshold is recalculated for the smoothed NDVI time series of all pixels in the study area to generate the optimal AUC raster at the pixel level for the study area. The expression for the optimal AUC raster is: ; in, Time interval The local optimal AUC value within, To preset the time step, , This represents the smoothed NDVI value for the corresponding time point.

[0012] Furthermore, the construction of the SSC-AUC function model includes acquiring ground sampling data of the 0-10cm topsoil layer in the study area and determining the conductivity of the soil leachate using the conductivity method. ,pass The measured SSC values ​​of ground sampling points were obtained by converting the SSC values ​​with the calibration relationship of SSC. The calibration relationship is based on the measured SSC values ​​of soil samples with different salinity gradients and their corresponding values. Values ​​were obtained through regression fitting; The optimal AUC value corresponding to each surface sample pixel is extracted. Using the measured SSC value of the sample points as the dependent variable and the optimal AUC value as the independent variable, an SSC-AUC function model is constructed. The model parameters are then calibrated using the nonlinear least squares method. The expression of the SSC-AUC function model is as follows: ; in, Optimal threshold The AUC value under the following conditions The baseline level of SSC under low salt stress conditions. Let AUC be the sensitivity coefficient to changes in salinity. This is the model bias term.

[0013] Furthermore, obtaining the SSC spatial distribution result and outputting pixel-level uncertainty labeling information includes: The calibrated SSC-AUC function model was applied pixel by pixel to the optimal AUC raster at the pixel level in the study area. The predicted value of soil salinity content for each effective pixel was calculated, and the SSC spatial distribution raster of the study area was generated. At the same time, a classification map was generated from the SSC spatial distribution raster according to the soil salinization classification standard. A pixel-level confidence index is constructed by combining the missing measurement ratio of each pixel, the smoothing residual, and the model fitting residual. Pixels with abnormal NDVI curves are marked as low-confidence pixels. A pixel-level uncertainty marker map is generated and output. The SSC spatial distribution results, the grading map, and the uncertainty marker map are exported into a general raster data format.

[0014] Secondly, a soil salinity inversion system based on the area under the thresholded NDVI curve includes: The data acquisition module is used to perform optical remote sensing image access, radiometric correction, atmospheric correction, quality control and NDVI calculation, and output pixel-level non-equidistant NDVI observation sequences; The time series construction module, connected to the data acquisition module, is used to perform time resampling, missing data interpolation and equal interval sequence generation, and outputs daily-scale equal interval NDVI sequences and missing data quality markers. The smoothing reconstruction module, connected to the timing construction module, is used to perform local least squares smoothing processing and output a smoothed NDVI sequence; The threshold AUC calculation module, connected to the smooth reconstruction module, is used to perform candidate threshold traversal, threshold constraint curve area under area calculation, and optimal threshold determination, and outputs the optimal threshold. and the corresponding AUC grid; The model calibration and inversion module is connected to the threshold AUC calculation module and is used to establish the SSC-AUC function model, calibrate the model parameters, and perform cell-level SSC inference. The mapping output module, connected to the model calibration and inversion module, is used to generate spatial distribution maps, hierarchical statistical maps, and uncertainty marker maps of soil salinity, and supports export in GeoTIFF or NetCDF format.

[0015] In summary, the present invention has the following beneficial technical effects: 1. This invention is based on the soil salinity inversion method of thresholded NDVI curve area. By standardizing and quality-controlling the surface reflectance data of multi-source optical satellites during the growing season of the target area and constructing NDVI sequences, the basic data for salt stress characterization are accurately extracted. Then, through time resampling and missing data interpolation, the non-uniform interval NDVI observation sequence is transformed into a diurnal scale equal interval sequence, which solves the integration bias problem caused by uneven optical satellite observation intervals and missing data, and provides a continuous and stable time series basis for subsequent curve integration. 2. This invention uses constraints on the sliding window length and polynomial order, and employs a combination of local least squares smoothing and multiple parameter combinations to denoise and reconstruct equally spaced NDVI sequences. While effectively suppressing the accumulation and amplification of noise such as cloud residue and atmospheric correction errors, it preserves the seasonal peak shape, amplitude, and phenological inflection point characteristics of the NDVI sequence, thus achieving the construction of a smooth NDVI curve that can be stably integrated.

[0016] 3. This invention constructs a set of NDVI candidate thresholds based on the vegetation ecological characteristics of the study area, and calculates the area under the curve (AUC) of the smooth NDVI sequence under each threshold by combining the trapezoidal method numerical integration traversal. The optimal threshold is determined by fitting the accuracy index, which realizes the accurate division between the effective growth stage and the non-growth background stage. This mechanism avoids the contamination of the cumulative characteristics by low NDVI backgrounds such as bare soil and sparse vegetation, making the AUC feature a core proxy indicator that can stably characterize seasonal salt stress. 4. This invention constructs an exponential SSC-AUC function model by combining the measured SSC values ​​of ground sampling points with the optimal AUC values. It then uses the nonlinear least squares method combined with a robust fitting strategy to calibrate the model parameters. The calibrated model is then applied to the pixel-level optimal AUC raster of the study area, achieving pixel-level accurate inversion of soil salinity content. At the same time, it integrates the pixel missing measurement ratio, smoothing residuals, and model fitting residuals to construct a confidence index and output pixel-level uncertainty labeling information, solving the problem of the lack of reliable expression in existing inversion results and achieving accurate quantification of the quality of inversion results.

[0017] 5. This invention uses NDVI time series as the core observation object throughout the process, without relying on a large number of external covariates such as topography and meteorology, which reduces the complexity of engineering implementation and promotion. Moreover, the optimal threshold and smoothing parameters are determined by data-driven methods, which makes the method consistent and transferable across regions and sensors. In complex surface scenarios where vegetation cover and mixed pixels are common, the method achieves stable inversion of soil salinity content based on the process response characteristics of salt stress during the growing season, effectively improving the spatial consistency and accuracy of the inversion results. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the overall process of a soil salinity inversion method based on the area under the thresholded NDVI curve according to an embodiment of the present invention.

[0019] Figure 2 This is a schematic diagram of NDVI smoothing and threshold AUC calculation in an embodiment of the present invention.

[0020] Figure 3 This is a schematic diagram of the SSC spatial distribution in an embodiment of the present invention.

[0021] Figure 4 This is an architecture diagram of a soil salinity inversion system based on the area under the thresholded NDVI curve in an embodiment of the present invention. Detailed Implementation

[0022] The present invention will be further described in detail below with reference to the accompanying drawings.

[0023] Example 1 Reference Figure 1 This embodiment of a soil salinity inversion method based on the area under the thresholded NDVI curve includes: S1. Acquire multi-source optical satellite surface reflectance data of the target growing season in the study area. After quality control, calculate the normalized vegetation index (NDVI) pixel by pixel to form a pixel-level non-equal interval NDVI observation sequence. S2. Perform time resampling and missing data interpolation on the non-equal interval NDVI observation sequences to generate equal interval NDVI sequences with a preset time step. S3. Local least squares smoothing is used to denoise and reconstruct the equally spaced NDVI sequences to obtain a smoothed NDVI time series. S4. Construct a set of NDVI candidate thresholds, iterate through and calculate the area under the curve (AUC) of the smoothed NDVI sequence under each threshold, combine the ground samples to determine the optimal threshold and generate the pixel-level optimal AUC raster for the study area. S5. Based on the measured soil salinity content (SSC) at ground sampling points and the corresponding optimal AUC value, construct the SSC-AUC function model to obtain the spatial distribution results of SSC and output pixel-level uncertainty labeling information.

[0024] Specifically, a soil salinity inversion method based on the area under the thresholded NDVI curve includes the following steps: like Figure 1 As shown, S1, acquire multi-source optical satellite surface reflectance data of the target growing season in the study area, and calculate the normalized vegetation index (NDVI) pixel by pixel after quality control to form a pixel-level non-equal interval NDVI observation sequence. This embodiment selects at least one type of optical satellite surface reflectance remote sensing data that can completely cover the entire study area where soil salinity is to be inverted. The time range for data selection is strictly limited to the target growing season of the study area. During the delineation of the growing season time window, months affected by snow cover in the study area are excluded through band feature identification and land cover interpretation of remote sensing images to avoid the abnormal surface reflectance caused by snow cover interfering with the temporal stability and feature authenticity of the subsequent NDVI sequence. The acquired optical satellite surface reflectance data is a standardized remote sensing product after radiometric and atmospheric correction. The errors of satellite sensor radiometric calibration, atmospheric scattering, atmospheric absorption and atmospheric refraction have been corrected to ensure the authenticity of the surface reflectance of the remote sensing data.

[0025] In a preferred embodiment of the present invention, the "target growing season" is not a fixed calendar period, but is dynamically determined based on the phenological characteristics of vegetation growth in the study area. The specific determination method is as follows: Based on NDVI time series data from the study area over many years (e.g., the past 5-10 years), a multi-year average NDVI sequence is calculated. By analyzing this sequence, two NDVI thresholds are set. and ,For example( Threshold and Based on the NDVI characteristics of bare soil and typical vegetation in the study area, representative bare soil sample point sets were selected within the study area based on the surface cover type of bare soil. Extract the annual NDVI time-series data for each sample point, calculate the mean and standard deviation of the NDVI background value of bare soil, and set the growing season initiation threshold as follows: To eliminate interference from bare soil background with high confidence, the bare soil sample point set refers to a set of spatial points within the study area that represent bare soil cover characteristics. These characteristics include bare soil surfaces without vegetation cover, bare surfaces with sparse vegetation cover (less than 5%), and areas covered by salt crust. These characteristics are automatically extracted using bare soil indices such as the Bare Soil Index (BSI) for threshold segmentation, and are used to statistically analyze the distribution characteristics of the bare soil background NDVI values. Simultaneously, typical vegetation sample point sets were selected based on the dominant vegetation types in the study area. Extract the NDVI time-series data of the peak period of the growing season for each sample point, and calculate the mean of the peak period of vegetation NDVI. with standard deviation Set the growing season termination threshold to To accurately identify the inflection point of vegetation aging, the dominant vegetation types include the crop types with the largest area in the study area and natural dominant vegetation communities. The sample points are automatically extracted with high NDVI stable pixels through NDVI time-series clustering analysis to statistically analyze the distribution characteristics of NDVI values ​​during the vigorous growth period of vegetation. Through the above-mentioned differential statistical threshold setting, the start and end time of the target growing season are dynamically adapted to the actual vegetation phenological patterns in the study area.

[0026] The start date of the target growing season is defined as: in the multi-year average NDVI sequence, starting from the beginning of the year, when the NDVI value first appears above the threshold for N consecutive days (e.g., N=5). The start date of the target growing season. Correspondingly, the end date of the target growing season is defined as: starting from the end of the growing season, when the NDVI value first appears below the threshold for N consecutive days. The end date. The growing season determined by this method can more accurately correspond to the actual active growth period of vegetation, thereby ensuring that the AUC (area below the threshold) calculated subsequently can more effectively capture the cumulative impact of salt stress on vegetation growth.

[0027] After data acquisition, pixel-level quality control is performed on the optical satellite surface reflectance data for each scene. This operation relies on the remote sensing image quality assessment dataset and pixel attribute recognition algorithm, employing pixel-by-pixel quality judgment and precise removal of abnormal pixels. Specifically: First, the Quality Assignment (QA) data or Scene Classification Layer (SCL) data accompanying the optical satellite surface reflectance data is retrieved. Bilinear interpolation is used to precisely spatially register these data with the surface reflectance data, controlling the registration error within a single pixel, thus establishing a one-to-one spatial matching relationship between the quality assessment information and the reflectance information of each pixel. Using the pre-defined feature codes for clouds, cirrus clouds, cloud shadows, snow cover, and water bodies in the quality assignment, the quality attributes of each pixel in the image are automatically identified and judged. Simultaneously, a reflectance threshold filtering method is used for secondary verification of abnormal pixels. This involves setting reflectance anomaly threshold ranges for the near-infrared and red bands based on the statistical characteristics of the surface reflectance in the study area and industry-standard thresholds. The near-infrared band reflectance threshold range is... The red light band reflectivity threshold range is ,in , These represent the minimum and maximum values ​​of normal surface reflectance in the near-infrared band of the study area, respectively. , These represent the minimum and maximum values ​​of normal surface reflectance in the red band of the study area, respectively, expressed by the formula: ; Identify and label anomalous pixels with high or low reflectance exceeding a threshold range to avoid leaving behind false pixels not identified by quality flags. After identifying and verifying the quality attributes of all pixels in the image, low-quality observation pixels identified as clouds, cirrus clouds, cloud shadows, snow cover, and water bodies are removed from the entire study area. Only valid observation pixels with surface cover types of soil, vegetation, or soil-vegetation mixture are retained. This process is achieved using binary space masking technology, constructing a mask matrix M(x,y). The value rules for the mask matrix are as follows: ; Where (x, y) are the spatial coordinates of the pixel, the effective reflectance data after masking is obtained by performing a pixel-by-pixel product operation between the mask matrix and the surface reflectance data, as shown in the formula: ; in , The effective reflectance data for the near-infrared and red bands after masking are extracted and integrated to form a non-equally spaced effective observation dataset after masking.

[0028] Subsequently, the normalized difference in vegetation index (NDVI) was calculated pixel-by-pixel for the quality-controlled, non-equal-interval effective observation data. Using each pixel as a calculation unit, the surface reflectance values ​​in the near-infrared and red bands were extracted. Following the calculation principle of the normalized difference in vegetation index, the NDVI was calculated using the formula... In the formula, For near-infrared band reflectivity, The red light band reflectance is used in this formula to quantitatively characterize the greenness of the vegetation canopy and the vegetation cover by using the difference and ratio of the reflectance between the near-infrared and red light bands.

[0029] After calculating the NDVI values ​​for all pixels in the entire study area, a time-series sorting operation was performed on the NDVI calculation results for each pixel within the study area. Using the imaging time of the optical satellite as the sorting criterion, the NDVI values ​​of each pixel at different imaging times were organized and sorted in time series, ultimately forming a non-equidistant NDVI observation sequence at the pixel level for the study area. ,in, For imaging time, As the NDVI value corresponding to the imaging time, this non-equal interval NDVI observation sequence realizes the temporal characterization of the vegetation growth status of each pixel in the study area during the target growing season.

[0030] S2. Perform time resampling and missing data interpolation on the non-equal interval NDVI observation sequences to generate equal interval NDVI sequences with a preset time step. Using pixel-level non-equidistant NDVI observation sequences in the study area To input the data, a continuous time axis was first constructed, starting with the start time of the target growing season in the study area. Starting from the timeline and ending at the target growing season. As the endpoint of the timeline, a continuous timeline covering the entire target growing season in the study area is constructed. Simultaneously, a preset time step is set and equally spaced time nodes are generated. In this embodiment, the preset time step is selected as a daily scale, i.e., the time step... Heaven, through formula Calculate and generate all time nodes on an equally spaced time axis, where k is a non-negative integer, taking values ​​of 0, 1, 2, ..., n, and n is determined by the total duration of the target growing season and the time step. Confirmed, satisfied This allows for the standardized and equidistant division of the target growing season time dimension in the study area, providing a unified time benchmark for the time resampling of non-equidistant NDVI observation sequences, and ensuring the consistency and comparability of the subsequent interpolated data in the time dimension.

[0031] After constructing the continuous equidistant time axis, for the non-equidistant NDVI observation sequence of each pixel in the study area, the NDVI value of the missing time nodes on the equidistant time axis is estimated using time interpolation, realizing the mapping from non-equidistant observations to the equidistant time axis. In this embodiment, an adaptive interpolation strategy for missing features is constructed. Based on the missing duration, time distribution characteristics, and sequence noise level of the non-equidistant NDVI observation sequence, linear interpolation, cubic spline interpolation, or weighted average interpolation is adaptively selected to perform the interpolation operation. All interpolation operations are performed on a single pixel as an independent processing unit to ensure the independence and integrity of the NDVI temporal features of each pixel. During the interpolation process, missing estimation is performed only based on the effective observation data within the equidistant time axis, strictly following the interpolation principle of the time dimension to avoid invalid interpolation across time ranges. The specific interpolation implementation method and corresponding formula algorithm are as follows: For the missing segment node where the time interval between adjacent effective observations is less than a preset threshold, the linear interpolation method is selected. A linear fitting relationship is constructed between the NDVI values ​​of two adjacent effective observation points and the time distance to calculate the estimated NDVI value of the missing node. Let the time of the missing node be denoted as . The time and NDVI value of the adjacent preceding and following effective observation points are respectively And satisfy The formula for calculating the estimated NDVI value of the missing node is: ; For nodes with long missing data segments exceeding a preset threshold and exhibiting significant seasonal fluctuations and continuous variations in the NDVI sequence, cubic spline interpolation is selected. Using valid observation points within equally spaced time axes as interpolation nodes, a smooth cubic spline function S(t) is constructed to fit the NDVI time-series characteristics of the valid observation points. This cubic spline function is applied to each sub-interval... The function is a cubic polynomial that satisfies the constraint that the function value, first derivative, and second derivative are continuous at the nodes. The missing nodes are obtained by solving the coefficient matrix of the cubic spline function. Corresponding NDVI estimate This enables accurate estimation of NDVI values ​​at missing nodes, ensuring the smoothness of the interpolated NDVI sequence curve and its consistency with seasonal variation characteristics. For non-equal-interval NDVI observation sequences with high observation noise and high NDVI value dispersion, a time-weighted average interpolation method is selected. Taking the missing node as the center, all valid observation points within its time neighborhood are selected as interpolation primitives. Weighting coefficients are set according to the time distance between the valid observation points and the missing node; the closer the time distance, the larger the weighting coefficient, and vice versa. Let the missing node be tx, and the valid observation points within its time neighborhood be... Then first use the formula Calculate the weighting coefficients for each valid observation point. ,in The minimum constant is used to avoid cases where the denominator is 0. The weight coefficients are then normalized to obtain the normalized weights. Finally, through the formula Calculate the estimated NDVI values ​​for the missing nodes, where, The normalized weight coefficient corresponding to the i-th valid observation point is... is the actual NDVI observation value of the i-th valid observation point, and n is the total number of data points.

[0032] When performing interpolation for missing data, no extrapolation interpolation is performed for consecutive missing data periods outside the effective observation range. That is, if there is no valid NDVI observation data for a certain time period on the equally spaced time axis as the basis for interpolation, then no NDVI value estimation is performed for any time nodes within that time period, and their missing status is retained. Simultaneously, the missing data characteristics of each pixel in the study area are quantitatively statistically analyzed, and the missing data ratio and the longest consecutive missing data length for each pixel are calculated. The missing data ratio is the ratio of the number of missing nodes for a single pixel on the equally spaced time axis to the total number of nodes, calculated using the formula... Calculate, where, The percentage of missing pixels. This represents the number of missing nodes for that pixel. The total number of nodes on the equally spaced time axis; the longest consecutive missing measurement length is the maximum number of consecutive missing measurement nodes for a single pixel on the equally spaced time axis, which is obtained by statistical analysis through a time-series traversal algorithm. The calculated missing measurement ratio and the longest consecutive missing measurement length are used as the quality markers of the pixel and associated with the corresponding pixel attribute information, providing a quantitative basis for the uncertainty expression and low-confidence pixel marking in subsequent steps.

[0033] After completing temporal resampling, missing data interpolation, and missing feature statistics for all pixels, the NDVI values ​​(including valid observations and interpolated estimates) of each pixel at all time points on the equally spaced time axis are integrated and sorted according to the chronological order of the time points. Finally, a pixel-level equally spaced NDVI sequence with a preset time step is generated for the study area. Each pixel in this sequence corresponds to an equally spaced NDVI time series data covering the target growing season, and the NDVI value at each time point corresponds to the value on the equally spaced time axis. Correspondingly, each pixel also carries quality marker information on the missing measurement ratio and the longest consecutive missing measurement length.

[0034] S3. Local least squares smoothing is used to denoise and reconstruct the equally spaced NDVI sequences to obtain a smoothed NDVI time series. like Figure 2 As shown, a local least squares smoothing algorithm is used to denoise and reconstruct equally spaced NDVI sequences. The core principle is to smooth the sequence through low-order polynomial fitting within a sliding window. First, the basic parameter constraints for local least squares smoothing are set: the sliding window length is denoted as m, and m takes the value of a positive odd number; the polynomial order is denoted as p, and the following conditions are met: The constraint conditions are such that the parameter constraint rules can avoid the distortion of time series features caused by polynomial overfitting, while ensuring that the amount of fitted data in the sliding window is sufficient to support the solution of polynomials of the corresponding order.

[0035] After setting the parameter constraints, using the equally spaced NDVI sequences of a single pixel within the study area as independent processing units, local least squares sliding fitting is performed sequentially on the equally spaced NDVI sequences of all pixels. Specifically, a sliding window of length m is constructed based on the time nodes on the equally spaced time axis, where the time nodes covered by the window are... , Let k be the half-window length of the sliding window, and k be the time node number corresponding to the center point of the sliding window. The sliding window is slid node by node along the time axis of the equally spaced NDVI sequence. For each sliding window position, the original equally spaced NDVI values ​​of all time nodes within the window are used as fitting samples. Let the time node corresponding to the j-th offset position within the window be k+j, and the original NDVI value be... , Fit a p-th order polynomial using the least squares criterion. The expression for the p-th order polynomial is: Where x is the offset of the time node within the window relative to the center point k, and y is the fitted NDVI value. These are the coefficients to be determined for a p-order polynomial. To solve for the coefficients of this polynomial, a least-squares objective function is constructed as follows: By taking the partial derivatives with respect to each coefficient and setting them to zero, we obtain the normal equation system. ,in, for 1-order matrix, coefficient vector Sample vector: ; Solving this system of normal equations yields the optimal solution for the polynomial coefficients. The NDVI value is obtained by fitting the time node k corresponding to the center point of the sliding window using the p-order polynomial obtained from the fitting. At this time, the relative offset x of the center point is 0. Substituting the value into the polynomial yields the fitted value. The fitted value is used as the smoothed NDVI value at time node k. After smoothing the position of a single window, the sliding window is then slid along the time axis for one time step. The above least squares fitting and center point value solving operations are repeated until the full sequence sliding fitting of the equally spaced NDVI sequence is completed, and finally the preliminary smoothed NDVI sequence corresponding to each pixel is obtained.

[0036] In this embodiment, the fitting result of local least squares smoothing can be equivalently expressed in the form of convolution, and its core calculation formula is as follows: Where h is the half-window length of the sliding window, calculated from the sliding window length m, using the following formula: , The convolution coefficients are determined by the sliding window length m and the polynomial order p. These coefficients are fixed values ​​and can be pre-calculated using the least squares fitting principle based on the values ​​of m and p. This represents the smoothed NDVI value at the k-th time node on an equally spaced time axis. Let j be the original equally spaced NDVI value at the (k+j)th time node on the equally spaced time axis, where j ranges from 1 to 1. h to h', covering all time points within the entire sliding window, the weighted summation operation within the sliding window can be directly completed using this convolutional calculation formula, quickly obtaining the smoothed NDVI value of each time point.

[0037] To avoid oversmoothing or undersmoothing caused by a single parameter combination and to ensure the optimal balance between denoising effect and temporal feature preservation, this step also sets multiple sets of (m,p) parameter combinations that satisfy the above parameter constraint rules. The above local least squares smoothing process is performed on the equally spaced NDVI sequences of the same pixel to obtain multiple sets of smoothed NDVI sequences corresponding to different parameter combinations. Each set of parameter combinations independently completes the full sequence sliding fitting and the generation of smoothed NDVI sequences, and the technical means, calculation formulas and operation procedures in the fitting process are consistent.

[0038] After smoothing all preset (m,p) parameter combinations, multiple smoothed NDVI sequences were comprehensively compared and screened. The integrity of the curve shape was the core screening criterion. Each smoothed NDVI sequence was verified to retain the seasonal peak shape, amplitude, phenological inflection point, and overall trend of the original equally spaced NDVI sequences. Simultaneously, the denoising effect of the sequences was examined. Smoothed NDVI sequences corresponding to parameter combinations that were over-smoothed (fuzzy temporal features, distorted peak shape and amplitude) or under-smoothed (ineffective suppression of residual noise, severe sequence fluctuations) were eliminated. The optimal (m,p) parameter combination was selected as the one that effectively suppressed residual noise while completely preserving the key features of the original temporal sequence. If multiple parameter combinations met the screening requirements, a secondary selection was performed based on the stability of subsequent AUC calculations. Finally, based on the smoothing result corresponding to the optimal (m,p) parameter combination, the final smoothed NDVI time series for each pixel in the study area was determined. All pixel smoothing processes followed a unified optimal (m,p) parameter combination rule to ensure the consistency of the smoothed NDVI time series across the entire study area.

[0039] S4. Construct a set of NDVI candidate thresholds, iterate through and calculate the area under the curve (AUC) of the smoothed NDVI sequence under each threshold, combine the ground samples to determine the optimal threshold and generate the pixel-level optimal AUC raster for the study area. Using the pixel-level smoothed NDVI time series of the study area as input data, the first step is to construct an NDVI candidate threshold set. This construction is based on the ecological characteristics of vegetation growth in the study area, specifically combining the growth greenness characteristics of the dominant vegetation types and the NDVI background values ​​of bare soil and sparse vegetation. A reasonable range and step size for the NDVI threshold are then set to form the NDVI candidate threshold set. The threshold range needs to consider the definition of the effective growth stage. The lower limit should avoid including low NDVI background values ​​of bare soil and sparse vegetation, while the upper limit should avoid truncating the normal growth NDVI values ​​of dominant vegetation. In this embodiment, the threshold step size is selected with equal intervals. Multiple candidate thresholds are obtained by dividing the value range with equal intervals. All candidate thresholds are denoted as follows: .

[0040] After constructing the candidate threshold set, it is used as the basic constraint for subsequent AUC calculations above the thresholds, providing a unified threshold benchmark for iteratively calculating AUC values ​​under each threshold. Following the completion of the candidate threshold set construction, smoothed NDVI time series data for each pixel within the study area are then processed. For each candidate threshold As independent constraints, iterate through and calculate the conditions that satisfy the growing season under each threshold. The area under the curve (AUC) above the threshold of the time period is theoretically defined as follows: In the formula, , The start and end times of the target growing season in the study area. This is a maximum value function, which retains only the portion of the smoothed NDVI value that is above the candidate threshold. At that time, take When participating in integral calculation, When the value is 0, it is used in the integral calculation to achieve the threshold constraint on the effective growth stage of vegetation in the study area, so as to ensure that the AUC value only quantifies the cumulative characteristics of vegetation growth in the effective growth stage, and eliminates the interference of low NDVI background such as bare soil and sparse vegetation on the characteristic quantity from the mechanism.

[0041] To enable the engineering calculation of this theoretical integral, this step employs the trapezoidal rule to perform numerical integration, dividing the continuous time axis of the target growing season into a preset time step. Divided into several consecutive time intervals ,in , These are two adjacent time nodes on an equally spaced time axis. For each preset time step set in step S2, For each time interval, calculate the area under the curve above the local threshold within that interval. Then, sum the local areas of all time intervals that meet the conditions to obtain the overall AUC value of the pixel under that candidate threshold. After completing the AUC calculation for a single candidate threshold, repeat the above operation with the remaining candidate thresholds as constraints until all candidate thresholds are calculated. The AUC value of each pixel in the study area is calculated by iteratively calculating the AUC value of each pixel under each candidate threshold, thereby realizing the multi-threshold quantization of salt stress features.

[0042] After completing the traversal calculation of pixel-level AUC values ​​under all candidate thresholds, the optimal NDVI threshold was determined by combining the measured soil salinity (SSC) data from ground sampling points in the study area. First, the AUC values ​​corresponding to all ground sampling pixels under each candidate threshold were extracted. Using the measured SSC values ​​of the ground sampling points as the dependent variable and the AUC values ​​of the corresponding sampling pixels under each candidate threshold as the independent variable, SSC-AUC fitting relationships were constructed for each candidate threshold. For each fitting relationship, the coefficient of determination was calculated. The two core accuracy metrics are the root mean square error (RMSE) and the coefficient of determination. This characterizes the degree of linear correlation between the AUC value and the measured SSC value. The closer the value is to 1, the stronger the correlation between the two. The root mean square error (RMSE) characterizes the degree of deviation between the fitted value and the measured value. The smaller the RMSE value, the higher the fitting accuracy. Subsequently, the optimal threshold is determined based on the accuracy index. The selection of the best option is determined by the following principle: The candidate threshold with the largest value and the smallest RMSE value is selected as the optimal NDVI threshold. If multiple candidate thresholds have similar accuracy indices, a secondary selection process is performed, taking into account the vegetation ecological characteristics of the study area, to ensure the optimal threshold. This approach satisfies both the fitting accuracy requirements and the ecological definition logic of the effective growth stage of vegetation in the study area, thus achieving the optimal threshold. Once determined, this threshold becomes the optimal dividing line between the effective growth stage and the non-growth background stage of vegetation in the study area.

[0043] Finally, the optimal AUC raster at the cell level for the study area is generated to determine the optimal threshold. To standardize the constraints, the area under the curve above the threshold was recalculated for the smoothed NDVI time series of all pixels within the study area. The trapezoidal rule was still used for numerical integration in the calculation process. For any time interval on the equally spaced time axis... The formula for calculating the local optimal AUC value within this interval is: ; in, Time interval The local optimal AUC value within, To preset the time step, , To calculate the smoothed NDVI value for the corresponding time node, the calculation first determines whether the time interval meets the requirements. The constraints are defined, and the local optimal AUC value is calculated for time intervals that meet the conditions. For time intervals that do not meet the conditions, the local optimal AUC value is recorded as 0. Then, the local optimal AUC values ​​of all time intervals within the target growing season are summed to obtain the optimal AUC value for each pixel. After calculating the optimal AUC values ​​for all pixels in the study area, the optimal AUC value of each pixel is associated with its corresponding spatial location. The values ​​are then stitched and integrated according to the spatial grid scale of the study area to generate the optimal AUC grid at the pixel level for the study area. The spatial resolution of this grid is consistent with the original optical satellite surface reflectance data. The value of each grid pixel is the optimal AUC value at that location. The value of this grid is negatively correlated with the soil salinity content in the study area and can be directly used as the core input feature for subsequent soil salinity content inversion.

[0044] S5. Based on the measured soil salinity content (SSC) at ground sampling points and the corresponding optimal AUC value, construct the SSC-AUC function model to obtain the spatial distribution results of SSC and output pixel-level uncertainty labeling information.

[0045] like Figure 3 As shown, this step first involves obtaining and calibrating the soil salinity (SSC) content at ground sampling points. For the 0-10cm surface soil sampling points in the study area, the conductivity (EC1:5) of the soil extract was determined using the conductivity method. Specifically, the soil extract was prepared by mixing soil and deionized water at a mass ratio of 1:5. After shaking, settling, and filtering, the conductivity value of the filtrate was measured. To achieve the conversion of electrical conductivity values ​​into soil salinity (SSC), a pre-constructed system was built. The calibration relationship with SSC was obtained through field measurements of soil samples at different salinity gradients. Specifically, representative soil samples with varying degrees of salinization from the study area were selected, and their actual soil salinity (SSC) values ​​were determined using either the drying and weighing method or the total ion method. Simultaneously, the corresponding values ​​for each sample were also measured. Value, in The values ​​are used as independent variables and the measured values ​​of SSC are used as dependent variables to perform regression fitting, thereby obtaining the quantitative calibration relationship between the two. Then, the values ​​of all ground sampling points are measured. Substituting the values ​​into the calibration relationship, the measured soil salinity (SSC) value for each ground sampling point is obtained, thus completing the calibration of the measured data of the ground sampling points and providing real and reliable dependent variable data for the subsequent construction of the SSC-AUC function model.

[0046] After completing the calibration of the measured SSC values ​​of the ground sampling points, the construction and parameter calibration of the SSC-AUC function model are performed. First, based on the spatial raster matching relationship of the study area, the optimal AUC values ​​of each ground sampling point pixel obtained in step S4 are extracted. This ensures that each ground sample point forms a set of matched data pairs of "optimal AUC value - measured SSC value"; the measured SSC value of the ground sample point is the dependent variable, and the optimal AUC value A(T) of the corresponding sample point pixel is the optimal AUC value. Using as the independent variable, an SSC-AUC function model is constructed. In this embodiment, an exponential function model is preferred to accurately characterize the nonlinear response of soil salinity content to changes in AUC value. The model expression is as follows: ,in, Optimal threshold The AUC value under the following conditions The baseline level of SSC under low salt stress conditions. The AUC is the sensitivity coefficient to changes in salinity, characterizing the degree to which changes in AUC affect soil salinity. This is a model bias term used to correct the model's fundamental biases. It employs a nonlinear least squares method to estimate the model's parameters. , , The core of calibration is to construct an objective function that represents the sum of squared residuals between the measured SSC values ​​and the model-predicted SSC values: ; Where n is the number of ground sampling points. Let SSC be the measured SSC value of the i-th ground sampling point. To find the optimal AUC value for the i-th ground sample pixel, the objective function is solved iteratively. Minimize to obtain model parameters , , The optimal estimate.

[0047] After completing the SSC-AUC function model construction, the spatial distribution results of soil salinity in the study area are generated. The calibrated SSC-AUC function model is then applied pixel-by-pixel to the optimal AUC raster of the study area output in step S4. That is, for each valid pixel in the study area, its optimal AUC value is calculated. Substituting the completed function model, the predicted value of soil salinity content for each pixel is calculated. The predicted salinity values ​​of all pixels are then associated with their corresponding spatial raster locations to generate a spatial distribution raster of SSCs for the entire study area. Simultaneously, based on industry standards and actual needs for salinization control and monitoring in the study area, soil salinization grading thresholds are set, and the SSC spatial distribution raster is divided into mild, moderate, and severe salinization levels according to salinity content, generating a soil salinization grading map for the study area, thus achieving a grading and visual representation of soil salinity content.

[0048] Finally, the construction and output of pixel-level uncertainty labeling information are performed. To quantify the reliability of the salinity inversion results for each pixel in the study area, a pixel-level confidence index is constructed by integrating the multi-dimensional characteristic indicators of each pixel in the previous steps. Specifically, this includes the pixel missing measurement ratio and the longest consecutive missing measurement length statistically obtained in step S2, the smoothing residual after local least squares smoothing in step S3, and the fitting residual from the model fitting process in this step. These indicators are quantified and integrated into a unified pixel-level confidence value through analytic hierarchy process (AHP) or weighted summation. The magnitude of the confidence value is positively correlated with the reliability of the pixel inversion results. Based on the confidence index, the study area... The reliability of pixels within the study area is assessed. Pixels with abnormal NDVI curves, excessively high missing measurement ratios, or smoothing residuals or fitting residuals exceeding preset thresholds are marked as low-confidence pixels. A pixel-level uncertainty marker map of the study area is generated, which can intuitively display the inversion reliability level of each pixel and the spatial distribution of low-confidence pixels. At the same time, the SSC spatial distribution raster of the study area, the salinization classification map, and the pixel-level uncertainty marker map are all exported to common raster data formats such as GeoTIFF and NetCDF to ensure that the results can be used in mainstream remote sensing and geographic information processing software, facilitating subsequent salinization monitoring, analysis, and engineering decision-making applications.

[0049] Example 2 The difference between this embodiment and Embodiment 1 is that this embodiment provides a specific simulation experiment; Taking a typical delta salinization area as the simulation study area, surface reflectance data from three mainstream optical satellites—Sentinel-2, Landsat-8 / 9, and MODIS—were selected as the simulation data source. The soil salinity content in this area exhibits significant spatial heterogeneity, and the vegetation types include salt-tolerant vegetation in saline-alkali land, cultivated crops, and forest land, which meets the inversion simulation requirements for complex salinization areas. The simulation data source selected is multi-source optical satellite surface reflectance data from Sentinel-2 (10m resolution), Landsat-8 / 9 (30m resolution), and MODIS (250m resolution) covering the target growing season of 2025 in the study area (excluding months affected by snow cover). All data are standardized products that have undergone radiometric and atmospheric correction. At the same time, the measured data from 108 surface soil ground sampling points of 0-10cm deployed in Example 1 were selected as simulation verification data, as shown in Table 1. Table 1. Calibration accuracy of sample point true values ​​(EC→SSC);

[0050] This simulation experiment obtains the optimal NDVI threshold for three optical satellite data sources through threshold iterative selection and multi-parameter combination screening. The optimal threshold for the Sentinel-2 data source was 0.16, for the Landsat-8 / 9 data source it was 0.15, and for the MODIS data source it was 0.14. The optimal thresholds for each data source were all within a reasonable range of 0.14 to 0.16, which is consistent with the ecological characteristics of the effective growth stage of vegetation in the study area. This indicates that the threshold traversal optimization strategy of the present invention is stable and can adaptively determine the optimal threshold according to the sequence characteristics of different data sources, avoiding the systematic bias of empirical thresholds.

[0051] This simulation experiment is based on the true values ​​of 108 ground sample points. The inversion accuracy index of three data sources, Sentinel-2, Landsat-8 / 9 and MODIS, is calculated. The specific results are shown in Table 2. Table 2. Optimal thresholds and inversion accuracy for different data sources;

[0052] This simulation experiment conducted a correlation analysis between the optimal AUC values ​​and the true SSC values ​​of the sample points from three data sources. The results showed that the AUC values ​​and SSC values ​​from all three data sources were significantly negatively correlated. That is, the higher the soil salinity, the lower the cumulative characteristic (AUC) above the NDVI threshold during the effective vegetation growth stage, which is highly consistent with the ecological response pattern of salt stress. This result proves that the threshold-constrained AUC feature constructed in this invention can stably characterize the degree of salt stress, is an effective proxy feature reflecting soil salinity, is not affected by the resolution of satellite data sources, and has good stability.

[0053] To make the present invention more applicable in different regions and application scenarios, several optional implementation methods are given below. These optional implementation methods are not mutually exclusive with Example 1 and can be combined as needed. (1) Adaptive integral window based on phenology: In addition to the fixed annual growing season window, the first derivative of smooth NDVI can be used to identify the turning point of greening and senescence, and the AUC calculation is limited to [SOS,EOS], thereby further reducing the influence of non-growing stages.

[0054] Dynamic threshold optimization for the target growing season: When determining the growing season window, the fixed thresholds T1 and T2 can be optimized into dynamic thresholds based on the statistical characteristics of NDVI sequences over many years.

[0055] (2) Land type grouping and labeling: When there are significant differences in land types in the study area, SSC-AUC models can be established for saline-alkali land, cultivated land, forest land and bare land respectively to improve the fitting consistency under different land types.

[0056] (3) Robustness and anomaly handling: Median filtering can be used to preprocess the spike anomalies in the NDVI sequence; quantile truncation or Winsorize can be used to process the extreme values ​​of AUC; and masking can be used to remove water bodies and permanent building areas.

[0057] (4) Uncertainty expression: The interval estimates of AUC and SSC can be formed based on the threshold traversal results. For example, the mean and standard deviation of SSC under several candidate thresholds can be output to characterize the parameter uncertainty. Alternatively, a pixel-level confidence index can be constructed based on the missing measurement ratio, smoothing residual and model residual.

[0058] (5) Multi-source fusion: While maintaining the consistency of the AUC definition, Sentinel can be used. 2. By fusing Landsat time series data with Landsat time series data on the time axis to form denser observations, the missing data can be reduced and the quality of curve reconstruction can be improved. Alternatively, weighted fusion can be used between multi-source results, with the weights determined by cloud coverage, number of observations and sensor noise assessment.

[0059] (6) Application extension: In addition to SSC, the threshold-constrained AUC framework of the present invention can also be extended to the inversion of other surface parameters that use vegetation process response as a proxy, such as the indirect monitoring of soil moisture stress, nutrient stress or heavy metal stress risk.

[0060] Example 3 This embodiment provides a soil salinity inversion system based on the area under the thresholded NDVI curve; like Figure 4 As shown, the system corresponding to this embodiment includes: a data acquisition module, a time series construction module, a smooth reconstruction module, a threshold AUC calculation module, a model calibration and inversion module, and a mapping output module, wherein: Data acquisition module: Accesses and manages Sentinel-2 remote sensing data and auxiliary data, and performs radiometric / atmospheric correction, geometric registration, cloud shadow removal and NDVI calculation; Time series construction module: Maps irregular observations into equally spaced sequences, and completes resampling, interpolation, and missing data quality labeling; Smoothing Reconstruction Module: Performs local least squares smoothing to generate smooth NDVI curves and can output derivative sequences and phenological inflection points; Threshold AUC Calculation Module: Calculates the AUC over a given threshold or a set of candidate thresholds, and outputs the threshold sensitivity curve and the optimal threshold. ; Model calibration and inversion module: Establishes the SSC-AUC function model and outputs the SSC raster at the pixel scale; Mapping output module: Generates SSC spatial distribution maps, hierarchical statistical maps, change detection maps, and uncertainty marker maps, and supports exporting to GeoTIFF, NetCDF, or raster database formats.

[0061] The system can be deployed on a local workstation, server, or cloud platform and can be integrated with common remote sensing processing toolchains. This embodiment also provides a computer-readable storage medium on which a computer program is stored. When the program runs on a processor, it executes steps S1–S5 of the above embodiment 1 to achieve automated SSC inversion and mapping.

[0062] The above are all preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Therefore, all equivalent changes made in accordance with the structure, shape and principle of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A method for soil salinity inversion based on the area under the thresholded NDVI curve, characterized in that, include: Multi-source optical satellite surface reflectance data of the target growing season in the study area were acquired. After quality control, the Normalized Difference Vegetation Index (NDVI) was calculated pixel by pixel to form a non-equal interval NDVI observation sequence at the pixel level. Time resampling and missing data interpolation are performed on non-equal interval NDVI observation sequences to generate equal interval NDVI sequences with a preset time step; Local least squares smoothing is used to denoise and reconstruct equally spaced NDVI sequences to obtain smoothed NDVI time series; Construct a candidate threshold set for NDVI, iterate through and calculate the area under the curve (AUC) of the smoothed NDVI sequence under each threshold, combine ground samples to determine the optimal threshold and generate a pixel-level optimal AUC raster for the study area. Based on the measured soil salinity content (SSC) at ground sampling points and the corresponding optimal AUC value, an SSC-AUC function model is constructed to obtain the spatial distribution results of SSC and output pixel-level uncertainty labeling information.

2. The method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 1, characterized in that, The calculation of the Normalized Difference Vegetation Index (NDVI) pixel by pixel after quality control includes: Acquire at least one optical satellite surface reflectance data covering the target growing season in the study area. After excluding months affected by snow cover in the study area, the surface reflectance data are standardized products that have undergone radiometric and atmospheric corrections. Pixel-level quality control is performed on the surface reflectance data to obtain quality-controlled non-equal interval effective observation data. Normalized Difference Vegetation Index (NDVI) was calculated pixel-by-pixel from the valid observation data. The NDVI values ​​of each pixel were then sorted chronologically by the satellite imaging date to form a non-equidistant NDVI observation sequence at the pixel level for the study area. ,in, For imaging time, The expression for the NDVI value corresponding to the imaging time is: ； In the formula, For near-infrared band reflectivity, This refers to the reflectivity in the red light band.

3. The method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 1, characterized in that, The generation of equally spaced NDVI sequences with a preset time step includes: A continuous time axis was constructed using the target growing season of the study area as the time range. A preset time step was set and equally spaced time nodes were generated. The preset time step was selected as a daily scale. The start and end times of the target growing season were calculated based on the historical NDVI time series of the study area to obtain the average NDVI sequence. The start date of the target growing season was the first time that the NDVI value in the multi-year average NDVI sequence continuously exceeded a preset threshold. The date on which the target growing season ends is the date on which the NDVI value last falls below a preset threshold for the last consecutive period. The date; For each pixel's non-equal interval NDVI observation sequence, time interpolation is used to estimate the NDVI value of the missing time nodes, realizing the mapping of non-equal interval observations to an equal interval time axis; For consecutive missing periods beyond the effective observation range, no extrapolation interpolation is performed. Simultaneously, the missing percentage and longest consecutive missing length for each pixel are statistically analyzed and used as quality markers. Finally, pixel-level NDVI sequences with pre-defined time steps and equal intervals for the study area are generated. ; Among them, threshold and The background NDVI values ​​were set based on the bare soil and typical vegetation in the study area.

4. The soil salinity inversion method based on the area under the thresholded NDVI curve according to claim 1, characterized in that, The method of using local least squares smoothing to denoise and reconstruct equally spaced NDVI sequences includes setting the sliding window length m to an odd number, and the polynomial order p being less than the sliding window length m. Based on the equally spaced NDVI sequences, a p-order polynomial is fitted within a sliding window of length m according to the least squares criterion. The polynomial fitting value at the center point of the sliding window is used as the smoothed output value. Full-sequence sliding fitting processing is then performed on the equally spaced NDVI sequences to obtain a smoothed NDVI time series. The expression for the smoothed NDVI time series is: ； in, The sliding window is half the window length. The convolution coefficients are determined by the sliding window length m and the polynomial order p. This represents the smoothed NDVI value at the k-th time node on an equally spaced time axis. This represents the original, equally spaced NDVI values ​​at the (k+j)th time node.

5. The method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 4, characterized in that, The method of using local least squares smoothing to denoise and reconstruct equally spaced NDVI sequences also includes setting multiple sets of conditions where m is an odd number and The (m,p) parameter combinations are used to perform local least squares smoothing on the equally spaced NDVI sequences. The curve shapes of the smoothed NDVI time series under different parameter combinations are compared, and the optimal (m,p) parameter combination is selected. The final smoothed NDVI time series is determined based on the smoothing result of the optimal (m,p) parameter combination.

6. The method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 1, characterized in that, The step of traversing and calculating the area under the curve (AUC) of the smoothed NDVI sequence at each threshold includes setting the range and step size of the NDVI threshold based on the ecological characteristics of vegetation growth in the study area, and constructing a set of candidate NDVI thresholds. For the smoothed NDVI time series, with candidate thresholds As constraints, calculate the conditions that must be met during the growing season. The area under the curve above the threshold for each time period is calculated using the trapezoidal rule to achieve numerical integration, yielding the pixel-level AUC value for each candidate threshold. The area under the threshold time integral is defined as: ； In the formula, , The start and end times of the target growing season in the study area. It is a function for maximizing the value.

7. A method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 6, characterized in that, The process of determining the optimal threshold by combining ground sampling points and generating the optimal AUC raster at the pixel level for the study area includes fitting the AUC values ​​of ground sampling point pixels under each candidate threshold with the measured soil salinity content (SSC) and calculating the coefficient of determination. The optimal threshold is determined by considering the root mean square error (RMSE). Based on the optimal threshold The area under the curve above the threshold is recalculated for the smoothed NDVI time series of all pixels in the study area to generate the optimal AUC raster at the pixel level for the study area. The expression for the optimal AUC raster is: ； in, Time interval The local optimal AUC value within, To preset the time step, , This represents the smoothed NDVI value for the corresponding time point.

8. The method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 1, characterized in that, The construction of the SSC-AUC function model includes obtaining ground sampling data of the 0-10cm topsoil layer in the study area and determining the conductivity of the soil leachate using the conductivity method. ,pass The measured SSC values ​​of ground sampling points were obtained by converting the SSC values ​​with the calibration relationship of SSC. The calibration relationship is based on the measured SSC values ​​of soil samples with different salinity gradients and their corresponding values. Values ​​were obtained through regression fitting; The optimal AUC value corresponding to each surface sample pixel is extracted. Using the measured SSC value of the sample points as the dependent variable and the optimal AUC value as the independent variable, an SSC-AUC function model is constructed. The model parameters are then calibrated using the nonlinear least squares method. The expression of the SSC-AUC function model is as follows: ； in, Optimal threshold The AUC value under the following conditions The baseline level of SSC under low salt stress conditions. Let AUC be the sensitivity coefficient to changes in salinity. This is the model bias term.

9. A method for soil salinity inversion based on the area under the thresholded NDVI curve according to claim 8, characterized in that, The process of obtaining the SSC spatial distribution result and outputting pixel-level uncertainty labeling information includes: The calibrated SSC-AUC function model was applied pixel by pixel to the optimal AUC raster at the pixel level in the study area. The predicted value of soil salinity content for each effective pixel was calculated, and the SSC spatial distribution raster of the study area was generated. At the same time, a classification map was generated from the SSC spatial distribution raster according to the soil salinization classification standard. A pixel-level confidence index is constructed by combining the missing measurement ratio of each pixel, the smoothing residual, and the model fitting residual. Pixels with abnormal NDVI curves are marked as low-confidence pixels. A pixel-level uncertainty marker map is generated and output. The SSC spatial distribution results, the grading map, and the uncertainty marker map are exported into a general raster data format.

10. A soil salinity inversion system based on the area under the thresholded NDVI curve, performed according to the method of claim 1, characterized in that, include: The data acquisition module is used to perform optical remote sensing image access, radiometric correction, atmospheric correction, quality control and NDVI calculation, and output pixel-level non-equidistant NDVI observation sequences; The time series construction module, connected to the data acquisition module, is used to perform time resampling, missing data interpolation and equal interval sequence generation, and outputs daily-scale equal interval NDVI sequences and missing data quality markers. The smoothing reconstruction module, connected to the timing construction module, is used to perform local least squares smoothing processing and output a smoothed NDVI sequence; The threshold AUC calculation module, connected to the smooth reconstruction module, is used to perform candidate threshold traversal, threshold constraint curve area under area calculation, and optimal threshold determination, and outputs the optimal threshold. and the corresponding AUC grid; The model calibration and inversion module is connected to the threshold AUC calculation module and is used to establish the SSC-AUC function model, calibrate the model parameters, and perform cell-level SSC inference. The mapping output module, connected to the model calibration and inversion module, is used to generate spatial distribution maps, hierarchical statistical maps, and uncertainty marker maps of soil salinity, and supports export in GeoTIFF or NetCDF format.

Images