A method for digital visualization of ocean topography based on 3D modeling

Abstract

The application discloses a kind of marine topography digital display methods based on three-dimensional modeling, comprising the following steps: step one: constructing multi-source marine data set;Step two: using the adaptive interpolation method based on position weight to carry out space fusion, generate three-dimensional space density field;Step three: using the improved Marching Cubes algorithm to carry out isosurface extraction, generate initial three-dimensional terrain model;Step four: using Delaunay triangulation reconstruction method to carry out local modeling, build complete three-dimensional terrain model;Step five: obtain texture terrain model;Step six: construct LOD tree structure and realize dynamic loading in combination with view frustum clipping and screen space error control strategy;Step seven: based on the dynamic loading result of LOD tree structure, through time controller drive node switching, carry out visual display.The application fuses three-dimensional modeling and LOD tree structure, realizes marine topography high-precision reconstruction and evolution display.

Description

Technical Field

[0001] This invention relates to the field of marine geographic information technology, and in particular to a method for digitally displaying marine terrain based on 3D modeling. Background Technology

[0002] With the increasing demand for marine resource development and seabed topography analysis, high-precision modeling and dynamic visualization technologies for deep-sea topography structures have received widespread attention. Existing 3D marine topography modeling methods mainly rely on a single data source (such as multibeam sonar) to construct surface models, or use traditional isosurface extraction methods to model regular voxel data. However, these methods generally suffer from the following problems in practical applications:

[0003] Multi-source ocean data suffers from significant differences in resolution, inconsistent coordinate systems, and uneven spatial distribution. Traditional interpolation methods exhibit poor reconstruction quality and large errors in sparse regions, making it difficult to effectively integrate heterogeneous multi-source data. The classic MarchingCubes algorithm suffers from topological anomalies and discontinuous normals, making it difficult to maintain the coherence and visual accuracy of complex seabed topography. Existing Delaunay reconstruction processes do not fully incorporate curvature or boundary features, resulting in the loss of details in high-curvature areas and a lack of effective local encryption modeling mechanisms. Texture mapping processes are mostly static texture processing, unable to overlay temporal information to display the terrain evolution process. The LOD hierarchical structure is simple to construct and lacks dynamic loading strategies based on error measurement and spatial organization, affecting model rendering efficiency and visual interactive experience.

[0004] Therefore, how to provide a digital visualization method for marine terrain based on 3D modeling is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0005] One objective of this invention is to propose a digital visualization method for marine topography based on 3D modeling. This invention integrates multi-source marine data interpolation modeling, an improved Marching Cubes algorithm, a Delaunay triangulation reconstruction method, texture mapping, and LOD multi-level dynamic loading technology to construct a high-precision temporal 3D topography model. It comprehensively realizes the efficient visualization expression of marine topography from static reconstruction to dynamic evolution process, and has the advantages of high modeling accuracy, strong timeliness, good system compatibility, and excellent display effect. It is suitable for marine scientific research, seabed topography evolution analysis, and underwater visualization simulation application scenarios.

[0006] A method for digitally displaying marine topography based on three-dimensional modeling according to an embodiment of the present invention includes the following steps:

[0007] Step 1: Acquire multibeam sonar detection data, satellite remote sensing images, underwater images, and hydrological profile images to construct a multi-source ocean dataset;

[0008] Step 2: Spatial fusion is performed on the multi-source ocean dataset using an adaptive interpolation method based on location weights to generate a three-dimensional spatial density field;

[0009] Step 3: Based on the three-dimensional spatial density field, the improved Marching Cubes algorithm is used to extract isosurfaces and generate an initial three-dimensional terrain model;

[0010] Step 4: Use the Delaunay triangulation reconstruction method to perform local modeling on the initial 3D terrain model to generate a local fine model, and then merge it with the initial 3D terrain model to construct a complete 3D terrain model;

[0011] Step 5: Perform texture mapping processing on the complete 3D terrain model to obtain a textured terrain model;

[0012] Step 6: Construct a LOD tree structure based on the textured terrain model, generate multi-level nodes according to spatial region division, and achieve dynamic loading by combining view frustum clipping and screen space error control strategies;

[0013] Step 7: Based on the dynamic loading results of the LOD tree structure, perform time-sharing scheduling and management of multiple texture terrain models from different periods, drive node switching through a time controller, and visualize the time series of terrain evolution process.

[0014] Optionally, step two specifically involves:

[0015] The spatial data in the multi-source ocean dataset are uniformly mapped to a regular cubic grid in the same three-dimensional Cartesian coordinate system, and the eight vertices of each regular cubic grid cell are set as the target points to be interpolated.

[0016] For each target point to be interpolated, a first search radius is set, and original observation data points are searched within the first search radius. The spatial Euclidean distance between each original observation data point and the target point to be interpolated is calculated.

[0017] When the number of original observation data points found is greater than or equal to the preset first number threshold, a first distance decay factor is set, and the weight value of each original observation data point is calculated.

[0018] ;

[0019] in, Indicates the first The weight values ​​of each original observation data point Indicates the first The spatial Euclidean distance between the original observation data points and the target point to be interpolated. Indicates the first distance attenuation factor;

[0020] When the number of original observation data points searched within the first search radius is less than the set first number threshold, the search radius is expanded to a second search radius, and a second distance attenuation factor is set. The second distance attenuation factor is less than the first distance attenuation factor. For the newly searched original observation data points in the second search radius, the weight value is calculated using the second distance attenuation factor.

[0021] Multiply the values ​​of all original observation data points by their corresponding weight values ​​and sum them up. Then divide by the sum of all weight values ​​to obtain the value of the target point to be interpolated.

[0022] The value of the target point to be interpolated is used as the density value corresponding to the regular cubic mesh cell;

[0023] Interpolation is performed on each regular cubic mesh element to generate a three-dimensional spatial density field.

[0024] Optionally, the improved Marching Cubes algorithm specifically includes:

[0025] Obtain regular cubic mesh elements in the three-dimensional spatial density field, and perform isosurface extraction processing on each regular cubic mesh element;

