Tunnel three-dimensional geologic body modeling method and system based on in-hole camera shooting and while-drilling lithology perception

Abstract

The invention provides a tunnel three-dimensional geologic body modeling method and system based on in-hole camera shooting and while-drilling lithology perception, and relates to the technical field of tunnel and underground engineering and intelligent geological information processing.The method comprises the steps that firstly, a drill hole wall optical image sequence and while-drilling engineering parameters are synchronously collected and preprocessed respectively; then establishing a depth-time mapping relation to realize accurate alignment of the multi-source data in a depth domain; then fracture intelligent identification and parameter extraction are realized through an improved U-Net network, and lithology classification is completed by using an XGBoost model; feature level and decision level fusion is carried out on the fracture features and the lithologic features, and rock mass quality grades are generated; and finally, constructing an integrated geologic model fusing lithologic distribution and a three-dimensional fracture network by adopting an implicit modeling and discrete fracture network technology. According to the method, the problem of splitting of in-hole visual information and while-drilling physical information is solved, automation and precision of geological recognition are improved, and reliable geological information support is provided for tunnel construction safety and support design.

Description

Technical Field

[0001] This invention relates to the fields of tunnel and underground engineering and intelligent geological information technology, and in particular to a method and system for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling-while-drilling lithology perception. Background Technology

[0002] The safe and efficient construction of tunnels and underground engineering projects depends on a precise understanding of the geological conditions of the rock mass ahead. The stability and deformation characteristics of the surrounding rock of a tunnel are jointly controlled by the rock mass fracture structure and lithological properties. Accurately identifying fracture characteristics, determining lithological distribution, and establishing a three-dimensional geological model to achieve integrated, refined identification and three-dimensional representation of these two aspects is of great significance for tunnel construction safety and support design.

[0003] While borehole imaging technology can provide continuous visual images of the borehole wall, current analysis heavily relies on manual interpretation, resulting in low efficiency and strong subjectivity. Although some studies have attempted to use image processing algorithms for automatic crack identification, there are still shortcomings in terms of accuracy, continuity, and geometric parameter extraction for scenarios with uneven lighting, blurred images, and complex crack morphologies.

[0004] Measurement while drilling (MWD) technology can collect parameters such as drill pressure, rotational speed, mechanical drilling rate, and torque in real time during drilling. These parameters are closely related to rock hardness and the degree of fracture development. Existing methods mostly rely on empirical formulas or statistical learning models to classify and predict lithology, but their accuracy is limited by the selection of sample features and the generalization ability of the model. Especially in strata with complex fracture development, lithology identification and geological interpretation are often ambiguous. For example, an increased drilling rate may indicate the encounter of a fracture zone or it may mean that the lithology has softened, so it is difficult to make a unique and accurate geological interpretation based on the parameters alone.

[0005] Currently, while borehole camera images can provide visual information about rock mass structure, they cannot directly reflect the physical properties of the rock mass, such as hardness, integrity, and strength. This makes it difficult to establish a quantitative relationship between lithology and mechanical properties based solely on images. Drilling parameters reflect indirect mechanical signals and lack a correspondence with image structural features, hindering the spatial fusion of geological information and creating a gap between "visual information" and "physical information." The lack of effective multi-source data fusion algorithms prevents the two types of data from complementing each other, hindering the construction of a unified and accurate geological interpretation model. Furthermore, existing 3D geological modeling methods rely solely on interpolation based on discrete core samples or geological profiles, failing to effectively integrate and reflect the 3D spatial distribution characteristics of key structural planes controlling rock mass stability. Summary of the Invention

[0006] This invention proposes a method and system for modeling three-dimensional geological bodies of tunnels based on borehole photography and drilling-while-drilling lithology perception. It aims to solve the problem that borehole photography images and drilling-while-drilling measurement data are difficult to integrate in traditional applications, resulting in the inability to complement the advantages of the two types of data and to build a unified and accurate geological interpretation model.

[0007] In a first aspect, the present invention provides a method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling-while-drilling lithology sensing, specifically including the following steps:

[0008] Step S1: Use an in-hole camera to acquire a sequence of optical images of the borehole wall and perform image preprocessing to generate a continuous unfolded image of the borehole wall;

[0009] Step S2: Use the measurement while drilling system to synchronously collect engineering parameters during the drilling process and perform data preprocessing;

[0010] Step S3: Establish a depth-time mapping relationship, resample the preprocessed drilling parameters from the time domain to the depth domain, so that they correspond one-to-one with the borehole wall unfolding diagram in the depth domain, and obtain depth-aligned multi-source data;

[0011] Step S4: Perform semantic segmentation of the hole wall unfolded diagram, and calculate the geometric parameters and three-dimensional spatial coordinates of each crack based on the segmentation results to form a crack feature dataset;

[0012] Step S5: Based on the depth-aligned multi-source data, extract time-domain and frequency-domain features, input them into a machine learning model for lithology classification, and output the lithology classification results along the borehole depth direction;

[0013] Step S6: The fracture feature dataset and the lithology classification results are fused at the feature level and the decision level to generate a rock mass quality grade along the borehole depth direction;

[0014] Step S7: Based on the fracture feature dataset, the lithology classification results, and the rock mass quality grading, an integrated model that includes both lithological entity distribution and a three-dimensional fracture network is constructed using implicit modeling technology and a discrete fracture network generation algorithm, and then visualized in three dimensions.

[0015] Furthermore, the specific method of step S1 includes:

[0016] S1-1: A high-resolution borehole imaging system is used. The probe is placed into the borehole at a constant speed of 2m / min and then pulled out at the same speed. During the pulling process, the probe continuously acquires a 360° annular image sequence of the borehole wall at fixed depth intervals. The acquisition software synchronously records the depth coordinates corresponding to each frame of the image.

[0017] S1-2: Convert the acquired raw image from the RGB color space to the HSV color space, keeping the chroma (H) and saturation (S) components unchanged, and only perform multi-scale Retinex processing on the luminance component (V).

[0018] S1-3: A bilateral filtering algorithm is used to eliminate residual noise in the image, preserve the edge structure information of the crack, and further improve the image quality;

[0019] S1-4: Between adjacent image frames, the SIFT algorithm is used for key point detection, and the kd-tree nearest neighbor search algorithm is used to quickly match the feature descriptors of the two frames.

[0020] S1-5: Based on the collocation set, use the RANSAC algorithm to estimate the homography matrix H between the two images;

[0021] S1-6: Using the homography matrix H, all images are projected onto a unified planar coordinate system. For overlapping areas, a multi-band mixing algorithm is used for fusion to preserve the image details to the maximum extent and generate an initial panoramic unfolded image.

[0022] S1-7: Perform contrast-limited adaptive histogram equalization (CLAHE) and adaptive gamma correction on the initial panoramic unfolded image to enhance image details, and perform standardized cropping on the image.

[0023] Furthermore, the specific method of step S2 includes:

[0024] S2-1: Engineering parameters during the drilling process are synchronously acquired using a measurement-while-drilling (MWD) system, including drill pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP). The acquisition frequency of these WWD parameters is preferably no less than 1 Hz, and the sampling interval is preferably 1 second. The sampling frequency can be adjusted according to the probe retraction rate. With the desired depth resolution According to the relation The calculation determines the frequency, and it can be adjusted within the range of 1-10Hz to adapt to different drilling speeds and resolution requirements;

[0025] S2-2: Adopt The standard deviation criterion is used to detect outliers in the data, and linear interpolation is used to repair the data, eliminating extreme anomalies caused by instrument vibration, signal loss, or sudden geological changes during drilling.

[0026] S2-3: Apply the SG filter to the time series of each parameter to smooth the signal fluctuation trend while suppressing high-frequency noise;

[0027] S2-4: To maintain the temporal synchronization and physical correlation among the smoothed parameters, the same SG filtering parameters are used to perform multi-parameter synchronous smoothing on the four sequences of drilling pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP), constructing a synchronous smoothing matrix that includes time coordinates and the physical values ​​of the smoothed parameters. ;

[0028] S2-5: Eliminate the differences in dimensions and numerical ranges between different parameters, and perform Z-score standardization on the smoothed data to obtain the standardized drilling parameter matrix.

[0029] Furthermore, the specific method of step S3 includes:

[0030] Based on the physical principles of drilling, the drilling depth is the integral of the mechanical drilling rate of power (ROP) with respect to time; therefore, a mapping model from the time domain to the depth domain is established:

[0031]

[0032] in, From start to time Total depth drilled, For mechanical rotation speed, For integration variables;

[0033] In discrete data processing, the trapezoidal numerical integration method is used for solving; for the th Each time point, and its corresponding depth value The calculation formula is as follows:

[0034]

[0035] in, For the first Drilling depth at each time step For the first Drilling depth at each time step For the first Mechanical drilling speed per time step For the first Mechanical drilling speed per time step The time interval between two adjacent time steps;

[0036] Then, linear interpolation is used to resample all parameters at equal intervals:

[0037]

[0038] in, For target depth Interpolation results of parameters at the location, For the first time to take pictures inside the hole A depth coordinate point, , Known adjacent depth points (satisfying) ), These correspond to the depth parameter values ​​respectively;

[0039] Finally, precise fusion and alignment of "visual information" and "physical information" in the depth domain were achieved, outputting a multi-source data matrix of depth and alignment. .

[0040] Furthermore, the specific method of step S4 includes:

[0041] S4-1: Introduce the Convolutional Attention Module (CBAM) that combines channel attention and spatial attention mechanisms to construct a U-Net convolutional neural network with an encoder-decoder structure. The encoder part extracts multi-scale features, and the decoder gradually restores the spatial resolution.

[0042] S4-2: A hybrid loss function combining binary cross-entropy loss and Dice loss is used to optimize the model and improve its generalization ability.

[0043] S4-3: Use the trained semantic segmentation model to infer the unfolded diagram of the hole wall and obtain the crack probability map. The optimal global threshold for the image is automatically calculated using the maximum inter-class variance method. The probabilistic map is transformed into a binary segmentation map, and the segmentation results are optimized by combining morphological processing.

[0044] S4-4: The high-quality binary image after morphological processing is refined using the structural surface averaging method to obtain the skeleton image. Then, the initial skeleton is optimized and denoised using the curve slope iteration method to remove abnormal points in the skeleton image and further enhance the smoothness of the fracture skeleton.

[0045] Furthermore, the specific method of step S5 includes:

[0046] S5-1: Extract time-domain and frequency-domain features from the aligned drill pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP) within a sliding window to construct a high-dimensional feature vector;

[0047] S5-2: The time-domain and frequency-domain features extracted from each drilling parameter are concatenated to construct a high-dimensional feature vector for each depth window. And it is standardized using Z-score;

[0048] S5-3: The XGBoost ensemble learning model is used to classify lithology based on the standard feature matrix.

[0049] Furthermore, the specific method of step S6 includes:

[0050] Output fracture characteristic data table Perform statistical aggregation along the borehole depth direction, sliding within the same depth window aligned with the drilling parameters, from... Extract all depth coordinates For cracks falling within the current window, data filtering is performed. Based on the filtered set of cracks, the statistical characteristics of the cracks in that window are calculated to form a crack feature vector. Furthermore, the fracture characteristic data table Transformed into a sequence distributed along the depth. ;

[0051] Fusing fracture feature vectors at the feature level and lithological characteristic vectors Construct joint feature vector ;

[0052] At the decision level, a random forest model is constructed to solve the complex multi-feature, nonlinear classification problem of rock mass quality grading, and the output is a multi-class probability vector of rock mass quality grades. The final rock mass quality grade is determined by the category with the highest probability. The rock mass quality is classified into five grades: Class I rock: hard rock; Class II rock: relatively hard rock; Class III rock: medium hard rock; Class IV rock: relatively soft rock; Class V rock: soft soil;

[0053] The fusion and identification results of all depth points are integrated to generate a continuous profile along the borehole depth direction, providing accurate borehole constraints for subsequent 3D geological modeling.

[0054] Furthermore, the specific method of step S7 includes:

[0055] The fused interpretation results from multiple boreholes are uniformly registered according to their spatial coordinates to form a three-dimensional structured dataset: ;

[0056] A three-dimensional lithological solid model is generated using an implicit modeling method, for each lithological category. Constructing implicit functions:

[0057]

[0058] in, For coordinate variables in space, For lithological categories The number of borehole constraint points, For the first The weight of each control point For radial basis functions RBF, No. The three-dimensional coordinates of each control point;

[0059] By implicit function voxelization and isosurface reconstruction, a continuous and smooth three-dimensional lithological entity distribution model is formed. ;

[0060] A deterministic fracture 3D reconstruction is performed. Using the 3D dip angle, 3D dip direction, 3D spatial trace length, and center point coordinates, the fracture is reconstructed as a deterministic fracture surface in 3D space and added to the deterministic fracture set. ;

[0061] Statistical simulations of stochastic fractures were conducted. Based on the attitude and trace length data of all deterministic fractures, their probability distribution model was statistically analyzed. According to the statistically obtained distribution model and fracture density field, a large number of simulated fractures were randomly generated using the Monte Carlo method, forming a set of stochastic fractures. ;

[0062] Deterministic fracture set and random fracture set The two components are merged to form a complete three-dimensional discrete fracture network model. ;

[0063] Finally, the lithological entity distribution model With three-dimensional discrete fracture network model Deep integration is achieved to realize a unified three-dimensional spatial expression of lithology, structural planes, and rock mass quality parameters, and to provide a visual representation.

[0064] Secondly, the present invention also provides a system for a method of modeling three-dimensional geological bodies of tunnels based on borehole photography and drilling-while-drilling lithology sensing, comprising:

[0065] The multi-source data acquisition and preprocessing module is used to simultaneously acquire and preprocess in-hole camera images and drilling parameters;

[0066] The depth domain data alignment and fusion module is used to achieve accurate alignment of in-hole camera images and drilling parameters in the depth domain;

[0067] The intelligent crack identification and parameter extraction module is used to automatically identify cracks and extract their parameters based on the borehole wall unfolding diagram.

[0068] The intelligent lithology sensing and classification module is used for automatic lithology classification based on drilling parameters.

[0069] The multi-source information fusion and rock mass quality classification module is used to fuse fracture and lithological information and perform rock mass quality classification.

[0070] The integrated 3D geological modeling and visualization module is used to construct and visualize an integrated 3D geological model that incorporates lithology and fracture networks.

[0071] Compared with the prior art, the present invention has the following advantages:

[0072] This invention overcomes the problem of the separation between "visual information" and "physical information" in traditional methods by precisely aligning and fusing the visual structural information provided by in-hole camera images with the physical and mechanical information reflected by drilling measurement data in the depth domain. This achieves efficient fusion and complementarity of multi-source data. It improves the fracture semantic segmentation network, realizes the automatic extraction of fracture geometric parameters, and uses an integrated learning model to classify the lithology of drilling parameters to solve the problem of multiple solutions in lithology interpretation, thereby improving the intelligence level of fracture identification and lithology perception. Through a dual fusion mechanism at the feature level and decision level, it achieves deep fusion and refined evaluation of rock mass information, constructs a high-precision integrated three-dimensional geological model, and promotes the intelligent and automated process of tunnel geological exploration.

[0073] Based on the implementation methods provided in the above aspects, this application can be further combined to provide more implementation methods. Attached Figure Description

[0074] The above and other objects, features, and advantages of exemplary embodiments of the present invention will become readily apparent upon reading the following detailed description with reference to the accompanying drawings. In the drawings, several embodiments of the invention are illustrated by way of example and not limitation, with the same or corresponding reference numerals denoteing the same or corresponding parts, wherein:

[0075] Figure 1 This is a flowchart of the method of the present invention;

[0076] Figure 2 This is a schematic diagram of the preprocessing of images captured by the camera inside the hole.

[0077] Figure 3 This is a schematic diagram of drilling parameter acquisition and preprocessing;

[0078] Figure 4 A schematic diagram of the Convolutional Attention Module (CBAM);

[0079] Figure 5 A schematic diagram of the improved U-Net semantic segmentation network structure;

[0080] Figure 6 Flowchart for XGBoost lithology classification;

[0081] Figure 7 This is a schematic diagram illustrating the deep integration of fractures and lithology.

[0082] Figure 8 This is a schematic diagram illustrating the construction of a three-dimensional geological model according to the present invention;

[0083] Figure 9This is a schematic diagram of the system of the present invention. Detailed Implementation

[0084] The exemplary embodiments disclosed in this application will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of this application are shown in the drawings, it should be understood that this application can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of this application and to fully convey the scope of this application to those skilled in the art. Unless otherwise specified, the technical means used in the embodiments are conventional means well known to those skilled in the art.

[0085] Combination Figure 1 As shown, a method for intelligent identification and 3D modeling of tunnel rock mass based on multi-source data fusion is described. The specific steps of the method include:

[0086] Step 1: Using an in-hole camera, a continuous sequence of optical images of the borehole wall is collected along the tunnel borehole. In the image preprocessing process, the Retinex algorithm is used for low-light enhancement, bilateral filtering is applied for noise suppression, and the crack edges are preserved. Based on SIFT feature matching and RANSAC algorithm, image geometric correction and panoramic stitching are performed to generate a continuous borehole wall unfolded image.

