A method and system for lunar-scale lake water body reconstruction based on multi-physics constraints and adaptive regularized spatiotemporal fusion

Abstract

This invention discloses a method and system for monthly-scale lake water body reconstruction based on multi-physics field constraints and spatiotemporal fusion, belonging to the field of water body reconstruction technology. The method includes the following steps: S1, collecting multi-source heterogeneous data and preprocessing the data; S2, generating an orthogonal time basis matrix based on the preprocessed multi-source heterogeneous data; S3, performing adaptive weighted ridge regression modeling and solving based on the orthogonal time basis matrix to generate an MNDWI long-term image set; S4, diagnosing and repairing abnormal months based on the MNDWI long-term image set, and extracting monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset. This invention not only eliminates the jagged effect of low-resolution data but also solves the spatiotemporal gap problem of high-resolution data, ensuring the continuity and logical consistency of monitoring data in the physical process.

Description

Technical Field

[0001] This invention belongs to the field of water body reconstruction technology, specifically relating to a method and system for monthly-scale lake water body reconstruction based on multi-physics field constraints and adaptive regularized spatiotemporal fusion. Background Technology

[0002] As a core component of the "Asian Water Tower," the lake system of the Qinghai-Tibet Plateau not only maintains the balance of the regional water cycle but also serves as a sensitive indicator of global climate change. The dynamic changes in the extent of lake water bodies are directly related to regional ecological security, water resource management, and disaster early warning. Therefore, acquiring long-term, high-frequency, and high-spatial-resolution dynamic monitoring data of lakes is of significant scientific importance and practical value for revealing the expansion and contraction patterns of lakes on an interannual scale and seasonal water level fluctuations on an intraannual scale (such as snowmelt replenishment, monsoon precipitation, and freeze-thaw processes).

[0003] Currently, most global and regional surface water datasets are constructed based on low to medium spatial resolution or single sensors. To address the contradiction between high temporal frequency and high spatial detail that a single sensor cannot simultaneously capture, spatiotemporal image fusion technology has emerged and received widespread attention. Early research proposed the Spatiotemporal Adaptive Reflectance Fusion Model (STARFM), which uses the spectral correlation between high and low resolution images for prediction, but its prediction accuracy is limited when land cover types undergo complex changes. Subsequently, the Enhanced Spatiotemporal Adaptive Reflectance Fusion Model (ESTARFM) was proposed, which improved the prediction ability for heterogeneous landscapes by introducing conversion coefficients. Further, the Flexible Spatiotemporal Fusion Algorithm based on unmixing (FSDAF) was proposed, which improved the ability to capture abrupt changes in land cover to some extent. However, these methods have several shortcomings: First, the lack of physical mechanism constraints leads to misidentification of land cover. Existing purely mathematical statistical models mainly rely on the spectral correlation between images, lacking an understanding of surface physical processes. During the complex "water-ice-snow" phase transition period of the Qinghai-Tibet Plateau, cloud shadows are easily misidentified as water expansion, or floating ice is misidentified as land. Second, the lack of time control points makes it impossible to reconstruct nonlinear processes. In extreme cases where there are no effective Landsat observations for several consecutive months during the rainy season or freezing period on plateaus, traditional methods cannot reconstruct nonlinear changes in lakes, such as rapid snowmelt expansion or rapid freezing in extreme cold. Third, the limitation of a single optical sensor results in weak anti-interference capabilities. Existing algorithms struggle to overcome the bottleneck of optical observation and cannot meet the monitoring needs of all-weather, long-term data collection. Fourth, the limitation of the linear mixture assumption leads to poor adaptability to complex and abrupt changes in land cover. When the land cover type undergoes complex changes, its prediction accuracy is severely limited. Therefore, there is an urgent need for a spatiotemporal reconstruction method for lake water bodies that can overcome the limitations of a single optical sensor, deeply fuse multi-source heterogeneous data, and introduce physical field strength constraints such as meteorological and microwave data to construct all-weather, long-term, and high-precision lake water body dataset products. Summary of the Invention

[0004] To address the problem of spatiotemporal gaps caused by extreme weather interference and scarce observation data in the monitoring of lakes on the Qinghai-Tibet Plateau, this invention proposes a method and system for monthly-scale lake water body reconstruction based on multi-physics field constraints and adaptive regularized spatiotemporal fusion.

[0005] The technical solution of this invention is: a method for reconstructing lunar-scale lake water bodies based on multi-physics constraints and adaptive regularized spatiotemporal fusion, comprising the following steps:

[0006] S1. Collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data;

[0007] S2. Based on the preprocessed multi-source heterogeneous data, generate an orthogonal time basis matrix;

[0008] S3. Based on the orthogonal time basis matrix, perform adaptive weighted ridge regression modeling and solution to generate MNDWI long time series image set;

[0009] S4. Based on the MNDWI long-term image set, perform abnormal month diagnosis and repair, and extract the monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset.

[0010] Furthermore, the multi-source heterogeneous data includes MODIS high-frequency imagery, Landsat high-resolution imagery, as well as ERA5 temperature data, Sentinel-1 SAR data, and SRTM topography.

[0011] Furthermore, in S1, missing values ​​are filled in for the MODIS high-frequency image, and the image is resampled according to the MNDWI index, vegetation index, and snow cover index to generate a time series matrix.

[0012] The MNDWI index was used as the main variable for fusion, and the vegetation index, snow cover index and quality label band were used as auxiliary variables to determine the number of historical valid observations of a pixel.

[0013] Convert the Kelvin temperature of the ERA5 air temperature to Celsius to generate a time series of lake air temperature;

[0014] Sentinel-1 SAR data is filtered, and the filtered image is binarized into water bodies and non-water bodies using the backscattering coefficient threshold. The proportion of water body pixels is counted to generate the SAR water body area ratio index.