[0026] Set an isosurface threshold, determine the relationship between the density value of the eight vertices of each regular cube grid cell and the isosurface threshold. If the density value is greater than the isosurface threshold, encode it as 1; otherwise, encode it as 0. Combine the codes of the eight vertices into an eight-bit binary number, search a preset topology table, and determine the isosurface boundary combination inside the regular cube grid cell.

[0027] For isosurface boundary combinations in the topology table that are not unique or discontinuous, topology correction is performed. The topology correction includes selecting a patch reconstruction method consistent with adjacent voxel combinations in the boundary topology structure with multiple solutions.

[0028] When calculating the location of the boundary intersection, the normal-guided interpolation method is used. Between a pair of adjacent corner points on the edge of a regular cubic mesh cell, the location of the intersection point on the boundary is calculated based on the density value and location of the two adjacent corner points, combined with the direction of the three-dimensional spatial density field, and is used as the vertex position of the triangular facet.

[0029] During the isosurface extraction process, the density gradient change rate is calculated for each regular cube grid cell. When the density gradient change rate exceeds a preset density change threshold, the current regular cube grid cell is subdivided. After subdivision, the side length of each regular cube grid cell is half the side length of the original regular cube grid cell.

[0030] Merge the triangular faces of all regular cube mesh units to generate an initial 3D terrain model.

[0031] Optionally, step four specifically involves:

[0032] Extract the vertex coordinates, corresponding normal vectors, and curvature information of all triangular patches in the initial three-dimensional terrain model;

[0033] For each triangular facet, calculate the angle between the normal vectors of adjacent triangular facets and the current triangular facet, and mark triangular facets whose normal vector angle is greater than a set angle threshold as high curvature regions;

[0034] In each high curvature region, a set of boundary points is extracted based on the envelope boundary, and the average side length between the boundary points is calculated. If the average side length is greater than the set thinning threshold, new sampling points are introduced between the boundary points, and the sampling density is one-third of the original side length.

[0035] The boundary points and the newly introduced sampling points are used together as the input point set. The Delaunay triangulation reconstruction method in three-dimensional space is used to reconstruct the input point set to generate a local triangular mesh model that satisfies the cavity maximization principle.

[0036] During the Delaunay reconstruction process, all non-convex elements that do not meet the three-dimensional Delaunay conditions are removed, and all valid elements that pass the tetrahedral circumscribed ball test are retained to form a local fine model.

[0037] Obtain the boundary region between the local fine model and the initial 3D terrain model, pair vertices in the boundary whose distance is less than a set fusion threshold, and generate connecting edges by interpolation along the shortest path;

[0038] The local fine model is inserted into the corresponding position in the initial three-dimensional terrain model to replace the original triangular facets, generating a fused complete three-dimensional terrain model.

[0039] Optionally, step five specifically includes:

[0040] Acquire satellite remote sensing images, underwater images, and hydrological profile images covering the area of ​​the complete three-dimensional terrain model, and unify them to the same spatial coordinate system as the complete three-dimensional terrain model;

[0041] Extract the vertex coordinates of each triangular facet in the complete 3D terrain model, and project each vertex coordinate onto the corresponding 2D plane of the image to establish the mapping relationship between the vertices of the complete 3D terrain model and the image pixels, thereby generating texture coordinates;

[0042] The satellite remote sensing images and the underwater images are processed by texture mapping. The images are selected as texture sources according to the preset priority, and the texture mapping of all triangular patches is completed.

[0043] After the texture mapping is completed, the contour images generated based on hydrological profile data are used as additional texture layers and overlaid on the surface of the already mapped complete 3D terrain model to form a textured terrain model.

[0044] Optionally, the construction of the LOD tree structure, which generates multi-level nodes according to spatial region division, specifically involves:

[0045] Based on the spatial range of the textured terrain model, the entire three-dimensional terrain area is spatially divided into layers, and a multi-level LOD tree structure is constructed in a top-down manner. Each level is refined according to the spatial resolution, with the top level being the lowest resolution model and the bottom level being the highest resolution model.

[0046] During the construction process, the spatial bounding box of the textured terrain model is used as the root node region. The region is divided according to the preset spatial size, and sub-regions in the form of octrees are generated. The resolution model data block corresponding to each sub-region is mounted as a node in the LOD tree structure.

[0047] Optionally, the node information of the LOD tree structure includes:

[0048] Resolution level information is used to identify the resolution accuracy level of the terrain data corresponding to the node;

[0049] Error factor, used to represent the projection error in screen space between the current node resolution model and the previous level resolution model;

[0050] Spatial index information is used to uniquely identify the spatial location of a node in the LOD tree structure, including the level number and spatial coordinates in the spatial partitioning system;

[0051] The child node pointer is used to record the list of the next level resolution model nodes of the node, and supports loading and releasing on demand;

[0052] After the LOD tree structure is constructed, all nodes are connected according to their hierarchical relationship to form a complete multi-resolution data structure.

[0053] Optionally, the dynamic loading achieved by combining view frustum clipping and screen space error control strategies specifically includes:

[0054] A view frustum structure is constructed, which consists of six clipping planes, namely the left clipping plane, right clipping plane, top clipping plane, bottom clipping plane, near clipping plane, and far clipping plane, which define the three-dimensional spatial region visible to the current camera.

[0055] Traverse the bounding box information of each node in the LOD tree structure, and sequentially determine the spatial relationship with the six clipping planes of the view frustum structure. If the bounding box does not intersect with any clipping plane at all, the node is determined to be invisible and the loading is skipped; if the bounding box intersects with all clipping planes or partially intersects with all clipping planes, it is determined to be a visible node.

[0056] For each visible node, based on the resolution level boundary to which the visible node belongs, the diagonal length of the bounding box, and the center distance between the camera and the node, the projection size of the visible node on the screen is determined, and the error estimate of the visible node in the screen space is obtained.

[0057] The error estimate is compared with a preset screen error tolerance. If the error estimate is greater than the preset screen error tolerance, it is determined that the current accuracy is insufficient, and the child nodes of the visible node are loaded to further refine the judgment. If the error estimate is less than or equal to the preset screen error tolerance, or the visible node has no child nodes, the resolution model data corresponding to the visible node is loaded into the rendering queue.

[0058] Optionally, the time controller specifically includes a timeline module, a keyframe index module, and a node scheduling management module;

