A model construction and water reserve estimation method based on turning point identification
By identifying water body turning points, removing outlier data, and performing regression fitting, the problems of water storage model accuracy and automated updating were solved, achieving high-precision, batch-based water storage estimation and meeting the dynamic monitoring needs of water resource management.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- JIANGSU PROVINCE SURVEYING & MAPPING ENG INST
- Filing Date
- 2026-06-11
- Publication Date
- 2026-07-14
Smart Images

Figure CN122391513A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of water resources monitoring and geographic information technology, specifically to a model construction and water storage estimation method based on inflection point identification. Background Technology
[0002] Water resources are a key element supporting regional economic and social development and ecological balance. Accurate and dynamic monitoring of water storage in lakes, reservoirs, and other water bodies is a crucial foundation for water resource allocation and management, flood and drought early warning, and ecological flow assurance. Traditional methods of obtaining water storage rely on manual cross-sectional measurements or water level-storage curves of single water bodies. These methods are not only inefficient but also difficult to update in batches across large areas and multiple water bodies, failing to meet the dynamic monitoring needs of basin-level water resource management departments.
[0003] With the development of geographic information technology, constructing a three-dimensional model of a body of water using underwater topographic survey data and then generating a water body area-water storage lookup table has become a common technique for estimating water storage. Specifically, discrete underwater elevation points of the water body are obtained using a single-beam or multi-beam echo sounder, and a three-dimensional model of the water body is generated through interpolation. Then, the contour line method or inundation analysis algorithm is used to calculate the water surface area and water storage corresponding to different water levels (depths), constructing an area-storage relationship curve. Finally, a mathematical model between area and water storage is established through regression fitting, realizing the estimation from water surface area to water storage.
[0004] However, the aforementioned method of directly fitting the area-water storage relationship based on full-scale 3D model data of water bodies has significant technical shortcomings in practical applications. Specifically, underwater topographic data of 248 lakes and 476 reservoirs in a certain area were measured using a single-beam echo sounder, and 3D models and area-water storage relationship curves for each water body were constructed based on the measured data. Verification results show that the mathematical model constructed using traditional methods, i.e., regression fitting using all measured data points, falls far short of the accuracy requirements for engineering applications: only about 6.5% of the 724 water bodies meet the accuracy standard of a relative error of less than 10%. This means that the mathematical model constructed using existing methods has significant estimation biases in more than 90% of the water bodies, making it difficult to support the decision-making needs of water resource management departments.
[0005] Further analysis revealed that the root cause of the aforementioned accuracy problem lies in the fact that traditional methods assume a strict monotonic functional mapping relationship between area and volume, and assume that all measured data points satisfy this mathematical premise. However, in real-world scenarios, especially for water bodies constrained by dams, artificial boundaries, or natural terrain, once the water level rises to a certain height, the water surface area, limited by the boundary, no longer increases significantly with rising water levels, while the water volume continues to grow. At this point, an anomaly of a one-to-many mapping appears on the area-water storage relationship curve, where the same water surface area value corresponds to multiple different water storage values. This violates the single-valued function premise required for regression fitting. Including these anomalous data points in the fitting will lead to distortion of the mathematical model.
[0006] More importantly, the aforementioned inherent defects are amplified dramatically under specific geographical conditions. The middle and lower reaches of the Yangtze River Plain, exemplified by a certain region, are dotted with numerous shallow lakes and reservoirs. In such environments, the relationship between water area and storage capacity exhibits extreme sensitivity: due to the extremely shallow water depth, even a tiny expansion or contraction of the water surface boundary (e.g., a change in the boundary of 1-2 pixels in a remote sensing image) can lead to a significant relative change in the calculated water storage capacity. In stark contrast, in some deep lakes, where depths can reach tens or even hundreds of meters, even a similar absolute change in surface area results in a much smaller relative rate of change in water volume compared to shallow water bodies. Therefore, while existing methods can achieve acceptable relative accuracy in deep lakes, in certain shallow plain areas, the relative error rises sharply, severely limiting the practicality of the mathematical models.
[0007] Furthermore, one of the core business needs of water resource management departments is to achieve annual dynamic updates of water storage to promptly grasp changes in the spatiotemporal distribution of water resources and support decision-making in water allocation, ecological water replenishment, and flood and drought control. However, existing methods are insufficient to support this high-frequency, large-scale, and routine update task due to problems such as the need for manual selection of effective data domains or adjustment of fitting parameters for each water body, insufficient accuracy, and inability to process data in batches. If traditional methods are used to model and update each of the 724 water bodies in a certain region one by one, not only will the workload be large and the cycle long, but the accumulation of errors from numerous low-precision mathematical models will also lead to a loss of credibility in the storage estimation results. Therefore, constructing a method that can build high-precision, batch, and automated mathematical models of water bodies is an urgent need to serve the annual operational update of water storage in a certain region.
[0008] No effective solutions have yet been proposed to address the problems in the relevant technologies. Summary of the Invention
[0009] In response to the problems in related technologies, this invention proposes a model construction and water storage estimation method based on inflection point identification to overcome the aforementioned technical problems existing in the existing related technologies.
[0010] Therefore, the specific technical solution adopted by the present invention is as follows: A method for model construction and water storage estimation based on inflection point identification includes: The discrete underwater elevation point cloud data of all water bodies to be measured within the target area is acquired, and the discrete underwater elevation point cloud data is preprocessed; based on the preprocessing results, a three-dimensional water body model of each water body is constructed. Based on the three-dimensional model of the water body, the water level elevation range and water level step size are determined, and an equally spaced water level sequence is generated. Based on the equally spaced water level sequence, the water surface area and water storage corresponding to each water level are calculated respectively, and the water surface area and water storage corresponding to each water level are paired to generate a water body area-water storage relationship dataset and a water body area-water storage relationship curve. Based on the water body area-water storage relationship dataset, a preset inflection point identification strategy is adopted to automatically identify the water body area-water storage relationship curve, determine the topographic inflection point on the area-water storage relationship curve that represents the boundary of the effective data domain, and obtain the water surface area and water storage corresponding to the topographic inflection point. Using the water surface area corresponding to the topographic inflection point as the threshold, outlier data points with water surface areas greater than the threshold are removed, and only valid data points within the valid data domain are retained. The number of valid data points is then checked to see if it meets the minimum fitting sample size requirement. Based on the valid data points within the retained valid data domain, a fitting function is selected from the preset candidate model library, and the least squares method is used for regression fitting to construct a high-precision area-water storage mathematical model of the water body. The optimal function form, fitting parameters and applicable domain of the high-precision area-water storage mathematical model are determined. For each water body to be measured within the target area, a corresponding area-water storage mathematical model is constructed to complete batch automated modeling. Based on the latest water surface area obtained from remote sensing interpretation, the corresponding high-precision area-water storage mathematical model is called to calculate the water storage of each water body in batches and generate a water storage update report.
[0011] The beneficial effects of this invention are as follows: 1. This invention identifies the topographic inflection point of the area-water storage relationship curve, uses this inflection point as the dividing criterion, eliminates abnormal data domains with mapping anomalies, and retains only valid data points that satisfy the single-value mapping relationship for modeling. This fundamentally solves the problem of mathematical mapping anomalies in full data fitting in the prior art, ensures the accuracy and validity of the modeling data, avoids estimation deviations caused by abnormal data interference, and significantly improves the model fitting accuracy.
[0012] 2. This invention abandons the drawbacks of traditional piecewise fitting modeling, and clearly defines the turning point as the boundary between the effective data domain and the outlier data domain. Its corresponding geographical significance is clear and it has clear physical interpretability. Compared with existing piecewise fitting methods, it does not require complex parameter settings, simplifies the modeling process, and ensures the uniformity and stability of the model, making it more practical for engineering.
[0013] 3. This invention includes a candidate model library with multiple function forms, which combines the least squares method to complete parameter estimation, adapts to different water body morphology characteristics, and provides multiple inflection point identification strategies that can be flexibly selected according to data quality and water body type, adapting to the modeling needs of different scenarios, avoiding the limitations of a single fitting method, and further improving the reliability of water storage estimation.
[0014] 4. This invention enables batch automated modeling of all water bodies to be measured within the target area. It uses a loop traversal method to execute the entire process sequentially without manual intervention. At the same time, it extracts water surface data through accurate remote sensing image interpretation and combines it with the constructed mathematical model to complete water storage estimation, efficiently outputting standardized results, meeting the needs of normalized and operational water resource management, and providing stable support for subsequent dynamic monitoring of water storage. Attached Figure Description
[0015] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0016] Figure 1 This is a flowchart of a method for model construction and water storage estimation based on inflection point identification according to an embodiment of the present invention; Figure 2 This is a schematic diagram comparing the fitting curves of three strategies according to an embodiment of the present invention; Figure 3 This is a schematic diagram illustrating the position markings of three strategic turning points according to embodiments of the present invention; Figure 4 This is a schematic diagram of a data simulation example using the second derivative method according to an embodiment of the present invention. Detailed Implementation
[0017] To further illustrate the various embodiments, the present invention provides accompanying drawings, which are part of the disclosure of the present invention. These drawings are mainly used to illustrate the embodiments and can be used in conjunction with the relevant descriptions in the specification to explain the operating principles of the embodiments. With reference to these drawings, those skilled in the art should be able to understand other possible implementation methods and the advantages of the present invention.
[0018] According to an embodiment of the present invention, a method for model construction and water storage estimation based on inflection point identification is provided.
[0019] The present invention will now be further described in conjunction with the accompanying drawings and specific embodiments, such as... Figure 1 As shown, the model construction and water storage estimation method based on inflection point identification according to an embodiment of the present invention includes: S1. Obtain discrete underwater elevation point cloud data of all water bodies to be measured within the target area, and preprocess the discrete underwater elevation point cloud data; based on the preprocessing results, construct a three-dimensional water body model for each water body. S2. Based on the three-dimensional model of the water body, determine the water level elevation range and water level step size, and generate an equally spaced water level sequence; based on the equally spaced water level sequence, calculate the water surface area and water storage corresponding to each water level, and pair the water surface area and water storage corresponding to each water level to generate a water body area-water storage relationship dataset and a water body area-water storage relationship curve. S3. Based on the water body area-water storage relationship dataset, a preset inflection point identification strategy is adopted to automatically identify the water body area-water storage relationship curve, determine the topographic inflection point on the area-water storage relationship curve that represents the boundary of the effective data domain, and obtain the water surface area and water storage corresponding to the topographic inflection point. S4. Using the water surface area corresponding to the topographic turning point as the threshold, remove abnormal data points with a water surface area greater than the threshold, retain only the valid data points within the valid data domain, and check whether the number of valid data points meets the minimum fitting sample size requirement. S5. Based on the valid data points within the retained valid data domain, select a fitting function from the preset candidate model library, use the least squares method for regression fitting, construct a high-precision area-water storage mathematical model of the water body, and determine the optimal function form, fitting parameters and applicable domain of the high-precision area-water storage mathematical model. S6. Construct corresponding area-water storage mathematical models for all water bodies to be measured within the target area to complete batch automated modeling; based on the latest water surface area obtained from remote sensing interpretation, call the corresponding high-precision area-water storage mathematical models to calculate the water storage of each water body in batches and generate a water storage update report.
[0020] In this optional embodiment, discrete underwater elevation point cloud data of all water bodies to be measured within the target area are acquired, and the discrete underwater elevation point cloud data are preprocessed. Based on the preprocessing results, a three-dimensional water body model of each water body is constructed, including: underwater topographic measurement of lakes and reservoirs using a single-beam echo sounder; setting up measurement cross-sections or survey line grids according to the water body morphology; collecting water surface elevation, bottom elevation, water depth, and plane coordinates to form discrete underwater elevation point cloud data containing plane coordinates and elevation values; removing outliers from the point cloud data based on the local neighborhood Laida criterion; and using a moving window weighted average method to smooth areas with slight noise. The system performs sliding processing and unifies the coordinates of all water bodies to a preset geodetic coordinate system (i.e., CGCS2000 National Geodetic Coordinate System) and a preset elevation datum (i.e., 1985 National Elevation Datum), resulting in preprocessed discrete point cloud data. The preprocessed discrete point cloud data is then used to identify and densify topographic feature points, generating encrypted feature points. The original underwater elevation measurement points, encrypted feature points, and boundary points from the water space survey results are then merged to form a unified discrete point set. Based on this discrete point set, an irregular triangular network is constructed using a triangulation algorithm. Finally, the irregular triangular network is interpolated into a regular grid three-dimensional water body model using the natural neighborhood method.
[0021] In this optional embodiment, constructing an irregular triangular mesh based on a set of discrete points using a triangulation algorithm includes: constructing a super triangle containing all discrete points as the initial mesh for triangulation; sequentially inserting each discrete point into the current triangular mesh, finding the triangle where the newly inserted point is located, and splitting the triangle into three new triangles; verifying the triangulation conditions edge by edge, and performing edge flipping on edges that do not satisfy the empty circumcircle criterion until the triangular mesh satisfies the empty circumcircle criterion; performing Laplace smoothing optimization on the generated triangular mesh, adjusting the vertex positions to improve the uniformity of the triangles, and completing the triangular mesh optimization according to a set number of iterations to obtain an irregular triangular mesh.
[0022] In this optional embodiment, the method of interpolating an irregular triangular mesh into a regular mesh 3D water body model using the natural neighborhood method includes: setting the target resolution of the 3D water body model and generating a regular mesh lattice on the horizontal projection plane; inserting each mesh point into the irregular triangular mesh, determining the Thiessen polygon cell where the mesh point is located, filtering to obtain the set of irregular triangular mesh vertices adjacent to the Thiessen polygon cell (i.e., the TIN vertex set), and using the irregular triangular mesh vertices as natural neighborhood points; calculating the Thiessen polygon area weight between the mesh point and each natural neighborhood point, and performing a weighted average of the elevation values of the natural neighborhood points according to the weights to obtain the interpolated elevation of the current mesh point; traversing all regular mesh points, sequentially completing the natural neighborhood filtering, weight calculation, and elevation interpolation to generate a complete regular elevation mesh; combining the generated regular elevation mesh with the boundary points in the water body spatial survey results, and uniformly storing them as a raster file format (i.e., a GeoTIFF format raster file) to obtain the 3D water body model for each water body.
[0023] In this optional embodiment, based on the three-dimensional model of the water body, the water level elevation range and water level step size are determined, and an equally spaced water level sequence is generated. Based on the equally spaced water level sequence, the water surface area and water storage corresponding to each water level are calculated, and the water surface area and water storage corresponding to each water level are paired to generate a water body area-water storage relationship dataset and a water body area-water storage relationship curve. This includes: statistically analyzing the elevation values of all grid cells in the three-dimensional model of the water body, obtaining the minimum and maximum elevation values, determining the water level elevation range, and setting a fixed water level step size within this range. The water level step size can be adaptively adjusted according to the actual conditions of different water bodies. Starting from the minimum elevation value, the water level step size is increased sequentially until the maximum elevation value is reached, generating an equally spaced water level sequence. For each water level value in the water level sequence, a flooding analysis algorithm is used to traverse all grid cells of the three-dimensional model of the water body. The algorithm compares the elevation values of raster cells with the current water level values. Based on the comparison results, the raster cells are labeled as water bodies and land areas respectively. The number of raster cells labeled as water bodies is counted and multiplied by the actual area represented by each raster cell labeled as a water body to calculate the water surface area of the corresponding water level. For each water level value in the water level sequence, the volume integral method is used to calculate the water column volume in each raster cell labeled as a water body. The water column volumes of all raster cells labeled as water bodies are summed to obtain the water storage of the corresponding water level. The water surface area and water storage of each water level are associated and paired one by one to construct a water body area-water storage relationship dataset. At the same time, a lookup table for water surface elevation, water depth, water surface area and water storage is generated. Scatter points are plotted with water surface area as the x-axis and water storage as the y-axis and connected to form a smooth curve to obtain the water body area-water storage relationship curve.
[0024] In this optional embodiment, based on the water body area-water storage relationship dataset, a preset inflection point identification strategy is used to automatically identify the water body area-water storage relationship curve, determine the topographic inflection points on the area-water storage relationship curve that represent the boundary of the effective data domain, and obtain the water surface area and water storage corresponding to the topographic inflection points. This includes: selecting a preset inflection point identification strategy to automatically identify the water body area-water storage relationship curve according to the data quality and water body morphology characteristics of the water body area-water storage relationship dataset; the inflection point identification strategy includes the slope change method, the cumulative growth rate method, and the second derivative method; when using the slope change method to identify inflection points, the discrete data points on the area-water storage relationship curve are sorted in ascending order of volume, and the central difference method combined with forward difference and backward difference methods is used to calculate the first derivative of the area with respect to volume for each data point; traversing all data points, confirming... To identify turning points, the method involves several steps. First, the point corresponding to the maximum value of the first derivative is determined. Then, points are searched sequentially along the direction of increasing volume. A slope descent threshold is used to filter out target turning points, and the corresponding water surface area and water storage are extracted. Second, the cumulative growth rate method is used to identify turning points. The maximum water surface area is determined centrally from the area-water storage relationship dataset. The cumulative growth rate of each data point's area relative to the maximum water surface area is calculated point by point. Data points are traversed in increasing volume order, and target turning points are filtered using a preset cumulative growth rate threshold. The corresponding water surface area and water storage are extracted. Third, the second derivative method is used to identify turning points. Based on the area-water storage relationship curve and the calculated first derivative, the second derivative of each data point is calculated. Discrete data points are traversed, and adjacent points where the second derivative changes from negative to positive are searched. Linear interpolation is used to locate the curve inflection point, which is then used as the target turning point. The corresponding water surface area and water storage are extracted.
[0025] In this optional embodiment, using the water surface area corresponding to the topographic inflection point as a threshold, outlier data points with water surface areas greater than the threshold are removed, retaining only valid data points within the valid data domain, and checking whether the number of valid data points meets the minimum fitting sample size requirement includes: defining three segments—free expansion zone, transition zone, and boundary constraint zone—based on the segmented characteristics of the area-water storage relationship curve; defining the interval with an area smaller than the water surface area corresponding to the topographic inflection point as the valid data domain with a strictly monotonically increasing function relationship, and defining the interval with an area greater than the water surface area corresponding to the topographic inflection point as the outlier data domain constrained by dams or natural shorelines and exhibiting one-to-many mapping anomalies; sorting the data points in the area-water storage relationship dataset in ascending order of volume, using the water surface area corresponding to the topographic inflection point as the discrimination threshold, comparing the water surface area of each data point with the threshold value, and identifying outlier data points with water surface areas greater than the threshold; removing all water surface areas greater than the threshold from the area-water storage relationship dataset. For outlier data points corresponding to the water surface area at the turning point, only valid data points with a water surface area less than or equal to the threshold and within the valid data domain are retained. Simultaneously, the number of removed data points, the total number of original data points, and the data removal ratio are calculated. The number of valid data points remaining after removing outlier data points is counted, and it is verified whether the number of valid data points meets the preset minimum fitting sample size requirement. If the number of valid data points meets the minimum fitting sample size requirement, subsequent modeling proceeds. If the number of valid data points does not meet the minimum fitting sample size requirement, a data insufficiency warning is issued, and the model is reverted to using the full dataset or a relaxed turning point threshold for fitting. The high-precision area-water storage mathematical model for the water body is marked as low confidence. After removing outlier data points, the remaining valid dataset, the area threshold and volume information corresponding to the turning point, data removal statistics, and the modeling domain corresponding to the valid dataset are output. This modeling domain is used as the applicable boundary for constructing the high-precision area-water storage mathematical model of the water body.
[0026] In this optional embodiment, based on valid data points within the retained valid data domain, a fitting function is selected from a pre-defined candidate model library, and regression fitting is performed using the least squares method to construct a high-precision area-water storage mathematical model of the water body. Determining the optimal function form, fitting parameters, and applicable domain of this high-precision area-water storage mathematical model includes: pre-defining a candidate model library containing multiple function forms and adapting it to different water body morphological characteristics; for each fitting function in the candidate model library, parameter estimation is performed using the least squares method based on the valid dataset; the fitting parameters are solved with the objective of minimizing the sum of squared residuals between the predicted values and actual observed values of the high-precision area-water storage mathematical model of the water body; linearization transformations are performed on power functions, exponential functions, and logarithmic functions to reconstruct the equation system and solve for the parameters; for polynomial functions, parameter estimation is completed by solving the normal equation system; and the fitting of each candidate fitting function is performed. The coefficients of determination were calculated for each result. The fitting function with the largest coefficient of determination was selected as the high-precision area-water storage mathematical model of the water body as the optimal high-precision area-water storage mathematical model, and the corresponding fitting parameters were determined. The water surface area corresponding to the topographic inflection point was used as the verification variable. The verification variable was substituted into the optimal high-precision area-water storage mathematical model to solve for the predicted water storage value. The true value of water storage for the corresponding area was retrieved from the original area-water storage lookup table. When the relative error between the predicted value and the true value does not exceed the preset value, the optimal high-precision area-water storage mathematical model is judged to meet the accuracy requirements. The constructed optimal high-precision area-water storage mathematical model was standardized and stored, and the water body identification information, the function form of the optimal high-precision area-water storage mathematical model, the fitting parameter values and the coefficient of determination were recorded simultaneously. The domain parameter information of the optimal high-precision area-water storage mathematical model was also retained.
[0027] In this optional embodiment, a corresponding area-water storage mathematical model is constructed for each water body to be measured within the target area to complete batch automated modeling. Based on the latest water surface area obtained from remote sensing interpretation, the corresponding high-precision area-water storage mathematical model is called to calculate the water storage of each water body in batches and generate a water storage update report. This includes: using programming software (preferably Python software), the operation of constructing an area-water storage mathematical model is performed for each water body to be measured within the target area through a loop to achieve fully automated batch modeling; after the area-water storage mathematical models corresponding to all water bodies within the target area are constructed, the result files are automatically generated in batches to form a complete area-water storage mathematical model library; using a normalized water body index, the latest water surface area data of each water body is obtained periodically, the corresponding area-water storage mathematical model of the water body is called to complete the water storage estimation in batches, the water storage data is updated synchronously, and a water storage update report is generated based on the updated water storage data.
[0028] In this optional embodiment, the use of a normalized water index to periodically acquire the latest water surface area data for each water body includes: periodically acquiring high-resolution optical remote sensing images of the target area, and sequentially performing radiometric calibration and atmospheric correction processing on the high-resolution optical remote sensing images; calculating the normalized water index based on the spectral reflectance of the green band and near-infrared band based on the radiometrically calibrated and atmospherically corrected high-resolution optical remote sensing images; determining the threshold of the normalized water index using a threshold method, and identifying areas where the threshold of the normalized water index is greater than or equal to a preset value as the boundaries of all water bodies in the study area; collecting vector layers of reservoirs and lakes from the latest annual land change survey data, associating and matching the identified water body boundaries with the land change survey vector boundaries in geographic information software, extracting the lake and reservoir water surface boundaries, and extracting the latest water surface area corresponding to each water body based on the lake and reservoir water surface boundaries.
[0029] According to another embodiment of the present invention, a model building and water storage estimation system based on inflection point identification is also provided, the system comprising: The water body 3D model construction module is used to acquire discrete underwater elevation point cloud data of all water bodies to be measured within the target area, and to preprocess the discrete underwater elevation point cloud data; based on the preprocessing results, a water body 3D model of each water body is constructed. The area-water storage relationship generation module is used to determine the water level elevation range and water level step size based on the three-dimensional model of the water body, and generate an equally spaced water level sequence; based on the equally spaced water level sequence, the water surface area and water storage corresponding to each water level are calculated respectively, and the water surface area and water storage corresponding to each water level are paired to generate the water body area-water storage relationship dataset and the water body area-water storage relationship curve. The inflection point identification module is used to automatically identify the area-water storage relationship curve of the water body based on the water body area-water storage relationship dataset, and determine the topographic inflection point on the area-water storage relationship curve that represents the boundary of the effective data domain, and obtain the water surface area and water storage corresponding to the topographic inflection point. The outlier removal module is used to remove outlier data points with a water surface area greater than the threshold corresponding to the terrain inflection point, retain only valid data points within the valid data domain, and check whether the number of valid data points meets the minimum fitting sample size requirement. The high-precision model building module is used to select fitting functions from a preset candidate model library based on valid data points within the retained valid data domain, perform regression fitting using the least squares method, construct a high-precision area-water storage mathematical model of the water body, and determine the optimal function form, fitting parameters, and applicable domain of the high-precision area-water storage mathematical model. The batch automation processing and result output module is used to build corresponding area-water storage mathematical models for all water bodies to be measured in the target area one by one, so as to complete batch automation modeling; based on the latest water surface area obtained by remote sensing interpretation, call the corresponding high-precision area-water storage mathematical model to calculate the water storage of each water body in batches and generate a water storage update report.
[0030] To facilitate the understanding of the above technical solutions of the present invention, the model construction and water storage estimation based on turning point recognition in the actual process of the present invention will be described in detail below.
[0031] Step 1: Acquisition of underwater terrain data and construction of a three-dimensional model of the water area.
[0032] 1) Underwater terrain data collection: Use a single-beam sounding instrument to conduct underwater terrain surveys on all water bodies to be measured in the target area, including various types of water bodies such as lakes and reservoirs. During the measurement, preset measurement sections or survey line grids according to the morphological characteristics of each water body. For narrow water bodies, adopt the cross-section layout method, and the section spacing is set to 100m - 1000m according to the bending degree of the river; for planar water bodies, adopt the grid layout method, and the grid spacing is set to 50m - 500m. Along the preset section or grid nodes, use a single-beam sounding instrument to collect water surface elevation, water bottom elevation, and water depth data point by point, and synchronously record the plane coordinates of each measurement point. Finally, obtain the discrete underwater elevation point cloud data of each water body, and the data recording format of each measurement point is (X, Y, Z), where X and Y are plane coordinates, and Z is the elevation value of this point.
[0033] 2) Preprocess the collected original point cloud data, including the following sub-steps: Outlier removal: Based on the 3σ principle of the local neighborhood (i.e., the Ralida criterion), remove obvious outliers caused by measurement errors or environmental interference. Specifically, for each measurement point, calculate its average distance from K adjacent points in the neighborhood (K takes 5 - 10). If this average distance exceeds 3 times the standard deviation of the global average distance, it is determined as an outlier and removed. [[ID=q15]]
[0034] Data smoothing: For areas with slight noise, use the moving window weighted average method for smoothing, and this step can be carried out in the data processing software supporting the single-beam sounding equipment. The specific implementation method is as follows: ① Window Size Setting: A square moving window with a 5×5 grid size is used. The determination of this size is based on the following: (a) According to previous experiments, the 3×3 window is insufficient to suppress random noise, and the smoothed data still has local abnormal fluctuations, affecting the stability of subsequent curvature calculations; (b) Although 7×7 and larger windows can effectively eliminate noise, they can lead to excessive smoothing of underwater micro-topographic features (such as ditches and steep slopes), resulting in the loss of topographic information; (c) The 5×5 window achieves the best balance between improving the signal-to-noise ratio and preserving the fidelity of topographic features. For sparsely populated areas (grid coverage <70%), the window is adaptively expanded to 7×7; for areas with complex topography (slope >15°), the window is adaptively reduced to 3×3. ② Weighting Function Design: A distance-elevation dual-factor weighting strategy is adopted. The weight of each grid cell within the window is determined by the planar distance weight. w d Weight of elevation similarity w z A joint decision.
[0035] ; in, m , n The row and column indices of the grid cells within the window are relative to the center grid; the planar distance weights use a Gaussian kernel function. ; in, d mn The distance is the Euclidean distance from the grid cell to the center of the window. σ d Take 1 / 3 of the window radius (approximately 1.67 times the grid spacing); the elevation similarity weight is: ; in, z mn The elevation values are those of the neighboring grid. z c The original elevation value of the center grid. δ z Take 0.5 times the standard deviation of the elevation of the validation set for this water body. ③ Smoothing calculation: The elevation values of all grid points within the window are weighted and averaged according to the normalized weights. The calculation result is... The elevation value after smoothing the central grid: ; in p=2 corresponds to half the width of a 5×5 window. The above calculation is performed sequentially on all grid cells in the 3D model of the water body to obtain the smoothed 3D model. ④ Boundary processing: When the moving window exceeds the effective data boundary of the 3D model of the water body, only the effective grid data within the window is used for calculation. If the number of effective grids is less than 50% of the total number of grids in the window (i.e., less than 13 grids), the original elevation value of the central grid remains unchanged. ⑤ Iteration and quality control: The above smoothing process can be iterated 1-3 times. After the first smoothing, the signal-to-noise ratio improvement rate is calculated; if it is less than 5%, the iteration is terminated. Coordinate system one: The coordinate data of all water bodies are unified to the CGCS2000 National Geodetic Coordinate System and the 1985 National Height Datum.
[0036] 3) For each water body, independently construct a 3D model of the water area based on the preprocessed discrete point cloud data. The specific steps are as follows: Topographic feature point identification and densification: To compensate for the sparse spatial distribution of single-beam bathymetry data, it is necessary to identify and densify topographic feature points based on the original bathymetry points. Specifically, this includes: ① Feature point type definition: Define the following four types of terrain feature points: The thallophora is the point in a river or lake basin where the water is deepest, reflecting the lowest point of the underwater topography.
[0037] Change point: The point where the underwater slope changes significantly, usually located at the boundary between the gentle area of the lake basin and the steep area of the lake shore.
[0038] Bifurcation point: A key node where underwater topography branches or converges, often found in irregular lake or reservoir inlets.
[0039] Interpolation points: Points that are inserted at a certain interval between feature points to control the smooth transition of terrain morphology.
[0040] ② Feature point identification method: For thallopath points, compare the water depth values of adjacent measuring points along the original sounding profile line, and select the point corresponding to the local maximum water depth value. For slope change points, calculate the slope change rate between adjacent measuring points; when the slope change rate exceeds a preset threshold (default 15%), mark the location as a slope change point. For bifurcation points, identify the points where three or more measuring lines intersect based on the spatial topological relationship of the measuring lines. For interpolation points, insert them at equal intervals of 2m in the water surface direction between adjacent feature points.
[0041] ③ Elevation assignment: The plane coordinates of all identified or generated densified points are determined based on the linear relationship between the original measuring points or feature points, and the elevation values are obtained through linear interpolation or cubic spline interpolation along the measuring line direction.
[0042] Multi-source data fusion and irregular triangulation network (TIN) construction combine original bathymetry points, generated encrypted feature points, and boundary points (dike lines, shorelines) from the water spatial survey results to form a unified discrete point set. Based on this point set, the Delaunay triangulation algorithm is used to construct an irregular triangulated network (TIN). This step can be implemented in geographic information software, and the specific steps are as follows: ① Super Triangle Initialization: Construct a super triangle containing all discrete points as the initial mesh for subdivision. ② Point-by-Point Insertion and Empty Circumcircle Criterion: Insert each discrete point into the current triangular mesh sequentially. For each newly inserted point, find its corresponding triangle and split that triangle into three new triangles. Then, check the Delaunay condition (i.e., no other point is contained within the circumcircle of any triangle) edge by edge. If it is not satisfied, flip the edges (swap diagonals) until the entire mesh satisfies the empty circumcircle criterion. ③ Triangle Optimization: Perform Laplacian smoothing optimization on the generated triangular mesh, adjusting vertex positions to improve the uniformity of the triangles. The number of optimization iterations is set to 3-5.
[0043] Based on the natural neighborhood method, regular grid interpolation is used to interpolate the TIN (Transcript Indicator) into a regular grid 3D water model on the basis of the constructed irregular triangular mesh. The natural neighborhood method has the advantages of good conformity preservation and insensitivity to data distribution, and is suitable for non-uniformly sampled data such as underwater topography. The specific steps are as follows: ① Mesh definition: Set the target resolution of the 3D water model (default 2m×2m), generate a regular grid of points on the horizontal projection plane, and denote each grid point as . P ( x , y ) ② Finding natural neighborhoods: For each grid point P ( x , y Insert it into the constructed TIN, determine its Voronoi cell (i.e., Thiessen polygon cell), and find the set of TIN vertices adjacent to that Voronoi cell. Q 1, Q 2,..., Q k These vertices are} P ③ Weight calculation: For each natural neighbor point Q j ,calculate P and Q j Voronoi area weights between λ j : ; Among them, the molecule is PVoronoi unit and its natural neighborhood points Q j The intersection area of the Voronoi units, with the denominator being P The total area of the Voronoi cells. The sum of the weights of all natural neighboring points is 1.
[0044] ④ Elevation interpolation: grid points P ( x , y elevation value Z ( P ) Calculated by weighted average of the elevation values of its natural neighboring points: ; in, k Represents the current interpolation point P The number of natural neighborhood points, k It's not a fixed value, but depends on the point. P The surrounding spatial distribution is dynamically determined. j =1,2,3,..., k , indicating counting sequentially from the first natural neighbor point to the second. k A natural neighborhood point. Z ( Q j ) represents the first j Known elevation values of natural neighboring points, λ j Representing the j The weight coefficients of each natural neighbor point.
[0045] ⑤ Grid traversal: Repeat steps ② to ④ above for all regular grid points to generate a complete elevation grid.
[0046] 3D model generation and quality control of water bodies: ① Generation of 3D Water Body Models: The generated regular elevation grid, combined with the adopted water body spatial survey boundaries (dike lines, shorelines), is uniformly stored as a GeoTIFF format raster file, i.e., an independent 3D water body model for each water body. ② Quality Control: 10% of the original sounding points that were not used in the modeling are randomly selected as a validation set. The planar coordinates of the validation points are substituted into the generated 3D water body model to extract interpolated elevations, and the deviation from the measured elevations is calculated. The quality acceptance criteria are: root mean square error (RMSE) ≤ 0.1m and maximum absolute error ≤ 0.3m. If these criteria are not met, the resolution of the 3D water body model is adjusted or the feature point density is increased, and the modeling process is repeated.
[0047] Step 2: Extract the area-water storage relationship curve.
[0048] 1) Determine the water level sequence: Elevation range determination: Statistically analyze the elevation values of all raster cells in the 3D model of the water area and obtain the minimum value. Z min (Elevation of the lowest point underwater) and maximum value Z max (The elevation of the highest point at the water body boundary usually corresponds to the elevation of the dam crest or the natural shoreline elevation). Water level step setting: within the water level range [ Z min, Z max The water level sequence is set at equal intervals. The water level step size Δ h set Set to 0.01m; the step size can be adjusted according to different water bodies. Water level sequence generation: From... Z min Initially, increase the step size Δ sequentially. h until reached Z max Generate water level sequence { h 1, h 2, h 3,..., h n},in h 1= Z min , h n = Z max , n This refers to the number of water level grades.
[0049] 2) For each water level value in the water level sequence h i Calculate the water surface area below this water level. A ( h i The calculation method employs a flooding analysis algorithm, specifically as follows: Traverse all raster cells in the 3D model of the water area; for each raster cell, compare its elevation values. Z grid Compared with the current water level h i :like Z grid ≤ h i If the grid cell is submerged in water, it is marked as a water area; if Z grid > h i If a grid cell is not submerged, it is marked as land. Area statistics: Count the number of all grid cells marked as water. N wet Multiply by the actual area represented by each grid cell A cellGet the current water level h i The area of the water surface below: ; 3) For each water level value in the water level sequence h i Calculate the water volume at this water level. V ( h i The calculation method uses the volume integral method, as detailed below: Single-cell volume calculation: For each submerged cell, calculate the volume of the water column within that cell. If the elevation of the cell is... Z grid The current water level is hi ( h i ≥ Z grid Then the height of the water column in the grid cell is ( h i Subtract Z grid The volume of the water column is: ; The current water level is obtained by summing the volumes of the water columns in all submerged grid cells. h i Total water volume below: ; in, i The index number of the submerged raster cell. For the first i Elevation values of individual grid cells A cell This represents the area of a single grid cell.
[0050] 4) Pair the calculated water surface area and water volume corresponding to each water level to generate a dataset of area-water storage relationship for each water body, and a lookup table of water surface elevation-water depth-area-water storage. Plot a scatter plot with water surface area A as the abscissa and water storage V as the ordinate, and connect them to form a curve to obtain the area-water storage relationship curve for the water body.
[0051] 5) After generating the area-water storage relationship dataset, to assess data quality and provide a reference for constructing subsequent mathematical models (i.e., high-precision area-water storage mathematical models), pre-analysis of the curves can be performed. This includes: using multiple regression methods to pre-fit the generated area-water storage relationship dataset to initially assess the data's fit and the optimal fitting function form. The fitting methods used include, but are not limited to, the following 11: Linear functions: Quadratic polynomial: cubic polynomial: ; Quadratic polynomial: Fifth-degree polynomial: Exponential function: Logarithmic function: Power function: Simple Gaussian function: Bigaussian function: ; Quadratic Gaussian function: .
[0052] Parameters in the above formulas a , b , c , d , e , f All are regression coefficients to be fitted, determined using the least squares method based on valid data points. Where: the parameters in the polynomial function represent the coefficients and constants of each power term; the parameters in the exponential function... a This is the scaling factor. b The growth rate coefficient; in the logarithmic function a The logarithmic gain coefficient, b The offset constant; in the power function a This is the scaling factor. b The exponent is the power; in the Gaussian function a , d , g Peak amplitude, b , e , h This is the peak position. c , f , i The standard deviation is given. This invention is not limited to the 11 function forms mentioned above; other suitable functions can be selected for fitting based on data characteristics.
[0053] Optimal fitting method selection: For the above 11 fitting methods, their coefficients of determination were calculated respectively. R 2 The calculation formula is: ; in, V i For the first i The actual volume value of each data point. This represents the predicted volume value of the current candidate mathematical model. This is the average value of the actual volume. n This represents the number of data points. R 2The value of is between 0 and 1; the closer it is to 1, the better the fit and the stronger the interpretability of the current candidate mathematical model on the data. Compare the various methods. R 2 Value, will R 2 The method with the largest value is marked as the candidate optimal fitting method for that water body, for reference in subsequent modeling. The function form and fitting parameters of this optimal method are also recorded. The generated area-water storage relationship dataset is stored in a structured format, with each record containing: water surface elevation, water depth, water surface area, and water storage. Additionally, the following metadata information is recorded: water body name, water body code, 3D model resolution of the water body, water level step size, number of data points, and whether the curve exhibits area saturation characteristics, for later reference.
[0054] Step 3: Automatically identify the boundaries of valid data domains (terrain inflection points).
[0055] Based on the generated area-water storage relationship curve, the effective data domain boundary points (i.e., inflection points) on the curve are automatically identified. These inflection points characterize the critical positions where the underwater topography of the water body transitions from a freely expanding zone to a boundary-constrained zone: before the inflection point (where the area is smaller than the area corresponding to the inflection point), there is a strict monotonic function mapping relationship between area and volume, which belongs to the effective data domain; after the inflection point (where the area is larger than the area corresponding to the inflection point), due to the constraints of dikes or natural shorelines, there is an anomaly of one-to-many mapping where the area tends to saturate while the volume continues to grow, which belongs to the invalid data domain.
[0056] This invention provides three independent inflection point identification strategies, which can be selected based on data quality (i.e., the data quality of the water body area-water storage relationship dataset) and water body morphology characteristics, or combined to improve robustness. The specific implementation methods of the three strategies are described below: Strategy A: Slope Variation Method. The slope variation method identifies inflection points based on the first derivative characteristics of the area-water storage curve. Its core idea is: within the effective data domain, the area monotonically increases with increasing volume, but the slope gradually increases; when entering the boundary constraint region, the area tends to saturate, and the slope drops sharply. The inflection point is the critical point where the slope changes from an upward trend to a sharp decline. The specific steps are as follows: Calculate the first derivative of area with respect to volume for discrete data points on the area-water storage curve ( A i , V i () i =1,2,..., n According to volume V (Sorted in ascending order), the first derivative of the area with respect to the volume at each point is calculated using the central difference method. k i : ; in, dA Representing area A The differential increment, dV Representing volume V The differential increment, dA / dV That is, the first derivative of area with respect to volume. i= 1,2,..., n The index number of the data point. n This represents the total number of data points. For endpoints ( i= 1 and i=n The derivative is calculated using either forward or backward differencing. k i The physical meaning of is the change in water surface area caused by a unit volume change, reflecting the sensitivity of the area to volume growth.
[0057] Identify the point of maximum slope, iterate through all data points, and find the first derivative. k i Take the data point with the maximum value P max Its corresponding derivative value is denoted as k max This point corresponds to the region where the area increases most rapidly with volume, typically located on steep slopes within a lake basin. From the point of maximum slope... P max Begin by searching point by point along the direction of increasing volume (backwards) until the first data point that satisfies the following conditions. P turn : ; in, α The slope descent threshold is used to control the sensitivity of inflection point identification. α The value range is 0.005-0.05, and the default value is 0.01 (i.e. 1% of the maximum slope).
[0058] Identify the turning point and use the points obtained from the search. P turn The corresponding area A turn and volume V turn This point is identified as the inflection point of the water body. It signifies that area growth has slowed significantly, and the water body has begun to enter the boundary constraint zone. The advantage of the slope change method is its simplicity and intuitiveness, making it suitable for regular water bodies with a clear unimodal slope characteristic in their area-water storage curves. However, its disadvantage is that the slope change may lead to misjudgments when the curve has local fluctuations or noise.
[0059] Strategy B: Cumulative Growth Rate Method. The cumulative growth rate method identifies inflection points based on the cumulative growth rate of area with increasing volume. Its core idea is: within the effective data domain, the area continuously increases with increasing volume, and the cumulative growth rate steadily rises; when entering the boundary constraint region, the area tends to saturate, and the cumulative growth rate approaches 1 (i.e., 100%). The inflection point is the position where the cumulative growth rate first reaches a preset threshold. This strategy is insensitive to data noise and has good robustness, making it the preferred strategy of this invention. The specific steps are as follows: Determine the maximum area value; find the maximum water surface area from the area-water storage relationship dataset. A max This typically corresponds to the water surface area at the highest water level (i.e., the crest elevation of the dam or the maximum inundation area). The cumulative growth rate is calculated for each data point ( A i , V i () i =1,2,..., n ), calculate its cumulative growth rate R i : ; in, R i The value range is 0-1, representing the proportion of the current area to the maximum area.
[0060] By volume V Traverse all data points in ascending order and find the first data point that satisfies the following condition. P turn : ,in, λ The cumulative growth rate threshold is set between 0.90 and 0.98, with a preferred value of 0.95.
[0061] Determine the turning point, and place the point P turn The corresponding area A turn and volume V turn This point was identified as the turning point for the water body. This point signifies that the area growth has reached its maximum. λ After this point, the area growth potential becomes extremely limited, and the water body enters the boundary constraint zone. Parameters λ The optimization determination of the above preferred values λ =0.95 was determined through the following experimental optimization process: representative water bodies within the study area (such as 248 lakes in a certain region) were selected, and candidate thresholds were set. λ∈{0.80,0.90,0.93,0.95,0.98,1.00}. For each candidate threshold, identify the inflection point and perform mathematical model construction. Calculate the fitting error of all water bodies under each threshold and the percentage of water bodies with an error less than 10%. Experimental results show that when... λ At a value of 0.95, the proportion of water bodies with an error of less than 10% reaches a peak of 98%, the average relative error is the lowest (3.24%), and there are no outlier error values, representing the optimal balance between accuracy and data utilization. Experimental results: Cumulative growth rate λ Parameter optimization experiments.
[0062] Measured underwater topographic data from 248 lakes in a certain region were selected. Seven candidate thresholds were set. λ ∈{0.80,0.90,0.93,0.95,0.98,1.00}. Where... λ =1.00 is equivalent to using the full dataset (without removing outliers) as a control group. For each candidate... λ Values are processed according to the cumulative growth rate method; different values... λ The batch fitting results for the given values are summarized in Table 1: Table 1 Cumulative Growth Rate λ Optimization results value MRE RE≤10% quantity RE≤10% Exception RE 0.80 5.89% 238 95.9% 5 0.90 3.0% 241 97.2% 1 0.93 3.28% 242 97.6% 1 0.95 3.24% 243 97.9% 0 0.98 5.61% 218 87.9% 0 1 20.98% 16 6.5% 0 In conclusion, λ A value of 0.95 strikes the best balance between fitting accuracy and data utilization, and is therefore determined as the preferred default parameter in this invention. Meanwhile, considering the differences in water body morphology, λ It can be adjusted within the range of 0.90-0.98.
[0063] Strategy C: Second Derivative Method. The second derivative method identifies inflection points based on the curvature change characteristics of the area-water storage curve. Its core idea is that the physical shape of the area-water storage curve indicates that the curve is concave within the effective data domain (second derivative is negative) and convex within the boundary constraint region (second derivative is positive). The inflection point is the point where the curve changes from concave to convex. The specific steps are as follows: Calculate the second derivative of area with respect to volume, based on the calculated first derivative. k i Further calculate the second derivative : ; Identify the zero-crossing points where the second derivative changes from negative to positive, and iterate through all data points to find the second derivative. The zero-crossing point where a negative value changes to a positive value. Since the data points are discrete values, the specific criterion is: there exist two adjacent points. i and i +1, satisfied. <0 and >0. The turning point is located at... i and i Between +1, linear interpolation can be used for precise positioning: ; ; Determine the inflection point, and interpolate the ( A turn , V turn This was identified as the turning point of the water body.
[0064] The advantages of the second derivative method are its theoretical rigor and ability to accurately locate inflection points of curves; its disadvantages are its reliance on the overall smoothness of the curve, lack of measurement noise, and the potential impact of measurement errors in actual data on the stability of the second derivative calculation. Regarding strategy selection and combination rules, the three strategies mentioned above can be used independently or in combination. This invention provides the following selection rules based on practical application scenarios: Single strategy selection: If the water body has a regular shape and high data quality, any strategy can be chosen; if the water body is a lake, reservoir, or other ise-shaped water body, strategy B (cumulative growth rate method) is preferred because it has the strongest robustness and highest accuracy; if the water body is a river or other linear water body, strategy A (slope change method) is preferred. Multi-strategy combination: To improve the stability and reliability of inflection point identification, a multi-strategy combination approach is recommended: simultaneously use strategy A, strategy B, and strategy C to identify inflection points, obtaining three candidate inflection points. A turn ( A ), V turn ( A )), ( A turn ( B ), V turn ( B )), ( A turn ( C ), V turn ( C If the three identification results are similar (area difference less than 10%), the average of the three is taken as the final inflection point; if there is a significant difference, the one with the smallest area among the three is taken as the final inflection point, to ensure that all potential outlier data points are eliminated and to prioritize the effectiveness of the fitted data domain. In a specific application example of this invention, for the batch processing of 724 water bodies in a certain area, strategy B (cumulative growth rate method) is adopted. λ With 0.95 as the default configuration, a 98% accuracy rate was achieved.
[0065] Step 4: Based on the identified inflection points, perform outlier data domain removal processing on the generated area-water storage relationship dataset.
[0066] 1) Definition and identification principle of abnormal data domain: The area-water storage relationship curve of a water body can be physically divided into three segments, as shown in Table 2: Table 2. Segmented curves showing the relationship between water body area and water storage. Section scope Physical characteristics Data properties Section I: Free Expansion Zone Area <Aturn The surface area of water increases continuously as the water level rises, and there is a strict monotonically increasing functional relationship between area and volume. Valid data field Section II: Transition Zone Area ≈ Aturn The growth of water surface area has begun to be constrained by boundaries, and the growth rate has slowed down, but it still maintains a basically monotonic relationship. Critical region Section III: Boundary Constraint Area Area > Aturn Water surface area tends to saturate due to the restriction of dikes or natural shorelines, resulting in multiple volume values for the same area, thus exhibiting a one-to-many mapping anomaly. Abnormal data fields (to be removed) 2) Outlier identification and removal: From the generated area-water storage relationship dataset, the data nature of each point is determined. Each data point in the dataset is recorded as ( A i ,V i ),in i =1,2,..., n According to volume V Sort in ascending order. The criteria are as follows: ; Based on the judgment results, perform the following removal operation: Removal execution: Remove all data that meet the criteria from the dataset. A i > A turn Data points. Removed records: Records the number of data points removed. n reject and the total number of original data points n total Calculate the rejection ratio: ; Data volume check after removal: Check the number of valid data points remaining after removal. n valid = n total - n reject : like n valid If the value is ≥5, the minimum fitted sample size requirement is met, and proceed to step five. n valid If the score is less than 5, a warning is issued, indicating that there are insufficient valid data points for the water body. In this case, the full data (or the threshold for the inflection point is relaxed) is used for fitting, and the mathematical model of the water body is marked as low confidence in the final result.
[0067] 3) After the removal operation is completed, the following information is output for building a high-precision area-water storage mathematical model: Valid dataset: The retained area is less than or equal to A turn Data point set {( Ai ,V i )∣ A i ≤ A turn Turning point information: A turn (Effective area threshold) V turn (Volume corresponding to the inflection point). Elimination statistics: number of eliminated data points, elimination ratio, number of valid data points. Modeling domain: area range of the valid dataset. A min , A turn ],in A min This is the minimum area corresponding to the lowest water level (usually close to 0). This domain will serve as the applicable boundary for the mathematical model constructed in step S5.
[0068] Step 5: Construct a high-precision area-water storage mathematical model.
[0069] 1) Model fitting and validation, based on the effective dataset obtained after removing outlier data domains {( A i , V i )∣ A i ≤ A turn This step constructs a high-precision area-water storage mathematical model for the water body. The modeling process references the process of generating the area-water storage relationship curve and the selection process for fitting the area-water storage relationship. It should be noted that the mathematical model in this invention (i.e., the high-precision area-water storage mathematical model) essentially refers to an area-volume regression fitting model (such as a power function, polynomial function, etc.), rather than a neural network model in deep learning, and there is no model training or optimization process. Candidate mathematical model library definition: This invention pre-defines a candidate mathematical model library containing 11 function forms. This model library covers various forms such as linear, polynomial, exponential, logarithmic, power functions, and Gaussian functions to adapt to different water body morphological characteristics. The mathematical model data source comes from valid datasets.
[0070] Mathematical model parameter estimation: For each candidate model, the least squares method is used to estimate the parameters of the effective dataset. The core objective of the least squares method is to find a set of model parameters that minimizes the sum of squared residuals between the mathematical model's predictions and the actual observations. This is achieved using a power function. V = a ⋅ A b For example, the parameter estimation process is as follows: Linearization transformation: Taking the logarithm of both sides of the original function transforms it into a linear form. ; Constructing a system of equations: To facilitate solving for the parameters of the power function using the least squares method, we take the natural logarithm of both sides of the original function and perform a linearization transformation: ; In the formula, Y This represents the transformed dependent variable, which is dimensionless. Y =ln( V ),in V Water storage capacity (unit: m³); X The independent variable after transformation is dimensionless. X =ln( A ),in A Water surface area (unit: m²); β 0 represents the transformed intercept, which is dimensionless. β 0 = ln( a ),in a The proportionality coefficient of the power function; β 1 represents the transformed slope, which is dimensionless. β 1= b ,in b is the exponent of the power function.
[0071] Solving parameters: based on valid data points ( X i , Y i The solution is obtained using the least squares method. ; in, , Represent , The sample mean. Restoration parameters: a = eβ 0, b = β 1.
[0072] For polynomial functions, the least squares method is used to solve the normal equations to obtain parameter estimates. For exponential and logarithmic functions, a similar linearization transformation is used for parameter estimation. The parameter estimation process for other functions (linear, polynomial, etc.) is similar to that for power functions, all using the least squares method. For optimal mathematical model selection, the coefficient of determination is calculated for each fitting method. R 2 ),Will R 2The largest mathematical model is taken as the optimal mathematical model. Mathematical model validation includes: validation point independent variable: inflection point area. Validation point predicted value: the water storage calculated by substituting the inflection point area into the optimal mathematical model. Validation point true value: the water storage corresponding to the inflection point area in the original area-water storage lookup table. Validation requirement: if the relative error between the validation point predicted value and the validation point true value does not exceed 10% (preset value), it meets the requirement. Validation results: Setting the water level (depth) step size to 0.01m, the water level gradually rises from the lowest to the highest level. The inundation analysis algorithm is used to calculate the water surface area and water storage corresponding to each water level value, generating an area-water storage relationship lookup table. Taking Xiaonan Lake and Shanzi Lake in a certain region as examples, traditional methods use all data points for fitting. Due to the inclusion of outliers constrained by boundaries, the mathematical model exhibits systematic bias in shallow water areas. The R² of Xiaonan Lake is only 0.6568, with a relative error (RE) of 16.13%. x ²、 x ³、 x 4 , x 5 These represent quadratic, cubic, quartic, and quintic terms, respectively. For example: x ² refers to x The square of the slope change method and the cumulative growth rate method are both effective strategies for identifying inflection points, effectively eliminating outlier data domains and showing significantly better fitting results than traditional methods. Among them, the cumulative growth rate method performs better. R The value reached 0.9991, and the average relative error decreased to 0.42%.
[0073] It needs to be explained that the fitting equations for the three strategies (Xiaonanhu) are as follows: Traditional method: RE (%) = 16.13; R 2 =0.9568; the fitted equation is .
[0074] Slope change method: RE (%) = 8.14; R 2 =0.9844; the fitted equation is .
[0075] Cumulative growth rate method: RE (%) = 0.42; R 2 =0.9991; the fitted equation is .
[0076] Three strategies for fitting the equation (Shanzihu): Traditional method: RE (%) = 23.37; R 2=0.8732; the fitted equation is .
[0077] Slope change method: RE (%) = 5.07; R 2 =0.9563; the fitted equation is .
[0078] Cumulative growth rate method: RE (%) = 2.28; R 2 =0.9969; the fitted equation is .
[0079] Figure 2 and Figure 3 The fitting curves of the three methods are compared, and the inflection point locations are marked. Figure 2 The subplots represent: (a) the curve fitted using the traditional method for Xiaonan Lake, (b) the curve fitted using the slope change method for Xiaonan Lake, (c) the curve fitted using the cumulative growth rate method for Xiaonan Lake, (d) the curve fitted using the traditional method for Shanzi Lake, (e) the curve fitted using the slope change method for Shanzi Lake, and (f) the curve fitted using the cumulative growth rate method for Shanzi Lake. Figure 3 The subplots represent: (a) finding the inflection point using the slope change method of Xiaonan Lake, (b) finding the inflection point using the cumulative growth rate method of Xiaonan Lake, (c) finding the inflection point using the slope change method of Shanzi Lake, and (d) finding the inflection point using the cumulative growth rate method of Shanzi Lake. The second derivative method is not applicable to these two cases. Reservoir simulation data is used to demonstrate the comparison of the fitted curves and the location of the inflection points using the second derivative method, as shown below. Figure 4 As shown. Figure 4 The subplots represent: (a) the location of the area inflection point, (b) the peak value of the second derivative, and (c) the area-volume relationship curve.
[0080] Based on the above modeling process, strategy 0 (traditional method) and strategy B (cumulative growth rate method, λ=0.95) were implemented in batches for 248 lakes in a certain region. The relative error (RE) of the fitted mathematical model for each lake was calculated, and the number of lakes was counted according to the error interval. The results are shown in Table 3. Table 3. Accuracy verification of the mathematical model for 248 lakes in a certain region. Fitting methods RE≤10% quantity RE>10% quantity Average RE Median RE Traditional methods 16 232 20.98% 20.62% Cumulative growth rate method 243 5 3.24% 3.24% The above results show that, for 248 lakes in a certain region, only 16 lakes meet the accuracy requirement of less than 10% error using the traditional method, while Strategy B of the present invention increases the compliance rate to about 98%, verifying the universality and high accuracy of the method of the present invention.
[0081] 2) Mathematical Model Storage and Domain Determination: The completed mathematical model is stored in a standardized format. Output includes: Water body identification information: water body name, water body code, and water body type. Optimal mathematical model: fitting function form, fitting parameter values, and coefficient of determination. R 2 Domain of the mathematical model: Lower limit of area A min Area limit A turn .
[0082] The above outputs will be summarized in the results files "Mathematical Model Information Table.xlsx" and "Fitting Equations and Parameters.txt" for subsequent water storage estimation and annual updates. The domain is determined, and the constructed mathematical model... V = f ( A It only has high-precision predictive ability within the valid data domain. Therefore, it is necessary to clearly define the applicable boundaries (domain) of the mathematical model: ; in, A min This is the lower limit of the area. A turn This is the upper limit of the area. When the water surface area to be estimated exceeds the above definition, the prediction results of the mathematical model may have a large deviation, which needs to be marked in the output results.
[0083] Step 6: Batch output of mathematical models and application of mathematical models.
[0084] 1) Batch processing workflow: Process all data within the target area. M For each water body, steps one through five are executed sequentially to construct a mathematical model of area-water storage for each body. The entire batch processing process requires no manual intervention; all parameters (such as water level step size, slope descent threshold, cumulative growth rate threshold, etc.) are determined using default configurations or adaptively, achieving fully automated modeling.
[0085] The batch processing using Python software employs a loop iteration method, and the specific process is as follows: For Water body number = 1 to M : Step 1: Obtain a 3D model of the water body.
[0086] Step 2: Extract the area-water storage relationship curve.
[0087] Step 3: Identify the turning point (valid data domain boundary).
[0088] Step 4: Remove abnormal data fields.
[0089] Step 5: Construct a mathematical model and select the best fitting method.
[0090] Store the mathematical model results of the water body in a temporary cache.
[0091] End For (This represents the end marker of a loop structure).
[0092] 2) Output of results files: After batch processing is completed, the following four results files will be automatically generated to form a complete mathematical model library of water bodies in a certain region: Fitting curve atlas (fitting curve plot.jpg); Mathematical model information table (mathematical model information table.xlsx); Fitting equation and parameter record (fitting equation and parameters of multiple fitting methods and the optimal fitting method.txt); Batch fitting result error analysis table (batch fitting result error analysis table.xlsx).
[0093] 3) Mathematical Model Application: The completed area-water storage mathematical model is used to support the annual dynamic update of water storage in a certain region. Specific application methods include: acquiring water surface area data; periodically (e.g., annually or quarterly) acquiring the latest water surface area data for all water bodies within the target area. The water surface area data primarily comes from remote sensing satellite image interpretation, and the normalized water index (NDWI) is used to extract the water surface area, obtaining the latest water surface area for each water body. A current The specific method is: Based on high-resolution optical remote sensing images (such as Gaofen series, Ziyuan series, Sentinel-2, etc.), after radiometric calibration and atmospheric correction, the Normalized Water Index (NDWI) is calculated. The threshold of NDWI (the threshold of the Normalized Water Index) is determined using the threshold method. When the threshold is ≥0.2 (preset value), it is identified as the boundary of all water bodies in the study area. The vector layers of reservoirs and lakes in the latest annual land change survey data are collected. The NDWI extracted boundaries and the vector boundaries of the annual land change survey are input into the geographic information software. After correlation, the water surface boundaries of lakes and reservoirs can be accurately extracted.
[0094] ; In the formula: Indicates spectral reflectance; Green Indicates the green light band; NIR Indicates the near-infrared band.
[0095] Bulk estimation of water storage capacity, for each water body, based on its latest water surface area. A current The system uses the constructed mathematical model and water storage calculation formula to calculate the current water storage volume. ; Based on the calculation results, an annual update report on water storage in a certain region is generated, including the total water storage volume and annual changes, as well as detailed changes in the storage volume of each water body. It should be noted that the method provided by this invention can be packaged into a software system and deployed on the server of a water resources management department to achieve the following functions: Data import interface: Supports batch import of underwater measured data and remote sensing water surface area data.
[0096] Automatic modeling engine: Calls steps one through five to achieve automated modeling of water bodies in a certain area.
[0097] Storage estimation module: Periodically obtain the latest water surface area and calculate water storage in batches.
[0098] Visualization: Displaying the distribution and changing trends of water storage in a certain area in the form of a map.
[0099] Report Export: Generate annual update reports and related deliverables with one click.
[0100] Specifically, this invention retains the flexibility of parameter adjustment: slope descent threshold α It can be adjusted within the range of 0.005-0.05; cumulative growth rate threshold. λ The value can be adjusted within the range of 0.90-0.98; users can fine-tune the parameters according to the water morphology characteristics of different regions (such as mountain reservoirs and plain lakes).
[0101] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.
Claims
1. A method for model construction and water storage estimation based on inflection point identification, characterized in that, include: Acquire discrete underwater elevation point cloud data of all water bodies to be measured within the target area, and preprocess the discrete underwater elevation point cloud data; Based on the preprocessing results, construct three-dimensional models of the water bodies for each water body; Based on the three-dimensional model of the water body, the water level elevation range and water level step size are determined, and an equally spaced water level sequence is generated. Based on the equally spaced water level sequence, the water surface area and water storage corresponding to each water level are calculated respectively, and the water surface area and water storage corresponding to each water level are paired to generate a water body area-water storage relationship dataset and a water body area-water storage relationship curve. Based on the water body area-water storage relationship dataset, a preset inflection point identification strategy is adopted to automatically identify the water body area-water storage relationship curve, determine the topographic inflection point on the area-water storage relationship curve that represents the boundary of the effective data domain, and obtain the water surface area and water storage corresponding to the topographic inflection point. Using the water surface area corresponding to the topographic inflection point as the threshold, outlier data points with water surface areas greater than the threshold are removed, and only valid data points within the valid data domain are retained. The number of valid data points is then checked to see if it meets the minimum fitting sample size requirement. Based on the valid data points within the retained valid data domain, a fitting function is selected from the preset candidate model library, and the least squares method is used for regression fitting to construct a high-precision area-water storage mathematical model of the water body. The optimal function form, fitting parameters and applicable domain of the high-precision area-water storage mathematical model are determined. For each water body to be measured within the target area, a corresponding area-water storage mathematical model is constructed to complete batch automated modeling. Based on the latest water surface area obtained from remote sensing interpretation, the corresponding high-precision area-water storage mathematical model is called to calculate the water storage of each water body in batches and generate a water storage update report.
2. The method for model construction and water storage estimation based on inflection point identification according to claim 1, characterized in that, The process involves acquiring discrete underwater elevation point cloud data of all water bodies to be measured within the target area and preprocessing the discrete underwater elevation point cloud data. Based on the preprocessing results, the three-dimensional models of each water body are constructed as follows: A single-beam echo sounder is used to conduct underwater topographic surveys of lakes and reservoirs. Measurement sections or survey line grids are laid out according to the water body morphology. Water surface elevation, bottom elevation, water depth and plane coordinates are collected to form discrete underwater elevation point cloud data containing plane coordinates and elevation values. Outliers in point cloud data are removed based on the local neighborhood Laida criterion. The moving window weighted average method is used to smooth areas with slight noise. The coordinates of all water bodies are unified to the preset geodetic coordinate system and preset elevation datum to obtain the preprocessed discrete point cloud data. The preprocessed discrete point cloud data is used to identify and encrypt terrain feature points, generate encrypted feature points, and integrate the original underwater elevation measurement points, encrypted feature points, and boundary points in the water spatial survey results to form a unified discrete point set. Based on the discrete point set, an irregular triangular network is constructed using a triangulation algorithm, and the irregular triangular network is interpolated into a regular grid water 3D model using the natural neighborhood method.
3. The method for model construction and water storage estimation based on inflection point identification according to claim 2, characterized in that, The construction of irregular triangular meshes based on discrete point sets and using triangulation algorithms includes: Construct a super triangle containing all discrete points as the initial mesh for partitioning; Insert each discrete point into the current triangulation in turn, find the triangle where the newly inserted point is located, and split the triangle into three new triangles. Check the triangulation condition edge by edge, and flip the edges that do not meet the empty circumcircle criterion until the triangulation meets the empty circumcircle criterion. The generated triangular mesh is optimized using Laplacian smoothing, and the vertex positions are adjusted to improve the uniformity of the triangles. The triangular mesh optimization is completed according to the set number of iterations to obtain an irregular triangular mesh.
4. The method for model construction and water storage estimation based on inflection point identification according to claim 2, characterized in that, The method of interpolating irregular triangular meshes into regular meshes to create a three-dimensional water area model includes: Set the target resolution of the 3D water model and generate a regular grid of points on the horizontal projection surface; Each grid point is inserted into the irregular triangular mesh, the Thiessen polygon cell in which the grid point is located is determined, the set of irregular triangular mesh vertices adjacent to the Thiessen polygon cell is obtained by filtering, and the irregular triangular mesh vertices are taken as natural neighborhood points. Calculate the Thiessen polygon area weights between grid points and their natural neighbors, and then calculate the weighted average of the elevation values of the natural neighbors based on these weights to obtain the interpolated elevation of the current grid point. Traverse all regular grid points, and sequentially complete natural neighborhood filtering, weight calculation and elevation interpolation to generate a complete regular elevation grid; The generated regular elevation grid is combined with the boundary points in the water spatial survey results and stored in a unified raster file format to obtain the three-dimensional water body model of each water body.
5. The method for model construction and water storage estimation based on inflection point identification according to claim 1, characterized in that, The method involves determining the water level elevation range and water level step size based on a three-dimensional water body model, and generating an equally spaced water level sequence. Based on the equally spaced water level sequence, the water surface area and water storage corresponding to each water level are calculated, and the water surface area and water storage corresponding to each water level are paired to generate a water body area-water storage relationship dataset and a water body area-water storage relationship curve, including: The elevation values of all grid cells in the three-dimensional model of the water body are statistically analyzed to obtain the minimum and maximum elevation values, determine the range of water level elevation, and set a fixed water level step within this range. The water level step can be adaptively adjusted according to the actual conditions of different water bodies. Starting from the minimum elevation, the water level is increased sequentially according to the set water level step size until the maximum elevation is reached, generating an equally spaced water level sequence. For each water level value in the water level sequence, the flooding analysis algorithm is used to traverse all grid cells of the three-dimensional water area model, compare the relationship between the grid cell elevation value and the current water level value, and based on the comparison results, the grid cells are marked as water area and land area respectively. The number of grid cells marked as water area is counted, and multiplied by the actual area represented by each grid cell marked as water area to calculate the water surface area of the corresponding water level. For each water level value in the water level sequence, the volume integral method is used to calculate the water column volume in each grid cell marked as a water area, and the water column volumes of all grid cells marked as water areas are summed to obtain the water storage at the corresponding water level. The water surface area and water storage corresponding to each water level are associated and paired one by one to construct a dataset of water body area-water storage relationship. At the same time, a lookup table of water surface elevation, water depth, water surface area and water storage relationship is generated. Scatter points are plotted with water surface area as the horizontal axis and water storage as the vertical axis and connected into a smooth curve to obtain the water body area-water storage relationship curve.
6. The method for model construction and water storage estimation based on inflection point identification according to claim 1, characterized in that, The data set based on the area-water storage relationship of water bodies employs a preset inflection point identification strategy to automatically identify the area-water storage relationship curve of water bodies, determine the topographic inflection points on the area-water storage relationship curve that represent the boundaries of the effective data domain, and obtain the water surface area and water storage corresponding to the topographic inflection points, including: Based on the data quality and water morphology characteristics of the water body area-water storage relationship dataset, a preset inflection point identification strategy is selected to automatically identify the water body area-water storage relationship curve. The inflection point identification strategies include the slope change method, the cumulative growth rate method, and the second derivative method. When identifying inflection points using the slope change method, the discrete data points on the area-water storage relationship curve are sorted in ascending order of volume. The central difference method combined with forward and backward difference methods is used to calculate the first derivative of the area with respect to the volume of each data point. All data points are traversed to determine the point corresponding to the maximum value of the first derivative. Then, the points are searched point by point along the direction of increasing volume. The target inflection point is selected by combining the slope descent threshold, and the water surface area and water storage corresponding to the inflection point are extracted. When using the cumulative growth rate method to identify inflection points, the maximum water surface area is determined from the area-water storage relationship dataset, and then the cumulative growth rate of the area of each data point relative to the maximum water surface area is calculated point by point. The data points are traversed in ascending order of volume, and the target inflection point is selected by combining the preset cumulative growth rate threshold, and the water surface area and water storage corresponding to the inflection point are extracted. When using the second derivative method to identify inflection points, the second derivative of each data point is calculated based on the area-water storage relationship curve and the calculated first derivative. The discrete data points are traversed, and the adjacent points where the second derivative changes from negative to positive are retrieved. The inflection point of the curve is located by linear interpolation and used as the target inflection point. The water surface area and water storage corresponding to the inflection point are then extracted.
7. The method for model construction and water storage estimation based on inflection point identification according to claim 1, characterized in that, The process of using the water surface area corresponding to the terrain inflection point as a threshold to remove outlier data points with a water surface area greater than the threshold, retaining only valid data points within the valid data domain, and checking whether the number of valid data points meets the minimum fitting sample size requirement includes: Based on the segmented characteristics of the area-water storage relationship curve, three segments are defined: the free expansion zone, the transition zone, and the boundary constraint zone. The interval with an area smaller than the water surface area corresponding to the topographic inflection point is defined as an effective data domain with a strictly monotonically increasing function relationship, while the interval with an area larger than the water surface area corresponding to the topographic inflection point is defined as an abnormal data domain constrained by dikes or natural shorelines and with one-to-many mapping anomalies. The data points in the area-water storage relationship dataset are sorted in ascending order of volume. The water surface area corresponding to the topographic turning point is used as the discrimination threshold. The water surface area of each data point is compared with the threshold to identify abnormal data points with a water surface area greater than the threshold. Remove all outlier data points from the area-water storage relationship dataset where the water surface area is greater than the water surface area corresponding to the topographic inflection point. Only retain valid data points where the water surface area is less than or equal to the threshold and is within the valid data domain. Simultaneously count the number of removed data points, the total number of original data points, and calculate the data removal ratio. The number of valid data points remaining after removing outlier data points is counted. It is then verified whether the number of valid data points meets the preset minimum fitting sample size requirement. If the number of valid data points meets the minimum fitting sample size requirement, the modeling process continues. If the number of valid data points does not meet the minimum fitting sample size requirement, an insufficient data warning is issued, and the model is reverted to using the full dataset or a relaxed inflection point threshold for fitting. The high-precision area-water storage mathematical model of the water body is then marked as having low confidence. After removing outlier data points, the system outputs the remaining valid dataset, the area threshold and volume information corresponding to the turning point, the data removal statistics, and the modeling domain corresponding to the valid dataset. The modeling domain is then used as the applicable boundary for constructing a high-precision area-water storage mathematical model of the water body.
8. The method for model construction and water storage estimation based on inflection point identification according to claim 1, characterized in that, Based on the valid data points within the retained effective data domain, a fitting function is selected from a preset candidate model library, and regression fitting is performed using the least squares method to construct a high-precision area-water storage mathematical model of the water body. The optimal function form, fitting parameters, and applicable domain of this high-precision area-water storage mathematical model are determined as follows: A predefined library of candidate models containing various function forms is provided, which is adapted to different water body morphological characteristics. For each fitting function in the candidate model library, the least squares method is used to estimate the parameters based on the effective dataset. The goal is to minimize the sum of squared residuals between the predicted values of the high-precision area-water storage mathematical model of the water body and the actual observed values. The linear transformation of the power function, exponential function, and logarithmic function is performed to reconstruct the system of equations to solve for the parameters. For polynomial functions, the parameter estimation is completed by solving the normal system of equations. The coefficient of determination is calculated for the fitting results of each candidate fitting function. The fitting function with the largest coefficient of determination is selected as the high-precision area-water storage mathematical model of the water body as the optimal high-precision area-water storage mathematical model, and the corresponding fitting parameters are determined. Using the water surface area corresponding to the topographic turning point as the verification independent variable, the verification independent variable is substituted into the optimal high-precision area-water storage mathematical model to solve the water storage prediction value. The true value of water storage for the corresponding area is retrieved from the original area-water storage lookup table. When the relative error between the predicted value and the true value does not exceed the preset value, it is determined that the optimal high-precision area-water storage mathematical model fits the accuracy requirements. The completed optimal high-precision area-water storage mathematical model is stored in a standardized manner, and the water body identification information, the function form of the optimal high-precision area-water storage mathematical model, the fitting parameter values and the coefficient of determination are recorded simultaneously. The domain parameter information of the optimal high-precision area-water storage mathematical model is also retained.
9. The method for model construction and water storage estimation based on inflection point identification according to claim 1, characterized in that, The process involves constructing corresponding area-water storage mathematical models for each water body within the target area to achieve automated batch modeling. Based on the latest water surface area obtained from remote sensing interpretation, the corresponding high-precision area-water storage mathematical models are used to calculate the water storage of each water body in batches, generating a water storage update report, including: Based on programming software, the operation of building an area-water storage mathematical model is performed on all water bodies to be tested in the target area by looping through them one by one, so as to achieve fully automated batch modeling. Once the area-water storage mathematical models for all water bodies within the target area are constructed, the resulting files will be automatically generated in batches, forming a complete area-water storage mathematical model library. The system employs a normalized water index to periodically obtain the latest water surface area data for each water body. It then calls the corresponding water body's area-water storage mathematical model to perform batch water storage estimation, synchronously updates the water storage data, and generates a water storage update report based on the updated water storage data.
10. The method for model construction and water storage estimation based on inflection point identification according to claim 9, characterized in that, The use of a normalized water index to periodically obtain the latest water surface area data for each water body includes: High-resolution optical remote sensing images of the target area are acquired periodically, and radiometric calibration and atmospheric correction are performed on the high-resolution optical remote sensing images in sequence. Based on high-resolution optical remote sensing images after radiometric calibration and atmospheric correction, the normalized water index is calculated using the spectral reflectance in the green and near-infrared bands. The threshold method is used to determine the threshold of the normalized water index, and the area where the threshold of the normalized water index is greater than or equal to the preset value is identified as the boundary of all water bodies in the study area. Collect vector layers of reservoirs and lakes from the latest annual land change survey data, match the identified water body boundaries with the land change survey vector boundaries in geographic information software, extract the lake and reservoir water surface boundaries, and extract the latest water surface area corresponding to each water body based on the lake and reservoir water surface boundaries.