[0015] Slope and distance from shore are calculated based on SRTM terrain.

[0016] Furthermore, S2 includes the following sub-steps:

[0017] S21. Using the lake temperature time series and SAR water area ratio index, the time series matrix is ​​physically cleaned to obtain the cleaned matrix.

[0018] S22. Perform singular value decomposition on the cleaning matrix after removing the mean.

[0019] S23. Based on the singular value decomposition results, select the top K right singular vectors with a cumulative variance contribution rate of 95% as the orthogonal time basis matrix.

[0020] Furthermore, physical rule cleaning includes removing cloud shadows and false snow, SAR penetration correction, and snow and ice value filling;

[0021] Specifically, removing cloud shadows and false snow means: when the monthly average temperature of the lake's temperature time series is greater than... Furthermore, when the MNDWI index is less than -0.1, it is judged as cloud shadow or false snow noise on the shore and is removed.

[0022] SAR penetration correction specifically refers to: when the SAR water area ratio index is greater than... Furthermore, the monthly average temperature of the lake's temperature time series is greater than If MODIS high-frequency images are missing or the MNDWI index is less than 0, it should be corrected to a typical value for water bodies.

[0023] Specifically, snow and ice value filling refers to: when the monthly average temperature of the lake's temperature time series is less than... Furthermore, when MODIS high-frequency images are missing, they are filled with typical ice and snow values.

[0024] Furthermore, the expression for singular value decomposition of the cleaning matrix after removing the mean is as follows:

[0025] ;

[0026] in, This represents the MNDWI time series matrix of MODIS after physical rule cleaning, missing value interpolation, and mean removal, with dimensions of [missing information]. N represents the total number of pixels in the lake area, and T represents the time step. This represents the left singular vector matrix, used to characterize the orthogonal transformation pattern of spatial pixel directions. Let represent a right singular vector matrix, whose column vectors are used to characterize the orthogonal transformation patterns in the time direction. This indicates transpose.

[0027] Furthermore, S3 includes the following sub-steps:

[0028] S31. Extract sparse observation vectors and time basis subsets based on orthogonal time basis matrices;

[0029] S32. Construct a regularized parameter field based on the preprocessed multi-source heterogeneous data;

[0030] S33. Construct the objective function based on the sparse observation vector, the time basis subset, and the regularization parameter field;

[0031] S34. Perform weighted ridge regression on the objective function to obtain the specific time coefficients;

[0032] S35. Traverse each cell, calculate the specific time coefficient of the cell based on the cell's sparse observation vector, diagonal weight matrix, slope, offshore distance, and historical valid observation count, and perform the calculation. This yields a one-dimensional time series vector of pixels, where... Let the orthogonal time basis matrix be composed of the first K right singular vectors, and its dimension be . , This represents the reconstructed MNDWI prediction sequence of the target pixel at T time steps, with dimension . , Indicates a specific time coefficient;

[0033] S36. Map and assemble the one-dimensional time series vector of the pixels according to the original geospatial coordinates to generate the MNDWI long time series image set.

[0034] Furthermore, the regularization parameter field The expression is:

[0035] ;

[0036] in, Indicates slope, Indicates the distance from the shore. This indicates the number of historical valid observations of a pixel. This represents the basic regularization coefficient, used to control the overall strength of the ridge regression penalty term. Indicates terrain adjustment. Indicates distance adjustment. Indicates sparse penalty;

[0037] objective function The expression is:

[0038] ;

[0039] in, Represents a sparse observation vector. Represents a time-based subset. Represents an unknown time coefficient vector. Indicates transpose. This represents the diagonal weight matrix that reflects the quality of the observations;

[0040] Specific time coefficient The expression is:

[0041] ;

[0042] in, Represents the identity matrix that matches the time base dimension.

[0043] Furthermore, S4 includes the following sub-steps:

[0044] S41. Calculate robust statistics based on the MNDWI long-term image set;

[0045] S42. Based on robust statistics, determine the disaster month and rectify the disaster month;

[0046] S43. Based on the repaired disaster month, extract dynamic thresholds and construct climatological prior constraints;

[0047] S44. Based on the dynamic threshold and climatic prior constraints, the maximum water surface area is taken as the target area for spatial masking.

[0048] S45. Use the spatial mask results to extract the polygonal boundaries of the lakes and generate a lake boundary dataset.

[0049] Furthermore, robust statistics The expression is:

[0050] ;

[0051] in, This represents a set of sequences of lake water areas. Indicates the first The lake water area corresponding to each time step. This indicates that the median operation is performed on the set of variables within the parentheses.

[0052] Based on the above method, this invention also proposes a monthly-scale lake water body reconstruction system based on multi-physics field constraints and adaptive regularized spatiotemporal fusion, comprising:

[0053] The data processing module is used to collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data.

[0054] The matrix generation module is used to generate orthogonal time basis matrices based on preprocessed multi-source heterogeneous data;

[0055] The image set generation module is used to perform adaptive weighted ridge regression modeling and solving based on the orthogonal time basis matrix to generate MNDWI long time series image sets;

[0056] The lake boundary dataset generation module is used to diagnose and repair abnormal months based on the MNDWI long-term image set, and extract the monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset.

[0057] The beneficial effects of this invention are:

[0058] (1) This invention breaks through the limitations of traditional fusion models that rely solely on optical band correlation, and innovatively introduces dual physical constraints of meteorological thermodynamics and microwave penetration. ERA5 temperature data is used to establish a thermodynamic threshold for water phase change, assisting in determining freezing and thawing states and solving the problem of ice-water confusion in winter. Sentinel-1 SAR data is used to construct a water body proportion index, providing penetration repair when there is full optical cloud coverage. This mechanism allows the global time base extracted through singular value decomposition to accurately reflect the true physical evolution trend of lake water, rather than cloud and snow meteorological interference, laying a clean data foundation for high-precision reconstruction.