[0059] The timeline module is used to construct a unified timeline, set the start and end times, time steps and playback frame rates, support time granularity control at the year, month, day and hour levels, and identify the terrain status at different times with a unified logical time sequence.

[0060] The keyframe index module is used to bind multiple texture terrain models from different periods to the timeline. It uses a key-value mapping method to record the relationship between each time point and the corresponding LOD tree node index. Each keyframe corresponds to a complete LOD tree node state set, including node level, texture path and resolution level information.

[0061] The node scheduling and management module extracts the LOD tree node status of the target time point from the key frame index module according to the current time axis pointer position, performs an unloading operation on the node that needs to be replaced, and performs a loading operation on the newly added node. During the node switching process, a caching mechanism is used to save the data blocks of the nearest key frames.

[0062] The time controller works in conjunction with the LOD tree structure to achieve orderly management and continuous display of multi-temporal texture terrain models, thus completing the time-series visualization of the seabed terrain evolution process.

[0063] The beneficial effects of this invention are:

[0064] This invention constructs a model based on improved Marching The Cubes algorithm and multi-source data fusion 3D modeling workflow addresses the problems of uneven data distribution, topological discontinuities, and weak dynamic visualization capabilities in existing marine terrain modeling. It employs an adaptive interpolation strategy based on spatial distance and attenuation factors to unify marine observation data from different sources into a regular 3D mesh structure, achieving balanced reconstruction of sparse and dense areas. During voxel isosurface extraction, topological consistency discrimination, normal-guided interpolation, and local density gradient subdivision strategies are introduced to effectively avoid topological artifacts and patch discontinuities, improving the accuracy of complex terrain reconstruction. For regions with abrupt curvature changes, high curvature regions are identified using normal vector angles, and local densification reconstruction is performed using a 3D Delaunay triangulation network, preserving terrain edge details and improving the overall model's local accuracy. After model construction, a hierarchical LOD tree structure is built, embedding resolution levels, error factors, and spatial index information. Combined with view frustum clipping and screen space error estimation strategies, visible model nodes are dynamically loaded based on camera position, improving rendering efficiency and interactive response speed. Finally, a time controller enables frame-by-frame scheduling and node state switching of the multi-temporal texture terrain model, continuously expressing the terrain evolution process. This invention achieves full-process technical integration of marine terrain modeling, from multi-source data fusion, detail-preserving reconstruction, dynamic multi-level loading to temporal display, significantly improving the accuracy, coherence, and spatiotemporal visualization of terrain modeling. Attached Figure Description

[0065] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings:

[0066] Figure 1 This is an overall flowchart of a method for digitally displaying marine terrain based on 3D modeling proposed in this invention;

[0067] Figure 2 This is a flowchart illustrating the process of extracting isosurfaces and generating an initial terrain model using an improved Marching Cubes algorithm, which is part of a digital visualization method for marine terrain based on 3D modeling proposed in this invention.

[0068] Figure 3 This is a flowchart illustrating the generation of a local fine model using the Delaunay triangulation reconstruction method, which is part of a digital visualization method for marine terrain based on 3D modeling proposed in this invention. Detailed Implementation

[0069] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0070] refer to Figures 1-3 A method for digitally displaying marine topography based on 3D modeling includes the following steps:

[0071] Step 1: Acquire multibeam sonar detection data, satellite remote sensing images, underwater images, and hydrological profile images to construct a multi-source ocean dataset;

[0072] Step 2: Spatial fusion is performed on the multi-source ocean dataset using an adaptive interpolation method based on location weights to generate a three-dimensional spatial density field;

[0073] Step 3: Based on the three-dimensional spatial density field, the improved Marching Cubes algorithm is used to extract isosurfaces and generate an initial three-dimensional terrain model;

[0074] Step 4: Use the Delaunay triangulation reconstruction method to perform local modeling on the initial 3D terrain model to generate a local fine model, and then merge it with the initial 3D terrain model to construct a complete 3D terrain model;

[0075] Step 5: Perform texture mapping processing on the complete 3D terrain model to obtain a textured terrain model;

[0076] Step 6: Construct a LOD tree structure based on the textured terrain model, generate multi-level nodes according to spatial region division, and achieve dynamic loading by combining view frustum clipping and screen space error control strategies;

[0077] Step 7: Based on the dynamic loading results of the LOD tree structure, perform time-sharing scheduling and management of multiple texture terrain models from different periods, drive node switching through a time controller, and visualize the time series of terrain evolution process.

[0078] In this embodiment, step two specifically includes:

[0079] The spatial data in the multi-source ocean dataset are uniformly mapped to a regular cubic grid in the same three-dimensional Cartesian coordinate system, and the eight vertices of each regular cubic grid cell are set as the target points to be interpolated.

[0080] For each target point to be interpolated, a first search radius is set, and original observation data points are searched within the first search radius. The spatial Euclidean distance between each original observation data point and the target point to be interpolated is calculated.

[0081] When the number of original observation data points found is greater than or equal to the preset first number threshold, a first distance decay factor is set, and the weight value of each original observation data point is calculated.

[0082] ;

[0083] in, Indicates the first The weight values ​​of each original observation data point Indicates the first The spatial Euclidean distance between the original observation data points and the target point to be interpolated. This represents the first distance decay factor, which is a positive real number and is used to control the decay rate of the spatial Euclidean distance on the weight value.

[0084] When the number of original observation data points searched within the first search radius is less than the set first number threshold, the search radius is expanded to a second search radius, and a second distance attenuation factor is set. The second distance attenuation factor is less than the first distance attenuation factor. For the newly searched original observation data points in the second search radius, the weight value is calculated using the second distance attenuation factor.

[0085] Multiply the values ​​of all original observation data points by their corresponding weight values ​​and sum them up. Then divide by the sum of all weight values ​​to obtain the value of the target point to be interpolated.

[0086] The value of the target point to be interpolated is used as the density value corresponding to the regular cubic mesh cell;

[0087] Interpolation is performed on each regular cubic mesh element to generate a three-dimensional spatial density field.

[0088] In this embodiment, the improved Marching Cubes algorithm is specifically as follows:

[0089] Obtain regular cubic mesh elements in the three-dimensional spatial density field, and perform isosurface extraction processing on each regular cubic mesh element;

