A method for karst cave system structure identification and visual modeling

Abstract

This invention belongs to the field of 3D modeling technology, specifically relating to a method for structural identification and visualization modeling of underground karst cave systems. The method includes the following steps: Step 1: Moving a data acquisition platform along a preset survey line within the underground karst cave, simultaneously acquiring laser scanning point cloud data and generating a point cloud coordinate system; Step 2: Generating cave structure evidence using an extended polarization reflection path weaving algorithm based on orthogonal dipole full polarization physical separation; Step 3: Generating the cave's implicit surface through variational implicit surface solving: establishing a voxel mesh and constructing an implicit distance field within a unified coordinate system; Step 4: Outputting 3D visualization modeling: performing iso-extraction on the cave's implicit surface to generate a 3D mesh model of the cave, and outputting a 3D geological profile and visualization model. This invention achieves stable identification and highly consistent modeling of complex spatial structures in underground karst caves, significantly improving the accuracy, repeatability, and engineering application value of cave structure interpretation.

Description

Technical Field

[0001] This invention belongs to the field of 3D modeling technology, specifically relating to a method for structural identification and visualization modeling of underground karst cave systems. Background Technology

[0002] Underground caves are widely found in soluble rock formations such as limestone and dolomite, often appearing as natural underground galleries, underground rivers, or large cavities. With the increasing demand for underground caves in tourism development, urban underground space utilization, and geological disaster prevention, accurate identification and visual modeling of their spatial structures has become a crucial technical challenge in engineering practice and scientific research. On the one hand, the internal spatial structure of underground caves is complex, with multi-scale, multi-branched, and highly irregular morphologies. On the other hand, the enclosed environment, poor lighting, high humidity, and presence of water in underground caves make conventional surface surveying techniques difficult to apply directly.

[0003] Current underground cave exploration and modeling technologies primarily rely on single sensing methods or weak fusion approaches. In non-contact detection, ground-penetrating radar (GPR) is widely used due to its ability to penetrate underground cavities and media interfaces. A common practice is to use single-polarization GPR to collect echo data along the survey line, identifying underground cavities based on reflected wave amplitude, arrival time, and hyperbolic morphology. However, single-polarization GPR only records echo information in one polarization direction, failing to distinguish echo variations caused by different scattering mechanisms. In underground cave environments, the scattering characteristics of electromagnetic waves differ significantly between the cave's surrounding rock and the water within the cave. Water often produces strong reflections, multiple reflections, and polarization rotation effects. These effects are highly superimposed on the reflections from the cave interface in the single-polarization echo, leading to a high misjudgment rate. Even with multi-line cross-verification methods, it is difficult to reliably distinguish cave boundaries from water multiples in space. Summary of the Invention

[0004] The main objective of this invention is to provide a method for identifying and visually modeling the structure of underground karst cave systems, which can achieve stable identification and highly consistent modeling of the complex spatial structure of underground karst caves, significantly improving the accuracy, repeatability, and engineering application value of cave structure interpretation.

[0005] To address the aforementioned technical problems, this invention provides a method for structural identification and visual modeling of underground karst cave systems, comprising the following steps: Step 1: Move the data acquisition platform along the pre-set survey line inside the underground cavern and use the fully polarized ground-penetrating radar system to collect fully polarized echo data of four polarization combinations: horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal. The fully polarized ground-penetrating radar system adopts an orthogonal dipole antenna design and sets the same trigger source for the four polarization combinations to generate a unified timestamp; simultaneously acquire laser scanning point cloud data and generate a point cloud coordinate system. Step 2: Generate cave structure evidence body based on the extended polarization reflection path weaving algorithm of orthogonal dipole full polarization physical separation. The cave structure evidence body includes cave evidence count, water evidence count, cave connectivity label and lithological interface label. Step 3: Variational implicit surface solution to generate the implicit surface of the cave: A voxel mesh is established in a unified coordinate system and an implicit distance field is constructed. A total driving force field is constructed, consisting of the consistency driving force of the cave structure evidence volume, the consistency driving force of the point cloud, the surface smoothness driving force, and the connectivity preservation driving force. The implicit distance field is updated using an alternating iterative process to output the implicit surface of the cave. Step 4: Output of 3D visualization modeling: Perform iso-value extraction on the implicit curved surface of the cave to generate a 3D mesh model of the cave, map the lithological interface labels and water evidence counts to the 3D geological attribute layer of the 3D mesh model of the cave, and output the 3D geological profile and visualization model.

[0006] Furthermore, by calibrating the target, an external parameter relationship between the coordinate system of the fully polarimetric ground-penetrating radar system and the point cloud coordinate system is established, and a unified coordinate system is generated. Based on the external parameter relationship, the fully polarimetric echo data and the laser scanning point cloud data are mapped to the unified coordinate system. In step 2, the fully polarimetric echo data is input into the polarization separation module. In the same polarization receiving link of the orthogonal dipole antenna, in-phase synthesis is performed to generate a common polarization synthesis channel. In the cross-polarization receiving link, orthogonal phase synthesis is performed to generate a polarization rotation channel. The common polarization coherent descriptor, the cross-polarization dominant descriptor, and the polarization rotation consistent descriptor are calculated on the common polarization synthesis channel and the polarization rotation channel, respectively. Based on the common polarization coherent descriptor, the cross-polarization dominant descriptor, and the polarization rotation consistent descriptor, the reflection event chain is generated and classified and filtered. Voxel evidence voting back projection is performed in combination with the line propagation speed.

[0007] Furthermore, in step two, the four-channel in-situ alignment and baseband construction process of the extended polarization reflection path weaving algorithm includes: performing sampling start point alignment at the same time stamp and generating channel sequences for the four polarization combinations of horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal; performing in-phase component extraction and orthogonal component extraction on each channel sequence to form a complex baseband channel sequence, and arranging the complex baseband channel sequences according to a unified time stamp to form a four-channel complex baseband set; performing in-phase synthesis on the horizontal-horizontal and vertical-vertical complex baseband sets in the in-polarization receiving link as input for the cave main reflection component data; and performing orthogonal phase synthesis on the horizontal-vertical and vertical-horizontal complex baseband sets in the cross-polarization receiving link as input for the water body multiple wave component data.

[0008] Furthermore, in step two, the polarization structure descriptor construction process includes: for each trace sequence of each survey line, a common-polarization coherent descriptor is calculated on the common-polarization synthesis channel using a sliding time window. The common-polarization coherent descriptor is obtained by performing point-by-point conjugate multiplication on the complex baseband sample values ​​of the common-polarization synthesis channel within the sliding time window, accumulating them within the sliding time window, and then performing amplitude normalization; a cross-polarization dominant descriptor is calculated on the polarization rotation channel using a sliding time window. The cross-polarization dominant descriptor is obtained by taking the amplitude of the complex baseband sample values ​​of the polarization rotation channel within the sliding time window and accumulating them within the sliding time window; a polarization rotation consistent descriptor is calculated on the horizontal-vertical channel and the vertical-horizontal channel using a sliding time window. The polarization rotation consistent descriptor is obtained by taking the phase of the complex baseband sample values ​​of the horizontal-vertical channel and the vertical-horizontal channel within the sliding time window and calculating the proportion of the phase difference within a preset stable range.

[0009] Furthermore, in step two, the process of extracting reflection events and generating reflection event chains includes: performing envelope generation on the co-polarization synthesis channel and detecting local peaks in each channel sequence to form a cave candidate event set; performing envelope generation on the polarization rotation channel and detecting local peaks in each channel sequence to form a water candidate event set; for the cave candidate event set, establishing a connection between two cave candidate events in adjacent channel sequences whose time positions are adjacent and whose envelope peak shapes meet a preset similarity threshold, and recording the continuously connected event sequences as cave reflection event chains; for the water candidate event set, generating water multiple event chains using the same connection rules as the cave candidate event set, and detecting peak clusters arranged in a repeating delay according to time position in each water multiple event chain to form water multiple cluster labels.

[0010] Furthermore, in step two, the event chain classification process based on descriptor and chain consistency includes: extracting continuous high-value segments of co-polarized coherent descriptors within each cave reflection event chain and generating cave interface consistency labels; extracting continuous high-value segments of cross-polarization dominant descriptors and polarization rotation consistent descriptors within each water multiple event chain, and combining them with water multiple cluster labels to generate water multiple consistency labels; selecting a cave interface candidate set from the cave reflection event chain based on the cave interface consistency labels, and selecting a water multiple candidate set from the water multiple event chain based on the water multiple consistency labels.

[0011] Furthermore, in step two, the process of determining the propagation speed includes: selecting event segments with hyperbolic trajectories from the candidate set of the cave interface, generating theoretical hyperbolas one by one using the candidate propagation speeds, calculating the cumulative time position deviation between the event segments and the theoretical hyperbolas, and selecting the candidate propagation speed with the smallest cumulative deviation as the propagation speed of the survey line.

