A water conservancy special terrain extraction method based on laser radar and satellite remote sensing

Abstract

The application discloses a special topography extraction method for water conservancy based on laser radar and satellite remote sensing, which comprises the following steps: collecting satellite radar data and satellite remote sensing data of a preset area; comparing and distinguishing the area according to the radar data and the remote sensing data to obtain a live area and a fuzzy area, extracting water conservancy boundaries and topography of the live area to obtain first data; cooperatively supplementing the fuzzy area to obtain water conservancy data, extracting water conservancy boundaries and topography of the water conservancy data by using closed contour lines to obtain second data; constructing a time-series water conservancy topology special topography extraction model according to the first data and the second data to extract water conservancy topography data DEM1 of the preset area; and constructing a two-dimensional water dynamic model based on the topography data and the water conservancy boundaries of the preset area to simulate flood inundation, extracting a water area according to an inundation range, and generating HDEM data of the preset area integrating water and land.

Description

Technical Field

[0001] This invention relates to the field of terrain extraction, and in particular to a method for terrain extraction in water conservancy based on lidar and satellite remote sensing. Background Technology

[0002] Topographic data is an important data support for water conservancy engineering design, watershed flood control simulation calculation, water resource management and allocation, soil erosion analysis for soil and water conservation, and is key to the construction of digital scenarios for watersheds, major rivers and major tributaries.

[0003] The precise design of water conservancy projects and the dynamic monitoring of watersheds heavily rely on the spatiotemporal accuracy of topographic data, especially the extraction of key water conservancy elements such as river boundaries, inundation simulations, and riverbed evolution. Traditional topographic mapping methods suffer from low efficiency, limited coverage, and an inability to capture dynamic changes, making it difficult to meet the demands of modern water conservancy informatization for real-time, comprehensive, and high-precision data. The development of lidar and satellite remote sensing technologies has provided new avenues for topographic data acquisition: lidar can accurately invert the three-dimensional morphology of the Earth's surface through point cloud data, while satellite remote sensing can achieve large-scale synchronous observation through radar, multispectral imagery, and other methods. The combination of these two technologies makes rapid and refined extraction of water conservancy topography possible.

[0004] However, existing technologies face the following challenges when fusing multi-source data: the imaging mechanisms of lidar data and satellite remote sensing data differ significantly; the former is sensitive to surface geometric features, while the latter is sensitive to spectral features, and direct fusion can easily lead to feature misalignment; in water conservancy scenarios, topographic data has strong time sensitivity, and existing single-source data model cannot take into account both spatial resolution and temporal continuity; water conservancy topology models need to reflect the evolution of topography over time, but traditional models are mostly built based on static data and lack the ability to dynamically fuse the "radar-remote sensing" time-series data chain, making it difficult to simulate dynamic processes such as water flow direction and flooding range expansion.

[0005] Furthermore, high-precision DEMs of buildings within water-affected areas are produced using technologies such as UAV aerial photography and LiDAR. Current methods for producing high-precision DEMs of water-affected areas suffer from the following problems: conventional DEMs do not consider the water-passing and water-blocking functions of water-affected structures, leading to distortions in topographic and water flow connectivity; additionally, the point cloud classification of structures is generally done manually, which is inefficient; conventional DEMs typically do not include underwater topographic information, resulting in unreasonable representations of water volume and flow in hydrological simulations, and failing to support accurate calculations of runoff generation and confluence, and precise simulations and extrapolations of watershed floods.

[0006] Therefore, there is an urgent need for a water conservancy topography extraction method based on multi-source data collaboration to fill the gap in watershed underlying surface topography production methods that take into account water flow connectivity, solve the problems of existing conventional topography that cannot truly express water flow connectivity and lack of underwater topography, improve the topography extraction accuracy in complex topography and dynamic scenarios, and thus improve the accuracy of runoff generation and confluence, flood simulation and other data on the watershed surface. Summary of the Invention

[0007] The purpose of this invention is to provide a method for extracting terrain specific to water conservancy based on lidar and satellite remote sensing.

[0008] To achieve the above objectives, the present invention is implemented according to the following technical solution:

[0009] This invention includes the following steps:

[0010] Collect satellite radar data and satellite remote sensing data of a preset area, and preprocess the radar data and the remote sensing data;

[0011] The radar data and the remote sensing data are compared and the regions are divided. The consistent regions are taken as the real regions and the inconsistent regions are taken as the fuzzy regions. The water conservancy boundaries and terrain are extracted from the real regions to obtain the first data.

[0012] Collaborative supplementary surveys are performed on the ambiguous area to obtain covered water conservancy data. Closed contour lines are used to extract water conservancy boundaries and terrain from the covered water conservancy data to obtain second data. The second data includes the terrain of the second region and the second water conservancy boundary.

[0013] Based on the first data and the second data, a time-series hydraulic topology-specific terrain extraction model is constructed to extract hydraulic topography data DEM1 for a preset area.

[0014] A two-dimensional hydrodynamic model is constructed based on the pre-defined terrain data and hydraulic boundaries to simulate flood inundation, and the water-affected areas are extracted according to the inundation range.

