Terrain real scene modeling method, device and equipment and storage medium thereof

Abstract

The invention relates to a terrain real scene modeling method, device and equipment and a storage medium thereof. According to the method, multi-source heterogeneous data related to the terrain in a target area is integrated, a multi-modal data cube is generated through noise filtering and data fusion under a unified coordinate system, then a data cavity area is identified based on density clustering and elevation gradient analysis, and geophysical features are extracted. Innovatively constructing a neural implicit inference model fused with physical constraints to perform joint modeling on the earth surface and underground structures, and finally converting implicit terrain representation into an explicit three-dimensional terrain model; seamless fusion modeling of the earth surface and the underground space under the complex shielding environment is achieved, the physical rationality of the model and the hidden area reconstruction precision are remarkably improved, meanwhile, it is ensured that the generated model strictly conforms to the geological physical law, and the geological engineering safety and the planning reliability are effectively guaranteed.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of real scene modeling, and particularly relates to a terrain real scene modeling method, device, equipment and storage medium thereof. BACKGROUND

[0002] In the field of terrain real scene modeling, traditional methods mainly rely on laser radar scanning and oblique photography technology to construct a digital model of the ground surface, but these means have significant limitations when facing vegetation-covered areas, building-sheltered areas and underground spaces. The problem of data voids caused by ground surface obstructions has long been difficult to solve, and although conventional interpolation algorithms can fill geometric gaps, they cannot restore the true geological structure and physical properties. Especially for negative obstacles such as underground caves and pipeline channels, existing technologies often can only obtain fragmented data through destructive means such as drilling exploration, which not only has high costs but also is difficult to form a continuous spatial expression.

[0003] In recent years, non-invasive detection technologies such as ground penetrating radar have provided a new way for underground modeling, but how to effectively integrate surface observation data with underground detection information still faces challenges. Purely data-driven methods can establish geometric correlations, but often violate physical laws such as rock-soil mechanics and gravity field distribution, resulting in insufficient reliability of the generated model in actual engineering applications. For example, in urban underground space planning, distorted geological models can cause safety hazards; in archaeological exploration, lack of physical constraints can lead to misjudgment of cultural heritage.

[0004] Current research attempts to introduce physical laws into neural networks, but mostly focuses on a single field (such as only considering the gravity field or only analyzing electromagnetic reflections), and fails to establish a multi-physical field collaborative constraint mechanism across geology, gravity, and electromagnetism. In addition, the intelligent inference of existing methods for data void areas often deviates from actual geological conditions, resulting in systematic deviations between the generated negative obstacle model and the real terrain. This technical defect seriously restricts the development of key fields such as digital twin cities and geological disaster warning. SUMMARY

[0005] Therefore, it is necessary to provide a terrain real scene modeling method, device, equipment and storage medium thereof in view of the above technical problems.

[0006] In a first aspect, the present application provides a terrain real scene modeling method, comprising:

[0007] S1, collecting multi-source heterogeneous data related to the terrain in a target area, the multi-source heterogeneous data comprising surface image data, surface point cloud data and electromagnetic wave reflection data;

[0008] S2, performing noise filtering and ground point classification on the surface point cloud data to obtain a digital surface model; performing time domain filtering and gain adjustment on the electromagnetic wave reflection data to obtain standardized radar reflection data;

[0009] S3, performing coordinate system alignment on the ground surface image data, the digital surface model, and the standardized radar reflection data based on a spatial coordinate conversion operation to obtain a multi-modal data cube;

[0010] S4, performing density clustering on the multi-modal data cube to obtain a data void region, and performing elevation gradient analysis on the data void region to obtain a potential negative obstacle region;

[0011] S5, performing reflection waveform feature extraction on the potential negative obstacle region to obtain a geophysical feature;

[0012] S6, constructing a neural network model fused with physical constraints based on the geophysical feature and the multi-modal data cube, training the neural network model so that a prediction result meets both the observation data of the multi-modal data cube and a physical law consistent with the geophysical feature representation, and obtaining a neural implicit inference model; wherein the neural network model takes spatial coordinates as input and takes a signed distance function and a material attribute as output;

[0013] S7, obtaining an implicit terrain representation based on the neural implicit inference model, converting the implicit terrain representation into an initial three-dimensional grid, and performing surface reconstruction on the initial three-dimensional grid to generate an explicit three-dimensional terrain model containing ground surface structures and underground negative obstacles.

[0014] In a second aspect, the present application also provides a terrain real scene modeling device for implementing the method described in the first aspect, which comprises:

[0015] A data acquisition and integration module is configured to acquire multi-source heterogeneous data related to terrain in a target area, wherein the multi-source heterogeneous data includes ground surface image data, ground surface point cloud data, and electromagnetic wave reflection data.

[0016] A data preprocessing module is configured to perform noise filtering and ground point classification on the ground surface point cloud data to obtain a digital surface model, and perform time domain filtering and gain adjustment on the electromagnetic wave reflection data to obtain standardized radar reflection data.

[0017] A coordinate alignment module is configured to perform coordinate system alignment on the ground surface image data, the digital surface model, and the standardized radar reflection data based on a spatial coordinate conversion operation to obtain a multi-modal data cube.

[0018] A void and obstacle identification module is configured to perform density clustering on the multi-modal data cube to obtain a data void region, and perform elevation gradient analysis on the data void region to obtain a potential negative obstacle region.

[0019] A geophysical feature extraction module is configured to perform reflection waveform feature extraction on the potential negative obstacle region to obtain a geophysical feature.

[0020] a model construction and training module configured to construct a neural network model fusing physical constraints based on the geophysical features and the multi-modal data cube, and train the neural network model to make the prediction results consistent with the observation data of the multi-modal data cube and comply with the physical laws represented by the geophysical features, to obtain a neural implicit inference model, wherein the neural network model takes spatial coordinates as input and outputs signed distance functions and material attributes;

[0021] a three-dimensional modeling module configured to obtain an implicit terrain representation based on the neural implicit inference model, and convert the implicit terrain representation into an initial three-dimensional grid, and perform surface reconstruction on the initial three-dimensional grid to generate an explicit three-dimensional terrain model containing surface structures and underground negative obstacles.

[0022] In a third aspect, the present application also provides a computer device, including a memory and a processor, the memory stores a computer program, and the processor implements the terrain real scene modeling method in the first aspect when executing the computer program.

[0023] In a fourth aspect, the present application also provides a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to implement the terrain real scene modeling method in the first aspect.

[0024] The terrain real scene modeling method, device, equipment and storage medium described above integrate multi-source heterogeneous data related to terrain in a target area, generate a multi-modal data cube through noise filtering and data fusion in a unified coordinate system, identify data hollow areas and extract geophysical features based on density clustering and elevation gradient analysis, innovatively construct a neural implicit inference model fusing physical constraints to jointly model surface and underground structures, and finally convert the implicit terrain representation into an explicit three-dimensional terrain model, thereby realizing seamless fusion modeling of surface and underground space in a complex shielding environment, significantly improving the physical rationality of the model and the reconstruction accuracy of hidden areas, while ensuring that the generated model strictly complies with the laws of geophysics, effectively guaranteeing the safety and reliability of geological engineering planning. BRIEF DESCRIPTION OF DRAWINGS

[0025] In order to more clearly illustrate the technical solutions in the embodiments or the related art, the following will briefly introduce the drawings needed to be used in the embodiments or the related art description. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor based on these drawings.

[0026] Figure 1 A flowchart of a terrain real scene modeling method provided by the present application;

[0027] Figure 2 A flowchart of a process for constructing and training a neural implicit inference model in an alternative embodiment of the present application;

[0028] Figure 3 A structural schematic diagram of a terrain real scene modeling device provided by the present application. DETAILED DESCRIPTION

[0029] In order to make the purposes, technical solutions and advantages of the present application clearer, the present application is further described in detail below in combination with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and do not limit the present application.

[0030] REFERENCE Figure 1 It shows a flowchart of a terrain real scene modeling method provided by the present application, which comprises the following steps:

[0031] S1, collecting multi-source heterogeneous data related to terrain in a target area, the multi-source heterogeneous data comprising surface image data, surface point cloud data and electromagnetic wave reflection data.

[0032] Specifically, the surface image data can be obtained by high-resolution aerial photography or satellite remote sensing, which has texture and ground object information and can reflect the topographic features, land cover types, etc. of the ground. The surface point cloud data can be collected by a laser radar system, which can provide high-precision three-dimensional spatial coordinates and accurately record the geometric shapes of the ground and its attached objects. In the collection process, the laser radar emits laser pulses and receives reflected signals to measure the coordinates of a large number of ground points, forming dense point cloud data. The electromagnetic wave reflection data comes from ground penetrating radar and other devices, which obtain information about underground geological structures by emitting electromagnetic waves into the ground and receiving reflected signals. When the ground penetrating radar is working, the antenna transmits high-frequency electromagnetic wave pulses into the underground medium, and when the electromagnetic wave encounters the interface between different media, it produces reflection. The receiving antenna captures these reflected signals, and after processing, the cross-section information of the underground strata can be obtained.

[0033] When collecting the above multi-source data, according to the geological complexity of the target area, the target depth, the accuracy requirements, etc., appropriate device parameters and collection schemes are selected to ensure the reliability and integrity of the data, and to provide sufficient data support for subsequent model construction.