[0059] (2) This invention designs a dual-driven spatiotemporal fusion strategy, which captures high-frequency temporal changes in MODIS through singular value decomposition and combines it with spatial texture information from Landsat. It abandons the commonly used global or block-based unified regularization methods in traditional algorithms and innovatively constructs an adaptive parameter field that includes terrain slope, distance from shore, and observation frequency. This strategy strengthens global trend constraints in the stable lake center region to suppress noise, and releases high-frequency change details in the sensitive lake shore region to restore the morphology, thus achieving accurate characterization of the complex dynamic processes of the lake.

[0060] (3) This invention successfully breaks through the inherent limitations of the spatiotemporal resolution of a single sensor, producing refined data with both 30-meter high spatial detail and high monthly temporal frequency through the aforementioned fusion method, and extracting lake water bodies based on this. In the event of severe data gaps caused by extreme climate on the plateau, a fully automated quality control system was further established. This system introduces absolute median difference to locate disaster months without physical abrupt changes, and designs neighborhood interpolation and trend backtracking strategies to fill the spatiotemporal blind spots by utilizing the high-frequency advantage of low-resolution data. The final 25-year lake boundary range dataset not only eliminates the jagged effect of low-resolution data, but also solves the spatiotemporal gap problem of high-resolution data, ensuring the continuity and logical consistency of monitoring data in physical processes, and filling the gap in the field of refined long-term time-series monitoring of lakes on the Qinghai-Tibet Plateau. Attached Figure Description

[0061] Figure 1 The flowchart shows a method for reconstructing lunar-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion.

[0062] Figure 2A schematic diagram of a method for reconstructing lake water bodies at the lunar scale based on multiphysics constraints and adaptive regularized spatiotemporal fusion;

[0063] Figure 3 This is a comparison chart showing the spatial detail reconstruction and monthly water morphology evolution of a typical year in 2010. Figure 3 (a) Figure 3 (b) A detailed comparison of typical areas in February 2010. Figure 3 (c) is a monthly evolution diagram from January to December. Detailed Implementation

[0064] The embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0065] like Figure 1 As shown, this invention provides a method for reconstructing monthly lake water bodies based on multi-physics constraints and adaptive regularized spatiotemporal fusion, comprising the following steps:

[0066] S1. Collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data;

[0067] S2. Based on the preprocessed multi-source heterogeneous data, generate an orthogonal time basis matrix;

[0068] S3. Based on the orthogonal time basis matrix, perform adaptive weighted ridge regression modeling and solution to generate MNDWI long time series image set;

[0069] S4. Based on the MNDWI long-term image set, perform abnormal month diagnosis and repair, and extract the monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset.

[0070] like Figure 2As shown, this invention proposes a long-term monthly scale lake boundary fine reconstruction method based on multi-physics field constraints and adaptive regularization. The technical route mainly consists of four core modules: (1) Data preparation and preprocessing, integrating MODIS high-frequency images, Landsat high-resolution images, and multi-source heterogeneous data such as ERA5 temperature, Sentinel-1 SAR, and SRTM topography; (2) Multi-physics field constraint time base construction, innovatively introducing meteorological thermodynamics and microwave penetration dual constraint mechanisms to clean high-frequency time series, and extracting a global time base reflecting the true evolution trend of lakes through singular value decomposition (SVD); (3) Adaptive regularization fusion modeling, constructing an adaptive regularization parameter field including topographic slope, offshore distance and observation frequency, and using a weighted ridge regression model to achieve deep fusion of high and low resolution spatiotemporal information; (4) Water body extraction and product generation, combining disaster monthly repair based on absolute median (MAD) and climatological dynamic threshold constraints, finally generating a continuous and seamless 30-meter monthly scale lake water body dataset from 2000 to 2025.

[0071] In this embodiment of the invention, the multi-source heterogeneous data includes MODIS high-frequency imagery, Landsat high-resolution imagery, ERA5 temperature data, Sentinel-1 SAR data, and SRTM terrain data.

[0072] In this embodiment of the invention, in S1, missing values ​​are filled in for the MODIS high-frequency image, and the image is resampled according to the MNDWI index, vegetation index and snow cover index to generate a time series matrix.

[0073] The MNDWI index was used as the main variable for fusion, and the vegetation index, snow cover index and quality label band were used as auxiliary variables to determine the number of historical valid observations of a pixel.

[0074] Convert the Kelvin temperature of the ERA5 air temperature to Celsius to generate a time series of lake air temperature;

[0075] Sentinel-1 SAR data is filtered, and the filtered image is binarized into water bodies and non-water bodies using the backscattering coefficient threshold. The proportion of water body pixels is counted to generate the SAR water body area ratio index.

[0076] Slope and distance from shore are calculated based on SRTM terrain.

[0077] MODIS processing: Monthly averages were synthesized using high-quality pixels free of clouds and aerosols; if missing values ​​were found, they were filled using a lenient condition of removing only the filler values. The MNDWI index, vegetation index (NDVI), and snow cover index (NDSI) were calculated and resampled to 30m, and exported as a time series matrix in CSV format.

[0078] Landsat processing: An optimal observation quality synthesis strategy is adopted, prioritizing the selection of pixels with the lowest cloud cover and best quality for mosaicking within the monthly window, and generating cloud, snow, and shadow masks using the QA_PIXEL band; the MNDWI index is calculated as the main fusion variable, while the vegetation index (NDVI), snow index (NDSI), and quality label band are output as auxiliary variables; the historical effective observation count (ObsCount) of each pixel is counted for subsequent weight allocation.