[0015] For details of water-related structures, LiDAR is used to acquire point cloud data of the water-related structure area. According to the functional characteristics, it is divided into water-blocking structures and water-passing structures. Graphical processing and a combination of addition and removal operations are performed, and it is fused with the underwater DEM. Weighted correction and Gaussian smoothing algorithms are used to embed it into the preset area water and land topographic data DEM1 to generate the preset area integrated water and land HDEM data.

[0016] Furthermore, the method for comparing and determining the region based on the radar data and the remote sensing data includes:

[0017] Edge features are extracted from optical images obtained from radar and satellite remote sensing data. Sensor registration is performed by fusing these edge features. Temporal weights are introduced to perform multi-temporal matching on the spatially registered radar and satellite remote sensing data. The expression is as follows:

[0018]

[0019] The time series weight at time point t is: Adjusted for seasonal variations, the overall similarity score is: The gradient component of the template image in the x-direction at point v is: The gradient component of the search image in the x-direction at point v is: The gradient component of the template image in the y-direction at point v is: The gradient component of the image in the y-direction at point v is: The total number of features is The number of timing templates is ;

[0020] Radar and satellite remote sensing data are grouped, with those having a matching degree greater than 0.673 grouped together. Based on the registered multi-source data, the spatiotemporal consistency of different groups is calculated:

[0021]

[0022] Where spatiotemporal consistency is The satellite remote sensing derived point cloud at time t is The radar point cloud at time t is The upper limit of observation time is The point cloud structure similarity is The point cloud is projected onto a 2.5D grid, the mean elevation and variance of intensity within the grid are calculated, and the SSIM comparison term is obtained based on the mean elevation and variance of intensity.

[0023] The temporal variance is obtained by summing the temporal standard deviations of spatiotemporal consistency. When the pixel-level consistency is greater than or equal to 0.786 and the temporal variance is less than 0.352, the elevation is stable, the intensity-vegetation index is strongly correlated, and the location of the group is the actual area. When the pixel-level consistency is less than 0.786 and the temporal variance is greater than 0.352, the elevation changes abruptly, the intensity-vegetation index is weakly correlated, and the location of the group is the ambiguous area.

[0024] Furthermore, the method for extracting the hydraulic boundaries and topography of the actual area includes:

[0025] The terrain of the first region is obtained by extracting terrain features from the actual area using SAR technology. Potential water bodies in the actual area are located by locating the slope threshold. The NDWI index is binarized in the potential water body area using the maximum inter-class variance method to extract the initial water body boundary. For the cloud-covered area of ​​the optical image, the water backscattering coefficient of the SAR image is used to fill the boundary missing. The two are then merged to generate the complete water body outline of the actual area, which is used as the first water boundary. The first water boundary and the terrain of the first region are output as the first data.

[0026] Furthermore, the method for obtaining topographic and overlay water conservancy data of the second region by collaborative supplementary surveying of the fuzzy region includes:

[0027] Construct an integrated "air-space-ground" observation network to achieve spatiotemporal data complementarity, acquire medium-density lidar point cloud data, radar back reflection coefficient and optical vegetation index, identify fuzzy region types through random forest model, and obtain high-priority regions, medium-priority regions and low-priority regions.

[0028] For high-priority areas, UAV data is introduced. The UAV altitude is equal to the maximum height of vegetation plus a safety margin. The overlap rate is proportional to the vegetation complexity. A mesh flight path is generated based on the LiDAR scanning angle.

[0029] The drone scans the same area from different directions, and the data is fused using a point cloud registration algorithm. The expression is:

[0030]

[0031] Where the rotation matrix is ​​A, the translation vector is w, and the u-th point in the source point cloud is... The source point cloud consists of point cloud data acquired from different scanning angles, which are to be registered to the target point cloud. The u-th point in the target point cloud is... The directional weighting coefficient is The normal vector is The merged data is G;

[0032] Spatiotemporal registration is performed on multi-temporal radar images and laser data of the same area. The real terrain is estimated using Kalman filtering and used as the second-region terrain data for the blurred area. For multi-echo lidar point cloud areas, the first echo is separated, the echo intensity difference is calculated, and potential river centerlines are identified by searching for local minimum points in areas with high echo intensity differences. For multi-source feature fusion areas, linear shadow features are extracted from optical images, and anomalous stripes of backscattering coefficients are analyzed from radar images. Combined with water conservancy planning data, spatial overlay analysis is used to locate the location information of water-related structures. The combined second-region terrain data, potential river centerlines, and location information of water-related structures are output as comprehensive water conservancy data.

[0033] Furthermore, the method for extracting the second data includes:

[0034] Based on the covered water conservancy data, the elevation gradient and flow direction of the river centerline neighborhood are obtained, and the robust fuzzy local information C-means clustering algorithm is used to initially separate water bodies and non-water bodies.

[0035] Given the objective function of fuzzy clustering, the expression is:

[0036]

[0037] Where the objective function is The prior probability that pixel a belongs to cluster c is The terrain dissimilarity measure between pixel a and cluster c is: The fuzzy local variables of pixel a and cluster c are: The ridge regression coefficient is The total number of pixels is The membership degree of pixel a to cluster c is The spectral dissimilarity measure between pixel a and cluster c is The balance parameter is Balancing spectral and topographic features, the elevation variation coefficient of pixel a is The elevation variation coefficient of cluster c is ;

[0038] Calculate the mean, covariance, and dissimilarity measures:

[0039]

[0040]

[0041]