[0090] Set an isosurface threshold, determine the relationship between the density value of the eight vertices of each regular cube grid cell and the isosurface threshold. If the density value is greater than the isosurface threshold, encode it as 1; otherwise, encode it as 0. Combine the codes of the eight vertices into an eight-bit binary number, search a preset topology table, and determine the isosurface boundary combination inside the regular cube grid cell.

[0091] For isosurface boundary combinations in the topology table that are not unique or discontinuous, topology correction is performed. The topology correction includes selecting a patch reconstruction method consistent with adjacent voxel combinations in the boundary topology structure with multiple solutions.

[0092] When calculating the location of the boundary intersection, the normal-guided interpolation method is used. Between a pair of adjacent corner points on the edge of a regular cubic mesh cell, the location of the intersection point on the boundary is calculated based on the density value and location of the two adjacent corner points, combined with the direction of the three-dimensional spatial density field, and is used as the vertex position of the triangular facet.

[0093] In specific operation, let the positions of the adjacent corner points be respectively and Its corresponding density value is and The isosurface threshold is .like and If the points are located on either side of the isosurface threshold, it can be determined that the isosurfaces intersect on the line connecting the two points. A linear interpolation method is then used to calculate the position coefficient of the intersection point on that edge. :

[0094] ;

[0095] The initial intersection point is :

[0096] ;

[0097] To further improve the geometric accuracy of the intersection points, three-dimensional density field gradient direction information is introduced. First, the local density gradient direction vector is calculated based on the voxels surrounding the cube, representing the direction of the fastest density change. The direction vector of the interpolated line segment is then weighted and fused with the density gradient direction vector, adjusting the directional offset of the interpolation result to better match the actual density change trend, thus obtaining the corrected intersection point coordinates. The final intersection point serves as the vertex position of the isosurface triangular patch, providing a more natural and continuous surface representation for the model.

[0098] During the isosurface extraction process, the density gradient change rate is calculated for each regular cube grid cell. When the density gradient change rate exceeds a preset density change threshold, the current regular cube grid cell is subdivided. After subdivision, the side length of each regular cube grid cell is half the side length of the original regular cube grid cell.

[0099] Merge the triangular faces of all regular cube mesh units to generate an initial 3D terrain model.

[0100] In this embodiment, the topology table is a core component of the improved Marching Cubes algorithm. Its function is to establish the correspondence between the eight vertex states of a regular cube mesh cell and the combinations of isosurface boundaries. The density values ​​of the eight vertices of each cube mesh cell are compared with a set isosurface threshold. Vertices above the threshold are encoded as 1, and those below the threshold are encoded as 0, forming an eight-bit binary number, resulting in 256 vertex states. The topology table predefines the boundary combinations crossed by the isosurface in these 256 states and their corresponding triangular facet connection methods, achieving efficient isosurface extraction based on table lookup.

[0101] In certain vertex states, there may be issues with the non-uniqueness or discontinuity of isosurface boundary combinations. Specifically, the same vertex state may correspond to multiple possible triangular facet stitching methods, or the proximity of isosurfaces to voxel boundaries may lead to discontinuities in connections between different voxels, resulting in cracks or geometric inconsistencies. These problems can easily cause structural defects in the extracted model, such as topological breaks and facet misalignment, affecting the continuity and realism of the final terrain model.

[0102] To address the aforementioned issues, topological correction strategies play a crucial role in the construction of isosurfaces. When multiple solutions exist at the isosurface boundary, the consistency between the boundary connection method and the patch combination relationship of surrounding voxels is compared. The optimal patch reconstruction method, consistent with adjacent voxel combinations, is selected, ensuring the isosurface extends continuously across multiple voxels and avoiding structural abrupt changes. Furthermore, for boundary regions with potential cracks, local normal vector direction and density gradient analysis optimize the patch interpolation positions and connection order, ensuring geometric continuity and topological closure among all isosurfaces.

[0103] The improved Marching Cubes algorithm, building upon traditional algorithms, introduces a topological consistency verification mechanism, a normal-guided interpolation strategy, and a local adaptive subdivision scheme, enhancing the geometric accuracy and topological continuity of isosurface extraction. In the topology table design, by supplementing and adjusting boundary combination data, it resolves the ambiguous patch splicing problem present in traditional algorithms, enhancing control over multiple boundary solutions. During patch construction, a consistency comparison strategy with adjacent voxel patch structures is introduced, enabling seamless connections of isosurfaces across voxel regions, effectively avoiding gaps and jumps. During boundary intersection calculation, interpolation guidance using vertex normals and density gradient directions optimizes the positioning of intersection points on voxel edges, making the extracted patches more closely resemble real terrain undulations. In regions with abrupt changes in spatial density gradients, the sampling resolution is dynamically adjusted through density change rate detection, automatically refining mesh cells in highly variable areas and enhancing the capture of detailed structures. The overall algorithm retains the lookup efficiency of the Marching Cubes framework while overcoming its shortcomings of topological non-closure and large interpolation errors, making it suitable for high-quality 3D model reconstruction in complex terrain environments.

[0104] In this embodiment, step four specifically includes:

[0105] Extract the vertex coordinates, corresponding normal vectors, and curvature information of all triangular patches in the initial three-dimensional terrain model;

[0106] For each triangular facet, calculate the angle between the normal vectors of adjacent triangular facets and the current triangular facet, and mark triangular facets whose normal vector angle is greater than a set angle threshold as high curvature regions;

[0107] In each high curvature region, a set of boundary points is extracted based on the envelope boundary, and the average side length between the boundary points is calculated. If the average side length is greater than the set thinning threshold, new sampling points are introduced between the boundary points, and the sampling density is one-third of the original side length.

[0108] The boundary points and the newly introduced sampling points are used together as the input point set. The Delaunay triangulation reconstruction method in three-dimensional space is used to reconstruct the input point set to generate a local triangular mesh model that satisfies the cavity maximization principle.

[0109] During the Delaunay reconstruction process, all non-convex elements that do not meet the three-dimensional Delaunay conditions are removed, and all valid elements that pass the tetrahedral circumscribed ball test are retained to form a local fine model.