[0079] ERA5 Processing: 2-meter air temperature data (Temperature_2m) from the ECMWF ERA5-Land reanalysis data were introduced. Using the maximum water body extent in the GSW as the region of interest, the daily average temperature within the target lake area was calculated using the reduceRegion function. Kelvin temperatures were converted to Celsius to generate daily lake temperature time series. This sequence will be used to subsequently determine the freezing and melting state of the lake.

[0080] SAR processing: Focal-median filtering with a radius of 30 meters was performed on the VV polarization band of the Sentinel-1 GRD product to remove speckle noise. A backscattering coefficient threshold was set to binarize the image into water bodies and non-water bodies. The proportion of water pixels within the lake ROI was statistically analyzed to generate a high-frequency SAR water area ratio index. .

[0081] SRTM Processing: Slope is calculated based on SRTM data. Using the Landsat reference image's Coordinate Reference System (CRS) and affine transformation matrix, the slope layer is forcibly reprojected and cropped to a 30m grid identical to Landsat's to eliminate spatial registration errors. A binary mask is constructed based on the maximum water surface area of ​​the GSW image. A fast distance transform is applied to calculate the Euclidean distance from each pixel to the nearest shore. Similarly, Landsat's CRS and Transform parameters are applied for resampling and alignment to generate the offshore distance matrix (Dist).

[0082] In this embodiment of the invention, S2 includes the following sub-steps:

[0083] S21. Using the lake temperature time series and SAR water area ratio index, the time series matrix is ​​physically cleaned to obtain the cleaned matrix.

[0084] S22. Perform singular value decomposition on the cleaning matrix after removing the mean.

[0085] S23. Based on the singular value decomposition results, select the top K right singular vectors with a cumulative variance contribution rate of 95% as the orthogonal time basis matrix.

[0086] In this embodiment of the invention, physical rule cleaning includes removing cloud shadows and false snow, SAR penetration correction, and snow and ice value filling.

[0087] Specifically, removing cloud shadows and false snow means: when the monthly average temperature of the lake's temperature time series is greater than... Furthermore, when the MNDWI index is less than -0.1, it is judged as cloud shadow or false snow noise on the shore and is removed.

[0088] SAR penetration correction specifically refers to: when the SAR water area ratio index is greater than... Furthermore, the monthly average temperature of the lake's temperature time series is greater than If MODIS high-frequency images are missing or the MNDWI index is less than 0, it should be corrected to a typical value for water bodies.

[0089] Specifically, snow and ice value filling refers to: when the monthly average temperature of the lake's temperature time series is less than... Furthermore, when MODIS high-frequency images are missing, they are filled with typical ice and snow values.

[0090] In this embodiment of the invention, the expression for singular value decomposition of the cleaning matrix after removing the mean is as follows:

[0091] ;

[0092] in, This represents the MNDWI time series matrix of MODIS after physical rule cleaning, missing value interpolation, and mean removal. Its dimension is N×T, where N represents the total number of pixels in the lake region and T represents the time step. This represents the left singular vector matrix, used to characterize the orthogonal transformation pattern of spatial pixel directions. Let represent a right singular vector matrix, whose column vectors are used to characterize the orthogonal transformation patterns in the time direction. This indicates transpose.

[0093] Constructing the initial time series matrix: Using the MODIS MNDWI time series matrix... (dimension) ,in The total number of pixels in the lake area. Import the time step into the local environment and perform physical cleaning and basis vector extraction.

[0094] Physical rule cleaning: for the original matrix The presence of cloud shadows, false snow, and missing data in the data led to the introduction of ERA5 temperature data. Ratio of SAR water bodies Perform monthly and pixel-by-pixel physical rule cleaning to generate a cleaning matrix. The cleaning rules are as follows:

[0095] Excluding cloud shadows and false snow: Monthly average temperature (Melting threshold) and At that time, noise identified as cloud shadows or false snow on the shore was removed and set to [missing information]. ;