[0042] The mean of cluster c is covariance is The transpose of the matrix is The spectral measure of pixel a is Data dimensions are The mean of pixel a in cluster c is ;

[0043] Calculate the local variation density and fuzzy local variables:

[0044]

[0045]

[0046] The local variation density of pixel a is The neighboring pixels of pixel a are The set of neighboring pixels is The mean of the spectral measurements of the neighboring pixels centered at pixel a is The variance of the spectral test of the neighborhood pixels centered at pixel a is Neighboring pixels The local variation density is ;

[0047] Calculate prior probability and membership degree:

[0048]

[0049]

[0050] The index clustering is as follows: The domain effect strength is The pixel a attribute is The number of clusters is The prior probability that pixel a belongs to cluster i is The spectral dissimilarity measure between pixel a and cluster i is The fuzzy local variables of pixel a and cluster i are ;

[0051] If the maximum error of the objective function from time point t to time point t+1 is less than the error threshold, then the membership matrix is ​​defuzzified to obtain the binary segmentation result. Morphological closing operation is used to fill the small holes in the binary segmentation result and extract the contour of the largest connected component.

[0052] Apply hydraulic set constraints to the contour of the largest connected region, forcing the hydraulic cluster centers to satisfy... Isolated clusters with areas smaller than a threshold are removed to obtain segmented data; the neighborhood elevation gradient is... The direction of water flow is ;

[0053] The gradient and gradient magnitude of the fuzzy water area are obtained from the segmented data. Anisotropic diffusion filtering is applied to the segmented data to preserve the sharpness of the channel boundaries. The expression is as follows:

[0054]

[0055]

[0056] The edge stopping function is: To control the diffusion intensity, the edge sensitivity coefficient in anisotropic diffusion filtering is... It adaptively adjusts based on the boundary gradient, with a gradient magnitude of s and an image intensity gradient of... ;

[0057] The segmented data boundary is used as the second water conservancy boundary data, and combined with the topographic data of the second region, the output is the second data.

[0058] Furthermore, the method for constructing a time-series hydraulic topology-specific terrain extraction model based on the first data and the second data includes:

[0059] A basic terrain grid is constructed based on the terrain data of the first and second datasets. Static parameters are extracted, and temporal difference analysis is performed on the hydraulic boundaries of the first and second datasets to identify boundary change areas. Vector-raster conversion is performed on the dynamic boundaries to generate a temporal raster layer and extract the dynamic topology node set. Static parameters include terrain elevation, slope, and aspect.

[0060] Dynamic nodes are matched with static nodes, and a spatiotemporal joint node set is generated through spatial overlay analysis. A node attribute table is established to record the spatial coordinates, temporal attributes, and hydraulic attributes of the nodes.

[0061] Based on the spatiotemporal joint node set, a dynamic topological edge set is constructed using the triangulation algorithm of the point set. Here, the adjacency relationship represents the direct connectivity between nodes, the temporal association represents the state change of the same node at different time steps, and the edge attributes include distance, flow direction, and water passage capacity.

[0062] Distance, water level difference, and time step are used as time-series weights for dynamic topological edges to construct a time-series hydraulic topology-specific terrain extraction model, and output the hydraulic topographic data DEM1 of the preset area.

[0063] Furthermore, the method for obtaining the water-affected area includes:

[0064] The first and second water conservancy boundaries are merged into a complete water conservancy boundary to form river network data;

[0065] The pre-defined regional topographic data DEM1 is divided into grids. The area adjacent to the water conservancy boundary and the area covering the water conservancy data are divided into grids with a grid size of 5m×5m, while other areas far from the river channel are divided into grids with a grid size of 20m×20m. The roughness, area coefficient, time step and other hydrodynamic parameters are assigned to the grid after division, and the grid model results are output.

[0066] Based on historical flow data of the preset area, a flood frequency curve is fitted to determine the design peak flow corresponding to the target frequency, which serves as the inflow boundary condition.

[0067] Based on the control volume method, the two-dimensional Saint-Venant shallow water equation is used to simulate flood inundation with the maximum design peak flow of the river in the preset area as the boundary condition, and the inundation range is output as the water-affected area.

[0068] The beneficial effects of this invention are:

[0069] This invention is a method for extracting terrain features specifically for water conservancy based on lidar and satellite remote sensing. Compared with existing technologies, this invention has the following technical advantages:

[0070] The proposed method for constructing a water conservancy-specific HDEM (High-Density Earth Model) considers the water flow and water retention characteristics of water-related structures such as dikes, reservoirs, dams, sluices, ponds, rubber dams, water-retaining roads, bridges, culverts, aqueducts, and inverted siphons on the basis of conventional DEMs. It refines and removes topographic features by adding or removing features, while also integrating underwater topographic data. This method can express the true connectivity of water flow and improve the accuracy of hydrological simulation analysis. It meets the needs of digital twin watershed applications such as the construction of three-dimensional scenes of the underlying surface of water conservancy basins, hydrological modeling and calculation, accurate calculation of runoff generation and runoff, accurate simulation and extrapolation of basin floods, assessment of flood inundation range and loss, and analysis of water balance relationships within the basin. It lays a solid foundation for the realization of precise scheduling of "one flow rate, one cubic meter of reservoir capacity, and one centimeter of water level" in watershed flood control and early warning, and greatly improves the level of digitalization, precision, and practicality of water conservancy.

