A method for automatically constructing a karst cave three-dimensional model based on multi-source point cloud fusion

Abstract

The application belongs to the technical field of detection, and particularly relates to a karst cave three-dimensional model automatic construction method based on multi-source point cloud fusion, which comprises the following steps: step S100, multi-hole multi-source sensor cooperative arrangement and original detection data acquisition; step S200, single-hole detection data preprocessing and karst cave space contour point cloud data extraction; step S300, multi-hole point cloud constrained by control points is registered and globally fused, the karst cave space contour point cloud data of all detection drill holes are registered and fused, and global karst cave point cloud and fusion quality signals are output; and step S400, a karst cave three-dimensional grid model or a fusion quality insufficient prompt information is output. The present application executes a three-level registration fusion strategy, and has high multi-hole registration precision. The present application effectively solves the error accumulation problem in multi-hole data fusion, ensures uniform distribution of overall fusion precision, and is suitable for large-scale chain-bead string-shaped karst cave detection scenes with many caves and many drill holes.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering investigation and underground space exploration technology, specifically a method for automatically constructing a three-dimensional model of a karst cave based on multi-source point cloud fusion. Background Technology

[0002] Karst geology is a special geological phenomenon formed by the long-term dissolution of soluble rocks (mainly carbonate rocks) by water, widely distributed in southwestern, southern, and eastern my country. In karst-developed areas, underground caves are a core geological risk source affecting the safety of engineering construction. The presence of caves significantly reduces the bearing capacity and stability of the overlying soil and rock mass. Triggered by external factors such as building loads, changes in groundwater dynamics, or seismic disturbances, they can easily lead to ground subsidence, uneven foundation settlement, and even structural damage to engineering structures, posing a serious threat to the construction and operation of major infrastructure such as roads, bridges, tunnels, subways, and high-rise buildings. Therefore, accurate exploration of underground caves and the establishment of high-precision three-dimensional geometric models in the early stages of engineering construction are important prerequisites for conducting engineering geological evaluation, cave treatment design, and construction safety control in karst areas.

[0003] Existing cave exploration and 3D modeling technologies have the following main shortcomings: First, the fusion of multi-site cloud data lacks effective constraints. Chain-like karst caves are distributed in series along their development direction, requiring data collection from multiple exploratory boreholes. Since the data from each borehole exists in an independent local coordinate system, without a unified coordinate reference constraint, the gradual accumulation of measurement errors when integrating the site clouds into the global coordinate system leads to significant translational and rotational deviations between different site clouds, ultimately causing overall distortion of the 3D model.

[0004] Secondly, the morphological characteristics of chain-like karst caves have not been fully utilized. Existing methods mostly treat karst caves as general irregular cavities for modeling, without taking into account the prior morphological characteristics of chain-like karst caves, such as the approximate ellipsoidal shape of the cave chambers and the interconnectedness of multiple chambers. This results in deficiencies in the geometric accuracy and structural integrity of the 3D reconstruction model, and it is also unable to automatically extract geometric parameters such as the cave volume and connecting passage parameters required for engineering.

[0005] To address the aforementioned issues, it is necessary to propose an automatic construction method for karst cave 3D models based on multi-source point cloud fusion. Summary of the Invention

[0006] To address the aforementioned shortcomings in existing technologies, this invention provides an automatic method for constructing a 3D model of a karst cave based on multi-source point cloud fusion, comprising the following steps: Step S100, Multi-source sensor deployment and raw detection data acquisition: The location data of the cave space is acquired by multi-source sensors deployed on the ground and in each detection borehole, and the raw dataset of the cave boundary is obtained. Step S200: Preprocessing of single-hole detection data and extraction of point cloud data of cave spatial contour: Step S201: Perform television image data processing and cave cavity segment identification on the image frame sequence data in the original dataset of the cave boundary to obtain the candidate location of the first cavity boundary; Step S202: The well diameter profile sequence and induction logging curve data in the original dataset of the cave boundary are subjected to noise filtering and cave spatial point cloud generation in sequence to obtain the cave spatial profile point cloud data and the second cave boundary candidate position in the local three-dimensional coordinate system with the borehole opening of each detection borehole as the origin. Step S203: Determine the final depth range of the cave segment by using the candidate positions of the first and second cave boundaries; extract and label the cave spatial contour point cloud data, add quality weight attributes to each point, and output the cave spatial contour point cloud data carrying quality weights. Step S300: Align and globally fuse the point cloud of the multi-hole constrained by the control point, align and fuse the point cloud data of the spatial contour of the cave in all the probe boreholes, eliminate the translational offset and rotational deviation between the local coordinate systems of each hole through the global binding adjustment algorithm, and output the global cave point cloud and fusion quality signal. Step S400: Reconstruction of the three-dimensional morphology of the cave based on the chain-like constraint: Select the processing path according to the fusion quality signal, automatically segment the cave chambers and fit the ellipsoid parameters on the global cave point cloud, perform model integrity repair and automatic calculation of engineering parameters, and output the three-dimensional mesh model of the cave or the prompt message of insufficient fusion quality.