[0110] The cavity maximization principle prioritizes selecting point sets that maximize the cavity structure within the reconstructed region when building a local triangular mesh model. This principle helps generate mesh structures with reasonable topological relationships and uniform element distribution, avoiding overly flat or degenerate triangular faces, and improving the geometric stability and numerical computation efficiency of the model. Based on this principle, when inserting a new point, all existing elements within the region affected by the current point are searched, and a new non-overlapping element combination that maximizes the connection between the new point and surrounding points is reconstructed.

[0111] The 3D Delaunay condition requires that any constructed tetrahedral element must have a circumscribed sphere (i.e., a sphere passing through all four vertices of the tetrahedron) that does not contain any points from the input point set other than the vertices of the tetrahedron. Tetrahedrons satisfying the 3D Delaunay condition are considered valid elements, which can avoid generating non-convex structures, excessively thin or highly skewed meshes, and improve the overall mesh quality and the stability of numerical simulations.

[0112] The tetrahedral circumsphere check is a procedure for verifying the 3D Delaunay conditions. For each candidate tetrahedral element, a unique circumsphere is first calculated based on the coordinates of its four vertices. Then, it is checked whether the sphere contains any other sampling points besides the element's vertices. If no other points fall inside the sphere, the tetrahedral element passes the check and can be retained in the final mesh; if other points exist, the element does not satisfy the 3D Delaunay conditions and must be discarded. Through this process, the final locally refined model possesses excellent geometric structure and data consistency.

[0113] Obtain the boundary region between the local fine model and the initial 3D terrain model, pair vertices in the boundary whose distance is less than a set fusion threshold, and generate connecting edges by interpolation along the shortest path;

[0114] The local fine model is inserted into the corresponding position in the initial three-dimensional terrain model to replace the original triangular facets, generating a fused complete three-dimensional terrain model.

[0115] In this embodiment, step five specifically includes:

[0116] Acquire satellite remote sensing images, underwater images, and hydrological profile images covering the area of ​​the complete three-dimensional terrain model, and unify them to the same spatial coordinate system as the complete three-dimensional terrain model;

[0117] Extract the vertex coordinates of each triangular facet in the complete 3D terrain model, and project each vertex coordinate onto the corresponding 2D plane of the image to establish the mapping relationship between the vertices of the complete 3D terrain model and the image pixels, thereby generating texture coordinates;

[0118] The satellite remote sensing images and the underwater images are processed by texture mapping. The images are selected as texture sources according to the preset priority, and the texture mapping of all triangular patches is completed.

[0119] After the texture mapping is completed, the contour images generated based on hydrological profile data are used as additional texture layers and overlaid on the surface of the already mapped complete 3D terrain model to form a textured terrain model.

[0120] In this embodiment, the construction of the LOD tree structure, which generates multi-level nodes according to spatial region division, specifically involves:

[0121] Based on the spatial range of the textured terrain model, the entire three-dimensional terrain area is spatially divided into layers, and a multi-level LOD tree structure is constructed in a top-down manner. Each level is refined according to the spatial resolution, with the top level being the lowest resolution model and the bottom level being the highest resolution model.

[0122] During the construction process, the spatial bounding box of the textured terrain model is used as the root node region. It is divided according to the preset spatial size to generate sub-regions in the form of an octree. The resolution model data block corresponding to each sub-region is mounted as a node in the LOD tree structure.

[0123] The bounding box is a minimum axis-aligned cubic envelope representing the area occupied by the textured terrain model in 3D space, typically defined as a minimum bounding cube. Its boundary is determined by the minimum and maximum coordinates of all vertices of the textured terrain model, forming a closed cubic structure with six faces, containing all the geometric information of the entire model. When constructing the LOD tree structure, this bounding box is used as the root node region, and the space is hierarchically divided accordingly. Based on a preset spatial size threshold, the bounding box is divided into eight equal-volume sub-regions, each serving as a child node of the root node, forming the first level of the octree. Subsequently, each sub-region is recursively divided according to the same rules until the preset minimum spatial unit scale or the resolution requirement of the terrain model is reached. Each branch node records the model data index and level number within the corresponding spatial region, enabling rapid spatial location and layered loading, thereby constructing an LOD octree structure that supports multi-resolution visualization.

[0124] In this embodiment, the node information of the LOD tree structure includes:

[0125] Resolution level information is used to identify the resolution accuracy level of the terrain data corresponding to the node;

[0126] Error factor, used to represent the projection error in screen space between the current node resolution model and the previous level resolution model;

[0127] Spatial index information is used to uniquely identify the spatial location of a node in the LOD tree structure, including the level number and spatial coordinates in the spatial partitioning system;

[0128] The child node pointer is used to record the list of the next level resolution model nodes of the node, and supports loading and releasing on demand;

[0129] After the LOD tree structure is constructed, all nodes are connected according to their hierarchical relationship to form a complete multi-resolution data structure.

[0130] In this embodiment, the dynamic loading achieved by combining view frustum clipping and screen space error control strategies specifically includes:

[0131] A view frustum structure is constructed that includes parameters such as camera position, orientation, field of view, near plane, and far plane. The view frustum structure consists of six clipping planes, namely the left clipping plane, right clipping plane, top clipping plane, bottom clipping plane, near clipping plane, and far clipping plane. The view frustum structure defines the three-dimensional spatial region visible to the current camera.

[0132] Traverse the bounding box information of each node in the LOD tree structure, and determine the spatial relationship with the six clipping planes of the view frustum structure in turn. If the bounding box does not intersect with any clipping plane at all, the node is determined to be invisible and the loading is skipped; if the bounding box intersects with all clipping planes or partially intersects with all clipping planes, it is determined to be a visible node.

[0133] For each visible node, based on the resolution level boundary to which the visible node belongs, the diagonal length of the bounding box, and the center distance between the camera and the node, the projection size of the visible node on the screen is determined, and the error estimate of the visible node in the screen space is obtained.

[0134] A bounding box is a simplified geometric structure used to quickly determine the extent and spatial relationships of objects in 3D space. It adopts an axis-aligned bounding box form, that is, a rectangular cube constructed with minimum and maximum coordinate values, used to enclose all the terrain model data contained in a LOD node. Each LOD node corresponds to a bounding box in the tree structure, defining the spatial extent covered by the node.