[0071] This invention improves the accuracy of water conservancy-specific topographic HDEMs by employing preprocessing, comparative regional division, extraction of first-region topographic data and first-region water conservancy boundary, multi-angle joint inversion, collaborative supplementary surveying, extraction of second-region topographic data and second-region water conservancy boundary, model construction, acquisition of water-related areas, and acquisition of HDEM results. Optimizing the extraction of water conservancy-specific topographic HDEMs significantly saves resources and improves work efficiency. It enables intelligent extraction of water conservancy-specific topographic HDEMs, real-time multi-angle joint inversion of the extracted data, and extraction of second-region topographic and water conservancy boundaries. This is of great significance for the extraction of water conservancy-specific topographic HDEMs, meeting the extraction requirements of different standards and possessing a certain degree of universality. It is also of great importance for the promotion and application of water conservancy-specific topographic HDEMs. Attached Figure Description

[0072] Figure 1 This is a flowchart illustrating the steps of a water conservancy-specific terrain extraction method based on lidar and satellite remote sensing according to the present invention.

[0073] Figure 2 This is a schematic diagram of the HDEM results for the bridge area in an embodiment of the present invention;

[0074] Figure 3 This is a schematic diagram of HDEM results for the reservoir dam area in an embodiment of the present invention;

[0075] Figure 4 This is a schematic diagram of the HDEM results of the sluice gate area in an embodiment of the present invention. Detailed Implementation

[0076] The present invention will be further described below through specific embodiments. The illustrative embodiments and descriptions herein are used to explain the present invention, but are not intended to limit the present invention.

[0077] The present invention discloses a method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing, comprising the following steps:

[0078] like Figure 1-2 As shown, this embodiment includes the following steps:

[0079] Collect satellite radar data and satellite remote sensing data of a preset area, and preprocess the radar data and the remote sensing data;

[0080] In the actual assessment, the XX water conservancy tributary was taken as the research object. Satellite radar data: Sentinel-1 spaceborne SAR, L-band, incident angle 35°, resolution 10m was used to acquire images of the flood season and the dry season. Airborne lidar points, cloud density 12 points / m², generated 0.5m resolution, including the three-dimensional structure of the dike and river channel; FMCW ground-based radar, resolution 0.1m, was used to perform local scanning of typical areas.

[0081] Satellite remote sensing data was used to extract the NDWI index, with a threshold of 0.4, to identify water bodies.

[0082] The radar data and the remote sensing data are compared and the regions are divided. The consistent regions are taken as the real regions and the inconsistent regions are taken as the fuzzy regions. The water conservancy boundaries and terrain are extracted from the real regions to obtain the first data.

[0083] In the actual assessment, 32% of the area was a vague area, while 68% of the area was a real area;

[0084] Collaborative supplementary surveys are performed on the ambiguous area to obtain covered water conservancy data. Closed contour lines are used to extract water conservancy boundaries and terrain from the covered water conservancy data to obtain second data. The second data includes the terrain of the second region and the second water conservancy boundary.

[0085] Based on the first data and the second data, a time-series hydraulic topology-specific terrain extraction model is constructed to extract hydraulic topography data DEM1 for a preset area.

[0086] A two-dimensional hydrodynamic model is constructed based on the pre-defined terrain data and hydraulic boundaries to simulate flood inundation, and the water-affected areas are extracted according to the inundation range.

[0087] For details of water-related structures, LiDAR is used to acquire point cloud data of the water-related structure area. According to the functional characteristics, it is divided into water-blocking structures and water-passing structures. Graphical processing and a combination of addition and removal operations are performed, and it is fused with the underwater DEM. Weighted correction and Gaussian smoothing algorithms are used to embed it into the preset area water and land topographic data DEM1 to generate the preset area integrated water and land HDEM data.

[0088] In this embodiment, the method for comparing and determining region division based on the radar data and the remote sensing data includes:

[0089] Edge features are extracted from optical images obtained from radar and satellite remote sensing data. Sensor registration is performed by fusing these edge features. Temporal weights are introduced to perform multi-temporal matching on the spatially registered radar and satellite remote sensing data. The expression is as follows:

[0090]

[0091] The time series weight at time point t is: Adjusted for seasonal variations, the overall similarity score is: The gradient component of the template image in the x-direction at point v is: The gradient component of the search image in the x-direction at point v is: The gradient component of the template image in the y-direction at point v is: The gradient component of the image in the y-direction at point v is: The total number of features is The number of timing templates is ;

[0092] Radar and satellite remote sensing data are grouped, with those having a matching degree greater than 0.673 grouped together. Based on the registered multi-source data, the spatiotemporal consistency of different groups is calculated:

[0093]

[0094] Where spatiotemporal consistency is The satellite remote sensing derived point cloud at time t is The radar point cloud at time t is The upper limit of observation time is The point cloud structure similarity is The point cloud is projected onto a 2.5D grid, the mean elevation and variance of intensity within the grid are calculated, and the SSIM comparison term is obtained based on the mean elevation and variance of intensity.

[0095] The temporal variance is obtained by summing the temporal standard deviations of spatiotemporal consistency. When the pixel-level consistency is greater than or equal to 0.786 and the temporal variance is less than 0.352, the elevation is stable, the intensity-vegetation index is strongly correlated, and the location of the group is the actual area. When the pixel-level consistency is less than 0.786 and the temporal variance is greater than 0.352, the elevation changes abruptly, the intensity-vegetation index is weakly correlated, and the location of the group is the ambiguous area.

