A three-dimensional terrain modeling method based on surveying data

Abstract

The present application relates to the technical field of three-dimensional data processing, and discloses a three-dimensional terrain modeling method based on surveying and mapping data, comprising the following steps: step one, obtaining multi-source heterogeneous surveying and mapping point cloud data of a region to be modeled, unifying the multi-source heterogeneous surveying and mapping point cloud data to a preset geographic coordinate system, and performing elevation normalization processing on underwater data by using synchronous tidal level data; step two, identifying homonymous feature elements of the processed multi-source heterogeneous surveying and mapping point cloud data in an overlapping area, calculating the vertical offset parameters between each load by using a least square adjustment algorithm, and performing elevation compensation on the multi-source heterogeneous surveying and mapping point cloud data by using the vertical offset parameters. The three-dimensional terrain modeling method adopts the overlapping area homonymous feature element matching and vertical offset parameter calculation technology, eliminates the systematic elevation difference of multi-source surveying and mapping data, achieves the unification of water and land benchmarks, realizes seamless smooth splicing of the terrain of the boundary area, and solves the problem of false steep cliffs and geometric misplacement at the water-land junction caused by traditional modeling.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional data processing technology, specifically a three-dimensional terrain modeling method based on surveying data. Background Technology

[0002] With the deepening development of concepts such as "Digital Earth," "Smart Cities," and "Digital Twin Watersheds," constructing high-precision, comprehensive, and seamlessly integrated 3D terrain models has become a crucial foundational task in fields such as water conservancy engineering planning, waterway dredging and maintenance, flood disaster simulation, and ecological environment monitoring. A realistic 3D terrain model needs to accurately represent the undulating morphology of the land surface and completely recreate the geometric features of the underwater riverbed. Especially in areas where land and water meet, the continuity and logical consistency of the terrain directly determine the accuracy of hydrodynamic numerical simulations and the reliability of engineering decisions.

[0003] To acquire high-precision geographic information data across the entire space, modern surveying and mapping technology has gradually shifted from a single operational mode to an integrated "land, sea, and air" collaborative operation mode. In actual operations, UAVs or airborne LiDAR equipped with high-resolution cameras are typically used to scan land and tidal flat areas to acquire high-density land point cloud data. Simultaneously, multibeam echo sounders mounted on survey vessels or unmanned surface vessels are used to conduct full-coverage underwater topographic surveys to acquire underwater sonar point cloud data. This method of acquiring multi-source heterogeneous data greatly expands the coverage of surveying and mapping operations, enabling the collection of all elements from land to underwater. However, it also brings enormous challenges to the fusion and post-processing of massive amounts of heterogeneous data. Especially in shallow water areas and intertidal zones with complex topographic environments, the differences in sensor physical characteristics, the lack of unified spatiotemporal references, and the interference of complex environmental noise make it a critical technical problem to be solved in the field of surveying and mapping geographic information.

[0004] However, existing 3D terrain modeling techniques still have the following significant defects and shortcomings when processing the aforementioned multi-source heterogeneous mapping data:

[0005] First, in the process of integrated land and water modeling, the surveying data from different sources have deviations in elevation benchmarks or differences in tide level correction models, which causes geometric misalignment between land topographic data and underwater bathymetry data at the junction. This results in false step or steep slope phenomena on the surface of the generated model, which seriously damages the continuity and authenticity of the topography in the land-water interface area.

[0006] Second, shallow water environments are affected by echo noise, water turbidity, and obstruction from terrain features, resulting in large areas of voids in the original survey data. Traditional techniques lack effective feature constraints when performing void interpolation and completion, which can easily introduce false depressions or false bulges into the model, causing severe distortion of underwater topography and making it impossible to restore the true landform undulations.

[0007] Third, existing terrain smoothing algorithms often focus on pursuing the visual continuity of curved surfaces, resulting in the removal of key polyline features at river cross-sections, shorelines, and river channel boundaries. Furthermore, the simplification of key landform structures means that the generated 3D model cannot retain the true geometric skeleton shape, reducing the accuracy of the terrain model in refined applications.

[0008] Fourth, current 3D terrain modeling results lack quantitative credibility assessments and cannot clearly indicate the reliability of data in different areas of the model. Due to the lack of evaluation criteria for data quality, operators find it difficult to determine which locations are high-risk error zones requiring supplementary surveys, resulting in a lack of evaluation standards for surveying results and increasing the risk of rework.

[0009] To address the aforementioned issues, this invention proposes a 3D terrain modeling method based on surveying data. Through vertical offset compensation, the construction of hydraulic monotonic constraints, and multi-objective joint optimization, seamless land-water integration, preservation of landform features, and quantitative evaluation of model quality are achieved. Summary of the Invention

[0010] To address the shortcomings of existing technologies, this invention provides a three-dimensional terrain modeling method based on surveying data, in order to solve the problems mentioned in the background section.

[0011] To achieve the above objectives, the present invention provides the following technical solution: a three-dimensional terrain modeling method based on surveying data, comprising:

[0012] Step 1: Obtain multi-source heterogeneous mapping point cloud data of the area to be modeled, unify the multi-source heterogeneous mapping point cloud data to the preset geographic coordinate system, and use synchronous tide data to perform elevation normalization processing on the underwater data.

[0013] Step 2: Identify the same-name feature elements in the overlapping area of ​​the processed multi-source heterogeneous mapping point cloud data, solve the vertical offset parameters between each load using the least squares adjustment algorithm, and use the vertical offset parameters to perform elevation compensation on the multi-source heterogeneous mapping point cloud data.

[0014] Step 3: Establish a vertical error model for the compensated survey points, and use the vertical error model to calculate the adaptive weights of each survey point.

[0015] Step 4: Extract structural polylines reflecting landform changes from the compensated multi-source heterogeneous mapping point cloud data, and construct a directed constraint graph with monotonically decreasing elevation logic by combining the preset water flow direction.

[0016] Step 5: Use the obtained adaptive weights as coefficients, and combine them with the extracted structural polyline and directed constraint graph to construct a joint optimization objective function that includes data fidelity terms, polyline preservation terms, and flow direction monotonic constraint terms;

