Water depth retrieval method based on SDW-LPT multi-temporal fusion remote sensing image

Abstract

The application discloses a water depth inversion method based on SDW-LPT multi-temporal fusion remote sensing image, and particularly relates to the field of water depth inversion. The original image of Sentinel-2 is subjected to radiation correction to eliminate system error, the filter window is adjusted according to the registration error, the filter weight is optimized in combination with the window size and different band characteristics, the image preprocessing is completed, and NDWI is calculated. Based on NDWI, the preprocessing image is subjected to cloud mask and land mask operation, ICESat-2 photon points are projected to the coordinate system consistent with the Sentinel-2 image, isolated noise photons are removed through DBSCAN density clustering, and underwater terrain photons are extracted, the standardized reference water depth aligned with each pixel is obtained through refraction and tide correction, and the training set and the verification set are divided according to the spatial distribution and the water depth range. The SDW-LPT method is used for multi-temporal image fusion, and multiple models are used for water depth inversion.

Description

Technical Field

[0001] This invention relates to the field of water depth inversion technology, specifically to a water depth inversion method based on SDW-LPT multi-temporal fusion remote sensing images. Background Technology

[0002] High-precision underwater topographic data of islands and reefs is a fundamental element supporting various marine engineering projects, including marine spatial planning, resource development, navigation safety, and ecological protection. Currently, mainstream depth acquisition methods mainly include active remote sensing technologies such as shipborne multibeam sonar and airborne lidar bathymetry systems. While these technologies offer high accuracy and active acquisition capabilities, their high cost and limited operating environments make them unsuitable for large-scale, high-frequency dynamic depth updates. In contrast, satellite depth inversion (SDB) technology based on optical remote sensing imagery offers advantages such as wide coverage, low cost, and short cycle time, becoming an effective supplement to traditional bathymetry methods and demonstrating great potential in improving the efficiency and spatial continuity of depth data acquisition. Among these, the Sentinel-2 series satellites, with their high spatial resolution, multispectral imaging capabilities, and 5-day revisit period, have become one of the most representative and widely used data sources for SDB.

[0003] However, the practical application of SDB still faces significant challenges. Single-temporal Sentinel-2 images are often affected by disturbances such as clouds and waves, leading to texture occlusion, noise, or outliers in water areas, which severely impacts the accuracy of SDB and the stability of water depth results. Furthermore, traditional strategies for manually selecting the best images are highly subjective and lack timeliness, making it difficult to support the needs of large-scale, high-frequency ocean bathymetry. Based on this, recent research has begun to focus on multi-temporal remote sensing image fusion strategies, utilizing the complementarity of time-series data in dimensions such as observation angle, wave state, and cloud cover to offset random noise and improve data quality, thereby enhancing the robustness of SDB. Existing techniques such as median synthesis and maximum outlier removal have improved inversion accuracy to some extent, providing new ideas for shallow water depth inversion. However, the effectiveness of these methods may be limited when image noise has spatial clustering characteristics or the overall quality of multi-temporal images is poor. Summary of the Invention

[0004] To address this, the present invention provides a water depth inversion method based on SDW-LPT multi-temporal fusion remote sensing images. This method integrates the image's multi-scale structural modeling capability with a local feature weighting mechanism, effectively suppressing noise interference while preserving spatial details, thereby improving SDB accuracy and solving the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a water depth inversion method based on SDW-LPT multi-temporal fusion remote sensing images, comprising the following steps:

[0006] The acquired Sentinel-2 raw images were radiometrically corrected to eliminate systematic errors. The filter window was then adjusted according to the registration error to adapt to the heterogeneity of the error. The filter weights were adjusted according to the characteristics of different bands based on the window size to balance noise reduction and preservation of ground feature details. This resulted in preprocessed Sentinel-2 raw images with consistent radiometric calibration, improved registration accuracy, and preservation of ground feature details, as well as Normalized Differential Water Index (NDWI) with high comparability in the same water body under multiple time phases and scenarios.

[0007] Combining the high reflectivity of clouds with the texture differences of land areas, cloud masking and land masking were performed on the preprocessed Sentinel-2 raw images based on the Normalized Differential Water Index (NDWI) to provide a geometrically consistent effective water range for subsequent multi-temporal fusion and ICESat-2 photon matching.

[0008] The ICESat-2 photon points are projected onto a coordinate system consistent with the preprocessed Sentinel-2 image. Each photon is associated with the water / non-water attributes and local texture information of the corresponding pixel. The DBSCAN density clustering method is used to perform spatial-height joint clustering of ICESat-2 photon heights. Based on the water connectivity and texture complexity provided by the preprocessed image (e.g., reflected by entropy and local variance), the neighborhood radius and minimum sample number of DBSCAN are adaptively set. This allows isolated noise photons to be preferentially removed in low-texture, low-noise water areas, while enhancing the retention of real underwater topographic photons in complex nearshore areas. The pixel-photon spatial correspondence brought by high registration accuracy is used to reduce pseudo-clustering caused by photons being misprojected onto land or the outside of the shoreline, thereby obtaining underwater topographic photons. Refraction correction and tidal level correction are then performed to obtain a reference water depth that is strictly aligned and standardized pixel by pixel with the preprocessed Sentinel-2 image. The reference water depth is divided into a model training set and a validation set according to the spatial distribution and water depth range.

[0009] The SDW-LPT method was used to fuse multi-temporal images to generate multi-temporal fused images. Based on the multi-temporal fused image data, water depth inversion was performed using a multi-band model, a logarithmic ratio model, and a quadratic polynomial ratio model (QPRM).

[0010] Using ICESat-2 validation data as a benchmark, the inversion results of the multi-band model, logarithmic ratio model, and quadratic polynomial ratio model (QPRM) are analyzed. Statistical indicators such as root mean square error (RMSE), systematic bias (Bias), and coefficient of determination (R²) are calculated to evaluate the performance of the fusion method.

[0011] Preferably, the steps for acquiring Sentinel-2 raw images are as follows: extracting four key band data (blue, green, red, and near-infrared) from the Sentinel-2 images, and simultaneously extracting image RPC parameters or ground control points (GCPs) for subsequent registration error calculation.

[0012] Preferably, the radiometric correction steps for the Sentinel-2 raw image are as follows:

[0013] System errors in the blue, green, red, and near-infrared bands of the image were corrected, and the corrected pixel values ​​are as follows:

[0014] Corrected pixel value = Original pixel value × Band gain coefficient (blue 1.05, green 1.02, red 1.01, near-infrared 0.98) - Band offset (obtained by statistically analyzing the black border pixels in the image);

[0015] For each pixel to be processed, the variance σ² is calculated in the 3×3 candidate window, and then compared with the adaptive noise threshold to classify the noise level; if it is classified as strong noise, then hybrid filtering (median + weighted mean) is applied; if it is weak noise, then only weighted mean filtering is required.

[0016] Preferably, the dynamic adjustment logic for the adaptive noise threshold is as follows:

[0017] Using the blue, green, red, and near-infrared bands of Sentinel-2 imagery as units, the pixel grayscale variance distribution of the entire image was statistically analyzed to determine the baseline threshold. The specific steps are as follows:

[0018] Calculate the global baseline threshold for Sentinel-2 images:

[0019] Divide the single-band image into blocks (e.g., 256×256 pixel blocks) and calculate the variance of each block. ;

[0020] variance of all blocks Perform statistical analysis to create a variance dataset and sort it from smallest to largest. Take the 75th percentile of the variance distribution (i.e., 75% of the block variances are less than or equal to this value, and 25% of the block variances are greater than or equal to this value) as the baseline threshold. ;

[0021] Threshold adjustment based on band characteristics: Introducing band weighting coefficients based on band noise sensitivity. Adjust the formula: The band weighting coefficients are obtained through experimental or statistical analysis, taking into account the noise characteristics of each band (blue / green / red / near-infrared) of Sentinel-2.

[0022] Within the same band, further spatial adaptive adjustments are made to address the differences in texture variance across different land cover areas (vegetation, water bodies, buildings) to obtain pixel-level local thresholds. The specific steps are as follows:

[0023] For each pixel's 3×3 neighborhood, calculate the texture complexity of the features within that neighborhood. ;