[0096] In this embodiment, the method for extracting the hydraulic boundaries and topography of the actual area includes:

[0097] The terrain of the first region is obtained by extracting terrain features from the actual area using SAR technology. Potential water bodies in the actual area are located by locating the slope threshold. The NDWI index is binarized in the potential water body area using the maximum inter-class variance method to extract the initial water body boundary. For the cloud-covered area of ​​the optical image, the water backscattering coefficient of the SAR image is used to fill the boundary missing. The two are then merged to generate the complete water body outline of the actual area, which is used as the first water boundary. The first water boundary and the terrain of the first region are output as the first data.

[0098] In this embodiment, the method for obtaining topographic and overlay water conservancy data of the second region by collaborative supplementary surveying of the fuzzy region includes:

[0099] Construct an integrated "air-space-ground" observation network to achieve spatiotemporal data complementarity, acquire medium-density lidar point cloud data, radar back reflection coefficient and optical vegetation index, identify fuzzy region types through random forest model, and obtain high-priority regions, medium-priority regions and low-priority regions.

[0100] For high-priority areas, UAV data is introduced. The UAV altitude is equal to the maximum height of vegetation plus a safety margin. The overlap rate is proportional to the vegetation complexity. A mesh flight path is generated based on the LiDAR scanning angle.

[0101] The drone scans the same area from different directions, and the data is fused using a point cloud registration algorithm. The expression is:

[0102]

[0103] Where the rotation matrix is ​​A, the translation vector is w, and the u-th point in the source point cloud is... The source point cloud consists of point cloud data acquired from different scanning angles, which are to be registered to the target point cloud. The u-th point in the target point cloud is... The directional weighting coefficient is The normal vector is The merged data is G;

[0104] Spatiotemporal registration is performed on multi-temporal radar images and laser data of the same area. The real terrain is estimated using Kalman filtering and used as the second-region terrain data for the blurred area. For multi-echo lidar point cloud areas, the first echo is separated, the echo intensity difference is calculated, and potential river centerlines are identified by searching for local minimum points in areas with high echo intensity differences. For multi-source feature fusion areas, linear shadow features are extracted from optical images, and anomalous stripes of backscattering coefficients are analyzed from radar images. Combined with water conservancy planning data, spatial overlay analysis is used to locate the location information of water-related structures. The combined second-region terrain data, potential river centerlines, and location information of water-related structures are output as comprehensive water conservancy data.

[0105] In this embodiment, the method for extracting the second data includes:

[0106] Based on the covered water conservancy data, the elevation gradient and flow direction of the river centerline neighborhood are obtained, and the robust fuzzy local information C-means clustering algorithm is used to initially separate water bodies and non-water bodies.

[0107] Given the objective function of fuzzy clustering, the expression is:

[0108]

[0109] Where the objective function is The prior probability that pixel a belongs to cluster c is The terrain dissimilarity measure between pixel a and cluster c is: The fuzzy local variables of pixel a and cluster c are: The ridge regression coefficient is The total number of pixels is The membership degree of pixel a to cluster c is The spectral dissimilarity measure between pixel a and cluster c is The balance parameter is Balancing spectral and topographic features, the elevation variation coefficient of pixel a is The elevation variation coefficient of cluster c is ;

[0110] Calculate the mean, covariance, and dissimilarity measures:

[0111]

[0112]

[0113]

[0114] The mean of cluster c is covariance is The transpose of the matrix is The spectral measure of pixel a is Data dimensions are The mean of pixel a in cluster c is ;

[0115] Calculate the local variation density and fuzzy local variables:

[0116]

[0117]

[0118] The local variation density of pixel a is The neighboring pixels of pixel a are The set of neighboring pixels is The mean of the spectral measurements of the neighboring pixels centered at pixel a is The variance of the spectral test of the neighborhood pixels centered at pixel a is Neighboring pixels The local variation density is ;

[0119] Calculate prior probability and membership degree:

[0120]

[0121]

[0122] The index clustering is as follows: The domain effect strength is The pixel a attribute is The number of clusters is The prior probability that pixel a belongs to cluster i is The spectral dissimilarity measure between pixel a and cluster i is The fuzzy local variables of pixel a and cluster i are ;

[0123] If the maximum error of the objective function from time point t to time point t+1 is less than the error threshold, then the membership matrix is ​​defuzzified to obtain the binary segmentation result. Morphological closing operation is used to fill the small holes in the binary segmentation result and extract the contour of the largest connected component.

[0124] Apply hydraulic set constraints to the contour of the largest connected region, forcing the hydraulic cluster centers to satisfy... Isolated clusters with areas smaller than a threshold are removed to obtain segmented data; the neighborhood elevation gradient is... The direction of water flow is ;

[0125] The gradient and gradient magnitude of the fuzzy water area are obtained from the segmented data. Anisotropic diffusion filtering is applied to the segmented data to preserve the sharpness of the channel boundaries. The expression is as follows:

[0126]

[0127]

[0128] The edge stopping function is: To control the diffusion intensity, the edge sensitivity coefficient in anisotropic diffusion filtering is... It adaptively adjusts based on the boundary gradient, with a gradient magnitude of s and an image intensity gradient of... ;