[0007] Preferably, step S201 includes: For the image frame sequence data of each borehole i, each frame is converted to the HSV color space, and the luminance channel V is extracted. Significant jumps in luminance values ​​at the boundary between the borehole wall and the cavity due to changes in reflection characteristics are detected in the image frame sequence data. Along the depth direction, the luminance values ​​at each depth layer Z are analyzed. k Luminance channel mean sequence v k The first-order difference is calculated as follows: ; Where k is the depth layer identifier; v k For depth layer Z k The brightness value extracted from the brightness channel; v k-1 For depth layer Z k-1 The brightness value extracted from the brightness channel; The depth layer k where the absolute value of the first-order difference exceeds the adaptive threshold is marked as the candidate location of the first hole boundary, as shown in the following formula: ; in An adaptive threshold; First-order difference The mean; First-order difference The standard deviation.

[0008] Preferably, step S202 includes: For the original dataset of cave boundaries after filtering out impulse noise, outliers, and low-frequency drift noise in the depth direction, it is divided into several depth layers according to the detection depth z. A consistency check of the radial distance sequence is performed between adjacent depth layers of each borehole i. The Z values ​​of each adjacent depth layer are then calculated. k and Z k+1 The average boundary recognition distance r(i, Z) at the same angle θ k ,θ) and r(i,Z) k+1 The difference between the contour abrupt change difference and the average boundary recognition distance r(i, Z) is used to obtain the contour abrupt change difference; if the contour abrupt change difference is less than or equal to the average boundary recognition distance r(i, Z), then the contour abrupt change difference is obtained. k ,θ) and r(i,Z) k+1 If the difference between the average values ​​of θ and Z is greater than 0.25, then the adjacent depth layers Z are considered to be adjacent. k and Z k+1 There is a significant abrupt change in angle θ between them, where k is the depth layer number symbol; If a significant abrupt change occurs simultaneously in a direction exceeding 30% of the circumference of the depth layer, the depth layer is determined to be a candidate location for the top or bottom interface of the cave, and the depth set of the candidate interfaces is recorded to obtain the candidate location for the second cavity boundary.

[0009] Preferably, step S203 includes: Add a quality weight attribute to each point in the point cloud and obtain the quality weight. As shown in the following formula: ; Where r(i, z, θ) is the average boundary identification distance of each detection borehole i at each detection depth z along each angle θ; in The median of the average boundary identification distance for each probe borehole i at all probe depths z along all angles θ; The numerical range of quality weight is 0 to 1. The smaller the deviation from the median value, the closer the weight is to 1, and the higher the reliability of the measurement point.

[0010] Preferably, step S300 includes: Step S301: Perform initial coordinate transformation based on the borehole control points. The coordinates of each point in the local coordinate system point cloud of each borehole location are initially mapped to the global geodetic coordinate system through a translation transformation based on the local coordinate system of the borehole location to the geodetic coordinate system, thus completing the initial transformation. Step S302: Point clouds of the same chamber or the same connecting channel detected by adjacent boreholes are spatially adjacent in the global coordinate system, but residual deviations exist due to the accumulation of measurement errors. These need to be further eliminated by a refinement transformation based on control point pairs with the same name. In the overlapping detection area of ​​adjacent boreholes i and i+1, a set of control point pairs with the same name is constructed. Step S303: Calculate the root mean square value of the residual coordinate difference between all control point pairs after the final transformation, and generate the fusion quality signal accordingly.

[0011] Preferably, step S302 includes solving the weighted least squares refined rigid transformation based on the control point pairs using the SVD method, as shown in the following equation: ; Where R is the orthogonal rotation matrix to be optimized; t is the translation vector to be optimized; Where w j Let be the confidence weight of the j-th control point pair, where j is the sequential number of the control point, j=1,2,...,m, m is the total number of control points, and m≥6; Among them (R) * , t * ) represents the optimal refined rigid transformation parameters, where R * Let t be a 3x3 orthogonal rotation matrix, where t * Let be the translation vector; where Let j be the three-dimensional geodetic coordinates of the same control point j in the i-th exploratory borehole; Let j be the three-dimensional coordinates in the geodetic coordinate system of the same control point j in the (i+1)th exploratory borehole. and Together they form a pair of control points with the same name.

[0012] Preferably, step S303 includes: When the root mean square value of the residual coordinate difference is less than the preset voxel side length of 0.05m, a high fusion quality signal S=1 is output; When the root mean square value of the residual coordinate difference is less than three times the preset voxel side length (0.15m) and greater than the preset voxel side length (0.05m), the output fusion quality signal S=2; When the root mean square value of the residual coordinate difference is greater than or equal to three times the preset voxel side length (0.15m), an insufficient fusion quality signal S=3 is output, and a hole-filling command is triggered.

[0013] Preferably, step S400, selecting the processing path based on the fusion quality signal, includes: If the fused quality signal is S=1, then the DBSCAN density clustering algorithm is used directly based on the global cave point cloud to automatically segment the cave. If the fused quality signal is S=2, then first perform point cloud smoothing processing based on moving least squares on the global cave point cloud, set the search radius to 0.10m, set the polynomial order to 2nd order, and then use the DBSCAN density clustering algorithm to automatically segment the cave. The DBSCAN core distance parameter is adaptively set to 0.3 times the average distance between adjacent boreholes, and the minimum number of points is set to 20. If the fusion quality signal is S=3, the cavern segmentation process is terminated, and a prompt message indicating that the borehole layout density and depth are unreasonable is output.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention employs a systematic, hierarchical noise processing approach, ensuring data quality. Targeting three typical types of noise in cave exploration data, this invention designs a serial, hierarchical processing flow that includes impulse noise and outlier filtering, depth-direction low-frequency drift noise filtering, and cross-depth layer consistency verification. Furthermore, a quantified quality weight attribute is added to each measurement point, enabling the effective transfer of data reliability information between subsequent processing steps, thus laying a solid data foundation for high-precision 3D reconstruction.

