Ground surface soil humidity multi-source data fusion and downscaling method and system

Abstract

The invention provides a surface soil humidity multi-source data fusion and downscaling method and system. The surface soil humidity multi-source data fusion and downscaling method comprises the following steps: performing quality evaluation and standardization processing on multi-source soil humidity data; carrying out anomaly detection and restoration by adopting a hierarchical spatial-temporal feature extractor; carrying out regional partitioning based on an adaptive boundary generation algorithm and a soft boundary transition model; data fusion and error control are realized through a regional cascade fusion network and an error chain tracking algorithm; and carrying out downscaling processing by utilizing a terrain adaptive downscaling framework and a multi-scale residual learning model to obtain a soil humidity data product with the resolution of 1km. The soil humidity data of complex terrains and high vegetation coverage areas can be effectively processed, the data fusion precision and downscaling effect are remarkably improved, and the technical problems of a traditional method in the aspects of data quality evaluation, anomaly detection, partition boundary processing, error control and the like are solved.

Description

Technical Field

[0001] The invention relates to a downscaling method, in particular to a method and system for fusion and downscaling of multi-source data of surface soil moisture. Background Art

[0002] As a key variable in terrestrial ecosystems, soil moisture plays an important role in the global energy cycle, water cycle and carbon cycle. Accurate soil moisture data is of great significance for climate change research, hydrological forecasting, agricultural production and ecosystem management. Especially in the context of global climate change, the spatiotemporal distribution characteristics of soil moisture have changed significantly, which directly affects the surface energy balance, vegetation growth and regional hydrological processes. Therefore, obtaining soil moisture data with high spatiotemporal resolution is of great scientific value for understanding regional ecohydrological processes and responding to climate change.

[0003] At present, soil moisture data mainly comes from ground station observations, remote sensing inversion and model simulation. Although ground station observations have high accuracy, their spatial representativeness is limited; passive microwave remote sensing (such as AMSR-E, SMOS, etc.) has a large spatial coverage, but its spatial resolution is usually 25-40 kilometers, which is difficult to meet the needs of regional scale applications; active microwave remote sensing (such as ASCAT) has a high temporal resolution, but the inversion accuracy in vegetation coverage areas is limited; land surface model simulation (such as CLM, NOAH, etc.) can provide continuous soil moisture data, but is subject to the uncertainty of input data and parameterization schemes. Traditional data fusion methods often use simple weighted average or Bayesian average models, and downscaling methods are mostly based on statistical regression or triangulation methods. These methods are not effective when dealing with complex terrain and high vegetation coverage areas.

[0004] The following technical problems exist in the current soil moisture data processing, including: the existing quality assessment methods are often based on a single indicator, which is difficult to fully reflect the quality characteristics of the data in terms of temporal continuity, spatial consistency and physical rationality; in addition, the outlier detection algorithm does not take environmental factors into consideration sufficiently, resulting in low anomaly recognition accuracy under extreme weather conditions; furthermore, the partitioning method mostly uses a fixed threshold, which is difficult to adapt to the spatial heterogeneity of terrain and vegetation cover, especially at the partition boundary, where mutations are prone to occur; in addition, the error transmission mechanism in the data fusion process has not been fully valued, and there is a lack of effective control over the error accumulation effect; in addition, the existing downscaling methods are difficult to maintain the spatial continuity and physical rationality of the data when dealing with transition zones and complex terrain areas, and lack effective assessment of the uncertainty of the downscaling results. These technical problems seriously restrict the accuracy and practicality of soil moisture data products. Summary of the invention

[0005] The purpose of the invention is to provide a method and system for fusion and downscaling of multi-source data of surface soil moisture, in order to solve the above problems.

[0006] Technical solution, according to one aspect of the present application, a method for fusion and downscaling of multi-source data of surface soil moisture includes the following steps:

[0007] Step S1, obtain six sets of original soil moisture data and unify the coordinate system and unit to obtain data sets Y1-Y6; calculate the data missing rate Mi and time continuity index Ti of data sets Y1-Y6 to generate the integrity index Ci; calculate the consistency index Li based on the data change rate Vi; use the physical constraint model to calculate the deviation Di between the theoretical value Pi and the actual value to generate the rationality index Ri; combine the integrity index Ci, consistency index Li, and rationality index Ri to obtain the quality index Qi; perform noise reduction and standardization on data sets Y1-Y6 to obtain data sets W1-W6 and their quality assessment matrix Q and reliability index R;

[0008] Step S2: Based on the data set W1-W6, the quality assessment matrix Q and the reliability index R, the time correlation matrix Tt and the space correlation matrix Ts are calculated, and the space-time feature set ST is generated by combining them; based on the space-time feature set ST, the time anomaly index At, the space anomaly index As and the environment anomaly index Ae are calculated to generate the anomaly labeling matrix A and the anomaly type matrix C; the data set W1-W6 is expanded into a tensor T(x, y, t), and the outliers are repaired according to the anomaly type matrix C to obtain the repaired data set W1'-W6';

[0009] Step S3, obtain the repair data set W1'-W6' and the pre-stored DEM data, calculate the slope S, slope aspect A and terrain undulation R of each spatial position, generate the terrain feature vector TF, and combine them to obtain the terrain feature set G; obtain NDVI and LAI data, combine with soil type data, calculate the vegetation cover index VI, generate the vegetation influence coefficient β, and form the vegetation feature set V; according to the terrain feature set G and the vegetation feature set V, calculate the terrain complexity and vegetation density, input the optimization function to obtain the optimal partition number N, and use the regional growing algorithm to generate the initial partition P; for each partition boundary in the initial partition P, extract the adjacent area features R1 and R2, calculate the transition zone feature B through the weight function ω, and generate the optimized partition D';

[0010] Step S4, obtain the repair data set W1'-W6', the optimized partition D' and the transition zone feature B, calculate the weight coefficient αi for the data in each partition, fuse the weighted data, and generate the initial fusion result F1; obtain the quality matrix Q, calculate the error transfer coefficient εj of each processing link, and build an error accumulation model; input the initial fusion result F1 into the error accumulation model, calculate the error estimation value e of each position; perform error compensation processing on the initial fusion result F1 according to the error estimation value e, and obtain the optimized fusion result F2;

[0011] Step S5, obtain the fused data set F2, terrain features G and vegetation features V, calculate the terrain correction coefficient K1 and the vegetation influence coefficient K2; input the fused data set F2, the terrain correction coefficient K1 and the vegetation influence coefficient K2 into the scale conversion function to generate preliminary 1km resolution data D1; obtain the pre-stored ground station measured data, and calculate the residual with the corresponding position data of the preliminary 1km resolution data D1; perform spatial distribution analysis on the residual, establish a residual correction model, correct the preliminary 1km resolution data D1, and obtain the final 1km resolution data product D2.

[0012] According to one aspect of the present application, a surface soil moisture multi-source data fusion and downscaling system is also provided, comprising:

[0013] at least one processor; and,

[0014] a memory communicatively connected to at least one of the processors; wherein,

[0015] The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the surface soil moisture multi-source data fusion and downscaling method described in any of the above technical solutions.

[0016] Beneficial effects: It can effectively process soil moisture data in complex terrain and high vegetation coverage areas, significantly improve data fusion accuracy and downscaling effect, and solve technical problems of traditional methods in data quality assessment, anomaly detection, partition boundary processing and error control. Related technical effects will be described in detail in conjunction with specific embodiments below. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a flow chart of the present invention.

[0018] Figure 2 It is a flow chart of step S1 of the present invention.

[0019] Figure 3 It is a flow chart of step S2 of the present invention.

[0020] Figure 4It is a flow chart of step S3 of the present invention.

[0021] Figure 5 It is a flow chart of step S4 of the present invention. DETAILED DESCRIPTION

[0022] According to one aspect of the present application, a method for fusion and downscaling of multi-source data of surface soil moisture is provided, comprising the following steps:

[0023] Step S1, obtain six sets of original soil moisture data and unify the coordinate system and unit to obtain data sets Y1-Y6; calculate the data missing rate Mi and time continuity index Ti of data sets Y1-Y6 to generate the integrity index Ci; calculate the consistency index Li based on the data change rate Vi; use the physical constraint model to calculate the deviation Di between the theoretical value Pi and the actual value to generate the rationality index Ri; combine the integrity index Ci, consistency index Li, and rationality index Ri to obtain the quality index Qi; perform noise reduction and standardization on data sets Y1-Y6 to obtain data sets W1-W6 and their quality assessment matrix Q and reliability index R;

[0024] Step S2: Based on the data set W1-W6, the quality assessment matrix Q and the reliability index R, the time correlation matrix Tt and the space correlation matrix Ts are calculated, and the space-time feature set ST is generated by combining them; based on the space-time feature set ST, the time anomaly index At, the space anomaly index As and the environment anomaly index Ae are calculated to generate the anomaly labeling matrix A and the anomaly type matrix C; the data set W1-W6 is expanded into a tensor T(x, y, t), and the outliers are repaired according to the anomaly type matrix C to obtain the repaired data set W1'-W6';

[0025] Step S3, obtain the repair data set W1'-W6' and the pre-stored DEM data, calculate the slope S, slope aspect A and terrain undulation R of each spatial position, generate the terrain feature vector TF, and combine them to obtain the terrain feature set G; obtain NDVI and LAI data, combine with soil type data, calculate the vegetation cover index VI, generate the vegetation influence coefficient β, and form the vegetation feature set V; according to the terrain feature set G and the vegetation feature set V, calculate the terrain complexity and vegetation density, input the optimization function to obtain the optimal partition number N, and use the regional growing algorithm to generate the initial partition P; for each partition boundary in the initial partition P, extract the adjacent area features R1 and R2, calculate the transition zone feature B through the weight function ω, and generate the optimized partition D';

[0026] Step S4, obtain the repair data set W1'-W6', the optimized partition D' and the transition zone feature B, calculate the weight coefficient αi for the data in each partition, fuse the weighted data, and generate the initial fusion result F1; obtain the quality matrix Q, calculate the error transfer coefficient εj of each processing link, and build an error accumulation model; input the initial fusion result F1 into the error accumulation model, calculate the error estimation value e of each position; perform error compensation processing on the initial fusion result F1 according to the error estimation value e, and obtain the optimized fusion result F2;