[0012] Furthermore, in step two, the voxel evidence voting backprojection process includes: establishing a voxel grid in a unified coordinate system and defining cave evidence counts and water evidence counts for the voxel grid; for each cave candidate event in the cave interface candidate set, converting the propagation distance of the cave candidate event into a spherical neighborhood based on the spatial location of the measuring point and the propagation speed of the measuring line, and performing cave evidence count accumulation within the voxels covered by the spherical neighborhood; for each water candidate event in the water multiple wave candidate set, performing water evidence count accumulation within the voxels covered by the spherical neighborhood using the same spherical neighborhood backprojection method as for the cave candidate events; the cave connectivity label and lithological interface label generation process includes: after completing the backprojection of all measuring lines, generating cave connectivity labels by performing three-dimensional connectivity domain marking based on the cave evidence count, and generating lithological interface labels in areas where the spatial gradient of the cave evidence count is greater than a preset gradient threshold.

[0013] Furthermore, in step three, the construction process of the total driving force field includes: generating the interface advancement direction based on the direction from high-count voxels to low-count voxels in the cave evidence count and writing it into the cave structure evidence volume consistency driving force; generating a point cloud distance field on the voxel grid based on laser scanning point cloud data and writing the equidistant surface direction of the point cloud distance field into the point cloud consistency driving force; generating a surface smoothing driving force based on the local second-order difference of the implicit distance field; performing symbolic reinitialization within the connected domain of the implicit distance field based on the cave connectivity label, and performing connection evolution at the contact boundary of the connected domain to write the connectivity maintenance driving force.

[0014] Furthermore, in step four, the three-dimensional visualization modeling output process also includes: automatically cutting the three-dimensional mesh model of the cave body according to preset cutting rules in a unified coordinate system to generate a three-dimensional geological profile, and outputting a visualization result containing the profile outline, profile attributes and water distribution.

[0015] The method for structural identification and visual modeling of underground karst cave systems proposed in this invention has the following beneficial effects: This invention introduces a fully polarized ground-penetrating radar acquisition method with orthogonal dipole antenna design. At the physical level, it simultaneously acquires fully polarized echo data of four polarization combinations: horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal. This allows the differences in electromagnetic scattering characteristics between the surrounding rock and water in underground caves to be fully preserved and utilized in subsequent processing.

[0016] Compared to existing methods that rely solely on single polarization or simple amplitude features, this invention enhances the information dimension at the data source, enabling the explicit differentiation of polarization rotation and multiple wave effects generated by water bodies through polarization physical separation. This significantly reduces misjudgments caused by water interference in cave structure identification. Through an extended polarization reflection path weaving algorithm, multi-polarization echo data is transformed into a cave structure evidence body comprising cave evidence counts, water body evidence counts, cave connectivity labels, and lithological interface labels. This elevates the representation of underground karst cave structures from a single interface interpretation to a multi-evidence, multi-attribute three-dimensional spatial representation, enhancing the stability and repeatability of the results.

[0017] During the modeling phase, this invention employs a variational implicit surface solution method, co-constraining the cave structure evidence volume and laser scanning point cloud data within a unified coordinate system. This ensures that the implicit surface of the cave satisfies geometric smoothness while strictly adhering to the cave evidence distribution and connectivity relationships, effectively avoiding spurious filling and topological breakage problems commonly found in traditional point cloud interpolation modeling. The alternating iterative update mechanism of the implicit distance field enables the model to gradually converge to a stable form in complex cave segments, demonstrating excellent adaptability to cave bifurcation, contraction, and arch changes.

[0018] In the visualization output stage, lithological interface labels and water body evidence counts are mapped to a three-dimensional geological attribute layer of the cave's three-dimensional mesh model. This allows the cave's geometry and geological attribute information to be presented intuitively within the same model, improving the model's interpretability and providing more intuitive and reliable technical support for tourism safety assessments, cave development planning, and public education. Overall, this invention achieves high-precision identification and highly consistent modeling of underground cave systems, possessing significant engineering application value and promotional significance. Attached Figure Description

[0019] Figure 1 A schematic diagram illustrating the automatic generation of a three-dimensional geological profile from a three-dimensional mesh model of a cave provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the variational implicit surface evolution principle based on evidence consistency driving force provided in an embodiment of the present invention; Figure 3 A schematic diagram illustrating the determination of medium propagation velocity based on hyperbolic trajectory fitting, provided in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the visualization effect of an interactive three-dimensional geological attribute model integrating multiple views and contour lines, provided in an embodiment of the present invention. Detailed Implementation

[0020] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0021] A method for structural identification and visualization modeling of underground karst cave systems, comprising the following steps: Step 1: Move the data acquisition platform along the pre-set survey line inside the underground cavern and use the fully polarized ground-penetrating radar system to collect fully polarized echo data of four polarization combinations: horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal. The fully polarized ground-penetrating radar system adopts an orthogonal dipole antenna design and sets the same trigger source for the four polarization combinations to generate a unified timestamp; simultaneously acquire laser scanning point cloud data and generate a point cloud coordinate system. Step 2: Generate cave structure evidence body based on the extended polarization reflection path weaving algorithm of orthogonal dipole full polarization physical separation. The cave structure evidence body includes cave evidence count, water evidence count, cave connectivity label and lithological interface label. Step 3: Variational implicit surface solution to generate the implicit surface of the cave: A voxel mesh is established in a unified coordinate system and an implicit distance field is constructed. A total driving force field is constructed, consisting of the consistency driving force of the cave structure evidence volume, the consistency driving force of the point cloud, the surface smoothness driving force, and the connectivity preservation driving force. The implicit distance field is updated using an alternating iterative process to output the implicit surface of the cave. Step 4: Output of 3D visualization modeling: Perform iso-value extraction on the implicit curved surface of the cave to generate a 3D mesh model of the cave, map the lithological interface labels and water evidence counts to the 3D geological attribute layer of the 3D mesh model of the cave, and output the 3D geological profile and visualization model.

[0022] In one implementation, the pre-planned survey lines are designed before entering the underground cave and implemented in a way that can be reproduced on-site. These pre-planned survey lines are typically laid out along visitor paths or the main corridors of the cave, ensuring that fully polarized echo data and laser scanning point cloud data cover the main geometric interfaces of the cave walls, ceiling, and floor. To enable repeated traversal of the survey lines within the cave, the pre-planned survey lines can be defined by ground markers and survey line direction markers. For example, a survey line marker can be set every 5 meters, with an additional direction marker at bends. This ensures that the data acquisition platform's path within the cave remains consistent with the survey line direction, reducing spatial errors caused by path differences between different data acquisitions.

[0023] Maintaining a constant or segmented constant speed while moving along the preset survey line facilitates the association of fully polarized echo data and laser scanning point cloud data from the same spatial location to similar unified timestamps. The moving speed of the acquisition platform can be selected from 0.2 m / s to 1.0 m / s to balance tunnel access safety and data coverage density. When the tunnel surface is uneven or has steps, the acquisition platform can use a short-distance stop-and-go approach to obtain denser fully polarized echo data and laser scanning point cloud data at the same location, thereby enhancing the ability to characterize the tunnel ceiling cavity boundaries and local concave and convex structures.

[0024] The fully polarimetric ground-penetrating radar system employs an orthogonal dipole antenna design to acquire fully polarimetric echo data for four polarization combinations: horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal. An orthogonal dipole antenna can be understood as two mutually orthogonal dipole radiating elements placed on the same mounting base. One dipole radiating element corresponds to the horizontal polarization direction, and the other corresponds to the vertical polarization direction. Selecting both the transmitter and receiver in the horizontal polarization direction forms a horizontal-horizontal polarization combination; selecting both the transmitter and receiver in the vertical polarization direction forms a vertical-vertical polarization combination; selecting both the transmitter and receiver in the horizontal polarization direction forms a horizontal-vertical polarization combination; and selecting both the transmitter and receiver in the vertical polarization direction forms a vertical-horizontal polarization combination. Through this physical arrangement of the orthogonal dipole antennas, the four polarization combinations can obtain comparable polarization responses at the same spatial measurement point, thus providing fundamental data for subsequent polarization separation and structure identification.

[0025] To ensure comparability among the four polarization combinations, the fully polarimetric ground-penetrating radar system uses the same trigger source to generate a unified timestamp for all four combinations. The same trigger source means that the transmission and reception sampling start times for the four polarization combinations (horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal) are all initiated by the same trigger signal. This trigger signal simultaneously drives the sampling gate of both the transmission pulse generation circuit and the reception sampling circuit, ensuring that the time origin of each echo record is consistent. The direct benefit of this is that the alignment of the echoes of the four polarization combinations on the time axis does not rely on post-hoc calculations to guess the alignment relationship; instead, it is guaranteed by hardware triggering. When strong reflective interfaces exist within the cavern, the arrival time difference between the echoes of the four polarization combinations mainly originates from the differences in polarization scattering by the medium, rather than from trigger drift or sampling start-point drift.

[0026] A unified timestamp can be generated using an incremental counting method. For example, each sampling triggered by the same trigger source generates one trigger count. The trigger count, together with a highly stable clock, generates a unified timestamp, ensuring that each echo record has a traceable acquisition sequence and time location. In subsequent processing, horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal fully polarized echo data from the same location can be grouped into a single set of fully polarized echo data based on the unified timestamp, avoiding mismatches caused by differences in acquisition sequence between different polarization combinations.