[0129] The segmented data boundary is used as the second water conservancy boundary data, and combined with the topographic data of the second region, the output is the second data.

[0130] In this embodiment, the method for constructing a time-series hydraulic topology-specific terrain extraction model based on the first data and the second data includes:

[0131] A basic terrain grid is constructed based on the terrain data of the first and second datasets. Static parameters are extracted, and temporal difference analysis is performed on the hydraulic boundaries of the first and second datasets to identify boundary change areas. Vector-raster conversion is performed on the dynamic boundaries to generate a temporal raster layer and extract the dynamic topology node set. Static parameters include terrain elevation, slope, and aspect.

[0132] Dynamic nodes are matched with static nodes, and a spatiotemporal joint node set is generated through spatial overlay analysis. A node attribute table is established to record the spatial coordinates, temporal attributes, and hydraulic attributes of the nodes.

[0133] Based on the spatiotemporal joint node set, a dynamic topological edge set is constructed using the triangulation algorithm of the point set. Here, the adjacency relationship represents the direct connectivity between nodes, the temporal association represents the state change of the same node at different time steps, and the edge attributes include distance, flow direction, and water passage capacity.

[0134] Distance, water level difference, and time step are used as time-series weights for dynamic topological edges to construct a time-series hydraulic topology-specific terrain extraction model, and output the hydraulic topographic data DEM1 of the preset area.

[0135] In this embodiment, the method for obtaining the water-affected area includes:

[0136] The first and second water conservancy boundaries are merged into a complete water conservancy boundary to form river network data;

[0137] The pre-defined regional topographic data DEM1 is divided into grids. The area adjacent to the water conservancy boundary and the area covering the water conservancy data are divided into grids with a grid size of 5m×5m, while other areas far from the river channel are divided into grids with a grid size of 20m×20m. The roughness, area coefficient, time step and other hydrodynamic parameters are assigned to the grid after division, and the grid model results are output.

[0138] Based on historical flow data of the preset area, a flood frequency curve is fitted to determine the design peak flow corresponding to the target frequency, which serves as the inflow boundary condition.

[0139] Based on the control volume method, the two-dimensional Saint-Venant shallow water equation is used to simulate flood inundation with the maximum design peak flow of the river in the preset area as the boundary condition, and the inundation range is output as the water-affected area.

[0140] In this embodiment, the method for obtaining the HDEM results includes:

[0141] Within the water-involved area, high-density lidar point cloud data of water-related structures is obtained by combining water conservancy planning data and coverage water conservancy data. A region-growing building classification and filtering algorithm based on planar fitting technology is used to extract water-related structures such as dikes, dams, culverts, bridges, sluices, flood walls, flood shelters, aqueducts, and inverted siphons. Within the water-involved area, a region-growing algorithm is used to segment the point cloud data into building footprints. The region-growing algorithm, based on planar fitting technology, is used to segment the ground feature footprints and extract building footprint blocks: ① An area set is created for adjacent building footprints. To begin creating a region, arbitrarily select a non-ground node as the initial seed point. Regions are formed by recursively connecting the eight neighboring non-ground nodes of the seed point. This process is repeated until every non-ground node is assigned to a region. ② Each region Divided into and i consists of two parts, among which Indicates the footer within the block. A series of edge points represent a block. If at least one laser point in a point's eight-neighborhood is a ground point, then that point is defined as an edge point; otherwise, it is defined as an interior point. ③ If a small area... If the area is less than the preset minimum area threshold or the number of laser footprints inside it is zero, then it is removed from the region set. Remove from ④ For the region Each small area Perform the following loop: Given an interior point , get and The best-fit plane of eight neighboring points, plane equation ,parameter By solving the deviation Minimum of the sum of squares Obtain; ⑤ Define As a small block obtained through the above process Let be a set. , ,in It's a small piece The first in If the represented small block is smaller than the preset minimum area threshold, then... If it is not completely contained by other blocks, then it will from Remove from middle; ⑥ merge Connect the blocks in the middle and get a new set. ⑦ Obtain the point cloud block of the roof surface of the water-related building, if If the building area is less than the preset minimum building area threshold, from Remove from ;final The remaining blocks are all used as point cloud blocks on the top surface of the water-related buildings; ⑧ Combined with manual correction, ensure that the topographic morphology of the water-related building area is completely and accurately expressed;

[0142] A TIN was constructed using point clouds on the top surfaces of water-related structures and supplementary topographic feature lines. A regular grid was interpolated to create a DEM2 file for the water-related structures. Based on the functional characteristics of the water-related structures, water-retaining structures such as dikes, dams, culverts, bridges, sluices, flood walls, water-repelling platforms, aqueducts, and inverted siphons were stored as DEM3 files as topographic features to be added to DEM1. Water-passing structures such as culverts and bridges were stored as DEM4 files as topographic features to be removed from DEM1.

[0143] The DEM1 and the DEM3 and DEM4 of the water-related structures are processed through image processing. The water-retaining structure DEM3 is added to the DEM1, and the water-passing structure DEM4 is removed. The resulting fusion yields the land-based water conservancy-specific topographic data DEM5, expressed as:

[0144]

[0145] The obtained terrestrial topographic data DEM5 is fused with the underwater DEM, and the expression is:

[0146]