[0027] Step S5, obtain the fused data set F2, terrain features G and vegetation features V, calculate the terrain correction coefficient K1 and the vegetation influence coefficient K2; input the fused data set F2, the terrain correction coefficient K1 and the vegetation influence coefficient K2 into the scale conversion function to generate preliminary 1km resolution data D1; obtain the pre-stored ground station measured data, and calculate the residual with the corresponding position data of the preliminary 1km resolution data D1; perform spatial distribution analysis on the residual, establish a residual correction model, correct the preliminary 1km resolution data D1, and obtain the final 1km resolution data product D2.

[0028] By building a complete multi-source soil moisture data processing system, the whole process of intelligent processing from data quality assessment to the final high-resolution product generation is realized. Specifically, in the data preprocessing stage, the three-level quality assessment framework (TQL) and the distribution adaptive normalization method (DANM) are used to accurately evaluate and standardize data quality; in the spatiotemporal analysis stage, the hierarchical spatiotemporal feature extractor (HSFE) and the environment-aware anomaly detector (EADD) are used to accurately identify and repair outliers; in the regional partitioning stage, the adaptive boundary generation algorithm (ABGA) and the soft boundary transition model (SBTM) are used to solve the partitioning problem of complex terrain areas; in the data fusion stage, the regional cascade fusion network and the error chain tracking algorithm (ECTA) are used to achieve high-precision data fusion and error control; in the downscaling stage, the terrain adaptive downscaling framework and the multiscale residual learning model (MRLM) are used to generate high-quality fine-scale data. This full-process intelligent processing method shows significant advantages in processing soil moisture data under complex terrain, high vegetation coverage and extreme weather conditions, which not only improves the automation level of data processing, but also significantly improves the accuracy and reliability of the final product. Especially in climate change research and hydrological model applications, the high-resolution soil moisture data generated by this method can more accurately reflect the characteristics of regional moisture changes and provide more reliable data support for ecological environment monitoring and agricultural production management.

[0029] like Figure 2 As shown, according to one aspect of the present application, step S1 is specifically:

[0030] Step S11, obtain six sets of original soil moisture data CCI, CLM, ERA_Interim, MOS, NOAH and VIC, and record them as data sets X1-X6 in sequence, read their latitude and longitude information and timestamp information; convert the projection coordinate system of X1-X6 to WGS84, generate a unified coordinate grid; unify the soil moisture values ​​in the data sets X1-X6 into m³ / m³, and generate standardized data sets Y1-Y6;

[0031] Step S12, according to the standardized data sets Y1-Y6, count the total number of pixels total_pixels and the number of missing pixels missing_pixels of each data set, and calculate the missing rate Mi; obtain the time series of the standardized data sets Y1-Y6, calculate the continuous observation days continuous_days and the total observation days total_days, and obtain the time continuity index Ti; input the missing rate Mi and the time continuity index Ti into the weighting function f(x) to generate the integrity index Ci; based on the time series of the standardized data sets Y1-Y6, calculate the data difference between adjacent moments divided by the time interval to obtain the change rate Vi; construct the data fluctuation range according to the change rate Vi and generate the envelope Ei; input the envelope Vi and the envelope Ei into the function h(x) to obtain the consistency index Li; obtain precipitation, evaporation and terrain data, input the physical constraint model k(x) to obtain the theoretical soil moisture value Pi; calculate the absolute difference between the standardized data sets Y1-Y6 and the theoretical soil moisture value Pi to obtain the deviation Di; input the deviation Di into the function m(x) to generate the rationality index Ri;

[0032] Step S13, construct sliding windows of different scales {t1, t2, ..., tn} for the time series of the standardized data sets Y1-Y6; calculate the data variance σ²(x, y, t) in each time window at each spatial position (x, y); select the time window ts with the smallest variance and perform weighted average filtering on the data at this position; repeat this process for all spatial positions to obtain the denoised data sets Z1-Z6;

[0033] Step S14: Calculate the probability distribution features Pi(x) of each of the denoised data sets Z1-Z6; construct a reference probability distribution Pr(x); map the denoised data sets Z1-Z6 to the reference distribution space through the cumulative distribution function Fi(x); apply the reverse mapping function F⁻¹r to generate the standardized data sets W1-W6;

[0034] Step S15, based on the quality index Qi, calculate the time dimension reliability index Rt for the time series of the standardized data set W1-W6; calculate the spatial dimension reliability index Rs for the spatial distribution of the standardized data set W1-W6; input Rt and Rs into the comprehensive evaluation function φ(x) to generate the final quality evaluation matrix Q and reliability index set R.

[0035] By introducing the three-level quality assessment framework (TQL) and the distributed adaptive normalization method (DANM), intelligent quality assessment and standardized processing of multi-source soil moisture data are realized. Specifically, in the data integrity assessment layer, by calculating the combination of the missing rate Mi and the temporal continuity index Ti, not only the data missing situation can be identified, but also the temporal continuity characteristics of the data can be evaluated; in the data consistency assessment layer, by calculating the data change rate Vi and the fluctuation envelope Ei, the abnormal fluctuation and mutation characteristics of the data can be effectively captured. Compared with the traditional single threshold method, the data that does not conform to the law of natural change can be more accurately identified; in the physical rationality assessment layer, by establishing a physical constraint model considering precipitation, evaporation and terrain factors, the quantitative evaluation of the physical desirability of the data is realized. This multi-level evaluation system can not only comprehensively evaluate the data quality, but also provide a reliable weight basis for subsequent data fusion. In addition, the adaptive multi-scale sliding window filter (AMWF) is used for data denoising. This method can adaptively adjust the time window size according to the local variance of the data, which not only ensures the smoothness of the data, but also avoids the information loss caused by over-smoothing. Finally, data standardization was achieved through the distribution adaptive normalization method. This method not only maintained the relative distribution characteristics of the data through probability distribution mapping, but also achieved the unification of the scales of data from different sources, laying the foundation for subsequent data fusion.

[0036] like Figure 3 As shown, according to one aspect of the present application, step S2 is specifically:

[0037] Step S21, according to the standardized data sets W1-W6 and the quality matrix Q outputted in step S1, set the time delay sequence τ={1,2,...,n}; for any two data sets Wi and Wj, calculate the Pearson correlation coefficient at each τ to obtain the time correlation matrix Tt(i,j); for each data set Wi, set the spatial displacements Δx and Δy, calculate the correlation coefficient between the position (x,y) and (x+Δx,y+Δy), and generate the spatial correlation matrix Ts(x,y); combine Tt and Ts to generate the spatiotemporal correlation feature set ST;

[0038] Step S22, using the spatiotemporal correlation feature set ST and the reliability index R, calculate the mean μt and standard deviation σt in the local time window at each time point t; calculate the difference between the data value Wi(t) and μt, and divide it by σt to obtain the time anomaly index At; calculate the mean μs and standard deviation σs of the local area for each spatial position (x, y), calculate the difference between the data value Wi(x, y) and μs, and divide it by σs to obtain the spatial anomaly index As; obtain the precipitation and temperature data of the corresponding position and time, input the environmental anomaly function f(x) to calculate the environmental anomaly index Ae; based on the combination of the temporal anomaly index At, the spatial anomaly index As and the environmental anomaly index Ae, generate the anomaly labeling matrix A, and generate the anomaly type matrix C according to the dominant anomaly type;

[0039] Step S23, obtain the standardized data set W1-W6, the outlier labeling matrix A and the outlier type matrix C, and construct a three-dimensional tensor T(x, y, t); for the data points marked as abnormal, select the time interpolation, space interpolation or time-space joint interpolation strategy according to the anomaly type in the anomaly type matrix C; for time anomalies, use the weighted average of adjacent time points for repair; for space anomalies, use local space interpolation for repair; for environmental anomalies, combine space-time information and environmental factors for comprehensive repair; after completing the repair of all abnormal points, obtain the repaired data set W1'-W6'.

[0040] By introducing a hierarchical spatiotemporal feature extractor (HSFE) and an environmentally aware anomaly detector (EADD), high-precision anomaly detection and repair of soil moisture data are achieved. Specifically, in the spatiotemporal correlation analysis, a complete spatiotemporal feature set ST is constructed by calculating the correlation coefficient matrix Tt under different time delays τ and the correlation matrix Ts under different spatial displacements. This method can not only capture the periodic change characteristics of the data, but also identify the abnormal distribution pattern in space. In the anomaly detection stage, the multi-dimensional identification of outliers is achieved by calculating the temporal anomaly index At and the spatial anomaly index As, combined with the environmental anomaly index Ae constructed by environmental factors. Compared with the traditional single threshold detection method, this comprehensive detection method considering environmental factors significantly improves the detection accuracy, especially in extreme weather events or complex terrain areas, and can more accurately distinguish between real outliers and natural fluctuations. In the data repair stage, by constructing a third-order tensor T(x,y,t) and an adaptive interpolation strategy based on the anomaly type, the spatiotemporal continuity of the data is maintained, and the over-smoothing of outliers by traditional interpolation methods is avoided. Especially in areas with high vegetation coverage, this method can better maintain the natural variation characteristics of soil moisture.

[0041] like Figure 4 As shown, according to one aspect of the present application, step S3 is specifically:

[0042] Step S31, obtain pre-stored DEM data, for each grid position (x, y), calculate the elevation difference between the central pixel and the adjacent pixel based on the 3×3 neighborhood, and obtain the slope data S; calculate the elevation gradient of each pixel in the four directions of east, west, south and north, and obtain the aspect data A; calculate the ratio of the maximum elevation difference to the horizontal distance in the 5×5 neighborhood, and obtain the terrain relief R; combine the slope data S, the aspect data A and the terrain relief R into a terrain feature vector TF; perform cluster analysis on the TF of all spatial positions to generate a terrain feature set G;

[0043] Step S32, obtaining pre-stored NDVI data and LAI data, and combining with pre-stored soil type data to calculate the vegetation coverage of each location; performing time series filtering on the NDVI value to remove abnormal fluctuations; calculating the thickness of the vegetation layer according to the LAI data; inputting the NDVI, LAI and soil type data into the vegetation impact model to calculate the vegetation water interception coefficient; calculating the vegetation impact coefficient β based on the interception coefficient and coverage; combining all vegetation related parameters to generate a vegetation feature set V;

[0044] Step S33, obtaining a terrain feature set G, a vegetation feature set V, and a restored data set W1'-W6'; calculating the standard deviation of the terrain features in the study area to obtain a terrain complexity index; statistically analyzing the spatial distribution density of vegetation features to obtain a vegetation density index; inputting the terrain complexity and vegetation density into the optimization function to calculate the optimal partition number N; selecting N initial seed points, performing regional growth based on feature similarity, and generating an initial partition result P, where N is a positive integer;