[0027] In the specific implementation of acquiring fully polarized echo data, either sequential polarization switching or parallel polarization reception can be used. In the sequential polarization switching method, the same trigger source controls the polarization selection of the transmitter and receiver according to a predetermined switching sequence. For example, a cyclic sequence of horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal can be used to acquire data, allowing the same spatial location to obtain fully polarized echo data of four polarization combinations in a very short time. When the acquisition platform moves at a low speed, this method allows the spatial locations corresponding to the four polarization combinations to be closer together. In the parallel polarization reception method, the transmitter is fixed in either the horizontal or vertical polarization direction, while the receiver simultaneously samples both directions. This allows for the acquisition of echo records with the same polarization and cross-polarization under a single trigger, and then the acquisition of the four polarization combinations is completed through switching at the transmitter. This method reduces the number of polarization switches, which is beneficial for improving acquisition efficiency and reducing the transient effects caused by mechanical or electrical switching.

[0028] Fully polarized ground-penetrating radar (GPR) data is typically stored as channel sequences, with each channel sequence corresponding to a measurement point location or sampling results within a short time period. To enhance the signal-to-noise ratio in karst cave environments, GPR systems can coherently or incoherently superimpose echoes from multiple transmissions under the same trigger count before storing them as a single channel sequence. For example, when electromagnetic interference is strong or medium loss is high within the cave, 16 superpositions can be selected to improve the visibility of stable reflections; when the cave structure changes rapidly or the acquisition platform moves quickly, 4 superpositions can be selected to balance resolution and efficiency.

[0029] There is a correlation between the time samples in fully polarized echo data and the propagation distance in the medium. If the round-trip time of the echo propagation is denoted as... The speed of electromagnetic wave propagation in the medium is denoted as . Then the one-way propagation distance can be calculated as follows: Perform the conversion. Among them, This represents the distance along the propagation path from the antenna phase center to the reflection point. This indicates the speed at which electromagnetic waves propagate through the cave medium. This indicates the time from the triggering moment of the same trigger source to the arrival of the echo at the receiver. The significance of using the same trigger source is also more intuitive here: the four polarization combinations are in the same... Recording echoes in a reference frame allows the same reflecting interface to reflect different polarization combinations. The anomaly can be interpreted as a difference in polarization scattering path or multiple wave path, rather than a pseudo-difference introduced by a change in the triggering origin.

[0030] refer to Figure 3This figure illustrates how the electromagnetic wave propagation velocity in the subsurface medium can be inverted by analyzing the diffraction patterns generated by point scatterers on a B-scan image (range-time profile) of fully polarimetric ground-penetrating radar data. The horizontal axis represents the distance traveled along the pre-defined survey line (in meters), and the vertical axis represents the two-way travel time of the radar echo (in nanoseconds). The alternating red and blue textures in the background represent the actual collected and pre-processed radar echo data. The obvious hyperbolic morphology is generated by local inhomogeneities in the subsurface (such as rock protrusions, cavity edges, or isolated rock blocks) acting as point diffraction sources. Because the distance between the radar antenna and the scattering point changes hyperbolically during the radar's movement, the echo arrival time also exhibits a hyperbolic trajectory. Figure 3 The core of this demonstration lies in showcasing the matching process of velocity analysis. The yellow asterisks marked as extracted event vertices in the diagram represent the peak positions of automatically detected scattering events.

[0031] To determine the medium velocity at this location, the system generated theoretical hyperbolic trajectories corresponding to different candidate velocities. The figure shows three curves of different colors: the green dashed line corresponds to the slower candidate velocity (v1), with a steep curve and a narrow hyperbola opening, significantly deviating from the energy ridge of the actual echo, indicating that the calculated dielectric constant is too high; the cyan dashed line corresponds to the faster candidate velocity (v3), with a gentler curve and an excessively wide opening, also failing to fit the data characteristics well; and the yellow solid line corresponds to the optimal matching velocity (v2), which perfectly coincides with the strong energy stripes in the background data. In the actual processing steps, the algorithm sets a velocity scanning range (e.g., 0.06 m / ns to 0.15 m / ns), generates a series of theoretical trajectories with a certain step size, and calculates the sum of the signal energy or correlation coefficient on each trajectory (i.e., the cumulative deviation or similarity shown in the figure). The velocity with the smallest cumulative deviation or the largest energy superposition value is determined as the representative propagation velocity of that measurement segment. This step is crucial for the subsequent time-depth conversion. Only by obtaining the accurate propagation speed can the time information of the radar echo be accurately converted into spatial depth information, thereby ensuring that the geometric position of the subsequently generated spherical shell evidence body in three-dimensional space is accurate.

[0032] Point cloud generation typically converts polar or spherical coordinate measurements into Cartesian coordinates. For example, if the distance measured by laser is denoted as... The horizontal scan angle is denoted as The vertical scanning angle is denoted as Then the coordinates of a point in the point cloud coordinate system can be expressed as: , , .in, This represents the coordinate components of a point along the forward axis in the point cloud coordinate system. This represents the coordinate components of a point along the lateral axis in the point cloud coordinate system. This represents the coordinate components of a point along the vertical axis in the point cloud coordinate system. Indicates the laser ranging value, and indicates the horizontal scanning angle. This indicates the vertical scanning angle. The advantage of this conversion method is that each laser echo can be directly mapped to a spatial point in the point cloud coordinate system, forming a geometric sample of the tunnel wall, tunnel top, and tunnel bottom. The point cloud coordinate system is naturally bound to the observation geometry of the laser scanning device, which facilitates the subsequent correlation of geometric information with the fully polarized echo data.

[0033] The acquisition frequency of laser scanning point cloud data can be selected according to the spatial scale inside the cave and the speed of the acquisition platform. For example, the laser scanning device can output one frame of point cloud at 10Hz, with each frame containing 200,000 points; when the acquisition platform moves at 0.5m / s, the interval between adjacent point cloud frames in the direction of travel is approximately 0.05m, ensuring sufficient overlap on the continuous surface of the cave and reducing geometric voids. If there is heavy dust or water mist inside the cave, an alternative implementation is to reduce the number of points per frame and increase the frame rate, for example, outputting point clouds at 20Hz to reduce the impact of motion distortion. At the same time, points with echo intensity below a preset intensity threshold are discarded during point cloud generation, thereby suppressing discrete noise caused by suspended particles.

[0034] In an optional implementation, to further reduce the impact of the acquisition platform's movement on data consistency, the acquisition platform can adopt a constant ground-elevation antenna installation method, keeping the distance between the orthogonal dipole antenna and the ground within the range of 0.05m to 0.30m. Maintaining the ground-elevation height is significant because changes in the propagation time of the air layer are directly reflected in the time samples of the fully polarized echo data. A stable ground-elevation height reduces this systematic time drift, making the echo times of the same reflecting interface more consistent across different measurement points, which is beneficial for the generation of subsequent reflection event chains. For laser scanning point cloud data, a stable ground-elevation height also reduces the pitch variation of the laser scanning device, lowers the geometric distortion of the point cloud coordinate system caused by attitude changes, and improves the consistency of point cloud stitching on continuous surfaces within the cave.

[0035] In this way, when the acquisition platform moves along the preset survey line in the underground cave, it can obtain fully polarized echo data with four polarization combinations (horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal) with a unified timestamp, and simultaneously obtain laser scanning point cloud data and point cloud coordinate system. This establishes a directly implementable data foundation for the subsequent extended polarization reflection path weaving algorithm and variational implicit surface solution based on orthogonal dipole fully polarized physical separation.

[0036] In one implementation, before the fully polarized echo data of four polarization combinations (horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal) enter the extended polarization reflection path weaving algorithm, sampling start point alignment under the same timestamp is performed. Alignment uses a unified timestamp as an index, grouping the fully polarized echo data of the four polarization combinations corresponding to the same unified timestamp into a group, and moving the sampling start points of the four channel sequences within each group to the same reference position. This concentrates the manifestation of subsequent polarization differences on medium scattering differences, avoiding the sampling start point offset from splitting the same reflection event into different time positions, and reducing the probability of false connections during reflection event chain generation. For example, sampling start point alignment can use the minimum energy criterion of the pre-triggered segment: finding the interval with the lowest continuous energy at the beginning of each channel sequence as a zero reference, and then aligning the four channel sequences to the same zero reference position.

[0037] To facilitate in-phase and quadrature phase synthesis, the four polarization combinations of fully polarized echo data are converted into complex-base bandpass sequences. In-phase and quadrature component extraction for each channel sequence can be achieved through quadrature mixing: the channel sequence is multiplied by the local oscillator cosine and sine sequences respectively, then low-pass filtered to obtain in-phase and quadrature sample sequences, which are then combined to form a complex-base bandpass sequence. The complex-base bandpass sequence makes the phase information explicit, allowing subsequent polarization rotation consistency descriptors to be directly obtained from the phase difference, reducing errors from indirectly inferring rotation characteristics using peak positions. For example, the sampling interval of the complex-base bandpass sequence can be 0.5 ns, with each channel sequence containing 2048 samples, thus covering the echo time range from near-field to deep regions.