[0015] This invention employs a three-level alignment fusion strategy, achieving high registration accuracy across multiple borehole locations. It organically combines three levels: initial translation transformation constrained by borehole geodetic coordinates, weighted SVD-refined rigid registration based on pairs of corresponding feature control points, and global bundled adjustment optimization. This approach utilizes borehole coordinates to provide global constraints, preventing errors in chain-like registration from concentrating at the end of the chain. Furthermore, feature registration eliminates rotational deviations introduced by measurement errors, effectively solving the error accumulation problem in multi-borehole data fusion. This ensures uniform distribution of overall fusion accuracy and is suitable for large-scale chain-like karst cave exploration scenarios with numerous chambers and boreholes. Attached Figure Description

[0016] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the point cloud of the karst cave space proposed in the embodiments of the present invention; Figure 3 This is a schematic diagram of the three-dimensional geometric mesh model of the karst cave proposed in the embodiments of the present invention. Detailed Implementation

[0017] To better understand the present invention, the following description, in conjunction with the accompanying drawings and examples, will further illustrate the content of the present invention.

[0018] like Figure 1As shown, this application relates to an automatic construction method for a 3D model of a karst cave based on multi-source point cloud fusion, comprising the following steps: Step S100, Multi-source sensor deployment and raw detection data acquisition: The location data of the cave space is acquired by multi-source sensors deployed on the ground and in each detection borehole, and the raw dataset of the cave boundary is obtained.

[0019] In this embodiment, based on the chain-like morphological characteristics of the karst cave, several exploratory boreholes are laid out around the cave at preset intervals along the cave's development direction. These boreholes are sequentially numbered with the identifier i, where i = 1, 2, ..., n; n is the total number of exploratory boreholes. The ground geodetic coordinates of each borehole are obtained, i.e., the borehole opening geodetic coordinates. Three types of detection equipment are placed inside each borehole: a borehole television camera, an ultrasonic caliper, and a three-component induction logging tool.

[0020] During the descent, the three types of detection equipment gradually collect raw data of the cave boundary at a preset step size of 0.5m. Each time the preset step size is reached, raw data of the cave boundary at that depth is collected in a 360-degree circumference on the horizontal plane, and the detection depth z corresponding to the raw data of the cave boundary is recorded simultaneously, until the bottom of the detection borehole is reached.

[0021] The borehole television camera uses 360° panoramic optical imaging to acquire images of the borehole wall and cave walls, obtaining image frame sequence data. The frame rate of the borehole television camera is no less than 25fps, and the resolution is no less than 1920×1080 pixels.

[0022] The ultrasonic wellbore measuring instrument measures the dimensions of the cave contour in various directions by reflecting sound pulses, thus obtaining a wellbore contour sequence. The ultrasonic wellbore measuring instrument has an angular resolution of 1° and a distance measurement accuracy better than ±3mm.

[0023] Among them, the three-component induction logging tool uses electromagnetic induction to detect the conductivity distribution in the cave space, assisting in judging the filling status and spatial boundary of the cave, and obtaining induction logging curve data.

[0024] Preferably, the estimated spacing of the exploratory boreholes is 10 meters.

[0025] Preferably, considering the beaded morphological characteristics of the karst cave, the preset spacing of the detection boreholes is 0.3 times the spacing between the chambers of the beaded karst cave.

[0026] In this embodiment, the acquisition of raw exploration data also includes: image frame sequence data, well diameter profile sequence data, and induction logging curve data collected by all devices in each exploration borehole are indexed by the exploration borehole sequence number i and the exploration depth z, and are uploaded in real time to the ground industrial control computer via a wired transmission bus. The data is then stored in a structured binary format on a solid-state drive array, forming a raw dataset of cave boundaries organized by borehole number. Each raw dataset of cave boundaries contains image frame sequence data, well diameter profile sequence data, and induction logging curve data from all exploration devices within the corresponding i-th borehole across the entire depth range.

[0027] It should be noted that step S100 uses the geodetic coordinates of the borehole opening as the reference control information throughout all subsequent steps. This design provides a reliable geometric constraint anchor point for the alignment and fusion of data from different borehole locations within a unified coordinate framework. The chambers of the chain-like karst cave are spatially discrete, and the spatial information of the chambers detected from different borehole locations belongs to their respective independent local coordinate systems. Without precise borehole opening control coordinates, the fusion of data from each borehole location will produce a significant systematic offset, leading to distortion in the final 3D model.

[0028] Step S200, preprocessing of single-hole detection data and extraction of cave spatial contour point cloud data, includes: Step S201: Perform television image data processing and cave cavity segment identification on the image frame sequence data in the original dataset of the cave boundary to obtain the candidate location of the first cavity boundary.