[0017] Step 6: Use numerical optimization algorithms to iteratively solve the constructed joint optimization objective function, calculate the elevation values ​​of each node of the terrain grid, and generate the surface of the three-dimensional terrain model.

[0018] Step 7: Using the residual matrix from the iterative solution process combined with the vertical error model, the elevation uncertainty of the terrain grid nodes is calculated through the error propagation law to generate a spatial reliability layer.

[0019] Step 8: Calculate the execution deviation of the 3D terrain model surface under the directed constraint graph, calculate the average violation, and output the 3D terrain model and the corresponding conformity certificate.

[0020] Preferably, step one further includes:

[0021] Sub-step The system acquires multi-source heterogeneous mapping point cloud data, including UAV image point clouds, multibeam sonar point clouds, and shore-based lidar point clouds. It then uses a spatial transformation matrix to uniformly transform the UAV image point clouds, multibeam sonar point clouds, and shore-based lidar point clouds to a preset geographic coordinate system. The transformation formula for the spatial transformation matrix is ​​expressed as: ,

[0022] in, The transformed geographic coordinate vector. The original coordinate vector of the sensor. The transformation matrix contains rotation and scaling parameters. It is the coordinate translation vector;

[0023] Sub-step Synchronous tide level data is introduced to perform elevation normalization processing on the converted multibeam sonar point cloud. A time-linear interpolation algorithm is used to calculate the instantaneous tide level at the sonar sampling time, and the underwater normalized elevation value is calculated based on the instantaneous tide level value. The formula for calculating the underwater normalized elevation value is expressed as follows: ,

[0024] in, The processed underwater naturalized elevation value. These are the original elevation observations of the multibeam sonar point cloud in the sensor coordinate system. The sonar sampling time corresponds to the multibeam sonar point cloud. and This refers to adjacent preceding and subsequent observation times in synchronized tide level data. and For a moment With time The corresponding observed tide height value;

[0025] Sub-step The processed underwater normalized elevation values, as well as the processed UAV image point cloud and shore-based lidar point cloud, are spatially voxelized and divided. Density criteria are used to remove sparse outlier noise points, generating a preprocessed multi-source heterogeneous mapping point cloud dataset. The formula for the density criteria is expressed as: ≥ ,

[0026] in, voxel grid Point cloud density values ​​within, This represents the total number of point clouds that fall within the voxel grid. The volume constant of the preset voxel lattice is... This is the preset minimum point cloud density threshold.

[0027] Preferably, step two further includes:

[0028] Sub-step For the generated preprocessed multi-source heterogeneous mapping point cloud dataset, spatial overlap region retrieval is performed. The preprocessed multi-source heterogeneous mapping point cloud dataset is divided into a reference point cloud and a point cloud to be adjusted. The centroid matching method is used to extract the corresponding feature elements of the reference point cloud and the point cloud to be adjusted within the computational grid. The formula for determining the horizontal Euclidean distance of the corresponding feature elements is expressed as:

[0029] ≤ ,

[0030] in, Points in the reference point cloud Points in the differential cloud to be leveled The horizontal Euclidean distance between them and For point horizontal coordinate value, and For point horizontal coordinate value, The preset threshold for searching the same-name feature;

[0031] Sub-step Using the obtained feature elements of the same name, an observation equation based on a linear adjustment model is constructed, and the vertical offset parameter of the point cloud to be adjusted relative to the reference point cloud is solved using the least squares adjustment algorithm; the formula for solving the vertical offset parameter is expressed as: ,

[0032] in, This is the vertical offset parameter. Feature points of the same name The benchmark elevation value, Feature points of the same name The elevation value to be corrected This represents the total number of identically named feature elements used in the calculation.

[0033] Sub-step The vertical offset parameters obtained from the solution are used to perform elevation compensation on the point cloud data to be adjusted in the multi-source heterogeneous mapping point cloud data. Gross observations are identified and removed using compensation residual judgment conditions to obtain multi-source heterogeneous mapping point cloud data with completed vertical offset correction. The formula for the compensation residual judgment conditions is expressed as: ≤ ,

[0034] in, Feature points of the same name The benchmark elevation value, Feature points of the same name The elevation value to be corrected This is the vertical offset parameter. The preset vertical mean square error constant of the surveying load.

[0035] Preferably, step three further includes:

[0036] Sub-step For multi-source heterogeneous mapping point cloud data with vertical offset correction, error estimation functions are established for each mapping load. These error estimation functions are then used to calculate the elevation observation uncertainty of the mapping points. The calculation formula for the error estimation function is as follows: ,

[0037] in, For surveying points Elevation observation uncertainty, The static systematic error constant of the mapping load, This refers to the dynamic proportional error coefficient of the sensor. For surveying points The measured distance value, To measure the scanning deflection angle when the payload emits a beam or sound beam;

[0038] Sub-step Using the calculated elevation observation uncertainty, an adaptive weight for each survey point is calculated using an inverse weighting rule; the formula for calculating the adaptive weight is as follows: ,

[0039] in, For surveying points Adaptive weights, For surveying points Elevation observation uncertainty, This is a pre-defined positive perturbation term to prevent numerical overflow;

[0040] Sub-step Based on the obtained adaptive weights, local statistical evaluation criteria are used to perform point reliability grading on multi-source heterogeneous mapping point cloud data that has undergone vertical offset correction, identifying and removing low-quality observation points with adaptive weights below a preset threshold; the determination formula for point reliability grading is expressed as: ≥ ,in, For surveying points Adaptive weights, This is the preset minimum reliability weight limit parameter.

[0041] Preferably, step four further includes:

[0042] Sub-step Feature tensor analysis is performed on a multi-source heterogeneous mapping point cloud dataset after processing, point reliability classification, and removal of low-quality observation points. Curvature determination criteria are used to identify and extract structural polyline point sets reflecting abrupt terrain changes. Structural polylines are then generated using a point-to-line algorithm. The calculation formula for the curvature determination criteria is expressed as follows: ,

[0043] in, For surveying points Surface variability value, , , For surveying points The eigenvalues ​​of the covariance matrix in the neighborhood and satisfying ≤ ≤ , Extract curvature thresholds for the preset structural polylines;