[0024] Set texture complexity coefficient ( Complex texture areas (For example, vegetation has a high entropy value) (lower the threshold to avoid misclassifying textures as noise); simple texture areas (such as water bodies) have a low entropy value. (Increase the threshold to enhance noise suppression), adjust the formula: Its texture complexity coefficient The value is determined by the complexity of the terrain texture.

[0025] Preferably, the dynamic adjustment steps for the filter window based on registration error are as follows:

[0026] Using RPC parameters or ground control points (GCP), a mapping relationship between image pixel coordinates (row and column numbers) and real geographic coordinates (latitude and longitude / UTM) is established. The geographic coordinates of the image pixel to be corrected are compared with the geographic coordinates of the reference base (reference image / real ground features). The difference is transformed by projection to calculate the registration error (Δx, Δy) of the current pixel. The error level is quantized based on the adaptive error threshold, and the window parameters are dynamically determined according to the error level.

[0027] Preferably, the adaptive error threshold adjustment logic is as follows:

[0028] First, the global registration error distribution of the image to be processed is statistically analyzed to determine the baseline range of the threshold:

[0029] Introducing scene weight coefficients Adjust the baseline threshold and correct the formula according to the application scenario: Among them, the scene weight coefficient The range to be set depends on the application scenario. Pixel-level local threshold Application scenarios include fine-scale ground feature classification / contour extraction, coarse-scale monitoring / coverage estimation, and routine image preprocessing;

[0030] To address the registration error sensitivity of different bands, a band correction coefficient is introduced. Fine-tuning yields the final adaptive error threshold: Among them, band correction coefficient The settings are based on the Sentinel-2 four-band characteristics.

[0031] Preferably, the steps for adjusting the filter weights based on the characteristics of different frequency bands are as follows:

[0032] Based on the window size, a preset basic weight matrix for the bands is defined as follows:

[0033] Based on the weight allocation principle, the sum of all elements in the weight matrix of each band is set to 1 (ensuring that the filtered pixel values ​​are within a reasonable dynamic range). The formula is: , Represents the original weights. This represents the sum of the original weights. The weight matrix (the sum of all elements is 1) represents the normalized weights and is ultimately used for filtering calculations. It is an unnormalized weight matrix preset based on band characteristics and window levels. The window levels are defined, and the rationality of the basic weight matrix is ​​preset and verified for different bands based on different window levels.

[0034] Filtering is performed according to the window that divides the noise level; the mean deviation of radiation, PSNR, and mean registration error before and after filtering are calculated. If any one of the following is true: mean deviation of radiation before and after filtering is greater than or equal to 5%, PSNR is less than 3dB, or the mean reduction of registration error is less than 30%, then the standard is not met, and the parameters are adjusted backtracked.

[0035] Preferably, a weighted gradient strategy is used in the boundary regions of different error levels and different bands to avoid image stitching artifacts; the specific settings of the weighted gradient strategy are as follows:

[0036] Determine the width of the boundary region: Transition region width = 1 / 2 of the window size difference, rounded up;

[0037] The basis for determining gradient weights: the weight matrix for each pixel within the transition region. Weight matrix on one side of the boundary and the weight matrix on the other side of the boundary Calculated using a distance-weighted average;

[0038] Quantization gradient: For each pixel in the transition region, calculate the gradient coefficient t, and then obtain the final weight using a linear gradient formula. The specific operation is as follows:

[0039] (1) Calculation of the gradual change coefficient t: ;

[0040] in, This represents the distance from the current pixel to the boundary of the reference region A, within the transition zone. From 0 to ; It is the width of the transition zone;

[0041] (2) The final weight matrix of pixels in the transition region, through... and Linear interpolation yields:

[0042] ;

[0043] in, These are the row and column indices of the weight matrix, calculated for... Renormalize to ensure that the sum of all elements is 1.

[0044] Preferably, the specific steps for fusing multi-temporal images using the SDW-LPT method are as follows:

[0045] Each preprocessed image is decomposed into a three-layer Laplacian pyramid and then filtered using a Gaussian filter. The low-frequency components of each layer are extracted, and the high-frequency detail map of the preprocessed image is obtained by removing the low-frequency components of each layer. The formula is as follows:

[0046] ;

[0047] In the formula, The high-frequency components of the Laplacian at layer n are represented by the high-frequency components of each layer of the pyramid. Since the pre-processed image is a denoised but detail-preserving multi-band image, the high-frequency components of each layer of the pyramid reflect the texture of the seabed and subtle changes in the light of the water, rather than random noise. This is the original input image for layer 𝑖; It is a Gaussian filter. , ;

[0048] To enhance the representation of complex textured regions, a local energy map is constructed in each image to reflect the degree of local variation. This energy map serves as a weighting factor for superimposing high-frequency information from multi-temporal images, including the signal-to-noise ratio obtained from radiometric correction and band-weighted filtering, local texture complexity, and consistency with the reference water depth. The formula is as follows:

[0049] ;

[0050] In the formula, This is the normalized energy weights graph. Indicated by The standard deviation of the 5×5 sliding window centered on the center. The maximum standard deviation for each graph;

[0051] Using the lowest-level low-frequency image of the Laplacian pyramid as the fusion base, high-frequency components are weighted and superimposed layer by layer to obtain the final fused image. The reconstruction formula is as follows:

[0052] ;

[0053] In the formula, For the final merged image; The lowest frequency image extracted from the Lth layer of the pyramid; For the high-frequency components of the nth layer; It is a normalized energy map; the constant 0.5 is a smoothing factor, which aims to improve the fusion contribution in low-energy regions and prevent the loss of details.

[0054] The preferred multi-band model is expressed as follows:

[0055] ;

[0056] In the formula, Because of the water depth, and For the parameter to be estimated, For the number of bands, Indicates the first The water reflectance corresponding to each waveband Indicates the deep water zone is in the 1st Reflectance in each wavelength band; , respectively, are the reflectance of the blue, green and red bands; m and n are constant terms, and the logarithmic adjustment factor is adjusted to two factors; this model is widely used in remote sensing water depth inversion by integrating the reflection of seabed information in different bands;

[0057] The logarithmic ratio model, mathematically expressed as follows:

[0058] ;

[0059] In the formula, Because of the water depth, and Here, n is the parameter to be estimated, and n is a fixed constant. In the logarithmic ratio model, n = 1000. Indicates the first The water reflectance corresponding to each band; since blue and green light have strong penetrating power in water, this model usually uses the blue and green band reflectance for calculation.

[0060] The quadratic polynomial ratio model is used for water depth estimation. First, the water depth index W is calculated using the following formula:

[0061] ;

[0062] In the formula, The depth index is used in the quadratic polynomial ratio model, where n is taken as 1500. , The water reflectance is defined for the blue and green bands, respectively; subsequently, a quadratic polynomial ratio model is established:

[0063] ;

[0064] In the formula, Because of the water depth, Let W be the median value, and a, b, and c be the coefficients to be estimated, which can be obtained by fitting using the least squares method.

[0065] This invention not only improves the quality of a single Sentinel-2 image through radiometric correction, registration error-driven window adjustment, and band-wise filter weight optimization, but also extends its advantages to the entire process of NDWI mask extraction, ICESat-2 photon denoising and depth correction, SDW-LPT multi-temporal fusion, and multi-model water depth inversion.

[0066] 1. Detail fidelity and radiation consistency improve the recognizability of the water-land boundary and shoal texture, directly improving the accuracy of NDWI mask and nearshore water depth inversion;

[0067] 2. Improved registration accuracy ensures spatial consistency between photons, pixels, and multiple time phases, reducing the "spurious dispersion" of training samples in the spectrum-depth relationship;

[0068] 3. The adaptive adjustment of the band filtering weights enables the blue and green water depth sensitive bands to retain the signals most sensitive to water depth changes while suppressing noise, providing a higher signal-to-noise ratio input for the SDW-LPT fusion weights and the subsequent three types of inversion models. Attached Figure Description

[0069] Figure 1 The inversion results of the fused SDW-LPT inversion method and multi-band model provided by this invention in three types of fused images;

[0070] Figure 2 The inversion results of the fusion median composite method and multi-band model provided by this invention in three types of fused images;

[0071] Figure 3 The inversion results of the fusion maximum outlier removal method and the multi-band model in three types of fused images provided by this invention;

[0072] Figure 4 The inversion results of the fused SDW-LPT inversion method and log ratio model provided by this invention in three types of fused images;

[0073] Figure 5 The inversion results of the fusion median synthesis method and log ratio model provided by this invention in three types of fused images;

[0074] Figure 6The inversion results of the maximum outlier removal method and the log ratio model provided by this invention in three types of fused images;

[0075] Figure 7 The inversion results of the fused SDW-LPT inversion method and quadratic polynomial ratio model provided by this invention in three types of fused images;

[0076] Figure 8 The inversion results of the fusion median synthesis method and the quadratic polynomial ratio model provided by this invention in three types of fused images;

[0077] Figure 9 The inversion results of the maximum outlier removal method and the quadratic polynomial ratio model provided by this invention in three types of fused images;

[0078] Figure 10 Comparison of RMSE for the three fusion methods provided by this invention under different water depth inversion models. Detailed Implementation

[0079] The following specific embodiments illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0080] To effectively improve the accuracy and stability of satellite remote sensing water depth inversion in typical shallow water areas of islands and reefs, this embodiment, based on the technical solution of the present invention, selects Sentinel-2 multispectral imagery and ICESat-2 ATL03 photon altimetry data (LI YU, ZHOU LEI, JIA DONGZHEN, et al. 2025. A denoising algorithm for ICESat-2 bathymetric photons based on point cloud gridding[J]. IEEE Transactions on Geoscience and Remote Sensing, 63: 1–22.) and selects a study area. The applicability, accuracy, and robustness of the SDW-LPT method of the present invention and existing methods under different inversion models are compared and analyzed.

