A reservoir sediment modeling method based on multi-sequential multi-beam sounding point cloud data

Abstract

This invention proposes a reservoir sediment modeling method based on multi-time-series multibeam bathymetry point cloud data, relating to the field of reservoir sediment monitoring and modeling. The method includes: constructing a spatial grid skeleton model of the reservoir area; collecting multi-time-series multibeam bathymetry point cloud data and corresponding hydrological auxiliary data; using an adaptive tile partitioning method for spatial division and a feature point downsampling method to generate multi-level point cloud tiles, establishing a three-dimensional composite index; processing the indexed multi-level point cloud tile set using an anisotropic ICP registration method to obtain a multi-time-series registered point cloud tile set; establishing an initial irregular triangular mesh sediment model, using historical hydrological auxiliary data for differentiated partitioning, and dynamically updating based on the partition status to obtain a dynamic reservoir sediment model; S5, calculating sedimentation parameters and using a three-dimensional engine to map the sedimentation to obtain an interactive three-dimensional dynamic sedimentation status map. This invention achieves high-precision dynamic visualization of the reservoir sedimentation process.

Description

Technical Field

[0001] This invention relates to the field of reservoir sediment monitoring and modeling technology, and in particular to a reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data. Background Technology

[0002] Sedimentation in reservoirs is one of the core issues affecting their long-term operational safety and water storage capacity. Inflowing water carries large amounts of suspended and bedload sediment, which continuously settles and accumulates in the reservoir area after the flow velocity decreases, gradually encroaching on the effective storage capacity and reducing the reservoir's flood control standards and water supply and storage capacity. As reservoirs age, sedimentation becomes increasingly severe. Failure to promptly grasp the spatial distribution and dynamic evolution of sedimentation will pose significant risks to reservoir safety operation and dredging decisions. Therefore, establishing an effective reservoir sedimentation monitoring and modeling system is a crucial technical requirement for reservoir operation and management.

[0003] Multibeam bathymetry, with its ability to acquire full-coverage, high-density underwater topographic data, has gradually become the mainstream method for underwater topographic surveying of reservoirs. Compared with traditional single-beam bathymetry and transect methods, multibeam bathymetry systems can acquire strip-shaped, high-density three-dimensional point cloud data in a single operation, significantly improving the accuracy of reservoir topographic reconstruction. However, the massive amount of multibeam bathymetry point cloud data, the accumulation of measurement errors between multiple time-series data, and the significant dynamic changes in reservoir topography with hydrological processes make efficient organization, accurate registration, and dynamic modeling of multi-time-series multibeam point cloud data a pressing technical challenge in current research and engineering applications.

[0004] Chinese invention patent application number 202411178622.3 discloses a method and system for short-term prediction of power generation of small hydropower stations. This method uses multibeam sonar to perform a full-coverage scan of the reservoir basin, acquiring regional sediment deposition distribution characteristics such as sediment thickness distribution, sediment volume distribution, and sediment morphology parameter distribution. Then, it uses the Delaunay triangular mesh algorithm to reconstruct the three-dimensional surface of the sediment deposits and combines this with a hybrid neural network model to make short-term predictions of the power generation of small hydropower stations. However, it cannot dynamically update the data based on the underwater topography of the reservoir area in a timely manner, leading to a decrease in data measurement accuracy and making it difficult to meet the engineering requirements for fine-grained modeling and dynamic evolution quantitative analysis of sediment deposition distribution in reservoirs. Summary of the Invention

[0005] In view of this, the present invention provides a reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data, in order to solve the problem of low dynamic update efficiency of sediment models caused by irregular terrain in the reservoir area in the existing technology, and realize high-precision dynamic visualization of the sediment deposition process in the reservoir.

[0006] The technical solution of this invention is implemented as follows:

[0007] On the one hand, this invention provides a method for reservoir sediment modeling based on multi-time-series multibeam bathymetry point cloud data, including:

[0008] S1. Construct a spatial grid skeleton model of the reservoir area based on reservoir mapping data, and collect multi-time series multi-beam point cloud data and corresponding hydrological auxiliary data of the reservoir according to a preset cycle;

[0009] S2. Preprocess the multi-time series multi-beam point cloud data to obtain processed point cloud data. Use an adaptive tile partitioning method to spatially divide the processed point cloud data. Use a feature point downsampling method that integrates hydrological scouring and sedimentation differential intensity index to generate multi-level point cloud tiles. Establish a three-dimensional composite index that integrates spatial coding, temporal labels and detail levels to obtain an indexed multi-level point cloud tile set.

[0010] S3. Using anisotropic ICP registration method on a tile-by-tile basis, the indexed multi-level point cloud tile set is processed, and constraints are applied based on hydrological auxiliary data to obtain a multi-temporal registration point cloud tile set.

[0011] S4. Based on the point cloud tile set with multi-time registration, the Delaunay triangulation method is used to establish an initial irregular triangular network sediment model. Historical hydrological auxiliary data is used to differentiate the spatial grid skeleton model of the reservoir area according to the distance of hydrological status to obtain the partition status. Based on the partition status, the initial irregular triangular network sediment model is dynamically updated to obtain a dynamic reservoir sediment model.

[0012] S5. Based on the elevation calculation of the vertices of adjacent temporal triangulation networks in the dynamic reservoir sediment model, the sedimentation parameters are obtained. The 3D engine is used to map the sedimentation parameters to the dynamic reservoir sediment model in the form of color gradients to obtain an interactive 3D dynamic sedimentation status map.

[0013] Based on the above technical solutions, preferably, step S1 specifically includes:

[0014] Acquire reservoir mapping data and historical hydrological auxiliary data, use the highest historical water level of the reservoir as a reference benchmark to determine the static constraint boundary of the reservoir area, and obtain the spatial range of the reservoir area; the hydrological auxiliary data includes water level, flow rate and sediment concentration;

[0015] The reservoir area is spatially divided into multiple grid units according to a preset size using a regular grid partitioning method. All grid units are then stored according to a quadtree coding rule to obtain the reservoir area spatial grid skeleton model.

[0016] Multibeam underwater sounding is performed on the entire reservoir area at a preset cycle to obtain sounding system parameters. These parameters are then associated and stored in the corresponding grid cells of the reservoir area spatial skeleton model, resulting in multi-time series multibeam point cloud data and corresponding hydrological auxiliary data. The sounding system parameters include the beam departure angles corresponding to each sounding point. Measuring water depth And sound velocity profile data.

[0017] Based on the above technical solutions, preferably, the step of using an adaptive tile partitioning method to spatially divide the processed point cloud data specifically includes:

[0018] Using the spatial grid skeleton model of the reservoir area as the initial tile division basis, the point cloud density of the processed point cloud data in each grid cell is statistically analyzed to obtain the point cloud density distribution of each grid cell.

[0019] Based on the point cloud density of each grid cell, a preset high-density threshold is applied. With preset low density threshold Perform hierarchical partitioning, where :