[0087] Step 2: Using a measurement-while-drilling system, synchronously acquire engineering parameters during the drilling process, including drill pressure, rotational speed, torque, and mechanical drilling speed. During data preprocessing, Savitzky-Golay (SG) filtering is used to smooth and denoise the time-series data, and Z-score normalization is used to eliminate the influence of dimensions.

[0088] Step 3: Establish a depth-time mapping relationship. Through interpolation algorithms, resample the preprocessed drilling parameters from the time domain to equally spaced depth coordinates, so that they correspond one-to-one with the borehole camera images in the depth domain, achieving accurate alignment of multi-source data.

[0089] Step 4: The preprocessed hole wall unfolded image is segmented into pixels by constructing and training a U-Net semantic segmentation network. Then, the dip angle, trace length, width, density and three-dimensional spatial coordinates of each crack are automatically extracted by skeletonization and connected component analysis morphological algorithms.

[0090] Step 5: Based on the aligned depth domain parameter data, extract time domain and frequency domain features; input the feature vectors into the XGBoost ensemble learning model for training and prediction, and output the lithology classification results along the borehole depth direction;

[0091] Step 6: Deeply fuse the identification results of fractures and lithology: Perform feature-level fusion, concatenating the fracture feature vector and the lithology feature vector into a joint feature vector; perform decision-level fusion, using a random forest model for comprehensive reasoning based on the joint feature vector to generate a refined rock mass quality classification.

[0092] Step 7: Integrate the fusion interpretation data of all boreholes, use implicit modeling technology and discrete fracture network generation algorithm to construct an integrated model that includes both lithological entity distribution and three-dimensional fracture network, and perform three-dimensional visualization.

[0093] Specifically, such as Figure 2 As shown, step 1 involves using an in-hole camera to acquire a continuous sequence of optical images of the borehole wall along the tunnel borehole. During image preprocessing, the Retinex algorithm is used for low-light enhancement, bilateral filtering is applied for noise suppression, and crack edges are preserved. Image geometric correction and panoramic stitching are performed based on SIFT feature matching and the RANSAC algorithm to generate a continuous unfolded borehole wall image. The specific method includes:

[0094] S1-1: A high-resolution intra-hole imaging system is used. The probe is lowered into the borehole at a constant speed of 2 m / min and then pulled out at the same speed. During the pulling process, the probe continuously acquires a 360° annular image sequence of the borehole wall at fixed depth intervals, forming a sequence frame. The acquisition software synchronously records the depth coordinates corresponding to each frame of the image: This establishes a "depth-image" index relationship, laying the foundation for subsequent alignment with the depth domain of the drilling data;

[0095] S1-2: Convert the acquired raw image from the RGB color space to the HSV color space, keeping the chroma (H) and saturation (S) components unchanged, and only process the luminance (Retinex).

[0096] The RGB model can be converted to the HSV model using the following formula:

[0097]

[0098]

[0099]

[0100] Where V is the luminance component, S is the saturation, H is the chroma, and R, G, B are the color space component values.

[0101] The chromaticity (H) and saturation (S) components are retained, while only the luminance (V) component undergoes MSR enhancement to reduce luminance unevenness and reflective interference.

[0102]

[0103] in, pixel coordinates The reflection component at that location, For scale quantity These are the weighting coefficients. For the original input image at pixel points The brightness value at that location, For the first Gaussian wrapping function of scale;

[0104] For the enhanced Linear contrast stretching is performed to map its grayscale range to [0,255], and then it is combined with the preserved chroma H and saturation S components to finally convert it back to the RGB color space, resulting in an enhanced image with uniform illumination.

[0105] S1-3: A bilateral filtering algorithm is used to eliminate residual noise in the image, preserve the edge structure information of the crack, and further improve the image quality.

[0106] The filter calculation formula is as follows:

[0107]

[0108] in, For target pixel The output grayscale value, To normalize the weights, For pixels The scope of the field For pixels in the neighborhood The original grayscale value, For spatial domain weight kernel, Kernel weights for the range;

[0109] S1-4: Keypoint detection is performed using the SIFT algorithm between adjacent image frames, and the feature descriptors of the two frames are quickly matched using the kd-tree nearest neighbor search algorithm.

[0110] The SIFT algorithm is used to detect scale- and rotation-invariant keypoints, and a 128-dimensional feature description vector is calculated for each keypoint.

[0111] A kd-tree nearest neighbor search algorithm is used for fast feature vector matching. To eliminate false matches, Lowe's ratio test is employed to retain matching pairs where the ratio of the nearest neighbor to the second nearest neighbor distance is less than 0.8, thereby obtaining a high-precision set of matching points. ;

[0112] S1-5: Based on the collocation set, estimate the homography matrix H between the two images using the RANSAC algorithm:

[0113] Based on the matching point set The RANSAC algorithm is used to estimate the homography matrix H between two frames of images.

[0114]

[0115]

[0116] Where H is the homography matrix, To find the matrix that minimizes the subsequent objective function, The origin point, for Corresponding to the transformed points;

[0117] S1-6: Using the homography matrix H, all images are projected onto a unified planar coordinate system. For overlapping areas, a multi-band mixing algorithm is used for fusion to preserve image detail information to the maximum extent, generating an initial panoramic unfolded image.

[0118] Based on the homography matrix H, all images are projected into a unified planar coordinate system. For overlapping regions, a multi-band mixing algorithm is used for fusion. This involves weighted fusion at different frequency band levels of the Laplacian pyramid, as shown in the following formula:

[0119]

[0120] in, The combined layers of the Laplace pyramid. The fusion weights for the Laplace pyramid layers. The two Laplace pyramid layers to be merged

[0121] Effectively eliminates stitching seams while preserving complete image details, generating an initial panoramic unfolded image;

[0122] S1-7: Perform contrast-limited adaptive histogram equalization (CLAHE) and adaptive gamma correction on the initial panoramic unfolded image to enhance image details, and perform normalized cropping on the image:

[0123] The initial panoramic unfolded image is subjected to contrast-limited adaptive histogram equalization (CLAHE) with a block size of 8×8 and a contrast limit threshold of 2.0. Bilinear interpolation is used to eliminate block artifacts.

[0124] The grayscale values ​​of all pixels in the CLAHE-processed image are normalized to the [0,1] interval, summed, and then divided by the total number of pixels in the image to obtain the overall average brightness of the image. The formula is as follows:

[0125]

[0126] in, This represents the total number of pixels in the image. For the first The grayscale value after normalization of each pixel

[0127] According to the average brightness Adaptive determination of gamma correction value Applied to each pixel of the image The corrected pixel values ​​are obtained. ;

[0128] The effective region of the image after CLAHE processing and adaptive gamma correction is located and cropped into a rectangular image of uniform width to form a standardized hole wall unfolded image.

[0129] Specifically, such as Figure 3 As shown, step 2 involves using a measurement-while-drilling system to synchronously acquire engineering parameters during the drilling process, including drill pressure, rotational speed, torque, and mechanical drilling speed. During data preprocessing, Savitzky-Golay (SG) filtering is used to smooth and denoise the time-series data, and Z-score normalization is used to eliminate the influence of dimensions. Specifically, this includes:

[0130] S2-1: Utilize a measurement-while-drilling (MWD) system to synchronously acquire multi-dimensional engineering parameters during the drilling process, including weight on bit (WOB), rotational speed (RPM), torque (TOR), and rate of penetration (ROP). The acquisition frequency of these WWD parameters is preferably no less than 1 Hz, and the sampling interval is preferably 1 second. The sampling frequency can be adjusted according to the probe retraction rate. With the desired depth resolution According to the relation The calculation determines the frequency, and it can be adjusted within the range of 1-10Hz to adapt to different drilling speeds and resolution requirements;

[0131] S2-2: Adopt The standard deviation criterion is used to detect outliers in the data, and linear interpolation is employed to repair the data, eliminating extreme anomalies caused by instrument vibration, signal loss, or sudden geological changes during drilling.

[0132] use Criterion: Calculate the mean of the parameter sequence. and standard deviation For each data point in the sequence Its standardized residual is ,like Then determine This is an outlier;

[0133] For data points marked as outliers, linear interpolation is used to repair them and ensure data continuity. The formula is:

[0134]

[0135] in, and It is the nearest normal data point before and after the outlier;

[0136] S2-3: Apply an SG filter to the time series of each parameter to smooth the signal fluctuation trend while suppressing high-frequency noise.