[0034] S2, noise filtering and ground point classification are performed on the surface point cloud data to obtain a digital surface model; time domain filtering and gain adjustment are performed on the electromagnetic wave reflection data to obtain standardized radar reflection data.

[0035] Specifically, for the ground point cloud data, the first step of noise filtering is a key step. The noise can come from various factors, such as measurement error, environmental interference, equipment precision limitation, etc. Statistical analysis-based methods can be used, such as calculating the neighborhood density of each point, and the points with density significantly lower than the surrounding area are identified as noise points for removal. Filtering algorithms such as Gaussian filtering, median filtering, etc. can also be used to smooth the point cloud data and remove isolated noise points. Based on the denoised point cloud data, ground point classification is performed to distinguish the ground part from the non-ground part (such as buildings, vegetation, etc.). This can be done with the help of a layer-by-layer classification algorithm, which first sets an initial ground seed point, and then gradually expands the classification outward based on the continuity of the ground points, slope changes, and other characteristics. Through this process, a digital surface model (DSM) can be constructed, which can truly reflect the relief form of the ground, including the top elevation information of all ground objects such as buildings and hills.

[0036] For electromagnetic wave reflection data, time domain filtering aims to remove useless high-frequency noise and low-frequency interference. For example, a Butterworth filter is used to filter the data, and the appropriate cutoff frequency is selected according to the spectral characteristics of the signal to retain the effective signal frequency band and remove noise components caused by poor device coupling, electromagnetic interference, etc. Gain adjustment is to compensate for the energy attenuation of electromagnetic waves during transmission, so that the reflection signal intensity at different depths is comparable. By analyzing the amplitude attenuation law of the reflection signal, a gain adjustment function is established to adjust the amplitude of the reflection signal at different times (corresponding to different depths), ensuring uniform signal intensity throughout the profile, thereby obtaining standardized radar reflection data to provide accurate data basis for subsequent geological structure analysis.

[0037] S3, based on the spatial coordinate conversion operation, aligning the coordinate systems of the ground image data, the digital surface model, and the standardized radar reflection data to obtain a multi-modal data cube.

[0038] Specifically, when aligning the coordinate systems, a unified reference coordinate system is first determined, which can be selected as the geodetic coordinate system (such as WGS-84 or CGCS2000) as the basis. For ground image data, its original coordinate system can be the row and column coordinates of the image itself. Through ground control points (GCPs) or known geographic coordinate information, geometric correction methods such as polynomial correction and affine transformation are used to convert the pixel coordinates of the image data to geographic coordinates. At the same time, combined with the elevation information of the digital surface model, the image is terrain corrected to eliminate the image distortion caused by the terrain undulation, so that the image data accurately reflects the true position and form of the ground.

[0039] The digital surface model itself has three-dimensional spatial coordinates in the generation process, but needs to be matched with image data and radar reflection data for precision and coordinate refinement. For standardized radar reflection data, the coordinate information is provided by the positioning system of ground penetrating radar (such as GPS), but there may be some errors and deviations. Therefore, the radar reflection data line position is matched with the corresponding position in the ground surface image and the DSM, and the two-dimensional line information of the radar reflection data is extended to three-dimensional spatial coordinates through spatial interpolation and other methods, so that it is aligned with the ground data in space. After completing the spatial coordinate conversion and alignment of each data, these multi-source data are integrated in a unified three-dimensional space to build a multi-modal data cube.

[0040] The data cube organically integrates the texture information of the ground image, the elevation information of the DSM, and the underground geological information of the radar reflection in the spatial dimension, forming a comprehensive data structure containing rich spatial and attribute information, providing a data support platform for subsequent data analysis and model construction.

[0041] S4, performing density clustering on the multi-modal data cube to obtain a data hollow area; performing elevation gradient analysis on the data hollow area to obtain a potential negative obstacle area.

[0042] Specifically, in the multi-modal data cube, the data distribution may be uneven, and some areas may form data hollows due to shielding, undetectable, and other reasons. The density clustering algorithm (such as the DBSCAN algorithm) divides clusters based on the local density difference of data points. First, set clustering parameters such as neighborhood radius and minimum number of contained points. In the space of the data cube, calculate the number of points contained in the neighborhood of each data point. If the number of points is greater than or equal to the minimum number of contained points, the point is marked as a core point, and from the core point, the neighborhood is continuously expanded to merge the points in the neighborhood that meet the conditions into a cluster. Through this process, the data dense area can be divided into multiple clusters, and the isolated points or low-density areas that are not clustered are identified as data hollow areas.

[0043] For these data void regions, further elevation gradient analysis is performed. The elevation gradient reflects the degree of fluctuation of the surface or underground structure in space. In the three-dimensional space of the data cube, the rate of change of the elevation value around the void region is calculated, that is, the gradient is obtained by calculating the ratio of the elevation difference between adjacent points to the horizontal distance. If the elevation gradient around a certain void region presents significant changes, such as a sudden increase in gradient or the presence of gradient anomalies, and such changes conform to the typical geological characteristics of negative obstacles (such as underground caves, pipeline channels, etc.) (for example, the elevation difference between the top of the cave and the surrounding rock-soil, the linear gradient characteristics of the pipeline channel, etc.), it can be judged that this region is a potential negative obstacle region. This analysis process comprehensively considers the geological structure law, the physical mechanism of the formation of negative obstacles, and combines with the existing geological prior knowledge to reasonably explain and judge the abnormal characteristics of the elevation gradient, so as to accurately identify the location where the negative obstacle may exist.

[0044] S5, performing reflection waveform feature extraction on the potential negative obstacle region to obtain geophysical features.

[0045] Specifically, the reflection waveform feature is the time-domain signal feature of the electromagnetic wave after reflection at the interface of different underground media, which contains geological information. In the potential negative obstacle region, first, the extracted radar reflection waveform is preprocessed, such as removing the direct current offset, performing zero-mean processing, etc., to eliminate the pseudo-signal components in the data and enhance the features of the true reflection signal. Then, feature extraction is performed from both time and frequency domains.

[0046] In terms of time domain features, the amplitude information of the reflection wave can be extracted, including peak amplitude, valley amplitude, and amplitude attenuation coefficient, etc. The size of the amplitude is closely related to the electrical difference of the underground medium, the fluctuation form of the reflection interface, and other factors. For example, when the electromagnetic wave reflects from a high-conductivity medium, a strong reflection signal may be generated; and when the reflection interface is uneven, the amplitude of the reflection wave will also show corresponding variation characteristics.

[0047] At the same time, the phase characteristics of the reflection waveform can be analyzed, such as zero phase, continuity of the phase axis, etc. The phase information reflects the propagation path of the reflection wave and the wave characteristics of the medium, and the discontinuity of the phase axis may indicate the existence of a negative obstacle (such as a cave that causes the reflection wave path to change, resulting in a broken phase axis).

[0048] In addition, the time delay feature of the reflection wave can also be extracted, that is, the reflection waveform is cross-correlated with the reference waveform to obtain the time difference of arrival of the reflection wave, which is related to the depth position of the negative obstacle, because the propagation speed of the electromagnetic wave in media at different depths is different, resulting in a change in the time delay of the reflection signal.

[0049] In the aspect of frequency domain feature extraction, the reflected waveform is subjected to fast Fourier transform (FFT) to obtain its spectral characteristics. The main frequency component, frequency bandwidth, and amplitude distribution of the reflected wave spectrum are analyzed. Different geological media have different absorption and scattering characteristics of electromagnetic waves, which leads to changes in the spectral characteristics of the reflected wave. For example, certain negative obstacles (such as cave fillings containing specific minerals) can produce strong absorption or scattering of electromagnetic waves of specific frequencies, thereby causing the reflected wave spectrum to exhibit characteristics such as amplitude reduction or frequency band broadening at the corresponding frequencies. By comprehensively analyzing the reflected waveform features in the time domain and the frequency domain, a comprehensive geophysical feature set is formed, which provides a key basis for subsequent geological modeling and obstacle identification.

[0050] S6, based on the geophysical features and the multi-modal data cube, a neural network model with physical constraints is constructed; by training the neural network model, the prediction result meets the observation data of the multi-modal data cube and the physical law represented by the geophysical features at the same time, and a neural implicit inference model is obtained; wherein the neural network model takes spatial coordinates as input and symbolic distance function and material properties as output.

[0051] Specifically, first, the basic architecture of the neural network is determined, which can adopt a multi-layer perceptron (MLP) structure, including an input layer, several hidden layers, and an output layer. The input layer receives spatial coordinate information (such as three-dimensional Cartesian coordinates), which covers position points at different depths from the surface to the underground. The hidden layer adopts a nonlinear activation function (such as ReLU or Tanh function), and learns the complex patterns in the data through the connection weights and biases between neurons. The output layer outputs the symbolic distance function (SDF) and the material properties. The symbolic distance function is used to represent the distance from any point in space to the boundary of the negative obstacle. When the point is inside the obstacle, the distance is negative; when it is outside, the distance is positive, and the value reflects the nearness or remoteness of the point to the boundary. The material properties are used to describe the geological material categories at the point (such as rock layers, soil layers, caves, pipelines, etc.).