[0038] The orthogonal dipole full polarization physical separation forms a common polarization combining channel and a polarization rotation channel on the same polarization receiving link and cross-polarization receiving link, respectively. When performing in-phase combining of the horizontal and vertical complex baseband sequences within the same polarization receiving link, the phase zero offset between the two channels is first eliminated, and then vector addition is performed. A directly implementable approach is to calculate the cross-correlation phase of the horizontal and vertical complex baseband sequence within the same time window, and use the cross-correlation phase to rotate the phase of one of the complex baseband sequences, making them as in-phase as possible in the main reflection section, and then adding them point by point to form a common polarization combining channel. Phase zero offset correction can suppress phase drift caused by differences in hardware link length or antenna feed line differences, making the in-phase combining closer to the superposition of the same scattering mechanism, thereby highlighting the weak polarization rotation and strong specular reflection components of the cavity interface.

[0039] To more effectively represent the polarization rotation and multiple scattering caused by water, the cross-polarized receiver link performs orthogonal phase synthesis on the complex baseband sets of the horizontal-vertical and vertical-horizontal channels to form a polarization rotation channel. One implementation involves first performing amplitude and phase equalization on the horizontal-vertical and vertical-horizontal channels to make their amplitudes statistically similar in the non-rotational scattering range, and then applying a polarization rotation filter to one of the channels. The phase shifted result is added to another channel to obtain a composite result that is more sensitive to polarization rotation. This composite result is used to characterize the water multiple component data because the strong polarization response of water bodies is more likely to introduce cross-polarization components and polarization direction changes. Orthogonal phase synthesis separates these changes from the same polarization background, making it easier to distinguish the cave main reflection component data from the water multiple component data in subsequent reflection event chain classification.

[0040] To transform the structural information in the co-polarization synthesis channel and the polarization rotation channel into a filterable basis, co-polarization coherent descriptors, cross-polarization dominant descriptors, and polarization rotation consistent descriptors are constructed on both channels. The descriptor calculation employs a sliding time window, with a sliding step size of one or two samples, allowing the descriptor to change continuously over time. This facilitates the location of stable segments of reflection events on each channel sequence. For example, the sliding time window length can be 64 samples, equivalent to a 32 ns time coverage. This length, in the context of a cave, covers both the main lobe and side lobes of a typical reflection wave packet without being too wide to mix wave packets from adjacent interfaces.

[0041] The copolarized coherent descriptor is calculated on the copolarized synthesis channel, aiming to measure the phase continuity and stability of the echo waveform over time. A feasible calculation method involves performing point-by-point conjugate multiplication and summation of adjacent samples within a sliding time window, followed by amplitude normalization. The corresponding calculation can be written as follows: .in, Indicates the first The co-polarization coherent descriptor corresponding to the sliding time window starting from each sample point; The copolarized synthesis channel indicates that the first... The complex basis band sample value at each sample point; express The complex conjugate; Indicates the sample length of the sliding time window; This represents the complex amplitude. By multiplying and summing the conjugates of adjacent samples, the components whose phase remains stable within a reflected wave packet can be accumulated into a large complex vector. This type of specular reflection at the cave interface typically makes... The wave packet segment exhibits continuously high values; multiple waves and scattering are more likely to cause phase perturbations, which cancel each other out after accumulation. The continuous high-value intervals are shorter.

[0042] The cross-polarization dominance descriptor is calculated over a polarization rotation channel, aiming to measure the concentration of cross-polarization energy within a time window. A direct implementation involves taking the amplitude of the complex basis band samples within the polarization rotation channel during the sliding time window and summing them, denoted as... .in, Indicates the cross-polarization dominant descriptor; Indicates the polarization rotating channel in the first... The complex basis band values ​​at each sample point; the meanings of the remaining symbols are consistent with the above. Water bodies and their multiple reflections more easily create sustained energy uplifts in the polarization rotation channel, leading to... The values ​​are significantly higher during the corresponding time period; the interface of the cavern is usually weaker in the polarization rotation channel. It is not easy to rise continuously, therefore This provides a stable basis for screening multiple wave event chains in water bodies.

[0043] The polarization rotation consistent descriptor is calculated on the horizontal-vertical channel and the vertical-horizontal channel, with the goal of determining whether the cross-polarization response exhibits a consistent rotational characteristic. First, the phase of the complex basis band samples in the horizontal-vertical channel and the vertical-horizontal channel is taken separately to form a phase sequence. Then, the phase difference is calculated, and the proportion of phase differences within a preset stable range is statistically analyzed. An implementable expression is: calculation within a sliding time window... , then calculate .in, Indicates the horizontal and vertical channels at the 1st Phase at each sample point; Indicates the vertical and horizontal channels at the 1st Phase at each sample point; Indicates phase difference; Indicates a polarization rotation consistent descriptor; This indicates an indicator function that returns 1 if the condition within the parentheses is true and 0 if it is false. This represents the center value of the phase difference within the sliding time window; the median of the phase difference can be taken. This represents half the width of the preset stable interval. For example, A value of 0.35 rad can be used to ensure the stable region covers common rotationally consistent scattering while still excluding random phase. Water multiples often originate from relatively fixed combinations of reflection paths, and the rotational characteristics of cross-polarization are more consistent across adjacent samples. This creates continuous high-value segments; scattering noise and unstable multipath propagation make phase difference jumps more likely. The value will decrease significantly.

[0044] After descriptor construction is complete, reflection event extraction and reflection event chain generation are performed. Reflection event extraction first generates envelopes on the common-polarization synthesis channel and the polarization rotation channel, then detects local peaks within each channel sequence. The envelope can be obtained using the Hilbert transform: for real signals... Calculate the Hilbert transform envelope Defined as .in, Indicates the channel signal at the 1st Real values ​​at each sample point; express The Hilbert transform result; The envelope represents the energy of the oscillating waveform, which is concentrated into a single peak shape. This makes the detection of local peaks more stable, especially when there is phase reversal or waveform distortion within the hole. Directly finding the peak of the original waveform can easily lead to confusion between the main lobe and the side lobes. The envelope peak is closer to the arrival time of the reflection event.

[0045] The detection of local peaks involves setting both a noise threshold and a peak shape constraint. For example, the noise threshold can be set as the mean of the envelope of the first 256 samples in the channel sequence plus six times its standard deviation. The peak shape constraint requires that the peak be a local maximum within 10 samples before and after it. The advantage of setting the threshold in this way is that system noise and environmental noise in the near-field of the cave are usually concentrated in the initial segment. Establishing the threshold using statistics from the initial segment can adapt to the noise level of the current cave segment, while suppressing the influence of individual spike noise and avoiding treating occasional interference as part of the cave's candidate event set. Local peaks detected in the co-polarization synthesis channel form the cave's candidate event set, while local peaks detected in the polarization rotation channel form the water body candidate event set.

[0046] The reflection event chain generation connects candidate events of adjacent trace sequences along the survey line. The connection rules use two conditions: temporal proximity and envelope peak shape satisfying a preset similarity threshold. Temporal proximity can be defined as the peak time position difference between adjacent trace sequences not exceeding 12 sample points, equivalent to a 6ns arrival time difference. This covers the geometric time difference caused by changes in survey point location and limits excessively wide connections that could lead to cross-interface misconnections. Envelope peak shape similarity is determined using the normalized peak inner product: the peak neighborhood centered on the peak of each candidate event is truncated into a peak sequence of length 21 sample points. The peak sequence is first divided by its own maximum value to normalize the amplitude, and then the inner product of two normalized peak sequences is calculated. .in, Indicates similarity; and These represent the normalized peaked sequences of the two candidate events at the 1st... The values ​​at each position; This represents the length of the peak sequence, with an example of 21. The similarity threshold can be set to 0.85, which makes the connection tend to link reflection events with similar wave packet shapes into a chain. The wave packet shapes at adjacent measuring points of the cave interface are usually more consistent, making it easier to form a stable cave reflection event chain. Water multiples are affected by multipath, and the wave packet shape is more likely to split or tail. The introduction of morphological similarity can prevent water multiples from being mistakenly linked to the cave reflection event chain.

[0047] After generating multiple wave event chains, it is necessary to detect peak clusters arranged with repetitive delays in time position within each multiple wave event chain to form multiple wave cluster tags. A direct approach is to arrange the peak time positions from earliest to latest on each channel sequence for the same multiple wave event chain, calculate the time position difference sequence between adjacent peaks, and mark the set of peaks whose time position differences fall within the same tolerance range for three or more consecutive times as peak clusters. The tolerance range can be four sample points, equivalent to a 2ns time tolerance. The repetitive delay arrangement corresponds to the superposition of multiple round-trip times of reflections along the same path. This regularity is more likely to occur when water interface reflections are strong and the number of reflections is high. The multiple wave cluster tags provide an executable basis for the subsequent generation of consistent multiple wave tags.

[0048] Subsequently, event chain classification based on descriptor and chain consistency is performed to generate cavern interface consistency tags and water multiple wave consistency tags. For each cavern reflection event chain, the copolarization coherent descriptor is read along the corresponding time position of the event chain, and continuous high-value segments are selected as cavern interface consistency tags. Continuous high-value segments can be determined by a copolarization coherent descriptor greater than or equal to 0.70 and the occurrence of at least 8 consecutive channel sequences. The phase continuity of the cavern interface at adjacent measurement points is stronger, and segments satisfying this condition are more concentrated; using a continuous length constraint can filter out short coherent burst noise. For each water multiple wave event chain, the cross-polarization dominant descriptor and polarization rotation consistency descriptor are read along the corresponding time position of the event chain, and continuous high-value segments are selected respectively. These are then combined with water multiple wave cluster tags to generate water multiple wave consistency tags. For example, the cross-polarization-dominant descriptor can be required to be greater than or equal to twice the median of its channel sequence, and the polarization-rotation-consistent descriptor can be required to be greater than or equal to 0.60, with at least 8 consecutive channel sequences, and also include a water multiple cluster label. The advantage of this combination is that the cross-polarization-dominant descriptor is responsible for picking out events with significant energy in the polarization rotation channel, the polarization-rotation-consistent descriptor is responsible for picking out rotational scattering with stable phase differences, and the water multiple cluster label adds the repetition delay characteristic of multiples to the screening condition. The combined constraint of these three factors can significantly reduce the probability of misclassifying cavern boundaries as water multiples.