[0029] Step S202: Noise filtering and spatial point cloud generation are performed on the well diameter profile sequence and induction logging curve data in the original dataset of the karst cave boundary to obtain the spatial profile point cloud data of the karst cave and the candidate position of the second cavity boundary in the local three-dimensional coordinate system with the borehole opening of each detection borehole as the origin.

[0030] Step S203: Determine the final depth range of the cave segment by using the candidate positions of the first and second cave boundaries; extract and label the cave spatial contour point cloud data, add quality weight attributes to each point, and output the cave spatial contour point cloud data carrying quality weights.

[0031] Specifically, the system receives the raw dataset of cave boundaries from each borehole acquired in step S100. It then processes the image frame sequence data in the raw dataset using television image data and identifies cave cavity segments to obtain the first candidate location for the cave boundary. Preprocessing is performed on the wellbore contour sequence and induction logging curve data in the raw dataset, including noise filtering and cave spatial point cloud generation, to obtain the cave spatial contour point cloud data and the second candidate location for the cave boundary in a local three-dimensional coordinate system with the borehole opening as the origin. The final cave cavity segment depth range is determined using the first and second candidate locations, resulting in the final cave boundary. Finally, the system extracts and labels the cave spatial contour point cloud data.

[0032] The process of obtaining the candidate location of the first cavity boundary includes: using a method combining threshold segmentation based on HSV color space and depth direction gradient detection on the image frame sequence data to identify the depth of the top and bottom interfaces of the cavity segment and obtain the candidate location of the first cavity boundary.

[0033] Specifically, the process of this application for processing television image data and identifying cave cavity segments in the original dataset of image frame sequences of cave boundaries is as follows: a method combining threshold segmentation based on HSV color space and depth direction gradient detection is used to identify the top and bottom interface depths of cave cavity segments from the image sequence.

[0034] Specifically, for the image frame sequence data of each borehole i, it is converted to the HSV color space frame by frame, and the luminance channel V is extracted. The image frame sequence data shows a significant jump in luminance values ​​at the boundary between the borehole wall and the cavity due to changes in reflection characteristics. Along the depth direction, the luminance values ​​at each depth layer Z are analyzed. k Luminance channel mean sequence v k The first-order difference is calculated as follows: ; Where k is the depth layer identifier, v k and v k-1 The depth layer Z is respectively k and Z k-1 The brightness value extracted from the brightness channel.

[0035] The depth layer k where the absolute value of the first-order difference exceeds the adaptive threshold is marked as the candidate location of the first hole boundary, as shown in the following formula: ; Among them An adaptive threshold; Among them and First-order difference The mean and standard deviation.

[0036] Among them, the candidate locations for the second cavity boundary include: For the original dataset of cave boundaries, after filtering out impulse noise, outliers, and low-frequency drift noise in the depth direction, it is divided into several depth layers according to the detection depth z. A consistency check of the radial distance sequence is performed between adjacent depth layers in each borehole i. The difference in the average boundary recognition distance at the same angle θ between adjacent depth layers is calculated to obtain the profile abrupt change difference. If the difference between the profile abrupt change difference and the average boundary recognition distance is greater than 0.25, it is determined that there is a significant abrupt change at that angle between the adjacent depths.

[0037] If a significant abrupt change is detected in more than 30% of the directions within a certain depth layer, then that depth layer is determined to be a candidate location for the top or bottom interface of the cave, and the depth set of the candidate interfaces is recorded to obtain the candidate location for the second cavity boundary.

[0038] Candidate locations of the top and bottom interfaces of the cave are identified, and candidate locations of the second cavity boundary are obtained.

[0039] Specifically, in this embodiment, the preprocessing process for the well diameter profile sequence and induction logging curve data in the original dataset of the karst cave boundary is as follows: Noise filtering and spatial point cloud generation are performed independently on the original dataset of the cave boundary for each exploratory borehole. Noise filtering specifically includes: for the original dataset of the cave boundary of the i-th exploratory borehole, impulse noise, outlier filtering, low-frequency drift noise filtering in the depth direction, and cross-depth layer consistency verification are performed respectively.

[0040] It should be noted that the original dataset of karst cave boundaries typically contains three types of noise: impulse noise caused by multipath reflection of sound waves, short-range anomalies caused by local fracturing of the rock wall, and low-frequency drift noise introduced by mechanical vibration of equipment. Corresponding filtering strategies are applied to each of these three types of noise, and they are processed sequentially.

[0041] The specific process for filtering out impulse noise and outliers is as follows: For each exploratory borehole i, the average boundary identification distance along the 360-degree circumference of the cavern boundary is calculated using the wellbore profile sequence data and induction logging curve data at each exploration depth z. This yields the average boundary identification distance r(i, z, θ) along each angle θ∈(0, 2π) at each exploration depth z. The median r(i, z) and the absolute median difference (MAD) of the average boundary identification distances along all angles θ at each exploration depth z for each exploratory borehole i are then calculated. r (i, z).

[0042] The average boundary recognition distance that meets the following conditions will be judged as impulse noise or an outlier and marked as invalid: ; All average boundary recognition distances marked as invalid are removed and replaced with linear interpolations of the average boundary recognition distances of adjacent angles θ to ensure the integrity of the full 360° contour data of each depth layer.