[0052] In constructing the network, the geophysical features and the observation data in the multi-modal data cube are integrated into the network structure as prior constraints. For example, the extracted reflection waveform features are taken as auxiliary input or feature fusion is performed in the intermediate layer of the network, so that the network can learn the geological patterns consistent with the geophysical laws. At the same time, in the loss function of the network, in addition to the error term between the predicted signed distance function and the actual observation data (such as known geological drilling data, partially interpreted negative obstacle boundaries, etc.), a physical constraint term is introduced. These physical constraint terms are based on physical laws such as geomechanics and electromagnetism, such as requiring the material property distribution output by the network to comply with the stratum deposition law (new on the old), and the correlation between the reflection wave features and the material properties to comply with the electromagnetic scattering principle. By reasonably designing the weight coefficients of the loss function, the influence of the data consistency term and the physical law term on the network training is balanced, so that the network can not only fit the observation data, but also comply with the geophysical laws in the learning process, thereby improving the reliability and generalization ability of the model.

[0053] In the training process, an optimization algorithm (such as Adam or SGD algorithm) is used to iteratively update the parameters of the network. In each iteration, a certain number of samples (spatial coordinate points and their corresponding true labels, such as known information obtained from multi-modal data cubes and geophysical features) are selected, the loss value between the forward propagation output of the network and the true label is calculated, and then the gradient is calculated through the back propagation algorithm to update the network weights. During the training process, the change curve of the loss value is continuously monitored, and when the loss value tends to be stable and the error on the validation set no longer significantly decreases, it is considered that the network training is completed, and the neural implicit inference model is obtained. This model can quickly and accurately predict the signed distance function and material properties at any given spatial coordinates, providing an efficient inference tool for subsequent three-dimensional terrain model construction.

[0054] S7, obtaining an implicit terrain representation based on the neural implicit inference model; converting the implicit terrain representation into an initial three-dimensional grid, and performing surface reconstruction on the initial three-dimensional grid to generate an explicit three-dimensional terrain model containing surface structures and underground negative obstacles.

[0055] Specifically, the implicit terrain representation output by the neural implicit inference model is a continuous function representation, which mathematically describes the terrain (including the surface and the underground part) of the entire target area based on the signed distance function and the material properties. The signed distance function defines the distance relationship between any point in space and the negative obstacle boundary, and the material properties explicitly specify the geological material category of each point. This implicit representation has high mathematical continuity and computability, which is convenient for geometric analysis and model conversion.

[0056] The process of converting the implicit terrain representation into an initial three-dimensional mesh mainly relies on the Marching Cubes algorithm or other isosurface extraction algorithms. First, according to the required model accuracy, the target area is divided into a regular three-dimensional voxel grid. Then, within each voxel unit, the signed distance function value of the boundary points is calculated by interpolation, and the isosurface with a signed distance function value of zero, i.e., the boundary of the negative obstacle, is found. According to the connection of the isosurface in different voxels, triangular or quadrilateral mesh patches are constructed, and the initial three-dimensional mesh model is formed. At this time, the initial three-dimensional mesh may contain a large number of mesh patches, and there are problems such as uneven mesh quality and complex topological structure, which need to be further optimized by surface reconstruction.

[0057] The purpose of surface reconstruction is to improve the quality and practicality of the mesh model and generate an explicit three-dimensional terrain model that meets geometric constraints and physical significance. First, mesh simplification can be performed, which can use a simplification algorithm based on Quadric error metric or a wavelet simplification method to remove redundant patches in the mesh, reduce the complexity of the model, and at the same time preserve key terrain features and obstacle boundary information. Then, mesh smoothing is performed, which uses Laplacian smoothing or Taubin smoothing algorithm to eliminate sharp protrusions and depressions on the mesh surface, and improve the visual effect and numerical stability of the mesh. Then, topological repair is performed to identify and repair topological errors such as non-manifold edges and hanging edges in the mesh, ensuring that the mesh has a reasonable topological structure and can accurately represent the connectivity of the ground surface and underground negative obstacles. Finally, combined with the ground image data and digital surface model, the model surface is texture mapped and material properties are assigned, so that the generated explicit three-dimensional terrain model is not only accurate in geometric shape, but also visually reflects the geological features and material distribution of the ground surface and underground.

[0058] Through the above conversion and optimization process, a high-quality explicit three-dimensional terrain model containing complete ground surface structure and underground negative obstacles is finally obtained, providing intuitive and accurate terrain reference for geological engineering, urban planning, archaeological exploration, and other fields.

[0059] The above-mentioned terrain real scene modeling method integrates multi-source heterogeneous data related to the terrain in the target area, generates a multi-modal data cube through noise filtering and data fusion in a unified coordinate system, and then identifies data void areas and extracts geophysical features based on density clustering and elevation gradient analysis. A neural implicit inference model that integrates physical constraints is innovatively constructed to jointly model the ground surface and underground structure, and finally converts the implicit terrain representation into an explicit three-dimensional terrain model. This method realizes seamless fusion modeling of the ground surface and underground space in a complex occlusion environment, significantly improves the physical reasonableness and reconstruction accuracy of hidden areas, and ensures that the generated model strictly conforms to the laws of geophysics, effectively ensuring the safety of geological engineering and the reliability of planning.

[0060] In an alternative embodiment, S4 comprises the following steps:

[0061] S41, performing density clustering on the ground point cloud data in the multi-modal data cube to generate a set of density anomaly regions.

[0062] Specifically, density clustering is an algorithm for grouping data based on the difference in data density, the core of which is to group points with similar density into the same group, while regions with large density differences are identified as anomalies. In this step, an improved algorithm based on DBSCAN (Density-Based Spatial Clustering of Applications with Noise) can be used. First, two key parameters are set: neighborhood radius ε and minimum neighborhood point number MinPts. The neighborhood radius ε determines the range of measuring the neighborhood of a point in a multi-dimensional space, while MinPts is used to determine whether a point is a core point. A core point is a point whose ε neighborhood contains at least MinPts points. For ground point cloud data, its three-dimensional coordinates (x, y, z) constitute the feature space for clustering analysis.

[0063] In the clustering process, the algorithm traverses each unmarked point and calculates the number of points in its ε neighborhood. If the number of points in the neighborhood of a point is greater than or equal to MinPts, the point is marked as a core point and the cluster is expanded with it as the center, and all points in the neighborhood are grouped into the same cluster. If the points in the neighborhood are also core points, their neighborhoods are further expanded until they cannot be expanded.

[0064] In this way, regions with high density are clustered into multiple clusters, while regions with low density (such as data voids or anomaly regions) are treated as noise points or small clusters. Finally, these noise points and small cluster regions that are not clustered into the main clusters are collected to generate a set of density anomaly regions. These density anomaly regions may correspond to special topographic structures on or near the ground surface, such as collapse areas, cave entrances, and areas above underground pipelines, providing a basis for further analysis.

[0065] S42, performing elevation gradient calculation on the set of density anomaly regions to generate an elevation gradient distribution map for each region.

[0066] Specifically, the elevation gradient reflects the rate of change of the terrain in space, which is an important indicator for identifying terrain mutations or potential geological structures. For each density anomaly region, an elevation gradient calculation method based on triangulation is used. First, the point cloud data in the region is Delaunay triangulated to form a series of irregular triangles. Each triangle is composed of three adjacent points, and the vertex contains the elevation information (z coordinate). For each triangle, calculate the elevation change rate of its three edges, that is, get the gradient value by calculating the ratio of the elevation difference between two points and the horizontal distance. Then, assign the gradient value of each edge to the midpoint position of the edge to form gradient sampling points.

[0067] To generate the elevation gradient distribution map, the inverse distance weighted (IDW) interpolation method is used. A regular grid is constructed in the region, and the grid cell size is determined according to the required accuracy. For each grid point, calculate the weighted average value of the gradient sampling points around it, and the weight is the inverse distance square of the sampling point to the grid point.

[0068] In this way, the discrete gradient sampling point values are interpolated to the entire region to generate a continuous elevation gradient distribution map. The distribution map represents the gradient size in color or grayscale, which intuitively shows the elevation change in the region. High gradient areas may correspond to steep slopes, cliffs, collapse edges, or the top boundaries of underground caves, providing key evidence for subsequent identification of concave regions.

[0069] S43, based on the elevation gradient distribution map, identify the concave region by calculating the normalized gradient difference.

[0070] Specifically, the normalized gradient difference is an indicator that measures the relative degree of change in the region's elevation gradient, used to identify concave regions that have significant differences from the surrounding environment. The specific calculation method is as follows: First, determine the size of the analysis window on the elevation gradient distribution map, and you can choose a 3x3 or 5x5 moving window. For each window center point, calculate the difference between its gradient value and the average value of all points in the window, that is, ΔG=G_center-G_mean. Then, divide this difference by the standard deviation σ_G of the gradient values in the window to get the normalized gradient difference NGD=ΔG / σ_G. By setting a threshold (such as NGD<-1.5), filter out points with significantly lower gradients than the surrounding area, which may belong to the boundaries or interiors of the concave region.

[0071] To further identify the depression regions, a region growing algorithm is adopted. Starting from the seed points whose NGD satisfies the threshold condition, the region is gradually expanded according to the predefined growing rules (e.g., the difference of NGD between adjacent points is less than the threshold, and the gradient directions are consistent). During the growing process, the statistical features (e.g., average gradient, area, etc.) of the region are updated in real time, and the maximum growing range is set to avoid over-expansion. Finally, multiple candidate depression regions are obtained. These depression regions may correspond to the low-lying areas on the ground surface, collapsed pits, partially collapsed areas of underground caves, etc., providing a basis for subsequent terrain continuity analysis.