[0081] This invention proposes a water depth inversion method based on SDW-LPT multi-temporal fusion remote sensing images, comprising the following steps:

[0082] Step 1: Extract the four key bands of blue, green, red, and near-infrared from the Sentinel-2 image, and simultaneously extract the image RPC parameters or ground control points (GCPs) for subsequent registration error calculation; eliminate systematic errors by radiometric correction on the acquired Sentinel-2 raw image, and then adjust the filter window according to the registration error to adapt to the heterogeneity of error. Combine the window size with the characteristics of different bands to adjust the filter weights to balance noise reduction and preservation of ground feature details. This results in a preprocessed Sentinel-2 raw image with consistent radiometric calibration, improved registration accuracy, and preservation of ground feature details, as well as a Normalized Differential Water Index (NDWI) with high comparability in the same water body under multiple time phases and scenarios.

[0083] Specifically, the radiometric correction steps for the Sentinel-2 raw images are as follows:

[0084] System errors in the blue, green, red, and near-infrared bands of the image were corrected, and the corrected pixel values ​​are as follows:

[0085] Corrected pixel value = Original pixel value × Band gain coefficient (blue 1.05, green 1.02, red 1.01, near-infrared 0.98) - Band offset (obtained by statistically analyzing the black border pixels in the image);

[0086] For each pixel to be processed, the variance σ² is calculated in the 3×3 candidate window, and then compared with the adaptive noise threshold to classify the noise level; if it is classified as strong noise, then hybrid filtering (median + weighted mean) is applied; if it is weak noise, then only weighted mean filtering is required.

[0087] Specifically, the dynamic adjustment logic for the adaptive noise threshold is as follows:

[0088] Using the blue, green, red, and near-infrared bands of Sentinel-2 imagery as units, the pixel grayscale variance distribution of the entire image was statistically analyzed to determine the baseline threshold. The specific steps are as follows:

[0089] Calculate the global baseline threshold for Sentinel-2 images:

[0090] Divide the single-band image into blocks (e.g., 256×256 pixel blocks) and calculate the variance of each block. ;

[0091] variance of all blocks Perform statistical analysis to create a variance dataset and sort it from smallest to largest. Take the 75th percentile of the variance distribution (i.e., 75% of the block variances are less than or equal to this value, and 25% of the block variances are greater than or equal to this value) as the baseline threshold. ;

[0092] Threshold adjustment based on band characteristics: Introducing band weighting coefficients based on band noise sensitivity. Adjust the formula: The band weighting coefficients are obtained through experimental or statistical analysis, taking into account the noise characteristics of each band (blue / green / red / near-infrared) of Sentinel-2.

[0093] Within the same band, further spatial adaptive adjustments are made to address the differences in texture variance across different land cover areas (vegetation, water bodies, buildings) to obtain pixel-level local thresholds. The specific steps are as follows:

[0094] For each pixel's 3×3 neighborhood, calculate the texture complexity of the features within that neighborhood. ;