[0135] During the view frustum clipping process, the visible nodes related to the field of view are quickly filtered out by determining the positional relationship between the bounding box and the six clipping planes of the camera's view frustum. If the bounding box is completely outside any clipping plane, it means that the node is invisible and can be skipped during loading; if the bounding box intersects with the view frustum or is completely inside it, it means that the node is within the field of view and can participate in subsequent loading judgments.

[0136] For visible nodes, their error estimate in screen space needs further evaluation to determine whether to load higher-precision child nodes. This error estimate is based on the node's resolution level boundary, the diagonal length of the bounding box, and the distance between the camera and the node's geometric center. Specifically, the geometric deviation is determined by the spatial resolution corresponding to the node's current level, resulting in the node resolution error. Based on the diagonal length of the bounding box, the actual size of the node in 3D space is estimated. Combining the actual size with the distance from the camera to the node's center, the projection size of the node on the screen is calculated using perspective projection relationships. The error estimate is the ratio between the number of pixels occupied by the projection size and the node resolution error, used to determine whether the current precision meets the display requirements.

[0137] The error estimate is compared with the preset screen error tolerance. If the error estimate is greater than the preset screen error tolerance, it is determined that the current accuracy is insufficient, and the child nodes of the visible node need to be loaded for further refinement to optimize the rendering quality and visual effect. If the error estimate is less than or equal to the preset screen error tolerance, or the visible node has no child nodes, the resolution model data corresponding to the visible node is loaded into the rendering queue.

[0138] View frustum clipping and screen space error control work together in the LOD tree node selection process to ensure that only terrain model data blocks within the current camera's field of view that meet visual accuracy requirements are loaded, thus achieving spatial and performance optimization.

[0139] In this embodiment, the time controller specifically includes a time axis module, a key frame index module, and a node scheduling management module;

[0140] The timeline module is used to construct a unified timeline, set the start and end times, time steps and playback frame rates, support time granularity control at the year, month, day and hour levels, and identify the terrain status at different times with a unified logical time sequence.

[0141] The keyframe index module is used to bind multiple texture terrain models from different periods to the timeline. It uses a key-value mapping method to record the relationship between each time point and the corresponding LOD tree node index. Each keyframe corresponds to a complete LOD tree node state set, including node level, texture path and resolution level information.

[0142] The node scheduling and management module extracts the LOD tree node status of the target time point from the keyframe index module based on the current timeline pointer position, performs an unloading operation on the nodes that need to be replaced, and performs a loading operation on the newly added nodes; during the node switching process, a caching mechanism is used to save the data blocks of the nearest keyframes;

[0143] The time controller works in conjunction with the LOD tree structure to achieve orderly management and continuous display of multi-temporal texture terrain models, thus completing the time-series visualization of the seabed terrain evolution process.

[0144] Example 1:

[0145] To verify the feasibility of this invention in practice, it was applied to the construction of a digital display system for a certain sea area. The goal is to leverage the method of this invention to perform high-precision terrain modeling, dynamic rendering, and multi-temporal evolution display of the target area, given that traditional two-dimensional topographic maps cannot accurately reflect complex seabed geomorphological structures. This comprehensively improves the visualization quality and interactive performance of marine topographic data. This scenario involves the support needs of marine geological research institutions for various aspects of work related to seabed tectonics, sedimentary evolution, and disaster assessment. In particular, the modeling and dynamic analysis of key structures such as seabed faults and sedimentary depressions place high demands on data accuracy and loading efficiency.

[0146] In this embodiment, the research team acquired raw depth point cloud data of the seabed using multibeam sonar equipment, with a data density of 12-15 points per square meter. This data was combined with satellite remote sensing imagery, underwater video images, and historical hydrological profile images of the corresponding sea area, resulting in a multi-source dataset of approximately 28 GB. Within the technical framework of this invention, all data was first unified into a regular cubic mesh structure in a three-dimensional Cartesian coordinate system, with a mesh resolution set to 5 meters.

[0147] The aforementioned multi-source data were spatially fused using an adaptive interpolation algorithm based on location weights. Using the eight vertices of each regular cube unit as target interpolation points, the number of original observation data points within search radii of 10 meters and 20 meters was counted, and weights were calculated based on actual distance and attenuation factors. The fused data yielded a smooth and internally coherent 3D density field. Based on this, an improved Marching Cubes algorithm was used for isosurface extraction. Normal-guided interpolation improved the geometric accuracy of isosurface boundaries, and automatic correction was performed for topological discontinuities. Finally, an initial 3D terrain model was generated, containing 2.6 million triangular faces with an average side length of approximately 2.3 meters and face normal deviation controlled within ±2.1°.

[0148] To further improve the local accuracy of the model, Delaunay triangulation was used for local reconstruction of high-curvature areas (such as seabed cliffs and trenches). During the reconstruction process, the included angle between adjacent angles was calculated for each high-curvature patch. If the angle exceeded 45 degrees, local refinement was triggered. The spacing between the triangular patch boundaries was controlled within 1.2 meters. After refinement, the number of patches in the region increased by approximately 4.8 times, significantly improving the model's detail representation. After fusion processing, a complete 3D terrain model was formed, with a total of 8.8 million patches.

[0149] Texture mapping was performed on the complete model, prioritizing remote sensing images with a resolution of 0.5 meters as the main texture source, supplemented by latent images for areas that were difficult to obtain. The final result was a 3D terrain model with realistic textures. To support browser loading efficiency, an eight-layer LOD tree structure was constructed based on the model, dividing spatial regions from top to bottom. The model contained a total of 17,952 nodes, with the lowest level resolution controlled within 1 meter and the highest level around 16 meters. During dynamic loading, view frustum clipping and screen space error control strategies were used to dynamically load and unload terrain nodes within the field of view, ensuring a stable average frame rate between 32 and 46 FPS on different devices.

[0150] To verify the system performance and the beneficial effects of the present invention, the performance indicators of the traditional single-layer high-precision model loading method and the LOD multi-resolution loading method of the present invention were compared, as shown in Table 1 below.

[0151] Table 1. Comparison of the method of the present invention and the traditional method under different operating scenarios.

[0152]