[0020] The point cloud density of the grid cells exceeds the preset high-density threshold. The grid cells are recursively subdivided into quadtrees to obtain the first level of detail tiles;

[0021] Set the point cloud density of the grid cells below a preset low-density threshold. The grid cells are reverse-merged into quadtrees and merged into the parent node tile level to obtain the third detail level tile.

[0022] Set the point cloud density of the grid cells to The grid cells of the interval retain the scale of the current level node to obtain the second level of detail tile;

[0023] An adaptive size set of tiles is formed based on the first level of detail tiles, the second level of detail tiles, and the third level of detail tiles;

[0024] For each tile in the adaptive-size tile set, the width is extended outwards on each of its four sides by [missing value]. The overlapping area is set inside the static constraint boundary of the reservoir area with a width of The boundary buffer is used to obtain a multi-layered tile structure.

[0025] Based on the above technical solutions, preferably, the method of generating multi-level point cloud tiles using the feature point downsampling method that integrates hydrological scouring and silting differential intensity indices specifically includes:

[0026] For the preprocessed point cloud in each tile of the multi-layer tile structure, with each point p as the core, search for k nearest neighbor points in a spherical neighborhood with radius r to construct the k-neighborhood point set of each point.

[0027] Calculate the mean curvature of each point based on the k-neighbor set. and normal rate of change ;

[0028] By querying the time-series dimension of a three-dimensional composite index, we can find the high-order sequence of point p in the first M historical time series. Calculate the time-series siltation differential intensity ;

[0029] The mean curvature, rate of change of normal, and temporal scouring and silting differential intensity of each point were normalized within the tile to which each point belonged, and a weighted summation method was used to calculate the comprehensive feature importance score of each point. The importance scores for each point are obtained.

[0030] Based on the target point cloud density of each detailed level tile, points within each tile are sorted from highest to lowest according to their comprehensive feature importance score, and then sequentially selected and retained. Points with a comprehensive feature importance score below a preset minimum importance threshold are excluded. In flat terrain areas, sampling is performed at preset spatial intervals. Perform uniform rasterization and crop the point cloud of each tile according to the upper limit of the point cloud density of each detail level to obtain multi-level point cloud tiles.

[0031] Based on the above technical solutions, preferably, the formula for calculating the comprehensive feature importance score W(p) is as follows:

[0032] ;

[0033] ;

[0034] in, Indicates the temporal scouring and silting differential intensity. For point In the Elevation values ​​in a historical time series This represents the total number of historical time series preceding the current time series. This is the hydrological similarity distance attenuation coefficient. For the current time series The hydrological state vector, For the first A historical time series of hydrological state vectors The distance represents the similarity of hydrological conditions. This represents the normalized mean curvature. Represents the rate of change of the normalized normal. This represents the normalized temporal scouring and silting difference intensity. The weighting coefficients represent the normalized mean curvature. The weighting coefficients represent the normalized rate of change of the normalized normal. The weighting coefficients represent the normalized temporal sludge differential intensity.

[0035] Based on the above technical solutions, preferably, step S3 specifically includes:

[0036] Using the sounding system parameters in the corresponding grid cells of the reservoir area spatial grid skeleton model, and based on the multibeam acoustic sounding error propagation relationship, each point in each tile of the multi-level point cloud tile set is... Build Diagonal anisotropic error covariance matrix And form a set of anisotropic error covariance matrices;

[0037] ;

[0038] in, The standard deviation of the horizontal position error in the direction perpendicular to the survey line. The standard deviation of the horizontal position error along the survey line direction. This represents the standard deviation of the elevation direction error.

[0039] In each temporal point cloud tile, the tile corresponding to the first acquisition time is used as the reference tile. If the mean of the determinant of the anisotropic error covariance matrix of tiles corresponding to other time times is less than that of the tile corresponding to the first acquisition time, then the tile corresponding to that time time is replaced as the reference tile. For the target tiles corresponding to the remaining time times, the reference tile is used as the registration reference, and the rotation matrix is ​​solved based on the anisotropic error covariance matrix and the objective function. With translation vector The process continues until the change in the objective function between two consecutive iterations is lower than a preset convergence threshold, at which point the converged rotation matrix is ​​output. With translation vector :

[0040] ;

[0041] in, This represents the three-dimensional coordinate vector of the i-th point in the target tile. Indicates the reference tile and The corresponding three-dimensional coordinate vector of the nearest neighbor matching point, Represents the i-th point in the target tile. The anisotropic error covariance matrix, Represents the i-th corresponding point in the reference tile. The anisotropic error covariance matrix;

[0042] Based on the converged rotation matrix With translation vector The target temporal point cloud tiles are transformed to obtain spatially aligned temporal tiles.

[0043] The spatially aligned point cloud tiles are fused by weighted average, and the spatial difference features of each point cloud are preserved based on the temporal labels in the 3D composite index, resulting in a set of point cloud tiles with multi-temporal registration.

[0044] Based on the above technical solutions, preferably, the step of using historical hydrological auxiliary data to differentiate the spatial grid skeleton model of the reservoir area according to the distance of hydrological status specifically includes:

[0045] Obtain hydrological auxiliary data for the current time series, including the current time series flow rate. Sand content and water level And construct the current time series hydrological state vector. Based on the historical hydrological auxiliary data of each grid cell in the reservoir area spatial grid skeleton model, the current time-series hydrological state vector is calculated. Historical hydrological state vectors at various time series Hydrological distance between :

[0046] ;

[0047] in, , , They are respectively the historical number The flow rate, sediment concentration, and water level at each time series, , , These represent the statistical standard deviations of flow, sediment concentration, and water level over all historical time series.

[0048] Will Below the preset similarity threshold Historical time series are identified as historically similar time series with similar hydrological conditions to the current time series, resulting in a set of historically similar time series. ;

[0049] In the statistical reservoir area spatial grid skeleton model, each tile is in Mean value of point cloud elevation variation across historical time series and compared with the preset elevation change threshold. Compare:

[0050] when Exceeding the preset elevation change threshold If the corresponding tile is identified as a tile with significant siltation changes, i.e., a highly active area, then the elevation change confirmation threshold is [value missing]. The area formed by extending a ring of adjacent tiles outward from the spatial distribution range of the high-activity area in the previous time series is defined as the boundary extension zone. The elevation change confirmation threshold of the boundary extension zone is used. ,in ;

[0051] when Below the preset elevation change threshold If the elevation change is minimal, the corresponding tile is classified as a low-activity zone, and the elevation change confirmation threshold is [value missing]. ;

[0052] Write the partition status of each tile and the corresponding elevation change confirmation threshold into the tile status flag of the three-dimensional composite index to obtain the partition status including high-activity area, low-activity area and boundary extension area.

[0053] Based on the above technical solutions, preferably, the dynamic updating of the initial irregular triangular mesh sediment model based on the partition status specifically includes:

[0054] The planar coordinates of the newly added point cloud data in the current time series are matched with the quadtree spatial index in the three-dimensional composite index to find the tile spatial code in which the newly added point cloud data falls, determine the set of tiles affected by the newly added data, and read the tile status flag in the three-dimensional composite index of each tile to obtain the partition status of each tile and the corresponding elevation change confirmation threshold, thus obtaining the set of tiles to be updated in the current time series and the elevation change confirmation threshold of each tile.

[0055] For each tile in the updated tile set, the elevation change of each point in the newly added point cloud data is compared with the elevation change confirmation threshold corresponding to that tile. The subset of point cloud data whose elevation change exceeds the corresponding elevation change confirmation threshold is retained as the current time series valid change point cloud subset, and the point cloud data whose elevation change does not exceed the corresponding elevation change confirmation threshold is determined as an invalid change subset.

[0056] The effective variation point cloud subsets within each tile are integrated into the existing point cloud of the corresponding tile and recursively split into left and right subsets until the number of point clouds in the subsets does not exceed the preset splitting threshold. Delaunay triangulation is performed independently on each subset to reconstruct the local irregular triangular mesh of the tile.

[0057] The initial irregular triangular mesh is reused for invalid variation subsets within each tile;

[0058] The reconstructed local irregular triangular mesh of each tile is merged with the reused initial irregular triangular mesh, and the current time series label is written into the time series dimension of the three-dimensional composite index of the corresponding tile to obtain the dynamic reservoir sediment model.

[0059] Based on the above technical solutions, preferably, step S5 specifically includes:

[0060] Based on the vertices of the triangular network corresponding to adjacent time series in the dynamic reservoir sediment model, the sedimentation parameters at each vertex are calculated, including sedimentation thickness. siltation rate and the amount of sediment in the area where each triangular facet is located. :

[0061] ;

[0062] ;

[0063] ;

[0064] in, The accumulation thickness between adjacent time sequences at the vertices of the triangular mesh. For the current time series The corresponding elevation value of the vertex The elevation value of the corresponding vertex in the adjacent preceding time series. The accumulation rate at that vertex. Adjacent time series and The time interval between The amount of sediment in the target area, For the first The area of ​​each triangular facet. For the first The thickness of the sediment at the three vertices of the triangular facet The mean;

[0065] By loading irregular triangular mesh model tiles corresponding to each time series and detail level of the dynamic reservoir sediment model using a 3D engine, and mapping each sedimentation parameter onto the surface of the dynamic reservoir sediment model in the form of color gradient, an interactive 3D dynamic sedimentation status map is obtained.

[0066] More preferably, the reservoir sediment modeling method further includes:

[0067] A time-series prediction model was constructed using a long short-term memory neural network. The model was trained by taking the historical sedimentation parameter sequences and corresponding hydrological auxiliary data sequences as inputs.

[0068] By inputting preset future working condition scenario parameters into the trained time series prediction model, the predicted results of sediment deposition volume, deposition range, and elevation changes within a specified time period under that working condition scenario are obtained. The preset future working condition scenario parameters are preset water level values ​​for a specified future time period set according to engineering requirements. Preset flow rate and preset values ​​of sand content .

[0069] The present invention has the following advantages over the prior art:

[0070] (1) By constructing a spatial grid skeleton model of the reservoir area, combined with anisotropic ICP registration, differentiated partitioning and Delaunay triangulation dynamic update, the limitations of the existing technology of static, low efficiency and disconnect from application are effectively overcome, and high-precision dynamic visualization and parameter quantification of the reservoir sediment deposition process are realized. The accuracy, efficiency and engineering practicality of sediment dynamic monitoring are significantly improved, and intuitive and reliable scientific decision-making basis is provided for reservoir operation management and dredging planning.

[0071] (2) The characteristic point downsampling method of integrating hydrological scour and deposition differential intensity index is adopted to effectively compress the data scale while taking into account the topographic geometric features and the spatial distribution of historical scour and deposition active areas, thereby reducing computational redundancy;

[0072] (3) Based on the propagation relationship of multi-beam acoustic sounding error, an anisotropic error covariance matrix is ​​constructed for each point in the tile, and the objective function is used to solve it to improve the spatial alignment accuracy of multi-time series point clouds.

[0073] (4) Construct the current time series hydrological state vector using historical hydrological auxiliary data, filter the historical similar time series set by hydrological state distance, and perform differentiated partitioning and threshold configuration based on the statistical results of the elevation change of each tile, so that the dynamic update process maintains a high sensitivity to change capture of active siltation areas, and reduces invalid local triangular network reconstruction calculations triggered by noise while ensuring the timeliness of the dynamic reservoir sediment model. Attached Figure Description

[0074] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0075] Figure 1 This is a flowchart of a reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data according to the present invention;

[0076] Figure 2 This is a block diagram of a reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data according to the present invention;

[0077] Figure 3This is a flowchart of the ICP registration process for a reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data according to the present invention. Detailed Implementation

[0078] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0079] like Figure 1 and Figure 2 As shown, this invention provides a reservoir sediment modeling method based on multi-time-series multibeam bathymetry point cloud data, comprising:

[0080] S1. Construct a spatial grid skeleton model of the reservoir area based on reservoir mapping data, and collect multi-time series multi-beam point cloud data and corresponding hydrological auxiliary data of the reservoir according to a preset cycle;

[0081] Specifically, step S1 includes:

[0082] Acquire reservoir mapping data and historical hydrological auxiliary data, use the highest historical water level of the reservoir as a reference benchmark to determine the static constraint boundary of the reservoir area, and obtain the spatial range of the reservoir area; the hydrological auxiliary data includes water level, flow rate and sediment concentration;

[0083] The reservoir area is spatially divided into multiple grid units according to a preset size using a regular grid partitioning method. All grid units are then stored according to a quadtree coding rule to obtain the reservoir area spatial grid skeleton model.

[0084] Multibeam underwater sounding is performed on the entire reservoir area at a preset cycle to obtain sounding system parameters. These parameters are then associated and stored in the corresponding grid cells of the reservoir area spatial skeleton model, resulting in multi-time series multibeam point cloud data and corresponding hydrological auxiliary data. The sounding system parameters include the beam departure angles corresponding to each sounding point. Measuring water depth And sound velocity profile data.

[0085] In one embodiment of the present invention, within the spatial range of the reservoir area, a regular grid partitioning method is used to spatially partition the area according to a preset size, resulting in multiple grid cells. The grid partitioning size is no larger than 10m × 10m × 10m, in order to balance computational efficiency and spatial resolution. The point cloud computing and interpolation resolution is less than 0.1m, the kriging interpolation resolution is 0.1~0.5m, the noise removal rate is ≥98%, and the temporal consistency error is ≤0.03m, ensuring the integrity of the single temporal point cloud and the adaptability to multiple temporal sequences.

