Method for simulating underwater topography of natural lake based on lakeshore topographic features

Abstract

The invention discloses a method for simulating underwater topography of a natural lake based on lakeshore topographic features, and relates to the technical field of earth science. According to the method, the natural process of simulating lake surface descending is taken as a main line, and lake surrounding terrain information and morphological parameters extracted from DEM data are fully utilized, so that reliable estimation of the water depth and the water reserve of the natural lake under the condition of no actual measurement data is realized; the method has the following main advantages: 1) the method can automatically adjust the calculation formula of the underwater topography according to the lake morphological parameters, thereby realizing the autonomous adjustment of the water depth prediction of different types of lakes, and effectively improving the precision and reliability of the simulation result; and 2) in the underwater terrain simulation process, the remodeling effect of the lake sediment on the underwater terrain is considered, and a correction formula is introduced into the gradient calculation module, so that the gradient calculation result is closer to the real lakebed form, and the accuracy of the simulation process is further improved.

Description

Technical Field

[0001] The invention relates to the field of earth science technology, and in particular to a method for simulating underwater topography of a natural lake based on lakeshore topography characteristics. Background Art

[0002] At present, the most accurate method for measuring lake water reserves is to use sonar equipment installed on ships or lidar carried by drones for field measurements to obtain the bathymetric point data of the measured lake, and then generate a bathymetric map in triangulated network or raster image format through the bathymetric points to estimate the lake water reserves. However, these field bathymetric methods are costly and difficult to be widely used in large-scale lake water resource surveys. In response to the above difficulties, some scholars have proposed and developed a series of lake depth and water reserve estimation methods. According to their theoretical mechanisms, these methods are mainly divided into two categories: one is to integrate optical images and satellite altimetry data to construct the lake area-elevation relationship, and further estimate the lake water reserves on this basis; the other is to reconstruct the underwater terrain based on the geomorphic features around the lake to generate a bathymetric map, and derive the lake water reserves from the bathymetric map.

[0003] The first type of method combines optical imagery and satellite altimetry data. First, based on optical images from different periods, the changes in the lake surface area are extracted and an area sequence is generated. Secondly, the satellite altimetry data is processed and corrected, and the altimetry data passing through the lake is extracted to calculate the water level elevation of the lake. Subsequently, the temporally adjacent area sequence is matched with the extracted water level elevation to obtain the correspondence between area and elevation, and then the most suitable geometric volume calculation formula for a specific lake is used to approximate the lake water storage. This method can estimate the lake water storage more accurately, but the premise is that the lake has obvious dry seasons and flood seasons, such as seasonal lakes or intermittent lakes. For most permanent lakes, this method cannot obtain a complete area-elevation correspondence. Li et al. combined altimetry data from multiple satellites and Landsat satellite image data to derive the area-elevation relationship of the reservoir, applied the obtained relationship to the surface water body incidence data to obtain the water depth value of the dynamic area of ​​the reservoir, and finally used an extrapolation method to obtain a complete water depth data set. [1] Although Li was able to obtain a complete bathymetric map after introducing the extrapolation method, the accuracy of the water depth value in the extrapolated area depends largely on the extrapolation method, and the maximum water depth needs to be introduced as a constraint.

[0004] The second type of method is usually based on Digital Elevation Model (DEM) data to extract topographic and geomorphic features around the lake for inferring underwater topography and estimating water depth and water reserves. Since the water surface area in the DEM data cannot obtain complete underwater topographic information, it is generally set to a fixed value, that is, the elevation of the water surface area. Therefore, DEM data cannot directly reflect the underwater topography of the lake area. Some scholars have conducted statistical analysis on a large number of lakes and reservoirs with known water depth information through geostatistical models, and constructed a nonlinear model based on surrounding terrain information, aiming to reveal the relationship between lake area and average water depth or maximum water depth. Subsequently, water reserves are estimated based on average water depth and lake area. The geostatistical model developed by Messager et al. estimates the volume of lakes with a surface area of ​​more than 10 hectares worldwide based on the surrounding terrain. [2] . Due to the offsetting effect of underestimation and overestimation, the water storage estimation results are acceptable for large-scale studies. However, the uncertainty and potential bias at the specific lake or local scale cannot be ignored. In the study of water storage estimation of specific lakes, some scholars believe that the underwater topography of the lake can be regarded as an extension of the surrounding topography and can be inferred using extrapolation or interpolation methods. For example, Tseng et al. used DEM data to reconstruct the underwater topography of Lake Mead by linearly extrapolating the lakeshore elevation. [3] Liu et al. combined the skeleton method to divide multiple lake profiles along the long axis of the lake, and performed interpolation operations on the profiles with the bathymetric data as the constraint condition to obtain the lake profile elevation. They then used these elevation data to generate lake bathymetric maps and extracted the water storage of the lake from the bathymetric maps. [4] . These methods based on extrapolation or interpolation can restore underwater topographic changes and estimate water reserves to a certain extent, but they usually rely on prior measured data such as the maximum water depth of the lake to constrain the interpolation results, otherwise the water reserves of the lake may be seriously overestimated. In addition, due to the gradual accumulation of historical sediments in lakes, most modern lakes usually have relatively flat bottom features. However, when estimating water depth, the current DEM-based water reserve estimation method generally roughly truncates it with the measured water depth when the predicted water depth is greater than the measured water depth. The water depth map obtained by this truncation method still has a flat lake bottom area, but it has no direct connection with the lake's sedimentation process, and ignores the reshaping of the underwater topography by the lake sedimentation process.