[0147] An adaptive weighted least squares method is used for smoothing, and unreasonable terrain and incomplete building features are checked, edited, and modified to generate an HDEM. Figure 2 This is a schematic diagram of the HDEM results for the bridge area. Figure 3 This is a schematic diagram of the HDEM results for the reservoir dam area. Figure 4 This is a schematic diagram of the HDEM results for the sluice gate area.

[0148] 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 extracting water conservancy-specific terrain based on lidar and satellite remote sensing, characterized in that, Includes the following steps: Collect satellite radar data and satellite remote sensing data of a preset area, and preprocess the radar data and the remote sensing data; The radar data and the remote sensing data are compared and the regions are divided. The consistent regions are taken as the real regions and the inconsistent regions are taken as the fuzzy regions. The water conservancy boundaries and terrain are extracted from the real regions to obtain the first data. Collaborative supplementary surveys are performed on the ambiguous area to obtain covered water conservancy data. Closed contour lines are used to extract water conservancy boundaries and terrain from the covered water conservancy data to obtain second data. The second data includes the terrain of the second region and the second water conservancy boundary. Based on the first data and the second data, a time-series hydraulic topology-specific terrain extraction model is constructed to extract hydraulic topography data DEM1 for a preset area. A two-dimensional hydrodynamic model is constructed based on the pre-defined terrain data and hydraulic boundaries to simulate flood inundation, and the water-affected areas are extracted according to the inundation range. For details of water-related structures, LiDAR is used to acquire point cloud data of the water-related structure area. According to the functional characteristics, it is divided into water-blocking structures and water-passing structures. Graphical processing and a combination of addition and removal operations are performed, and it is fused with the underwater DEM. Weighted correction and Gaussian smoothing algorithms are used to embed it into the preset area water and land topographic data DEM1 to generate the preset area integrated water and land HDEM data.

2. The method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing according to claim 1, characterized in that, The method for determining region division based on comparison of the radar data and the remote sensing data includes: Edge features are extracted from optical images obtained from radar and satellite remote sensing data. Sensor registration is performed by fusing these edge features. Temporal weights are introduced to perform multi-temporal matching on the spatially registered radar and satellite remote sensing data. The expression is as follows: The time series weight at time point t is: Adjusted for seasonal variations, the overall similarity score is: The gradient component of the template image in the x-direction at point v is: The gradient component of the search image in the x-direction at point v is: The gradient component of the template image in the y-direction at point v is: The gradient component of the image in the y-direction at point v is: The total number of features is The number of timing templates is ; Radar and satellite remote sensing data are grouped, with those having a matching degree greater than 0.673 grouped together. Based on the registered multi-source data, the spatiotemporal consistency of different groups is calculated: Where spatiotemporal consistency is The satellite remote sensing derived point cloud at time t is The radar point cloud at time t is The upper limit of observation time is The point cloud structure similarity is The point cloud is projected onto a 2.5D grid, the mean elevation and variance of intensity within the grid are calculated, and the SSIM comparison term is obtained based on the mean elevation and variance of intensity. The temporal variance is obtained by summing the temporal standard deviations of spatiotemporal consistency. When the pixel-level consistency is greater than or equal to 0.786 and the temporal variance is less than 0.352, the elevation is stable, the intensity-vegetation index is strongly correlated, and the location of the group is the actual area. When the pixel-level consistency is less than 0.786 and the temporal variance is greater than 0.352, the elevation changes abruptly, the intensity-vegetation index is weakly correlated, and the location of the group is the ambiguous area.

3. The method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing according to claim 1, characterized in that, The method for extracting hydraulic boundaries and topography of the actual area includes: The terrain of the first region is obtained by extracting terrain features from the actual area using SAR technology. Potential water bodies in the actual area are located by locating the slope threshold. The NDWI index is binarized in the potential water body area using the maximum inter-class variance method to extract the initial water body boundary. For the cloud-covered area of ​​the optical image, the water backscattering coefficient of the SAR image is used to fill the boundary missing. The two are then merged to generate the complete water body outline of the actual area, which is used as the first water boundary. The first water boundary and the terrain of the first region are output as the first data.

4. The method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing according to claim 1, characterized in that, A method for obtaining topographic and overlay water conservancy data of a second region by collaborative supplementary surveying of the aforementioned ambiguous region includes: Construct an integrated "air-space-ground" observation network to achieve spatiotemporal data complementarity, acquire medium-density lidar point cloud data, radar back reflection coefficient and optical vegetation index, identify fuzzy region types through random forest model, and obtain high-priority regions, medium-priority regions and low-priority regions. For high-priority areas, UAV data is introduced. The UAV altitude is equal to the maximum height of vegetation plus a safety margin. The overlap rate is proportional to the vegetation complexity. A mesh flight path is generated based on the LiDAR scanning angle. The drone scans the same area from different directions, and the data is fused using a point cloud registration algorithm. The expression is: Where the rotation matrix is ​​A, the translation vector is w, and the u-th point in the source point cloud is... The source point cloud consists of point cloud data acquired from different scanning angles, which are to be registered to the target point cloud. The u-th point in the target point cloud is... The directional weighting coefficient is The normal vector is The merged data is G; Spatiotemporal registration is performed on multi-temporal radar images and laser data of the same area. The real terrain is estimated using Kalman filtering and used as the second-region terrain data for the blurred area. For multi-echo lidar point cloud areas, the first echo is separated, the echo intensity difference is calculated, and potential river centerlines are identified by searching for local minimum points in areas with high echo intensity differences. For multi-source feature fusion areas, linear shadow features are extracted from optical images, and anomalous stripes of backscattering coefficients are analyzed from radar images. Combined with water conservancy planning data, spatial overlay analysis is used to locate the location information of water-related structures. The combined second-region terrain data, potential river centerlines, and location information of water-related structures are output as comprehensive water conservancy data.