[0086] In one embodiment of the present invention, a multibeam echo sounder system is used to conduct underwater echo sounding operations on the entire reservoir area or key siltation areas according to a preset time series period (such as monthly, quarterly, annual, or rainy season) to acquire multi-time series point cloud data. The point cloud data includes three-dimensional coordinates (X, Y, Z), measurement timestamps, sediment particle size correlation information, and echo sounder system parameters. The echo sounder system parameters include the beam number and beam exit angle corresponding to each echo point. and measuring water depth During the data collection process, hydrological auxiliary data such as reservoir water level, flow rate, and sediment concentration were recorded simultaneously to provide data support for subsequent sedimentation analysis and prediction. The time-series bathymetry data were correlated and stored in the corresponding grid cells of the reservoir area spatial grid skeleton model to obtain multi-time-series multibeam point cloud data and corresponding hydrological auxiliary data.

[0087] S2. Preprocess the multi-time series multi-beam point cloud data to obtain processed point cloud data. Use an adaptive tile partitioning method to spatially divide the processed point cloud data. Use a feature point downsampling method that integrates hydrological scouring and sedimentation differential intensity index to generate multi-level point cloud tiles. Establish a three-dimensional composite index that integrates spatial coding, temporal labels and detail levels to obtain an indexed multi-level point cloud tile set.

[0088] In one embodiment of the present invention, preprocessing includes noise filtering, outlier removal, and data completion to ensure that the integrity and accuracy of a single temporal point cloud meet the requirements of subsequent modeling.

[0089] In one embodiment of the present invention, the step of spatially partitioning the processed point cloud data using an adaptive tile partitioning method specifically includes:

[0090] Using the spatial grid skeleton model of the reservoir area as the initial tile division basis, the point cloud density of the processed point cloud data in each grid cell is statistically analyzed to obtain the point cloud density distribution of each grid cell.

[0091] Based on the point cloud density of each grid cell, a preset high-density threshold is applied. With preset low density threshold Perform hierarchical partitioning, where :

[0092] The point cloud density of the grid cells exceeds the preset high-density threshold. The grid cells are recursively subdivided into quadtrees to obtain the first level of detail tiles;

[0093] Set the point cloud density of the grid cells below a preset low-density threshold. The grid cells are reverse-merged into quadtrees and merged into the parent node tile level to obtain the third detail level tile.

[0094] Set the point cloud density of the grid cells to The grid cells of the interval retain the scale of the current level node to obtain the second level of detail tile;

[0095] An adaptive size set of tiles is formed based on the first level of detail tiles, the second level of detail tiles, and the third level of detail tiles;

[0096] For each tile in the adaptive-size tile set, the width is extended outwards on each of its four sides by [missing value]. The overlapping area is set inside the static constraint boundary of the reservoir area with a width of The boundary buffer is used to obtain a multi-layered tile structure.

[0097] In one embodiment of the present invention, a feature point downsampling method integrating hydrological scouring and sedimentation differential intensity indices is used to generate multi-level point cloud tiles, specifically including:

[0098] For the preprocessed point cloud in each tile of the multi-layer tile structure, with each point p as the core, search for k nearest neighbor points in a spherical neighborhood with radius r to construct the k-neighborhood point set of each point.

[0099] Calculate the mean curvature of each point based on the k-neighbor set. and normal rate of change ;

[0100] ;

[0101] in, , These represent the principal curvatures in the two principal directions at that point;

[0102] By querying the time-series dimension of a three-dimensional composite index, we can find the high-order sequence of point p in the first M historical time series. Calculate the time-series siltation differential intensity ;

[0103] The mean curvature, rate of change of normal, and temporal scouring and silting differential intensity of each point were normalized within the tile to which each point belonged, and a weighted summation method was used to calculate the comprehensive feature importance score of each point. The importance scores for each point are obtained.

[0104] Based on the target point cloud density of each detailed level tile, points within each tile are sorted from highest to lowest according to their comprehensive feature importance score, and then sequentially selected and retained. Points with a comprehensive feature importance score below a preset minimum importance threshold are excluded. In flat terrain areas, sampling is performed at preset spatial intervals. Perform uniform rasterization and crop the point cloud of each tile according to the upper limit of the point cloud density of each detail level to obtain multi-level point cloud tiles.

[0105] Furthermore, the formula for calculating the comprehensive feature importance score W(p) is as follows:

[0106] ;

[0107] ;

[0108] in, Indicates the temporal scouring and silting differential intensity. For point In the Elevation values ​​in a historical time series This represents the total number of historical time series preceding the current time series. This is the hydrological similarity distance attenuation coefficient. For the current time series The hydrological state vector, For the first A historical time series of hydrological state vectors The distance represents the similarity of hydrological conditions. This represents the normalized mean curvature. Represents the rate of change of the normalized normal. This represents the normalized temporal scouring and silting difference intensity. The weighting coefficients represent the normalized mean curvature. The weighting coefficients represent the normalized rate of change of the normalized normal. The weighting coefficients represent the normalized temporal sludge differential intensity.

[0109] This invention employs a feature point downsampling method that integrates hydrological scour and deposition differential intensity indices. This method effectively compresses the data size while taking into account the topographic geometric features and the spatial distribution of historically active scour and deposition areas, thereby reducing computational redundancy.

[0110] In one embodiment of the present invention, a three-dimensional composite index integrating spatial coding, temporal tags, and detail levels is established, specifically including:

[0111] In the spatial grid skeleton model of the reservoir area, the point cloud is adaptively tiled using a quadtree pyramid partitioning algorithm, and the spatial scale of the tiles is dynamically adjusted according to the point cloud density.

[0112] A feature point downsampling method based on comprehensive feature importance scores is used to generate point cloud tiles with multiple detail levels;

[0113] Point cloud data from different time series use the same tile size and are distinguished by time series labels;

[0114] A three-dimensional composite index integrating spatial coding, temporal labels, and detail level identifiers is constructed, and tile data is stored in a compressed format.

[0115] In one embodiment of the present invention, the initial tile side length is set to 50m×50m, and the tile is adaptively subdivided according to the point cloud density within a single tile. When the point cloud density within a tile exceeds a preset threshold (e.g., 100 points / m), the tile is further subdivided. 2 For areas requiring high precision in multibeam bathymetry, the system automatically subdivides the data down to sub-tiles (minimum side length 10m×10m). In areas with gentle sedimentation changes and sparse point clouds (such as deep water areas in front of dams), the parent tile scale is preserved to avoid redundant computational overhead. For detailed dimensions, feature point downsampling technology is used to generate multi-level point cloud tiles, adapting to different visualization and analysis needs.

[0116] Specifically, the point cloud density of LOD0 (first level of detail tile) is ≥100 points / m. 2 The original high-precision point cloud core features are preserved for detailed depiction of siltation and accurate parameter calculation; LOD1 (second level of detail tile) is downsampled to 50 points / m through feature point downsampling. 2 Key topographic features are preserved for routine sedimentation analysis; LOD2 (third level of detail tiles) are downsampled to 20 points / m. 2 Used for rapid rendering of the entire reservoir area and observation of overall trends.