[0072] S44, performing terrain continuity analysis on the depression region by calculating the curvature variation rate of adjacent grids in the depression region, screening abnormal isolated regions as a candidate negative obstacle region set.

[0073] Specifically, the curvature variation rate reflects the change of the bending degree of the terrain in a local range, and is an important indicator for evaluating terrain continuity and identifying abnormal terrain structures. For each depression region, it is first divided into regular grids, and the grid size is determined according to the required accuracy. Then, the curvature value of each grid cell is calculated. The curvature calculation is based on the method of fitting a quadratic surface: in the neighborhood of each grid point (such as a 3x3 grid window), a quadratic surface z=ax 2 +by 2 +cxy+dx+ey+f is fitted, and the coefficients a, b, c, d, e, f are solved by the least squares method. The curvature value is calculated from the quadratic coefficients a and b, and the specific formula is The greater the curvature value, the higher the bending degree of the terrain at that point.

[0074] The calculation method of the curvature variation rate is: for each grid point, the relative difference between its curvature value and the average curvature value in the surrounding neighborhood (such as a 3x3 window) is calculated, i.e., CVR=(κ_center-κ_mean) / κ_meanx100%. Wherein, CVR represents the curvature variation rate, κ_center represents the curvature value at the center grid point, and κ_mean represents the average curvature value in the neighborhood (such as a 3x3 grid window) around the center grid point.

[0075] By setting a threshold (e.g. CVR > 50%), regions with sharp curvature changes are identified, which can correspond to abrupt points of the terrain, such as sudden collapse areas of cave ceilings, local uplifts of underground pipelines, etc. Based on the screening results of the curvature change rate, combined with the prior knowledge of the continuity of the terrain (e.g. the surface should have a relatively smooth transition, and the underground structure should conform to the geological layer distribution rule), morphological analysis methods (such as opening operation, closing operation) are used to remove regions with good continuity with the surrounding terrain, and to retain regions that are abnormally isolated and have sharp curvature changes in the terrain, forming a candidate negative obstacle region set. These candidate regions can include the entrance of an underground cave, a collapse area not covered by vegetation, an exposed part of an underground pipeline, etc., providing a basis for further shelter type verification.

[0076] S45, based on the optical texture data in the multi-modal data cube, excluding vegetation covered areas by HSV color space segmentation, performing shelter type verification on the candidate negative obstacle region set to obtain verified potential negative obstacle regions.

[0077] Specifically, the optical texture data contains surface color and texture information, which can effectively distinguish different ground object types. In this step, HSV (Hue, Saturation, Value) color space is used for vegetation covered area segmentation. First, convert the RGB optical image data to HSV color space. In the HSV space, the hue value of vegetation is usually in the green range (about 80°-140°), the saturation is high (>0.4), and the value is moderate (0.2-0.8). According to these characteristics, set the threshold range, and extract the vegetation covered area by color segmentation algorithm. In the segmentation process, adaptive threshold optimization based on Otsu method is used to improve the segmentation accuracy.

[0078] After excluding the vegetation covered area, shelter type verification is performed on the candidate negative obstacle region set. The shelter type includes non-vegetation shelter caused by buildings, artificial structures, shadows, etc. For each candidate region, analyze its features in the optical texture data: building shelter areas usually have regular geometric shapes, uniform textures, and specific colors (such as red, gray, etc.); artificial structures (such as roads, bridges) have linear or block texture characteristics; shadow areas are characterized by low brightness, low saturation, and no obvious texture features. By constructing a corresponding feature discrimination model (such as a support vector machine-based classifier), the candidate regions are classified to exclude non-negative obstacle shelter areas. Finally, the verified potential negative obstacle regions are obtained, which are obvious depressions, discontinuities or abnormal structures on the surface or near-surface, and are not covered by vegetation or other non-geological factors, providing accurate input for further geophysical feature extraction and modeling.

[0079] In an alternative embodiment, S5 comprises the following steps:

[0080] S51, performing a ground penetrating radar B-scan profile extraction operation on the potential negative obstacle region to obtain a reflection waveform dataset.

[0081] Specifically, first determine the working parameters of the ground penetrating radar, including center frequency, antenna spacing, sampling frequency, etc. The selection of the center frequency can be determined according to the target depth and the electromagnetic properties of the geological medium, high frequency (such as 1GHz) is suitable for shallow high resolution detection, and low frequency (such as 400MHz) is suitable for deep detection. The antenna spacing should be adjusted according to the geological structure density of the target area, which can be between 0.5-2 meters, to ensure the spatial continuity of the data.

[0082] In the potential negative obstacle region, grid-like survey lines are laid out according to the preset line spacing (which can be 0.5-1 meters). The ground penetrating radar antenna moves uniformly along the survey line, emits electromagnetic wave pulses and receives reflected signals. The reflection signals at each survey line position are digitized to form an A-scan (reflection waveform perpendicular to the time axis). A plurality of adjacent A-scans are arranged and combined according to the survey line position to form a B-scan profile, each B-scan profile contains underground reflection information, reflecting the reflection characteristics of different depths of geological interfaces and target bodies.

[0083] The extraction process of the reflection waveform dataset includes preprocessing and feature extraction of the B-scan profile. First, denoising the original B-scan data, using wavelet transform or adaptive filtering algorithm to remove high frequency noise and power frequency interference. Then, time zero correction is performed to align the start time of the reflection signal with the electromagnetic wave emission time, ensuring the accuracy of the depth calculation. Next, through the amplitude recovery algorithm to compensate the energy attenuation of electromagnetic wave in the propagation process, enhance the visibility of weak reflection signal. Finally, the processed B-scan profile is organized according to the spatial position and depth information to form the reflection waveform dataset, which provides the basis for the subsequent reflection feature parameter extraction.

[0084] S52, performing adaptive threshold segmentation on the reflection waveform dataset to obtain reflection feature parameters.

[0085] Specifically, adaptive threshold segmentation is a method of dynamically determining the segmentation threshold according to the local statistical characteristics of the data, which can effectively distinguish between valid reflection signals and background noise. First, the reflection waveform dataset is divided into blocks, and the entire dataset is divided into multiple local blocks, each block containing a certain number of reflection waveforms (such as a 10x10 waveform matrix). In each local block, the local mean and standard deviation of the reflection waveform are calculated, and the segmentation threshold is dynamically determined by the formula T=μ+kσ, where μ is the local mean, σ is the standard deviation, and k is an empirical coefficient, which can be taken as 2-3.

[0086] For each reflected waveform, compare its amplitude value with the corresponding local threshold T. If the amplitude value is greater than T, it is judged as an effective reflection signal; otherwise, it is regarded as background noise. In this way, the reflected waveform is divided into an effective signal segment and a noise segment. The start point and end point of the effective signal segment correspond to the arrival time and duration of the reflected wave, respectively, reflecting the burial depth and horizontal size of the underground target body. At the same time, the peak amplitude, waveform width, rise time and other parameters of the effective signal segment are calculated, which can represent the morphological characteristics of the reflected wave and the physical properties of the underground geological structure. For example, a high-amplitude reflection may indicate a high-conductivity medium interface or a sharp boundary of a geological body; a wide waveform may correspond to a thick-layer medium or multiple reflections of a rough interface. By statistically analyzing these reflection characteristic parameters, a set of reflection characteristic parameters is formed, which provides data support for subsequent geophysical feature fusion.

[0087] S53, collecting gravity field distribution data of the point positions corresponding to the potential negative obstacle area by the gravity measuring instrument; performing Bouguer anomaly calculation on the gravity field distribution data to obtain an underground density anomaly distribution map.

[0088] Specifically, first, a suitable gravity measuring instrument is selected, such as a proton magnetometer or a superconducting gravity meter. In the potential negative obstacle area, gravity measurement points are arranged according to a preset grid spacing (for example, 5-20 meters). At each point, the instrument is stably placed and multiple repeated measurements are performed to eliminate the influence of instrument drift and environmental interference. The spatiotemporal information of the measured gravity acceleration data is recorded, including the point position coordinates, measurement time and environmental parameters (such as air temperature, air pressure, etc.).

[0089] Bouguer anomaly calculation is a process of topographic and density correction on gravity field measurement data, aiming to eliminate the interference of surface topography and known density medium on the gravity field, and extract underground density anomaly information. The specific steps are as follows: first, the original gravity acceleration data is corrected for latitude to eliminate the gravity changes caused by earth rotation and latitude changes; then, the gravity value of the measurement point is reduced to the same reference height (such as sea level) by performing altitude correction; then, the gravity influence of the surrounding topography on the measurement point is calculated and deducted by performing topographic correction; finally, the gravity contribution of known density medium (such as surface loose deposits) is eliminated by performing density correction according to known geological data and surface rock density.

[0090] The gravity anomaly data after the above correction is the Bouguer anomaly value. The Bouguer anomaly value is interpolated to a regular grid according to the point position coordinates to generate an underground density anomaly distribution map. The distribution map displays the spatial distribution of underground density anomalies in the form of contour lines or pseudo-color images, reflecting the density differences of underground geological structures, such as low-density anomalies that may correspond to caves, fractures or low-density rock layers, and high-density anomalies that may indicate ore bodies, dikes or dense rock layers.