[0043] in This is the preset impulse noise removal factor, with a default value of 3.5.

[0044] It should be noted that when the impulse noise rejection factor is set to 3.5, it can effectively remove strong outliers with a deviation exceeding 3.5 times the absolute median, while retaining reasonable radial distance abrupt changes caused by geometric changes in the cave wall.

[0045] The specific process for filtering out low-frequency drift noise in the depth direction is as follows: The average boundary identification distance of each borehole i is obtained along the detection depth z at the same angle θ. The Savitzky-Golay filter is used to filter the data sequence for each angle θ. The filter window length is set to 11 depth sampling layers, the polynomial fitting order is set to 3, and the filter output is used as the drift estimate of the angle θ at the depth z.

[0046] It should be noted that, as the ultrasonic caliper and three-component induction logging tool descend along the depth direction, the data sequence of the cave boundary identification distance at the same angle θ shows a gradual change in the continuous rock mass. However, the low-frequency drift introduced by mechanical vibration will be superimposed on the true signal of the cave boundary identification distance, forming a trend deviation.

[0047] The specific process of cross-depth layer consistency verification is as follows: For the original dataset of cave boundaries, after filtering out impulse noise, outliers, and low-frequency drift noise in the depth direction, it is divided into several depth layers according to the detection depth z. A consistency check of the radial distance sequence is performed between adjacent depth layers of each borehole i. The Z values ​​for each adjacent depth layer are calculated. k and Z k +1 is the average boundary recognition distance r(i,Z) at the same angle θ. k ,θ)and r(i,Z) k+1 The difference between the contour abrupt change difference and the average boundary recognition distance r(i, Z) is obtained. k ,θ)and r(i,Z) k+1 If the difference between the average values ​​of (θ) and (θ) is greater than 0.25, then the adjacent depth layer Z is determined to be... k and Z k+1There is a significant abrupt change at this angle θ, where k is the depth layer number.

[0048] If a significant abrupt change is detected in more than 30% of the directions within a certain depth layer, then that depth layer is determined to be a candidate location for the top or bottom interface of the cave, and the depth set of the candidate interfaces is recorded to obtain the candidate location for the second cavity boundary.

[0049] As one implementation of this embodiment, the intersection of the first cavity boundary candidate position and the second cavity boundary candidate position is taken to obtain the final cavity boundary, thereby determining the final depth range of the karst cavity segment.

[0050] In this embodiment, the specific process of extracting and quality labeling the spatial contour point cloud data of the karst cave is as follows: like Figure 2 As shown, in this embodiment, after filtering out impulse noise, outliers, low-frequency drift noise in the depth direction, and cross-depth layer consistency of the original dataset of the cave boundary of a certain borehole i, the subsequent spatial point cloud generation of the cave is performed, and the points of each depth layer Z are then processed. k The cleaned full-circumference radial distance sequence is converted into a local three-dimensional coordinate point set. A local coordinate system is established with the borehole opening as the origin, north as the X-axis, east as the Y-axis, and vertical downward as the Z-axis. All local three-dimensional coordinate points of the borehole are then input into this local coordinate system.

[0051] Based on the determined depth range of the cave cavity segment, a subset of the point cloud within the cave cavity segment is extracted from the original contour point cloud to form the cave spatial contour point cloud. Simultaneously, a quality weight attribute is assigned to each point in the point cloud, calculated based on the reciprocal of the deviation of each measurement value from the median in the first-stage impulse noise detection. The calculation formula is as follows: ; Where r(i,z,θ) is the average boundary identification distance of each detection borehole i at each detection depth z along each angle θ∈(0,2π); in The median of the average boundary identification distance for each probe borehole i at all probe depths z along all angles θ; The quality weight ranges from 0 to 1. The smaller the deviation from the median, the closer the weight is to 1, reflecting a higher reliability of the measurement point. The quality weight attribute will be passed to step S300 along with the point cloud data and will be directly used in the subsequent calculation of the reliability weight of the control points, realizing the effective transfer of single-well point cloud quality information between steps.

[0052] Finally, the point cloud of the spatial contour of the karst cave at each borehole location is stored in .pcd format, with the borehole location number i, the geodetic coordinates of the borehole opening, the final cavity boundary, the depth range of the karst cavity section, and the mass weight of each point all transferred to step S300.

[0053] Step S300: Alignment and global fusion of multi-hole point cloud data based on control point constraints. Alignment and data fusion are performed on the spatial contour point cloud data of all exploration boreholes. Translational offset and rotational deviation between local coordinate systems of each borehole are eliminated by a global binding adjustment algorithm. The global cave point cloud and fusion quality signal are output.

[0054] Specifically, the point cloud data of the spatial contour of the karst caves from all the boreholes are aligned and fused. A global binding adjustment algorithm is used to eliminate translational and rotational deviations between the local coordinate systems of each borehole to obtain the overall spatial contour point cloud data of the karst caves. The global karst cave point cloud and fusion quality signal are then output.