[0005] Therefore, a new solution to the above problems needs to be proposed. Summary of the invention

[0006] The object of the present invention is to provide a method for simulating the underwater topography of a natural lake based on lakeshore topography characteristics, so as to solve the technical problems raised in the background technology.

[0007] To achieve the above object, the present invention provides the following technical solution: a method for simulating the underwater topography of a natural lake based on lakeshore topography characteristics, comprising at least the following steps: S1: Collect 30m resolution DEM data around the lake as input variables; S2: constructing an underwater terrain simulation algorithm based on the terrain continuity principle, wherein the underwater terrain simulation algorithm consists of a parameter initialization module and an underwater elevation calculation module; S3: The DEM data of underwater terrain elevation including the lake range is obtained through underwater terrain simulation algorithm. Combined with the given lake surface elevation or lake surface boundary, the water depth and water storage of the lake can be estimated.

[0008] Furthermore, the parameter initialization module in S2 is used to initialize parameters, and the parameters that need to be initialized include the water mask matrix and the lakeshore slope. The parameter initialization module needs to calculate intermediate parameters before initializing the parameters, and the intermediate parameters include at least the lake boundary, the buffer range and the slope map.

[0009] Furthermore, the calculation of lake boundary and water mask matrix comprises at least the following steps: Since the water pixel values ​​in the DEM data used are reassigned to the water surface elevation, the DEM values ​​of the water surface area need to be uniformly set to the water surface elevation, and the water-land binary classification is generated by screening the swb band values ​​in the DEM to distinguish between water bodies and land areas; in: 255: land, that is, in the result of land-water classification, a pixel value of 255 indicates that the location is a "land" area; 0: water body, that is, a pixel value of 0 indicates that the location is a "water body" area; The results of the water-land binary classification usually contain fragmented water patches, so a connected area detection algorithm is needed to retain the largest connected area as the lake body. After determining the main body of the lake, the lake boundary pixels are identified based on the edge detection function; The lake boundary pixels are assigned a value of 2, and the lake body pixels are assigned a value of 1 and stored in a new matrix with the same spatial coordinate system as the original data as the water mask matrix containing the lake boundary. The remaining pixels in the matrix are assigned a value of 0, that is, in the water mask matrix, the pixel value of 2 represents the lake boundary, the pixel value of 1 represents the water body, and the pixel value of 0 represents the land.

[0010] Furthermore, the calculation of the buffer range includes at least the following steps: The buffer range is determined by the comprehensive multi-level buffer, and the maximum buffer range of 21 pixels is first delineated; Secondly, multiple levels of buffers are delineated within the range at a distance of every three pixels, and the average elevation value within the range of each level of buffer is calculated; Therefore, each level of buffer is matched with an elevation value, and the turning point of the multi-level buffer elevation value is searched from the lake boundary to the outside, that is, the point where the elevation value changes from rising to falling. The area from the outside of the buffer to the lake shore boundary corresponding to this turning elevation value is the final buffer range; The elevation values ​​within this buffer range will be involved in subsequent calculations. If there is no turning point, 21 pixels will be used as the buffer range.

[0011] Furthermore, the calculation of the lakeshore slope involves three sets of data in the DEM: the elevation band, the lake boundary array, and the water mask matrix. In the water mask matrix, for each lake boundary pixel, that is, the pixel value is 2, the eight-neighborhood detection method is used to identify the land area in its domain, that is, the adjacent pixel value is 0; Extract the elevation value within the buffer range from the elevation band in the direction identified as land, and calculate the slope in that direction based on linear or polynomial fitting; After all calculations in the land direction are completed, the average of all fitted slope calculations is taken as the starting value for the underwater terrain slope calculation; Repeat the above steps for all lakeshore boundary pixels, map the lakeshore boundary pixel value domain to the slope value domain within the buffer zone, and reduce the impact of local extreme terrain on the calculation results.

[0012] Furthermore, the calculation of the slope map and the initial slope stretching at least comprises the following steps: Using the elevation band within the buffer area, the slope map of the buffer area is output based on the slope calculation method of ArcGIS. That is, the slope calculated using a 3×3 window can better represent the slope situation in the entire buffer area. The characteristics of slope map calculation are combined with those of lakeshore slope calculation, and the average value of lakeshore slope is selected as the criterion to decide whether to stretch the entire lakeshore slope. If the mean value of the lakeshore boundary pixels is greater than , indicating that the lakeshore area is relatively steep. Directly using the lakeshore slope calculation may lead to a serious overestimation of the water depth. At this time, the lakeshore slope will be stretched and transformed, and it will be reassigned according to the value range in the slope map.