[0044] Sub-step The extracted structural polylines are used as geometric constraint boundaries, and the river centerline nodes are extracted based on the preset water flow direction. A distance-weighted interpolation algorithm is used to construct the basic mesh topology of a directed constraint graph covering the region to be modeled. The calculation formula of the distance-weighted interpolation algorithm is expressed as:

[0045] ,

[0046] in, The initial elevation values ​​of the basic grid topology nodes. Mapping points within the neighborhood Elevation observation values, The corresponding inverse distance weights, This represents the number of observation points within the neighborhood.

[0047] Sub-step Based on the constructed basic grid topology, the grid edges are weighted according to the preset water flow direction to establish a directed constraint graph that satisfies the monotonically decreasing elevation property; the determination condition for the monotonically decreasing elevation property is expressed as follows: ≥ ,

[0048] in, For upstream nodes in a directed constraint graph Elevation value, For downstream nodes in a directed constraint graph Elevation value, The preset minimum riverbed gradient threshold, For nodes With nodes The horizontal projection distance between them.

[0049] Preferably, step five further includes:

[0050] Sub-step Using the obtained adaptive weights and the generated structural polylines, an elevation data fidelity term function based on discrete grid points is constructed. A bilinear interpolation operator is used to establish the mapping relationship between the grid node elevations and the observed elevations of the survey points. The calculation formula for the data fidelity term function is expressed as follows:

[0051] ,

[0052] in, The total energy value of the data fidelity function. This represents the total number of survey points. For surveying points Adaptive weights, For surveying points The set of nodes in the grid cell. For grid nodes At the location of the survey point The bilinear interpolation basis function values ​​at the location, For the grid nodes to be solved The elevation is unknown. For surveying points The actual elevation observation value;

[0053] Sub-step An anisotropic diffusion tensor is established based on the extracted structural polylines. The polyline blocking coefficient is calculated for each connecting edge in the mesh topology, and a polyline preservation term function that integrates gradient smoothing and feature preservation is constructed. The calculation formula of the polyline preservation term function is expressed as follows:

[0054] ,

[0055] in, To preserve the total energy value of the piecewise linear function, The set of all connecting edges in the mesh topology. and For the unknown elevation of the edge grid nodes, For connecting edges The broken-line blocking coefficient;

[0056] When connecting edges When it intersects with any structural polygonal line on the horizontal projection plane Values When connecting edges When it does not intersect any structural polygonal line on the horizontal projection plane Values ;

[0057] Sub-step Using the constructed directed constraint graph, an asymmetric penalty mechanism is introduced for grid node pairs that violate monotonicity logic. A flow monotonic constraint term function is constructed, and a joint optimization objective function is generated by combining the data fidelity term function and the piecewise linear preservation term function. The calculation formula of the joint optimization objective function is expressed as follows:

[0058] ,

[0059] in, To jointly optimize the overall objective value of the objective function, , , The preset weighting coefficients for data items, line graph items, and flow direction items are: Let be the set of directed edges in a directed constraint graph. For the upstream node Pointing to downstream nodes The directed edge, and The corresponding grid node elevation is unknown. The preset minimum forced elevation difference threshold, It is a linear rectification function.

[0060] Preferably, step six further includes:

[0061] Sub-step For the constructed joint optimization objective function, the first-order partial derivatives with respect to the unknown elevations of each node in the terrain grid are calculated to obtain the global gradient vector at the current iteration step; the formula for calculating the global gradient vector is expressed as: ,

[0062] in, For the first The global gradient vector at the nth iteration For the first Global grid node elevation vector at the next iteration For Hamiltonian operators, , , The preset weighting coefficients for data items, line graph items, and flow direction items are: For the gradient components of the data fidelity term, For the gradient components of the piecewise linear preservation term, The gradient components of the flow direction towards the monotonic constraint term;

[0063] Sub-step Using the calculated global gradient vector, a reverse gradient update is performed on the global grid node elevation vector using a preset step size factor to generate an updated global grid node elevation vector; the calculation formula for the reverse gradient update is expressed as: ,

[0064] in, For the first The global grid node elevation vector after the next iteration For the first Global grid node elevation vector at the next iteration The preset learning rate step size factor, For the first The global gradient vector at the next iteration;

[0065] Sub-step The convergence of the updated global grid node elevation vectors is checked, the Euclidean distance norm between two adjacent iterations is calculated, the iteration termination state is determined using the truncation error criterion, and the surface of the 3D terrain model is output; the formula for the truncation error criterion is expressed as:

[0066] ≤ ,

[0067] in, For vectors Norm notation, This represents the total number of terrain grid nodes. For grid nodes In the Elevation value at the next iteration For grid nodes In the Elevation value at the next iteration This is a preset iterative convergence threshold constant.

[0068] Preferably, step seven further includes:

[0069] Sub-step Posterior error statistics are performed on the output 3D terrain model surface. The unit weight variance factor is calculated using the sum of squared residuals from the data fidelity term to evaluate the overall fit of the grid model with respect to the original observation data. The formula for calculating the unit weight variance factor is as follows:

[0070] ,

[0071] in, Unit weighted variance factor This represents the total number of survey points. The total number of grid nodes. For surveying points Adaptive weights, For the surface of the 3D terrain model at the survey points Interpolated elevation at the location, For surveying points The actual elevation observation value;

[0072] Sub-step The constructed joint optimization objective function is used to calculate the second-order partial derivative of the objective function with respect to the grid node elevations to construct a sparse Hessian matrix. The elevation uncertainty of each grid node is then calculated using the approximate values ​​of the diagonal elements of the sparse Hessian matrix combined with the unit weight variance factor. The formula for calculating the elevation uncertainty is as follows: ,

[0073] in, For grid nodes Elevation uncertainty, Unit weighted variance factor To jointly optimize the sparse Hessian matrix corresponding to the objective function in the th... Line number The diagonal element values ​​of the column. This is a preset regularization damping coefficient to prevent the denominator from being zero;

[0074] Sub-step Using the calculated elevation uncertainty, a nonlinear mapping relationship between the uncertainty value and spatial reliability is established using a negative exponential decay function, generating a spatial reliability layer covering all grid nodes; the mapping formula for the spatial reliability is expressed as: ,