[0045] Step S34, obtain the initial partition P of step S33, and extract each partition boundary line; for each boundary line, obtain the terrain characteristics R1 and soil moisture characteristics R2 of the adjacent areas on both sides; calculate the weight coefficients ω1 and ω2, input the terrain characteristics R1 and the soil moisture characteristics R2 into the weighting function, and obtain the boundary transition characteristics; expand a specific width to both sides along the boundary to construct a transition zone, and record the characteristic value B of each point in the transition zone; combine the original partition P with the transition zone feature B to generate an optimized partition D' and a complete transition zone feature set B.

[0046] By introducing the adaptive boundary generation algorithm (ABGA) and the soft boundary transition model (SBTM), the intelligent partitioning and boundary optimization of the study area are realized. Specifically, in the stage of terrain feature extraction, the complete terrain feature vector TF is constructed by calculating the slope S, the slope aspect A and the terrain relief R. This method can fully reflect the complexity of the terrain. In the vegetation impact assessment, the vegetation-soil moisture interaction model is established and the vegetation impact coefficient β is calculated, which effectively quantifies the impact of vegetation on soil moisture. In the process of dynamic partition generation, the optimal number of partitions N is adaptively determined by calculating the terrain complexity and vegetation density, avoiding the limitations of the traditional fixed number of partitions method. Especially in terms of boundary processing, the transition function B(x,y) constructed by the soft boundary transition model realizes a smooth transition between partitions, effectively solving the mutation problem at the boundary of the traditional hard partitioning method. This partitioning method that comprehensively considers the characteristics of terrain, vegetation and soil moisture shows obvious advantages in complex terrain areas and areas with significant changes in vegetation coverage. It can more accurately reflect the spatial distribution characteristics of soil moisture and provide a reliable spatial reference framework for subsequent data fusion and downscaling.

[0047] like Figure 5 As shown, according to one aspect of the present application, step S4 is specifically:

[0048] Step S41, obtain the repaired data set W1'-W6', partition D' and transition zone feature B; for each partition, extract the time series features and spatial distribution features of each data source in the partition; calculate the similarity between each data source and the overall feature of the region to obtain the initial weight; adjust the initial weight based on the quality matrix Q to generate the final weight coefficient αi; apply the weight coefficient αi to the data in each partition for weighted averaging; in the transition zone area, calculate the dynamic weight according to the transition zone feature B, and perform a smooth transition on the fusion results of adjacent partitions; combine the fusion results of all partitions to generate the initial fusion result data set F1;

[0049] Step S42, obtain the initial fusion result data set F1 and the quality matrix Q; calculate the error transfer coefficient εj for each link in the data processing process, including standardization, anomaly detection, partition optimization and fusion; construct an error accumulation function based on the error transfer coefficient, and analyze the error transmission path in the processing process; for each data point in the initial fusion result data set F1, calculate its cumulative error estimate e; design an error compensation strategy based on the size and distribution characteristics of the error estimate e; apply the compensation strategy to the initial fusion result data set F1 to obtain the optimized fusion result F2; in the compensation process, adopt local re-fusion for high error areas, and keep the original fusion results for low error areas; integrate the compensated results to generate the final fusion data set F2.

[0050] By introducing the regional cascade fusion network and the error chain tracking algorithm (ECTA), high-precision fusion and error control of multi-source data are achieved. Specifically, in the regional adaptive fusion stage, the dynamic adjustment of data fusion weights is achieved by constructing the time dimension weight matrix Wt and the space dimension weight matrix Ws. This method can adaptively allocate fusion weights according to the spatiotemporal characteristics and quality indicators of the data, significantly improving the accuracy of the fusion results. In terms of transition zone processing, by constructing the distance decay function d(x) and the dynamic change coefficient γ, a smooth transition between partitions is achieved, effectively avoiding the discontinuity problem at the partition boundary of the traditional fusion method. In terms of error transmission control, by constructing the error transmission graph G(V,E) and the error accumulation function A(p), the whole process of error tracking and control is achieved. This method can not only identify the main source of error, but also effectively compensate for the error, especially in areas with large differences in data quality. This method can significantly improve the reliability of the fusion results. In addition, through the adaptive adjustment of the error compensation kernel function K(e), differentiated processing of errors of different degrees is achieved, avoiding the limitations of the traditional unified compensation method.

[0051] According to one aspect of the present application, step S5 is specifically:

[0052] Step S51, obtain the final fused data set F2, terrain features G and vegetation features V; calculate the terrain relief, slope change rate and slope aspect variability in each 1km grid, and generate the terrain correction coefficient K1; count the vegetation coverage change and vegetation type distribution in each 1km grid, and calculate the vegetation influence coefficient K2; construct a spatial interpolation grid from 25km to 1km; for each 25km pixel, extract the terrain correction coefficient K1 and vegetation influence coefficient K2 within its corresponding range, and establish a local terrain-vegetation weight model; distribute the values ​​of the final fused data set F2 to the 1km grid through the local weight model to generate preliminary high-resolution data D1;

[0053] Step S52, obtaining the preliminary downscaling result D1 and the pre-stored measured data of the ground site; for each measured site location, extracting the data of the corresponding location in the preliminary downscaling result D1 and the 3×3 grids around it; calculating the difference between the measured value and the preliminary downscaling result D1 to obtain a residual sequence; performing spatial autocorrelation analysis on the residual sequence, and calculating the spatial distribution characteristics of the residuals; based on the spatial distribution characteristics of the residuals, constructing a Kriging interpolation model to generate a residual distribution map for the entire region; correcting the preliminary downscaling result D1 according to the residual distribution map, and applying the corresponding residual correction value to each 1km grid; combining the corrected data to generate the final 1km resolution soil moisture data product D2; and generating a correction coefficient and uncertainty estimate for each grid point as a quality assessment indicator.

[0054] By introducing a terrain adaptive downscaling framework and a multiscale residual learning model (MRLM), high-precision downscaling and correction of soil moisture data are achieved. Specifically, in the multi-factor downscaling model, by constructing the terrain correction factor K1 and the vegetation influence factor K2, the precise quantification of the terrain and vegetation effects is achieved. This method can effectively handle the scale conversion problem under complex terrain conditions. In the residual analysis stage, by constructing the local weighted regression equation and the spatial variation function γ(h), the accurate modeling of the residuals and the capture of the spatial distribution characteristics are achieved. Compared with the traditional linear interpolation method, the spatial variation characteristics of the residuals can be more accurately described. In the correction process, the nonlinear correction of the residuals is achieved through the dynamic adjustment of the adaptive correction function f(rc). This method can adaptively adjust the correction strength according to the residual amplitude and reliability, avoiding the limitations of the traditional fixed correction coefficient method. Especially in mountainous areas and areas with complex vegetation coverage, this downscaling method that considers the influence of multiple factors shows significant advantages, not only improving the spatial resolution of the downscaling results, but also maintaining the physical rationality of the data.

[0055] According to one aspect of the present application, step S15 is specifically:

[0056] Step S151, obtain the quality index Qi, extract the time series data for each data set W1-W6; calculate the autocorrelation coefficient rt(k) of each time series, where k is the time lag; construct the time continuity feature vector Vt based on the autocorrelation coefficient rt(k) and apply principal component analysis to extract the main time change mode Mt; compare the main time change mode Mt with the preset time change threshold, and calculate the time dimension deviation δt; calculate the time dimension reliability index Rt based on the time dimension deviation δt and the quality index Qi;

[0057] Step S152, obtain the standardized data set W1-W6, and construct a local area matrix L(x, y) for each spatial position; calculate the spatial autocorrelation coefficient rs(h) in the local area matrix L(x, y), where h is the spatial distance; construct a spatial continuity feature vector Vs based on the spatial autocorrelation coefficient rs(h); perform a variogram analysis on the spatial continuity feature vector Vs to obtain the spatial structure feature Gs; compare the spatial structure feature Gs with the regional geographical features to calculate the spatial dimension deviation δs; calculate the spatial dimension reliability index Rs based on the spatial dimension deviation δs and the quality index Qi;

[0058] Step S153, obtain the time dimension reliability index Rt and the space dimension reliability index Rs; construct the spatiotemporal weight matrix W(t,s); for each data point, extract the time dimension reliability index Rt and the space dimension reliability index Rs values ​​in its time and space neighborhood; apply the spatiotemporal weight matrix W(t,s) to perform weighted combination of the time dimension reliability index Rt and the space dimension reliability index Rs to generate an initial comprehensive index I; apply local smoothing filtering to the initial comprehensive index I to eliminate mutation points; standardize the filtered results to generate the final quality assessment matrix Q; calculate the reliability level of each position based on the quality assessment matrix Q to form a reliability index set R.

[0059] By introducing a joint evaluation mechanism of time and space dimensions, the accurate quantification of soil moisture data quality is achieved. Specifically, in the time dimension evaluation, by calculating the autocorrelation coefficient rt(k) and constructing the time continuity feature vector Vt, the periodic change characteristics and mutation points of the data can be effectively captured. Compared with the traditional single time series analysis method, it can more comprehensively reflect the time evolution characteristics of the data; in the spatial dimension evaluation, by constructing the local area matrix L(x,y) and calculating the spatial autocorrelation coefficient rs(h), it can not only identify spatial anomalies, but also evaluate the spatial continuity of the data. Especially in areas with complex terrain, this method can more accurately reflect the spatial distribution characteristics of the data. In the final comprehensive evaluation, the dynamic weight allocation of the time and space dimensions is realized by constructing the spatiotemporal weight matrix W(t,s). This method can adaptively adjust the evaluation weight according to the local characteristics of the data, avoiding the limitations of the traditional fixed weight method. At the same time, by performing local smoothing filtering on the comprehensive index I, the noise in the evaluation results is effectively eliminated, and the stability and reliability of the evaluation results are improved. This multi-dimensional quality assessment method can not only comprehensively reflect the quality characteristics of the data, but also provide a reliable weight basis for subsequent data fusion. Especially in climate change research, it can more accurately assess the credibility of data from different sources.

[0060] According to one aspect of the present application, step S23 is specifically:

[0061] Step S231, obtain data sets W1-W6, construct a time-space-observation value tensor T(x, y, t); calculate the sliding window statistics St of the time-space-observation value tensor T in the time dimension, including the mean, standard deviation and trend coefficient; calculate the local statistics Ss of the time-space-observation value tensor T in the spatial dimension, including spatial gradient and curvature characteristics; combine the sliding window statistics St and the local statistics Ss to form a spatiotemporal feature tensor TSF; construct a local correlation matrix LC based on the spatiotemporal feature tensor TSF for subsequent anomaly repair;