5. The method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing according to claim 1, characterized in that, The method for extracting the second data includes: Based on the covered water conservancy data, the elevation gradient and flow direction of the river centerline neighborhood are obtained, and the robust fuzzy local information C-means clustering algorithm is used to initially separate water bodies and non-water bodies. Given the objective function of fuzzy clustering, the expression is: Where the objective function is The prior probability that pixel a belongs to cluster c is The terrain dissimilarity measure between pixel a and cluster c is: The fuzzy local variables of pixel a and cluster c are: The ridge regression coefficient is The total number of pixels is The membership degree of pixel a to cluster c is The spectral dissimilarity measure between pixel a and cluster c is The balance parameter is Balancing spectral and topographic features, the elevation variation coefficient of pixel a is The elevation variation coefficient of cluster c is ; Calculate the mean, covariance, and dissimilarity measures: The mean of cluster c is covariance is The transpose of the matrix is The spectral measure of pixel a is Data dimensions are The mean of pixel a in cluster c is ; Calculate the local variation density and fuzzy local variables: The local variation density of pixel a is The neighboring pixels of pixel a are The set of neighboring pixels is The mean of the spectral measurements of the neighboring pixels centered at pixel a is The variance of the spectral test of the neighborhood pixels centered at pixel a is Neighboring pixels The local variation density is ; Calculate prior probability and membership degree: The index clustering is as follows: The domain effect strength is The pixel a attribute is The number of clusters is The prior probability that pixel a belongs to cluster i is The spectral dissimilarity measure between pixel a and cluster i is The fuzzy local variables of pixel a and cluster i are ; If the maximum error of the objective function from time point t to time point t+1 is less than the error threshold, then the membership matrix is ​​defuzzified to obtain the binary segmentation result. Morphological closing operation is used to fill the small holes in the binary segmentation result and extract the contour of the largest connected component. Apply hydraulic set constraints to the contour of the largest connected region, forcing the hydraulic cluster centers to satisfy... Isolated clusters with areas smaller than a threshold are removed to obtain segmented data; the neighborhood elevation gradient is... The direction of water flow is ; The gradient and gradient magnitude of the fuzzy water area are obtained from the segmented data. Anisotropic diffusion filtering is applied to the segmented data to preserve the sharpness of the channel boundaries. The expression is as follows: The edge stopping function is: To control the diffusion intensity, the edge sensitivity coefficient in anisotropic diffusion filtering is... It adaptively adjusts based on the boundary gradient, with a gradient magnitude of s and an image intensity gradient of... ; The segmented data boundary is used as the second water conservancy boundary data, and combined with the topographic data of the second region, the output is the second data.

6. The method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing according to claim 1, characterized in that, A method for constructing a time-series hydraulic topology-specific terrain extraction model based on the first data and the second data includes: A basic terrain grid is constructed based on the terrain data of the first and second datasets. Static parameters are extracted, and temporal difference analysis is performed on the hydraulic boundaries of the first and second datasets to identify boundary change areas. Vector-raster conversion is performed on the dynamic boundaries to generate a temporal raster layer and extract the dynamic topology node set. Static parameters include terrain elevation, slope, and aspect. Dynamic nodes are matched with static nodes, and a spatiotemporal joint node set is generated through spatial overlay analysis. A node attribute table is established to record the spatial coordinates, temporal attributes, and hydraulic attributes of the nodes. Based on the spatiotemporal joint node set, a dynamic topological edge set is constructed using the triangulation algorithm of the point set. Here, the adjacency relationship represents the direct connectivity between nodes, the temporal association represents the state change of the same node at different time steps, and the edge attributes include distance, flow direction, and water passage capacity. Distance, water level difference, and time step are used as time-series weights for dynamic topological edges to construct a time-series hydraulic topology-specific terrain extraction model, and output the hydraulic topographic data DEM1 of the preset area.

7. The method for extracting water conservancy-specific terrain based on lidar and satellite remote sensing according to claim 1, characterized in that, The method for obtaining the water-affected area includes: The first and second water conservancy boundaries are merged into a complete water conservancy boundary to form river network data; The pre-defined regional topographic data DEM1 is divided into grids. The area adjacent to the water conservancy boundary and the area covering the water conservancy data are divided into grids with a grid size of 5m×5m, while other areas far from the river channel are divided into grids with a grid size of 20m×20m. Roughness, area coefficient and time step are assigned to the grid after division, and the grid model results are output. Based on historical flow data of the preset area, a flood frequency curve is fitted to determine the design peak flow corresponding to the target frequency, which serves as the inflow boundary condition. Based on the control volume method, the two-dimensional Saint-Venant shallow water equation is used to simulate flood inundation with the maximum design peak flow of the river in the preset area as the boundary condition, and the inundation range is output as the water-affected area.

Images