[0075] in, For grid nodes Spatial credibility score For elevation uncertainty, This is the preset standard deviation constant for the project's target accuracy.

[0076] Preferably, step eight further includes:

[0077] Sub-step Using the generated 3D terrain model surface and the constructed directed constraint graph, all directed edges in the graph are traversed, and the hydraulic logic violation value of each directed edge is calculated using a monotonicity test function; the formula for calculating the hydraulic logic violation value is expressed as:

[0078] ,

[0079] in, For directed edges Hydraulic logic violation value, As an upstream node, For downstream nodes, and For the surface of the 3D terrain model at the nodes With nodes The final elevation value at the location, The preset minimum riverbed gradient threshold, For nodes With nodes The horizontal projected distance between them;

[0080] When the calculation result is less than The time value is , indicating that the constraint is satisfied when the calculation result is greater than Retaining this value at this time indicates a violation of the constraint;

[0081] Sub-step All hydraulic logic violation values ​​obtained through statistical calculation are used to calculate the average violation amount using a global averaging algorithm, and the logic violation rate is calculated using a counting statistics method; the calculation formulas for the average violation amount and the logic violation rate are expressed as follows:

[0082] , ,

[0083] In the formula, For the average violation amount, Let be the set of directed edges of a directed constraint graph. For directed edges Hydraulic logic violation value, This represents the total number of directed edges. For logical violation rate, For hydraulic logic violation values ​​greater than The total number of directed edges;

[0084] Sub-step The average violation rate and logical violation rate are used to generate a consistency score using a weighted scoring model. A consistency certificate is then constructed, comprising the average violation rate, logical violation rate, and consistency score. The consistency certificate is then packaged and output along with the 3D terrain model surface. The formula for calculating the consistency score is as follows: ,

[0085] in, For consistency scoring, The pre-set penalty weighting coefficient for violations. The default violation rate penalty weighting coefficient is set. For the average violation amount, This represents the logical violation rate.

[0086] Step 1 involves unifying multi-source heterogeneous point cloud data into a preset geographic coordinate system and using synchronized tidal data to perform elevation normalization processing on underwater data, thereby eliminating spatial reference differences between different sensors and dynamic water level errors caused by tidal fluctuations.

[0087] Step two involves identifying the same-named feature elements in the overlapping area of ​​the processed multi-source heterogeneous mapping point cloud data and calculating the vertical offset parameters to accurately correct the elevation misalignment caused by system errors under different loads. This effectively eliminates the false step phenomenon in the water-land boundary zone and achieves seamless splicing and precise alignment of multi-source mapping data in the vertical direction.

[0088] Step 3: By establishing a vertical error model for the compensated survey points and calculating the adaptive weights of each survey point, the reliability of the survey data under different observation conditions is quantified. This ensures that the modeling process can automatically reduce the influence of low-quality data in areas with blurred edges or large incident angles, and improve the fitting weight of the terrain model to high-precision observations.

[0089] Step four involves extracting structural polylines from the compensated multi-source heterogeneous mapping point cloud data and constructing a directed constraint graph by combining it with the direction of water flow. This transforms the discrete mapping point data into a geometric network with topological relationships, providing clear morphological framework constraints and hydraulic logic guidance for terrain surface reconstruction, and preventing terrain distortions that violate physical laws from occurring in data void areas.

[0090] Step 5: By constructing a joint optimization objective function that includes data fidelity terms, polyline preservation terms, and flow direction monotonicity constraints, the survey observations, geomorphic feature lines, and hydraulic monotonicity are transformed into a multi-constraint mechanism at the mathematical level, achieving a multi-objective balance that ensures observation accuracy while maintaining topographic geometric features and physical logic rationality.

[0091] Step six involves iteratively solving the constructed joint optimization objective function using numerical optimization algorithms to approximate the optimal terrain surface solution under multiple constraints. This automatically fills in the spatial missing areas of the survey data and smooths out random measurement noise, generating a high-precision three-dimensional terrain model surface that combines geometric continuity, distinctive features, and physical realism.

[0092] Step 7: By using the residual matrix from the iterative solution process in conjunction with the vertical error model to perform error propagation calculation, the implicit mathematical statistical error is transformed into the explicit terrain grid node elevation uncertainty, which intuitively shows the accuracy distribution of each region of the terrain model and provides a quantitative scientific basis for operators to identify low-quality data areas and formulate supplementary measurement plans.

[0093] Step 8: By statistically analyzing the execution deviation of the 3D terrain model surface under the directed constraint graph and calculating the average violation, a quality acceptance standard for the 3D terrain model is established from the hydraulic logic level to ensure that the geometric accuracy of the delivered results meets the standards and conforms to the water flow physics laws of natural river channels, effectively reducing the model applicability risk in engineering applications.

[0094] An electronic device, characterized in that it comprises: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to, when executing the instructions, implement the steps of the three-dimensional terrain modeling method based on surveying data.

[0095] This invention provides a three-dimensional terrain modeling method based on surveying data. It has the following beneficial effects:

[0096] 1. This invention employs overlapping region homonymous feature element matching and vertical offset parameter calculation technology to eliminate systematic elevation differences in multi-source surveying data, achieve the unification of land and water benchmarks, realize seamless and smooth splicing of terrain in boundary areas, and solve the shortcomings of traditional modeling that results in false steep slopes and geometric misalignments at the junction of land and water.

[0097] 2. This invention employs a directed constraint graph construction based on water flow direction and a physical constraint technique of monotonically decreasing elevation. In the void interpolation, it forces the model to follow the hydraulic logic, thereby suppressing abnormal terrain undulations, restoring the true riverbed morphology, and solving the problem of false depressions or uplifts caused by missing data and noise.

[0098] 3. This invention employs geomorphic structure polyline extraction and multi-objective joint optimization solution technology, introducing the geometric skeleton as a hard constraint into the smoothing process, achieving the effect of retaining steep slopes and riverbank features while removing noise, thus realizing the complete preservation of the topographic geometric skeleton and solving the shortcomings of conventional smoothing algorithms in erasing key structural features.