[0062] Step S232, obtain the anomaly type matrix C, extract the time anomaly point set At; for each anomaly point in the time anomaly point set At, obtain its TSF feature vector and LC correlation data; construct an adaptive time weight function wt, and the weight coefficient decays with the time distance and correlation; based on the adaptive time weight function wt, weight the effective values ​​of the n moments before and after the anomaly point to generate a time dimension repair value Rt; apply physical constraints to the repair result to make the repair value within a reasonable range;

[0063] Step S233, obtain a set of spatial outlier points As; for each outlier point, construct a spatial influence domain D(x, y); extract terrain gradient g and vegetation coverage v in the spatial influence domain D(x, y); construct a spatial weight function ws based on the terrain gradient g and vegetation coverage v, and the weight varies with the terrain similarity and vegetation consistency; apply the spatial weight function ws to weight the valid values ​​in the spatial influence domain D(x, y) to generate a spatial dimension repair value Rs; perform a consistency check on the spatial dimension repair value Rs and the surrounding grid points to ensure spatial continuity;

[0064] Step S234, obtain the time repair value Rt and the space repair value Rs; construct the repair optimization function f(Rt, Rs), and calculate the credibility index of the time and space repair based on the TSF feature tensor; select the optimal repair value according to the credibility index, or perform weighted combination according to the credibility index; perform localized iterative optimization on the repair result until the first-order and second-order derivatives of adjacent time and space positions meet the continuity condition; reorganize the repaired data into the original data format to generate the repaired data set W1'-W6'.

[0065] By introducing the spatiotemporal feature tensor TSF and the adaptive repair strategy, high-precision identification and repair of outliers are achieved. Specifically, in the tensor construction stage, the complete spatiotemporal feature tensor TSF is constructed by calculating the sliding window statistics St in the time dimension and the local statistics Ss in the space dimension. This method can simultaneously capture the temporal continuity and spatial correlation of the data. In the temporal anomaly repair, the accurate repair of outliers is achieved through the dynamic adjustment of the adaptive temporal weight function wt. This method can adaptively assign weights according to the temporal correlation of the data, avoiding the limitations of the traditional fixed weight method. In terms of spatial anomaly repair, by constructing the spatial influence domain D(x,y) and the spatial weight function ws based on the terrain gradient g and the vegetation coverage v, spatial interpolation considering the influence of terrain and vegetation is achieved. This method is particularly suitable for data repair in complex terrain and high vegetation coverage areas. In the repair optimization process, the intelligent selection of the optimal repair value is achieved through the repair optimization function f(Rt,Rs) and the credibility evaluation based on the TSF feature, and the spatiotemporal continuity of the repair result is ensured through localized iterative optimization. This repair method that comprehensively considers multi-dimensional features significantly improves the accuracy and reliability of outlier repair.

[0066] According to one aspect of the present application, step S33 is specifically:

[0067] Step S331, obtain the terrain feature set G, calculate the coefficients of variation cv1 and cv2 of the slope and slope aspect in the study area; obtain the vegetation feature set V, calculate the spatial heterogeneity index h1 of the vegetation index; construct the complexity evaluation function f(cv1, cv2, h1), generate the regional complexity matrix CM; apply wavelet transform to the regional complexity matrix CM, extract the characteristic frequencies at different scales, and obtain the scale characteristic vector SF; calculate the optimal partition number estimate N1 based on the scale characteristic vector SF;

[0068] Step S332, obtain the complexity matrix CM, construct the spatial distance matrix D; calculate the local density ρ of each location point and the minimum distance δ to the high-density point; identify the density peak point on the ρ-δ graph, and select the first N1 points as the initial seed point set S; extract the feature vector FV of each seed point, the feature vector FV contains the terrain, vegetation and soil moisture characteristics; calculate the Mahalanobis distance between the seed points, and generate the seed point optimization set S'; N1 is a preset natural number;

[0069] Step S333, obtain the optimized seed point set S', construct the feature space distance function d(x); for each non-seed point p, calculate its feature distance d(p, s) with all seed points; set a dynamic growth threshold θ, the initial value of which is determined based on the local feature variance; when the feature distance d(p, s) is less than θ, assign the non-seed point p to the nearest seed point; θ is adaptively adjusted along with the regional growth process, and the adjustment amplitude is inversely proportional to the feature variance in the current partition; generate the initial partition result P1;

[0070] Step S334, obtain the initial partition P1, and construct a partition boundary set B; for each boundary line in the partition boundary set B, extract the feature sequence of the strip areas on both sides thereof; calculate the Jensen-Shannon divergence JS of the feature sequence as a boundary significance index; stratify the boundaries based on the Jensen-Shannon divergence JS value and the set JS threshold, keep the high JS value boundaries sharp, and smooth the low JS value boundaries; perform a region merging test inside the partition, and merge adjacent partitions when the feature similarity exceeds the threshold; generate the final partition result P.

[0071] By innovatively introducing the complexity evaluation function and density peak clustering method, the intelligent partitioning and optimization of the study area are realized. Specifically, in the complexity evaluation stage, by calculating the coefficients of variation cv1 and cv2 of slope and aspect and the spatial heterogeneity index h1 of vegetation index, a complete complexity evaluation function f(cv1,cv2,h1) is constructed. This method can fully reflect the terrain and vegetation complexity of the region. In the seed point selection process, the automatic identification of the optimal seed point is realized by calculating the local density ρ and the density peak spacing δ, avoiding the subjectivity of traditional manual selection of seed points. In the regional growth stage, the precise division of the partition boundary is achieved through the adaptive adjustment of the dynamic growth threshold θ. This method can automatically adjust the growth range according to the local feature variance, which is particularly suitable for the processing of feature transition zones. In the boundary optimization link, by calculating the Jensen-Shannon divergence JS of the feature sequence and the boundary significance stratification, the intelligent optimization of the partition boundary is realized, which effectively solves the mutation problem of the traditional partition method at the boundary and makes the partition result more in line with the actual geographical characteristics.

[0072] According to one aspect of the present application, step S41 is specifically:

[0073] Step S411, obtain the data set W1'-W6' output by step S2 and the partition data D' output by step S3. Extract the time series characteristics of the data in each partition Di and calculate the Fourier transform coefficient F(ω). Apply singular spectrum analysis to F(ω) to extract the main periodic component PC. Calculate the projection coefficient α of each data source on PC. Construct the time dimension weight matrix Wt based on the α value. Calculate the consistency index Ct of each data source with the overall time change pattern of the region.

[0074] Step S412: Obtain the time dimension weight matrix Wt output in step S411, and extract the spatial distribution characteristics in each partition Di. Construct a generalized form of the local spatial autocorrelation index Moran's I, and calculate the spatial structure characteristics of different directions and scales. Calculate the matching degree β with the typical spatial structure of the region for each data source. Construct the spatial dimension weight matrix Ws based on the β value. Combine the spatial consistency index Cs of the data source with the time consistency index Ct of step S411 to generate the initial comprehensive weight W0.

[0075] Step S413, obtain the initial comprehensive weight W0 of step S412 and the quality matrix Q output by step S1. Construct a weight optimization function L(W), and the objective function includes data quality terms, spatial smoothness terms, and time continuity terms. For each partition Di, extract the local quality index qi. Based on qi, modify W0 to generate a weight matrix W1 that considers quality factors. Apply the Lagrange multiplier method to solve the optimal solution of L(W) and obtain the final weight coefficient αi.

[0076] Step S414, obtain the final weight coefficient αi of step S413 and the transition zone feature B output by step S3. Apply αi to the data in each partition Di for weighted fusion to obtain the partition fusion result Fi. Construct the distance attenuation function d(x) at the transition zone position. Calculate the dynamic change coefficient γ of the adjacent partition weights based on d(x). For the data in the transition zone area, use the bilateral filtering algorithm considering γ for smoothing. Combine the fusion results of all partitions and the transition zone processing results to generate a complete fusion data set F1.

[0077] By introducing a multi-dimensional weight optimization system and a dynamic transition zone processing mechanism, high-precision fusion of partitioned data is achieved. Specifically, in the temporal feature extraction stage, the periodic characteristics of the data are accurately quantified by applying singular spectrum analysis to the Fourier transform coefficient F(ω) and calculating the projection coefficient α. This method can effectively identify the change patterns at different time scales. In the spatial feature evaluation, the multi-directional and multi-scale spatial structural feature analysis is achieved by constructing the generalized Moran's I index, which can capture spatial correlation more comprehensively than the traditional single spatial statistical method. In the weight optimization process, the global optimization of the weight is achieved by constructing an optimization function L(W) containing data quality terms, spatial smoothing terms, and temporal continuity terms. This method not only considers the impact of data quality, but also ensures the spatiotemporal continuity of the fusion results. Especially in the processing of transition zones, by constructing the distance attenuation function d(x) and the dynamic change coefficient γ, a smooth transition between adjacent partitions is achieved, which effectively solves the discontinuity problem of the traditional fusion method at the partition boundary. This multi-dimensional fusion optimization method shows obvious advantages in processing areas with complex terrain and significant changes in vegetation coverage, and significantly improves the spatial consistency and temporal continuity of the fusion results.

[0078] According to one aspect of the present application, step S42 is specifically:

[0079] Step S421, obtain the intermediate results of each processing link from step S1 to S41 and the quality matrix Q output from step S1. Construct an error transfer graph G(V,E), where the node V represents the processing link and the edge E represents the data flow. For each processing link, extract the quality index qi of its input and output data. Calculate the error change rate δe between adjacent processing links. Construct an error transfer matrix M based on δe, where the matrix element mij represents the error transfer coefficient εj from link i to link j.

[0080] Step S422, obtain the error transfer matrix M output by step S421 and the fusion result F1 output by step S41. For each spatial position p(x,y), construct the error accumulation function A(p). Extract the error transfer coefficient εj of all upstream processing links at point p. Calculate the cumulative effect λt of the error in the time dimension and the cumulative effect λs in the spatial dimension. Construct an error accumulation model based on λt and λs to obtain the error estimate e(p) at each position. Generate a global error distribution map E.

[0081] Step S423, obtain the error distribution map E output by step S422, and perform hierarchical clustering on the errors. Calculate the spatial autocorrelation feature r(h) of each error level, where h is the spatial lag distance. Construct the error compensation kernel function K(e), and the kernel function parameters are dynamically adjusted with the error amplitude and spatial distribution characteristics. For high error areas, extract local area data to construct tensor T, and apply tensor decomposition to obtain the main error mode. Generate the error compensation strategy matrix C.