[0049] A candidate set of cave interfaces is obtained by filtering from the cave reflection event chain based on the cave interface consistency label, and a candidate set of water multiples is obtained by filtering from the water multiples event chain based on the water multiples consistency label. The cave interface candidate set is used for subsequent propagation velocity determination and cave evidence count accumulation, while the water multiples candidate set is used for water evidence count accumulation. The two form complementary evidence information within the voxel grid.

[0050] The propagation velocity is determined by selecting event segments exhibiting hyperbolic trajectories from the candidate set of cavern interfaces. A theoretical hyperbola is generated for each candidate propagation velocity, and the cumulative time position deviation between the event segment and the theoretical hyperbola is calculated. The candidate propagation velocity with the smallest cumulative deviation is taken as the propagation velocity of the survey line. The reason for selecting hyperbolic trajectory event segments is that point scatterers or locally protruding interfaces often exhibit a hyperbolic shape along the survey line. The curvature of the hyperbola is directly affected by the propagation velocity. Hyperbolic matching allows the propagation velocity to be estimated from its geometric shape without relying on prior external dielectric parameters. The theoretical hyperbola can be... Expression, among which, Indicates a horizontal offset of Round-trip transmission time; Indicates the candidate propagation speed; This represents the vertical distance from the scatterer to the antenna phase center; This represents the horizontal offset distance of the measuring point relative to the vertex of the hyperbola. For each The vertex time position of the event segment can be used as the initial value for fitting. Then, the theoretical time position of each measuring point is calculated and the difference is made with the actual event time position. The absolute differences are accumulated to obtain the cumulative time position deviation. For example, the candidate propagation velocity can be taken from 0.06 m / ns to 0.15 m / ns, the step is 0.005 m / ns, and the event segment length can be taken as 25 channel sequences, which can find a stable minimum deviation point within the range of common lithology in the tunnel.

[0051] After determining the propagation velocity of the survey line, a voxel evidence voting back projection is performed to generate cave evidence counts and water evidence counts, and further generate cave connectivity labels and lithological interface labels. The voxel grid is established in a unified coordinate system, with voxel side lengths set to 0.20m to match the voxel scale with the geometric scale of the cave corridor. This effectively represents the curvature changes of the cave boundary without resulting in an excessively large number of voxels leading to excessively long back projection times. For each cave candidate event, the propagation distance is first converted from the event's time and location, then converted to a spherical neighborhood, and cave evidence counts are accumulated within the voxels covered by the spherical neighborhood. The propagation distance conversion uses... ,in, Indicates the one-way propagation distance; Indicates the propagation speed of the measuring line; This represents the round-trip propagation time of the candidate event within the cavern. The round-trip propagation time is calculated by converting the time position of the candidate event within the cavern and the sampling interval. For example, if the sampling interval is 0.5 ns and the time position of the candidate event within the cavern is the 800th sample point, Equal to 400 ns. The realization of the spherical shell neighborhood can be achieved using voxel center distance determination: taking the spatial position of the measurement point corresponding to the candidate event in the cavity as the sphere center, and... For each voxel in the voxel grid that falls within the thickness range of the spherical shell (where the shell radius is used), the cavity evidence count is incremented by 1. The shell thickness can be taken as the diagonal length of the voxel, ensuring continuous coverage without discontinuities during voxel discretization. For each water body candidate event, the same shell neighborhood back projection method as for cavity candidate events is used to accumulate water body evidence counts within the voxels covered by the shell neighborhood, thus forming two types of evidence distributions—cavities and water bodies—on the same voxel grid.

[0052] To ensure that back projection can still be performed directly on the long-distance measurement line of the cave, the search for the overburden voxels in the vicinity of the spherical shell can be limited to a local bounding box centered on the center of the sphere. For example, the side length of the bounding box is taken as... m, which makes the bounding box completely cover the spherical shell and leaves a discrete margin, and then only makes distance determination within the range of the voxel index corresponding to the bounding box, can reduce the number of voxel traversals in a single event back projection from the number of global voxels to the number of local voxels, and the back projection time is more stable.

[0053] After completing the back projection of all survey lines, 3D connectivity labeling is performed based on the cave evidence count to generate cave connectivity labels. A direct approach is to select voxels with a cave evidence count greater than or equal to 3 as the cave candidate voxel set, and then traverse the candidate voxel set using the 26-neighborhood connectivity rule, assigning connected voxel sets the same connectivity domain number. This connectivity domain number is the cave connectivity label. The advantage of setting the cave evidence count threshold to 3 is that if a voxel is repeatedly covered by the spherical neighbor of multiple cave candidate events, it indicates that the voxel is in a stable cave interface back projection intersection region; single or double coverage is more likely to originate from sporadic peaks or local noise, and the threshold can filter out such unstable contributions. The 26-neighborhood connectivity rule can also consider obliquely connected voxels as connected, which is beneficial for maintaining the connectivity of cave arches or inclined cave walls.

[0054] Lithological interface labels are generated in regions where the spatial gradient of cavern evidence counts exceeds a preset gradient threshold. The spatial gradient can be approximated using three-dimensional difference: for each voxel, the difference in cavern evidence counts with adjacent voxels in three coordinate directions is calculated to form the gradient magnitude. .in, Indicates the gradient magnitude; , , These represent the differences in cave evidence counts across the three coordinate directions. A preset gradient threshold of 4 is used, indicating that regions with gradient magnitudes greater than or equal to 4 are labeled as lithological interfaces. The reason for using a gradient threshold is that the evidence distribution near the cave interface changes more sharply, and the difference in cave evidence counts between the cave interior and the surrounding rock area is more pronounced. Gradient filtering can extract the interface location from the transition zone between high and low evidence areas, serving as an interface guide for subsequent variational implicit surface solutions.

[0055] The final cave structure evidence body is constructed using a voxel grid in a unified coordinate system. Each voxel stores cave evidence counts and water evidence counts, as well as cave connectivity labels and lithological interface labels. Cave evidence counts and water evidence counts provide dense three-dimensional statistical evidence, cave connectivity labels provide clues to the connectivity structure of the cave network, and lithological interface labels provide location clues to the transition zones between cave boundaries and surrounding rock. These four types of information together constitute the input that can be directly used for subsequent variational implicit surface solving.

[0056] In another optional implementation, the voxel mesh side length is set to 0.10m for fine modeling. After the cavity evidence count and water evidence count are accumulated, a three-dimensional morphological closure operation is performed to fill the small holes into a continuous set of voxels. Then, three-dimensional connected component labeling is performed to generate cavity connectivity labels. This can improve the integrity of connected components at fine-scale voxels and reduce connected component breaks caused by discrete sampling.

[0057] Through the above implementation process, the orthogonal dipole fully polarized physical separation separates the water body multiple wave component data from the cave body main reflection component data. The extended polarization reflection path weaving algorithm transforms the four-channel polarization information into a co-polarized coherent descriptor, a cross-polarization dominant descriptor, and a polarization rotation consistent descriptor. After reflection event extraction and reflection event chain generation, event chain classification, propagation speed determination, and voxel evidence voting back projection, a cave body structure evidence body containing cave body evidence count, water body evidence count, cave body connectivity label, and lithological interface label is obtained.

[0058] In one implementation, the voxel grid uses a unified coordinate system to define its spatial extent, enabling it to cover the survey corridor corresponding to the four polarization combinations of horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal fully polarized echo data, as well as the cave wall area covered by laser scanning point cloud data. The voxel edge length can be 0.20m, facilitating the inclusion of cave evidence counts, water evidence counts, cave connectivity labels, and lithological interface labels within the same spatial bearing structure, while ensuring continuous surface detail in the subsequently extracted 3D cave mesh model. When the cave morphology changes drastically or finer cave wall textures are required, the voxel edge length can be 0.10m to enhance the representation of cave ceiling arches and cave wall undulations.

[0059] The implicit range field is constructed on a voxel mesh, and the implicit range field is denoted as . ,in To represent a three-dimensional spatial position within a unified coordinate system. This indicates the signed distance from this location to the zero isosurface corresponding to the implicit surface of the cavity. The sign symbol allows the implicit surface of the cavity to be accessed via... The point set description shows that the interior and exterior of the cavity correspond to different symbol intervals. The initial construction of the implicit distance field can directly utilize the cavity connectivity labels: set the voxels within the cavity's connected domains indicated by the connectivity labels to negative values, and set the voxels outside the labels to positive values. Then, perform distance field reinitialization on the symbol field, making the value of each voxel close to its geometric distance to the zero isosurface. This makes the initial implicit cavity surface closely resemble the cavity's connected structure, which is already determined by voxel evidence. The alternating iterative process converges from a geometrically reasonable starting point, reducing the chance of leading the implicit cavity surface to incorrect connected domains. One feasible reinitialization method is to use the fast marching method: take the voxels at the sign change points as the front set, propagate the distance from the front inwards and outwards, and fill the voxel grid, so that the implicit distance field satisfies distance monotonicity in space.