[0099] 4. This invention employs uncertainty calculation and spatial reliability layer generation technology based on the error propagation law to quantitatively evaluate the numerical reliability of each node in the model, achieving a visual early warning effect for high-risk error areas, realizing refined quality acceptance of surveying and mapping results, and solving the problem of terrain models lacking quantitative evaluation standards. Attached Figure Description

[0100] Figure 1 This is a flowchart of a three-dimensional terrain modeling method based on surveying data according to the present invention. Detailed Implementation

[0101] To enable those skilled in the art to understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. Other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort should fall within the scope of protection of the present invention.

[0102] The present invention will now be described in detail with reference to the accompanying drawings:

[0103] Example 1: Implementation of Refined Topographic Modeling Based on Inland River Management Projects

[0104] Please see the appendix Figure 1 This embodiment is applied to a digital survey project for river management in a mountainous area. The survey area is approximately 2 kilometers long and includes undulating riverbank slopes, shallow water tidal flats, and deep-water main navigation channels.

[0105] Step one involves using a multi-rotor UAV equipped with a five-lens tilting camera to conduct low-altitude photogrammetry of the riverbank and tidal flat areas, generating high-density land image point cloud data. Simultaneously, an unmanned surface vessel (USV) equipped with a multibeam echo sounder performs a full-coverage scan of the underwater river channel, acquiring underwater sonar point cloud data. The data is then uniformly converted to the CGCS2000 coordinate system using a pre-defined parameter transformation model. Synchronous water level data is obtained from water level gauges deployed upstream and downstream of the survey area during the measurement period. The instantaneous water level value at the moment of USV operation is calculated using a time-linear interpolation method, thus normalizing the elevation of the underwater sonar point cloud data from the instantaneous water surface to a unified geodetic datum.

[0106] Step two: The system automatically retrieves the spatial overlap between UAV imagery point clouds and underwater sonar point clouds in the shallow water tidal flat area, dividing this overlap region into a regular computational grid with sides of 1 meter. Using the centroid matching method, corresponding feature elements within each grid are extracted. Calculations reveal a vertical systematic bias of 15 centimeters on average between the two sets of data in the overlap region. The vertical bias parameter is calculated using the least squares adjustment algorithm and applied to the underwater sonar point cloud data. A uniform elevation addition compensation is performed on all underwater measurement points to eliminate the systematic stratification caused by equipment installation errors.

[0107] Step 3: For the compensated mapping points, the system establishes a vertical error model based on the physical characteristics of the sensors. For multibeam sonar point clouds, the error value is calculated based on the beam incidence angle, and higher elevation observation uncertainty is assigned to measurement points with large incident angles at the edges. For UAV image point clouds, the uncertainty is calculated based on the matching pixel reprojection error. An adaptive weight for each measurement point is calculated using an inverse weighting rule, and sparse noise points with weight values ​​below 0.1 are marked as low-quality observation points and removed, retaining high-reliability core data for subsequent calculations.

[0108] Step four involves using a curvature tensor analysis algorithm to automatically extract structural polylines reflecting the ridgeline of the riverbank, the toe of the slope, and the edge of the deep channel from the land and underwater point clouds. Combined with the natural north-to-south flow direction of the river, the centerline nodes of the river channel are extracted to construct a directed constraint graph covering the entire region. The system automatically identifies node pairs in the grid topology that violate the hydraulic logic of "high upstream, low downstream" and establishes corresponding mathematical penalty terms to ensure that the generated riverbed model conforms to the natural laws of gravity flow.

[0109] Step 5: Construct the joint optimization objective function: The first component is the data fidelity term, using the adaptive weights obtained in Step 3 to ensure the model closely adheres to the high-precision measurement points; the second component is the polyline preservation term, utilizing the structural polylines extracted in Step 4 to block the smooth diffusion of the mesh and prevent the riverbank edges from being blurred; the third component is the flow direction monotonic constraint term, applying a quadratic penalty to all mesh edges that violate the monotonically decreasing logic. The weight coefficients for the data term are set to 1.0, the polyline term to 0.8, and the flow direction term to 1.5.

[0110] Step six involves numerically iterating the joint optimization objective function using the conjugate gradient method. After 50 iterations, the global gradient vector approaches zero, and the truncation error is less than the preset convergence threshold of 1 mm. The system outputs the final elevation values ​​of each node in the terrain mesh, generating a 3D terrain mesh model with a resolution of 0.5 meters. This model preserves the steep features of the riverbanks while automatically filling in shallow water blind spots, and the riverbed elevation exhibits a continuous and natural downward trend.

[0111] Step seven involves calculating the posterior elevation uncertainty of each grid node using the diagonal elements of the Hessian matrix and the unit weight variance factor obtained during the solution process. The results show that the elevation uncertainty is approximately 5 cm in the data overlap area and the flat riverbed area; however, it rises to 20 cm in the data hole filling area. Based on this, the system generates a red-yellow-green three-color confidence layer to visually display the model's accuracy distribution.

[0112] Step 8: The system traverses all edges of the 3D terrain model under the directed constraint graph and statistically finds that 0.5% of the edges have minor hydraulic logic violations, mainly in local deep pool areas, with an average violation amount of less than 2 cm. The overall consistency score is calculated to be 98 points. The system automatically generates a consistency certificate containing this score and violation rate statistics, and packages it with the 3D model file for delivery to the engineering unit.

[0113] Example 2: Implementation of Complex Terrain Modeling Based on Coastal Island and Reef Engineering

[0114] Please see the appendix Figure 1 This embodiment is applied to a preliminary survey project for the construction of a wharf around an island. The survey area has a complex environment, including exposed rocks, intertidal mudflats, and underwater reefs.

[0115] Step one involves using an airborne lidar to scan the island and surrounding exposed reefs to obtain land point clouds; and using a multibeam echo sounder to measure the deep water area to obtain water depth point clouds. Due to the dramatic tidal changes, high-precision tide gauge data for this sea area is introduced, and a time-linear interpolation algorithm is used to calculate the instantaneous tide height at the time of each sonar beam emission, accurately reducing the water depth data to the local theoretical depth datum.