[0096] SAR penetration correction: when and If MODIS data is missing or indicates non-water ( The value was forcibly corrected to the typical value for water bodies (0.4), and the optical omissions were corrected by utilizing microwave penetration.

[0097] Ice and snow value fill: when When the extreme cold threshold is missing and MODIS data is missing, it is forcibly filled with the typical value for ice and snow (-0.2) to maintain the continuity of the winter series.

[0098] In the above threshold settings, although the theoretical freezing point of pure water... However, considering that ERA5 data represents 2-meter air temperature rather than surface temperature, and that lake water has a large specific heat capacity, its freezing / thawing process lags behind air temperature changes. Therefore, the following settings are made: The safety buffer, Ensure the water body is in a completely liquid or dissolved state to eliminate the risk of misjudgment; To ensure the water body enters a stable freezing period, consistent with the thermodynamic freezing characteristics of lakes on the Qinghai-Tibet Plateau. Due to inherent speckle noise in SAR imagery, single-pixel interpretation is unreliable. Cross-validation with cloud-free Landsat imagery revealed that when the proportion of low-scattering (water) SAR pixels within the lake's ROI exceeds 60%, the confidence level for the area to be considered water under optical cloud cover is [not specified]. This method effectively distinguishes between scenes with full cloud cover and those with full snow cover. By analyzing the MNDWI histogram distribution of historical images of the study area, the peak values ​​for typical water bodies are between 0.4 and 0.6, while the peak values ​​for land / snow are between -0.4 and -0.2. Although cloud shadows have lower MNDWI values, they are typically distributed in the transition zone between -0.1 and 0. A cutoff threshold of -0.1 was set to maximize the removal of noise interference from cloud shadows and wet shorelines, even at the cost of sacrificing a small number of mixed pixels.

[0099] SVD time base extraction: This involves processing the matrix after physical rule cleaning and interpolation. After removing the mean, singular value decomposition (SVD) is performed.

[0100] Construction of the global time basis matrix B: Select the top K right singular vectors (K is usually 2-5 in this embodiment) with a cumulative variance contribution rate of 95% as the orthogonal time basis matrix B (dimension). This matrix captures the dominant trends in lake water volume changes with seasons and interannual variations. Each column vector in matrix B represents a dominant pattern of lake change during that period, such as seasonal fluctuations or interannual expansion trends; it is a globally known quantity in the subsequent reconstruction model.

[0101] In this embodiment of the invention, S3 includes the following sub-steps:

[0102] S31. Extract sparse observation vectors and time basis subsets based on orthogonal time basis matrices;

[0103] S32. Construct a regularized parameter field based on the preprocessed multi-source heterogeneous data;

[0104] S33. Construct the objective function based on the sparse observation vector, the time basis subset, and the regularization parameter field;

[0105] S34. Perform weighted ridge regression on the objective function to obtain the specific time coefficients;

[0106] S35. Traverse each cell, calculate the specific time coefficient of the cell based on the cell's sparse observation vector, diagonal weight matrix, slope, offshore distance, and historical valid observation count, and perform the calculation. This yields a one-dimensional time series vector of pixels, where... express, This represents the reconstructed MNDWI prediction sequence of the target pixel at T time steps, with dimension . , Indicates a specific time coefficient;

[0107] S36. Map and assemble the one-dimensional time series vector of the pixels according to the original geospatial coordinates to generate the MNDWI long time series image set.

[0108] In this embodiment of the invention, the regularization parameter field The expression is:

[0109] ;

[0110] in, Indicates slope, Indicates the distance from the shore. This indicates the number of historical valid observations of a pixel. Denotes the basic regularization coefficient, used to control the overall strength of the ridge regression penalty term. Indicates terrain adjustment. Indicates distance adjustment. Indicates sparse penalty;

[0111] objective function The expression is:

[0112] ;

[0113] in, Represents a sparse observation vector. Represents a time-based subset. Represents an unknown time coefficient vector. Indicates transpose. This represents the diagonal weight matrix that reflects the quality of the observations;

[0114] Specific time coefficient The expression is:

[0115] ;

[0116] in, Represents the identity matrix that matches the time base dimension.

[0117] This section defines the core mathematical model used to reconstruct 30-meter high-resolution imagery. This model is a pixel-by-pixel solver that does not directly generate data but defines how to calibrate the global time base using sparse Landsat observations.

[0118] Basic assumptions and sparse observation vectors ( )extract:

[0119] Assuming any 30-meter pixel Full-time variation (dimension is) It can be represented by a linear combination of the global time base matrix B:

[0120] ;

[0121] in, It is a problem to be solved. A time coefficient vector determines how a pixel responds to a global trend. However, in real-world scenarios, limited by Landsat's long revisit period and frequent cloud and snow obstruction on the Tibetan Plateau, pixel... There are only m high-quality valid observations in T months (typically) To solve for the above coefficients. We extract the MNDWI values ​​of these m valid observations from the Landsat time series data and construct a sparse observation vector in chronological order. (dimension is) Simultaneously, rows that strictly correspond in time to these m valid observations are extracted from the global time base matrix B, forming a time base subset. (dimension is) By extracting and This method successfully transforms the problem of solving the full-time evolution problem into a mathematical solution process that uses extremely sparse local known observations to fit the global trend.

[0122] Constructing an adaptive regularized parameter field:

[0123] To address the ill-conditioned equations caused by the sparsity of Landsat observations, a regularized parameter field incorporating information on topography, spatial location, and observation abundance is constructed. .

[0124] In areas with significant topographic relief, the geometric correction residuals at 30-meter resolution typically increase exponentially with increasing slope. Experience shows that slopes greater than... The probability of sub-pixel-level misalignment increases significantly in certain regions, therefore high regularization weights are assigned to suppress spatial errors; while Registration accuracy is high in flat areas, and weights are reduced to preserve details. Therefore, in terrain adjustment... When the slope At that time, the coefficient was set to 2.0 to suppress terrain registration errors; slope The time factor is set to 0.7. Distance adjustment... At that time, when the distance from the shore (In deep water or on land) the coefficient is set to 2.0, forcibly following the time base trend; distance (For lakeshore variation areas), the coefficient is set to 0.5 to preserve high-frequency details. Sparsity penalty. The number of valid historical observations for Landsat is extremely low ( The number of pixels increases nonlinearly. To prevent overfitting.

[0125] Solving weighted ridge regression:

[0126] To solve for the unknown time coefficient vector This invention constructs a weighted ridge regression objective function. This function aims to minimize the weighted residuals between Landsat sparse observations and model predictions, while introducing a regularization penalty term to constrain model complexity and prevent overfitting. Its formula is defined as follows:

[0127] ;

[0128] in It is the time basis subset corresponding to the observation time. It is a diagonal weight matrix that reflects the observation quality (high-quality effective observation pixels have a weight of 1, pixels with slight pollution from clouds and snow have a correspondingly reduced weight, and severely occluded pixels have a weight of 0). This is the adaptive regularization parameter specifically designed for this pixel, calculated in the previous step.

[0129] Solve for the time coefficient a:

[0130] In order to obtain the objective function The global optimal solution for the above objective function with respect to the time coefficient vector Find the partial derivative and set it to zero. After algebraic simplification, derive the analytical solution (closed-form solution) to the weighted ridge regression problem. Using this closed-form solution, the algorithm can efficiently and stably compute any effective pixel within the study area. Dedicated time coefficient vector This allows for the calibration of parameters for global temporal variation trends based on local sparse observations.

[0131] Full-time reconstruction:

[0132] Traverse each cell within the reconstructed region and read the Landsat observation vector at that location. Weight And the corresponding physical parameters (slope, distance, observation count). Substitute these into the closed-form formula in step 4 to calculate the specific time coefficient for that pixel. Perform the operation. Through this matrix multiplication operation, a continuous MNDWI prediction value for the pixel over a period of 300 months (i.e., covering every month within the study period from 2000 to 2025) is obtained. This step successfully transforms a single pixel from incomplete, locally sparse observations into a complete and continuous one-dimensional time series vector covering the entire study period.

[0133] Generation of 30m MNDWI long-term imagery:

[0134] After predicting all valid pixels, the spatial dimensions are assembled and output. The system assembles and outputs massive amounts of one-dimensional continuous time series vectors according to their original geospatial coordinates. The images are then remapped and reassembled. The predicted values ​​of all pixels are reorganized according to their spatial location, and the final output is a continuous, seamless monthly 30-meter spatial resolution MNDWI long-term image set from 2000 to 2025.

[0135] In this embodiment of the invention, S4 includes the following sub-steps:

[0136] S41. Calculate robust statistics based on the MNDWI long-term image set;

[0137] S42. Based on robust statistics, determine the disaster month and rectify the disaster month;

[0138] S43. Based on the repaired disaster month, extract dynamic thresholds and construct climatological prior constraints;

[0139] S44. Based on the dynamic threshold and climatic prior constraints, the maximum water surface area is taken as the target area for spatial masking.

[0140] S45. Use the spatial mask results to extract the polygonal boundaries of the lakes and generate a lake boundary dataset.

[0141] In this embodiment of the invention, robust statistics The expression is:

[0142] ;

[0143] in, This represents a set of sequences of lake water areas. Indicates the first The lake water area corresponding to each time step. This indicates that the median operation is performed on the set of variables within the parentheses.

[0144] Based on the above method, this invention also proposes a monthly-scale lake water body reconstruction system based on multi-physics field constraints and adaptive regularized spatiotemporal fusion, comprising:

[0145] The data processing module is used to collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data.

[0146] The matrix generation module is used to generate orthogonal time basis matrices based on preprocessed multi-source heterogeneous data;

[0147] The image set generation module is used to perform adaptive weighted ridge regression modeling and solving based on the orthogonal time basis matrix to generate MNDWI long time series image sets;

[0148] The lake boundary dataset generation module is used to diagnose and repair abnormal months based on the MNDWI long-term image set, and extract the monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset.

[0149] This step describes how to apply the long-term temporal images generated by the above model to the final vector product extraction and fully automated quality control.

[0150] Disaster Month Diagnosis:

[0151] For disaster monthly diagnosis of reconstructed images, if the water area in a certain month shows a sudden jump compared to the adjacent month, the median absolute difference (MAD) is introduced as a robust statistic, and its calculation formula is as follows:

[0152] ;

[0153] in Given a set of lake water area sequences, when the lake area changes at adjacent time steps satisfy... or effective observation ratio Furthermore, severe snow pollution triggers abnormal mutation conditions, thus classifying it as a disaster month; otherwise, it is classified as normal imagery.

[0154] Handling of abnormal mutations:

[0155] For disaster months identified as "abnormal mutations" by MAD (Modular Image Deposition) diagnostics, an automated restoration strategy is implemented. Priority is given to using neighborhood interpolation results from adjacent normal months; if continuous temporal gaps prevent interpolation, model backtracking is performed, using pure prediction results (PRED) reconstructed solely from the MODIS time base to fill in the gaps, ensuring the reconstructed image sequence does not violate physical laws. The restored images are then incorporated into the subsequent water body extraction process along with the normal images.

[0156] Dynamic threshold extraction:

[0157] For both time-series and restored images, pixel-by-pixel water segmentation is performed. For each single image, the algorithm dynamically calculates its local optimal segmentation threshold (dynamic threshold extraction) to adapt to changes in the spectral characteristics of water bodies under different seasons and lighting conditions.

[0158] Climatic threshold constraint (Otsu):

[0159] To prevent the dynamic threshold from failing due to local noise, a climatological prior constraint is introduced. High-quality observation months (observation ratio) from 2000–2021 are utilized. And the integration is efficient ) Calculate the median of the optimal Otsu threshold for each month to establish a monthly threshold climatology. For a single-period image, calculate the local Otsu threshold. and forcibly restrict to Within the range.

[0160] GSW spatial mask:

[0161] After threshold segmentation, the maximum water surface area of ​​the Global Surface Water Data Set (GSW) is used as the target area for spatial masking. This step effectively eliminates the interference of mountain shadows, topographic relief, and other non-water features around the lake on the water extraction results.

[0162] Vectorization and Product Generation:

[0163] Vectorization was performed on the binary image after spatial masking to extract smooth and continuous lake polygon boundaries. The final output generated a "Monthly 30m Lake Boundary Dataset from 2000 to 2025", filling the gap in the field of refined long-term time-series monitoring of lakes on the Qinghai-Tibet Plateau.

[0164] Through all the steps described above, this embodiment successfully generated a monthly 30-meter lake water body range vector dataset for Namtso from 2000 to 2025. The reconstruction results accurately reflect the dynamic changes of the lake during the dry season, wet season, and freezing season. In the cloudless months, it is highly consistent with the original Landsat observations, and in the fully clouded months, it provides reliable estimation results through physical constraints.

[0165] The accuracy of the reconstruction results was verified using the retained validation pixels. The root mean square error (RMSE) and coefficient of determination (R²) for the entire time series from 2000 to 2025 were statistically analyzed. 2 The results are presented every five years. As shown in Table 1, due to the complex meteorological conditions of the Tibetan Plateau and the Landsat sensor stripe loss (SLC-off), the average annual effective observation rate in the study area is only about 23.8%, and in some years (such as 2005) it is even less than 15%. Despite the extremely scarce data, in the adaptive regularization model constructed in this invention, the R of the Namtso region... 2 The values ​​generally remained between 0.781 and 0.878, indicating a good model fit. The RMSE values ​​ranged from 0.113 to 0.157, with relatively small fluctuations, indicating a generally low error level and no significant accuracy decay in years with extremely missing data. Overall, the model exhibits good stability across different years, with high overall fit, small error, and strong temporal consistency and reliability.

[0166] Table 1

[0167] time Valid observation ratio (%) <![CDATA[R 2 ]]> RMSE 2000 20.6 0.845 0.135 2005 13.6 0.781 0.157 2010 17.2 0.834 0.145 2015 29.2 0.878 0.113 2020 25.3 0.861 0.121 2025 28.1 0.869 0.119

[0168] To further verify the model's spatial reconstruction capabilities under different seasons and meteorological backgrounds, this study selected 2010 as a typical test year to comprehensively assess the spatial morphological details and intra-year evolution of Namtso Lake. Figure 3 (a) and Figure 3 (b) Shows a comparison of the spatial detail reconstruction results of a typical area in February 2010. Figure 3 (c) shows the complete monthly evolution of water morphology from January to December.

[0169] In terms of spatial detail reconstruction, due to the limitations of sensor spatial resolution, the original MODIS image (500m) exhibits severe pixel mixing at the lake edge, and the extracted water boundary (golden solid line) shows obvious "jaggedness," making it difficult to accurately depict the tortuous details of the lake shoreline. Figure 3(a)). In contrast, our method (MPAR-STF) fully utilizes the 30m spatial texture information of Landsat to effectively downscale and reconstruct the low-resolution signal. The reconstructed water body boundary (red solid line) is smooth and continuous, accurately conforming to the complex shoreline morphology of Namtso Lake, and significantly eliminating the spatial discretization error of the low-resolution data. Figure 3 (b)). Regarding the spatiotemporal evolution over the year ( Figure 3 (c) During the freezing periods from January to March and November to December, MODIS observations based solely on optical methods are prone to confusing floating ice, snow cover, and water bodies, leading to drastic, non-physical fluctuations in the extracted range. This method introduces a multi-physics constraint mechanism using air temperature (ERA5) and microwave (SAR) data, effectively eliminating false snow noise and ensuring that the reconstructed boundary strictly follows the thermodynamic and structural characteristics of water ice, maintaining the physical authenticity of the lake outline in winter. During the rainy season from June to September, cloud cover in the plateau region is extremely frequent, causing boundary breaks in optical observations due to shadows or thick cloud cover. Comparative results show that even with cloud contamination or missing bands in the Landsat base imagery for some months, the red boundary calculated by this method still maintains extremely high integrity. This is thanks to the global trend information extracted from the SVD time base, enabling the model to penetrate local noise and extrapolate missing boundaries.

[0170] In terms of spatial accuracy, compared to the jagged boundaries of MODIS which suffer from severe mixed pixel effects, this invention makes full use of 30 meters of spatial texture information for reconstruction, significantly reducing the boundary displacement error from 199.80 meters to 39.79 meters, and achieving high-precision and fine restoration of lake shore details.

[0171] In terms of anti-interference, during the winter freezing period, the invention introduces a multi-physics field constraint mechanism of air temperature and microwave to effectively eliminate false snow noise, overcome the defect of easy confusion between ice and water by single optical observation, and maintain the authenticity of the water body boundary contour.

[0172] In terms of temporal continuity, when faced with image gaps caused by high-frequency thick clouds during the rainy season or missing sensor strips, such as in March and August, this method relies on the global time base trend to penetrate local noise, accurately deduce and repair the spatiotemporal gaps, and ensure the high continuity and integrity of the long-term evolution process.

[0173] Quantitative evaluation shows that, based on the complete 30m Landsat original images from February, April, November, and December, compared with the original MODIS extraction results, the reconstruction results of this method have improved the annual effective intersection-over-union (IoU) from 86.37% to 94.54%, and the average boundary displacement error (BDE) has significantly decreased from 199.80 meters (approximately 0.5 MODIS pixels) to 39.79 meters (approximately 1 Landsat pixel), effectively eliminating the spatial discrepancy error of low-resolution data, as shown in Table 2.

[0174] Table 2

[0175] month MODIS_IoU(%) Ours_IoU(%) MODIS_BDE(m) Ours_BDE(m) 2 86.24 95.33 198.72 39.82 4 85.04 92.75 212.38 49.58 11 86.97 94.79 191.30 35.06 12 87.22 95.27 196.77 34.68 average value 86.37 94.54 199.80 39.79

[0176] Compared with existing technologies, this invention, based on the GEE cloud platform and local high-performance computing environment, deeply integrates heterogeneous data from multiple sources such as optical, meteorological, and microwave data. It constructs a multi-physics-constrained time-based cleaning model and an adaptive regularized spatiotemporal fusion algorithm (MPAR-STF), and combines disaster monthly restoration and climatological threshold constraints for post-processing, achieving seamless reconstruction of lake water bodies at a 30-meter monthly scale from 2000 to 2025. The high spatiotemporal resolution reconstructed dataset effectively overcomes the spatiotemporal gaps caused by cloud and snow obstruction and scarce observations in plateau regions. Especially under extreme observation conditions such as ice-water mixing, cloud shadow interference, and missing sensor bands, it maintains the consistency of physical laws and accurately restores lake shoreline details. Furthermore, it is of great significance in revealing the response mechanism of lakes to climate change, assessing water resource security in the service area, and supporting the simulation of hydrological processes in high-altitude and cold environments.

[0177] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A method for reconstructing lunar-scale lake water bodies based on multi-physics constraints and adaptive regularized spatiotemporal fusion, characterized in that, Includes the following steps: S1. Collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data; S2. Based on the preprocessed multi-source heterogeneous data, generate an orthogonal time basis matrix; S3. Based on the orthogonal time basis matrix, perform adaptive weighted ridge regression modeling and solution to generate MNDWI long time series image set; S4. Based on the MNDWI long-term image set, perform abnormal month diagnosis and repair, and extract the monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset.

2. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 1, characterized in that, The multi-source heterogeneous data includes MODIS high-frequency imagery, Landsat high-resolution imagery, ERA5 temperature data, Sentinel-1 SAR data, and SRTM topography.

3. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 1, characterized in that, In step S1, missing values ​​are filled in for the MODIS high-frequency image, and the image is resampled according to the MNDWI index, vegetation index, and snow cover index to generate a time series matrix. The MNDWI index was used as the main variable for fusion, and the vegetation index, snow cover index and quality label band were used as auxiliary variables to determine the number of historical valid observations of a pixel. Convert the Kelvin temperature of the ERA5 air temperature to Celsius to generate a time series of lake air temperature; Sentinel-1 SAR data is filtered, and the filtered image is binarized into water bodies and non-water bodies using the backscattering coefficient threshold. The proportion of water body pixels is counted to generate the SAR water body area ratio index. Slope and distance from shore are calculated based on SRTM terrain.

4. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 1, characterized in that, S2 includes the following sub-steps: S21. Using the lake temperature time series and SAR water area ratio index, the time series matrix is ​​physically cleaned to obtain the cleaned matrix. S22. Perform singular value decomposition on the cleaning matrix after removing the mean. S23. Based on the singular value decomposition results, select the top K right singular vectors with a cumulative variance contribution rate of 95% as the orthogonal time basis matrix.

5. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 4, characterized in that, The physical rule cleaning includes removing cloud shadows and false snow, SAR penetration correction, and snow and ice value filling; The removal of cloud shadows and false snow specifically refers to: when the monthly average temperature of the lake's temperature time series is greater than... Furthermore, when the MNDWI index is less than -0.1, it is judged as cloud shadow or false snow noise on the shore and is removed. The SAR penetration correction specifically refers to: when the SAR water area ratio index is greater than... Furthermore, the monthly average temperature of the lake's temperature time series is greater than If MODIS high-frequency images are missing or the MNDWI index is less than 0, it should be corrected to a typical value for water bodies. The specific method for filling in the ice and snow value is as follows: when the monthly average temperature of the lake's temperature time series is less than... Furthermore, when MODIS high-frequency images are missing, they are filled with typical ice and snow values.

6. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 4, characterized in that, The expression for performing singular value decomposition on the cleaning matrix after removing the mean is as follows: ； in, This represents the MNDWI time series matrix of MODIS after physical rule cleaning, missing value interpolation, and mean removal. Describes a left singular vector matrix. Describes a right singular vector matrix. This indicates transpose.

7. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 1, characterized in that, S3 includes the following sub-steps: S31. Extract sparse observation vectors and time base subsets based on the orthogonal time base matrix; S32. Construct a regularized parameter field based on the preprocessed multi-source heterogeneous data; S33. Construct the objective function based on the sparse observation vector, the time basis subset, and the regularization parameter field; S34. Perform weighted ridge regression on the objective function to obtain the specific time coefficients; S35. Traverse each cell, calculate the specific time coefficient of the cell based on the cell's sparse observation vector, diagonal weight matrix, slope, offshore distance, and historical valid observation count, and perform the calculation. This yields a one-dimensional time series vector of pixels, where... Let represent the orthogonal time basis matrix formed by the first K right singular vectors. This represents the reconstructed MNDWI prediction sequence of the target pixel at T time steps. Indicates a specific time coefficient; S36. Map and assemble the one-dimensional time series vector of the pixels according to the original geospatial coordinates to generate the MNDWI long time series image set.

8. The method for reconstructing monthly-scale lake water bodies based on multiphysics constraints and adaptive regularized spatiotemporal fusion according to claim 7, characterized in that, The regularization parameter field The expression is: ； in, Indicates slope, Indicates the distance from the shore. This indicates the number of historical valid observations of a pixel. Represents the basic regularization coefficient. Indicates terrain adjustment. Indicates distance adjustment. Indicates sparse penalty; The objective function The expression is: ； in, Represents a sparse observation vector. Represents a time-based subset. Represents an unknown time coefficient vector. Indicates transpose. This represents the diagonal weight matrix that reflects the quality of the observations; The specific time coefficient The expression is: ； in, Represents the identity matrix that matches the time base dimension.

9. The method for reconstructing monthly-scale lake water bodies based on multi-physics constraints and adaptive regularized spatiotemporal fusion according to claim 1, characterized in that, S4 includes the following sub-steps: S41. Calculate robust statistics based on the MNDWI long-term image set; S42. Based on robust statistics, determine the disaster month and rectify the disaster month; S43. Based on the repaired disaster month, extract dynamic thresholds and construct climatological prior constraints; S44. Based on the dynamic threshold and climatic prior constraints, the maximum water surface area is taken as the target area for spatial masking. S45. Use the spatial mask results to extract the polygonal boundaries of the lakes and generate a lake boundary dataset; The robust statistics The expression is: ； in, This represents a set of sequences of lake water areas. Indicates the first The lake water area corresponding to each time step. This indicates that the median operation is performed on the set of variables within the parentheses.

10. A lunar-scale lake water body reconstruction system based on multi-physics constraints and adaptive regularized spatiotemporal fusion, characterized in that, include: The data processing module is used to collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data. The matrix generation module is used to generate orthogonal time basis matrices based on preprocessed multi-source heterogeneous data; The image set generation module is used to perform adaptive weighted ridge regression modeling and solving based on the orthogonal time basis matrix to generate MNDWI long time series image sets; The lake boundary dataset generation module is used to diagnose and repair abnormal months based on the MNDWI long-term image set, and extract the monthly lake water body boundaries based on the repaired MNDWI long-term image set to generate a lake boundary dataset.

Images