[0082] Step S424, obtain the error compensation strategy matrix C of step S423 and the fusion result F1 of step S41. For high error areas, extract data in the local spatiotemporal window. Construct a local weighted reconstruction model based on heteroscedastic kernel. Calculate the difference δ between the reconstructed data and the original fusion result. When δ exceeds the threshold, replace the original result with the reconstructed result; when δ is small, correct it by the error compensation kernel function K(e). For low error areas, keep the original fusion result. Combine the processing results of all areas to generate the final fusion data set F2.

[0083] By introducing the error transfer graph G(V,E) and the adaptive error compensation mechanism, high-precision error control in the fusion process is achieved. Specifically, in the error transfer analysis stage, by constructing a complete error transfer graph and calculating the error change rate δe, the error transfer path in the processing flow is accurately characterized. This method can effectively identify the main sources and accumulation nodes of the error. In the error accumulation modeling, by constructing an error accumulation model that considers the cumulative effect λt of the time dimension and the cumulative effect λs of the space dimension, all-round tracking of the error is achieved. Compared with the traditional single error evaluation method, the spatial distribution and time evolution of the error can be estimated more accurately. In the error compensation stage, a differentiated error compensation strategy is implemented through the adaptive kernel function K(e) based on the error amplitude and spatial distribution characteristics. In particular, for high error areas, the main error mode is extracted by tensor decomposition, and the local weighted reconstruction model based on the heteroscedastic kernel is used for accurate correction. This full-process error control method not only significantly improves the reliability of the fusion results, but also provides an important basis for the continuous improvement of data quality. In particular, when the quality of multi-source data varies greatly, the accuracy and stability of the fusion results can be better guaranteed.

[0084] According to one aspect of the present application, step S52 is specifically:

[0085] Step S521, obtain the preliminary downscaling result D1 output by step S51 and the pre-stored ground site measured data S. For each site location (x, y), extract the corresponding location in D1 and its 3×3 neighborhood data matrix N(x, y). Calculate the difference between the measured value s(x, y) and the center value of N(x, y) to obtain the initial residual r(x, y). Construct a local weighted regression equation and calculate the partial correlation coefficient ρi between the residual and the terrain factor and vegetation factor. Based on ρi, screen the set of significant influencing factors F and generate the residual feature vector RV.

[0086] Step S522, obtain the residual feature vector RV of step S521, and construct the spatial variation function γ(h). Calculate the variation function values ​​in different directions to identify the anisotropic characteristics of the residual spatial distribution. Based on the variation function fitting results, determine the optimal search radius R and anisotropy parameter θ. Perform detrending analysis on the residual data to separate the global trend term GT and the local fluctuation term LT. Construct the residual spatial distribution model M(x,y) = GT + LT to generate the trend surface prediction residual RT.

[0087] Step S523, obtain the residual trend surface RT and local fluctuation term LT output by step S522. Construct an adaptive kernel function K(d), where d is the spatial distance. For each predicted position p, calculate the influence weight w of the surrounding sites based on K(d). Extract the local fluctuation feature LTp within the influence range. Apply the local polynomial interpolation method to calculate the local residual prediction value rp of point p. Combine RT and rp to generate the residual distribution map RD for the entire region.

[0088] Step S524, obtain the residual distribution map RD of step S523 and the downscaling result D1 of step S51. For each 1km grid point, extract its residual correction value rc. Construct an adaptive correction function f(rc), and the correction coefficient varies nonlinearly with the residual amplitude and reliability. Perform a local consistency check on the correction results to ensure the spatial continuity of the corrected data. Calculate the local statistical characteristics of each grid point before and after correction, and generate the correction coefficient matrix K and uncertainty estimate U. Combine the corrected data, correction coefficient and uncertainty estimate to generate the final 1km resolution soil moisture data product D2.

[0089] In this embodiment, the data processing flow is:

[0090] S1. Multi-source data preprocessing and quality assessment

[0091] S11. Data reading and basic preprocessing: read six sets of original soil moisture data: CCI (X1), CLM (X2), ERA_Interim (X3), MOS (X4), NOAH (X5), VIC (X6); perform spatiotemporal registration on each set of data, unify the projection coordinate system to WGS84; unify the data unit to m³ / m³; output the standardized six sets of data sets Y1-Y6;

[0092] S12, realize dynamic weight quality assessment system through three-layer quality assessment framework (TQL, Triple-layer Quality Level)

[0093] The first layer: data integrity assessment, calculate the missing rate Mi = missing_pixels / total_pixels of each data set; evaluate time continuity Ti = continuous_days / total_days; generate integrity index Ci = f(Mi, Ti), where f is a weighted function; the second layer: data consistency assessment; calculate the data change rate at adjacent moments Vi = |SMt - SMt-1| / Δt; construct the data fluctuation envelope Ei = g(Vi, t); generate consistency index Li = h(Vi,Ei); the third layer: physical rationality assessment; introduce the "soil moisture physical constraint model" Pi = k(precipitation, evaporation, terrain); o calculate the deviation from the physical model Di = |SMi - Pi|; generate rationality index Ri = m(Di); output the comprehensive quality index Qi = w1Ci + w2Li + w3Ri for each data set;

[0094] S13, filtering the time-varying data noise by adaptive multiscale sliding window filter (AMWF): construct multiscale time windows {t1, t2, ..., tn}; for each spatial position (x,y), calculate the local variance σ²(x,y,t); determine the optimal time window scale ts = argmin(σ²); apply adaptive weight filtering; finally output the denoised data sets Z1-Z6;

[0095] S14. Standardization and assimilation of multi-source data through the distribution adaptive normalization method (DANM): Calculate the

[0096] The probability distribution characteristics Pi(x); construct the reference distribution Pr(x); realize data conversion through probability mapping, Yi=F⁻¹r(Fi(Xi)), where F is the cumulative distribution function; output the standardized data set W1-W6;

[0097] S15, quality assessment index generation, including: calculating the time dimension reliability index Rt, the space dimension reliability index Rs, and the comprehensive reliability index R = φ(Rt, Rs) based on the quality index Qi of S12, and outputting the final quality assessment matrix Q and the reliability index set R.

[0098] Where: SMt represents the soil moisture value at time t; Δt is the time interval; w1, w2, w3 are weight coefficients, and w1 +w2 + w3 = 1; φ is the comprehensive evaluation function; F and F⁻¹ represent the forward and reverse probability distribution mapping functions respectively; the final output includes: six standardized data sets (W1-W6), quality assessment matrix (Q) and reliability index set (R).

[0099] S2, spatiotemporal correlation analysis and anomaly detection,

[0100] S21: Multidimensional spatiotemporal correlation modeling is performed through the hierarchical spatiotemporal feature extractor (HSFE), inputting the standardized data set {W1-W6} and the quality matrix Q; constructing the temporal correlation matrix, Tt(i,j) = corr(Wi(t), Wj(t+τ)), where τ is the time lag; constructing the spatial correlation matrix: Ts(x,y) = corr(Wi(x,y), Wi(x+Δx,y+Δy)); outputting the spatiotemporal correlation feature set ST;

[0101] S22. Outlier identification and classification through environment-aware anomaly detector (EADD):

[0102] Read ST and reliability index R, calculate multi-scale anomaly index: At = |Wi(t) - μt| / σt (time anomaly), As = |Wi(x,y) - μs| / σs (spatial anomaly), Ae = f(precipitation, temperature) (environmental anomaly); output the outlier labeling matrix A and the anomaly type matrix C;

[0103] S23, using adaptive tensor completion algorithm to repair spatiotemporal continuity: input {W1-W6}, A, C:

[0104] Construct a third-order tensor T(x,y,t), select a repair strategy based on the anomaly type C, and output the repaired data set {W1'-W6'}

[0105] S3, Adaptive region segmentation and boundary optimization

[0106] S31, terrain feature extraction, read DEM data, {W1'-W6'}, calculate terrain indicators: slope S, slope aspect A, terrain undulation R; construct terrain feature vector TF, and output terrain feature set G.

[0107] S32, vegetation impact assessment, read NDVI data, {W1'-W6'}, build "vegetation-soil moisture interaction model", VI = h(NDVI, LAI, soil_type), calculate vegetation impact coefficient β, and output vegetation feature set V.

[0108] S33. Generate dynamic partitions according to the adaptive boundary generation algorithm (ABGA): read G, V, {W1'-W6'}, calculate the optimal number of partitions: N = f(terrain complexity, vegetation density), generate initial partitions P, and output the partition result D.

[0109] S34. Optimize the boundary transition zone through the soft boundary transition model (SBTM): read D, {W1'-W6'}, and construct the transition function: B(x,y) = ω1R1 + ω2R2; Ri is the adjacent area feature; output the optimized partition D' and transition zone feature B.

[0110] S4. Multi-source data fusion and error propagation control

[0111] S41, Regional Adaptive Fusion

[0112] Read {W1'-W6'}, D', B, and build a "partition cascade fusion network": apply different weights to each partition, consider the transition zone feature B; output the initial fusion result F1.

[0113] S42, suppress error propagation through error chain tracking algorithm (ECTA): read F1, quality matrix Q, build error propagation model, apply error compensation, and output optimized fusion result F2;

[0114] S5. Terrain-aware downscaling and accuracy assessment

[0115] S51, multi-factor downscaling model construction, read F2, terrain characteristics G, vegetation characteristics V, build a "terrain adaptive downscaling framework": terrain correction factor K1, vegetation influence factor K2 and scale conversion function S(x), and output the preliminary downscaling result D1.

[0116] S52, optimize and verify accuracy through multi-scale residual learning model (MRLM): read D1, measured data, calculate residual distribution, apply residual correction, and output the final 1km resolution data product D2.

[0117] The character descriptions are as follows: Wi: standardized data of the ith data source; μt, σt: mean and standard deviation of the time dimension;

[0118] μs, σs: mean and standard deviation of spatial dimension; ω1, ω2: weight coefficients; τ: time lag parameter; Δx, Δy: spatial offset; the final output is soil moisture data product D2 with a resolution of 1 km, including: soil moisture value, quality assessment indicators, and uncertainty estimation.

[0119] According to another aspect of the present application, a surface soil moisture multi-source data fusion and downscaling system is also provided, which is characterized by comprising:

[0120] at least one processor; and,