[0137] To preserve the physical trend of the drilling signal and suppress high-frequency noise, SG filtering is applied to the time series data, for a window width of... The order of the polynomial is The SG filter, which is located at any position within the window. The fitted value is expressed as:

[0138]

[0139] in, This refers to the local coordinate index within the current sliding window. For index The value obtained by polynomial fitting at that point. The coefficients to be determined are obtained by minimizing the fitting error of all data points within the window. To determine;

[0140] The smooth output of the window center point is:

[0141]

[0142] in, For global location index, For local coordinate indices within the sliding window, For the window radius, These are the filter weight coefficients. For global position The original parameter values ​​collected at the location;

[0143] S2-4: To maintain the temporal synchronization and physical correlation among the smoothed parameters, the same SG filtering parameters are used to perform multi-parameter synchronous smoothing on the four sequences of drilling pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP), constructing a synchronous smoothing matrix that includes time coordinates and the physical values ​​of the smoothed parameters:

[0144] The four signals—whole drill weight (WOB), rotational speed (RPM), torque (TOR), and rate of drilling (ROP)—are subjected to synchronous SG filtering using the same filtering window and polynomial order, and a synchronous smoothing matrix is ​​applied. :

[0145]

[0146] S2-5: Eliminate the differences in dimensions and numerical ranges between different parameters, and perform Z-score standardization on the smoothed data to obtain the standardized drilling parameter matrix:

[0147] Due to the significant differences in the physical dimensions and numerical ranges of the drilling parameters, Z-score normalization was used to smooth the matrix to ensure numerical comparability in subsequent modeling and feature extraction. Normalize each parameter column:

[0148]

[0149] in, The standardized dimensionless value. For parameters at time points The smoothed physical value, This is the sample mean of the parameter sequence. This is the sample standard deviation of the parameter sequence;

[0150] After standardization, a standardized parameter matrix is ​​formed. :

[0151]

[0152] Specifically, step 3, establishing the depth-time mapping relationship, involves resampling the preprocessed drilling parameters from the time domain to equally spaced depth coordinates using an interpolation algorithm. This ensures a one-to-one correspondence between the preprocessed drilling parameters and the borehole camera images in the depth domain, achieving precise alignment of multi-source data. This includes:

[0153] Based on the physical principles of drilling, borehole depth is the integral of the rate of penetration (ROP) with respect to time. Therefore, a mapping model from the time domain to the depth domain is established:

[0154]

[0155] in, From start to time Total depth drilled, For mechanical rotation speed, For integration variables;

[0156] In discrete data processing, the trapezoidal numerical integration method is used for solving the problem. For the th Each time point, and its corresponding depth value The calculation formula is as follows:

[0157]

[0158] in, For the first Drilling depth at each time step For the first Drilling depth at each time step For the first Mechanical drilling speed per time step For the first Mechanical drilling speed per time step The time interval between two adjacent time steps;

[0159] Then, linear interpolation is used to resample all parameters at equal intervals:

[0160]

[0161] in, For target depth Interpolation results of parameters at the location, For the first time to take pictures inside the hole A depth coordinate point, , Known adjacent depth points (satisfying) ), These correspond to the depth parameter values ​​respectively;

[0162] Finally, precise fusion and alignment of "visual information" and "physical information" in the depth domain were achieved, outputting a multi-source data matrix of depth and alignment. :

[0163]

[0164] Specifically, such as Figure 4-5 As shown, step 4 involves performing pixel-level crack segmentation on the preprocessed hole wall unfolded image by constructing and training a U-Net semantic segmentation network. Subsequently, morphological algorithms such as skeletonization and connected component analysis are used to automatically extract the dip angle, trace length, width, density, and three-dimensional spatial coordinates of each crack. Specifically, this includes:

[0165] S4-1: Introducing a convolutional attention module CBAM that combines channel attention and spatial attention mechanisms to construct a U-Net convolutional neural network with an encoder-decoder structure. The encoder extracts multi-scale features, and the decoder gradually restores the spatial resolution.

[0166] The U-Net network employs a symmetrical encoder-decoder structure. The encoder extracts multi-level, multi-scale features from the input image, while the decoder progressively restores the abstract features extracted by the encoder to the original image resolution and performs pixel-level classification. After each downsampling stage in the encoder path, a CBAM attention module is integrated.

[0167] The Convolutional Attention Module (CBAM) consists of a channel attention submodule and a spatial attention submodule, which are expressed as follows:

[0168] Channel attention :

[0169]

[0170] in, Given the input feature map, For global average pooling, For global max pooling, It is a multilayer perceptron. Use the Sigmoid activation function;

[0171] Spatial attention :

[0172]

[0173] in, These are intermediate features after processing by the channel attention module. For channel-dimensional average pooling, Max pooling for the channel dimension. For standard convolutional layers, Use the Sigmoid activation function;

[0174] Finally, the output features after fusion attention are obtained. ;

[0175] S4-2: A hybrid loss function combining binary cross-entropy loss and Dice loss is used to optimize the model and improve its generalization ability.

[0176] During model training, a hybrid loss function combining binary cross-entropy loss and Dice loss is used for model optimization, as detailed below:

[0177] The hybrid loss function Defined as:

[0178]

[0179] Binary cross-entropy loss :

[0180]

[0181] in, The total number of samples, For the first The true labels of each sample For the first The predicted probability that a sample belongs to a crack;

[0182] Dice loss :

[0183]

[0184] in, The total number of samples, For the first The predicted probability that a sample belongs to a crack. For the first The true labels of each sample To ensure numerical stability, a smoothing term is used to avoid denominators of 0.

[0185] The model was trained using the Adam optimizer with an initial learning rate of 1×10⁻⁶. -4 The training batch size was set to 8, and the training was iterated for 200 rounds until the loss function fully converged. During training, data augmentation operations such as random rotation, flipping, and brightness changes were applied to the input image to improve the model's robustness to different imaging conditions.

[0186] S4-3: Use the trained semantic segmentation model to infer the unfolded diagram of the hole wall and obtain the crack probability map. The optimal global threshold for the image is automatically calculated using the maximum inter-class variance method. The probabilistic map is transformed into a binary segmentation map, and the segmentation results are optimized by combining morphological processing.

[0187] The trained semantic segmentation model is used to infer the new pore wall unfolded diagram and output a crack probability map. The Otsu algorithm is used to automatically calculate the optimal global threshold for the image. Thus, the probability graph Convert to binary image

[0188]

[0189] in, The coordinates of the pixels. For the original image in Pixel value at coordinates This is the binarization threshold;

[0190] The following morphological operations are performed sequentially on the initial binary image to optimize segmentation quality. First, a morphological opening operation is performed to remove unwanted details and noise. The opening operation is a process of erosion followed by dilation, and its expression is:

[0191]

[0192] in, The image after the opening operation. For the original image, As a structural element, This represents the coordinate offset within the structuring element;

[0193] Then, the area threshold filtering method is used to calculate the area of ​​each connected component in the binary graph. Set an area threshold For all that satisfy The connected components are identified as residual noise, and their pixel values ​​are set as the background color. This effectively filters out large discrete noise points that were not cleared in the opening operation, providing a high-quality binary image for subsequent skeleton extraction.

[0194] S4-4: The high-quality binary image after morphological processing is refined using the structural surface averaging method to obtain the skeleton image. Then, the initial skeleton is optimized and denoised using the curve slope iteration method to remove outliers in the skeleton image and further enhance the smoothness of the fracture skeleton.

[0195] The high-quality binary image after morphological processing is refined using the structural surface averaging method to extract the fracture skeleton. This method locates the upper and lower edges of non-zero pixel regions column by column and calculates their midpoint coordinates, as shown in the following formula:

[0196]

[0197] in, With the center point, The upper and lower boundaries of the non-zero pixels in this column.

[0198] all Construct the initial skeleton point set ,

[0199] Then, the curve slope iteration method is used to optimize and denoise the initial skeleton. For an ordered skeleton point set... Calculate the slope between adjacent points in sequence. ,when If a point is identified as a noise-induced abrupt change, it is removed from the skeleton point set. This iterative process is repeated until no more slope abrupt changes satisfying the above conditions appear in the skeleton point set, thus obtaining a smooth and continuous final fracture skeleton. .

[0200] S4-5: For each optimized fracture skeleton, perform quantitative extraction of geometric parameters and spatial coordinates, including: dip angle, trace length, fracture width, fracture density, and three-dimensional spatial coordinates.

[0201] For each marked single crack, based on its optimized skeleton and binary region, the following geometric parameters and spatial coordinates are quantitatively extracted:

[0202] inclination :

[0203]

[0204] in, This represents the absolute value of the slope of the fracture skeleton.