[0091] S54, fuse the reflection feature parameters and the underground density anomaly distribution map to obtain geophysical features.

[0092] Specifically, first, the reflection feature parameters (such as amplitude, waveform width, arrival time, etc.) are aligned in space with the underground density anomaly distribution map. Through a geographic registration algorithm, the reflection feature parameter set and the density anomaly distribution map are converted to the same coordinate system and spatial resolution, ensuring the correspondence of the two in position and depth.

[0093] The fusion method can use a linear combination based on weights or a nonlinear mapping based on machine learning. In the linear combination method, according to prior geological knowledge and data sensitivity analysis, weight coefficients are assigned to the reflection feature parameters and the density anomaly. For example, if the target area is dominated by electromagnetic reflection features, the reflection feature parameters are given a higher weight; if the density anomaly is more critical for identifying the target body, the weight of the density anomaly is increased. Through weighted summation, a comprehensive geophysical feature map is generated.

[0094] In the machine learning method, a training data set is constructed, including reflection feature parameters, density anomaly values, and their corresponding geological attribute labels of known geological structures. Using algorithms such as support vector machines (SVM) and random forests (RF), a classification or regression model is trained to learn the complex mapping relationship between reflection features and density anomalies and geological attributes. Through model prediction, the reflection features and density anomalies of unknown areas are converted into geophysical features, achieving comprehensive representation of underground geological structures.

[0095] The fused geophysical features not only contain electromagnetic reflection information but also incorporate density anomaly constraints, which can more accurately reflect the physical properties and spatial distribution of underground geological structures. For example, reflection feature parameters can indicate the position and shape of geological interfaces, while density anomaly distribution maps provide density information of geological bodies, and the combination of the two can effectively distinguish different types of negative obstacles (such as cavities, karst, ore bodies, etc.), providing more reliable basis for subsequent modeling and inference. In this way, the obtained geophysical feature set can comprehensively and accurately describe the underground geological characteristics of potential negative obstacle regions, laying a solid foundation for constructing high-precision three-dimensional terrain models.

[0096] Reference Figure 2 In an optional embodiment, S6 includes the following steps:

[0097] S61, performing spatial coordinate resolution on the multi-modal data cube to obtain a discrete spatial coordinate set; performing parameter deconstruction on the geophysical features to obtain a geological constraint matrix, a gravity difference threshold, and an electromagnetic similarity threshold; the geological constraint matrix represents the mapping relationship of lithology physical parameters.

[0098] Specifically, in the multi-modal data cube, the surface image data, the surface point cloud data, and the normalized radar reflection data have been aligned to a unified coordinate system. The purpose of spatial coordinate resolution is to extract the spatial information in these data and form a discrete spatial coordinate set. Each coordinate point (x i ,y i ,z i ) represents a specific geographic location and elevation information. Specifically, the surface image data provides two-dimensional geographic coordinates (x, y), while the surface point cloud data and the radar reflection data provide elevation information z. By interpolation and sampling, these data are integrated into a regular three-dimensional grid, and each grid point is a discrete spatial coordinate point.

[0099] The geophysical features include the reflected waveform characteristics and the underground density anomaly distribution map. The extraction of the reflected waveform characteristic parameters is obtained by analyzing the ground penetrating radar B-scan profile, including amplitude, waveform width, rise time, etc. The underground density anomaly distribution map is obtained by gravity measurement and Bouguer anomaly calculation. The geological constraint matrix is a multi-dimensional array, and its elements represent the mapping relationship between different lithologies, such as density, electromagnetic characteristics, elastic parameters, etc. Specifically, the geological constraint matrix can be represented as:

[0100]

[0101] where ρ i represents the density of the i-th lithology, μ i represents the electromagnetic characteristics of the i-th lithology, and so on. The gravity difference threshold and the electromagnetic similarity threshold are set according to the statistical analysis and experience of the actual measurement data, and are used to evaluate the matching degree between the simulation results and the measured data.

[0102] S62, a neural network architecture containing a symbolic distance prediction branch and a material attribute prediction branch is constructed, and the discrete spatial coordinate set is input into a position encoder for feature enhancement to obtain a high-dimensional feature vector set; the discrete spatial coordinate set is input into the position encoder to obtain a high-dimensional feature vector set.

[0103] Specifically, the neural network architecture includes two main branches, namely the symbolic distance prediction branch and the material attribute prediction branch. The symbolic distance prediction branch is used to predict the distance of any point in space to the negative obstacle boundary, and the material attribute prediction branch is used to predict the geological material category of the point. The input of the neural network is the discrete spatial coordinate set, which is enhanced by the position encoder to generate a high-dimensional feature vector set.

[0104] The position encoder is used to convert spatial coordinates into high-dimensional feature vectors to enhance the model's ability to learn spatial information. The position encoder can adopt a multi-layer perceptron (MLP) structure to map low-dimensional coordinates to high-dimensional space through nonlinear transformation.

[0105] The purpose of feature enhancement is to enrich the feature vector by introducing additional feature information such as local density, elevation gradient, etc. For example, local density information and elevation gradient information can be input as additional inputs along with spatial coordinates into the position encoder. This approach can improve the model's ability to model complex terrain.

[0106] S63, input the high-dimensional feature vector set into the neural network to perform forward propagation calculation to generate predicted symbol distance values and predicted material property values.

[0107] Specifically, forward propagation calculation is the core process of neural networks, which generates the final prediction result through layer-by-layer calculation of the input layer, hidden layer, and output layer. The output layer of the symbol distance prediction branch can be a neuron for predicting the signed distance of each spatial point to the negative obstacle boundary. The output layer of the material property prediction branch can be a multi-classification layer for predicting the geological material category of each spatial point. The activation function of the output layer can use the softmax function to convert the output value into a probability distribution.

[0108] S64, reconstruct the three-dimensional terrain geometry based on the predicted symbol distance values, and calculate the geological stability error by using the lithology physical parameters in the geological constraint matrix and the rock-soil mechanics equilibrium criterion.

[0109] Specifically, the predicted symbol distance values can be used to reconstruct the three-dimensional terrain geometry, and the Marching Cubes algorithm can be used to extract the zero level surface of the signed distance function to obtain the surface model of the terrain. The basic idea of the Marching Cubes algorithm is to divide the three-dimensional space into regular cubic grids, and then determine whether the terrain surface passes through each cubic according to the signed distance value of each cubic vertex. If it does, generate the corresponding triangular patches.

[0110] The geological stability error is calculated by the rock-soil mechanics equilibrium criterion to evaluate the error of the reconstructed three-dimensional terrain model in terms of geological stability. The rock-soil mechanics equilibrium criterion includes the Mohr-Coulomb criterion and the Drucker-Prager criterion. The specific calculation formula can be expressed as:

[0111]

[0112] where σ i is the calculated stress of the i-th point, σ critical,i is the critical stress value obtained from the geological constraint matrix, and N is the total number of points.

[0113] S65, reconstruct the density distribution of the subsurface medium based on the predicted material property values, obtain a simulated gravity field; when the deviation value between the measured gravity field and the simulated gravity field exceeds the gravity difference threshold, the deviation value between the measured gravity field and the simulated gravity field is taken as the gravity matching error.

[0114] Specifically, according to the density information in the geological constraint matrix, the material property of each spatial point is converted into a density value. The simulated gravity field can be calculated by numerical methods such as finite difference method or finite element method. According to Newton's law of gravity, the calculation formula of the gravity field is:

[0115]

[0116] where G is the gravitational constant, ρ i is the density of the i-th point, ΔV i is the volume element of the i-th point, r is the position of the calculation point, r i is the position of the i-th point.

[0117] The gravity matching error is obtained by comparing the simulated gravity field and the measured gravity field. The specific calculation formula is:

[0118]

[0119] where g sim,i is the simulated gravity field value, g meas,i is the measured gravity field value, and M is the number of measurement points. When E gravity exceeds the set gravity difference threshold, it is taken as the gravity matching error.

[0120] S66, reconstruct the electromagnetic parameter distribution of the medium based on the predicted material property values, obtain a simulated wave form; when the similarity between the simulated wave form and the measured wave form is lower than the electromagnetic similarity threshold, calculate the wave form difference between the simulated wave form and the measured wave form as the electromagnetic reflection error.

[0121] Specifically, according to the electromagnetic parameter information in the geological constraint matrix, the material property of each spatial point is converted into an electromagnetic parameter value. The reconstruction of electromagnetic parameter distribution can improve the accuracy of reconstruction by introducing electromagnetic tomography technology and inversion algorithm.

[0122] In actual operation, finite difference time domain method (FDTD) or finite element method (FEM) can be used to reconstruct the electromagnetic parameter distribution. FDTD method can effectively simulate the propagation and reflection of electromagnetic wave in medium by discretizing the derivatives of Maxwell equations in time and space. FEM method can handle complex geometric structures and boundary conditions by dividing the calculation region into a finite number of elements and approximately solving Maxwell equations in each element.

[0123] The simulation waveform can be calculated by numerical methods such as FDTD or FEM. According to Maxwell's equations, the propagation and reflection of electromagnetic waves are affected by the electromagnetic parameters of the medium. During the calculation, appropriate boundary conditions and initial conditions are set to simulate the real geological environment and detection conditions.