[0055] In this embodiment, the point cloud data of the karst cave spatial contour of each borehole, the geodetic coordinates of each borehole opening, and the mass weight of each point are received from step S200. Through control point constraint transformation and a global binding adjustment algorithm, the translational offset and rotational deviation between the local coordinate systems of each borehole are resolved. All borehole point clouds are then fused into a unified global geodetic coordinate system, and the global karst cave point cloud and fused quality signal are output. The specific process is as follows: Step S301: Perform initial coordinate transformation based on the borehole control points. The coordinates of each point in the local coordinate system point cloud of each borehole location are initially mapped to the global geodetic coordinate system through a translation transformation based on the local coordinate system of the borehole location to the geodetic coordinate system, thus completing the initial transformation.

[0056] Step S302: Point clouds detected by adjacent boreholes in the same cave or connecting passage are spatially adjacent in the global coordinate system, but residual deviations exist due to accumulated measurement errors. These deviations need to be further eliminated through a refinement transformation based on control point pairs with the same name features. In the overlapping detection area of ​​adjacent boreholes i and i+1, i.e., in the connected sections of the karst caves or adjacent caves jointly detected by the two boreholes, a set of control point pairs with the same name is constructed based on the spatial intersection point of the centroid trajectory line of the connecting passage cross-section.

[0057] In this embodiment, the number of effective control point pairs is required to be greater than or equal to 6, and the confidence weight of each control point is directly taken as the average value of the quality weights of each measurement point in the neighborhood of the control point output in step S200.

[0058] Based on control point pairs, the SVD method is used to solve the weighted least squares refined rigid transformation: ; Where R and t are the orthogonal rotation matrix and translation vector to be optimized, respectively; Where wj is the confidence weight of the j-th control point pair. The specific value is directly derived from the quality weight in step S200. j is the sequential number of the control point, j=1,2,...,m, m is the total number of control points, and m≥6.

[0059] Among them (R) * , t * ) represents the optimal refined rigid transformation parameters, which is the combination of the rotation matrix and translation vector that minimizes the objective function, where R * A 3x3 orthogonal rotation matrix is ​​used to eliminate rotational deviations between adjacent pore point clouds, where t * It is a translation vector used to eliminate translational deviations between adjacent pore point clouds.

[0060] in Let j be the three-dimensional geodetic coordinates of the same control point j in the i-th exploratory borehole; Let j be the three-dimensional coordinates in the geodetic coordinate system of the same control point j in the (i+1)th exploratory borehole. and Together, they form a pair of control points with the same name. Theoretically, they should describe the same physical feature point in the cave space (such as the centroid of the cross section of the connecting passage), and the difference in their coordinates is the residual registration deviation that needs to be eliminated through refinement transformation.

[0061] As one implementation method of this embodiment, after iteratively processing all adjacent borehole pairs to form a chain registration, to avoid the accumulation of errors along the chain direction, a global bundled adjustment optimization model containing the constraints of all borehole control points is further constructed. The model uses the weighted sum of squared residuals of all control point pairs as the objective function to simultaneously optimize the refinement transformation parameters of all boreholes, distributing the error evenly across all boreholes. After the global bundled adjustment converges, the point clouds of each borehole after refinement transformation are merged, and overlapping redundant points are removed by voxel filtering. The voxel side length is set to 0.05m to form the global cave point cloud P.

[0062] Step S303: Calculate the root mean square value of the residual coordinate difference between all control point pairs after the final transformation, and generate the fusion quality signal accordingly.

[0063] When the root mean square value of the residual coordinate difference is less than the preset voxel side length of 0.05m, a high fusion quality signal S=1 is output; When the root mean square value of the residual coordinate difference is less than three times the preset voxel side length (0.15m) and greater than the preset voxel side length (0.05m), the output fusion quality signal S=2; When the root mean square value of the residual coordinate difference is greater than or equal to three times the preset voxel side length (0.15m), an insufficient fusion quality signal S=3 is output, and a hole-filling command is triggered. The triggering of the hole replenishment command indicates that the density and depth of the detection boreholes are unreasonable, resulting in distortions, mainly translation and rotation, in the data fusion process. The number of boreholes should be increased or the detection equipment should be monitored to ensure data accuracy.

[0064] It should be noted that this step combines the initial translation transformation constrained by the borehole geodetic coordinates with three levels of alignment processing: refined SVD registration based on corresponding feature control point pairs, and global bundled adjustment. This constitutes a complete solution for multi-borehole data alignment and correction. Simply relying on the translation transformation of the borehole coordinates cannot eliminate the rotational deviation introduced by the accumulation of measurement errors, while relying solely on feature registration faces the risk of mismatch when the geometries of chain-like karst caves are highly similar. Combining the two allows for the organic integration of the global constraints provided by the borehole coordinates and the local accuracy advantages of feature registration. Global bundled adjustment further solves the error accumulation problem of chain registration, which is particularly crucial for large-scale chain-like karst cave exploration scenarios with many caves and boreholes, ensuring that the overall fusion accuracy is evenly distributed rather than concentrated at the end of the chain.

[0065] Step S400: Reconstruction of the three-dimensional morphology of the karst cave based on the chain-like constraint: Select the processing path according to the fusion quality signal, automatically segment the cave chambers and fit the ellipsoid parameters on the global karst cave point cloud, perform model integrity repair and automatic calculation of engineering parameters, and output the three-dimensional mesh model of the karst cave.