[0121] a memory communicatively connected to at least one of the processors; wherein,

[0122] The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the surface soil moisture multi-source data fusion and downscaling method described in any of the above technical solutions.

[0123] In another embodiment of the present application, steps S1 and S2 are as in the previous embodiment. Steps S3 to S5 are further as follows:

[0124] Step S3, obtain the restored data set W1'-W6' output by step S2 and the pre-stored DEM data, calculate the elevation gradient in the x-direction and y-direction in each grid, and obtain the slope data S. Use the elevation gradient to calculate the slope data A of each grid. Calculate the elevation difference between the central grid and the surrounding grids in each 3×3 window to obtain the terrain roughness index R. Calculate the maximum elevation difference in the 5×5 window to obtain the terrain undulation H. Combine S, A, R, and H to generate a terrain feature vector TF, and combine the terrain type data to generate a terrain feature set G. Obtain NDVI time series data and soil type data, calculate the maximum and minimum NDVI values ​​of each grid, and obtain the vegetation density index VCI. Calculate the vegetation impact coefficient β based on VCI, soil texture data and slope S. Combine β with the seasonal factor to generate the vegetation impact matrix VM, and combine VCI, β and VM to obtain the vegetation feature set V. Calculate the terrain complexity and vegetation density based on G and V to generate the optimal number of partitions N. Combine G, V and W1'-W6' in the feature space, apply the watershed algorithm, and get the initial partition D. Extract the boundary pixel B of D, and for each boundary pixel, obtain the adjacent area features R1 and R2, combine the local terrain feature G, construct the transition function T, and generate the optimized partition D' and transition zone feature B.

[0125] Step S31, obtain the pre-stored DEM data, for each grid position, extract the elevation difference of the adjacent grids in the x direction and divide it by the horizontal distance to obtain the x-direction gradient data Gx. Extract the elevation difference of the adjacent grids in the y direction and divide it by the horizontal distance to obtain the y-direction gradient data Gy. Take the square root of the square sum of Gx and Gy to obtain the slope data S. Use the inverse tangent value of Gx and Gy to calculate the slope data A. Construct a 3×3 window with each grid as the center, calculate the average of the sum of the squares of the elevation differences between all grids in the window and the central grid, and then take the square root to obtain the terrain roughness index R. Construct a 5×5 window, calculate the difference between the maximum elevation value and the minimum elevation value in the window, and obtain the terrain undulation H. Combine S, A, R, and H at each position into a vector to generate a terrain feature vector TF. Obtain the pre-stored terrain type data, combine it with TF, and generate a terrain feature set G.

[0126] Step S32: Obtain the pre-stored NDVI time series data, extract the annual NDVI sequence for each grid position, and calculate its maximum NDVImax and minimum NDVImin. Subtract NDVImin from the current NDVI value and divide it by the difference between the maximum and minimum values ​​to obtain the vegetation density index VCI. Obtain the pre-stored soil type data and extract the soil texture information at each location.

[0127] Obtain the slope data S of step S31, multiply VCI, soil texture and slope S by preset weight coefficients k1, k2, k3 respectively, and sum them to obtain vegetation influence coefficient β. Obtain seasonal factor data, perform function operation with β, and generate vegetation influence matrix VM. Combine VCI, β and VM to generate vegetation feature set V.

[0128] Step S33, obtain the terrain feature set G of step S31 and the vegetation feature set V of step S32, calculate the terrain complexity of the study area, and multiply the complexity value by the coefficient α. Calculate the vegetation density value and multiply it by the coefficient β. Add the weighted terrain complexity and vegetation density to obtain the optimal number of partitions N. Obtain the data set W1'-W6' output by step S2, and combine G, V and soil moisture data to construct the feature space F. Apply the watershed algorithm in the F space, introduce the terrain gradient as a segmentation constraint, and consider the spatial continuity of the vegetation to generate the initial partition D with the area number and boundary point coordinates.

[0129] Step S34, obtain the initial partition D of step S33, and extract the set B of all boundary pixels. For each pixel in B, obtain the features R1 and R2 of the adjacent areas on both sides. Obtain the value G(x, y) of the terrain feature G of step S31 at the pixel position. Multiply R1, R2 and G(x, y) by the weight coefficients w1, w2, w3 respectively, and sum them to construct the transition function value T(x, y) of the pixel. Apply tensor voting method to smooth all boundary pixels while maintaining the continuity of terrain features. Generate the optimized partition boundary, obtain the optimized partition D', and record the transition zone feature B.

[0130] Step S4, obtain the restored data set W1'-W6' output by step S2, the optimized partition D' and transition zone feature B, terrain feature set G, and vegetation feature set V output by step S3, and construct a regional feature matrix RF for each partition, including soil moisture values, terrain features, and vegetation features. Obtain the quality matrix Q output by step S1, and calculate the fusion weight w of each data source in each partition in combination with the temporal stability of the data. Apply the corresponding weights to the soil moisture data in each partition for weighted fusion to generate the initial fusion result F1. Obtain the error estimate E of each processing link, construct the error transfer equation, and calculate the error E(t+1) at the next moment. Calculate the error compensation value C based on E and spatial correlation, subtract the compensation value C from F1, and obtain the optimized fusion result F2.

[0131] Step S41, obtain the data set W1'-W6' output by step S2, the partition D' output by step S3, the transition zone feature B, the terrain feature set G and the vegetation feature set V. For each partition, the soil moisture data, terrain features and vegetation features are combined to construct the regional feature matrix RF. Obtain the quality matrix Q output by step S1, and calculate the stability index of each data source at different time points. Input the quality index and time stability into the weight calculation function to obtain the fusion weight w of each data source in each partition. For the soil moisture data in each partition, apply the corresponding weight w for weighted averaging to generate the fusion value of the partition. In the transition zone area at the boundary of the partition, adjust the fusion weight according to the transition zone feature B. Combine the fusion results of all partitions to generate the initial fused data set F1.

[0132] Step S42: Obtain the initial fusion result F1 of step S41 and the quality matrix Q output by step S1. Obtain the pre-stored error estimation parameters α and β. For each time point t, multiply the current error E(t) by the coefficient α, and add the soil moisture change ΔSM(t) multiplied by the coefficient β to obtain the error estimate E(t+1) at the next moment. Calculate the spatial correlation in the study area, combine it with the error estimate E, and generate the error compensation coefficient C. Subtract the corresponding compensation value C from each data point in F1 to generate the optimized fusion result F2.

[0133] Step S5, obtain the fusion result F2 output by step S4, the terrain feature set G output by step S3, the vegetation feature set V and the pre-stored surface temperature LST data, and construct the scale conversion function S(x). Calculate the terrain correction coefficient K1 using the slope, aspect and elevation data. Calculate the vegetation influence coefficient K2 based on the NDVI and LAI data. Multiply F2 with S(x), K1 and K2 to obtain the preliminary downscaling result D1. Obtain the ground measured data M, calculate the difference with D1 to obtain the residual R. Construct a residual model RM that includes terrain, vegetation and soil types, correct D1, and obtain the final 1km resolution data product D2.

[0134] Step S51, obtain the fusion result F2 output by step S4, the terrain feature set G output by step S3, the vegetation feature set V and the pre-stored surface temperature LST data. Input the terrain index and vegetation index into function h to construct the scale conversion function S(x). Obtain the slope, aspect and elevation data, input function f to calculate the terrain correction coefficient K1. Obtain NDVI and LAI data, input function g to calculate the vegetation influence coefficient K2. Multiply each pixel value of F2 by S(x), K1, and K2 at the corresponding position to generate preliminary 1km resolution data D1.

[0135] Step S52, obtain the downscaling result D1 of step S51 and the pre-stored ground measured data M. Calculate the difference between the corresponding position values ​​of M and D1 to obtain the residual sequence R. Calculate the standard deviation σ_r of R. Obtain terrain characteristics, vegetation characteristics and soil type data, input the residual model function to obtain the residual prediction value RM. Add D1 to RM to obtain the corrected data D2. Calculate the quality index and uncertainty estimate of each pixel, and combine them with D2 to form the final 1km resolution soil moisture data product.

[0136] The specific data flow of Example 2 is as follows:

[0137] S3, Adaptive region segmentation and boundary optimization

[0138] S31. Terrain feature extraction

[0139] Read DEM data D and the restored soil moisture dataset {W1'-W6'}; calculate terrain gradient: slope S = ∂z / ∂x + ∂z / ∂y, where z is elevation; slope aspect A = arctan(∂z / ∂y, ∂z / ∂x);

[0140] Calculate the terrain roughness index: R = √(1 / n∑(zi-z̄)²), zi is the local elevation value;

[0141] Construct terrain feature vector: TF = [S, A, R, H], H is terrain relief; output terrain feature set G = {TF, R, terrain_type};

[0142] S32, vegetation impact assessment, read NDVI time series data N, soil type data ST, G;

[0143] Calculate the vegetation density index: VCI = (NDVI - NDVImin) / (NDVImax - NDVImin);

[0144] Construct vegetation-soil moisture relationship: β = k1·VCI + k2·soil_texture + k3·slope; k1, k2, k3 are weight coefficients

[0145] Generate vegetation impact matrix: VM = f(β, seasonal_factor); output vegetation feature set V = {VCI,β, VM}

[0146] S33: Generate dynamic partitions using terrain-vegetation collaborative partitioning algorithm (TVCA): read G, V, {W1'-W6'};

[0147] Calculate the optimal number of partitions: N = α·terrain_complexity + β·vegetation_density; α, β are adaptive weights

[0148] Construct feature space: F = [G, V, soil_moisture]; apply watershed segmentation: introduce terrain gradient constraints, consider vegetation continuity, and output initial partition D = {region_id, boundary_points}.

[0149] S34. Optimize the boundary transition zone through the multi-dimensional boundary smoothing algorithm (MBSA); read D, {W1'-W6'}, G, V; identify the boundary pixel B = {p | p∈boundary(D)}; construct the transition function: T(x,y) = w1·R1 + w2·R2 +w3·G(x,y); Ri is the adjacent area feature.

[0150] Optimize the boundary: Apply the tensor voting method to smooth the boundary, maintain the continuity of terrain features, and output the optimized partition D' and transition zone features B

[0151] S4, Multi-source data fusion and error propagation control,

[0152] S41. Regional adaptive fusion through multi-scale feature fusion network (MSFN): read {W1'-W6'}, D', B, G, V; construct regional feature matrix: RF = [soil_moisture, terrain, vegetation]; calculate fusion weight:

[0153] w(i,j) = f(quality_index, temporal_stability); perform partition fusion SM_fused =∑(wi·SMi); output the initial fusion result F1

[0154] S42. Suppress error propagation through hierarchical error compensation model (HECM): read F1, quality matrix Q, error estimate E; construct error propagation model: E(t+1) = α·E(t) + β·ΔSM(t); calculate error compensation: C = g(E,spatial_correlation); apply error correction: F2 = F1 – C; output optimized fusion result F2.

[0155] S5. Terrain-aware downscaling and accuracy assessment

[0156] S51. Construct a multi-factor downscaling model through the adaptive multiscale decomposition algorithm (AMDA): read F2, G, V, and land surface temperature LST; construct a scale conversion function: S(x) = h(terrain_factor, vegetation_index); calculate terrain correction: K1 = f(slope, aspect, elevation); apply vegetation influence: K2 = g(NDVI, LAI); integrate downscaling:

[0157] D1 = F2·S(x)·K1·K2; output the preliminary downscaling result D1.

[0158] S52: Optimize and verify accuracy through the residual learning optimization framework (RLOF): read D1, measured data M, quality index Q; calculate residual distribution: R = M - D1; σ_r = std(R); build residual model: RM = f(terrain,vegetation, soil_type); apply correction: D2 = D1 + RM; evaluate accuracy: calculate RMSE, MAE, R²; output final product D2={SM_1km, quality_index, uncertainty}. Where: ∂z / ∂x, ∂z / ∂y are elevation gradients

[0159] z̄ is the mean elevation; σ_r is the residual standard deviation; SM_1km is the soil moisture at 1 km resolution; quality_index is the quality index; uncertainty is the uncertainty estimate.

[0160] It should be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction. In order to avoid unnecessary repetition, the present invention will not further describe various possible combinations.

Claims

1. The multi-source data fusion and downscaling method of surface soil moisture is characterized by: The steps include: Step S1, obtain six sets of original soil moisture data and unify the coordinate system and unit to obtain data sets Y1-Y6; calculate the data missing rate Mi and time continuity index Ti of data sets Y1-Y6 to generate the integrity index Ci; calculate the consistency index Li based on the data change rate Vi; use the physical constraint model to calculate the deviation Di between the theoretical value Pi and the actual value to generate the rationality index Ri; combine the integrity index Ci, consistency index Li, and rationality index Ri to obtain the quality index Qi; perform noise reduction and standardization on data sets Y1-Y6 to obtain data sets W1-W6 and their quality assessment matrix Q and reliability index R; Step S2: Based on the data set W1-W6, the quality assessment matrix Q and the reliability index R, the time correlation matrix Tt and the space correlation matrix Ts are calculated, and the space-time feature set ST is generated by combining them; based on the space-time feature set ST, the time anomaly index At, the space anomaly index As and the environment anomaly index Ae are calculated to generate the anomaly labeling matrix A and the anomaly type matrix C; the data set W1-W6 is expanded into a tensor T(x, y, t), and the outliers are repaired according to the anomaly type matrix C to obtain the repaired data set W1'-W6'; Step S3, obtain the restoration data set W1'-W6' and the pre-stored DEM data, calculate the slope S, slope aspect A and terrain undulation R of each spatial position, generate the terrain feature vector TF, and combine them to obtain the terrain feature set G; obtain NDVI and LAI data, combine with soil type data, calculate the vegetation cover index VI, generate the vegetation influence coefficient β, and form the vegetation feature set V; according to the terrain feature set G and the vegetation feature set V, calculate the terrain complexity and vegetation density, input the optimization function to obtain the optimal partition number N, and use the regional growing algorithm to generate the initial partition P; For each partition boundary in the initial partition P, extract the adjacent region features R1 and R2, calculate the transition zone feature B through the weight function ω, and generate the optimized partition D'; Step S4, obtaining the repaired data set W1'-W6', the optimized partition D' and the transition zone feature B, calculating the weight coefficient αi for the data in each partition, fusing the weighted data, and generating an initial fusion result F1; Obtain the quality matrix Q, calculate the error transfer coefficient εj of each processing link, and build an error accumulation model; input the initial fusion result F1 into the error accumulation model, and calculate the error estimation value e of each position; perform error compensation processing on the initial fusion result F1 according to the error estimation value e to obtain the optimized fusion result F2; Step S5, obtaining the fused data set F2, the terrain features G and the vegetation features V, and calculating the terrain correction coefficient K1 and the vegetation influence coefficient K2; The fused data set F2, terrain correction coefficient K1 and vegetation influence coefficient K2 are input into the scale conversion function to generate preliminary 1km resolution data D1; Obtain the pre-stored measured data of ground stations and calculate the residuals with the corresponding position data of the preliminary 1km resolution data D1; perform spatial distribution analysis on the residuals, establish a residual correction model, correct the preliminary 1km resolution data D1, and obtain the final 1km resolution data product D2.

2. The method for fusion and downscaling of multi-source data of surface soil moisture according to claim 1, characterized in that: The step S1 is specifically as follows: Step S11, obtain six sets of original soil moisture data CCI, CLM, ERA_Interim, MOS, NOAH and VIC, and record them as data sets X1-X6 in sequence, read their latitude and longitude information and timestamp information; convert the projection coordinate system of X1-X6 to WGS84 to generate a unified coordinate grid; unify the units of the soil moisture values ​​in the data sets X1-X6 to generate standardized data sets Y1-Y6; Step S12, according to the standardized data sets Y1-Y6, count the total number of pixels total_pixels and the number of missing pixels missing_pixels of each data set, and calculate the missing rate Mi; obtain the time series of the standardized data sets Y1-Y6, calculate the continuous observation days continuous_days and the total observation days total_days, and obtain the time continuity index Ti; input the missing rate Mi and the time continuity index Ti into the weighting function f(x) to generate the integrity index Ci; based on the time series of the standardized data sets Y1-Y6, calculate the data difference between adjacent moments divided by the time interval to obtain the change rate Vi; Construct the data fluctuation range according to the change rate Vi and generate the envelope Ei; input the envelope Vi and the envelope Ei into the function h(x) to obtain the consistency index Li; obtain precipitation, evaporation and terrain data, input the physical constraint model k(x) to obtain the theoretical soil moisture value Pi; calculate the absolute difference between the standardized data set Y1-Y6 and the theoretical soil moisture value Pi to obtain the deviation Di; input the deviation Di into the function m(x) to generate the rationality index Ri; Step S13, construct sliding windows of different scales {t1, t2, ..., tn} for the time series of the standardized data sets Y1-Y6; calculate the data variance σ²(x, y, t) in each time window at each spatial position (x, y); select the time window ts with the smallest variance, and perform weighted average filtering on the data at this position; Repeat this process for all spatial locations to obtain the denoised data sets Z1-Z6; Step S14: Calculate the probability distribution features Pi(x) of each of the denoised data sets Z1-Z6; construct a reference probability distribution Pr(x); map the denoised data sets Z1-Z6 to the reference distribution space through the cumulative distribution function Fi(x); apply the reverse mapping function F⁻¹r to generate the standardized data sets W1-W6; Step S15, based on the quality index Qi, calculate the time dimension reliability index Rt for the time series of the standardized data set W1-W6; calculate the spatial dimension reliability index Rs for the spatial distribution of the standardized data set W1-W6; input Rt and Rs into the comprehensive evaluation function φ(x) to generate the final quality evaluation matrix Q and reliability index set R.

3. The method for fusion and downscaling of multi-source data of surface soil moisture according to claim 2, characterized in that: The step S2 is specifically as follows: Step S21, according to the standardized data sets W1-W6 and the quality matrix Q outputted in step S1, set the time delay sequence τ={1,2,...,n}; for any two data sets Wi and Wj, calculate the Pearson correlation coefficient at each τ to obtain the time correlation matrix Tt(i,j); For each data set Wi, set the spatial displacement Δx and Δy, calculate the correlation coefficient between the position (x, y) and (x+Δx, y+Δy), and generate the spatial correlation matrix Ts(x, y); combine Tt and Ts to generate the spatiotemporal correlation feature set ST; Step S22, using the spatiotemporal correlation feature set ST and the reliability index R, calculate the mean μt and standard deviation σt in the local time window at each time point t; calculate the difference between the data value Wi(t) and μt, and divide it by σt to obtain the time anomaly index At; calculate the mean μs and standard deviation σs of the local area for each spatial position (x, y), calculate the difference between the data value Wi(x, y) and μs, and divide it by σs to obtain the spatial anomaly index As; obtain the precipitation and temperature data of the corresponding position and time, input the environmental anomaly function f(x) to calculate the environmental anomaly index Ae; based on the combination of the temporal anomaly index At, the spatial anomaly index As and the environmental anomaly index Ae, generate the anomaly labeling matrix A, and generate the anomaly type matrix C according to the dominant anomaly type; Step S23, obtain the standardized data set W1-W6, the outlier labeling matrix A and the outlier type matrix C, and construct a three-dimensional tensor T(x, y, t); for the data points marked as abnormal, select the time interpolation, space interpolation or time-space joint interpolation strategy according to the anomaly type in the anomaly type matrix C; for time anomalies, use the weighted average of adjacent time points for repair; for space anomalies, use local space interpolation for repair; for environmental anomalies, combine space-time information and environmental factors for comprehensive repair; after completing the repair of all abnormal points, obtain the repaired data set W1'-W6'.

4. The method for fusion and downscaling of multi-source data of surface soil moisture as described in claim 3, characterized in that: The step S3 is specifically as follows: Step S31, obtaining pre-stored DEM data, and for each grid position (x, y), calculating the elevation difference between the central pixel and the adjacent pixel based on a 3×3 neighborhood to obtain slope data S; Calculate the elevation gradient of each pixel in the four directions of east, west, south and north to obtain the aspect data A; calculate the ratio of the maximum elevation difference to the horizontal distance in a 5×5 neighborhood to obtain the terrain relief R; combine the slope data S, the aspect data A and the terrain relief R into a terrain feature vector TF; perform cluster analysis on the TF of all spatial positions to generate a terrain feature set G; Step S32, obtaining pre-stored NDVI data and LAI data, and combining with pre-stored soil type data to calculate the vegetation coverage of each location; performing time series filtering on the NDVI value to remove abnormal fluctuations; calculating the thickness of the vegetation layer according to the LAI data; inputting the NDVI, LAI and soil type data into the vegetation impact model to calculate the vegetation water interception coefficient; calculating the vegetation impact coefficient β based on the interception coefficient and coverage; combining all vegetation related parameters to generate a vegetation feature set V; Step S33, obtaining a terrain feature set G, a vegetation feature set V, and a restored data set W1'-W6'; calculating the standard deviation of the terrain features in the study area to obtain a terrain complexity index; statistically analyzing the spatial distribution density of vegetation features to obtain a vegetation density index; inputting the terrain complexity and vegetation density into the optimization function to calculate the optimal partition number N; selecting N initial seed points, performing regional growth based on feature similarity, and generating an initial partition result P, where N is a positive integer; Step S34, obtain the initial partition P of step S33, and extract each partition boundary line; for each boundary line, obtain the terrain characteristics R1 and soil moisture characteristics R2 of the adjacent areas on both sides; calculate the weight coefficients ω1 and ω2, input the terrain characteristics R1 and the soil moisture characteristics R2 into the weighting function, and obtain the boundary transition characteristics; expand a specific width to both sides along the boundary to construct a transition zone, and record the characteristic value B of each point in the transition zone; combine the original partition P with the transition zone feature B to generate an optimized partition D' and a complete transition zone feature set B.