[0095] Set texture complexity coefficient ( Complex texture areas (For example, vegetation has a high entropy value) (lower the threshold to avoid misclassifying textures as noise); simple texture areas (such as water bodies) have a low entropy value. (Increase the threshold to enhance noise suppression), adjust the formula: Its texture complexity coefficient The value is determined by the complexity of the terrain texture.

[0096] Specifically, the steps for dynamically adjusting the filter window based on registration error are as follows:

[0097] Using RPC parameters or ground control points (GCP), a mapping relationship between image pixel coordinates (row and column numbers) and real geographic coordinates (latitude and longitude / UTM) is established. The geographic coordinates of the image pixel to be corrected are compared with the geographic coordinates of the reference base (reference image / real ground features). The difference is transformed by projection to calculate the registration error (Δx, Δy) of the current pixel. The error level is quantized based on the adaptive error threshold, and the window parameters are dynamically determined according to the error level.

[0098] Specifically, the adaptive error threshold adjustment logic is as follows:

[0099] First, the global registration error distribution of the image to be processed is statistically analyzed to determine the baseline range of the threshold:

[0100] Introducing scene weight coefficients Adjust the baseline threshold and correct the formula according to the application scenario: Among them, the scene weight coefficient The range to be set depends on the application scenario. Pixel-level local threshold Application scenarios include fine-scale ground feature classification / contour extraction, coarse-scale monitoring / coverage estimation, and routine image preprocessing;

[0101] To address the registration error sensitivity of different bands, a band correction coefficient is introduced. Fine-tuning yields the final adaptive error threshold: Among them, band correction coefficient The settings are based on the Sentinel-2 four-band characteristics.

[0102] Specifically, the steps for adjusting the filter weights based on the characteristics of different frequency bands are as follows:

[0103] Based on the window size, a preset basic weight matrix for the bands is defined as follows:

[0104] Based on the weight allocation principle, the sum of all elements in the weight matrix of each band is set to 1 (ensuring that the filtered pixel values ​​are within a reasonable dynamic range). The formula is: , Represents the original weights. This represents the normalized weights, which are ultimately used in the filtering calculation as a weight matrix (the sum of all elements is 1). The sum of the original weights is an unnormalized weight matrix preset based on band characteristics and window levels; and the window levels are defined, and the rationality of the basic weight matrix is ​​preset and verified for different bands based on different window levels.

[0105] Filtering is performed according to the window that divides the noise level; the mean deviation of radiation, PSNR, and mean registration error before and after filtering are calculated. If any one of the following is true: mean deviation of radiation before and after filtering is greater than or equal to 5%, PSNR is less than 3dB, or the mean reduction of registration error is less than 30%, then the standard is not met, and the parameters are adjusted backtracked.

[0106] Specifically, a weighted gradient strategy is used in the boundary regions of different error levels and different bands to avoid image stitching artifacts; the specific settings of the weighted gradient strategy are as follows:

[0107] Determine the width of the boundary region: Transition region width = 1 / 2 of the window size difference, rounded up;

[0108] The basis for determining gradient weights: the weight matrix for each pixel within the transition region. Weight matrix on one side of the boundary and the weight matrix on the other side of the boundary Calculated using a distance-weighted average;

[0109] Quantization gradient: For each pixel in the transition region, calculate the gradient coefficient t, and then obtain the final weight using a linear gradient formula. The specific operation is as follows:

[0110] (1) Calculation of the gradual change coefficient t: ;

[0111] in, This represents the distance from the current pixel to the boundary of the reference region A, within the transition zone. From 0 to ; It is the width of the transition zone;

[0112] (2) The final weight matrix of pixels in the transition region, through... and Linear interpolation yields:

[0113] ;

[0114] in, These are the row and column indices of the weight matrix, calculated for... Renormalize to ensure that the sum of all elements is 1.

[0115] Step 2: Combining the high reflectivity characteristics of clouds with the texture differences of land areas, cloud masking and land masking are performed on the preprocessed Sentinel-2 raw images based on the Normalized Differential Water Index (NDWI) to provide a geometrically consistent effective water body range for subsequent multi-temporal fusion and ICESat-2 photon matching.

[0116] The ICESat-2 photon points are projected onto a coordinate system consistent with the preprocessed Sentinel-2 image. Each photon is associated with the water / non-water attributes and local texture information of the corresponding pixel. The DBSCAN density clustering method is used to perform spatial-height joint clustering of ICESat-2 photon heights. Based on the water connectivity and texture complexity provided by the preprocessed image (e.g., reflected by entropy and local variance), the neighborhood radius and minimum sample number of DBSCAN are adaptively set. This allows isolated noise photons to be preferentially removed in low-texture, low-noise water areas, while enhancing the retention of real underwater topographic photons in complex nearshore areas. The pixel-photon spatial correspondence brought by high registration accuracy is used to reduce pseudo-clustering caused by photons being misprojected onto land or the outside of the shoreline, thereby obtaining underwater topographic photons. Refraction correction and tidal level correction are then performed to obtain a reference water depth that is strictly aligned and standardized pixel by pixel with the preprocessed Sentinel-2 image. The reference water depth is divided into a model training set and a validation set according to the spatial distribution and water depth range.

[0117] Because the ATL03 product contains a large number of background noise photons, and relying solely on the confidence labels provided by the manufacturer is insufficient to accurately separate signal from noise, this invention employs the DBSCAN density clustering method to replace traditional classification methods, in order to extract high-confidence photon points representing real underwater terrain. The specific processing steps are as follows:

[0118] First, all photon points are statistically analyzed in layers according to elevation values ​​to identify water surface locations. Water surface echoes are typically characterized by concentrated density, high energy, and large quantity; therefore, the bottom boundary of the elevation layer with the highest photon count is selected as the boundary between surface and underwater photons. Second, DBSCAN clustering is performed on surface and underwater photons separately to remove noise points and non-topographical photons, obtaining underwater topographically reflected photons. To enhance the local adaptability of the clustering algorithm, the entire photon track is divided into multiple sub-segments at 50m intervals along the ICESat-2 orbit direction, and DBSCAN processing is performed independently on each sub-segment, with adaptive clustering parameters set for different sub-segments. When the number of photons in a neighborhood of a certain segment is less than 3, the minimum neighborhood number is uniformly set to 3 to ensure the basic executability of the clustering algorithm.

[0119] Specifically, in the traditional ICESat-2 data processing workflow, laser pulses are typically simplified to propagate in a straight line through the air, without considering the refraction effect that occurs when the beam crosses the air-water interface. Under clear water conditions, this assumption leads to an underestimation of the underwater photon elevation. Therefore, this invention introduces a refraction correction model based on Snell's law to correct the spatial position of underwater photons in terms of both propagation direction and velocity. The corrected photon position is spatially closer to the actual seabed reflection point, providing more reliable underwater topographic reference data for subsequent depth inversion and accuracy verification.

[0120] Specifically, since ICESat-2 data from different observation dates are used, the extracted water depth information is inevitably affected by tidal variations. Therefore, this invention uses the mean sea surface height (MSSH) provided by the global mean sea surface model DTU21MSS as a reference surface, subtracting the corresponding MSSH value from the ellipsoidal height of each photon point to estimate the relative water depth. Essentially, this process maps water depths from different periods to a single standard sea surface reference system, avoiding systematic biases caused by tidal changes.

[0121] Step 3: The accuracy of remote sensing water depth inversion is closely related to the quality of the input image. However, due to the susceptibility of optical remote sensing images to cloud cover, sea waves, and shoreline disturbances, single-temporal images often contain random noise and discontinuous areas, limiting the stability and generalization ability of the shallow sea SDB model. Therefore, a method that fuses multi-scale features of images with a regional statistical weighting mechanism (SDW-LPT method) is proposed. This method preserves the image structure while suppressing noise, generating multi-temporal fused images. Based on the multi-temporal fused image data, water depth inversion is performed using a multi-band model, a logarithmic ratio model, and a quadratic polynomial ratio model (QPRM).

[0122] Specifically, the steps for fusing multi-temporal images using the SDW-LPT method are as follows:

[0123] Each preprocessed image is decomposed into a three-layer Laplacian pyramid and then filtered using a Gaussian filter. The low-frequency components of each layer are extracted, and the high-frequency detail map of the preprocessed image is obtained by removing the low-frequency components of each layer. The formula is as follows:

[0124] ;

[0125] In the formula, The high-frequency components of the Laplacian at layer n are represented by the high-frequency components of each layer of the pyramid. Since the pre-processed image is a denoised but detail-preserving multi-band image, the high-frequency components of each layer of the pyramid reflect the texture of the seabed and subtle changes in the light of the water, rather than random noise. This is the original input image for layer 𝑖; It is a Gaussian filter. , ;

[0126] To enhance the representation of complex textured regions, a local energy map is constructed in each image to reflect the degree of local variation. This energy map serves as a weighting factor for superimposing high-frequency information from multi-temporal images, including the signal-to-noise ratio obtained from radiometric correction and band-weighted filtering, local texture complexity, and consistency with the reference water depth. The formula is as follows:

[0127] ;

[0128] In the formula, This is the normalized energy weights graph. Indicated by The standard deviation of the 5×5 sliding window centered on the center. The maximum standard deviation for each graph;

[0129] Using the lowest-level low-frequency image of the Laplacian pyramid as the fusion base, high-frequency components are weighted and superimposed layer by layer to obtain the final fused image. The reconstruction formula is as follows:

[0130] ;

[0131] In the formula, For the final merged image; The lowest frequency image extracted from the Lth layer of the pyramid; For the high-frequency components of the nth layer; It is a normalized energy map; the constant 0.5 is a smoothing factor, which aims to improve the fusion contribution in low-energy regions and prevent the loss of details.

[0132] Specifically, the multi-band model is represented as:

[0133] ;

[0134] In the formula, Because of the water depth, and For the parameter to be estimated, For the number of bands, Indicates the first The water reflectance corresponding to each waveband Indicates the reflectivity of deep water areas; These represent the reflectance of the blue, green, and red bands, respectively. is the reflectance of the deep water area; m and n are constant terms, and the logarithmic adjustment factor is adjusted to two factors; this model is widely used in remote sensing water depth inversion by integrating the reflection of seabed information in different wavebands.

[0135] The logarithmic ratio model, mathematically expressed as follows:

[0136] ;

[0137] In the formula, Because of the water depth, and Here, n is the parameter to be estimated, and n is a fixed constant. In the logarithmic ratio model, n = 1000. Indicates the first The water reflectance corresponding to each band; since blue and green light have strong penetrating power in water, this model usually uses the blue and green band reflectance for calculation.

[0138] The quadratic polynomial ratio model is used for water depth estimation. First, the water depth index W is calculated using the following formula:

[0139] ;

[0140] In the formula, The depth index is used in the quadratic polynomial ratio model, where n is taken as 1500. , The water reflectance is defined for the blue and green bands, respectively; subsequently, a quadratic polynomial ratio model is established:

[0141] ;

[0142] In the formula, Because of the water depth, Let W be the median value, and a, b, and c be the coefficients to be estimated, which can be obtained by fitting using the least squares method.

[0143] Step 4: Using ICESat-2 validation data as a benchmark, analyze the inversion results of the multi-band model, logarithmic ratio model, and quadratic polynomial ratio model (QPRM), calculate statistical indicators such as root mean square error (RMSE), systematic bias (Bias), and coefficient of determination (R²), and evaluate the performance of the fusion method.

[0144] To evaluate the water depth inversion results driven by multi-temporal fusion images, this invention employs RMSE, Bias, and R... 2 The formula for measuring the consistency and reliability between the retrieved water depth value and the measured water depth data is as follows:

[0145] ;

[0146] ;

[0147] ;

[0148] In the formula, To invert water depth, To represent the actual water depth, To verify the number of water depth points. and These are the average values ​​of the inverted value and the measured value, respectively.

[0149] 1. Study Area

[0150] This embodiment selects Atoll as the experimental area. Atoll A consists of multiple reef platforms, lagoons, and offshore cliffs of varying sizes, exhibiting typical coral reef geomorphological characteristics such as dramatic topographic relief, diverse seabed types, and significant water depth gradients. The water in this area is clear, with exposed seabed and good optical permeability. However, Atoll A is located in a tropical marine climate zone, experiencing frequent fog and significant fluctuations in aerosol concentrations year-round. This can easily lead to problems such as occlusion, uneven brightness, and edge perturbation in single-temporal images, limiting the accuracy and stability of the traditional single-image SDB method.

[0151] 2. Data and Preprocessing

[0152] 2.1 Sentinel-2 Image Data

[0153] The optical remote sensing images used in this study are from the Sentinel-2 satellite system. The system's multispectral instrument (MSI) operates in a 786km sun-synchronous orbit, providing a wide swath imaging capability of 290km. Through the coordinated operation of the Sentinel-2A and Sentinel-2B twin satellites, high-frequency coverage of global land and coastal areas can be achieved every 5 days, offering high temporal resolution and observational stability. Sentinel-2 provides multispectral observation data in 13 bands, with the blue, green, red, and near-infrared bands having a spatial resolution of 10m. All image data is available through the Copernicus Open Data Platform (https: / / browser.dataspace.copernicus.eu / ).

[0154] This embodiment uses Sentinel-2 Level-2A images as the data source. This level of imagery has undergone radiometric and atmospheric correction and can be directly used for water reflectance extraction and water depth inversion modeling. To eliminate land interference, an improved Normalized Difference Water Index (NDWI) was used to mask non-water areas. The imagery used spanned from January 1, 2019 to January 1, 2022, with a cloud cover of no more than 10% as the selection criterion. A total of 36 Sentinel-2 images meeting the criteria were ultimately obtained.

[0155] 2.2. ICESat-2 laser altimetry data

[0156] (1) Data Overview

[0157] This embodiment uses ATL03 photon point cloud data provided by the ICESat-2 satellite as the training sample and validation benchmark for the water depth inversion model. The Advanced Topographic Laser Altimeter System (ATLAS) onboard ICESat-2 is an important space observation platform for current shallow sea bathymetry research. This system operates at a high pulse repetition frequency of 10 kHz, splitting a single green laser pulse into six independent beams arranged in a dual-beam structure. The spacing within each beam group is approximately 90 m, and the spacing between groups is approximately 3 km. The ATL03 product is the foundational photon dataset for ICESat-2, providing information such as the longitude, latitude, corresponding ellipsoidal elevation, reflection intensity, and classification label for each detected photon, possessing rich three-dimensional structure and spatiotemporal attributes. The ATL03 data used in this embodiment covers the Atoll region at location A and can be downloaded from the NASA EarthData open platform (https: / / search.earthdata.nasa.gov / search).

[0158] (2) ICESat-2 data processing

[0159] Table 1 lists the transit time and beam information of the ICESat-2 ATL03 data covering the Atoll area used in this embodiment.

[0160] Table 1 Overview of ICESat-2 data for the study area

[0161]

[0162] Approximately 90,000 reliable underwater topographic photon data points were extracted from the study area, distributed at depths ranging from 0 to 20 meters. The ICESat-2 data still require preprocessing, including background noise removal, refraction correction of underwater photon paths, tidal level correction, and division of the training and validation sets.

[0163] To ensure the water depth inversion model has good physical representativeness during the modeling phase, this embodiment divides the preprocessed ICESat-2 underwater photon data into training and validation sets. Considering the non-uniformity of water depth distribution in shallow sea areas, random partitioning might lead to insufficient representativeness of the dataset at specific depth ranges, ultimately affecting the model's generalization ability. Therefore, this embodiment adopts a stratified sampling strategy based on depth distribution. First, all valid photon points are sorted in ascending order according to water depth values. Then, samples are drawn at 10%, 20%, and 30% quantile intervals to construct the validation set, ensuring representative sample distribution across shallow, medium, and deep water ranges. The remaining 90% of the data is used for model training. This partitioning method balances the breadth of sample distribution with modeling depth sensitivity, not only improving the model's adaptability to complex water depth gradients but also providing reliable data for subsequent accuracy evaluation.

[0164] 3. Experimental Design

[0165] 3.1. Multi-temporal image fusion experiment

[0166] To systematically evaluate the impact of multi-temporal fusion methods on the accuracy of water depth inversion, this embodiment designed a comparative experiment of three fusion strategies: median synthesis, maximum outlier removal, and SDW-LPT. The median synthesis method represents a robust statistical approach, focusing on outlier suppression; the maximum outlier removal method is based on spectral dispersion screening, emphasizing the removal of systematic noise. These two methods differ significantly in their mechanism and processing logic, providing an effective benchmark for the SDW-LPT method.

[0167] Median composite is a classic temporal fusion method. Its core idea is to calculate the median reflectance of each band in multi-temporal images at the pixel level to construct a fused image. This method requires no complex modeling and has advantages such as simple implementation, high computational efficiency, and strong noise suppression capabilities. However, when there is cloud cover or drastic edge changes between temporal images, this method is prone to problems such as blurred fusion results, leading to loss of image details.

[0168] The maximum outlier removal method, based on the principle of spectral consistency, first calculates the spectral deviation of each pixel in the full-time sequence. Then, it progressively removes the most inconsistent pixels from high to low outlier until a set threshold is met or the number of retained pixels reaches a lower limit. Finally, the filtered pixel set is averaged to generate the final fused image. Compared to the median fusion method, this method is more capable of identifying extreme outliers, significantly improving the spectral consistency and physical plausibility of the fused image, and mitigating systematic errors caused by atmospheric disturbances, sea surface changes, or observation conditions to some extent. However, its fusion result is sensitive to the threshold and sample size, and may lead to loss of detail when there are significant differences in the quality of time-series images or uneven acquisition intervals.

[0169] The experimental area was selected as Atoll A. The image data came from 36 Sentinel-2 Level-2A images with cloud cover not exceeding 10% from January 1, 2019 to January 1, 2022. To ensure fair comparison, a 3×3 mean filter was uniformly applied to the Sentinel-2 images before processing by all three fusion methods to reduce random noise and suppress geometric registration error. The parameter settings in the maximum outlier removal method were based on the recommended values ​​in existing literature (CHU S, CHENG L, RUAN X, et al. Technical framework for shallow-water bathymetry with high reliability and Nomissing data based on time-series sentinel-2 images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(11): 8745-8763.).

[0170] To further explore the impact of the number of fused images on the accuracy of water depth inversion, this embodiment sets a minimum of 3 fused images to meet the basic requirements of the median composite method and the maximum outlier removal method for time series length. Based on this, multiple fusion groups are constructed by adding one image at a time. For each image quantity, five different temporal combinations are randomly selected for fusion experiments, and the results are input into each fusion strategy for water depth inversion. The mean and standard deviation of accuracy indicators are obtained through multiple experiments to evaluate the performance stability and inversion reliability of each fusion method under different temporal coverage, providing a quantitative basis for reasonably determining the fusion time series length in practical applications.

[0171] 3.2. Water depth inversion and accuracy assessment

[0172] The fused images generated by combining different fusion strategies and the number of fused images in section 3.1 were then input into the multi-band model, the logarithmic ratio model, and the QPRM for depth inversion. Based on the same image input, the applicability of different inversion models was evaluated, which helps to analyze the interaction between the fusion strategy and the inversion model, and further verifies the effectiveness of multi-temporal fusion in improving inversion accuracy.

[0173] To objectively evaluate the accuracy performance of different fusion strategies and inversion models, this embodiment uses photon depth data from the ICESat-2 ATL03 product as a reference. The validation data underwent DBSCAN denoising, refraction correction based on Snell's law, and tidal level correction using the DTU21MSS model in the previous processing stage, ultimately obtaining standardized measured depth values. The inversion results were matched pixel-level with measured photon points using spatial registration, and accuracy evaluation metrics were calculated for each fusion and model combination. Accuracy evaluation adopted the aforementioned evaluation system, including RMSE, Bias, and R... 2 By comparing the performance of different fusion strategies and model combinations under various indicators, the applicability and effectiveness of the multi-temporal fusion method in water depth inversion are systematically analyzed.

[0174] 4. Results and Analysis

[0175] 4.1. Comparison of water depth inversion accuracy using different fusion strategies

[0176] To systematically evaluate the impact of multi-temporal image fusion strategies on the accuracy of shallow water depth inversion, this embodiment uses 36 Sentinel-2 images of the Atoll area at location A as a basis, and constructs fused images using SDW-LPT, median composite method, and maximum outlier removal method, respectively. Based on this, nine fusion-modeling combination schemes are formed by combining three typical SDB modeling methods: multi-band model, logarithmic ratio model, and quadratic polynomial ratio model, and comprehensive experimental comparisons are conducted. The experiment uses the coefficient of determination R0. 2The root mean square error (RMSE) and systematic bias (Bias) were used as accuracy evaluation indicators, and the analysis was conducted from two dimensions: overall fitting accuracy and segmented water depth performance. The results are as follows: Figures 1-9 As shown.

[0177] Figures 1-9 The scatter plot relationship between the fused image-driven water depth inversion results and the ICESat-2 measured water depth is shown. Figures 1-9 The horizontal axis represents the measured water depth by photons, and the vertical axis represents the water depth retrieved by remote sensing. The solid blue line represents the straight line obtained by fitting the data points, and the dashed red line represents the ideal fitted line y=x, with accuracy indicators labeled. The horizontal subplots represent the performance of the same water depth retrieval model under different image fusion strategies, used to quantitatively compare the differences in retrieval accuracy of fusion methods; the vertical subplots reflect the regression performance of different modeling methods under the same fusion strategy, used to evaluate the robustness of fused images under each retrieval model.

[0178] Figure 1 , Figure 4 , Figure 7 This is the inversion result of SDW-LPT; Figure 2 , Figure 5 , Figure 8 This is the inversion result of the median synthesis method; Figure 3 , Figure 6 , Figure 9 This is the inversion result of the maximum outlier removal method. When using a multi-band model for modeling, such as... Figure 1 The SDW-LPT fused image shown achieved the best inversion performance, R 2 The Pt / R² value reached 0.93, and the RMSE was 0.45m, demonstrating the best fit. In contrast, Figure 2 Median synthesis method and Figure 3 The R² values ​​of the maximum outlier removal method were 0.91 and 0.92, and the RMSE values ​​were 0.49m and 0.50m, respectively. Although they also had high fitting accuracy, they were weaker than the SDW-LPT method in terms of both accuracy and stability.

[0179] Figures 1-3 The inversion results of the multi-band model under three types of fused imagery; Figures 4-6 The results are the inversion results of the log-ratio model under three types of fused images; Figures 7-9 This presents the inversion results of the quadratic polynomial ratio model under three types of fused imagery. Under the same fusion strategy, the multi-band model generally exhibits the best fitting effect, with its R² and RMSE indices being superior to the other two modeling methods. Taking SDW-LPT fused imagery as an example... Figure 1 The multi-band model achieved an R² of 0.93 and an RMSE of 0.45m, while Figure 4 The log-ratio model in Figure 7The R² values ​​of the quadratic polynomial ratio model in the two methods are 0.84 and 0.90, respectively, and the RMSE values ​​are 0.68m and 0.52m, respectively, both of which are weaker than the multi-band model. Similar trends exist in the median composite method and the maximum outlier removal method, further verifying that the multi-band model has better inversion performance and robustness under various fused image driving conditions.

[0180] To further reveal the adaptability of different image fusion strategies in various water depth ranges, the RMSE and Bias indices of each fusion method in four typical depth ranges of 0~5m, 5~10m, 10~15m and 15~20m were statistically analyzed. The results are shown in Table 2.

[0181] Table 2. Comparison of accuracy of different fusion strategies and water depth inversion models at various water depths.

[0182]

[0183] Table 2 shows that the SDW-LPT fusion strategy exhibits good inversion performance across all water depth ranges, with a more pronounced effect in the shallow to medium water depth region. In the 0–5 m water depth range, the RMSE of SDW-LPT combined with the multi-band model and the quadratic polynomial ratio model are 0.38 m and 0.35 m, respectively, significantly better than other fusion methods. Even in the logarithmic ratio model with relatively large errors, the RMSE of SDW-LPT is still 0.58 m, still superior to other methods. In the 5–10 m water depth range, SDW-LPT and the median synthesis method perform similarly under the multi-band model, with an RMSE of 0.52 m for both. However, in the logarithmic ratio model and the quadratic polynomial ratio model, the RMSE of SDW-LPT is superior to other fusion strategies, demonstrating stronger adaptability and stability in shallow to medium water depths. Entering the 10-15m range, the reduced sample size led to an increase in error. The median synthesis method showed slightly better RMSE in some models, but SDW-LPT performed more robustly in terms of systematic bias control, with the lowest bias of only -0.74m. In the 15-20m water depth range, the sample size was severely insufficient, and the errors of all models increased significantly. However, the RMSE of the combination of SDW-LPT and the quadratic polynomial ratio model remained at 4.80m, still the lowest, demonstrating good noise resistance and robustness.

[0184] Figures 1-9The results in Table 2 show that the SDW-LPT fusion strategy outperforms traditional methods in overall accuracy, adaptability to different inversion models, and performance across various water depths. Particularly in the shallow and mid-shallow water regions (0–10 m), the multi-band model and quadratic polynomial ratio model driven by SDW-LPT fusion imagery achieved the best RMSE index. Even in deep water where training samples are scarce, the strategy maintains strong stability and noise resistance. In summary, the SDW-LPT fusion method demonstrates excellent inversion capabilities across various models due to its ability to enhance detailed features and suppress outliers.

[0185] Furthermore, under different fusion strategies, the error of the same inversion model increases overall with water depth, and the negative bias is more significant in deep water. There are two main influencing factors: First, with increasing depth, the two-way attenuation of the water column increases, causing the contribution of the bottom to surface reflection to decrease exponentially, resulting in a decrease in effective SNR; second, with increasing depth, the number of usable depth points for ICESat-2 decreases significantly, and single-point uncertainty increases (affected by receptor scattering, refraction correction residuals, etc.), leading to a scarcity of training samples and exacerbating heteroscedasticity. The combination of these two factors causes RMSE to increase with depth, exhibiting a stronger systematic underestimation in deep water.

[0186] 4.2 Impact of Image Quantity on Water Depth Inversion Accuracy

[0187] To systematically evaluate the impact of the number of multi-temporal images on the accuracy of water depth inversion, the number of input images was used as a control variable in the experiment while keeping the fusion strategy and inversion model combination unchanged. Specifically, 3 to 36 images were selected as sample inputs, increasing progressively. For each set of image numbers, 5 sets of combined samples were randomly generated, and the fusion and inversion process was executed. The corresponding root mean square error was then calculated. Figure 2 The study demonstrates the trend of root mean square error (RMSE) as a function of the number of images and its error fluctuation range under three modeling frameworks: multi-band model, log-ratio model, and quadratic polynomial ratio model.

[0188] Figure 2The results show that different fusion strategies exhibit highly consistent trends during the depth inversion process. First, as the number of fused images gradually increases, the RMSE of the three depth inversion models generally decreases, and the error fluctuations gradually converge, fully validating the significant advantages of multi-temporal fusion in improving inversion accuracy and result stability. Second, under all model conditions, the SDW-LPT fusion strategy consistently achieves the lowest RMSE value and the smallest error variance, demonstrating superior stability and broad adaptability. Further observation reveals that in the low-sample stage with 3 to 10 images, SDW-LPT maintains stable error control, while the median synthesis method and the maximum outlier removal method exhibit significant error fluctuations in this stage, indicating their sensitivity to sample size and limited generalization ability. Overall, SDW-LPT demonstrates particularly strong robustness under small sample conditions.

[0189] Further analysis from the perspective of different modeling methods reveals that the SDW-LPT fusion strategy exhibits good convergence characteristics under both the multi-band model and the quadratic ratio model. When the number of fused images increases to approximately 20, the RMSE curves of both models tend to stabilize, and the error standard deviation converges significantly, indicating that under sufficient sample size, the model performance can quickly stabilize and reach a high level of accuracy. In contrast, in the log-ratio model, although the RMSE also gradually decreases with the increase in the number of images, the overall convergence speed is slower, and the error fluctuation is larger, showing that this model is more sensitive to the quality of the fused images and input noise.

[0190] In summary, the number of images significantly affects the accuracy of water depth inversion. As the number of images increases, the accuracy of all models generally improves, and error fluctuations tend to converge. Among them, the SDW-LPT fusion strategy consistently performs best in terms of accuracy and stability, especially showing greater robustness when there are insufficient samples or significant noise interference. In contrast, the median synthesis method and the maximum outlier removal method have larger errors and poorer stability with a low number of images.

[0191] This embodiment addresses the accuracy and stability issues in shallow water depth remote sensing inversion by proposing a Laplacian pyramid fusion strategy based on a local standard deviation weighting mechanism. Using 36 Sentinel-2 images of Atoll A as experimental samples, experiments were conducted under three frameworks: multi-band model, logarithmic ratio model, and quadratic polynomial ratio model. Results show that the SDW-LPT method has significant advantages in improving the quality of fused images, enhancing model robustness, and reducing water depth inversion errors. Compared with the median synthesis method and the maximum outlier removal method, SDW-LPT achieves R² and RMSE of 0.93 and 0.45m respectively in the multi-band model, and also maintains leading performance in the logarithmic ratio and quadratic polynomial ratio models. Furthermore, the SDW-LPT method exhibits higher accuracy in the 0–10m range, with the highest accuracy at 0.35m in the 0–5m range. Even in the 15–20m depth range, where samples are scarce and noise interference is significant, the SDW-LPT method still demonstrates good robustness and adaptability. Further experiments with incremental image fusion showed that the model converged when the number of images reached approximately 20, maintaining stable performance even with low sample counts. Therefore, this study validated the comprehensive performance of the SDW-LPT fusion strategy across multiple dimensions, including model accuracy improvement, error control, and sample adaptation. However, it should be noted that factors such as water optical extinction, the scarcity of ICESat-2 depth points in deep water sections, and increased uncertainty at single points can all cause heteroscedasticity in errors with depth and introduce systematic biases, constituting sources of uncertainty in topographic inversion accuracy.

[0192] This invention, after radiometric correction, adaptive window adjustment based on registration error, and optimization of band-specific filter weights, obtains preprocessed Sentinel-2 multispectral images with consistent radiometric calibration, improved registration accuracy, and preserved ground feature details. Based on this, a depth-oriented inversion method is proposed, including:

[0193] (1) NDWI masking and effective water extraction based on high-quality preprocessed images

[0194] Normalized Differential Water Index (NDWI) was calculated using the preprocessed Sentinel-2 blue, green, red, and near-infrared reflectance bands. Because radiometric correction eliminates system gain and offset errors, and the band-specific filtering weights suppress noise while preserving details such as shorelines and shoals, the NDWI numerical distribution is more stable, resulting in higher comparability of the same water body across multiple time phases and scenes. Based on this, combining the high reflectance characteristics of clouds and the texture differences in land areas, a unified or adaptive threshold strategy was used to perform cloud and land masking on the NDWI images, extracting continuous and clearly defined effective water body regions. Simultaneously, the registration error field was used to perform sub-pixel-level correction on the mask boundaries, avoiding water-land boundary jitter caused by multi-temporal registration errors, thus providing a geometrically consistent water body range for subsequent multi-temporal fusion and ICESat-2 photon matching.

[0195] (2) ICESat-2 photonic denoising and normalized reference depth construction based on water mask and texture prior

[0196] ICESat-2 photon points are projected onto a coordinate system consistent with the preprocessed Sentinel-2 image. Utilizing the high-precision registration relationship obtained during preprocessing, each photon is associated with the water / non-water attributes and local texture information of its corresponding pixel. While retaining only valid water areas, the DBSCAN density clustering method is used for spatial-height joint clustering of photon heights. Based on the water connectivity and texture complexity provided by the preprocessed image, the neighborhood radius and minimum sample number of DBSCAN are adaptively set. This prioritizes the removal of isolated noise photons in low-texture, low-noise water areas, while enhancing the retention of photons representing real underwater topography in complex near-shore areas. The pixel-photon spatial correspondence resulting from high registration accuracy reduces pseudo-clustering caused by photons being misprojected onto land or outside the shoreline.

[0197] Refraction and tidal level corrections were applied to the underwater topographic photons obtained from clustering. A standardized reference depth dataset, strictly aligned pixel-by-pixel with the preprocessed Sentinel-2 image, was constructed and divided into a model training set and a validation set based on spatial distribution and depth range. Due to the improved radiometric consistency and registration accuracy in the early stages, the dispersion of multi-temporal reflectance scatter points corresponding to the same depth was reduced, and the physical consistency of the training samples was significantly enhanced.

[0198] (3) Enhancement of multi-temporal image fusion and water depth-sensitive information using SDW-LPT

[0199] For preprocessing Sentinel-2 images of multiple time phases in the same area, the SDW-LPT multi-time phase fusion method is used within the water mask area.

[0200] In the Laplacian pyramid decomposition stage, the denoised but detail-preserving multi-band imagery directly affects the preprocessed output, enabling the high-frequency components of each pyramid layer to more accurately reflect the shallow seabed texture and subtle changes in water illumination, rather than random noise. In the multi-temporal weighting stage, the signal-to-noise ratio (SNR) obtained from radiometric correction and band-weighted filtering, local texture complexity, and consistency with the reference depth are all incorporated into the SDW-LPT weighting design. This ensures that depth-sensitive phases with smaller registration errors and lower noise receive higher weights during fusion. In areas with larger errors, the registration error level and window parameters from preprocessing are used to constrain the spatial interpolation scale during fusion, preventing the amplification of pseudo-textures introduced by incorrect registration. The final result is a multi-temporal fused image that outperforms single-temporal images in spectral consistency, spatial registration, and detail fidelity, providing high-quality spectral input for subsequent depth inversion.

[0201] (4) Multi-model collaborative water depth inversion based on fused images

[0202] Using the multi-temporal water reflectance fused by SDW-LPT as input and the standardized reference water depth, after refraction and tidal correction and strictly corresponding to the pixel, as output, a multi-model collaborative water depth inversion system is constructed: The multi-band model utilizes the absolute reflectance and their combination characteristics of the blue, green, red, and near-infrared bands after radiometric correction, and improves the response sensitivity to different water depth ranges by leveraging the bottom sediment and water body detail textures preserved in preprocessing; The logarithmic ratio model constructs a logarithmic ratio water depth index based on high-quality blue-green bands, and the robustness of this index under multi-temporal and multi-water body conditions is ensured by noise suppression and radiometric consistency in the early stage; The quadratic polynomial ratio model QPRM further introduces quadratic and cross terms on the basis of logarithmic ratio, fully exploring the nonlinear relationship between multiple bands. Its parameter fitting relies on a "cleaner" spectrum-water depth correspondence in the sample, which significantly improves the model's fitting ability in medium-deep water and complex bottom sediment areas.

[0203] The three models work collaboratively under the same fused image and unified reference depth system: the multi-band model provides an overall trend characterization, the log-ratio model enhances sensitivity to transparent water bodies in shallow areas, and QPRM corrects nonlinear residuals based on this. By weighted fusion or partition selection of the inversion results from the three models, high-precision inversion of shallow water depth is achieved.

[0204] Therefore, this invention achieves high-precision and stable acquisition of shallow water depth information through an integrated technology chain from "active adaptation + precise quality improvement" preprocessing to "multi-source collaboration + multi-model joint inversion".

[0205] Although the present invention has been described in detail above with general descriptions and specific embodiments, modifications or improvements can be made to it, which will be obvious to those skilled in the art. Therefore, all such modifications or improvements made without departing from the spirit of the present invention fall within the scope of protection claimed by the present invention.

Claims

1. A water depth inversion method based on SDW-LPT multi-temporal fusion remote sensing image, characterized in that: The method comprises the following steps: The acquired Sentinel-2 original image is subjected to radiation correction to eliminate system errors, and then a filter window is adjusted according to a registration error, a filter weight is adjusted according to characteristics of different bands in combination with a window size, and a preprocessed Sentinel-2 original image is obtained and a normalized difference water index NDWI is calculated, so as to balance noise reduction and detail preservation of ground objects; The Sentinel-2 original image acquisition step is as follows: four key band data of blue, green, red and near-infrared are extracted from the Sentinel-2 image, and image RPC parameters or ground control points GCP are extracted; The radiation correction step of the Sentinel-2 original image is as follows: System errors of the image blue, green, red and near-infrared bands are corrected, and the corrected pixel value is as follows: Corrected pixel value = original pixel value * band gain coefficient - band offset; The variance σ 2 of a 3*3 candidate window of each pixel to be processed is calculated, and then compared with an adaptive noise threshold to divide the noise level; if it is strong noise, mixed filtering is performed, and the mixed filtering includes median and weighted mean filtering; if it is weak noise, only weighted mean filtering is performed; The dynamic adjustment logic of the adaptive noise threshold is as follows: The Sentinel-2 image is taken as a unit of four bands of blue, green, red and near infrared, the pixel gray scale variance distribution of the whole image is counted respectively to determine the reference threshold The specific operation is as follows: The global reference threshold of the Sentinel-2 image is calculated: The single-band image is divided into blocks, and the variance of each block is calculated ; variance of all blocks do statistical analysis, form a variance dataset and sort from small to large, take the 75 percentile of variance distribution as the benchmark threshold ; Band weight coefficient is introduced based on band noise sensitivity , adjustment formula: adjusted threshold ; band weight coefficient is combined with the noise characteristics of each band of Sentinel-2, obtained through experimental or statistical analysis; In the same band, for different texture variance of ground object area, further spatial adaptive adjustment, get pixel level local threshold The specific operation is as follows: For each pixel, a texture complexity coefficient is set for its 3x3 neighborhood , the adjustment formula is: ; determined by the degree of texture fluctuation of the pixel's local neighborhood The filter window dynamic adjustment step based on the registration error is as follows: The mapping relationship between the image pixel coordinates and the real geographical coordinates is established by using the RPC parameters or the ground control points GCP, the geographical coordinates of the pixel to be corrected are compared with the geographical coordinates of the reference benchmark, the difference is converted by projection to calculate the registration error of the current pixel; the adaptive error threshold is quantized into an error level, and the window parameters are dynamically determined according to the error level; The adaptive error threshold adjustment logic is as follows: First, the global registration error distribution of the image to be processed is counted to determine the reference interval of the threshold, and the steps are as follows: Introducing scene weight coefficient , adjusting the reference threshold according to the application scene, and revising the formula: ; wherein the scene weight coefficient The setting range depends on the application scene, is the pixel-level local threshold The application scene includes fine feature classification / contour extraction, coarse scale monitoring / coverage estimation, and regular image preprocessing; A band correction coefficient is introduced for the sensitivity of the registration error of different bands , and a final adaptive error threshold is obtained by fine tuning: ; wherein the setting of the band correction coefficient is based on the characteristics of the four bands of Sentinel-2; The filter weight adjustment step according to the characteristics of different bands is as follows: The window size is combined, and a preset band basic weight matrix is set, which is as follows: Based on the weight distribution principle, the sum of all elements of the weight matrix of each band is equal to 1, formula: , represents the original weight, represents the sum of the original weight, represents the normalized weight, the final weight matrix used for filtering calculation is the un-normalized weight matrix preset according to the band characteristics and window level; and the level of the window is defined, based on different window levels, the base weight matrix is preset and verified for reasonableness by bands; The filter is executed according to the window of the divided noise level; the radiation mean deviation before and after filtering, PSNR and registration error mean are calculated, and if any one of the following conditions is met, it is not up to standard, and the parameters are adjusted back: In the boundary area of different error levels and different bands, a weight gradual change strategy is adopted to avoid image splicing marks; the specific setting of the weight gradual change strategy is as follows: The width of the boundary area is determined: transition zone width = 1 / 2 of the window size difference, rounded up; The calculation basis for determining the gradient weight: the weight matrix of each pixel in the transition zone The weight matrix on one side of the boundary The weight matrix on the other side of the boundary Weighted calculation according to distance proportion Quantitative gradient: for each pixel in the transition zone, the gradient coefficient t is calculated, and the maximum weight is obtained through a linear gradient formula, and the specific operation is as follows: (1) Gradation coefficient t calculation: ; wherein, represents the distance of the current pixel to the border of the reference A region within the transition zone from 0 to ; is the width of the transition zone; (2) The most right weight matrix of the pixels in the transition zone is obtained by and linear interpolation. ； wherein, are row and column indices of the weight matrix, computed after renormalized to ensure all elements sum to 1; The preprocessed Sentinel-2 original image is subjected to cloud mask and land mask processing based on the normalized difference water index NDWI, so as to provide an effective water body range with geometric consistency for subsequent multi-temporal fusion and ICESat-2 photon matching. ICESat-2 photon points are projected into the coordinate system consistent with the preprocessed Sentinel-2 image, and the DBSCAN density clustering method is used for spatial-height joint clustering of photon height. The neighborhood radius and minimum sample number of DBSCAN are set adaptively, and isolated noise photons are removed to enhance the retention capability of real underwater terrain photons. Underwater terrain photons are extracted, and refraction correction and tide correction are performed to obtain the standardized reference water depth aligned with the preprocessed Sentinel-2 image pixel by pixel. According to the spatial distribution and water depth range division model, the training set and the verification set are divided; SDW-LPT method is used to fuse multi-temporal images to generate multi-temporal fusion images. Based on the multi-temporal fusion image data, multi-band model, logarithmic ratio model and quadratic polynomial ratio model QPRM are used for water depth inversion.

2. The water depth retrieval method based on SDW-LPT multi-temporal fusion remote sensing image according to claim 1, characterized in that: The specific steps of SDW-LPT method for fusing multi-temporal images are as follows: The pre-processed image is decomposed into three layers of Laplacian pyramid, and a Gaussian filter is used The low-frequency component of each layer is extracted, and the high-frequency detail image is obtained by removing the low-frequency component of each layer, and the formula is as follows: ； In the formula, represents the Laplacian high-frequency component of the i-th layer; is the original input image of the i-th layer; is a Gaussian filter, , ; A local energy map is constructed in each image to reflect the local change degree of the region and serve as a weighting factor to realize the superposition of high-frequency information of multi-temporal images. The formula is as follows: ； wherein is the normalized energy weight map, denotes the standard deviation of a 5x5 sliding window centered at is the maximum standard deviation for each image.​ The lowest frequency image of the Laplacian pyramid is used as the fusion base, and the high-frequency components are weighted and superimposed layer by layer to obtain the fusion image. The reconstruction formula is as follows: ； wherein, is the final fused image; is the lowest frequency image extracted from the pyramid Lth layer; is the high frequency component of the ith layer; is the normalized energy map; the constant 0.5 is a smoothing factor.

3. The water depth retrieval method based on SDW-LPT multi-temporal fusion remote sensing image according to claim 1, characterized in that: The multi-band model is represented as: ； wherein, is the water depth, and is the parameter to be estimated, is the number of wave bands, represents the reflectivity of the water body corresponding to the th wave band, represents the reflectivity of the deep water area under the th wave band; are the reflectivities of the blue, green and red wave bands, respectively. m, n are constant terms, and the logarithmic adjustment factor is adjusted to two factors; The logarithmic ratio model is as follows: ； In the formula, Because of the water depth, and Here, n is a fixed constant and is the parameter to be estimated. Indicates the first Water reflectance corresponding to each waveband; The quadratic polynomial ratio model is used for water depth estimation. First, the water depth index W is calculated as follows: ； In the formula, is the water depth index, , The water reflectance of two bands selected from the four key bands of blue, green, red and near-infrared; then, a quadratic polynomial ratio model is established: ； wherein is the water depth, is the median value of W, a, b, c are coefficients to be estimated, obtained by least square fitting.

Images