Water volume calculation method based on spaceborne lidar data and geometric space constraint
By using a method based on spaceborne lidar data and geometric spatial constraints, the problem of poor morphological adaptability in lake underwater topography estimation under sparse satellite trajectory data was solved, achieving high-precision lake water volume calculation, which is suitable for water volume monitoring in complex environments.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHINA UNIV OF GEOSCIENCES (WUHAN)
- Filing Date
- 2026-04-03
- Publication Date
- 2026-07-10
AI Technical Summary
Existing technologies, when using sparse satellite trajectory data to extrapolate the underwater topography of the entire lake, have poor model morphology adaptability, are prone to generating unreasonable extrapolation values in areas without data, and are difficult to efficiently monitor lake water volume in complex environments.
A method based on spaceborne lidar data and geometric spatial constraints is adopted. The lake boundary is extracted through multispectral remote sensing images, the photon point set is separated by kernel density estimation, the water depth is corrected by Snell's law, a geometric scaling factor and a polynomial fitting model are constructed, and the water volume is calculated by volume integral method. Geometric scaling constraints and an adaptive polynomial fitting model are introduced.
It achieves high-precision lake water volume estimation under sparse data conditions, prevents unreasonable extrapolation, adapts to various lake morphologies, and is suitable for efficient water volume monitoring in complex environments.
Smart Images
Figure CN122362371A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the fields of remote sensing geographic information processing and hydrological monitoring technology, and in particular to a method for calculating water volume based on spaceborne lidar data and geometric spatial constraints. Background Technology
[0002] Lake bathymetry and water volume are key parameters for global water cycle research, water resource management, and regional climate change monitoring. Traditional lake bathymetry mainly relies on sonar ship surveys or manual vertical surveys. These methods are costly, inefficient, and difficult to implement on a large scale in data-scarce areas with complex terrain and harsh environments, such as the Qinghai-Tibet Plateau.
[0003] With the development of satellite remote sensing technology, it has become possible to use spaceborne photon counting lidar such as ICESat-2 for underwater topographic surveying. ICESat-2 can acquire high-precision underwater photon point clouds, but its ground trajectory is distributed in a discrete linear pattern and is limited by orbital intervals, often only covering local areas of lakes, such as only the edges or parts of the bays.
[0004] Existing techniques for extrapolating the topography of a whole lake using sparse trajectories often employ parabolic fitting or conventional interpolation methods. However, these methods suffer from significant drawbacks, such as poor morphological adaptability. For instance, simple linear interpolation or logarithmic fitting struggles to simultaneously adapt to various complex lake sedimentary morphologies, including V-shaped valley lakes and U-shaped basin lakes. Therefore, there is an urgent need for modeling methods that can ensure the rationality and stability of topographic extrapolation by introducing physical geometric constraints under conditions of data sparsity. Summary of the Invention
[0005] This invention proposes a water volume calculation method based on spaceborne lidar data and geometric spatial constraints, aiming to solve the technical problems in the existing technology of poor model morphology adaptability and unreasonable extrapolation values in areas without data when using sparse satellite trajectory data to extrapolate the underwater topography of the whole lake.
[0006] This invention provides a method for calculating water volume based on spaceborne lidar data and geometric spatial constraints, including:
[0007] S1. Acquire multispectral optical remote sensing images covering the target lake, calculate the normalized differential water index, extract the lake's land-water boundary through threshold segmentation, and generate a closed lake vector polygon.
[0008] S2, acquire ICESat-2 ATL03 photon point cloud data covering the target lake, separate background noise in the vertical elevation direction using kernel density estimation, and extract the photon point set on the lake surface and the photon point set on the lake bottom;
[0009] S3, based on Snell's law, performs two-way refraction correction on the extracted lake bottom photon point set to obtain trajectory data containing geographic coordinates and measured water depth;
[0010] S4. Perform a full-range scan of the geometric centroid of the lake vector polygon to obtain the maximum geometric radius and corresponding water depth of the lake, as well as the maximum coverage radius and corresponding water depth of the measured trajectory; construct a geometric scale factor, and combine it with the measured maximum water depth to calculate the physical upper limit of the maximum depth of the entire lake.
[0011] S5. Calculate the normalized distance from each measured water depth point to the geometric centroid, establish a polynomial fitting model with the normalized distance as the independent variable and the water depth as the dependent variable, optimize the parameters using the least squares method, and use the physical upper limit value as the inequality constraint condition of the solution space.
[0012] S6. The interior of the lake vector polygon is rasterized, and a normalized distance field of the whole lake is generated based on Euclidean distance transformation. The optimal parameters obtained in step S5 are substituted into the polynomial fitting model, and the water depth is calculated grid by grid to construct the underwater digital elevation model of the lake.
[0013] S7. Calculate the average lake surface elevation based on the lake surface photon point set extracted in step S2. Combine this with the underwater digital elevation model of the lake constructed in step S6, and use the volume integration method to calculate the total water volume of the lake.
[0014] The technical effect of the water volume calculation method based on spaceborne lidar data and geometric space constraints disclosed in this invention is that by introducing "geometric proportion constraints" and "adaptive polynomial fitting model", the overfitting and divergence problems of traditional methods when extrapolating sparse data are solved, and high-precision estimation of lake water volume without actual measured water depth data is achieved.
[0015] Furthermore, the specific formula for calculating the normalized difference water index in S1 is as follows:
[0016] ;
[0017] In the formula, For green band reflectivity, This refers to the reflectivity in the near-infrared band.
[0018] Furthermore, in step S2, the specific steps for extracting photon points using the kernel density estimation method include: constructing the probability density function of the photon point cloud in the vertical elevation direction. :
[0019] ;
[0020] in, This is the elevation value. For the first The elevation of a photon The total number of photons, For bandwidth, The Gaussian kernel function;
[0021] By using peak detection, the elevation point set corresponding to the first peak of the probability density function is extracted as the lake surface photon point set, and the elevation point set corresponding to the second peak located below the lake surface is extracted as the lake bottom photon point set.
[0022] Furthermore, the calculation formula for the two-way refraction correction based on Snell's law in S3 is as follows:
[0023] ;
[0024] in, To correct for the true water depth, For satellite observation of depth, The laser incident angle, For the angle of refraction, and These are the refractive indices of air and water, respectively.
[0025] Furthermore, S4 specifically includes:
[0026] S41, Calculate the geometric centroid C of the lake vector polygon;
[0027] S42, perform a 360-degree rotation scan with the centroid C as the origin, calculate the Euclidean distance from all points on the lake boundary to the centroid C, and define the maximum value as the geometric radius of the lake. ;
[0028] S43, calculate the Euclidean distance from all points in the measured water depth trajectory data to the centroid C, and define the maximum value as the trajectory coverage radius. ;
[0029] S44, Constructing Geometric Scale Factors ;
[0030] S45, obtain the maximum water depth in the measured trajectory. Combined with safety redundancy coefficient Calculate the upper limit of the maximum depth during the fitting process. :
[0031] ;
[0032] Wherein, the security redundancy coefficient The value is determined based on the spatial distribution characteristics of the remote sensing water depth index of the corresponding lake.
[0033] Furthermore, the geometric scaling factor in S44 Based on the topographic features of the lake, linear proportions, exponential functions, or power functions are used for fitting and adjustment.
[0034] Furthermore, the safety redundancy factor in S45 The determination method includes: separately statistically analyzing the set of remote sensing water depth indices within the coverage area of the ICESat-2 measured trajectory and the set of remote sensing water depth indices within the vector boundary of the entire lake, and extracting the feature statistics of both, wherein the feature statistics are the 95th quantiles; the safety redundancy coefficient The value of is determined by the ratio of the relative depth trend of the entire lake to the measured depth trend of the local area, after being normalized by a geometric scaling factor.
[0035] Furthermore, the polynomial fitting model in S5 is:
[0036] ;
[0037] in, To predict water depth, For distance from the center of mass Normalized distance, The maximum depth of the lake center is to be determined. Let be the morphological index to be solved, and be . Coefficients to be solved;
[0038] The least squares method is used to solve for the parameters of the polynomial fitting model. and When, the following inequality constraint is applied:
[0039] .
[0040] Furthermore, S6 specifically includes:
[0041] S61, rasterize the internal region of the lake vector polygon and determine the area of the raster unit;
[0042] S62, calculate the Euclidean distance from each grid cell to the nearest lake shore boundary, and generate the whole lake distance field matrix;
[0043] S63, Obtain the maximum value in the distance field matrix. Divide all values in the matrix by Normalization is performed to obtain the normalized range field matrix. ;
[0044] S64, Fit the optimal parameters obtained in step S5 and Substituting into the polynomial fitting model, for the matrix The water depth distribution matrix of the entire lake was obtained by performing element-by-element calculations.
[0045] Furthermore, in step S7, the specific formula for calculating the total water volume of the lake using the volume integral method is as follows:
[0046] ;
[0047] in, For the lake's water volume, The total number of grid cells within the lake area. The elevation of the lake surface. For the first The lakebed elevation of each grid, This refers to the area of a grid cell; only statistics are presented. Effective water grid. Attached Figure Description
[0048] Figure 1 An illustration of the water depth extraction effect from ICESat-2 photon data is provided for an embodiment of the present invention.
[0049] Figure 2 This is a schematic diagram illustrating the principle of geometric constraints and maximum water depth estimation based on trajectory range.
[0050] Figure 3 A schematic diagram of lake water depth data fitting;
[0051] Figure 4 A schematic diagram of the 3D underwater topography of a lake. Detailed Implementation
[0052] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. 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.
[0053] To address the problems of poor model adaptability and the tendency to generate unreasonable extrapolation values in areas without data when using sparse satellite trajectory data to extrapolate the underwater topography of a whole lake in existing technologies, this invention provides a water volume calculation method based on spaceborne lidar data and geometric spatial constraints. Figures 1 to 4 As shown, the specific steps include:
[0054] Step S1: Obtaining the Lake Boundary. This step aims to obtain a high-precision two-dimensional lake boundary, providing a foundation for subsequent geometric constraints and spatial mapping.
[0055] First, select multispectral optical remote sensing images that cover the target lake and are cloudless or have few clouds, such as Sentinel-2 MSI images or Landsat 8 OLI images. Taking Sentinel-2 images as an example, its spatial resolution is 10 meters, which can provide fine boundary information.
[0056] Then, the Normalized Difference Water Index (NDWI) is calculated. The formula is as follows:
[0057] ;
[0058] In the formula, For green band reflectivity, This refers to near-infrared reflectance. NDWI can effectively enhance water body information while suppressing vegetation and soil information.
[0059] Finally, a suitable threshold (e.g., 0) is set to binarize the NDWI image. Pixels with NDWI values greater than the threshold are identified as water bodies, while others are considered non-water bodies. Morphological processing (such as dilation and erosion) is then applied to the binary image to fill in tiny voids within the water bodies and smooth the boundaries. The processed water body mask is converted to Shapefile vector polygon format, generating a closed, high-precision lake vector polygon. This vector polygon precisely defines the lake's land-water boundary, providing spatial constraints for subsequent calculations of the geometric centroid, distance field, etc., ensuring that all subsequent calculations are performed within the lake region and avoiding the introduction of external noise.
[0060] Step S2: Photon point cloud processing. This step aims to accurately separate the effective signal photons representing the lake surface and lake bottom from the raw photon point cloud data of ICESat-2.
[0061] First, acquire ICESat-2 ATL03 level data covering the target lake. ATL03 data provides geographic coordinates and elevation information for each photon, but it contains a significant amount of atmospheric noise, solar background noise, etc.
[0062] Then, in order to extract the effective signal from the noise, this invention employs kernel density estimation (KDE). Specifically, along the satellite orbit direction, a probability density function is constructed in the vertical elevation direction for the photon point cloud within a certain spatial window. :
[0063] ;
[0064] in, This is the elevation value. For the first The elevation of a photon The total number of photons, For bandwidth (which can be determined through adaptive algorithms, such as those based on Scott's rules). This is the Gaussian kernel function.
[0065] In the vertical elevation direction, the water surface typically reflects strongly, forming a dense cloud of photons, corresponding to the first peak of the probability density function. After the laser pulse penetrates the water, it reflects at the lake bottom, forming a weaker second peak (located below the lake surface). Peak detection algorithms can accurately identify the precise elevations of these two peaks. The elevation point set corresponding to the first peak is extracted as the lake surface photon point set, and the elevation point set corresponding to the second peak is extracted as the lake bottom photon point set. Using the KDE method, background noise can be effectively filtered, and even under low signal-to-noise ratio conditions, the weak underwater signal can be stably extracted, providing a reliable data source for subsequent acquisition of high-precision measured water depth data.
[0066] from Figure 2 This is a two-dimensional scatter plot, with the horizontal axis representing the distance (or latitude and longitude) traveled by the satellite along its orbit, and the vertical axis representing elevation. The plot shows the raw photon point cloud distribution acquired by the ICESat-2 satellite as it traverses a lake. Two main photon aggregation zones are clearly visible in the plot:
[0067] The dense bands at the top (A) correspond to the set of photons on the lake surface. Due to the strong reflection from the water surface, the photon density is highest here, forming the first peak of the probability density function.
[0068] The diffuse stripe (B) below corresponds to the photon set at the lake bottom. After the laser pulse penetrates the water, it undergoes diffuse reflection at the lake bottom, resulting in a weak signal, low photon density, and some fluctuations, corresponding to the second peak of the probability density function.
[0069] In addition, the figure also contains a large number of isolated noise points (C), which are interference signals such as atmospheric scattering and solar background noise. Using the kernel density estimation (KDE) method described in step S2 of this invention, these two main peaks can be effectively separated in the vertical elevation direction, thereby accurately extracting the lake surface point set and the lake bottom point set. Different colors or markers are used in the figure to distinguish between the extracted signal photons and noise photons. Figure 1 This intuitively demonstrates the effectiveness of the present invention in extracting weak underwater signals from complex background noise, laying the foundation for obtaining accurate measured depth data in the future.
[0070] Step S3: Water Depth Geometric Correction. This step aims to eliminate the refraction effect of laser light in the water, converting the slant range apparent depth observed by satellite into the true vertical water depth.
[0071] Laser pulses propagate in a straight line in the atmosphere, but their propagation direction changes upon entering water because water has a higher refractive index than air. The "water depth" recorded by ICESat-2 is actually the slant path length of the laser in the water. Therefore, geometric correction must be performed according to Snell's law. The correction formula is:
[0072] ;
[0073] in, To correct for the true water depth, (For satellite observation depth, i.e., the slant distance along the orbit from the photon on the water surface to the photon on the bottom). The laser incident angle, For the angle of refraction, and These are the refractive indices of air and water, respectively.
[0074] After this correction, each lakebed photon point obtains a physically accurate vertical water depth value corresponding to its geographical location, generating a trajectory dataset containing latitude, longitude, and measured water depth. This step is a crucial prerequisite for ensuring the absolute accuracy of subsequent terrain modeling and water depth fitting, avoiding systematic biases caused by refraction.
[0075] Step S4: Geometric Constraint Boundary Construction. This step is one of the core innovations of this invention. It aims to utilize the geometric morphology information of the lake to provide a physically reasonable upper limit for the maximum water depth for subsequent model fitting, preventing the model from over-extrapolating in sparse data regions (especially the center of the lake) and generating unrealistic "abyss values".
[0076] The specific implementation process is as follows:
[0077] S41. Using GIS spatial analysis algorithms, calculate the geometric centroid C of the lake's vector polygon. The geometric centroid is the geometric center of the lake's planar shape; for regular lakes, it is close to the physical center of gravity.
[0078] S42, with the centroid C as the origin, performs a 360-degree omnidirectional scan to calculate the Euclidean distance from all vertices on the lake boundary to the centroid C. The maximum value of these distances is defined as the geometric radius of the lake. This radius represents the maximum extent of the lake's planar shape in the radial direction.
[0079] S43, calculate the Euclidean distance from all points in the measured water depth trajectory data (i.e., the data points processed in step S3) to the centroid C. Take the maximum value among these distances and define it as the trajectory coverage radius. This radius represents the radial extent of the lake area actually detected by the satellite trajectory.
[0080] S44, Constructing Geometric Scale Factors This scaling factor quantifies the relationship between the overall size of the lake and the coverage of known data. If the trajectory only covers the lake's edge, this factor will be much greater than 1. Depending on the lake's topographic features, this geometric scaling factor can be fitted and adjusted using linear scaling, exponential functions, or power functions to adapt to the morphological characteristics of different lake basins.
[0081] S45, obtain the maximum water depth in the measured trajectory. Combined with safety redundancy coefficient Calculate the upper limit of the maximum depth during the fitting process. :
[0082] ;
[0083] Wherein, the security redundancy coefficient The value of is determined based on the spatial distribution characteristics of the remote sensing depth index (SDB, such as relative depth retrieved from satellite imagery) of the corresponding lake. Specifically, the SDB index set within the coverage area of the ICESat-2 measured trajectory and the SDB index set within the vector boundary of the entire lake are statistically analyzed separately, and the characteristic statistics of both are extracted (e.g., the 95th percentile is selected to exclude extreme noise). The value of SDB is determined by normalizing the ratio of the overall lake's relative depth trend to the local measured depth trend using a geometric scaling factor. For example, if the overall lake's SDB distribution shows that the depth in the central region is much greater than that at the edges, then... It will be greater than 1, and vice versa. In this way, It is not a simple mathematical extrapolation value, but a physical upper limit that integrates the planar geometry of the lake and the known depth distribution trend, providing a reliable boundary constraint for subsequent fitting.
[0084] Figure 2 The geometric principle of the core innovative step S4 (geometric constraint boundary construction) of this invention is illustrated in the figure. Much larger In this case, the satellite trajectory does not cover the deepest part of the lake (the center). According to this invention, a geometric scaling factor quantifies the degree of this information gap. This is combined with the measured maximum water depth along the trajectory. This allows us to calculate the physical upper limit of the lake's maximum depth. The diagram visually illustrates how the planar geometry of the lake can be used to constrain and infer the water depth in unknown areas (especially the center of the lake), effectively preventing excessive divergence during model extrapolation.
[0085] Step S5: Water Depth Cross-Section Fitting. This step aims to construct a mathematical function with good generalization ability that can describe the cross-sectional morphology of the lake based on discrete measured water depth points.
[0086] First, calculate the Euclidean distance from all measured water depth points to the geometric centroid C, and divide it by the geometric radius of the lake. After normalization, the normalized distance d is obtained, and its value range is: . Represents the position of the center of mass. Represents the boundary of a lake.
[0087] Then, a system is established to predict water depth using the normalized distance d as the independent variable. This is a multinomial fitting model for the dependent variable. The model used in this invention is:
[0088] ;
[0089] in, To predict water depth, For distance from the center of mass Normalized distance, The maximum depth of the lake center is to be determined. Let be the morphological index to be solved, and be . The coefficients to be solved (controlling the steepness of the cross-sectional curve).
[0090] Finally, the least squares method is used to solve for the parameters of the polynomial fitting model. and When, the following inequality constraint is applied:
[0091] .
[0092] This constraint ensures that the calculated maximum depth of the lake center does not exceed physical geometric limits, effectively avoiding the problem of over-extrapolation caused by sparse measured data (e.g., the trajectory does not pass through the deepest point). Through this step, the optimal combination of model parameters that best fits the measured data and conforms to physical laws is obtained.
[0093] Figure 3 The process and effect of step S5 (water depth profile fitting) of the present invention are shown. The figure is a two-dimensional coordinate system, with the horizontal axis representing the normalized distance d from the centroid C (ranging from 0 to 1) and the vertical axis representing the water depth D (positive values indicate depth). Figure 3 The fitting results with and without constraints are visually compared. The fitted curve with constraints (C) is confined within a reasonable physical upper limit in the lake center region, exhibiting a smooth overall shape that conforms to the physical law of the lake cross-section decreasing from the lake center to the shore, without exhibiting the illogical "abyss value" extrapolation seen in the unconstrained curve (B). This figure vividly demonstrates the crucial role of the constraint mechanism proposed in this invention in ensuring the physical rationality of the model.
[0094] Step S6: Reconstructing the 3D Topography of the Entire Lake. This step aims to extend the 2D cross-sectional model obtained in Step S5 to the 2D planar space of the entire lake, constructing a continuous underwater topography of the entire lake. Specifically, it includes the following sub-steps:
[0095] S61, the internal region of the lake vector polygon is rasterized to determine the area of each raster cell; each raster cell has a fixed area. For example, the resolution of the original remote sensing image can be set to 10 meters × 10 meters.
[0096] S62, calculate the Euclidean distance from each grid cell to the nearest lake shore boundary, and generate the whole lake distance field matrix.
[0097] S63, Obtain the maximum value in the distance field matrix. This represents the maximum distance from the center of the lake to the shore. Divide all values in the matrix by... Normalization is performed to obtain the normalized range field matrix. Each value in this matrix also falls within... Within the interval, it has the exact same physical meaning as the normalized distance d in step S5 (offshore normalized distance).
[0098] S64, Fit the optimal parameters obtained in step S5 and Substituting into the polynomial fitting model, for the matrix A stepwise calculation is performed to obtain the water depth distribution matrix of the entire lake. By mapping the one-dimensional cross-sectional curve to a two-dimensional plane through a distance field, this method can generate a smooth three-dimensional terrain that conforms to the geometry of the lake, avoiding the "bull's-eye" or "step" effects that may be produced by traditional interpolation methods.
[0099] Figure 4 This image showcases the final result of step S6 (full lake 3D terrain reconstruction) of the present invention. It is a 3D perspective rendering that vividly presents the three-dimensional shape of the lake. Figure 4 The final output of this invention is a continuous, smooth, and geometrically constrained three-dimensional underwater topographic model of the entire lake. Based on this model, the total water volume of the lake can be accurately calculated using the volume integration method described in step S7. The figure visually illustrates the complete technical path from one-dimensional sparse trajectory data to two-dimensional cross-section fitting, and then to three-dimensional full-lake topographic reconstruction, highlighting the invention's ability to perform high-precision three-dimensional topographic modeling under sparse data conditions.
[0100] Step S7: Lake Water Volume Estimation. This step uses the constructed underwater topographic model (DEM) to calculate the total water volume of the lake. The volume integration method is employed, which involves accumulating the volume of all effective water cells within the lake. The calculation formula is:
[0101] ;
[0102] in, For the lake's water volume, The total number of grid cells within the lake area. The elevation of the lake surface. For the first The lakebed elevation of each grid, This refers to the area of a grid cell; only statistics are presented. The effective water grid is determined. During the calculation, only those satisfying the criteria are counted. The effective water grid automatically eliminates invalid areas where the lake bottom is higher than the lake surface, which may be caused by model errors, thus ensuring the physical accuracy and stability of water volume estimation.
[0103] The beneficial effects of this invention are:
[0104] 1. By using geometric scale constraints, the algorithm is effectively prevented from calculating "abyss values" or "overfitting values" in areas without data (especially in the center of a lake), thus ensuring the numerical stability of the terrain model.
[0105] 2. Polynomial fitting models use morphological indices Automatic optimization can adaptively simulate a cone ( From frying pan type ( The various lake morphologies are better than parabolic fitting.
[0106] 3. This method provides confidence intervals based on physical boundaries while ensuring the accuracy of macroscopic water volume estimation, making it particularly suitable for monitoring the storage of lakes in areas without actual measurement data, such as the Qinghai-Tibet Plateau.
[0107] Example embodiments have been disclosed herein, and while specific terminology has been used, it is for illustrative purposes only and should be construed as such, and is not intended to be limiting. In some instances, it will be apparent to those skilled in the art that features, characteristics, and / or elements described in conjunction with particular embodiments may be used alone, or in combination with features, characteristics, and / or elements described in conjunction with other embodiments, unless otherwise expressly indicated. Therefore, those skilled in the art will understand that various changes in form and detail may be made without departing from the scope of the invention as set forth in the appended claims.
Claims
1. A method for calculating water volume based on spaceborne lidar data and geometric spatial constraints, characterized in that, include: S1. Acquire multispectral optical remote sensing images covering the target lake, calculate the normalized differential water index, extract the lake's land-water boundary through threshold segmentation, and generate a closed lake vector polygon. S2, acquire ICESat-2 ATL03 photon point cloud data covering the target lake, separate background noise in the vertical elevation direction using kernel density estimation, and extract the photon point set on the lake surface and the photon point set on the lake bottom; S3, based on Snell's law, performs two-way refraction correction on the extracted lake bottom photon point set to obtain trajectory data containing geographic coordinates and measured water depth; S4. Using the geometric centroid of the vector polygon of the lake as the origin, perform an all-round scan to obtain the maximum geometric radius and corresponding water depth of the lake, as well as the maximum coverage radius and corresponding water depth of the measured trajectory; construct a geometric scale factor, and combine it with the measured maximum water depth to calculate the physical limit value of the maximum depth of the entire lake; S5. Calculate the normalized distance from each measured water depth point to the geometric centroid, establish a polynomial fitting model with the normalized distance as the independent variable and the water depth as the dependent variable, optimize the parameters using the least squares method, and use the physical constraint value as the inequality constraint condition of the solution space. S6. The interior of the lake vector polygon is rasterized, and a normalized distance field of the whole lake is generated based on Euclidean distance transformation. The optimal parameters obtained in step S5 are substituted into the polynomial fitting model, and the water depth is calculated grid by grid to construct the underwater digital elevation model of the lake. S7. Calculate the average lake surface elevation based on the lake surface photon point set extracted in step S2. Combine this with the underwater digital elevation model of the lake constructed in step S6, and use the volume integration method to calculate the total water volume of the lake.
2. The method according to claim 1, characterized in that, The specific formula for calculating the Normalized Differential Water Index (NDWI) in S1 is as follows: ； In the formula, For green band reflectivity, This refers to the reflectivity in the near-infrared band.
3. The method according to claim 1, characterized in that, In step S2, the specific steps for extracting photon points using the kernel density estimation method include: constructing the probability density function of the photon point cloud in the vertical elevation direction. : ； in, This is the elevation value. For the first The elevation of a photon The total number of photons, For bandwidth, The Gaussian kernel function; By using peak detection, the elevation point set corresponding to the first peak of the probability density function is extracted as the lake surface photon point set, and the elevation point set corresponding to the second peak located below the lake surface is extracted as the lake bottom photon point set.
4. The method according to claim 1, characterized in that, The calculation formula for two-way refraction correction based on Snell's law in S3 is as follows: ； in, To correct for the true water depth, For satellite observation of depth, The laser incident angle, For the angle of refraction, and These are the refractive indices of air and water, respectively.
5. The method according to claim 1, characterized in that, S4 specifically includes: S41, Calculate the geometric centroid C of the lake vector polygon; S42, perform a 360-degree rotation scan with the centroid C as the origin, calculate the Euclidean distance from all points on the lake boundary to the centroid C, and define the maximum value as the geometric radius of the lake. ; S43, calculate the Euclidean distance from all points in the measured water depth trajectory data to the centroid C, and define the maximum value as the trajectory coverage radius. ; S44, Constructing Geometric Scale Factors ; S45, obtain the maximum water depth in the measured trajectory. Combined with safety redundancy coefficient Calculate the upper limit of the maximum depth during the fitting process. : ； Wherein, the security redundancy coefficient The value is determined based on the spatial distribution characteristics of the remote sensing water depth index of the corresponding lake.
6. The method according to claim 5, characterized in that, The geometric scaling factor in S44 Based on the topographic features of the lake, linear proportions, exponential functions, or power functions are used for fitting and adjustment.
7. The method according to claim 5, characterized in that, The safety redundancy factor in S45 The determination method includes: separately statistically analyzing the set of remote sensing water depth indices within the coverage area of the ICESat-2 measured trajectory and the set of remote sensing water depth indices within the vector boundary of the entire lake, and extracting the feature statistics of both, wherein the feature statistics are the 95th quantiles; the safety redundancy coefficient The value of is determined by the ratio of the relative depth trend of the entire lake to the measured depth trend of the local area, after being normalized by a geometric scaling factor.
8. The method according to claim 1, characterized in that, The polynomial fitting model in S5 is as follows: ； in, To predict water depth, For distance from the center of mass Normalized distance, The maximum depth of the lake center is to be determined. Let be the morphological index to be solved, and be . Coefficients to be solved; The least squares method is used to solve for the parameters of the polynomial fitting model. and When, apply the following inequality constraint: 。 9. The method according to claim 1, characterized in that, S6 specifically includes: S61, rasterize the internal region of the lake vector polygon and determine the area of the raster unit; S62, calculate the Euclidean distance from each grid cell to the nearest lake shore boundary, and generate the whole lake distance field matrix; S63, Obtain the maximum value in the distance field matrix. Divide all values in the matrix by Normalization is performed to obtain the normalized range field matrix. ; S64, Fit the optimal parameters obtained in step S5 and Substituting into the polynomial fitting model, for the matrix The water depth distribution matrix of the entire lake was obtained by performing element-by-element calculations.
10. The method according to claim 1, characterized in that, In S7, the specific formula for calculating the total water volume of the lake using the volume integral method is as follows: ； in, For the lake's water volume, The total number of grid cells within the lake area. The elevation of the lake surface. For the first The lakebed elevation of each grid, This refers to the area of a grid cell; only statistics are presented. Effective water grid.