5. The method for fusion and downscaling of multi-source data of surface soil moisture according to claim 4, characterized in that: The step S4 is specifically as follows: Step S41, obtaining the repaired data set W1'-W6', partition D' and transition zone feature B; for each partition, extracting the time series features and spatial distribution features of each data source in the partition; Calculate the similarity between each data source and the overall characteristics of the region to obtain the initial weight; Adjust the initial weight based on the mass matrix Q to generate the final weight coefficient αi; For the data in each partition, the weight coefficient αi is applied for weighted averaging; In the transition zone area, the dynamic weight is calculated according to the transition zone feature B, and the fusion results of adjacent partitions are smoothly transitioned; the fusion results of all partitions are combined to generate the initial fusion result dataset F1; Step S42, obtaining the initial fusion result data set F1 and the quality matrix Q; calculating the error transfer coefficient εj for each link in the data processing flow, including standardization, anomaly detection, partition optimization and fusion; Based on the error transfer coefficient, an error accumulation function is constructed to analyze the error transfer path in the processing flow; for each data point in the initial fusion result data set F1, its cumulative error estimate e is calculated; according to the size and distribution characteristics of the error estimate e, an error compensation strategy is designed; the compensation strategy is applied to the initial fusion result data set F1 to obtain the optimized fusion result F2; In the compensation process, local re-fusion is adopted for high error areas, and the original fusion results are kept for low error areas; The compensated results are integrated to generate the final fused data set F2.

6. The method for multi-source data fusion and downscaling of surface soil moisture according to claim 5, characterized in that: The step S5 is specifically as follows: Step S51, obtain the final fused data set F2, terrain features G and vegetation features V; calculate the terrain relief, slope change rate and slope aspect variability in each 1km grid, and generate the terrain correction coefficient K1; count the vegetation coverage change and vegetation type distribution in each 1km grid, and calculate the vegetation influence coefficient K2; construct a spatial interpolation grid from 25km to 1km; for each 25km pixel, extract the terrain correction coefficient K1 and vegetation influence coefficient K2 within its corresponding range, and establish a local terrain-vegetation weight model; distribute the values ​​of the final fused data set F2 to the 1km grid through the local weight model to generate preliminary high-resolution data D1; Step S52, obtaining the preliminary downscaling result D1 and the pre-stored measured data of the ground station; for each measured station location, extracting the data of the corresponding location in the preliminary downscaling result D1 and the surrounding 3×3 grid; calculating the difference between the measured value and the preliminary downscaling result D1 to obtain a residual sequence; performing spatial autocorrelation analysis on the residual sequence to calculate the spatial distribution characteristics of the residual; Based on the spatial distribution characteristics of residuals, a Kriging interpolation model is constructed to generate a residual distribution map for the entire region. The preliminary downscaling result D1 is corrected according to the residual distribution map, and the corresponding residual correction value is applied to each 1km grid. The corrected data are combined to generate the final 1km resolution soil moisture data product D2. At the same time, the correction coefficient and uncertainty estimate of each grid point are generated as quality assessment indicators.

7. The method for multi-source data fusion and downscaling of surface soil moisture according to claim 6, characterized in that: The step S15 is specifically as follows: Step S151, obtain the quality index Qi, extract the time series data for each data set W1-W6; calculate the autocorrelation coefficient rt(k) of each time series, where k is the time lag; construct the time continuity feature vector Vt based on the autocorrelation coefficient rt(k) and apply principal component analysis to extract the main time change mode Mt; compare the main time change mode Mt with the preset time change threshold, and calculate the time dimension deviation δt; calculate the time dimension reliability index Rt based on the time dimension deviation δt and the quality index Qi; Step S152, obtaining standardized data sets W1-W6, and constructing a local area matrix L(x, y) for each spatial position; Calculate the spatial autocorrelation coefficient rs(h) in the local area matrix L(x,y), where h is the spatial distance; construct the spatial continuity feature vector Vs based on the spatial autocorrelation coefficient rs(h); perform variogram analysis on the spatial continuity feature vector Vs to obtain the spatial structure feature Gs; compare the spatial structure feature Gs with the regional geographical features to calculate the spatial dimension deviation δs; calculate the spatial dimension reliability index Rs based on the spatial dimension deviation δs and the quality index Qi; Step S153, obtain the time dimension reliability index Rt and the space dimension reliability index Rs; construct the spatiotemporal weight matrix W(t,s); for each data point, extract the time dimension reliability index Rt and the space dimension reliability index Rs values ​​in its time and space neighborhood; apply the spatiotemporal weight matrix W(t,s) to perform weighted combination of the time dimension reliability index Rt and the space dimension reliability index Rs to generate an initial comprehensive index I; apply local smoothing filtering to the initial comprehensive index I to eliminate mutation points; standardize the filtered results to generate the final quality assessment matrix Q; calculate the reliability level of each position based on the quality assessment matrix Q to form a reliability index set R.

8. The method for fusion and downscaling of multi-source surface soil moisture data according to claim 6, characterized in that: The step S23 is specifically as follows: Step S231, obtain data sets W1-W6, construct a time-space-observation value tensor T(x, y, t); calculate the sliding window statistics St of the time-space-observation value tensor T in the time dimension, including the mean, standard deviation and trend coefficient; calculate the local statistics Ss of the time-space-observation value tensor T in the spatial dimension, including spatial gradient and curvature characteristics; The sliding window statistics St and the local statistics Ss are combined to form the spatiotemporal feature tensor TSF; based on the spatiotemporal feature tensor TSF, the local correlation matrix LC is constructed for subsequent anomaly repair; Step S232, obtain the anomaly type matrix C, extract the time anomaly point set At; for each anomaly point in the time anomaly point set At, obtain its TSF feature vector and LC correlation data; construct an adaptive time weight function wt, and the weight coefficient decays with the time distance and correlation; based on the adaptive time weight function wt, weight the effective values ​​of the n moments before and after the anomaly point to generate a time dimension repair value Rt; apply physical constraints to the repair result to make the repair value within a reasonable range; Step S233, obtain a set of spatial outlier points As; for each outlier point, construct a spatial influence domain D(x, y); extract terrain gradient g and vegetation coverage v in the spatial influence domain D(x, y); construct a spatial weight function ws based on the terrain gradient g and vegetation coverage v, and the weight varies with the terrain similarity and vegetation consistency; apply the spatial weight function ws to weight the valid values ​​in the spatial influence domain D(x, y) to generate a spatial dimension repair value Rs; perform a consistency check on the spatial dimension repair value Rs and the surrounding grid points to ensure spatial continuity; Step S234, obtaining a temporal repair value Rt and a spatial repair value Rs; Construct a repair optimization function f(Rt, Rs) to calculate the credibility index of time and space repair based on the TSF feature tensor; select the optimal repair value according to the credibility index, or perform a weighted combination according to the credibility index; Perform localized iterative optimization on the repair results until the first-order and second-order derivatives of adjacent spatiotemporal positions meet the continuity conditions; The repaired data are reorganized into the original data format to generate the repaired data sets W1'-W6'.

9. The method for multi-source data fusion and downscaling of surface soil moisture according to claim 6, characterized in that: The step S33 is specifically as follows: Step S331, obtain the terrain feature set G, calculate the coefficients of variation cv1 and cv2 of the slope and slope aspect in the study area; obtain the vegetation feature set V, calculate the spatial heterogeneity index h1 of the vegetation index; construct the complexity evaluation function f(cv1, cv2, h1), generate the regional complexity matrix CM; apply wavelet transform to the regional complexity matrix CM, extract the characteristic frequencies at different scales, and obtain the scale characteristic vector SF; Calculate the optimal partition number estimate N1 based on the scale feature vector SF; Step S332, obtain the complexity matrix CM, construct the spatial distance matrix D; calculate the local density ρ of each location point and the minimum distance δ to the high-density point; identify the density peak point on the ρ-δ graph, and select the first N1 points as the initial seed point set S; extract the feature vector FV of each seed point, the feature vector FV contains the terrain, vegetation and soil moisture characteristics; calculate the Mahalanobis distance between the seed points, and generate the seed point optimization set S'; N1 is a preset natural number; Step S333, obtain the optimized seed point set S', construct the feature space distance function d(x); for each non-seed point p, calculate its feature distance d(p,s) with all seed points; set the dynamic growth threshold θ, the initial value of which is determined based on the local feature variance; When the feature distance d(p,s) is less than θ, the non-seed point p is assigned to the nearest seed point; θ is adaptively adjusted during the region growth process, and the adjustment amplitude is inversely proportional to the feature variance in the current partition; the initial partition result P1 is generated; Step S334, obtaining the initial partition P1, constructing a partition boundary set B; for each boundary line in the partition boundary set B, extracting the feature sequence of the strip areas on both sides thereof; calculating the Jensen-Shannon divergence JS of the feature sequence as a boundary significance index; stratifying the boundaries based on the Jensen-Shannon divergence JS value and the set JS threshold, keeping the boundaries with high JS values ​​sharp, and smoothing the boundaries with low JS values; Perform a regional merging test within the partition, and merge adjacent partitions when the feature similarity exceeds a threshold; Generate the final partition result P.

10. A surface soil moisture multi-source data fusion and downscaling system, characterized in that: include: at least one processor; as well as, a memory communicatively connected to at least one of the processors; wherein, The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the surface soil moisture multi-source data fusion and downscaling method described in any one of claims 1 to 9.

Images