[0124] The electromagnetic reflection error is obtained by comparing the simulation waveform and the measured waveform. The waveform similarity can be evaluated by correlation coefficient or mutual information and other indicators. Specifically, first, align the simulation waveform and the measured waveform on the same time axis, and then calculate the similarity between them. If the similarity is lower than the set electromagnetic similarity threshold, it means that there is a significant difference between the simulation result and the measured data, and the model needs to be adjusted and optimized.

[0125] S67, weightedly fuse the geological stability error, the gravity matching error and the electromagnetic reflection error to obtain a physical constraint loss term; and construct an optimization objective function based on the physical constraint loss term.

[0126] Specifically, the purpose of weighted fusion is to balance the influence of different error terms, so that the optimization objective function can consider the physical constraints of geological stability, gravity field and electromagnetic field. In actual operation, the selection of weight coefficients can be adjusted according to the specific application scene and data characteristics. For example, in geological disaster warning, the weight of geological stability error can be set higher to highlight the constraint of geological stability; in resource exploration, the weights of gravity matching error and electromagnetic reflection error can be set higher to highlight the constraints of gravity field and electromagnetic field.

[0127] The optimization objective function combines the physical constraint loss term and the data-driven loss term. The data-driven loss term is used to evaluate the difference between the prediction result and the observation data, while the physical constraint loss term is used to ensure that the prediction result conforms to the physical law. By minimizing the optimization objective function, the optimal neural network parameters can be obtained, thereby realizing high-precision modeling of underground geological structure.

[0128] The design of the optimization objective function can also consider how to effectively balance the influence of data-driven terms and physical constraint terms. For example, an adaptive weight adjustment mechanism can be introduced to dynamically adjust the weight coefficients according to the error changes during training, so as to improve the convergence speed and generalization ability of the model.

[0129] S68, update the parameters of the neural network architecture through the backpropagation algorithm to reduce the value of the optimization objective function; repeat the forward propagation to parameter update operation until convergence to obtain a trained neural implicit inference model.

[0130] Specifically, the gradient is calculated layer by layer from the output layer to the input layer by the chain rule, and the network parameters are updated using the gradient descent method until the optimization objective function converges. The convergence criterion can be that the change rate of the loss function value is lower than a certain threshold, or the error on the validation set no longer decreases significantly. To prevent overfitting, an early stopping strategy can be used during training, that is, when the error on the validation set no longer decreases for a certain number of iterations, the training is stopped.

[0131] In an alternative embodiment, the implicit terrain representation is converted into an initial three-dimensional mesh, surface reconstruction is performed on the initial three-dimensional mesh, and an explicit three-dimensional terrain model containing surface structures and underground negative obstacles is generated, including the following steps:

[0132] S71, spatial sampling is performed on the implicit terrain representation to obtain a discrete signed distance field.

[0133] Specifically, the implicit terrain representation is a continuous mathematical description that defines the distance of any point in space to the terrain surface through a signed distance function (SDF). To convert it into a discrete form that can be processed by a computer, spatial sampling is required.

[0134] In terms of sampling resolution, higher resolution can capture more details, but will result in a sharp increase in data volume, increasing the burden of subsequent processing. Therefore, the resolution can be reasonably determined according to the application scenario. For example, in city-level terrain modeling, the resolution can be set to the meter level; while in fine archaeological site modeling, the resolution may need to reach the centimeter level. The sampling mode can use regular grid sampling, as it can ensure uniform distribution of sampling points, facilitating subsequent three-dimensional reconstruction processing.

[0135] The construction of the signed distance field (SDF) is the core goal of spatial sampling. During the sampling process, for each sampling point, the distance to the nearest terrain surface is calculated, and a positive or negative sign is assigned to represent the relative position relationship between the point and the surface. There are mainly two methods to calculate the signed distance: one is accurate calculation based on geometric models, and the other is approximate calculation based on numerical methods.

[0136] The accurate calculation based on geometric models is suitable for cases where high-precision geometric models already exist. At this time, geometric algorithms can be used, such as drawing a perpendicular line from the sampling point to the surface of the geometric model, and calculating the length of the perpendicular line as the signed distance. However, in most real-world scenarios, accurate geometric models do not exist, so approximate calculation can be performed with the help of numerical methods. Common numerical methods include ray casting and gradient descent-based optimization methods.

[0137] The ray casting method shoots rays from the sampling points in multiple directions, detects the intersection with the existing point cloud data or image data, and determines the signed distance according to the distance and direction of the nearest intersection point.

[0138] The gradient descent-based optimization method converts the signed distance calculation into an optimization problem. A local search is performed in the neighborhood of the sampling point, and the search direction is adjusted iteratively to gradually approach the terrain surface, and finally the signed distance is obtained.

[0139] To balance the computational efficiency and accuracy, a multi-level sampling strategy can be used. First, coarse sampling with low resolution is performed to quickly obtain a rough signed distance field. Then, according to the needs, fine sampling with high resolution is performed in the area with rich details to supplement more detailed information. In addition, parallel computing technology can be used to distribute the sampling task to multiple processors or computing units to improve the sampling efficiency.

[0140] During the sampling process, the signed distance field can also be smoothed to reduce noise and artifacts in the sampling process. The smoothing methods include Gaussian filtering and median filtering, etc. Gaussian filtering can effectively reduce high-frequency noise by convolution operation. Median filtering has better suppression effect on outliers and is suitable for processing terrain data containing sharp features.

[0141] S72, performing moving cube isosurface extraction on the discrete signed distance field to obtain an initial triangular mesh.

[0142] Specifically, the Marching Cubes algorithm is an isosurface extraction algorithm widely used in medical imaging, geological modeling and computer graphics, etc. Its basic principle is to divide the discrete three-dimensional data into a series of cubic cells, and then determine the intersection of the isosurface and the cubic cell according to the signed distance value of each cubic vertex. For each cubic, according to the relationship between the vertex value and the isosurface threshold, the intersection points of the isosurface and the cubic edge are determined, and the positions of the intersection points are calculated by interpolation. Finally, according to the pre-defined lookup table, the intersection points are connected into triangular patches, thereby constructing the triangular mesh model of the isosurface.

[0143] In practical applications, the discrete signed distance field may have data discontinuity or noise, which will affect the extraction effect of the Marching Cubes algorithm. To improve the robustness of isosurface extraction, some preprocessing and post-processing measures can be taken.

[0144] In the preprocessing stage, data interpolation and filtering techniques can be used to smooth and repair the signed distance field. Data interpolation methods such as bicubic interpolation can increase the continuity of the data, while filtering techniques such as anisotropic diffusion filtering can reduce noise while preserving details.

[0145] In the post-processing stage, it is necessary to clean up and optimize the extracted initial triangular mesh, including removing isolated triangles, repairing the topology of the mesh, and smoothing the mesh surface. For example, by calculating the area and normal vector of each triangle, identify and remove triangles with too small area or abnormal normal vector. At the same time, use mesh smoothing algorithms such as Laplacian smoothing to adjust the mesh vertices and improve the overall quality of the mesh.

[0146] To further improve the efficiency and accuracy of isosurface extraction, adaptive subdivision technology can be introduced. Adaptive subdivision dynamically adjusts the size and subdivision level of the cube according to the local variation of the signed distance field. In areas with sharp changes in signed distance (such as steep slopes or depressions in the terrain), finer subdivision is performed to capture more details; while in areas with gentle changes in signed distance (such as flat areas), coarser subdivision is used to reduce the amount of calculation.

[0147] S73, performing an edge collapse simplification operation on the initial triangular mesh to obtain a simplified mesh model.

[0148] Specifically, edge collapse is a mesh simplification algorithm, the core of which is to reduce the number of triangles by iteratively collapsing edges in the mesh, while preserving the main geometric features and topology of the mesh as much as possible.

[0149] The basic process of edge collapse operation is: evaluate the importance of each edge in the mesh, according to the length of the edge, the geometric variation of the surrounding triangles, and the curvature of the region where the edge is located, etc. The importance of the edge reflects the influence of the edge on the overall shape of the mesh after the edge is collapsed. According to the importance of the edge from small to large, the edge with the lowest importance is preferentially selected for collapse. The collapse operation merges the two vertices of an edge into a new vertex, and updates the triangles adjacent to the edge. After each edge collapse operation, the importance of the surrounding edges affected by the operation is re-evaluated, and the sorted list is updated. Repeat the above process until the desired simplification target is reached, such as reducing the number of triangles to a specified number or achieving the desired simplification rate.

[0150] In the mesh simplification process, the main features and details of the terrain are preserved. To achieve this goal, feature preservation mechanisms and error control strategies can be introduced into the edge collapse algorithm.

[0151] The feature preservation mechanism focuses on identifying and protecting key features of the terrain, such as peaks, valleys, rivers, etc. This can be achieved by giving higher weights to edges in feature areas during edge importance evaluation, making them less likely to be collapsed. In addition, a feature-based multi-resolution representation method can be used to divide the terrain into different feature areas and simplify each area separately to better preserve local features.

[0152] Error control strategies aim to limit the extent of shape deviation from the original mesh due to simplification. Geometric error metrics can be employed to evaluate the shape changes induced by simplification operations. Geometric error metrics include vertex position error, normal vector variation error, and volume variation error, etc. In the edge collapse process, the geometric error between the new vertex position and the original mesh surface is calculated whenever an edge is collapsed, and compared with a pre-set error threshold. If the error exceeds the threshold, the edge is not allowed to collapse, or the position of the new vertex can be adjusted to reduce the error.