[0013] Furthermore, the underwater elevation calculation module in S2 at least includes the following steps: The calculation formula for the pixel elevation in the underwater area is as follows:

[0014] In the formula, is the elevation value of the adjacent land pixel, is the distance between two pixels, The slope of the current calculated pixel; Because the The calculation of pixel elevation involves the elevation values ​​of neighboring pixels, that is, Therefore, the calculation order and calculation mode of underwater area pixels are particularly important; Based on the natural law that the lake surface elevation remains basically the same when the lake level drops, a cycle model is used to simulate and calculate the underwater terrain elevation; Slope calculation, in the formula, is a key input parameter. On the vertical section where the pixel is located, the underwater terrain is approximated as a quadratic function:

[0015] In the formula, Represents the elevation value, represents the independent variable related to the distance between the pixel and the lake shore; In order to simplify the calculation, it is necessary to assume a lowest point on the profile as the lake bottom point and establish a coordinate system on the profile. According to the real terrain change law, on the vertical profile of the lake, the lake bottom point is usually closer to the shore with steeper terrain. Therefore, the lake shore slopes on both sides of the profile are extracted, and the assumed lake bottom point in the slope calculation process is determined based on the slope values ​​on both sides. The horizontal distance between the lake bottom point and the left bank is calculated using the following formula (3):

[0016] In the formula, is the distance between the two banks of the lake on the current profile, and Represent the lakeshore slopes of the left and right banks respectively; After determining the horizontal position of the hypothetical lake bottom point, assuming that the slope value of the lake bottom point is infinitely close to 0, the slope of the target pixel is the derivative of the formula at this point, and the calculation formula is as follows:

[0017] In the formula, is the distance between the target pixel and the left bank; Since the sedimentation process in lakes will reshape the lake bottom topography and make it flat, it is necessary to correct the calculation results in the formula to make the simulation results closer to the actual underwater topography. The specific correction formula is as follows:

[0018] In the formula, is the corrected slope of the current calculated pixel, is the uncorrected current pixel slope, is the distance between the target pixel and the opposite bank in the current land boundary, It is a power index determined according to the characteristics of the lake, and its value range is as follows:

[0019] In the formula, The average value of the slope map within the buffer range in the parameter initialization module.

[0020] Furthermore, the specific process of calculating the underwater area pixel elevation includes at least the following steps: Perform multiple cycles, which imitate the process of lake level drop. Starting from the lake boundary, the lake level gradually drops. At this time, the underwater pixels will be gradually exposed along the shore to the lake bottom as the lake level drops. The areas with steep terrain have fewer underwater pixels exposed, while the areas with flat terrain have more underwater pixels exposed. According to the exposure order of the pixels, the elevation values ​​are calculated one by one. In the first cycle, starting from the current lake boundary, the elevations of all pixels adjacent to the boundary are calculated, and then the minimum elevation value is used as the stopping threshold of the first lake level drop process; It should be noted that once the elevation value of the underwater pixel is calculated, the pixel value at the corresponding position in the water mask matrix will change from 1 to 0, that is, it will change from water to land; For the current land boundary, continue to extract its inner boundary and calculate the elevation value of the inner boundary, select pixels that are not less than the stop threshold, and repeat the above operation until the elevation values ​​of all pixels in a certain round of calculation are not greater than the stop threshold, at which point the loop ends and enters the next loop; At the end of each cycle, the elevations of the lake boundary pixels tend to be consistent, and the simulation results are consistent with the lake level drop process. When there are no water pixels in the water mask matrix, that is, all values ​​are 0, the elevation calculation ends.

[0021] Furthermore, the specific estimation formula of the water reserve in S3 is as follows:

[0022] in is the volume of lake water, For the The volume of water pixels, For the The area of ​​a water pixel is numerically equal to the product of the length and width of the pixel. For the The simulated water depth of each water pixel.

[0023] Compared with the prior art, the present invention has the following beneficial effects: 1. The present invention takes "simulating the natural process of lake level drop" as the main line, makes full use of the terrain information and morphological parameters around the lake extracted from the DEM data, and realizes the reliable estimation of the water depth and water storage of natural lakes without measured data. The main advantages are as follows: 1) This method can automatically adjust the calculation formula of the underwater topography according to the lake morphological parameters, thereby realizing the autonomous adjustment of the water depth prediction of different types of lakes, and effectively improving the accuracy and reliability of the simulation results; 2) In the process of underwater topography simulation, the reshaping effect of lake sediments on the underwater topography is taken into account, and a correction formula is introduced in the slope calculation module to make the slope calculation result closer to the actual lake bottom morphology, further improving the accuracy of the simulation process.