[0116] Step two: The system identifies a spatial overlap between the intertidal zone exposed by the airborne lidar and the multibeam echo sounder at low tide. By extracting rock feature points within the overlapping area as corresponding elements, a 20-centimeter vertical deviation is calculated, primarily caused by residual errors in the sound velocity profile correction. The calculated vertical offset parameters are then used to perform overall correction of the underwater point cloud, achieving precise unification of the land-sea benchmark.

[0117] Step 3: Establish an error model based on the scanning distance of the lidar and the beam angle of the sonar. For lidar data, the farther the distance from the scanning center, the higher the uncertainty; for sonar data, the larger the beam incident angle, the higher the uncertainty. Calculate the adaptive weights for each measurement point, and remove outlier noise points generated by fish schools and bubbles in the water body through local density statistics.

[0118] Step four involves analyzing the rate of change of the normal vector on the point cloud surface to extract the cliff lines at the island's edge and the top contours of the underwater reefs as structural polygons. A pre-defined water flow direction is set, radiating outwards from the island's center. The ebb tide direction is simulated, and a directed constraint graph is constructed to force the terrain elevation to generally decrease along the radial direction, preventing the appearance of ring-shaped depressions that violate physical laws in the interpolation area.

[0119] Step 5: Construct a joint optimization objective function, focusing on strengthening the weight of the polyline preservation term to protect the sharp edge features of reefs and cliffs. Specifically, the weight coefficient for the data fidelity term is set to 1.0, the weight coefficient for the polyline preservation term is set to 2.0 to handle complex and fragmented terrain, and the weight coefficient for the flow direction monotonicity constraint term is set to 0.5 to allow for some local undulations, such as seafloor sand waves.

[0120] Step six involves iterative optimization using the steepest descent method. A small learning rate step size factor is set to avoid oscillations due to the complex terrain surrounding the island. After 120 iterations, the objective function converges. The generated model successfully reconstructs the layered structure of the underwater reef, and in shallow, fractured waters inaccessible to lidar and survey vessels, it naturally generates reasonable connection surfaces through optimization algorithms.

[0121] Step 7: Based on the optimized residual matrix, calculate the elevation uncertainty of the grid nodes. The results show that in rock and seabed areas covered by measured data, the confidence score is above 90; in shallow-water fracture zones generated by interpolation, the confidence score automatically decreases to around 60. The system generates a spatial confidence layer rendered in grayscale gradients, clearly identifying low-confidence areas that require further focused measurements.

[0122] Step eight involves performing a hydraulic logic check on the final model. The calculated average elevation violation in the low tide direction is 3 cm, with a logic violation rate of 1.2%. Based on the weighted scoring model, the consistency score of this result is 95 points. The system outputs a 3D digital elevation model file and a corresponding consistency certificate, providing basic geographic data with quantifiable quality commitments for wharf site selection design.

[0123] Embodiments of the present invention have been presented and described. It will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to the embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A three-dimensional terrain modeling method based on surveying data, characterized in that, include: Step 1: Obtain multi-source heterogeneous mapping point cloud data of the area to be modeled, unify the multi-source heterogeneous mapping point cloud data to the preset geographic coordinate system, and use synchronous tide data to perform elevation normalization processing on the underwater data. Step 2: Identify the same-name feature elements in the overlapping area of ​​the processed multi-source heterogeneous mapping point cloud data, solve the vertical offset parameters between each load using the least squares adjustment algorithm, and use the vertical offset parameters to perform elevation compensation on the multi-source heterogeneous mapping point cloud data. Step 3: Establish a vertical error model for the compensated survey points, and use the vertical error model to calculate the adaptive weights of each survey point. Step 4: Extract structural polylines reflecting landform changes from the compensated multi-source heterogeneous mapping point cloud data, and construct a directed constraint graph with monotonically decreasing elevation logic by combining the preset water flow direction. Step 5: Use the obtained adaptive weights as coefficients, and combine them with the extracted structural polyline and directed constraint graph to construct a joint optimization objective function that includes data fidelity terms, polyline preservation terms, and flow direction monotonic constraint terms; Step 6: Use numerical optimization algorithms to iteratively solve the constructed joint optimization objective function, calculate the elevation values ​​of each node of the terrain grid, and generate the surface of the three-dimensional terrain model. Step 7: Using the residual matrix from the iterative solution process combined with the vertical error model, the elevation uncertainty of the terrain grid nodes is calculated through the error propagation law to generate a spatial reliability layer. Step 8: Calculate the execution deviation of the 3D terrain model surface under the directed constraint graph, calculate the average violation, and output the 3D terrain model and the corresponding conformity certificate.

2. The three-dimensional terrain modeling method based on surveying data according to claim 1, characterized in that, Step one further includes: Sub-step The system acquires multi-source heterogeneous mapping point cloud data, including UAV image point clouds, multibeam sonar point clouds, and shore-based lidar point clouds. It then uses a spatial transformation matrix to uniformly transform the UAV image point clouds, multibeam sonar point clouds, and shore-based lidar point clouds to a preset geographic coordinate system. The transformation formula for the spatial transformation matrix is ​​expressed as: , in, The transformed geographic coordinate vector. The original coordinate vector of the sensor. The transformation matrix contains rotation and scaling parameters. It is the coordinate translation vector; Sub-step Synchronous tide level data is introduced to perform elevation normalization processing on the converted multibeam sonar point cloud. A time-linear interpolation algorithm is used to calculate the instantaneous tide level at the sonar sampling time, and the underwater normalized elevation value is calculated based on the instantaneous tide level value. The formula for calculating the underwater normalized elevation value is expressed as follows: , in, The processed underwater naturalized elevation value. These are the original elevation observations of the multibeam sonar point cloud in the sensor coordinate system. The sonar sampling time corresponds to the multibeam sonar point cloud. and This refers to adjacent preceding and subsequent observation times in synchronized tide level data. and For a moment With time The corresponding observed tide height value; Sub-step The processed underwater normalized elevation values, as well as the processed UAV image point cloud and shore-based lidar point cloud, are spatially voxelized and divided. Density criteria are used to remove sparse outlier noise points, generating a preprocessed multi-source heterogeneous mapping point cloud dataset. The formula for the density criteria is expressed as: ≥ , in, voxel grid Point cloud density values ​​within, This represents the total number of point clouds that fall within the voxel grid. The volume constant of the preset voxel lattice is... This is the preset minimum point cloud density threshold.