[0066] Specifically, step S400 involves the reconstruction of the three-dimensional morphology of the cave based on the chain-like constraint; and the completion of missing surfaces, removal of rough edges, fitting of ellipsoidal parameters, rendering, and marking of geometric parameters for the overall cave spatial contour point cloud data based on the global cave point cloud and the fusion quality signal.

[0067] like Figure 3 As shown, in this embodiment, the global cave point cloud and fusion quality signal output in step S300 are received. Combining the prior morphological constraints of the beaded cave system, automatic cave segmentation, ellipsoidal parameter fitting, surface reconstruction, and global mesh stitching are performed on a GPU computing platform to generate a three-dimensional geometric mesh model of the cave.

[0068] As one implementation method of this embodiment, the processing path is first selected based on the fusion quality signal. If the fusion quality signal is S=1, then the cavern segmentation is performed directly based on the global cavern point cloud. If the fusion quality signal is S=2, firstly, point cloud smoothing based on moving least squares is performed on the global cave point cloud to eliminate the influence of residual errors introduced by moderate fusion quality on surface reconstruction. Then, the DBSCAN density clustering algorithm is used to automatically segment the global point cloud into several cave sub-point clouds and several connected channel point cloud segments. For each cave sub-point cloud, based on the prior constraint of the beaded cave's approximate ellipsoid, a weighted least squares ellipsoid fitting method is used to estimate the cave's geometric center, three-axis radii, and attitude rotation matrix. The weights of the weighted least squares ellipsoid fitting method are directly used from the output of step S200 and passed through step S300 to retain the quality weights. After completing the above operations, cave segmentation is performed.

[0069] In one implementation of this embodiment, the search radius for point cloud smoothing based on the moving least squares method is set to 0.10m, and the polynomial order is 2. The core distance parameter of the DBSCAN density clustering algorithm is adaptively set to 0.3 times the average distance between adjacent boreholes, and the minimum number of points is set to 20. If the fusion quality signal is S=3, then no cavity segmentation will be performed, and the following information will be output: The hole layout density and layout depth are unreasonable.

[0070] In this embodiment, after the cavern segmentation is completed, model integrity repair, comprehensive accuracy evaluation, automatic calculation of engineering parameters, and final output are performed to form a complete modeling result package for engineering applications, outputting a karst cave network model. The karst cave mesh model undergoes topological integrity detection, automatically locating and repairing geometric defects such as hanging edges, non-manifold vertices, and cavities. Cavity repair employs a surface completion method based on radial basis function (RBF) interpolation, using the cavity boundary normal vector as a constraint to ensure the continuity of the repaired surface with the surrounding geometry. After repair, the volume of each cavern and the overall engineering parameters of the karst cave are calculated based on the cavern ellipsoidal geometric descriptor.

[0071] The final results are uploaded to the engineering database server in three types of files: 3D models in obj and ply formats, engineering parameter reports in XML format (including the volume of each cavern, spatial coordinates, series direction, and connecting passage parameters), and comprehensive evaluation reports in PDF format, thus completing the entire process of automatically constructing the 3D model of the karst cave.

[0072] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0073] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0074] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0075] These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable apparatus for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0076] The above are merely embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention are included within the scope of the claims of the present invention pending approval.

Claims

1. A method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion, characterized in that, The method includes the following steps: Step S100, Multi-source sensor deployment and raw detection data acquisition: The location data of the cave space is acquired by multi-source sensors deployed on the ground and in each detection borehole, and the raw dataset of the cave boundary is obtained. Step S200: Preprocessing of single-hole detection data and extraction of point cloud data of cave spatial contour: Step S201: Perform television image data processing and cave cavity segment identification on the image frame sequence data in the original dataset of the cave boundary to obtain the candidate location of the first cavity boundary; Step S202: The well diameter profile sequence and induction logging curve data in the original dataset of the cave boundary are subjected to noise filtering and cave spatial point cloud generation in sequence to obtain the cave spatial profile point cloud data and the second cave boundary candidate position in the local three-dimensional coordinate system with the borehole opening of each detection borehole as the origin. Step S203: Determine the final depth range of the cave segment by using the candidate positions of the first and second cave boundaries; extract and label the cave spatial contour point cloud data, add quality weight attributes to each point, and output the cave spatial contour point cloud data carrying quality weights. Step S300: Align and globally fuse the point cloud of the multi-hole constrained by the control point, align and fuse the point cloud data of the spatial contour of the cave in all the probe boreholes, eliminate the translational offset and rotational deviation between the local coordinate systems of each hole through the global binding adjustment algorithm, and output the global cave point cloud and fusion quality signal. Step S400: Reconstruction of the three-dimensional morphology of the cave based on the chain-like constraint: Select the processing path according to the fusion quality signal, automatically segment the cave chambers and fit the ellipsoid parameters on the global cave point cloud, perform model integrity repair and automatic calculation of engineering parameters, and output the three-dimensional mesh model of the cave or the prompt message of insufficient fusion quality.

2. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 1, characterized in that, Step S201 includes: For the image frame sequence data of each borehole i, each frame is converted to the HSV color space, and the luminance channel V is extracted. Significant jumps in luminance values ​​at the boundary between the borehole wall and the cavity due to changes in reflection characteristics are detected in the image frame sequence data. Along the depth direction, the luminance values ​​at each depth layer Z are analyzed. k Luminance channel mean sequence v k The first-order difference is calculated as follows: ； Where k is the depth layer identifier; v k For depth layer Z k The brightness value extracted from the brightness channel; v k-1 For depth layer Z k-1 The brightness value extracted from the brightness channel; The depth layer k where the absolute value of the first-order difference exceeds the adaptive threshold is marked as the candidate location of the first hole boundary, as shown in the following formula: ； in An adaptive threshold; First-order difference The mean; First-order difference The standard deviation.

3. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 1, characterized in that, Step S202 includes: For the original dataset of cave boundaries after filtering out impulse noise, outliers, and low-frequency drift noise in the depth direction, it is divided into several depth layers according to the detection depth z. A consistency check of the radial distance sequence is performed between adjacent depth layers of each borehole i. The Z values ​​of each adjacent depth layer are then calculated. k and Z k+1 The average boundary recognition distance r(i, Z) at the same angle θ k ,θ) and r(i,Z) k+1 The difference between the contour abrupt change difference and the average boundary recognition distance r(i, Z) is used to obtain the contour abrupt change difference; if the contour abrupt change difference is less than or equal to the average boundary recognition distance r(i, Z), then the contour abrupt change difference is obtained. k ,θ) and r(i,Z) k+1 If the difference between the average values ​​of θ and Z is greater than 0.25, then the adjacent depth layers Z are considered to be adjacent. k and Z k+1 There is a significant abrupt change in angle θ between them, where k is the depth layer number; If a significant abrupt change occurs simultaneously in a direction exceeding 30% of the circumference of the depth layer, the depth layer is determined to be a candidate location for the top or bottom interface of the cave, and the depth set of the candidate interfaces is recorded to obtain the candidate location for the second cavity boundary.

4. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 1, characterized in that, Step S203 includes: Add a quality weight attribute to each point in the point cloud and obtain the quality weight. As shown in the following formula: ； Where r(i, z, θ) is the average boundary identification distance of each detection borehole i at each detection depth z along each angle θ; in The median of the average boundary identification distance for each probe borehole i at all probe depths z along all angles θ; The numerical range of quality weight is 0 to 1. The smaller the deviation from the median value, the closer the weight is to 1, and the higher the reliability of the measurement point.

5. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 1, characterized in that, Step S300 includes: Step S301: Perform initial coordinate transformation based on the borehole control points. The coordinates of each point in the local coordinate system point cloud of each borehole location are initially mapped to the global geodetic coordinate system through a translation transformation based on the local coordinate system of the borehole location to the geodetic coordinate system, thus completing the initial transformation. Step S302: Point clouds of the same chamber or the same connecting channel detected by adjacent boreholes are spatially adjacent in the global coordinate system, but residual deviations exist due to the accumulation of measurement errors. These need to be further eliminated by a refinement transformation based on control point pairs with the same name. In the overlapping detection area of ​​adjacent boreholes i and i+1, a set of control point pairs with the same name is constructed. Step S303: Calculate the root mean square value of the residual coordinate difference between all control point pairs after the final transformation, and generate the fusion quality signal accordingly.

6. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 5, characterized in that, Step S302 includes solving the weighted least squares refined rigid transformation based on the control point pairs using the SVD method, as shown in the following equation: ； Where R is the orthogonal rotation matrix to be optimized; t is the translation vector to be optimized; Where w j Let be the confidence weight of the j-th control point pair, where j is the sequential number of the control point, j=1,2,...,m, m is the total number of control points, and m≥6; Among them (R) * , t * ) represents the optimal refined rigid transformation parameters, where R * Let t be a 3x3 orthogonal rotation matrix, where t * Let be the translation vector; where Let j be the three-dimensional geodetic coordinates of the same control point j in the i-th exploratory borehole; Let j be the three-dimensional coordinates in the geodetic coordinate system of the same control point j in the (i+1)th exploratory borehole. and Together they form a pair of control points with the same name.

7. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 5, characterized in that, Step S303 includes: When the root mean square value of the residual coordinate difference is less than the preset voxel side length of 0.05m, a high fusion quality signal S=1 is output; When the root mean square value of the residual coordinate difference is less than three times the preset voxel side length (0.15m) and greater than the preset voxel side length (0.05m), the output fusion quality signal S=2; When the root mean square value of the residual coordinate difference is greater than or equal to three times the preset voxel side length (0.15m), an insufficient fusion quality signal S=3 is output, and a hole-filling command is triggered.

8. The method for automatically constructing a 3D model of a karst cave based on multi-source point cloud fusion as described in claim 1, characterized in that, Step S400, selecting the processing path based on the fusion quality signal includes: If the fused quality signal is S=1, then the DBSCAN density clustering algorithm is used directly based on the global cave point cloud to automatically segment the cave. If the fused quality signal is S=2, then first perform point cloud smoothing processing based on moving least squares on the global cave point cloud, set the search radius to 0.10m, set the polynomial order to 2nd order, and then use the DBSCAN density clustering algorithm to automatically segment the cave. The DBSCAN core distance parameter is adaptively set to 0.3 times the average distance between adjacent boreholes, and the minimum number of points is set to 20. If the fusion quality signal is S=3, the cavern segmentation process is terminated, and a prompt message indicating that the borehole layout density and depth are unreasonable is output.

Images