[0024] 2. Compared with the existing lake water reserve prediction methods, the method of the present invention can actively adjust the calculation method through the extracted lake morphological parameters. The simulation results of this method are good for lakes of different areas and geometric shapes. In the absence of any measured prior data, the model only relies on the DEM data around the lake as input to achieve low-cost and high-efficiency water depth and water reserve estimation. At the same time, the simulation results are highly consistent with the field measurement results of radar depth sounders, etc., which can provide reliable data support for regional water resources management. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for describing the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0026] Figure 1 It is the technical roadmap of the present invention; Figure 2 is a composite schematic diagram of the present invention, wherein Figure 2 (a) is the determination diagram of the buffer in the parameter initialization module, where Figure 2 (b) and (c) are the calculation process of lakeshore slope; Figure 3 It is a schematic diagram of the simulated lake level drop process in the underwater elevation calculation link of the present invention; Figure 4 It is a schematic diagram of the parameters of the slope calculation process of the present invention and the difference of underwater terrain simulation before and after the slope correction; Figure 5 The slope calculation of the present invention and Schematic diagram of Figure 6 A three-dimensional water depth map of a lake simulated by the present invention; Figure 7A three-dimensional lake depth map generated based on measured data according to the present invention; Figure 8 It is a scatter plot of the measured water depth map and the simulated water depth map of the present invention. DETAILED DESCRIPTION

[0027] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0028] See also Figure 1 - Figure 8 , a method for simulating underwater topography of a natural lake based on lakeshore topography features, comprising at least the following steps: Step 1: Collect 30m resolution DEM data around the lake as input variables. See the table below for details:

[0029] Step 2: Construct an underwater terrain simulation algorithm based on the principle of terrain continuity. The algorithm consists of a parameter initialization module and an underwater elevation calculation module. The specific details are as follows: (1) Parameter initialization module The parameters that need to be initialized include the water mask matrix and the lakeshore slope. In initializing the above parameters, intermediate parameters such as the lake boundary, buffer range, and slope map must be calculated first.

[0030] Lake boundary and water mask matrix: Since the water pixel values ​​in the DEM data used this time are reassigned to the water surface elevation, the DEM values ​​of the water surface area need to be uniformly set to the water surface elevation, and the water-land binary classification is generated by screening the swb band values ​​in the DEM (255: land; 0: water body) to distinguish water and land areas. The results of the water-land binary classification usually contain fragmented water patches, and the connected area detection algorithm needs to be used to retain the largest connected area as the lake body. After determining the lake body, the lake boundary pixels are identified based on the edge detection function. The lake boundary pixels are assigned a value of 2, and the lake body pixels are assigned a value of 1 and stored in a new matrix with the same spatial coordinate system as the original data as the water mask matrix containing the lake boundary. The remaining pixels in the matrix are assigned a value of 0. That is, in the water mask matrix, the pixel value of 2 represents the lake boundary, the pixel value of 1 represents the water body, and the pixel value of 0 represents the land.

[0031] Buffer range: The buffer range is determined by multi-level buffers. First, the maximum buffer range of 21 pixels is delineated. Then, multi-level buffers are delineated within the range at a distance of every 3 pixels. The average elevation value within each level of buffer is calculated. Thus, each level of buffer matches an elevation value. The turning point of the multi-level buffer elevation value (i.e., the point where the elevation value changes from rising to falling) is searched from the lake boundary to the outside. The area from the outside of the buffer to the lakeshore boundary corresponding to this turning elevation value is the final buffer range. The elevation values ​​within this buffer range will participate in subsequent calculations. If there is no turning point, 21 pixels are used as the buffer range.

[0032] Lakeshore slope: This step involves three sets of data: the elevation band in the DEM, the lake boundary array, and the water mask matrix. Figure 2 As shown in (b), for each lake boundary pixel (pixel value is 2) in the water mask matrix, the eight-neighborhood detection method is used to identify the land area within its boundary (i.e., the adjacent pixel value is 0). Figure 2 As shown in (c), in the direction identified as land, the elevation value within the buffer range is extracted from the elevation band, and the slope in this direction is calculated based on linear or polynomial fitting. After all land directions are calculated, the average of all fitted slope calculations is taken as the starting value for underwater terrain slope calculation. Repeat the above steps for all lakeshore boundary pixels, and map their value range to the slope value range within the buffer range to reduce the impact of local extreme terrain on the local lakeshore on the calculation results.

[0033] Slope map and initial slope stretching: First, use the elevation band within the buffer range and output the slope map of the buffer range based on the slope calculation method of ArcGIS. The slope calculated by this method using a 3×3 window is more representative of the slope situation in the entire buffer area. The lakeshore slope calculated in c) is more representative of the slope from land to water. In order to comprehensively consider the characteristics of the two, we choose to use the average value of the lakeshore slope as the criterion to decide whether to stretch the entire lakeshore slope. If the average value of the lakeshore boundary pixel is greater than , indicating that the lakeshore area is relatively steep. Directly using the lakeshore slope calculation may lead to a serious overestimation of the water depth. At this time, the lakeshore slope will be stretched and transformed, and it will be reassigned according to the value range in the slope map.