[0153] The comparative data in Table 1 clearly shows that the present invention has significant performance advantages over traditional single-layer modeling methods in the process of 3D terrain modeling and display in different operating scenarios.

[0154] In submarine fault regions, the traditional single-layer model method requires 875MB of data to be loaded in a single operation, with an initial loading time of 12.8 seconds, an average frame rate of only 18 FPS, and a maximum error of 5.6 pixels. In contrast, the method of this invention only requires loading 218MB of data, reducing the initial loading time to 4.4 seconds, increasing the average frame rate to 42 FPS, and reducing the maximum error to 2.1 pixels, achieving a 4.3-fold improvement in local accuracy. This demonstrates superior lightweight and high-precision capabilities.

[0155] In applications involving deposited depressions, traditional methods require 722MB of data, take 10.7 seconds to load, operate at 20 FPS, and have an error of 4.9 pixels. In contrast, this invention requires only 203MB of data and 3.8 seconds to load, increases the frame rate to 44 FPS, reduces the error to 1.8 pixels, and improves local accuracy by 4.7 times, greatly enhancing display smoothness and detail reproduction capabilities.

[0156] In flat, deep-sea areas, traditional methods require loading 641MB of data, achieving an average frame rate of 26 FPS and a maximum error of 3.2 pixels. This invention, however, only requires loading 156MB of data to achieve a high frame rate of 46 FPS, reducing the error to 1.2 pixels and improving local accuracy by 2.9 times, demonstrating its advantage of maintaining high performance even in low-complexity regions.

[0157] This embodiment demonstrates the significant advantages of the present invention in the visualization of complex terrain. It effectively solves the problems of slow loading, low frame rate, large error and insufficient local precision in traditional methods. It not only improves the expression accuracy and local detail representation of 3D terrain models, but also significantly enhances the visualization efficiency of large-scale terrain data. Furthermore, it adaptively adjusts the loaded content according to the user's perspective, achieving a unity of efficient rendering and precision control, and provides stable and reliable technical support for the visualization of multi-temporal seafloor terrain evolution processes.

[0158] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for digital presentation of marine topography based on three-dimensional modeling, characterized by, The method comprises the following steps: Step 1: Obtain multi-beam sonar detection data, satellite remote sensing images, underwater images and hydrographic profile images, and construct a multi-source marine data set; Step 2: Perform spatial fusion on the multi-source marine data set using an adaptive interpolation method based on position weight to generate a three-dimensional spatial density field; Step 3: Based on the three-dimensional spatial density field, an improved Marching Cubes algorithm is used to extract isosurface to generate an initial three-dimensional terrain model; Step 4: The initial three-dimensional terrain model is reconstructed using a Delaunay triangular network method for local modeling to generate a local fine model, and the initial three-dimensional terrain model is fused to construct a complete three-dimensional terrain model; Step 5: Perform texture mapping processing on the complete three-dimensional terrain model to obtain a texture terrain model; Step 6: Based on the texture terrain model, construct a LOD tree structure, generate multiple hierarchical nodes according to spatial region division, and realize dynamic loading combined with view frustum culling and screen space error control strategy; Step 7: Based on the dynamic loading result of the LOD tree structure, perform time scheduling and management on multiple texture terrain models of different periods, drive node switching through a time controller, and perform visual time sequence display of the terrain evolution process; The improved Marching Cubes algorithm specifically comprises: Obtain regular cubic grid cells in the three-dimensional spatial density field, and perform isosurface extraction processing on each regular cubic grid cell; Set an isosurface threshold value, judge the size relationship between the density values of the eight vertices of each regular cubic grid cell and the isosurface threshold value, if the density value is greater than the isosurface threshold value, encode it as 1, otherwise encode it as 0, and combine the eight vertex encodings into an eight-bit binary number, look up a preset topological table, and determine the isosurface boundary combination inside the regular cubic grid cell; For the isosurface boundary combinations with non-uniqueness or discontinuity in the topological table, topological correction is performed, and the topological correction includes selecting a face sheet reconstruction method consistent with the adjacent voxel combination in the boundary topological structure with multiple solutions; In the calculation of the boundary intersection position, a normal-guided interpolation method is used, and the position of the intersection on the boundary is calculated between a pair of adjacent corner points of the edge of the regular cubic grid cell according to the density values and positions of the two adjacent corner points and the direction of the three-dimensional spatial density field, and is taken as the position of the triangular face sheet vertex; In the isosurface extraction process, the density gradient change rate of each regular cubic grid cell is calculated, and when the density gradient change rate exceeds a preset density change threshold, the current regular cubic grid cell is subdivided, and the edge length of each regular cubic grid cell after subdivision is half of the edge length of the original regular cubic grid cell; Merge the triangular face sheets of all regular cubic grid cells to generate an initial three-dimensional terrain model.

2. The method for digital presentation of marine topography based on three-dimensional modeling according to claim 1, characterized in that, The step 2 specifically comprises: Map the spatial data in the multi-source marine data set to regular cubic grids in the same three-dimensional Cartesian coordinate system, and set the eight vertices of each regular cubic grid cell as the interpolation target points; For each target point to be interpolated, a first search radius is set, and original observation data points are searched within the first search radius, and a spatial Euclidean distance between each original observation data point and the target point to be interpolated is calculated; When the number of searched original observation data points is greater than or equal to a preset first number threshold, a first distance decay factor is set, and a weight value of each original observation data point is calculated; ； wherein, represents a weight value of the th original observation data point, represents a spatial Euclidean distance between the th original observation data point and the target point to be interpolated, represents a first distance decay factor; When the number of searched original observation data points within the first search radius is less than the preset first number threshold, the search radius is expanded to a second search radius, a second distance decay factor is set, the second distance decay factor is less than the first distance decay factor, and a weight value of a newly searched original observation data point in the second search radius is calculated using the second distance decay factor; After multiplying the values of all original observation data points by corresponding weight values respectively and summing, the value of the target point to be interpolated is obtained by dividing the sum of all weight values. The value of the target point to be interpolated is taken as a density value corresponding to the regular cubic grid unit. All regular cubic grid units are sequentially subjected to interpolation processing to generate a three-dimensional space density field.