[0117] In one embodiment of the present invention, the k value is dynamically adapted according to the tile level (k=10 for the first detail level tile, k=8 for the second detail level tile, and k=5 for the third detail level tile), and the neighborhood search radius is consistent with the radius filtering radius (0.5m) to ensure that the neighborhood point set can reflect the local terrain undulations.

[0118] In one embodiment of this invention, the S3M format is used to store the data of each tile at different detail levels, while simultaneously converting it to the Cesium-compatible 3DTiles format (adapting to Cesium's spatiotemporal data loading standard). The Draco compression algorithm is used to perform differentiated compression of vertex coordinates at different levels—the compression error for LOD0 level is ≤0.005m, and the compression error for LOD1~LOD3 levels is ≤0.01m, ensuring accuracy requirements at each level while reducing storage capacity. A three-dimensional composite index is established, integrating spatial encoding, temporal tags, and detail levels. The three-dimensional composite index uses quadtree encoding to quickly locate the spatial position of tiles, the temporal index associates multiple temporal data of the same spatial tile, and the detail level index supports the Cesium engine in dynamically calling tiles at different detail levels as needed (e.g., calling the first detail level tile for detailed analysis, and calling the third detail level tile for global preview). This embodiment adds a temporal index dimension to the tiled data structure, adapting to Cesium's water conservancy visualization scenarios and improving multi-terminal lightweight display capabilities.

[0119] S3. Using anisotropic ICP registration method on a tile-by-tile basis, the indexed multi-level point cloud tile set is processed, and constraints are applied based on hydrological auxiliary data to obtain a multi-temporal registration point cloud tile set.

[0120] like Figure 3 As shown, specifically, step S3 includes:

[0121] Using the sounding system parameters in the corresponding grid cells of the reservoir area spatial grid skeleton model, and based on the multibeam acoustic sounding error propagation relationship, each point in each tile of the multi-level point cloud tile set is... Build Diagonal anisotropic error covariance matrix And form a set of anisotropic error covariance matrices:

[0122] ;

[0123] in, The standard deviation of the horizontal position error in the direction perpendicular to the survey line. The standard deviation of the horizontal position error along the survey line direction. This represents the standard deviation of the elevation direction error.

[0124] In each temporal point cloud tile, the tile corresponding to the first acquisition time is used as the reference tile. If the mean of the determinant of the anisotropic error covariance matrix of tiles corresponding to other time times is less than that of the tile corresponding to the first acquisition time, then the tile corresponding to that time time is replaced as the reference tile. For the target tiles corresponding to the remaining time times, the reference tile is used as the registration reference, and the rotation matrix is ​​solved based on the anisotropic error covariance matrix and the objective function. With translation vector The process continues until the change in the objective function between two consecutive iterations is lower than a preset convergence threshold, at which point the converged rotation matrix is ​​output. With translation vector :

[0125] ;

[0126] in, This represents the three-dimensional coordinate vector of the i-th point in the target tile. Indicates the reference tile and The corresponding three-dimensional coordinate vector of the nearest neighbor matching point, Represents the i-th point in the target tile. The anisotropic error covariance matrix, Represents the i-th corresponding point in the reference tile. The anisotropic error covariance matrix;

[0127] Based on the converged rotation matrix With translation vector The target temporal point cloud tiles are transformed to obtain spatially aligned temporal tiles.

[0128] The spatially aligned point cloud tiles are fused by weighted average, and the spatial difference features of each point cloud are preserved based on the temporal labels in the 3D composite index, resulting in a set of point cloud tiles with multi-temporal registration.

[0129] Understandable. It is the beamout angle. With measuring water depth The function, as Increased significantly; The accuracy of GNSS / INS integrated navigation is determined by the fact that it is approximately constant within a single tile. It measures water depth. The function of the sound velocity profile correction residual reflects the elevation uncertainty caused by the curvature of the sound ray.

[0130] In the iterative solution of anisotropic ICP, the registration point pair established in the current iteration step is used as the basis. As the basic unit of computation, where For a point in the target tile, For the nearest point in the reference tile, the inverse of the joint anisotropic error covariance matrix is ​​obtained by considering the points. The objective function is constructed using the Mahalanobis distance between the registered point pairs as the accuracy-weighted matrix.

[0131] Based on the propagation relationship of multi-beam acoustic sounding errors, an anisotropic error covariance matrix is ​​constructed for each point within the tile, and solved using an objective function to improve the spatial alignment accuracy of multi-temporal point clouds.

[0132] In one embodiment of the present invention, the weighted average fusion includes: constructing a multi-temporal point cloud fusion dataset, labeling the temporal attributes of each point cloud data according to the time dimension; for multi-temporal point clouds at the same spatial location, using a weighted average method to fuse coordinate data, preserving temporal difference characteristics, and providing accurate data support for dynamic congestion analysis. To improve processing speed, a fine grid with a spacing of 0.1M is further used within the tile, projecting the point cloud onto the grid node constraint framework, and performing point search and registration according to the node framework to improve registration speed.

[0133] In one embodiment of the present invention, the horizontal error of the registration using a subdivided mesh model is ≤ ±0.05m, the elevation error is ≤ ±0.05m, and the single-tile registration processing time is ≤ 5ms.

[0134] S4. Based on the point cloud tile set with multi-time registration, the Delaunay triangulation method is used to establish an initial irregular triangular network sediment model. Historical hydrological auxiliary data is used to differentiate the spatial grid skeleton model of the reservoir area according to the distance of hydrological status to obtain the partition status. Based on the partition status, the initial irregular triangular network sediment model is dynamically updated to obtain a dynamic reservoir sediment model.

[0135] In one embodiment of the present invention, an initial irregular triangular mesh sediment model is established based on a point cloud tile set with multi-temporal registration using the Delaunay triangulation method, specifically including:

[0136] Constructing the outer contour of tile point cloud based on Alpha Shape algorithm;

[0137] During the construction of the triangulation network, a complete tile-level irregular model is generated by recursively splitting the point cloud subset, independently constructing sub-triangulation networks and merging them.

[0138] Triangular mesh removal is performed based on the outer contour to obtain the final irregular model, and the triangle index number is stored.

[0139] In one embodiment of the present invention, the spatial grid skeleton model of the reservoir area is differentiated and partitioned based on hydrological state distance using historical hydrological auxiliary data, specifically including:

[0140] Obtain hydrological auxiliary data for the current time series, including the current time series flow rate. Sand content and water level And construct the current time series hydrological state vector. Based on the historical hydrological auxiliary data of each grid cell in the reservoir area spatial grid skeleton model, the current time-series hydrological state vector is calculated. Historical hydrological state vectors at various time series Hydrological distance between :

[0141] ;

[0142] in, , , They are respectively the historical number The flow rate, sediment concentration, and water level at each time series, , , These represent the statistical standard deviations of flow, sediment concentration, and water level over all historical time series.

[0143] Will Below the preset similarity threshold Historical time series are identified as historically similar time series with similar hydrological conditions to the current time series, resulting in a set of historically similar time series. ;

[0144] In the statistical reservoir area spatial grid skeleton model, each tile is in Mean value of point cloud elevation variation across historical time series and compared with the preset elevation change threshold. Compare:

[0145] when Exceeding the preset elevation change threshold If the corresponding tile is identified as a tile with significant siltation changes, i.e., a highly active area, then the elevation change confirmation threshold is [value missing]. The area formed by extending a ring of adjacent tiles outward from the spatial distribution range of the high-activity area in the previous time series is defined as the boundary extension zone. The elevation change confirmation threshold of the boundary extension zone is used. ,in ;

[0146] when Below the preset elevation change threshold If the elevation change is minimal, the corresponding tile is classified as a low-activity zone, and the elevation change confirmation threshold is [value missing]. ;

[0147] Write the partition status of each tile and the corresponding elevation change confirmation threshold into the tile status flag of the three-dimensional composite index to obtain the partition status including high-activity area, low-activity area and boundary extension area.

[0148] Understandably, a preset threshold for elevation change is used. The threshold for confirming elevation changes is: , as well as All settings can be configured according to actual usage, and this invention does not impose specific limitations on them.

[0149] Furthermore, the initial irregular triangular mesh sediment model is dynamically updated based on the partition status, specifically including:

[0150] The planar coordinates of the newly added point cloud data in the current time series are matched with the quadtree spatial index in the three-dimensional composite index to find the tile spatial code in which the newly added point cloud data falls, determine the set of tiles affected by the newly added data, and read the tile status flag in the three-dimensional composite index of each tile to obtain the partition status of each tile and the corresponding elevation change confirmation threshold, thus obtaining the set of tiles to be updated in the current time series and the elevation change confirmation threshold of each tile.

[0151] For each tile in the updated tile set, the elevation change of each point in the newly added point cloud data is compared with the elevation change confirmation threshold corresponding to that tile. The subset of point cloud data whose elevation change exceeds the corresponding elevation change confirmation threshold is retained as the current time series valid change point cloud subset, and the point cloud data whose elevation change does not exceed the corresponding elevation change confirmation threshold is determined as an invalid change subset.

[0152] The effective variation point cloud subsets within each tile are integrated into the existing point cloud of the corresponding tile and recursively split into left and right subsets until the number of point clouds in the subsets does not exceed the preset splitting threshold. Delaunay triangulation is performed independently on each subset to reconstruct the local irregular triangular mesh of the tile.

[0153] The initial irregular triangular mesh is reused for invalid variation subsets within each tile;

[0154] The reconstructed local irregular triangular mesh of each tile is merged with the reused initial irregular triangular mesh, and the current time series label is written into the time series dimension of the three-dimensional composite index of the corresponding tile to obtain the dynamic reservoir sediment model.

[0155] This invention utilizes historical hydrological auxiliary data to construct the current time-series hydrological state vector, filters historical similar time-series sets by hydrological state distance, and performs differentiated partitioning and threshold configuration based on the statistical results of the elevation change amplitude of each tile. This ensures a high sensitivity to change capture in active sedimentation areas during dynamic updates, thereby reducing invalid local triangulation reconstruction calculations triggered by noise while ensuring the timeliness of the dynamic reservoir sediment model.

[0156] S5. Based on the elevation calculation of the vertices of adjacent temporal triangulation networks in the dynamic reservoir sediment model, the sedimentation parameters are obtained. The 3D engine is used to map the sedimentation parameters to the dynamic reservoir sediment model in the form of color gradients to obtain an interactive 3D dynamic sedimentation status map.

[0157] Specifically, step S5 includes:

[0158] Based on the vertices of the triangular network corresponding to adjacent time series in the dynamic reservoir sediment model, the sedimentation parameters at each vertex are calculated, including sedimentation thickness. siltation rate and the amount of sediment in the area where each triangular facet is located. :

[0159] ;

[0160] ;

[0161] ;

[0162] in, The accumulation thickness between adjacent time sequences at the vertices of the triangular mesh. For the current time series The corresponding elevation value of the vertex The elevation value of the corresponding vertex in the adjacent preceding time series. The accumulation rate at that vertex. Adjacent time series and The time interval between The amount of sediment in the target area, For the first The area of ​​each triangular facet. For the first The thickness of the sediment at the three vertices of the triangular facet The mean;

[0163] By loading irregular triangular mesh model tiles corresponding to each time series and detail level of the dynamic reservoir sediment model using a 3D engine, and mapping each sedimentation parameter onto the surface of the dynamic reservoir sediment model in the form of color gradient, an interactive 3D dynamic sedimentation status map is obtained.

[0164] This invention, by constructing a spatial grid skeleton model of the reservoir area and combining anisotropic ICP registration, differentiated zoning, and Delaunay triangulation dynamic updates, effectively overcomes the limitations of existing technologies, such as static nature, low efficiency, and disconnect from application. It achieves high-precision dynamic visualization and parameter quantification of the reservoir sedimentation process, significantly improving the accuracy, efficiency, and engineering practicality of sediment dynamic monitoring, and providing an intuitive and reliable scientific decision-making basis for reservoir operation management and dredging planning.

[0165] In one embodiment of the present invention, dynamic sludge visualization is achieved by using the Cesium engine to realize three-dimensional dynamic visualization. Relying on its 3DTiles loading capability, irregular model tiles at various time sequences and detail levels are loaded. Through Cesium's material mapping interface, parameters such as sludge thickness and rate are superimposed on the surface of the three-dimensional model in the form of color gradient (red-yellow-green gradient, red ≥ 0.8m, yellow 0.3~0.8m, green < 0.3m), supporting material transparency adjustment and layered rendering. With the help of Cesium's time axis control, the historical sludge process can be traced back along the time dimension to realize a smooth dynamic demonstration of the migration and expansion of sludge areas.

[0166] Furthermore, the interactive 3D dynamic siltation status map now features a time-series dynamic playback function, adapting to reservoir cross-section linkage analysis scenarios.

[0167] In one embodiment of the present invention, the interactive 3D dynamic siltation status map adds an interactive rapid cross-section analysis function. Based on the Cesium interactive interface, operation is adapted to allow users to draw cross-section lines arbitrarily in the Cesium 3D scene or reservoir plan map (generated through two-point / multi-point positioning and freehand drawing, with cross-section line accuracy ≤0.1m), enabling rapid extraction and temporal evolution analysis of cross-section siltation data. This embodiment optimizes data extraction efficiency and 3D linkage logic, reducing analysis time to within 3 seconds, providing precise support for local engineering planning.

[0168] In one embodiment of the present invention, the interactive 3D dynamic siltation situation map supports fast rendering and spatial query analysis based on Cesium. The real-time modeling response speed of a single frame is <200ms, and the dynamic rendering frame rate is ≥30fps when high-speed caching is enabled. It can achieve smooth switching and dynamic playback of at least 6 time-series siltation models, and supports adjusting the playback speed according to the time step (month / quarter / year), intuitively presenting the year-by-year expansion trend of the siltation area at the tail of the reservoir. It supports the spatial query function built into Cesium. Clicking on any position of the model can display the siltation thickness, rate and time series information of the corresponding area in real time, solving the problem that the existing technology cannot visualize dynamic siltation and interactive query.

[0169] In one embodiment of the present invention, no less than 100 borehole sampling points or underwater robot detection data are used, and the relative error of the siltation parameter calculation is strictly controlled within 5%, which meets the accuracy requirements for reservoir engineering applications.

[0170] In one embodiment of the present invention, the reservoir sediment modeling method further includes:

[0171] A time-series prediction model was constructed using a long short-term memory neural network. The model was trained by taking the historical sedimentation parameter sequences and corresponding hydrological auxiliary data sequences as inputs.

[0172] By inputting preset future working condition scenario parameters into the trained time series prediction model, the predicted results of sediment deposition volume, deposition range, and elevation changes within a specified time period under that working condition scenario are obtained. The preset future working condition scenario parameters are preset water level values ​​for a specified future time period set according to engineering requirements. Preset flow rate and preset values ​​of sand content .

[0173] Understandably, the irregular triangular network is reconstructed using the predicted results as vertex elevation data to obtain a predicted irregular triangular network sediment model. Based on the future sedimentation distribution, sedimentation range expansion trend, and elevation change reflected by the predicted irregular triangular network sediment model, a sedimentation risk assessment report is output for reservoir scheduling optimization, dredging project planning, and reservoir capacity assessment.

[0174] 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 reservoir sediment modeling method based on multi-sequential multi-beam sounding point cloud data, characterized in that, include: S1. Construct a spatial grid skeleton model of the reservoir area based on reservoir mapping data, and collect multi-time series multi-beam point cloud data and corresponding hydrological auxiliary data of the reservoir according to a preset cycle; S2. Preprocess the multi-temporal multi-beam point cloud data to obtain processed point cloud data. Use an adaptive tile partitioning method to spatially divide the processed point cloud data to obtain a multi-level tile structure. Use a feature point downsampling method that integrates hydrological scouring and sedimentation differential intensity index to generate multi-level point cloud tiles. Establish a three-dimensional composite index that integrates spatial coding, temporal labels and detail levels to obtain an indexed multi-level point cloud tile set. The method of generating multi-level point cloud tiles using feature point downsampling with integrated hydrological scouring and silting differential intensity indices specifically includes: For the preprocessed point cloud in each tile of the multi-layer tile structure, with each point p as the core, search for k nearest neighbor points in a spherical neighborhood with radius r to construct the k-neighborhood point set of each point. Calculate the mean curvature of each point based on the k-neighbor set. and normal rate of change ; By querying the time-series dimension of a three-dimensional composite index, we can find the high-order sequence of point p in the first M historical time series. Calculate the time-series siltation differential intensity ; The mean curvature, rate of change of normal, and temporal scouring and silting differential intensity of each point were normalized within the tile to which each point belonged, and a weighted summation method was used to calculate the comprehensive feature importance score of each point. The importance scores for each point are obtained. Based on the target point cloud density of each detailed level tile, points within each tile are sorted from highest to lowest according to their comprehensive feature importance score, and then sequentially selected and retained. Points with a comprehensive feature importance score below a preset minimum importance threshold are excluded. In flat terrain areas, sampling is performed at preset spatial intervals. Perform uniform rasterization and crop the point cloud of each tile according to the upper limit of the point cloud density of each detail level to obtain multi-level point cloud tiles. S3. Using anisotropic ICP registration method on a tile-by-tile basis, the indexed multi-level point cloud tile set is processed, and constraints are applied based on hydrological auxiliary data to obtain a multi-temporal registration point cloud tile set. S4. Based on the point cloud tile set with multi-time registration, the Delaunay triangulation method is used to establish an initial irregular triangular network sediment model. Historical hydrological auxiliary data is used to differentiate the spatial grid skeleton model of the reservoir area according to the distance of hydrological status to obtain the partition status. Based on the partition status, the initial irregular triangular network sediment model is dynamically updated to obtain a dynamic reservoir sediment model. S5. Based on the elevation calculation of the vertices of adjacent temporal triangulation networks in the dynamic reservoir sediment model, the sedimentation parameters are obtained. The 3D engine is used to map the sedimentation parameters to the dynamic reservoir sediment model in the form of color gradients to obtain an interactive 3D dynamic sedimentation status map.

2. The reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 1, characterized in that: Step S1 specifically includes: Acquire reservoir mapping data and historical hydrological auxiliary data, use the highest historical water level of the reservoir as a reference benchmark to determine the static constraint boundary of the reservoir area, and obtain the spatial range of the reservoir area; the hydrological auxiliary data includes water level, flow rate and sediment concentration; The reservoir area is spatially divided into multiple grid units according to a preset size using a regular grid partitioning method. All grid units are then stored according to a quadtree coding rule to obtain the reservoir area spatial grid skeleton model. Multibeam underwater sounding is performed on the entire reservoir area at a preset cycle to obtain sounding system parameters. These parameters are then associated and stored in the corresponding grid cells of the reservoir area spatial skeleton model, resulting in multi-time series multibeam point cloud data and corresponding hydrological auxiliary data. The sounding system parameters include the beam departure angles corresponding to each sounding point. Measuring water depth And sound velocity profile data.

3. The reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 2, characterized in that: The process of spatially partitioning the processed point cloud data using an adaptive tile partitioning method specifically includes: Using the spatial grid skeleton model of the reservoir area as the initial tile division basis, the point cloud density of the processed point cloud data in each grid cell is statistically analyzed to obtain the point cloud density distribution of each grid cell. Based on the point cloud density of each grid cell, a preset high-density threshold is applied. With preset low density threshold Perform hierarchical partitioning, where : The point cloud density of the grid cells exceeds the preset high-density threshold. The grid cells are recursively subdivided into quadtrees to obtain the first level of detail tiles; Set the point cloud density of the grid cells below a preset low-density threshold. The grid cells are reverse-merged into quadtrees and merged into the parent node tile level to obtain the third detail level tile. Set the point cloud density of the grid cells to The grid cells of the interval retain the scale of the current level node to obtain the second level of detail tile; An adaptive size set of tiles is formed based on the first level of detail tiles, the second level of detail tiles, and the third level of detail tiles; For each tile in the adaptive-size tile set, the width is extended outwards on each of its four sides by [missing value]. The overlapping area is set inside the static constraint boundary of the reservoir area with a width of The boundary buffer is used to obtain a multi-layered tile structure.

4. The reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 3, characterized in that: The formula for calculating the comprehensive feature importance score W(p) is as follows: ； ； in, Indicates the temporal scouring and silting differential intensity. For point In the Elevation values ​​in a historical time series This represents the total number of historical time series preceding the current time series. This is the hydrological similarity distance attenuation coefficient. For the current time series The hydrological state vector, For the first A historical time series of hydrological state vectors The distance represents the similarity of hydrological conditions. This represents the normalized mean curvature. Represents the rate of change of the normalized normal. This represents the normalized temporal scouring and silting difference intensity. The weighting coefficients represent the normalized mean curvature. The weighting coefficients represent the normalized rate of change of the normalized normal. The weighting coefficients represent the normalized temporal sludge differential intensity.

5. A reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 2, characterized in that: Step S3 specifically includes: Using the sounding system parameters in the corresponding grid cells of the reservoir area spatial grid skeleton model, and based on the multibeam acoustic sounding error propagation relationship, each point in each tile of the multi-level point cloud tile set is... Build Diagonal anisotropic error covariance matrix And form a set of anisotropic error covariance matrices; ； in, The standard deviation of the horizontal position error in the direction perpendicular to the survey line. The standard deviation of the horizontal position error along the survey line direction. This represents the standard deviation of the elevation direction error. In each temporal point cloud tile, the tile corresponding to the first acquisition time is used as the reference tile. If the mean of the determinant of the anisotropic error covariance matrix of tiles corresponding to other time times is less than that of the tile corresponding to the first acquisition time, then the tile corresponding to that time time is replaced as the reference tile. For the target tiles corresponding to the remaining time times, the reference tile is used as the registration reference, and the rotation matrix is ​​solved based on the anisotropic error covariance matrix and the objective function. With translation vector The process continues until the change in the objective function between two consecutive iterations is lower than a preset convergence threshold, at which point the converged rotation matrix is ​​output. With translation vector : ； in, This represents the three-dimensional coordinate vector of the i-th point in the target tile. Indicates the reference tile and The corresponding three-dimensional coordinate vector of the nearest neighbor matching point, Represents the i-th point in the target tile. The anisotropic error covariance matrix, Represents the i-th corresponding point in the reference tile. The anisotropic error covariance matrix; Based on the converged rotation matrix With translation vector The target temporal point cloud tiles are transformed to obtain spatially aligned temporal tiles. The spatially aligned point cloud tiles are fused by weighted average, and the spatial difference features of each point cloud are preserved based on the temporal labels in the 3D composite index, resulting in a set of point cloud tiles with multi-temporal registration.

6. The reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 1, characterized in that: The method of using historical hydrological auxiliary data to differentiate the spatial grid skeleton model of the reservoir area based on hydrological state distance specifically includes: Obtain hydrological auxiliary data for the current time series, including the current time series flow rate. Sand content and water level And construct the current time series hydrological state vector. Based on the historical hydrological auxiliary data of each grid cell in the reservoir area spatial grid skeleton model, the current time-series hydrological state vector is calculated. Historical hydrological state vectors at various time series Hydrological distance between : ； in, , , They are respectively the historical number The flow rate, sediment concentration, and water level at each time series, , , These represent the statistical standard deviations of flow, sediment concentration, and water level over all historical time series. Will Below the preset similarity threshold Historical time series are identified as historically similar time series with similar hydrological conditions to the current time series, resulting in a set of historically similar time series. ; In the statistical reservoir area spatial grid skeleton model, each tile is in Mean value of point cloud elevation variation across historical time series and compared with the preset elevation change threshold. Compare: when Exceeding the preset elevation change threshold If the corresponding tile is identified as a tile with significant siltation changes, i.e., a highly active area, then the elevation change confirmation threshold is [value missing]. The area formed by extending a ring of adjacent tiles outward from the spatial distribution range of the high-activity area in the previous time series is defined as the boundary extension zone. The elevation change confirmation threshold of the boundary extension zone is used. ,in ; when Below the preset elevation change threshold If the elevation change is minimal, the corresponding tile is classified as a low-activity zone, and the elevation change confirmation threshold is [value missing]. ; Write the partition status of each tile and the corresponding elevation change confirmation threshold into the tile status flag of the three-dimensional composite index to obtain the partition status including high-activity area, low-activity area and boundary extension area.

7. The reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 6, characterized in that: The dynamic updating of the initial irregular triangular mesh sediment model based on the partition status specifically includes: The planar coordinates of the newly added point cloud data in the current time series are matched with the quadtree spatial index in the three-dimensional composite index to find the tile spatial code in which the newly added point cloud data falls, determine the set of tiles affected by the newly added data, and read the tile status flag in the three-dimensional composite index of each tile to obtain the partition status of each tile and the corresponding elevation change confirmation threshold, thus obtaining the set of tiles to be updated in the current time series and the elevation change confirmation threshold of each tile. For each tile in the updated tile set, the elevation change of each point in the newly added point cloud data is compared with the elevation change confirmation threshold corresponding to that tile. The subset of point cloud data whose elevation change exceeds the corresponding elevation change confirmation threshold is retained as the current time series valid change point cloud subset, and the point cloud data whose elevation change does not exceed the corresponding elevation change confirmation threshold is determined as an invalid change subset. The effective variation point cloud subsets within each tile are integrated into the existing point cloud of the corresponding tile and recursively split into left and right subsets until the number of point clouds in the subsets does not exceed the preset splitting threshold. Delaunay triangulation is performed independently on each subset to reconstruct the local irregular triangular mesh of the tile. The initial irregular triangular mesh is reused for invalid variation subsets within each tile; The reconstructed local irregular triangular mesh of each tile is merged with the reused initial irregular triangular mesh, and the current time series label is written into the time series dimension of the three-dimensional composite index of the corresponding tile to obtain the dynamic reservoir sediment model.

8. The reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 1, characterized in that: Step S5 specifically includes: Based on the vertices of the triangular network corresponding to adjacent time series in the dynamic reservoir sediment model, the sedimentation parameters at each vertex are calculated, including sedimentation thickness. siltation rate and the amount of sediment in the area where each triangular facet is located. : ； ； ； in, The accumulation thickness between adjacent time sequences at the vertices of the triangular mesh. For the current time series The corresponding elevation value of the vertex The elevation value of the corresponding vertex in the adjacent preceding time series. The accumulation rate at that vertex. Adjacent time series and The time interval between The amount of sediment in the target area, For the first The area of ​​each triangular facet. For the first The thickness of the sediment at the three vertices of the triangular facet The mean; By loading irregular triangular mesh model tiles corresponding to each time series and detail level of the dynamic reservoir sediment model using a 3D engine, and mapping each sedimentation parameter onto the surface of the dynamic reservoir sediment model in the form of color gradient, an interactive 3D dynamic sedimentation status map is obtained.

9. A reservoir sediment modeling method based on multi-time series multibeam bathymetry point cloud data as described in claim 8, characterized in that: The reservoir sediment modeling method also includes: A time-series prediction model was constructed using a long short-term memory neural network. The model was trained by taking the historical sedimentation parameter sequences and corresponding hydrological auxiliary data sequences as inputs. By inputting preset future working condition scenario parameters into the trained time series prediction model, the predicted results of sediment deposition volume, deposition range, and elevation changes within a specified time period under that working condition scenario are obtained. The preset future working condition scenario parameters are preset water level values ​​for a specified future time period set according to engineering requirements. Preset flow rate and preset values ​​of sand content .

Images