[0205] Long :

[0206]

[0207] in, This represents the total number of skeleton pixels in the crack. For the index of the point, These are two adjacent coordinate points;

[0208] Crack width :

[0209]

[0210] in, This represents the number of sampling points along the skeleton. For the index of the point, No. The width of the crack at each sampling point;

[0211] Crack density :

[0212]

[0213] in, To calculate the area, It is the sum of the lengths of all the crack traces;

[0214] By utilizing the geometric mapping relationship between the unfolded diagram of the hole wall and the three-dimensional cylindrical surface, the two-dimensional pixel coordinates on the skeleton are... Inverse calculation to three-dimensional coordinates in the borehole coordinate system Generate a three-dimensional spatial skeleton point set :

[0215]

[0216]

[0217] in, For circumferential angle, Where is the borehole radius, This represents the total width of the panoramic unfolded image. The coordinates are the hole depth coordinates;

[0218] Based on the three-dimensional spatial skeleton point set, the three-dimensional point set of the fracture skeleton is regarded as a structural surface. The least squares method is used for plane fitting, and the plane equation is set as follows: By minimizing Obtain the normal vector ;

[0219] Based on the fitted plane normal vector Calculate the three-dimensional dip angle of the fracture : Based on the projection of the normal vector onto the horizontal plane Calculate the three-dimensional dip of the fracture : ; Three-dimensional trace length : ;

[0220] Final output fracture feature data table , This serves as the basic input for subsequent 3D fracture modeling and lithological fusion.

[0221] Specifically, further, step 5 involves extracting time-domain and frequency-domain features based on the aligned depth-domain parameter data; inputting the feature vectors into the XGBoost ensemble learning model for training and prediction; and outputting lithology classification results along the borehole depth direction, specifically including:

[0222] S5-1: Extract time-domain and frequency-domain features from the aligned drill weight (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling rate (ROP) within a sliding window to construct a high-dimensional feature vector:

[0223] In the temporal features, each parameter is within the depth window. sequence within Statistical characteristics are defined as including: mean. Standard deviation skewness kurtosis coefficient of variation and the correlation coefficient between drilling pressure and torque These factors collectively reflect the average stress state, fluctuation characteristics, and nonlinear mechanical response of the rock strata during the drilling process.

[0224] In the frequency domain, the time-domain signal is transformed to the frequency domain using the Fast Fourier Transform (FFT), and the transformation formula is as follows:

[0225]

[0226] in, For frequency domain signals, For frequency domain variables, The length of the timing signal. The first time domain signal One sample, The imaginary unit;

[0227] Frequency domain characteristics can effectively distinguish the unique mechanical response modes corresponding to different lithologies such as hard rock, soft rock, and fracture zones. These mainly include: dominant frequency. Power spectral energy Spectral center Spectral entropy ;

[0228] For each depth window, the time-domain and frequency-domain features of the four types of parameters are concatenated in a fixed order to form the high-dimensional feature vector of that window:

[0229]

[0230] in, For time-domain eigenvalues, For frequency domain eigenvalues, WOB, RPM, TOR, and ROP correspond to drilling pressure, rotational speed, torque, and mechanical drilling speed, respectively.

[0231] S5-2: The time-domain and frequency-domain features extracted from each drilling parameter are concatenated to construct a high-dimensional feature vector for each depth window. And it is standardized using Z-score:

[0232] The time-domain and frequency-domain features are concatenated to construct a high-dimensional feature vector for each depth window. To avoid the impact of differences in the dimensions of different features on the model convergence, the feature vectors are... Z-score standardization is adopted.

[0233]

[0234] in, and These are the mean and standard deviation of the feature across all training samples, respectively. and These are the normalized eigenvector and the original eigenvector, respectively.

[0235] The standardized feature matrix is ​​obtained by continuously distributing it along the borehole depth direction. ;

[0236] S5-3: Lithology classification is performed using the XGBoost ensemble learning model on the standard feature matrix.

[0237] The XGBoost ensemble learning model was used for lithology classification. This model combines... Decision Tree The prediction is performed, and its output is:

[0238]

[0239] in, For the first The predicted output for each sample, The number of base learners, For the first The base learner for the first... Feature vector of each sample The output,

[0240] The model training objective function consists of a loss function and a regularization term:

[0241]

[0242] in, The total number of samples, For the first The loss term for each sample, For the first The complexity of the regularization term for each basis learner;

[0243] After the model is trained, the probability distribution of each lithology category is output using the Softmax function. :

[0244]

[0245] in, No. Each sample in the input features Below, it belongs to the lithological category. The probability, Model for the first Each sample belongs to category The original prediction score, The sum of the original predicted scores for all lithological categories;

[0246] The final lithological classification result is determined by the category with the highest probability. By performing sliding window prediction on a continuous depth domain, a sequence of lithological profiles along the borehole direction can be obtained. .

[0247] Specifically, such as Figure 6 As shown, step 6 involves deep fusion of the fracture and lithology identification results: feature-level fusion is performed, concatenating the fracture feature vector and the lithology feature vector into a joint feature vector; decision-level fusion is performed, using a random forest model for comprehensive reasoning based on the joint feature vector to generate a refined rock mass quality classification, specifically including:

[0248] Output fracture characteristic data table Perform statistical aggregation along the borehole depth direction, sliding within the same depth window aligned with the drilling parameters, from... Extract all depth coordinates For cracks falling within the current window, data filtering is performed. Based on the filtered set of cracks, the statistical characteristics of the cracks in that window are calculated to form a crack feature vector. , Furthermore, the fracture characteristic data table Transformed into a sequence distributed along the depth. ;

[0249] Lithological characteristic vector Includes time-domain features, frequency-domain features, and XGBoost pre-classification output.

[0250] Fusing fracture feature vectors at the feature level and lithological characteristic vectors Construct joint feature vector ;

[0251] At the decision level, a random forest model is constructed to solve the complex multi-feature, nonlinear classification problem of rock mass quality grading. After model training, the joint feature vector is used... The input is the rock mass quality grade, and the output is a multi-class probability vector. The final rock mass quality grade is determined by the category with the highest probability. Among them, rock mass quality grade Rocks are classified into five categories: Category I: hard rock; Category II: relatively hard rock; Category III: medium-hard rock; Category IV: relatively soft rock; Category V: soft soil.

[0252] The fusion and identification results of all depth points are integrated to generate a continuous profile along the borehole depth direction. This provides precise borehole constraints for subsequent three-dimensional geological modeling.

[0253] Specifically, such as Figure 7-8 As shown, step 7, which integrates the fusion interpretation data from all boreholes, employs implicit modeling techniques and a discrete fracture network generation algorithm to construct an integrated model that simultaneously includes lithological entity distribution and a three-dimensional fracture network, and then provides a three-dimensional visualization. Specifically, this includes:

[0254] The fused identification results at each depth location are integrated with the extracted three-dimensional coordinate information of the fracture and uniformly packaged into a standardized borehole data unit: ,in Including three-dimensional tilt angle Three-dimensional tendency Three-dimensional trace length and the set of three-dimensional coordinates of the skeleton points ;

[0255] The fused interpretation results from multiple boreholes are uniformly registered according to their spatial coordinates to form a three-dimensional structured dataset: ;

[0256] A three-dimensional lithological solid model is generated using an implicit modeling method, for each lithological category. Constructing implicit functions:

[0257]

[0258] in, For coordinate variables in space, For lithological categories The number of borehole constraint points, For the first The weight of each control point For radial basis functions RBF, No. The three-dimensional coordinates of each control point;

[0259] Establish a regular 3D voxel mesh in the modeling space. For any voxel... Its final lithological category Determined by the following formula:

[0260]

[0261] in, For a voxel position in three-dimensional space, An index for lithological categories, This represents the total number of lithological categories. No. Lithology in location The implicit function value at that location;

[0262] In two types of lithology and Construct a smooth interface between them and define its isosurface function. And solve for its zero isosurface:

[0263]

[0264] in, A point in three-dimensional space. No. Lithology at point The implicit function value at that location. No. Lithology at point The implicit function value at that location;

[0265] This isosurface represents the boundary between the two lithologies. The moving cube algorithm is used to extract the isosurface from the scalar field. Within each voxel, this algorithm calculates the vertex positions of the triangular facets using linear interpolation.

[0266]

[0267] in, Let the spatial coordinate vector of the first sampling point be... The spatial coordinate vector of the second sampling point. isosurface function At point The function value at that point, isosurface function At point The function value at that point, The constant value of the target isosurface;

[0268] By traversing all voxels, the calculated set of triangular facets is integrated and optimized to form a continuous and smooth three-dimensional lithological entity distribution model. ;

[0269] Deterministic 3D reconstruction of fractures is performed for each fracture identified and extracted from borehole camera images. Using its quantized geometric parameters, a precise geometric reconstruction can be performed on it in three-dimensional space:

[0270] Based on the tendency of the crack With tilt angle Calculate and determine the unit normal vector of the crack plane. :

[0271]

[0272] Using the coordinates of the fracture center point as the base point, the normal vector To determine the direction, the spatial orientation of the fracture is established. Then, a disk model is used to geometrically simplify the fracture, and its radius is determined. From three-dimensional trace length The disk crack surface can be represented as a set of points satisfying the following conditions:

[0273]

[0274] in, For any point in three-dimensional space, For the first Coordinates of the center point of the cracked surface of the disk For the first The unit normal vector of the cracked surface of a disk;

[0275] Add all reconstructed deterministic fracture patches to the deterministic fracture set.

[0276] Statistical simulation of stochastic fractures was conducted. Based on the dip and dip angle data of all deterministic fractures, the normal vector was fitted using a Fisher distribution, and its probability density function was:

[0277]

[0278] in, is a unit direction vector in three-dimensional space. The average direction vector of the distribution. For concentration parameters, It is a hyperbolic sine function. is the normalization constant of the distribution;

[0279] Based on the three-dimensional trace lengths of all deterministic cracks, a probability density function that follows a log-normal distribution is fitted:

[0280]

[0281] in, Let be a random variable representing the length of the crack trace. Mean parameter on a logarithmic scale The standard deviation parameter represents the logarithmic scale. is the normalization constant of the probability density function;

[0282] Statistically analyze its probability distribution model, and based on the statistically obtained distribution model and fracture density field A large number of simulated fractures are randomly generated using the Monte Carlo method, forming a set of random fractures. ;

[0283] Deterministic fracture set and random fracture set The two components are merged to form a complete three-dimensional discrete fracture network model. ;

[0284] To achieve precise geometric fusion, a three-dimensional discrete fracture network model is used. As a cutting tool, it is used to analyze the distribution model of lithological entities. Perform three-dimensional Boolean operations:

[0285]

[0286] Geometrically equivalent to using the fracture surface to cut a complete rock mass, generating a collection of independent rock blocks divided by the fracture surface, thus reproducing the spatial cutting relationship of the fracture to the rock mass;

[0287] Based on geometric fusion, each geometric element in the model is assigned corresponding geological attributes to construct a unified three-dimensional attribute field. :

[0288]

[0289] in, Lithological category Rock mass quality grade It is three-dimensional inclination. For three-dimensional tilt angle, For the fracture density field;

[0290] This joint attribute model realizes a unified digital expression of lithology, structural planes, and rock mass quality parameters in three-dimensional space;

[0291] Finally, the lithological entity distribution model With three-dimensional discrete fracture network model Deep integration is achieved to realize a unified three-dimensional spatial expression of lithology, structural planes, and rock mass quality parameters, and to provide a visual representation.

[0292] like Figure 9 As shown, the present invention also discloses a tunnel three-dimensional geological body modeling system based on borehole photography and drilling lithology perception, including: a multi-source data acquisition and preprocessing module, a depth domain data alignment and fusion module, a fracture intelligent identification and parameter extraction module, a lithology intelligent perception and classification module, a multi-source information fusion and rock mass quality grading module, and a three-dimensional geological body integrated modeling and visualization module.

[0293] Multi-source data acquisition and preprocessing module: used for the synchronous acquisition and preprocessing of in-hole camera images and drilling parameters. It integrates image acquisition and drilling parameter acquisition units, can control the probe to move at a constant speed and record images and depth synchronously, and simultaneously acquire drilling pressure, rotation speed, torque and mechanical drilling speed in real time.

[0294] The depth domain data alignment and fusion module is used to achieve accurate alignment of in-hole camera images and drilling parameters in the depth domain. It establishes a depth-time mapping model based on the time integral of the mechanical drilling rate, transforms the drilling parameters from the time domain to the depth domain, and resamples them to the same depth coordinates as the image through a linear interpolation algorithm. Finally, it outputs a depth-aligned multi-source data matrix.

[0295] Intelligent fracture identification and parameter extraction module: used for automatic fracture identification and parameter extraction. Based on the improved U-Net semantic segmentation network, the module automatically extracts the dip angle, trace length, width, density, and three-dimensional spatial coordinates of the segmented fractures through algorithms such as skeletonization and connected component analysis, and outputs a fracture feature data table;

[0296] Lithology Intelligent Perception and Classification Module: This module is used to achieve automatic lithology classification. It extracts time-domain and frequency-domain features from aligned drilling parameters to construct feature vectors, and uses a trained XGBoost model for prediction. Finally, it outputs a lithology classification profile along the borehole depth direction.

[0297] Multi-source information fusion and rock mass quality classification module: used for comprehensive analysis of fracture and lithology information, concatenating fracture feature vectors and lithology feature vectors, inputting them into a random forest model for inference, and finally outputting a rock mass quality classification profile (Class I to V) along the borehole depth direction.

[0298] The integrated 3D geological body modeling and visualization module employs implicit modeling techniques to generate lithological entity models and utilizes a discrete fracture network algorithm to generate deterministic and stochastic fractures. Geometric fusion of the rock mass and fracture network is achieved through 3D Boolean operations, ultimately generating an interactive 3D visualized geological model.

[0299] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling-while-drilling lithology sensing, characterized in that, Specifically, the steps include the following: Step S1: Use an in-hole camera to acquire a sequence of optical images of the borehole wall and perform image preprocessing to generate a continuous unfolded image of the borehole wall; Step S2: Use the measurement while drilling system to synchronously collect engineering parameters during the drilling process and perform data preprocessing; Step S3: Establish a depth-time mapping relationship, resample the preprocessed drilling parameters from the time domain to the depth domain, so that they correspond one-to-one with the borehole wall unfolding diagram in the depth domain, and obtain depth-aligned multi-source data; Step S4: Perform semantic segmentation of the hole wall unfolded diagram, and calculate the geometric parameters and three-dimensional spatial coordinates of each crack based on the segmentation results to form a crack feature dataset; Step S5: Based on the depth-aligned multi-source data, extract time-domain and frequency-domain features, input them into a machine learning model for lithology classification, and output the lithology classification results along the borehole depth direction; Step S6: The fracture feature dataset and the lithology classification results are fused at the feature level and the decision level to generate a rock mass quality grade along the borehole depth direction; Step S7: Based on the fracture feature dataset, the lithology classification results, and the rock mass quality grading, an integrated model that includes both lithological entity distribution and a three-dimensional fracture network is constructed using implicit modeling technology and a discrete fracture network generation algorithm, and then visualized in three dimensions.

2. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology sensing according to claim 1, characterized in that, The specific method for acquiring a sequence of optical images of the borehole wall using an in-hole camera device and performing image preprocessing to generate a continuous borehole wall unfolded image in step S1 includes: S1-1: A high-resolution borehole imaging system is used. The probe is placed into the borehole at a constant speed of 2m / min and then pulled out at the same speed. During the pulling process, the probe continuously acquires a 360° annular image sequence of the borehole wall at fixed depth intervals. The acquisition software synchronously records the depth coordinates corresponding to each frame of the image. S1-2: Convert the acquired raw image from the RGB color space to the HSV color space, keeping the chroma (H) and saturation (S) components unchanged, and only perform multi-scale Retinex processing on the luminance component (V). The RGB model can be converted to the HSV model using the following formula: ； ； ； Where V is the luminance component, S is the saturation, H is the chroma, and R, G, B are the color space component values. The chromaticity (H) and saturation (S) components are retained, while only the luminance (V) component undergoes MSR enhancement to reduce luminance unevenness and reflective interference. ； in, pixel coordinates The reflection component at that location, For scale quantity These are the weighting coefficients. For the original input image at pixel points The brightness value at that location, For the first Gaussian wrapping function of scale; For the enhanced Linear contrast stretching is performed, and the result is combined with the preserved chroma (H) and saturation (S) components, and finally converted back to the RGB color space to obtain an enhanced image with uniform illumination. S1-3: A bilateral filtering algorithm is used to eliminate residual noise in the image, preserve the edge structure information of the crack, and further improve the image quality; The filter calculation formula is as follows: ； in, For target pixel The output grayscale value, To normalize the weights, For pixels The scope of the field For pixels in the neighborhood The original grayscale value, For spatial domain weight kernel, Kernel weights for the range; S1-4: Between adjacent image frames, the SIFT algorithm is used for key point detection, and the kd-tree nearest neighbor search algorithm is used to quickly match the feature descriptors of the two frames. The SIFT algorithm is used to detect scale- and rotation-invariant keypoints, and a 128-dimensional feature description vector is calculated for each keypoint. A kd-tree nearest neighbor search algorithm is used for fast feature vector matching. To eliminate false matches, Lowe's ratio test is employed to retain matching pairs where the ratio of the nearest neighbor to the second nearest neighbor distance is less than 0.8, thereby obtaining a high-precision set of matching points. ; S1-5: Based on the collocation set, use the RANSAC algorithm to estimate the homography matrix H between the two images; Based on the matching point set The RANSAC algorithm is used to estimate the homography matrix H between two frames of images; ； ； Where H is the homography matrix, To find the matrix that minimizes the subsequent objective function, The origin point, for Corresponding to the transformed points; S1-6: Using the homography matrix H, all images are projected onto a unified planar coordinate system. For overlapping areas, a multi-band mixing algorithm is used for fusion to preserve the image details to the maximum extent and generate an initial panoramic unfolded image. S1-7: Perform contrast-limited adaptive histogram equalization (CLAHE) and adaptive gamma correction on the initial panoramic unfolded image to enhance image details, and perform standardized cropping on the image.

3. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology perception according to claim 1, characterized in that, The specific method for synchronously acquiring engineering parameters during the drilling process and performing data preprocessing using the measurement-while-drilling system described in step S2 includes: S2-1: Engineering parameters during the drilling process are synchronously acquired using a measurement-while-drilling (MWD) system, including drill pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP). The acquisition frequency of these WWD parameters is preferably no less than 1 Hz, and the sampling interval is preferably 1 second. The sampling frequency can be adjusted according to the probe retraction rate. With the desired depth resolution According to the relation The calculation determines the frequency, and it can be adjusted within the range of 1-10Hz to adapt to different drilling speeds and resolution requirements; S2-2: Adopt The standard deviation criterion is used to detect outliers in the data, and linear interpolation is used to repair the data, eliminating extreme anomalies caused by instrument vibration, signal loss, or sudden geological changes during drilling. use Criterion: Calculate the mean of the parameter sequence. and standard deviation For each data point in the sequence Its standardized residual is ,like Then determine This is an outlier; For data points marked as outliers, linear interpolation is used to repair them to ensure data continuity. The formula is as follows: ； in, and It is the nearest normal data point before and after the outlier; S2-3: Apply the SG filter to the time series of each parameter to smooth the signal fluctuation trend while suppressing high-frequency noise; To preserve the physical trend of the drilling signal and suppress high-frequency noise, SG filtering is applied to the time series data, for a window width of... The order of the polynomial is The SG filter, which is located at any position within the window. The fitted value is expressed as: ； in, This refers to the local coordinate index within the current sliding window. For index The value obtained by polynomial fitting at that point. The coefficients to be determined are obtained by minimizing the fitting error of all data points within the window. To determine; The smooth output of the window center point is: ； in, For global location index, For local coordinate indices within the sliding window, For the window radius, These are the filter weight coefficients. For global position The original parameter values ​​collected at the location; S2-4: To maintain the temporal synchronization and physical correlation among the smoothed parameters, the same SG filtering parameters are used to perform multi-parameter synchronous smoothing on the four sequences of drilling pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP), constructing a synchronous smoothing matrix that includes time coordinates and the physical values ​​of the smoothed parameters. ; S2-5: Eliminate the differences in dimensions and numerical ranges between different parameters, and perform Z-score standardization on the smoothed data to obtain the standardized drilling parameter matrix; Due to the significant differences in the physical dimensions and numerical ranges of the drilling parameters, Z-score normalization was used to smooth the matrix to ensure numerical comparability in subsequent modeling and feature extraction. Normalize each parameter column: ； in, The standardized dimensionless value. For parameters at time points The smoothed physical value, This is the sample mean of the parameter sequence. This is the sample standard deviation of the parameter sequence; After standardization, a standardized parameter matrix is ​​formed. : 。 4. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology sensing according to claim 1, characterized in that, The specific method for establishing the depth-time mapping relationship in step S3, which involves resampling the preprocessed drilling parameters from the time domain to the depth domain so that they correspond one-to-one with the borehole wall unfolding diagram in the depth domain, to obtain depth-aligned drilling parameter data, includes: Based on the physical principles of drilling, the drilling depth is the integral of the mechanical drilling rate of power (ROP) with respect to time; therefore, a mapping model from the time domain to the depth domain is established: ； in, From start to time Total depth drilled, For mechanical rotation speed, For integration variables; In discrete data processing, the trapezoidal numerical integration method is used for solving; for the th Each time point, and its corresponding depth value The calculation formula is as follows: ； in, For the first Drilling depth at each time step For the first Drilling depth at each time step For the first Mechanical drilling speed per time step For the first Mechanical drilling speed per time step The time interval between two adjacent time steps; Then, linear interpolation is used to resample all parameters at equal intervals: ； in, For target depth Interpolation results of parameters at the location, For the first time to take pictures inside the hole A depth coordinate point, , Known adjacent depth points (satisfying) ), These correspond to the depth parameter values ​​respectively; Finally, precise fusion and alignment of "visual information" and "physical information" in the depth domain were achieved, outputting a multi-source data matrix of depth and alignment. .

5. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology sensing according to claim 1, characterized in that, The specific method for performing crack semantic segmentation on the unfolded diagram of the borehole wall in step S4, and calculating the geometric parameters and three-dimensional spatial coordinates of each crack based on the segmentation results to form a crack feature dataset includes: S4-1: Introduce the Convolutional Attention Module (CBAM) that combines channel attention and spatial attention mechanisms to construct a U-Net convolutional neural network with an encoder-decoder structure. The encoder part extracts multi-scale features, and the decoder gradually restores the spatial resolution. The Convolutional Attention Module (CBAM) consists of a channel attention submodule and a spatial attention submodule, which are expressed as follows: Channel attention : ； in, Given the input feature map, For global average pooling, For global max pooling, It is a multilayer perceptron. Use the Sigmoid activation function; Spatial attention : ； in, These are intermediate features after processing by the channel attention module. For channel-dimensional average pooling, Max pooling for the channel dimension. For standard convolutional layers, Use the Sigmoid activation function; Finally, the output features after fusion attention are obtained. ; S4-2: A hybrid loss function combining binary cross-entropy loss and Dice loss is used to optimize the model and improve its generalization ability. During model training, a hybrid loss function combining binary cross-entropy loss and Dice loss is used for model optimization, as detailed below: The hybrid loss function Defined as: ； Binary cross-entropy loss : ； in, The total number of samples, For the first The true labels of each sample For the first The predicted probability that a sample belongs to a crack; Dice loss : ； in, The total number of samples, For the first The predicted probability that a sample belongs to a crack. For the first indivual Sample true labels To ensure numerical stability, a smoothing term is used to avoid denominators of 0. The model was trained using the Adam optimizer with an initial learning rate of 1×10⁻⁶. -4 The training batch size is set to 8, and the training is iterated for 200 rounds until the loss function fully converges. During the training process, data augmentation operations such as random rotation, flipping, and brightness changes are applied to the input image to improve the robustness of the model to different imaging conditions. S4-3: Use the trained semantic segmentation model to infer the unfolded diagram of the hole wall and obtain the crack probability map. The optimal global threshold for the image is automatically calculated using the maximum inter-class variance method. The probabilistic map is transformed into a binary segmentation map, and the segmentation results are optimized by combining morphological processing. The trained semantic segmentation model is used to infer the new pore wall unfolded diagram and output a crack probability map. The Otsu algorithm is used to automatically calculate the optimal global threshold for the image. Thus, the probability graph Convert to binary image ; ； in, The coordinates of the pixels. For the original image in Pixel value at coordinates This is the binarization threshold; The following morphological operations are performed sequentially on the initial binary image to optimize segmentation quality: first, morphological opening is performed to remove unwanted details and noise; then, area thresholding is used to calculate the area of ​​each connected component in the binary image. Set an area threshold For all that satisfy The connected components are identified as residual noise, and their pixel values ​​are set as the background color. This effectively filters out large discrete noise points that were not cleared in the opening operation, providing a high-quality binary image for subsequent skeleton extraction. S4-4: The high-quality binary image after morphological processing is refined using the structural surface averaging method to obtain the skeleton image. Then, the initial skeleton is optimized and denoised using the curve slope iteration method to remove abnormal points in the skeleton image and further enhance the smoothness of the fracture skeleton. For each optimized fracture skeleton, quantitative extraction of geometric parameters and spatial coordinates is performed, including: dip angle, trace length, fracture width, fracture density, and three-dimensional spatial coordinates; The high-quality binary image after morphological processing is thinned using the structural surface averaging method to extract the fracture skeleton. This method generates a single-pixel-width skeleton that maintains the topological connectivity of the fracture by locating the upper and lower edges of non-zero pixel regions column by column and calculating their midpoint coordinates. Then, the curve slope iteration method is used to optimize and denoise the initial skeleton. For an ordered skeleton point set... Calculate the slope between adjacent points in sequence. ,when If the point is determined to be a sudden change point caused by noise, it is removed from the skeleton point set. The above iterative process is repeated until no slope sudden change points satisfying the above conditions appear in the skeleton point set, thereby obtaining a smooth and continuous final fracture skeleton. For each marked single crack, based on its optimized skeleton and binary region, the following geometric parameters and spatial coordinates are quantitatively extracted: inclination : ； in, This represents the absolute value of the slope of the fracture skeleton. Long : ； in, This represents the total number of skeleton pixels in the crack. For the index of the point, These are two adjacent coordinate points; Crack width : ； in, This represents the number of sampling points along the skeleton. For the index of the point, No. The width of the crack at each sampling point; Crack density : ； in, To calculate the area, It is the sum of the lengths of all the crack traces; By utilizing the geometric mapping relationship between the unfolded diagram of the hole wall and the three-dimensional cylindrical surface, the two-dimensional pixel coordinates on the skeleton are... Inverse calculation to three-dimensional coordinates in the borehole coordinate system Generate a three-dimensional spatial skeleton point set : ； ； in, For circumferential angle, Where is the borehole radius, This represents the total width of the panoramic unfolded image. The coordinates are the hole depth coordinates; Based on the three-dimensional spatial skeleton point set, the three-dimensional point set of the fracture skeleton is regarded as a structural surface. The least squares method is used for plane fitting, and the plane equation is set as follows: By minimizing Obtain the normal vector ; Based on the fitted plane normal vector Calculate the three-dimensional dip angle of the fracture : Based on the projection of the normal vector onto the horizontal plane Calculate the three-dimensional dip of the fracture : ; Three-dimensional trace length : ; Final output fracture feature data table , This serves as the basic input for subsequent 3D fracture modeling and lithological fusion.

6. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology perception according to claim 1, characterized in that, The specific method for extracting time-domain and frequency-domain features from the depth-aligned drilling parameter data in step S5, inputting them into a machine learning model for lithology classification, and outputting lithology classification results along the borehole depth direction includes: S5-1: Extract time-domain and frequency-domain features from the aligned drill pressure (WOB), rotational speed (RPM), torque (TOR), and mechanical drilling speed (ROP) within a sliding window to construct a high-dimensional feature vector; In the temporal features, each parameter is within the depth window. sequence within The defined statistical characteristics include: mean, standard deviation, skewness, kurtosis, coefficient of variation, and correlation coefficient between drill pressure and torque, which together reflect the average stress state, fluctuation characteristics, and nonlinear mechanical response of the rock formation during the drilling process. In the frequency domain, the time-domain signal is converted to the frequency domain using the Fast Fourier Transform (FFT), and the conversion formula is as follows: ； in, For frequency domain signals, For frequency domain variables, The length of the timing signal. The first time domain signal One sample, The imaginary unit; Frequency domain characteristics can effectively distinguish the unique mechanical response modes corresponding to different lithologies such as hard rock, soft rock, and fracture zones. These mainly include: dominant frequency, power spectral energy, spectral centroid, and spectral entropy. S5-2: The time-domain and frequency-domain features extracted from each drilling parameter are concatenated to construct a high-dimensional feature vector for each depth window. And it is standardized using Z-score; The time-domain and frequency-domain features are concatenated to construct a high-dimensional feature vector for each depth window. To avoid the impact of differences in the dimensions of different features on the model convergence, the feature vectors are... Z-score standardization is adopted. ； in, and These are the mean and standard deviation of the feature across all training samples, respectively. and The normalized eigenvectors and the original eigenvectors are used to obtain the standardized feature matrix that is continuously distributed along the borehole depth direction. ; S5-3: Lithology classification is performed using the XGBoost ensemble learning model on the standard feature matrix; Lithology classification is performed using the XGBoost ensemble learning model; this model combines... Decision Tree The prediction is performed, and its output is: ； in, For the first The predicted output for each sample, The number of base learners, For the first The base learner for the first... Feature vector of each sample The output; The model training objective function consists of a loss function and a regularization term: ； in, The total number of samples, For the first The loss term for each sample, For the first The complexity of the regularization term for each basis learner; After the model is trained, the probability distribution of each lithology category is output using the Softmax function. The final lithological classification result is determined by the category with the highest probability. By performing sliding window prediction on a continuous depth domain, a sequence of lithological profiles along the borehole direction can be obtained. .

7. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology perception according to claim 1, characterized in that, The specific method for fusing the fracture feature dataset with the lithology classification results at the feature level and decision level in step S6 to generate rock mass quality grading along the borehole depth direction includes: Output fracture characteristic data table Perform statistical aggregation along the borehole depth direction, sliding within the same depth window aligned with the drilling parameters, from... Extract all depth coordinates For cracks falling within the current window, data filtering is performed. Based on the filtered set of cracks, the statistical characteristics of the cracks in that window are calculated to form a crack feature vector. Furthermore, the fracture characteristic data table Transformed into a sequence distributed along the depth. ; Fusing fracture feature vectors at the feature level and lithological characteristic vectors Construct joint feature vector ; At the decision level, a random forest model is constructed to solve the complex multi-feature, nonlinear classification problem of rock mass quality grading, and the output is a multi-class probability vector of rock mass quality grades. The final rock mass quality grade is determined by the category with the highest probability. The rock mass quality is classified into five grades: Class I rock: hard rock; Class II rock: relatively hard rock; Class III rock: medium hard rock; Class IV rock: relatively soft rock; Class V rock: soft soil; The fusion and identification results of all depth points are integrated to generate a continuous profile along the borehole depth direction, providing accurate borehole constraints for subsequent 3D geological modeling.

8. The method for modeling three-dimensional geological bodies of tunnels based on in-hole photography and drilling lithology perception according to claim 1, characterized in that, The specific method for constructing an integrated model that simultaneously includes lithological entity distribution and a three-dimensional fracture network based on the fracture feature dataset, the lithology classification results, and the rock mass quality grading in step S7, and for performing three-dimensional visualization, includes: The fused interpretation results from multiple boreholes are uniformly registered according to their spatial coordinates to form a three-dimensional structured dataset: ; A three-dimensional lithological solid model is generated using an implicit modeling method, for each lithological category. Constructing implicit functions: ； in, For coordinate variables in space, For lithological categories The number of borehole constraint points, For the first The weight of each control point For radial basis functions RBF, No. The three-dimensional coordinates of each control point; By implicit function voxelization and isosurface reconstruction, a continuous and smooth three-dimensional lithological entity distribution model is formed. ; A deterministic fracture 3D reconstruction is performed. Using the 3D dip angle, 3D dip direction, 3D spatial trace length, and center point coordinates, the fracture is reconstructed as a deterministic fracture surface in 3D space and added to the deterministic fracture set. ; Statistical simulations of stochastic fractures were conducted. Based on the attitude and trace length data of all deterministic fractures, their probability distribution model was statistically analyzed. According to the statistically obtained distribution model and fracture density field, a large number of simulated fractures were randomly generated using the Monte Carlo method, forming a set of stochastic fractures. ; Deterministic fracture set and random fracture set The two components are merged to form a complete three-dimensional discrete fracture network model. ; Finally, the lithological entity distribution model With three-dimensional discrete fracture network model Deep integration is achieved to realize a unified three-dimensional spatial expression of lithology, structural planes, and rock mass quality parameters, and to provide a visual representation.

9. A system used in the tunnel three-dimensional geological body modeling method based on borehole photography and drilling lithology perception as described in claim 1, characterized in that, include: The multi-source data acquisition and preprocessing module is used to simultaneously acquire and preprocess in-hole camera images and drilling parameters; The depth domain data alignment and fusion module is used to achieve accurate alignment of in-hole camera images and drilling parameters in the depth domain; The intelligent crack identification and parameter extraction module is used to automatically identify cracks and extract their parameters based on the borehole wall unfolding diagram. The intelligent lithology sensing and classification module is used for automatic lithology classification based on drilling parameters. The multi-source information fusion and rock mass quality classification module is used to fuse fracture and lithological information and perform rock mass quality classification. The integrated 3D geological modeling and visualization module is used to construct and visualize an integrated 3D geological model that incorporates lithology and fracture networks.

Images