[0153] To ensure that the simplified mesh model meets the application requirements, quality evaluation and iterative optimization can be performed on the simplification process. Quality evaluation indicators can be defined from multiple aspects such as geometric accuracy, topological correctness, and visual fidelity. Geometric accuracy evaluation mainly focuses on the positional deviation and shape similarity between the simplified mesh and the original mesh; topological correctness evaluation checks whether the simplified mesh has topological abnormalities, such as non-manifold edges or holes; visual fidelity evaluation focuses on measuring the appearance difference between the simplified mesh and the original mesh from the perspective of human visual perception.

[0154] Based on the results of quality evaluation, the parameters of the edge collapse algorithm can be adjusted, such as modifying the edge importance evaluation function, adjusting the error threshold, etc., and the simplification operation is performed again. Through multiple iterations of optimization, the quality of the simplified mesh is gradually improved, so that it meets the simplification target while being as close as possible to the characteristics and appearance of the original mesh.

[0155] S74, perform texture mapping on the simplified mesh model based on the surface image data to obtain a textured three-dimensional model.

[0156] Specifically, the corresponding texture coordinates are calculated for each vertex of the simplified mesh model. Texture coordinates define the projected position of a vertex on a two-dimensional texture image, which can be represented as a normalized two-dimensional vector (u, v), where u and v correspond to the horizontal and vertical coordinates of the texture image, respectively. There are multiple methods to calculate texture coordinates, including planar projection, cylindrical projection, and spherical projection, etc. For terrain models, planar projection is the most commonly used method, as it can better maintain the planar structure of the terrain and the integrity of the texture.

[0157] According to the texture coordinates of the vertices, the corresponding pixel values in the surface image data are sampled. To improve the quality of texture mapping, methods such as bilinear interpolation or cubic spline interpolation can be used to obtain smoother pixel value transition effects on the texture image. Bilinear interpolation obtains the pixel value at the target position by calculating the weighted average of the four adjacent pixels, which can effectively reduce the mosaic effect in texture mapping; cubic spline interpolation is based on higher-order polynomial functions for interpolation, which can generate smoother texture transitions.

[0158] The pixel values obtained by sampling are assigned to the corresponding vertices of the simplified mesh model, and the texture color values are filled to the entire triangular surface through an interpolation algorithm. During the texture mapping process, attention should be paid to the seamless splicing of the texture and the fusion of multiple texture images. For terrain models spanning multiple images, texture fusion techniques such as weighted average fusion or multi-resolution pyramid fusion can be used to eliminate obvious seams and inconsistencies at the texture boundaries.

[0159] To improve the effect and quality of texture mapping, it is necessary to preprocess and enhance the surface image data. The preprocessing steps include radiation correction, geometric correction, and atmospheric correction of the image, etc. to eliminate various distortions and errors in the imaging process and obtain accurate texture information.

[0160] Radiometric correction is mainly used to correct the radiometric brightness differences in the image caused by sensor response characteristics, sunlight conditions, and terrain shadows, etc. Methods include statistical-based radiometric correction, physical model-based radiometric correction, etc. Geometric correction is to convert the image data from the original sensor coordinate system to a unified geographic coordinate system and correct the geometric distortion in the image, such as perspective distortion, lens distortion, etc. Geometric correction usually uses a polynomial correction model or an affine transformation method based on control points. Atmospheric correction aims to eliminate the effects of atmospheric scattering and absorption on image data, restore the true value of surface reflectance, and improve the accuracy of texture.

[0161] Texture enhancement techniques can further improve the visual effect of texture mapping. For example, by enhancing the contrast, sharpening the edges, adjusting the color balance, etc. of the image, the texture becomes clearer and more realistic. In addition, texture synthesis techniques can also be used to synthesize and repair the texture information in areas where the image data is missing or of poor quality, ensuring the integrity and consistency of the texture mapping.

[0162] S75, performing a Poisson equation solving operation on the textured three-dimensional model to obtain a topologically optimized explicit three-dimensional terrain model.

[0163] Specifically, the basic principle of Poisson equation solving is: assuming that the input data (such as point cloud or mesh model) provides a direction field or normal vector field, find a scalar function f, so that its gradient field is as consistent as possible with the input direction field d. At the same time, in order to ensure the uniqueness of the solution, appropriate boundary conditions need to be imposed. In terrain modeling, the homogeneous Dirichlet boundary condition is usually used, i.e. setting the value of f to zero on the boundary of the model.

[0164] The algorithm steps of Poisson equation solving include data preprocessing, voxelization, sparse matrix construction and solving, and isosurface extraction, etc.

[0165] In the data preprocessing stage, the input textured 3D model is preprocessed, and the normal vector information of each vertex is extracted. If there is no normal vector data in the model, the normal vector can be estimated by fitting or calculating the gradient of the local neighborhood. At the same time, the model is normalized and mapped into a unit cube to facilitate subsequent calculation and solution.

[0166] In the voxelization stage, the normalized model is voxelized, i.e. the continuous three-dimensional space is divided into a regular voxel grid. Each voxel records the point cloud density and normal vector statistical information in the region. The resolution of voxelization directly affects the accuracy and computational complexity of the Poisson equation solution, and can be reasonably selected according to actual needs.

[0167] In the sparse matrix construction and solution stage, according to the voxelized data, the sparse matrix system corresponding to the Poisson equation is constructed. The matrix system contains a large number of unknowns (each voxel corresponds to an unknown scalar value), but has sparsity, so efficient sparse matrix solving algorithms such as conjugate gradient method or multigrid method can be used for fast solution. In the solving process, through iterative optimization, the approximate solution of the scalar function f that satisfies the Poisson equation and the boundary condition is gradually obtained.

[0168] In the isosurface extraction stage, the obtained scalar function f is used to generate the final topologically optimized 3D terrain model through isosurface extraction algorithms such as MarchingCubes algorithm. The zero isosurface of f can be selected as the extraction target because it can well balance the integrity and simplicity of the model.

[0169] The 3D model obtained by solving the Poisson equation is generally reasonable in topological structure, but may still need further topological optimization and post-processing. The purpose of topological optimization is to eliminate possible topological defects in the model, such as non-manifold edges and internal closed holes, to ensure the topological consistency of the model.

[0170] The post-processing steps include smoothing of the model, recalculation of the normal vector, and updating of the texture mapping, etc. Smoothing can eliminate the step effect and noise on the model surface by applying surface smoothing algorithms such as Taubin smoothing or bilateral filtering. Recalculation of the normal vector is based on the smoothed model surface to provide accurate normal vector information for subsequent lighting rendering and visual display. At the same time, according to the changes of the model after topological optimization, the texture mapping can be adjusted and updated to ensure the correct alignment of the texture with the model surface and the consistency of the visual effect.

[0171] The terrain real scene modeling method integrates multi-source heterogeneous data related to the terrain in the target area, generates a multi-modal data cube through noise filtering and data fusion in a unified coordinate system, identifies data hollow areas and extracts geophysical features based on density clustering and elevation gradient analysis, innovatively constructs a neural implicit inference model that fuses physical constraints to jointly model the surface and underground structure, and finally converts the implicit terrain representation into an explicit three-dimensional terrain model. The seamless fusion modeling of the surface and underground space in a complex shielding environment is realized, the physical rationality of the model is significantly improved, the reconstruction accuracy of the hidden area is improved, the generated model strictly conforms to the geological physical law, and the safety and reliability of the geological engineering are effectively ensured.

[0172] It should be understood that, although each step in the flowchart involved in each embodiment as described above is displayed in sequence according to the arrow, these steps are not necessarily executed in the order indicated by the arrow. Unless otherwise specified herein, the execution of these steps is not strictly limited in sequence, and these steps can be executed in other orders. Moreover, at least part of the steps in the flowchart involved in each embodiment as described above can include multiple steps or stages, which are not necessarily executed at the same time, but can be executed at different times, and the execution order of these steps or stages is not necessarily sequential, but can be executed alternately or alternately with at least part of other steps or steps or stages in other steps.

[0173] Based on the same inventive concept, the embodiments of the present application also provide a device for implementing the terrain real scene modeling method described above. The implementation scheme for solving the problem provided by the device is similar to the implementation scheme described in the above method, so the specific limitations in one or more terrain real scene modeling device embodiments provided below can refer to the limitations of the terrain real scene modeling method described above, which will not be repeated here.

[0174] In one exemplary embodiment, as shown in Figure 3 A terrain real scene modeling device 30 is provided for implementing the method in each of the method embodiments described above, and the device includes:

[0175] The data acquisition and integration module 31 is configured to acquire multi-source heterogeneous data related to the terrain in the target area, and the multi-source heterogeneous data includes surface image data, surface point cloud data, and electromagnetic wave reflection data.

[0176] The data preprocessing module 32 is configured to perform noise filtering and ground point classification on the surface point cloud data to obtain a digital surface model, and perform time domain filtering and gain adjustment on the electromagnetic wave reflection data to obtain standardized radar reflection data.

[0177] The coordinate alignment module 33 is configured to perform coordinate system alignment on the ground surface image data, the digital surface model and the standardized radar reflection data based on a spatial coordinate conversion operation to obtain a multi-modal data cube.