3. The three-dimensional terrain modeling method based on surveying data according to claim 2, characterized in that, Step two further includes: Sub-step For the generated preprocessed multi-source heterogeneous mapping point cloud dataset, spatial overlap region retrieval is performed. The preprocessed multi-source heterogeneous mapping point cloud dataset is divided into a reference point cloud and a point cloud to be adjusted. The centroid matching method is used to extract the corresponding feature elements of the reference point cloud and the point cloud to be adjusted within the computational grid. The formula for determining the horizontal Euclidean distance of the corresponding feature elements is expressed as: ≤ ， in, Points in the reference point cloud Points in the differential cloud to be leveled The horizontal Euclidean distance between them and For point horizontal coordinate value, and For point horizontal coordinate value, The preset threshold for searching the same-name feature; Sub-step Using the obtained feature elements of the same name, an observation equation based on a linear adjustment model is constructed, and the vertical offset parameter of the point cloud to be adjusted relative to the reference point cloud is solved using the least squares adjustment algorithm; the formula for solving the vertical offset parameter is expressed as: , in, This is the vertical offset parameter. Feature points of the same name The benchmark elevation value, Feature points of the same name The elevation value to be corrected This represents the total number of identically named feature elements used in the calculation. Sub-step The vertical offset parameters obtained from the solution are used to perform elevation compensation on the point cloud data to be adjusted in the multi-source heterogeneous mapping point cloud data. Gross observations are identified and removed using compensation residual judgment conditions to obtain multi-source heterogeneous mapping point cloud data with completed vertical offset correction. The formula for the compensation residual judgment conditions is expressed as: ≤ , in, Feature points of the same name The benchmark elevation value, Feature points of the same name The elevation value to be corrected This is the vertical offset parameter. The preset vertical mean square error constant of the surveying load.

4. The three-dimensional terrain modeling method based on surveying data according to claim 3, characterized in that, Step three further includes: Sub-step For multi-source heterogeneous mapping point cloud data with vertical offset correction, error estimation functions are established for each mapping load. These error estimation functions are then used to calculate the elevation observation uncertainty of the mapping points. The calculation formula for the error estimation function is as follows: , in, For surveying points Elevation observation uncertainty, The static systematic error constant of the mapping load, This refers to the dynamic proportional error coefficient of the sensor. For surveying points The measured distance value, To measure the scanning deflection angle when the payload emits a beam or sound beam; Sub-step Using the calculated elevation observation uncertainty, an adaptive weight for each survey point is calculated using an inverse weighting rule; the formula for calculating the adaptive weight is as follows: , in, For surveying points Adaptive weights, For surveying points Elevation observation uncertainty, This is a pre-defined positive perturbation term to prevent numerical overflow; Sub-step Based on the obtained adaptive weights, local statistical evaluation criteria are used to perform point reliability grading on multi-source heterogeneous mapping point cloud data that has undergone vertical offset correction, identifying and removing low-quality observation points with adaptive weights below a preset threshold; the determination formula for point reliability grading is expressed as: ≥ ,in, For surveying points Adaptive weights, This is the preset minimum reliability weight limit parameter.

5. The three-dimensional terrain modeling method based on surveying data according to claim 4, characterized in that, Step four further includes: Sub-step Feature tensor analysis is performed on a multi-source heterogeneous mapping point cloud dataset after processing, point reliability classification, and removal of low-quality observation points. Curvature determination criteria are used to identify and extract structural polyline point sets reflecting abrupt terrain changes. Structural polylines are then generated using a point-to-line algorithm. The calculation formula for the curvature determination criteria is expressed as follows: , in, For surveying points Surface variability value, , , For surveying points The eigenvalues ​​of the covariance matrix in the neighborhood and satisfying ≤ ≤ , Extract curvature thresholds for the preset structural polylines; Sub-step The extracted structural polylines are used as geometric constraint boundaries, and the river centerline nodes are extracted based on the preset water flow direction. A distance-weighted interpolation algorithm is used to construct the basic mesh topology of a directed constraint graph covering the region to be modeled. The calculation formula of the distance-weighted interpolation algorithm is expressed as: ， in, The initial elevation values ​​of the basic grid topology nodes. Mapping points within the neighborhood Elevation observation values, The corresponding inverse distance weights, This represents the number of observation points within the neighborhood. Sub-step Based on the constructed basic grid topology, the grid edges are weighted according to the preset water flow direction to establish a directed constraint graph that satisfies the monotonically decreasing elevation property; the determination condition for the monotonically decreasing elevation property is expressed as follows: ≥ , in, For upstream nodes in a directed constraint graph Elevation value, For downstream nodes in a directed constraint graph Elevation value, The preset minimum riverbed gradient threshold, For nodes With nodes The horizontal projection distance between them.

6. The three-dimensional terrain modeling method based on surveying data according to claim 5, characterized in that, Step five further includes: Sub-step Using the obtained adaptive weights and the generated structural polylines, an elevation data fidelity term function based on discrete grid points is constructed. A bilinear interpolation operator is used to establish the mapping relationship between the grid node elevations and the observed elevations of the survey points. The calculation formula for the data fidelity term function is expressed as follows: ， in, The total energy value of the data fidelity function. This represents the total number of survey points. For surveying points Adaptive weights, For surveying points The set of nodes in the grid cell. For grid nodes At the location of the survey point The bilinear interpolation basis function values ​​at the location, For the grid nodes to be solved The elevation is unknown. For surveying points The actual elevation observation value; Sub-step An anisotropic diffusion tensor is established based on the extracted structural polylines. The polyline blocking coefficient is calculated for each connecting edge in the mesh topology, and a polyline preservation term function that integrates gradient smoothing and feature preservation is constructed. The calculation formula of the polyline preservation term function is expressed as follows: ， in, To preserve the total energy value of the piecewise linear function, The set of all connecting edges in the mesh topology. and For the unknown elevation of the edge grid nodes, For connecting edges The broken-line blocking coefficient; When connecting edges When it intersects with any structural polygonal line on the horizontal projection plane Values When connecting edges When it does not intersect any structural polygonal line on the horizontal projection plane Values ; Sub-step Using the constructed directed constraint graph, an asymmetric penalty mechanism is introduced for grid node pairs that violate monotonicity logic. A flow monotonic constraint term function is constructed, and a joint optimization objective function is generated by combining the data fidelity term function and the piecewise linear preservation term function. The calculation formula of the joint optimization objective function is expressed as follows: ， in, To jointly optimize the overall objective value of the objective function, , , The preset weighting coefficients for data items, line graph items, and flow direction items are: Let be the set of directed edges in a directed constraint graph. For the upstream node Pointing to downstream nodes The directed edge, and The corresponding grid node elevation is unknown. The preset minimum forced elevation difference threshold, It is a linear rectifier function.