[0034] (2) Underwater elevation calculation module Calculation of underwater pixel elevation: The calculation formula for the elevation of underwater area pixels is as follows:

[0035] In the formula, is the elevation value of the adjacent land pixel, is the distance between two pixels, is the slope of the current calculated pixel. The calculation of the pixel elevation involves the elevation values ​​of neighboring pixels (i.e. ), so the calculation order and calculation mode of underwater area pixels are particularly important. Based on the natural law that the lake surface elevation is basically the same when the lake level drops, a cycle mode is used to simulate the underwater terrain elevation. The specific process is as follows: Figure 3 As shown in the figure, this loop imitates the process of lake level drop. Starting from the lake boundary, the lake level gradually drops. At this time, the underwater pixels will be gradually exposed along the shore to the bottom of the lake as the lake level drops. There are fewer underwater pixels exposed in the steep area, and more underwater pixels exposed in the flat area. According to the exposure order of the pixels, the elevation values ​​are calculated one by one. In the first loop, starting from the current lake boundary, the elevations of all pixels adjacent to the boundary are calculated, and then the minimum elevation value is used as the stop threshold of the first lake level drop process. It should be noted that once the elevation value of the underwater pixel is calculated, the pixel value of the corresponding position in the water mask matrix will change from 1 to 0, that is, from water to land. For the current land boundary, continue to extract its inner boundary and calculate the elevation value of the inner boundary, select pixels that are not less than the stop threshold, and repeat the above operation until the elevation values ​​of all pixels in a certain round of calculation are not greater than the stop threshold. At this time, the loop ends and enters the next loop. At the end of each cycle, the elevation of the lake boundary pixels tends to be consistent, and the simulation results are consistent with the lake level drop process. When there is no water pixel in the water mask matrix (that is, the value is 1), the elevation calculation ends.

[0036] Slope calculation: In the formula, is the key input parameter. On the vertical section where the pixel is located, the underwater terrain is approximated as a quadratic function:

[0037] In the formula, Represents the elevation value, represents the independent variable related to the distance between the pixel and the lake shore. In order to simplify the calculation, we need to assume a lowest point on the profile as the lake bottom point and establish the following Figure 4 The coordinate system shown in Figure 1. According to the law of real terrain changes, in the vertical section of the lake, the lake bottom point is usually closer to the steeper shore. Therefore, we extract the lake shore slopes on both sides of the section and determine the assumed lake bottom point in the slope calculation process based on the slope values ​​on both sides. The horizontal distance between the lake bottom point and the left bank can be calculated using the following formula:

[0038] In the formula, is the distance between the two banks of the lake on the current profile, and Represent the lakeshore slopes of the left and right banks respectively. After determining the horizontal position of the hypothetical lake bottom point, we assume that the slope value of the lake bottom point is infinitely close to 0, then the slope of the target pixel is the derivative of the formula at this point, and the calculation formula is as follows:

[0039] In the formula, is the distance between the target pixel and the left bank. Indicates As in formula (2) As in formula (2) Enter, then the formula Derivation Then we get formula (4). Since the sedimentation process in lakes will reshape the lake bottom topography and make it flat, we need to correct the calculation results in the formula to make the simulation results closer to the real underwater topography. The specific correction formula is as follows:

[0040] In the formula, is the corrected slope of the current calculated pixel, is the uncorrected current pixel slope, is the distance between the target pixel and the opposite bank in the current land boundary, and The relationship is as Figure 5 As shown, It is a power index determined according to the characteristics of the lake, and its value range is as follows:

[0041] In the formula, The average value of the slope map within the buffer range in the parameter initialization module. Figure 4 As shown, the red curve represents the underwater terrain simulation result without slope correction, and the blue curve represents the underwater terrain simulation result after slope correction.

[0042] Step 3: Through the underwater terrain simulation algorithm, the DEM data of the underwater terrain elevation including the lake range is obtained. Combined with the given lake surface elevation or lake surface boundary, the water depth and water storage of the lake can be estimated. The specific estimation formula of water storage is as follows:

[0043] in is the volume of lake water, For the The volume of water pixels, For the The area of ​​a water pixel is numerically equal to the product of the length and width of the pixel. For the The simulated water depth of each water pixel.

[0044] To prove the advantages of this case over the prior art, the following comparison is made: The existing underwater terrain reconstruction methods based on DEM lack holistic considerations for the calculation of lakeshore slope or oversimplify the underwater elevation calculation process. [5] , obtain the original lake basin DEM, and obtain the lake basin DEM based on the hydrological analysis rules; then interpolate the underwater area in the lake basin DEM to be processed, and the interpolation method is simple linear interpolation, that is, assume that the slope of the underwater terrain is consistent with the slope on the shore, and the underwater terrain elevation decreases linearly at a fixed slope. Then merge the interpolation area with the surface area to obtain the lake basin DEM. Although the above method does not rely on measured water depth measurement data, linear interpolation obviously ignores the law of change of the actual underwater terrain, resulting in a large deviation between the interpolation result and the actual situation.