[0060] The overall driving force field consists of the consistency driving force of the cave structure evidence, the consistency driving force of the point cloud, the surface smoothness driving force, and the connectivity preservation driving force. In each round of alternating iterations, the same voxel mesh is used as the computational domain for updating and superimposing. The consistency driving force of the cave structure evidence comes from the cave evidence count and the water body evidence count. Its purpose is to make the implicit surface of the cave tend to surround regions with high cave evidence counts and low water body evidence counts within the cave, while excluding regions with low cave evidence counts or high water body evidence counts from the cave. In actual calculations, a scalar field of cave evidence counts is first formed on the voxel mesh, denoted as... ,in Indicates position The corresponding voxel-based cave evidence count; simultaneously forming a scalar field for water body evidence count, denoted as... ,in Indicates position The water body evidence count corresponds to the voxels. The interface advancement direction is determined by pointing from high-count voxels of cave body evidence count to low-count voxels; this can be calculated using discrete gradients. The opposite direction is used as the direction of propulsion, where Represents the spatial gradient operator, Indicates the count of evidence related to the cave in the location. The gradient vector. Using a direction from high-count voxels to low-count voxels can push the implicit surface of the cavern into the transition zone of cavern evidence counting, making the zero isosurface closer to the boundary region of evidence from high to low. This is consistent with the spatial distribution of the cavern interface after the voxel evidence voting back projection, which shows a higher inner side and a lower outer side. To constrain the advancement direction of water evidence counting, water suppression conditions can be applied to the cavern evidence counting when calculating the advancement direction, for example, only when the conditions are met... voxels are used The reverse direction is used as the direction of interface advancement, while satisfying voxels are used The positive direction is used as the interface advancement direction, which makes the implicit curved surface of the cave more inclined to retreat in the area where the water evidence count is dominant, thus maintaining consistency with the water multiple wave component data obtained by orthogonal dipole full polarization physical separation.

[0061] The driving force for point cloud consistency comes from laser-scanned point cloud data, aiming to ensure that the implicit curved surfaces of the cave body are closely aligned with the geometric observations of the cave walls. This is achieved by first generating a point cloud distance field on a voxel mesh, denoted as... ,in This represents the Euclidean distance from the location to the nearest point in the laser-scanned point cloud data. The calculation of the point cloud distance field can be accelerated using spatial indexing structures, such as establishing a 3D nearest neighbor search structure for the laser-scanned point cloud data and then querying the nearest point distance for each voxel center. To avoid unnecessary interference from the spatial domain far from the cave walls on point cloud consistency, the nearest point search radius can be set to 2.0m. When a voxel center cannot find a nearest neighbor within this radius, the point cloud distance field for that voxel remains at 2.0m. The driving force for point cloud consistency uses the equidistant surface direction of the point cloud distance field, which is also the gradient direction of the point cloud distance field. ,in The direction of increasing distance is indicated. Incorporating the equidistant surface direction of the point cloud distance field into the point cloud consistency driving force allows the implicit surface of the cave to tend towards a position with smaller point cloud distances along the normal direction. This causes the zero isosurface to gradually approach the dense point cloud region on the cave wall, thus locking the geometric position of the implicit surface of the cave to actual cave wall observations. The benefit of this approach is that when cave evidence counts are sparse or widely distributed, the point cloud consistency driving force provides a stable geometric reference, preventing the implicit surface of the cave from drifting in space solely based on back projection evidence.

[0062] The driving force for surface smoothing originates from the local second-order difference of the implicit range field, aiming to suppress voxel-level steps or jagged edges on the implicit surface of the cavity. The local second-order difference can be implemented using a three-dimensional discrete Laplacian operator, denoted as... ,in Represents the discrete Laplace operator. Indicates the location The second-order difference combination of the implicit range field values ​​of the neighboring voxels. A direct implementation is to compute on the six neighborhoods of the voxel grid. That is, to put the position The implicit range field values ​​of adjacent voxels in the three coordinate directions are summed and then subtracted by a multiple of the implicit range field value at the current position. The smooth surface driving force is... The sign and amplitude guide the implicit distance field to evolve towards a smoother shape, reducing local spikes and isolated depressions on the implicit surface of the cavity while maintaining its overall shape. This facilitates a more regular distribution of triangular facets when generating a 3D mesh model of the cavity through subsequent iso-evaluation. To avoid excessively smoothing out the real edges of the cavity walls, the surface smoothing driving force can be activated only in the neighborhood of voxels where the cavity evidence count is below a certain preset threshold. For example, regions with a cavity evidence count less than 3 are smoothed first, allowing more detail to be preserved at interfaces with strong evidence.

[0063] The connectivity maintenance mechanism originates from cavern connectivity tags, aiming to stabilize the cavern connectivity structure during alternating iterations and address breaks or unexpected mergers in cavern connected domains during numerical evolution. Specifically, it involves performing a consistent reinitialization of the implicit distance field within each connected domain based on the cavern connectivity tags: for each cavern connected domain, the implicit distance field symbols of all voxels within the domain are unified to cavern-internal symbols, and voxels outside the domain are unified to cavern-external symbols. Then, a distance field reinitialization is performed on the symbol boundaries, ensuring the implicit distance field re-satisfies its morphology. This prevents symbol flip islands from appearing inside the cavern due to local advancements during alternating iterations, and also prevents symbol leakage from outside the cavern into the connected domain. The connection evolution method for connected domains at contact boundaries addresses situations where narrow channels exist between connected domains but are disconnected due to voxel discretization or insufficient local evidence. After each iteration, boundary voxel pairs between different cavern connected domains are searched. When the distance between boundary voxel pairs is less than or equal to 0.30m, the voxel path between the two boundary voxel pairs is set as a zero-equivalence channel. This means adjusting the implicit distance field value of the voxels along the path to 0, and then re-initializing the distance field in the neighborhood of that channel, ensuring continuous connection of the cavern's implicit surface at that location. Voxel paths between boundary voxel pairs can be generated using a 3D line rasterization algorithm, ensuring continuous path coverage on the voxel mesh. This method is simple to execute and computationally stable.

[0064] refer to Figure 2 This diagram visually illustrates how the cave boundary gradually converges from an initial hypothetical state to the actual geological interface. Figure 2In the computational domain shown, the background is divided into regular two-dimensional grids, representing a slice of a three-dimensional voxel grid. Each grid cell (voxel) stores the cave structure evidence data calculated through previous steps, including cave evidence counts and water body evidence counts. The shades of color in the figure represent the intensity of the cave evidence counts; darker colors indicate a greater likelihood that the location belongs to the cave interior, i.e., stronger co-polarization coherence of the fully polarized radar echo and higher laser point cloud density. The dashed circles in the figure represent the initial state (t=0) of the implicit surface evolution, typically initialized as a simple geometry enclosing the survey trajectory or a coarse envelope generated based on connectivity labels. The solid outline represents the final cave implicit surface (zero isosurface) that converges after multiple iterations.

[0065] The evolution process is controlled by a total driving force field, with the arrows in the figure indicating the direction of the driving force. This driving force field consists of three parts: first, the evidence consistency driving force, which uses the spatial gradient information of the cavity evidence count to move the surface from the transition zone of the low evidence zone to the high evidence zone, i.e., the arrow in the figure points to the high gradient boundary; second, the point cloud consistency driving force, which, based on the distance field constructed from the laser point cloud, pulls the surface toward the actual laser observation point position to correct for possible spatial positioning errors in the radar data; and finally, the surface smoothing driving force, which, based on local second-order difference operators (such as the Laplacian operator), suppresses non-physical noise and spikes generated on the surface during the evolution process. During the iteration process, the computational unit uses the level set method or a similar numerical scheme to update the signed distance function value on the grid node according to the total driving force vector at the current position. As the number of iterations t increases, the implicit surface will automatically deform, split, or merge, eventually locking at the position where the cavity evidence count gradient is the largest and best fits the point cloud, i.e., the convergence state shown by the solid line in the figure. This variational evolution method effectively solves the problem that traditional explicit modeling is unable to handle complex changes in topology (such as branches and loops), and can generate closed, smooth and physically meaningful geological interfaces of caves based on multi-source data fusion.

[0066] An alternating iterative process is used to update the implicit distance field to output the implicit surface of the cavity. In one iteration, the implicit distance field is first updated based on the cavity structure evidence consistency driving force and the point cloud consistency driving force, then smoothed based on the surface smoothing driving force, and finally re-initialized and evolved based on the connectivity preservation driving force. The rationale behind this arrangement is that the cavity structure evidence consistency driving force and the point cloud consistency driving force are responsible for pushing the implicit surface of the cavity to the correct spatial position; the surface smoothing driving force improves the local geometric quality after the position is basically correct; and the connectivity preservation driving force pulls the topology back to the connected framework indicated by the cavity connectivity label after geometric evolution, reducing topological drift. The update of the implicit distance field can adopt a normal-progressive approach: calculating the gradient of the implicit distance field. The normal direction is obtained by normalization. ,in Indicates position The distance field normal direction, Represent the vector norm; then update the propulsion by using the component of the total driving force field in the normal direction. An implementable update syntax is as follows: ,in Indicates the first The implicit distance field values ​​during round iterations Indicates the first The implicit distance field values ​​during round iterations Indicates the step size, Indicates position The total driving force field vector, This represents the vector dot product. The step size can be set to 0.05m, ensuring that the magnitude of each advance is less than the voxel edge length, thereby gradually approximating the target surface while maintaining convergence stability. The number of iterations can be set to 60. The iteration ends when the position change of the zero isosurface is less than 1 voxel edge length in 5 consecutive iterations, and the implicit surface of the cavity is output.