7. The three-dimensional terrain modeling method based on surveying data according to claim 6, characterized in that, Step six further includes: Sub-step For the constructed joint optimization objective function, the first-order partial derivatives with respect to the unknown elevations of each node in the terrain grid are calculated to obtain the global gradient vector at the current iteration step; the formula for calculating the global gradient vector is expressed as: , in, For the first The global gradient vector at the nth iteration For the first Global grid node elevation vector at the next iteration For Hamiltonian operators, , , The preset weighting coefficients for data items, line graph items, and flow direction items are: For the gradient components of the data fidelity term, For the gradient components of the piecewise linear preservation term, The gradient components of the flow towards the monotonic constraint term; Sub-step Using the calculated global gradient vector, a reverse gradient update is performed on the global grid node elevation vector using a preset step size factor to generate an updated global grid node elevation vector; the calculation formula for the reverse gradient update is expressed as: , in, For the first The global grid node elevation vector after the next iteration For the first Global grid node elevation vector at the next iteration The preset learning rate step size factor, For the first The global gradient vector at the next iteration; Sub-step The convergence of the updated global grid node elevation vectors is checked, the Euclidean distance norm between two adjacent iterations is calculated, the iteration termination state is determined using the truncation error criterion, and the surface of the 3D terrain model is output; the formula for the truncation error criterion is expressed as: ≤ ， in, For vectors Norm notation, This represents the total number of terrain grid nodes. For grid nodes In the Elevation value at the next iteration For grid nodes In the Elevation value at the next iteration This is a preset iterative convergence threshold constant.

8. A three-dimensional terrain modeling method based on surveying data according to claim 7, characterized in that, Step seven further includes: Sub-step Posterior error statistics are performed on the output 3D terrain model surface. The unit weight variance factor is calculated using the sum of squared residuals from the data fidelity term to evaluate the overall fit of the grid model with respect to the original observation data. The formula for calculating the unit weight variance factor is as follows: ， in, Unit weighted variance factor This represents the total number of survey points. The total number of grid nodes. For surveying points Adaptive weights, For the surface of the 3D terrain model at the survey points Interpolated elevation at the location, For surveying points The actual elevation observation value; Sub-step The constructed joint optimization objective function is used to calculate the second-order partial derivative of the objective function with respect to the grid node elevations to construct a sparse Hessian matrix. The elevation uncertainty of each grid node is then calculated using the approximate values ​​of the diagonal elements of the sparse Hessian matrix combined with the unit weight variance factor. The formula for calculating the elevation uncertainty is as follows: , in, For grid nodes Elevation uncertainty, Unit weighted variance factor To jointly optimize the sparse Hessian matrix corresponding to the objective function in the th... Line number The diagonal element values ​​of the column. This is a preset regularization damping coefficient to prevent the denominator from being zero; Sub-step Using the calculated elevation uncertainty, a nonlinear mapping relationship between the uncertainty value and spatial reliability is established using a negative exponential decay function, generating a spatial reliability layer covering all grid nodes; the mapping formula for the spatial reliability is expressed as: , in, For grid nodes Spatial credibility score For elevation uncertainty, This is the preset standard deviation constant for the project's target accuracy.

9. A three-dimensional terrain modeling method based on surveying data according to claim 8, characterized in that, Step eight further includes: Sub-step Using the generated 3D terrain model surface and the constructed directed constraint graph, all directed edges in the graph are traversed, and the hydraulic logic violation value of each directed edge is calculated using a monotonicity test function; the formula for calculating the hydraulic logic violation value is expressed as: ， in, For directed edges Hydraulic logic violation value, As an upstream node, For downstream nodes, and For the surface of the 3D terrain model at the nodes With nodes The final elevation value at the location, The preset minimum riverbed gradient threshold, For nodes With nodes The horizontal projected distance between them; When the calculation result is less than The time value is , indicating that the constraint is satisfied when the calculation result is greater than Retaining this value at this time indicates a violation of the constraint; Sub-step All hydraulic logic violation values ​​obtained through statistical calculation are used to calculate the average violation amount using a global averaging algorithm, and the logic violation rate is calculated using a counting statistics method; the calculation formulas for the average violation amount and the logic violation rate are expressed as follows: ， ， In the formula, For the average violation amount, Let be the set of directed edges of a directed constraint graph. For directed edges Hydraulic logic violation value, This represents the total number of directed edges. For logical violation rate, For hydraulic logic violation values ​​greater than The total number of directed edges; Sub-step The average violation rate and logical violation rate are used to generate a consistency score using a weighted scoring model. A consistency certificate is then constructed, comprising the average violation rate, logical violation rate, and consistency score. The consistency certificate is then packaged and output along with the 3D terrain model surface. The formula for calculating the consistency score is as follows: , in, For consistency scoring, The pre-set penalty weighting coefficient for violations. The default violation rate penalty weighting coefficient is set. For the average violation amount, This represents the logical violation rate.

10. An electronic device, characterized in that, include: processor; A memory for storing processor-executable instructions; wherein the processor is configured to, when executing the instructions, implement the steps of the three-dimensional terrain modeling method based on mapping data as described in any one of claims 1 to 9.

Images