[0045] As the literature [6] . This method obtains the water boundary of the lake and reservoir from the lake and reservoir water and land binary matrix, and determines the slope of the water boundary pixel of the lake and reservoir by polynomial fitting according to the slope of the pixel within the preset range of the reverse slope of the water boundary pixel of the lake and reservoir, and obtains the slope matrix; according to the surrounding slope matrix and the current water surface elevation, by assuming that the water surface height is constantly decreasing, the simulated elevation of the pixels exposed to the water surface is iteratively calculated until all the water body pixels are calculated, and the simulated lake and reservoir underwater terrain matrix is ​​obtained; the simulated lake and reservoir underwater terrain matrix is ​​corrected by using the above-water part to obtain the final lake and reservoir underwater terrain matrix. In this method, the slope matrix is ​​calculated only by polynomial fitting, ignoring the impact of local terrain anomalies. Secondly, in the process of underwater terrain calculation, although the calculation is performed according to the pattern of water level drop, the calculation formula is single, and the reshaping of the terrain by lake sediments is not considered, and the simulation results are relatively fixed. In order to highlight the advantages of this method, this paper uses the literature [6] The method in the literature [7] The lake water reserves in the simulation are compared with the maximum water depth simulation results. The comparison results are shown in Table 2. The average percentage error of the water reserve estimation results of the method of the present application is about 8% lower than that of the comparison method, and the estimation accuracy of the maximum water depth is also higher. Because the water levels in different years are different, that is, the water reserves values ​​are also different in different years, the comparison results in Table 2 are the absolute percentage errors between the estimated water reserves and the measured water reserves.

[0046] Table 1 shows the water storage results derived from the simulated lake underwater topography

[0047] Table 2 shows the present invention and application literature [2] Comparison results of papers

[0048] In addition, this application shows scatter plots of simulated and measured bathymetric maps, e.g. Figure 8 As shown, the scattered points are randomly and evenly taken from the lake range, with a number of 2500, which can cover the entire lake. At the same time, the correlation coefficient is calculated , root mean square error , mean absolute error , the calculation formula is as follows:

[0049] In the formula, Indicates The true value of the data point, that is, the actual measured or observed value, Indicates The predicted value of a data point, that is, the prediction result generated by the model or algorithm based on the input data, Represents the mean of the true value, that is, the arithmetic mean of all true values, which is used to measure the central tendency of the data. Indicates the total number of data points, that is, the number of true values ​​and predicted values ​​is the same. In the above formula, It measures the model's ability to explain data variation, that is, the quality of the model's prediction. The closer it is to 1, the better the model's prediction effect. By taking the average and square root of the squared errors, the impact of larger errors is emphasized. It measures the average absolute difference between the predicted value and the true value, reflecting the average level of model error. Figure 8 According to the results in The average 9.26m, average It is 8.07m. Considering that the maximum water depth of the lake in this example ranges from 18 to 130m, it can be considered that the simulation results of this example perform well and can more realistically reflect the changing trend of underwater terrain.

[0050] References: LI, Y., GAO, H., JASINSKI, M. F., et al. Deriving high-resolutionreservoir bathymetry from ICESat-2 prototype photon-counting lidar andLandsat imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019,57(10): 7883-7893. DOI: 10.1109 / TGRS.2019.2917012. MESSAGER M L, LEHNER B, GRILL G, et al. Estimating the volume andage of water stored in global lakes using a geo-statistical approach[J].Nature Communications, 2016, 7(1): 13603. DOI:10.1038 / ncomms13603. TSENG, K. H., SHUM, C. K., KIM, J. W., et al. Integrating Landsatimageries and digital elevation models to infer water level change in HooverDam[J]. IEEE Journal of Selected Topics in Applied Earth Observations andRemote Sensing, 2016, 9(4): 1696-1709. DOI: 10.1109 / JSTARS.2015.2500599. LIU, K., SONG, C., ZHAN, P., et al. A low-cost approach for lakevolume estimation on the Tibetan Plateau: coupling the lake hypsometric curve and bottom elevation[J]. Frontiers in Earth Science, 2022, 10: 925944. DOI:10.3389 / feart.2022.925944. Ma Xiaoqi, Jing Haoxin, Zhu Bin. A method, equipment, device and medium for reconstructing underwater terrain of lakes and reservoirs based on DEM[P]. Shandong Province: CN202310777707.2, 2024-03-26. Lu Shanlong, Fang Chun, Li Mingyang, et al. Method and device for underwater three-dimensional terrain simulation of lakes and reservoirs based on digital elevation model[P]. Beijing: CN202110702775.3, 2021-10-19. FANG C, LU S, LI M, et al. Lake water storage estimation methodbased on similar characteristics of above-water and underwater topography[J]. Journal of Hydrology, 2023, 618: 129146. DOI:10.1016 / j.jhydrol.2023.129146. In summary: In view of the problems of water depth uncertainty and oversimplification of interpolation methods involved in the prediction of lake and reservoir water reserves, the present invention has designed for the first time an underwater terrain prediction idea of ​​"considering the influence of different lake sedimentation", making full use of the lake geometric feature information extracted from DEM, and realizing underwater terrain trend prediction and water reserve calculation. First, the initial parameters such as buffer range, slope, current lake boundary and water mask matrix are initialized based on the original DEM data; secondly, the slope calculation and correction formula are used to simulate the underwater terrain slope in combination with the reshaping effect of sediments on underwater terrain; then, based on the principle of elevation consistency when the lake level is lowered, the underwater terrain elevation change is cyclically calculated; finally, the lake water depth and water reserve are extracted from the simulated underwater terrain results. This method makes full use of DEM terrain information and takes into account the influence of sedimentation, thereby realizing fully automated estimation of lake water depth and water reserve, significantly improving the accuracy of lake and reservoir water reserve prediction in areas without measured data, and is suitable for hydrological and meteorological research, water resources management and other fields.