[0067] In an optional implementation, the alternating iterative process can divide the consistency driving force of the cave structure evidence body and the consistency driving force of the point cloud into two sub-loops: first, execute 10 iterations primarily driven by the consistency driving force of the cave structure evidence body, so that the implicit surface of the cave enters the high-gradient transition zone of the cave evidence count; then execute 10 iterations primarily driven by the point cloud consistency driving force, so that the implicit surface of the cave conforms to the point cloud of the cave wall; and then enter a complete loop containing the surface smoothing driving force and the connectivity maintenance driving force. This makes it easier to lock the implicit surface of the cave to the correct cave wall position when the cave evidence count is sparse but the point cloud coverage is sufficient, while reducing the chance of the point cloud consistency driving force pulling the surface toward the point cloud of obstacles outside the cave.

[0068] After outputting the implicit surface of the cavity, a 3D visualization model is generated. When generating the 3D mesh model of the cavity through isosurface extraction, the zero isosurface of the implicit distance field is selected as the extraction target. A 3D isosurface extraction algorithm is used to convert the implicit distance field on the voxel mesh into a triangular mesh. A common implementation involves checking the sign combination of the implicit distance field values ​​of the eight vertices of each voxel cube, generating triangular patches based on the sign combination by looking up a table, and then performing linear interpolation on the vertex positions of the triangular patches to obtain more accurate zero-intersection positions. Linear interpolation can be written as... ,in and This represents the coordinates of the two endpoints of a voxel edge. This represents the coordinates of the intersection point between the zero isosurface and the edge of the voxel. This represents the interpolation coefficients. Interpolation coefficients can be derived from... Received, among which Endpoints The implicit distance field values, Endpoints The implicit distance field values ​​are used. This results in a smoother boundary for the 3D mesh model of the cave, and the triangular faces are closer to the actual curved surface. For example, when the voxel side length is 0.20m, the number of triangular faces in the 3D mesh model of the cave can be between 300,000 and 1,200,000, depending on the spatial scale of the cave and the complexity of the surface.

[0069] When mapping lithological interface labels and water evidence counts to the 3D geological attribute layers of the cave's 3D mesh model, a consistent coordinate reference is maintained during the mapping, and sampling is performed directly from the voxel mesh to the mesh vertices. The mapping of water evidence counts uses trilinear interpolation: for each mesh vertex location, its corresponding cubic voxel cell is found in the voxel mesh, the water evidence counts of the eight vertices of that cube are read, and interpolation is used to obtain the vertex water evidence value. When used for visualization, the water evidence value can be mapped to color, transparency, or texture intensity, creating a visual indication of water distribution in space. To ensure comparability across different cave sections, the water evidence value can be normalized to the maximum water evidence count across the entire field; for example, the maximum water evidence count corresponds to a display intensity of 1.0, and 0 corresponds to a display intensity of 0.0. The mapping of lithological interface labels can use neighborhood determination: when a lithological interface label exists in the voxel containing a mesh vertex or in its 26 neighboring voxels, the interface attribute of that mesh vertex is set as an interface point; otherwise, it is set as a non-interface point. This allows the voxel-level interface markers to be extended onto the curved mesh, making it easier to display the location of lithological interfaces on the surface of the 3D mesh model of the cave in the form of outlines or contour lines, which helps to interpret the transition zone between the cave boundary and the surrounding rock.

[0070] The 3D geological profile output uses preset cutting rules to automatically cut the 3D mesh model of the cave. These preset cutting rules can generate a cutting plane every 10m along a preset survey line direction, or every 15m along the principal axis of the cave's 3D mesh model. The intersection of the cutting plane and the cave's 3D mesh model yields the profile outline, and 3D geological attribute layers are simultaneously sampled along this outline, forming a visualization of profile attributes and water distribution. Setting the cutting plane orthogonal to the preset survey line direction ensures the profile aligns with the survey line interpretation perspective of the fully polarized echo data. This facilitates the identification of changes in the cave roof height, cave wall contraction, and concentrated areas of water evidence within the profile outline, corresponding to different segments of the echo event chain, thus forming a verifiable interpretation loop. For situations where it is necessary to highlight the risk points of a certain tunnel section, the preset cutting rules can also include an optional method of densifying the cutting based on dense areas of lithological interface labels: in sections where lithological interface labels are continuously distributed in a unified coordinate system, the cutting interval is adjusted from 10m to 3m, so that the profile can reveal the spatial location of the interface transition zone in more detail.

[0071] refer to Figure 1 ,exist Figure 1 In the illustrated embodiment, the computer system or data processing terminal has completed the construction of the implicit surface of the cave based on the aforementioned variational implicit surface solution steps, and performed iso-value extraction operations (e.g., the moving cube algorithm), thereby generating a three-dimensional mesh model of the cave composed of a large number of triangular facets. This model is placed in a globally unified coordinate system, which is determined based on the alignment relationship between the survey line trajectory of the fully polarimetric ground-penetrating radar system and the laser scanning point cloud data. Figure 1 This demonstrates in detail how the system performs automatic sectioning to obtain information about the internal geological structure of a specific location. The three-dimensional tubular structure shown in the figure is the generated three-dimensional mesh model of the cave. Its surface color represents different geological properties. For example, dark blue areas correspond to areas with high water body evidence counts, indicating the presence of water accumulation or highly water-bearing media, while brown or gray areas correspond to dry lithological cave walls or surrounding rock matrix. The red plane P1 in the figure represents the preset sectioning plane.

[0072] The definition of the cutting plane P1 is based on the orientation of a preset survey line or a spatial location specified interactively by the user. It is typically configured to be perpendicular to the tangent direction of the survey line or perpendicular to the main axis of the tunnel to obtain a view that best reflects the cross-sectional morphology of the tunnel. During processing, the computational unit first determines the spatial equation of the cutting plane P1, then traverses all triangular faces of the tunnel's 3D mesh model, calculating the intersections of plane P1 with the edges of the triangular faces. Connecting these intersections forms a closed wireframe, i.e., the red profile outline shown in the figure. Importantly, this process involves not only geometric interception but also simultaneous attribute sampling. While calculating the intersections, the system calculates the attribute values ​​of each point on the profile outline using linear interpolation based on the geological attribute values ​​of the mesh vertices (such as lithological interface labels and water evidence counts). The generated profile outline not only contains the geometric clearance information of the tunnel at that location but also carries attribute information on water distribution and lithological variations. Figure 1 The two-dimensional projection shown on the right is a geological profile map after mapping the three-dimensional interception results onto a two-dimensional plane. It clearly marks the water distribution area and the lithological cave wall area. This automated cutting and attribute mapping mechanism enables technicians to quickly extract key cross-sectional data that can be used for engineering drawing and safety assessment from complex full three-dimensional models, realizing the dimensionality reduction analysis from three-dimensional spatial field to two-dimensional engineering drawing.

[0073] The visualization model output includes a 3D mesh model of the cave, a 3D geological attribute layer, and a 3D geological profile. The 3D mesh model of the cave is used for overall spatial display, the 3D geological attribute layer is used to express the spatial distribution of lithological interface labels and water evidence counts, and the 3D geological profile provides a more readable cross-sectional interpretation view. For ease of engineering application, the visualization model can be output as 3D interactive scene data and profile image data. The 3D interactive scene data retains the vertex attributes of the 3D mesh model and the 3D geological attribute layer, while the profile image data retains the overlay rendering result of the profile outline and profile attributes. For example, the profile image data can be output at a resolution of 1920 x 1080, with one profile image output for each cutting plane, and vector data of the profile outline is output simultaneously for direct overlay in planning software.

[0074] refer to Figure 4 This figure illustrates the core 3D view portion of the final graphical user interface delivered to the user. The main element shown is a 3D mesh model of the cavity generated after variational implicit surface solving and isometry extraction. Figure 1 They focus on different cutting processes. Figure 4 The focus is on the comprehensive expression of overall attributes and the presentation of spatial measurement information. The surface texture of this 3D model is not a simple lighting rendering, but a synthetic material generated by mapping multidimensional information from the cave structure evidence. Specifically, the color distribution on the model surface directly corresponds to the physical properties of the underground medium: dark blue or cyan areas are mapped to water evidence counts in the voxel grid through trilinear interpolation, intuitively revealing water-rich areas at the bottom of the cave or in fissures, which serves as a warning for the drainage design of underground engineering; brown or yellowish-brown areas correspond to lithological interface labels or low water count areas, representing stable surrounding rock structures.