3. The method for digital presentation of marine topography based on three-dimensional modeling according to claim 1, characterized in that, The step four is specifically: The vertex coordinates, corresponding normal vectors and curvature information of all triangular facets in the initial three-dimensional terrain model are extracted; For each triangular facet, an angle between normal vectors of adjacent triangular facets and the current triangular facet is calculated, and a triangular facet with an angle greater than a preset angle threshold is marked as a high curvature area; In each high curvature area, a boundary point set is extracted according to an envelope boundary, an average edge length between boundary points is calculated, and if the average edge length is greater than a preset refinement threshold, a new sampling point is introduced between the boundary points, and a sampling density is one-third of an original edge length; The boundary points and the newly introduced sampling points are collectively taken as an input point set, a Delaunay triangular net reconstruction method in a three-dimensional space is used to reconstruct the input point set, and a local triangular grid model satisfying a cavity maximization principle is generated; In the Delaunay reconstruction process, all non-convex units not satisfying a three-dimensional Delaunay condition are removed, and all effective units passing a tetrahedron circumscribed sphere test are retained to form a local fine model; The boundary regions of the local fine model and the initial three-dimensional terrain model are obtained, vertices in the boundary regions with a distance less than a preset fusion threshold are paired, and a connecting edge is generated by interpolation along a shortest path; The local fine model is inserted into the initial three-dimensional terrain model at a position corresponding to the original triangular facet, the original triangular facet is replaced, and a fused complete three-dimensional terrain model is generated.

4. The method for digital presentation of marine topography based on three-dimensional modeling according to claim 1, characterized in that, The step five is specifically: Satellite remote sensing images, hidden shooting images and hydrological profile images covering a region of the complete three-dimensional terrain model are obtained, and are unified to a same spatial coordinate system as the complete three-dimensional terrain model; Vertex coordinates of each triangular facet in the complete three-dimensional terrain model are extracted, and each vertex coordinate is projected to a two-dimensional plane of a corresponding image to establish a mapping relationship between the vertex of the complete three-dimensional terrain model and an image pixel, and a texture coordinate is generated. The satellite remote sensing image and the candid image are respectively subjected to texture mapping processing, image is selected as a texture source according to a preset priority, and texture mapping of all triangular facets is completed; After the texture mapping is completed, the contour image generated based on the hydrological profile data is taken as an additional texture layer, is superimposed on the surface of the complete three-dimensional terrain model which has been mapped, and a texture terrain model is formed.

5. The method for digital presentation of marine topography based on three-dimensional modeling according to claim 1, characterized in that, The LOD tree structure is constructed, and multiple levels of nodes are generated according to spatial region division, specifically as follows: Based on the spatial range of the texture terrain model, the entire three-dimensional terrain region is subjected to spatial hierarchical division, and a multi-level LOD tree structure is constructed in a top-down manner, wherein each level is gradually refined according to spatial resolution, the topmost level is a lowest-resolution model, and the bottommost level is a highest-resolution model; In the construction process, the spatial bounding box of the texture terrain model is taken as a root node region, the root node region is divided according to a preset spatial size, eight sub-regions in the form of an octree are respectively generated, and each sub-region corresponding resolution model data block is hung as a node in the LOD tree structure.

6. The method for digital presentation of marine topography based on three-dimensional modeling according to claim 1, characterized in that, The node information of the LOD tree structure includes: resolution level information for identifying the resolution accuracy level of the terrain data corresponding to the node; an error factor for representing the projection error of the current node resolution model and the last level resolution model in the screen space; spatial index information for uniquely identifying the spatial position of the node in the LOD tree structure, including the level number and spatial coordinates in the spatial division system; a sub-node pointer for recording the next level resolution model node list of the node, supporting on-demand loading and release; After the LOD tree structure is constructed, all nodes are connected according to the hierarchical relationship to form a complete multi-resolution data structure.

7. The method for digital presentation of marine topography based on three-dimensional modeling according to claim 1, characterized in that, The dynamic loading is realized by combining the view frustum clipping and the screen space error control strategy, specifically as follows: a view frustum structure is constructed, the view frustum structure is composed of six clipping planes, in sequence, a left clipping plane, a right clipping plane, an upper clipping plane, a lower clipping plane, a near clipping plane and a far clipping plane, and the three-dimensional space region visible by the current camera is limited through the view frustum structure; the bounding box information of each node in the LOD tree structure is traversed, and the spatial relationship with the six clipping planes of the view frustum structure is sequentially judged, if the bounding box is completely disjointed with any clipping plane, the node is judged as invisible and skipped for loading, if the bounding box intersects or partially intersects with all the clipping planes, the node is judged as a visible node; for each visible node, the projection size of the visible node on the screen is determined based on the resolution level boundary of the visible node, the diagonal length of the bounding box and the distance between the camera and the center of the node, and an error estimation value of the visible node in the screen space is obtained; the error estimation value is compared with a preset screen error tolerance, if the error estimation value is greater than the preset screen error tolerance, it is judged that the current accuracy is insufficient, and the child nodes of the visible node are loaded for further refinement and judgment, if the error estimation value is less than or equal to the preset screen error tolerance, or the visible node has no child node, the resolution model data corresponding to the visible node is loaded into a rendering queue.

8. The method for digital presentation of oceanographic relief based on three-dimensional modeling according to claim 1, characterized in that, The time controller specifically comprises a time axis module, a key frame index module and a node scheduling management module; The time axis module is used for constructing a unified time progress line, setting start and end time, time step and playing frame rate, supporting year, month, day and hour level time granularity control, and identifying terrain states in different periods with unified logical time sequence; The key frame index module is used for corresponding binding of multiple texture terrain models in different periods and the time axis, and recording the index relationship between each time point and the corresponding LOD tree node in a key value mapping mode, wherein each key frame corresponds to a complete LOD tree node state set, including node level, texture path and resolution level information; The node scheduling management module extracts the LOD tree node state of the target time point from the key frame index module according to the position of the current time axis pointer, performs unloading operation on the nodes needing to be replaced, performs loading operation on the added nodes, and adopts a cache mechanism to save the data blocks of the adjacent key frames in the node switching process; The time controller cooperates with the LOD tree structure to realize the ordered management and continuous display of the multi-time phase texture terrain model, and completes the time sequence visualization of the seafloor terrain evolution process.

Images