[0178] The cavity and obstacle identification module 34 is configured to perform density clustering on the multi-modal data cube to obtain a data cavity region, and perform elevation gradient analysis on the data cavity region to obtain a potential negative obstacle region.

[0179] The geophysical feature extraction module 35 is configured to perform reflection waveform feature extraction on the potential negative obstacle region to obtain a geophysical feature.

[0180] The model construction and training module 36 is configured to construct a neural network model with physical constraints based on the geophysical feature and the multi-modal data cube, and train the neural network model so that a prediction result meets both observation data of the multi-modal data cube and a physical law represented by the geophysical feature, to obtain a neural implicit inference model, wherein the neural network model takes spatial coordinates as input and outputs a signed distance function and a material attribute.

[0181] The three-dimensional modeling module 37 is configured to obtain an implicit terrain representation based on the neural implicit inference model, convert the implicit terrain representation into an initial three-dimensional grid, and perform surface reconstruction on the initial three-dimensional grid to generate an explicit three-dimensional terrain model containing ground surface structures and underground negative obstacles.

[0182] Embodiments of the present application also provide a computer device including a memory and a processor, the memory storing a computer program, and the processor implementing steps in each method embodiment as described above when executing the computer program.

[0183] Embodiments of the present application also provide a computer readable storage medium storing a computer program, the computer program being executed by a processor to implement steps in each method embodiment as described above.

[0184] For the device embodiment, since it basically corresponds to the method embodiment, the related parts are described in the part of the method embodiment. The device embodiments described above are only illustrative, and the components described as separate components can be or can not be physically separated, and the components displayed as units can be or can not be physical units, i.e., they can be located in one place or distributed on multiple network units. Some or all of the modules can be selected to achieve the purpose of the present disclosure according to actual needs. Those skilled in the art can understand and implement it without creative labor.

[0185] The above-described embodiments only express several implementation manners of the application, the description is more specific and detailed, but it cannot be understood as the limitation of the patent scope of the application. It should be pointed out that for ordinary skilled in the art, without departing from the concept of the application, several modifications and improvements can be made, which are within the protection scope of the application.

Claims

1. A terrain real scene modeling method, characterized in that: The method comprises: S1. Collecting multi-source heterogeneous data related to the terrain in the target area, wherein the multi-source heterogeneous data includes surface image data, surface point cloud data, and electromagnetic wave reflection data; S2. performing noise filtering and ground point classification on the surface point cloud data to obtain a digital surface model; performing time domain filtering and gain adjustment on the electromagnetic wave reflection data to obtain standardized radar reflection data; S3. Based on a spatial coordinate conversion operation, aligning the coordinate systems of the surface image data, the digital surface model, and the standardized radar reflection data to obtain a multimodal data cube; S4. Performing density clustering on the multimodal data cube to obtain a data hole area; performing elevation gradient analysis on the data hole area to obtain a potential negative obstacle area; S5. Extracting reflection waveform features from the potential negative obstacle area to obtain geophysical features. S6. Constructing a neural network model that integrates physical constraints based on the geophysical features and the multimodal data cube; training the neural network model so that prediction results are consistent with the observed data of the multimodal data cube and conform to the physical laws represented by the geophysical features, thereby obtaining a neural implicit inference model; wherein the neural network model takes spatial coordinates as input and uses a signed distance function and material properties as output; S7. Obtain an implicit terrain representation based on the neural implicit inference model; convert the implicit terrain representation into an initial three-dimensional mesh, perform surface reconstruction on the initial three-dimensional mesh, and generate an explicit three-dimensional terrain model including surface structures and underground negative obstacles.

2. The method according to claim 1, characterized in that The S4 includes: S41, performing density clustering on the surface point cloud data in the multimodal data cube to generate a set of density anomaly regions; S42, performing elevation gradient calculation on the set of density anomaly areas to generate an elevation gradient distribution map of each area; S43, identifying a concave area by calculating a normalized gradient difference based on the elevation gradient distribution map; S44, performing terrain continuity analysis on the concave area by calculating the curvature change rate of adjacent grids in the concave area, and screening abnormal isolated areas as a candidate negative obstacle area set; S45. Based on the optical texture data in the multimodal data cube, exclude vegetation-covered areas through HSV color space segmentation, perform occlusion type verification on the candidate negative obstacle area set, and obtain the verified potential negative obstacle area.

3. The method according to claim 1, characterized in that The S5 includes: S51, performing a ground penetrating radar B-scan profile extraction operation on the potential negative obstacle area to obtain a reflection waveform dataset; S52, performing adaptive threshold segmentation on the reflection waveform data set to obtain reflection feature parameters; S53, collecting gravity field distribution data of points corresponding to the potential negative obstacle area using a gravity measuring instrument; performing Bouguer anomaly calculation on the gravity field distribution data to obtain an underground density anomaly distribution map; S54: Fusing the reflection characteristic parameters and the underground density anomaly distribution map to obtain the geophysical characteristics.

4. The method according to any one of claims 1 to 3, characterized in that The S6 includes: S61, performing spatial coordinate analysis on the multimodal data cube to obtain a discrete spatial coordinate set; performing parameter deconstruction on the geophysical features to obtain a geological constraint matrix, a gravity difference threshold, and an electromagnetic similarity threshold; the geological constraint matrix represents a mapping relationship of lithologic physical parameters; S62, constructing a neural network architecture including a signed distance prediction branch and a material attribute prediction branch, inputting the discrete space coordinate set into a position encoder for feature enhancement to obtain a high-dimensional feature vector set; inputting the discrete space coordinate set into a position encoder to obtain a high-dimensional feature vector set; S63, inputting the high-dimensional feature vector set into the neural network to perform forward propagation calculation to generate predicted symbol distance values ​​and predicted material attribute values; S64, reconstructing a three-dimensional terrain geometry based on the predicted signed distance value, and calculating a geological stability error using the rock physical parameters in the geological constraint matrix and the geotechnical equilibrium criterion; S65. Reconstructing the underground medium density distribution based on the predicted material property value to obtain a simulated gravity field; when the deviation value between the measured gravity field and the simulated gravity field exceeds the gravity difference threshold, using the deviation value between the measured gravity field and the simulated gravity field as a gravity matching error; S66. Reconstructing a medium electromagnetic parameter distribution based on the predicted material property value to obtain a simulated waveform; when the similarity between the simulated waveform and the measured waveform is lower than the electromagnetic similarity threshold, calculating a waveform difference between the simulated waveform and the measured waveform as an electromagnetic reflection error; S67: Weightedly fuse the geological stability error, the gravity matching error, and the electromagnetic reflection error to obtain a physical constraint loss term; and construct an optimization objective function based on the physical constraint loss term. S68. Update the parameters of the neural network architecture through the back-propagation algorithm to reduce the value of the optimization objective function; repeat the forward propagation to parameter update operation until convergence to obtain the trained neural implicit inference model.

5. The method according to claim 4, characterized in that The step of converting the implicit terrain representation into an initial three-dimensional mesh, performing surface reconstruction on the initial three-dimensional mesh, and generating an explicit three-dimensional terrain model including surface structures and underground negative obstacles comprises: S71. spatially sampling the implicit terrain representation to obtain a discrete signed distance field; S72, performing marching cube isosurface extraction on the discrete signed distance field to obtain an initial triangular mesh; S73, performing an edge collapse simplification operation on the initial triangular mesh to obtain a simplified mesh model; S74, performing texture mapping on the simplified mesh model based on the surface image data to obtain a textured three-dimensional model; S75 , performing a Poisson equation solving operation on the textured three-dimensional model to obtain the topology-optimized explicit three-dimensional terrain model.

6. A terrain real scene modeling device, used to implement the method according to any one of claims 1 to 5, characterized in that: The device comprises: A data acquisition and integration module is used to collect multi-source heterogeneous data related to the terrain in the target area, wherein the multi-source heterogeneous data includes surface image data, surface point cloud data, and electromagnetic wave reflection data; A data preprocessing module is used to perform noise filtering and ground point classification on the surface point cloud data to obtain a digital surface model; and perform time domain filtering and gain adjustment on the electromagnetic wave reflection data to obtain standardized radar reflection data; A coordinate alignment module is used to align the surface image data, the digital surface model and the standardized radar reflection data based on a spatial coordinate conversion operation to obtain a multimodal data cube; A hole and obstacle identification module is used to perform density clustering on the multimodal data cube to obtain data hole areas; perform elevation gradient analysis on the data hole areas to obtain potential negative obstacle areas; A geophysical feature extraction module is used to extract reflection waveform features of the potential negative obstacle area to obtain geophysical features; a model building and training module for constructing a neural network model that incorporates physical constraints based on the geophysical features and the multimodal data cube; training the neural network model so that prediction results are consistent with the observed data of the multimodal data cube and conform to the physical laws representing the geophysical features, thereby obtaining a neural implicit inference model; wherein the neural network model takes spatial coordinates as input and uses signed distance functions and material properties as output; A three-dimensional modeling module is used to obtain an implicit terrain representation based on the neural implicit inference model; convert the implicit terrain representation into an initial three-dimensional mesh, perform surface reconstruction on the initial three-dimensional mesh, and generate an explicit three-dimensional terrain model including surface structures and underground negative obstacles.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the method according to any one of claims 1 to 5 is implemented.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 5 is implemented.

Images