[0075] To enhance the measurability of the model, Figure 4The illustrated embodiment displays white contour lines overlaid on a 3D curved surface. These contour lines are vector lines automatically generated by the system through horizontal cropping of the 3D mesh model at preset elevation intervals (e.g., every 1 or 2 meters). They curve with the surface's undulations, allowing users to intuitively judge the slope changes of the cave walls by observing the density of the contour lines—dense contour lines indicate steep cave walls, while sparse lines indicate gentle slopes. Furthermore, the figure also shows a red marker line and the text "current profile position P1," indicating that the visualization system supports interactive operation. Users can dynamically adjust the position of the cutting plane in this 3D view by clicking or dragging with the mouse, and the system will update the 2D profile in the auxiliary view in real time. This visualization method transforms the implicit mathematical model into an intuitive engineering picture, not only displaying the geometry of the cave but also integrating physical parameters retrieved from polarimetric radar. It provides geological analysts with a WYSIWYG digital exploration platform, enabling a comprehensive assessment of the structural stability and water hazard risks of underground cave systems.

[0076] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these specific embodiments are merely illustrative. Those skilled in the art can omit, substitute, and modify the details of the above methods and systems in various ways without departing from the principles and essence of the present invention. For example, combining the above method steps to perform substantially the same function and achieve substantially the same result according to substantially the same method falls within the scope of the present invention. Therefore, the scope of the present invention is defined only by the appended claims.

Claims

1. A method for structural identification and visual modeling of underground karst cave systems, characterized in that, The method includes the following steps: Step 1: Move the data acquisition platform along the pre-set survey line inside the underground cavern and use the fully polarized ground-penetrating radar system to collect fully polarized echo data of four polarization combinations: horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal. The fully polarized ground-penetrating radar system adopts an orthogonal dipole antenna design and sets the same trigger source for the four polarization combinations to generate a unified timestamp; simultaneously acquire laser scanning point cloud data and generate a point cloud coordinate system. Step 2: Generate cave structure evidence body based on the extended polarization reflection path weaving algorithm of orthogonal dipole full polarization physical separation. The cave structure evidence body includes cave evidence count, water evidence count, cave connectivity label and lithological interface label. Step 3: Variational implicit surface solution to generate the implicit surface of the cave: A voxel mesh is established in a unified coordinate system and an implicit distance field is constructed. A total driving force field is constructed, consisting of the consistency driving force of the cave structure evidence volume, the consistency driving force of the point cloud, the surface smoothness driving force, and the connectivity preservation driving force. The implicit distance field is updated using an alternating iterative process to output the implicit surface of the cave. Step 4: Output of 3D visualization modeling: Perform iso-extraction on the implicit curved surface of the cave to generate a 3D mesh model of the cave, map the lithological interface labels and water evidence counts to the 3D geological attribute layer of the 3D mesh model of the cave, and output the 3D geological profile and visualization model.

2. The method as described in claim 1, characterized in that, By calibrating the target, the external parameter relationship between the coordinate system of the fully polarized ground-penetrating radar system and the point cloud coordinate system is established and a unified coordinate system is generated. Based on the external parameter relationship, the fully polarized echo data and the laser scanning point cloud data are mapped to the unified coordinate system. In step 2, the fully polarized echo data is input into the polarization separation module. In the same polarization receiving link of the orthogonal dipole antenna, in-phase synthesis is performed to generate a common polarization synthesis channel. In the cross polarization receiving link, orthogonal phase synthesis is performed to generate a polarization rotation channel. The co-polarization coherent descriptor, the cross-polarization dominant descriptor, and the polarization rotation consistent descriptor are calculated on the co-polarization synthesis channel and the polarization rotation channel, respectively. Based on the co-polarization coherent descriptor, the cross-polarization dominant descriptor, and the polarization rotation consistent descriptor, reflection event chain generation and classification screening are performed, and voxel evidence voting back projection is performed in combination with the line propagation speed.

3. The method according to claim 2, characterized in that, In step two, the four-channel in-situ alignment and baseband construction process of the extended polarization reflection path weaving algorithm includes: performing sampling start point alignment at the same time stamp and generating channel sequences for the four polarization combinations of horizontal-horizontal, vertical-vertical, horizontal-vertical, and vertical-horizontal; performing in-phase component extraction and orthogonal component extraction on each channel sequence to form a complex baseband channel sequence, and arranging the complex baseband channel sequences according to a unified time stamp to form a four-channel complex baseband set; performing in-phase synthesis on the horizontal-horizontal and vertical-vertical complex baseband sets in the in-polarization receiving link as input for the cave main reflection component data; and performing orthogonal phase synthesis on the horizontal-vertical and vertical-horizontal complex baseband sets in the cross-polarization receiving link as input for the water body multiple wave component data.

4. The method according to claim 3, characterized in that, In step two, the polarization structure descriptor construction process includes: for each trace sequence of each survey line, a common-polarization coherent descriptor is calculated on the common-polarization synthesis channel using a sliding time window. The common-polarization coherent descriptor is obtained by performing point-by-point conjugate multiplication on the complex baseband sample values ​​of the common-polarization synthesis channel within the sliding time window, accumulating them within the sliding time window, and then performing amplitude normalization; a cross-polarization dominant descriptor is calculated on the polarization rotation channel using a sliding time window. The cross-polarization dominant descriptor is obtained by taking the amplitude of the complex baseband sample values ​​of the polarization rotation channel within the sliding time window and accumulating them within the sliding time window; a polarization rotation consistent descriptor is calculated on the horizontal-vertical channel and the vertical-horizontal channel using a sliding time window. The polarization rotation consistent descriptor is obtained by taking the phase of the complex baseband sample values ​​of the horizontal-vertical channel and the vertical-horizontal channel within the sliding time window and calculating the proportion of the phase difference within a preset stable range.

5. The method according to claim 4, characterized in that, In step two, the process of extracting reflection events and generating reflection event chains includes: performing envelope generation on the co-polarization synthesis channel and detecting local peaks in each channel sequence to form a cave candidate event set; performing envelope generation on the polarization rotation channel and detecting local peaks in each channel sequence to form a water candidate event set; for the cave candidate event set, establishing a connection between two cave candidate events with adjacent time positions and envelope peak shapes that meet a preset similarity threshold in adjacent channel sequences, and recording the continuously connected event sequences as cave reflection event chains; for the water candidate event set, generating water multiple event chains using the same connection rules as the cave candidate event set, and detecting peak clusters arranged with repeated delays in time position in each water multiple event chain to form water multiple cluster labels.

6. The method according to claim 5, characterized in that, In step two, the event chain classification process based on descriptor and chain consistency includes: extracting continuous high-value segments of co-polarized coherent descriptors within each cave reflection event chain and generating cave interface consistency labels; extracting continuous high-value segments of cross-polarization dominant descriptors and polarization rotation consistent descriptors within each water multiple event chain, and combining them with water multiple cluster labels to generate water multiple consistency labels; selecting a cave interface candidate set from the cave reflection event chain based on the cave interface consistency labels, and selecting a water multiple candidate set from the water multiple event chain based on the water multiple consistency labels.

7. The method according to claim 6, characterized in that, In step two, the process of determining the propagation speed includes: selecting event segments with hyperbolic trajectories from the candidate set of the cave interface, generating theoretical hyperbolas one by one using the candidate propagation speeds, calculating the cumulative time position deviation between the event segments and the theoretical hyperbolas, and selecting the candidate propagation speed with the smallest cumulative deviation as the propagation speed of the survey line.

8. The method according to claim 7, characterized in that, In step two, the voxel evidence voting backprojection process includes: establishing a voxel grid in a unified coordinate system and defining cave evidence counts and water evidence counts for the voxel grid; for each cave candidate event in the cave interface candidate set, converting the propagation distance of the cave candidate event into a spherical neighborhood based on the spatial location of the measuring point and the propagation speed of the measuring line, and performing cave evidence count accumulation within the voxels covered by the spherical neighborhood; for each water candidate event in the water multiples candidate set, performing water evidence count accumulation within the voxels covered by the spherical neighborhood using the same spherical neighborhood backprojection method as for cave candidate events; the cave connectivity label and lithological interface label generation process includes: after completing the backprojection of all measuring lines, generating cave connectivity labels by performing three-dimensional connectivity domain marking based on the cave evidence count, and generating lithological interface labels in areas where the spatial gradient of the cave evidence count is greater than a preset gradient threshold.

9. The method according to claim 1, characterized in that, In step three, the construction process of the total driving force field includes: generating the interface advancement direction based on the direction from high-count voxels to low-count voxels in the cave evidence count and writing the cave structure evidence volume consistency driving force; generating the point cloud distance field on the voxel grid based on the laser scanning point cloud data and writing the equidistant surface direction of the point cloud distance field into the point cloud consistency driving force; generating the surface smoothing driving force based on the local second-order difference of the implicit distance field; performing symbolic reinitialization within the connected domain of the implicit distance field based on the cave connectivity label, and performing connection evolution at the contact boundary of the connected domain to write the connectivity maintenance driving force.

10. The method according to claim 1, characterized in that, In step four, the three-dimensional visualization modeling output process also includes: automatically cutting the three-dimensional mesh model of the cave body in a unified coordinate system according to the preset cutting rules to generate a three-dimensional geological profile, and outputting a visualization result containing the profile outline, profile attributes and water distribution.

Images