[0051] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other specific forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations falling within the meaning and scope of the equivalent elements of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

Claims

1. A method for simulating underwater topography of a natural lake based on lakeshore topography characteristics, characterized in that: At least the following steps are included: S1: Collect 30m resolution DEM data around the lake as input variables; S2: constructing an underwater terrain simulation algorithm based on the terrain continuity principle, wherein the underwater terrain simulation algorithm consists of a parameter initialization module and an underwater elevation calculation module; S3: Obtain DEM data of underwater terrain elevation including the lake range through underwater terrain simulation algorithm, and estimate the water depth and water storage corresponding to the lake in combination with the given lake surface elevation or lake surface boundary.

2. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 1, characterized in that: The parameter initialization module in S2 is used to initialize parameters. The parameters that need to be initialized include the water mask matrix and the lakeshore slope. The parameter initialization module needs to calculate intermediate parameters before initializing the parameters. The intermediate parameters include at least the lake boundary, the buffer range and the slope map.

3. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 2, characterized in that: The calculation of lake boundary and water mask matrix comprises at least the following steps: Since the water pixel values ​​in the DEM data used are reassigned to the water surface elevation, the DEM values ​​of the water surface area need to be uniformly set to the water surface elevation, and the water-land binary classification is generated by screening the swb band values ​​in the DEM to distinguish between water bodies and land areas; in: 255: land, that is, in the result of land-water classification, a pixel value of 255 indicates that the location is a "land" area; 0: water body, that is, a pixel value of 0 indicates that the location is a "water body" area; The results of the water-land binary classification usually contain fragmented water patches, so a connected area detection algorithm is needed to retain the largest connected area as the lake body. After determining the main body of the lake, the lake boundary pixels are identified based on the edge detection function; The lake boundary pixels are assigned a value of 2, and the lake body pixels are assigned a value of 1 and stored in a new matrix with the same spatial coordinate system as the original data as the water mask matrix containing the lake boundary. The remaining pixels in the matrix are assigned a value of 0, that is, in the water mask matrix, the pixel value of 2 represents the lake boundary, the pixel value of 1 represents the water body, and the pixel value of 0 represents the land.

4. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 3, characterized in that: The calculation of the buffer range comprises at least the following steps: The buffer range is determined by the comprehensive multi-level buffer, and the maximum buffer range of 21 pixels is first delineated; Secondly, multiple levels of buffers are delineated within the range at a distance of every three pixels, and the average elevation value within the range of each level of buffer is calculated; Therefore, each level of buffer is matched with an elevation value, and the turning point of the multi-level buffer elevation value is searched from the lake boundary to the outside, that is, the point where the elevation value changes from rising to falling. The area from the outside of the buffer to the lake shore boundary corresponding to this turning elevation value is the final buffer range; The elevation values ​​within this buffer range will be involved in subsequent calculations. If there is no turning point, 21 pixels will be used as the buffer range.

5. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 4, characterized in that: The calculation of lakeshore slope involves three sets of data: elevation band in DEM, lake boundary array and water mask matrix. In the water mask matrix, for each lake boundary pixel, that is, the pixel value is 2, the eight-neighborhood detection method is used to identify the land area in its domain, that is, the adjacent pixel value is 0; Extract the elevation value within the buffer range from the elevation band in the direction identified as land, and calculate the slope in that direction based on linear or polynomial fitting; After all calculations in the land direction are completed, the average of all fitted slope calculations is taken as the starting value for the underwater terrain slope calculation; Repeat the above steps for all lakeshore boundary pixels, map the lakeshore boundary pixel value domain to the slope value domain within the buffer zone, and reduce the impact of local extreme terrain on the calculation results.

6. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 5, characterized in that: The calculation of the slope map and the initial slope stretching at least include the following steps: Using the elevation band within the buffer area, the slope map of the buffer area is output based on the slope calculation method of ArcGIS. That is, the slope calculated using a 3×3 window can better represent the slope situation in the entire buffer area. The characteristics of slope map calculation are combined with those of lakeshore slope calculation, and the average value of lakeshore slope is selected as the criterion to decide whether to stretch the entire lakeshore slope. If the average value of the lakeshore boundary pixels is greater than tan(5°), it means that the lakeshore area is relatively steep. Directly using the lakeshore slope calculation will seriously overestimate the water depth. At this time, the lakeshore slope will be stretched and reassigned according to the value range in the slope map.

7. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 1, characterized in that: The underwater elevation calculation module in S2 at least comprises the following steps: The calculation formula for the pixel elevation in the underwater area is as follows: Z i =Z-d·k i (1) Where Z is the elevation value of the adjacent land pixel, d is the distance between two pixels, and k i The slope of the current calculated pixel; Since the calculation of the elevation of the i-th pixel involves the elevation values ​​of the neighboring pixels, that is, the value of Z, the calculation order and calculation mode of the pixels in the underwater area are particularly important; Based on the natural law that the lake surface elevation remains basically the same when the lake level drops, a cycle model is used to simulate and calculate the underwater terrain elevation; Slope calculation, in formula (1), k i is the key input parameter, which approximates the underwater terrain as a quadratic function on the vertical profile of the i-th pixel: Z=ax 2 +bx+c(2) In the formula, Z is the elevation value of the adjacent land pixel, and x represents the independent variable related to the distance between the pixel and the lake shore; In order to simplify the calculation, it is necessary to assume a lowest point on the profile as the lake bottom point and establish a coordinate system on the profile. According to the real terrain change law, on the vertical profile of the lake, the lake bottom point is usually closer to the steeper shore. Therefore, the lake shore slopes on both sides of the profile are extracted, and the assumed lake bottom point in the slope calculation process is determined based on the slope values ​​on both sides. The horizontal distance between the lake bottom point and the left bank is calculated using the following formula (3): Where D i is the distance between the two banks of the lake on the current profile, s l and r Represent the lakeshore slopes of the left and right banks respectively; After determining the horizontal position of the assumed lake bottom point, assuming that the slope value of the lake bottom point is infinitely close to 0, the slope of the target pixel is the derivative of formula (1) at this point, and the calculation formula is as follows: Where l is the distance between the target pixel and the left bank; Since the sedimentation process in lakes will reshape the lake bottom topography and make it flat, it is necessary to correct the calculation results in formula (4) to make the simulation results closer to the actual underwater topography. The specific correction formula is as follows: In the formula, k icorrection is the corrected slope of the current calculated pixel, k i is the uncorrected current pixel slope, d i is the distance between the target pixel and the opposite bank in the current land boundary, and n is the power exponent determined according to the characteristics of the lake, with the value range as follows: In the formula, S avg The average value of the slope map within the buffer range in the parameter initialization module.

8. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 7, characterized in that: The specific process of calculating the underwater area pixel elevation includes at least the following steps: Perform multiple cycles, which imitate the process of lake level drop. Starting from the lake boundary, the lake level gradually drops. At this time, the underwater pixels will be gradually exposed along the shore to the lake bottom as the lake level drops. The areas with steep terrain have fewer underwater pixels exposed, while the areas with flat terrain have more underwater pixels exposed. According to the exposure order of the pixels, the elevation values ​​are calculated one by one. In the first cycle, starting from the current lake boundary, the elevations of all pixels adjacent to the boundary are calculated, and then the minimum elevation value is used as the stopping threshold of the first lake level drop process; It should be noted that once the elevation value of the underwater pixel is calculated, the pixel value at the corresponding position in the water mask matrix will change from 1 to 0, that is, it will change from water to land; For the current land boundary, continue to extract its inner boundary and calculate the elevation value of the inner boundary, select pixels that are not less than the stop threshold, and repeat the above operation until the elevation values ​​of all pixels in a certain round of calculation are not greater than the stop threshold, at which point the loop ends and enters the next loop; At the end of each cycle, the elevations of the lake boundary pixels tend to be consistent, and the simulation results are consistent with the lake level drop process. When there are no water pixels in the water mask matrix, that is, all values ​​are 0, the elevation calculation ends.

9. The method for simulating underwater topography of a natural lake based on lakeshore topography characteristics according to claim 1, characterized in that: The specific estimation formula for the water reserves in S3 is as follows: V=∑v i =∑s i ×depth i (7) Where V is the volume of the lake, v i is the volume of the ith water pixel, s i is the area of ​​the ith water pixel, which is numerically equal to the product of the length and width of the pixel. i is the simulated water depth of the i-th water pixel.

Images