Three-dimensional resistivity-based comprehensive survey method for tunnel engineering in limestone interlayer section

Abstract

The present application provides a comprehensive investigation method for tunnel engineering in limestone interlayer section based on three-dimensional resistivity, relates to the technical field of engineering geological investigation, and comprises the following steps: S1: collecting regional geological data of a tunnel site area to obtain a karst development region; S2: establishing a resistivity value database of different rock-soil bodies; S3: designing a survey line arrangement scheme according to the tunnel orientation of the tunnel site area, testing the two-dimensional profile resistivity of the tunnel site area based on the survey line arrangement scheme, and obtaining two-dimensional profile resistivity data; S4: constructing a three-dimensional geological model, extracting a three-dimensional abnormal body and a soluble rock three-dimensional geological body according to the three-dimensional geological model and the karst development region; S5: drilling verification is performed on the three-dimensional abnormal body and the soluble rock three-dimensional geological body, the karst development region is corrected according to the verification result, and the three-dimensional geological model is updated; and S6: the updated three-dimensional geological model is combined to guide engineering geological line selection. The present application provides a precise basis for tunnel engineering geological line selection and site construction.

Description

Technical Field

[0001] The present invention relates to the technical field of engineering geological survey, and in particular to a comprehensive survey method for tunnel engineering in limestone interlayer sections based on three-dimensional resistivity. Background Art

[0002] Limestone areas have complex geological and geomorphological conditions and are susceptible to factors such as geological structure and groundwater. This creates significant potential risks for railway construction in limestone areas. The occurrence and development of railway construction risks in limestone areas are particularly unique, especially when tunneling through limestone interlayers. Limited by the limitations of surveying and drilling techniques, drilling alone is difficult to determine the spatial location and extent of karst development within these interlayers. Consequently, the inability to accurately ascertain geological conditions makes it difficult to proactively avoid these interlayers or develop targeted protective measures, leading to significant engineering risks. Currently, the exploration of limestone interlayers and their karst development in tunneling projects typically relies on drilling combined with two-dimensional geophysical profile analysis. However, this approach struggles to accurately determine the spatial distribution of interlayers and the accuracy of karst development is poor. Inaccuracies in this analysis can easily lead to ineffective design measures, potentially causing engineering accidents and impacting construction. Therefore, developing a comprehensive survey method with three-dimensional properties that can accurately reflect the patterns of limestone interlayers and reveal the extent of karst development is an urgent need. Summary of the Invention

[0003] In view of this, the present invention proposes a comprehensive survey method for tunnel engineering in limestone interlayers based on three-dimensional resistivity. In the survey work of limestone interlayers in tunnel engineering, it solves the shortcomings of poor accuracy of conventional drilling combined with two-dimensional geophysical profile analysis of karst spatial distribution and development degree, and provides an accurate basis for geological line selection and on-site construction of tunnel engineering.

[0004] The technical solution of the present invention is achieved as follows:

[0005] The present invention provides a comprehensive survey method for tunnel engineering in limestone interlayer sections based on three-dimensional resistivity, comprising the following steps:

[0006] S1: Collect regional geological data of the tunnel site, analyze and preliminarily determine the regional geological data, and determine the karst development area;

[0007] S2: Select typical areas in the tunnel site to conduct high-density resistivity testing and establish a database of resistivity values ​​for different rock and soil bodies;

[0008] S3: Design a survey line layout plan based on the tunnel direction of the tunnel site, measure the two-dimensional profile resistivity of the tunnel site based on the survey line layout plan, and obtain two-dimensional profile resistivity data;

[0009] S4: Preprocess the 2D profile resistivity data based on the resistivity database of different rock and soil bodies, and construct a 3D geological model. Extract 3D anomalies and soluble rock 3D geological bodies based on the 3D geological model and karst development areas;

[0010] S5: Drill and verify the 3D anomaly bodies and soluble rock 3D geological bodies, revise the karst development area based on the verification results, and update the 3D geological model;

[0011] S6: Combine the updated 3D geological model to guide engineering geological line selection.

[0012] Based on the above technical solution, preferably, step S1 includes:

[0013] S11: Collect regional geological data of the tunnel site, including regional geological maps, geological reports, drilling data and auxiliary information;

[0014] S12: Develop a survey plan based on regional geological data and conduct on-site geological surveys based on the survey plan, including field investigations, topographic and geomorphological surveys, and establishing an image database;

[0015] S13: Organize and standardize the drilling data and draw a drilling histogram. Based on the drilling histogram, analyze the characteristics and distribution of the limestone layer and preliminarily determine the degree of karst development.

[0016] S14: Establish a GIS project based on the regional geological map and drill hole data, outline the distribution range of the limestone layer in the GIS layer, and use the interpolation method to estimate the thickness distribution of the limestone layer;

[0017] S15: Based on topographic features, drilling data and field geological survey results, mark potential karst development areas and use different symbols or colors to indicate the degree of karst development;

[0018] S16: According to step S14 and step S15, a limestone distribution range map is generated, and a karst development area distribution map is produced.

[0019] Based on the above technical solution, preferably, step S2 includes:

[0020] S21: Select a typical area within the tunnel site based on regional geological data;

[0021] S22: Sampling of rock and soil in typical areas, where the sampling points cover all rock and soil types within the tunnel site;

[0022] S23: Conduct high-density resistivity testing on the collected samples, calculate the resistivity values ​​of various rock and soil bodies, and draw box plots to show the resistivity distribution characteristics;

[0023] S24: Obtain the resistivity range reference standard of different rock-soil bodies, determine the resistivity range of various rock-soil bodies according to the test result and the reference standard, and establish the resistivity value database of different rock-soil bodies.

[0024] Based on the technical scheme, preferably, the step S3 comprises:

[0025] S31: Design a survey line arrangement scheme, the geophysical survey lines are parallel to the tunnel axis direction, the interval between the geophysical survey lines is D, n geophysical survey lines are arranged according to the interval D, and the survey line arrangement scheme covers the karst development area;

[0026] S32: Position the geophysical survey line by using GPS or a total station, arrange the electrodes at the designed geophysical survey line positions according to the interval D, and record the coordinates of each electrode;

[0027] S33: Perform multi-round resistivity tests by using a high-density resistivity meter, and take the average value of the multi-round test results as the resistivity data;

[0028] S34: Form the profile data of the tunnel site area based on the coordinates of the electrodes, and generate two-dimensional profile resistivity data by combining the resistivity data and the profile data.

[0029] Based on the technical scheme, preferably, the interval D is 10-20 m.

[0030] Based on the technical scheme, preferably, in the resistivity test process, if the resistivity value of the current round appears to be abnormal, the abnormal section is positioned, a geophysical survey line is added, and the interval D between the geophysical survey lines of the section is set to 3-5 m.

[0031] Based on the technical scheme, preferably, the step S4 comprises:

[0032] S41: Import the survey line terrain data, correct the two-dimensional profile resistivity data according to a terrain correction algorithm, and remove noise by using a filtering algorithm to obtain the two-dimensional profile resistivity data after preliminary processing;

[0033] S42: Perform inversion calculation on the two-dimensional profile resistivity data after preliminary processing by using an inversion algorithm, and output the two-dimensional profile resistivity data after inversion;

[0034] S43: Compare the two-dimensional profile resistivity data after inversion with the resistivity value database of different rock-soil bodies, quantify the difference by using a statistical method, and mark the regions with significant differences;

[0035] S44: Correct the regions with significant differences, and return to the step S42 until the inversion result is consistent with the resistivity value database, and output the two-dimensional profile resistivity data after processing;

[0036] S45: Input the processed two-dimensional profile resistivity data into MATLAB and create a structure array, where each structure corresponds to a two-dimensional profile;

[0037] S46: Determine the grid range in the X, Y, and Z directions according to the range of the structure, and set the grid resolution. Merge all processed two-dimensional profile resistivity data into a three-dimensional point set. Each point in the three-dimensional point set includes its X, Y, and Z coordinates and the corresponding resistivity value. Use the interpolation algorithm to perform linear interpolation on the three-dimensional point set to obtain three-dimensional voxel data.

[0038] S47: setting a resistivity threshold of the anomaly, extracting an isosurface from the three-dimensional voxel data according to the resistivity threshold of the anomaly, and taking the area surrounded by the isosurface as a preliminary anomaly, wherein the isosurface is a collection of points with the same resistivity value; performing connectivity analysis on the preliminary anomaly, classifying the connected areas as one anomaly, and obtaining a three-dimensional anomaly;

[0039] S48: Set the resistivity range of the limestone body, preliminarily extract the limestone body according to the isosurface method, optimize the extraction results using the 3D morphology method, and perform continuity analysis on the optimized extraction results to identify the limestone interlayer, i.e., the 3D soluble rock geological body.

[0040] Based on the above technical solution, preferably, in step S47, the connectivity analysis process is:

[0041] S471: Convert the 3D voxel data into binary data, assign the data points on the preliminary abnormal volume a value of 1, and assign the other data points in the 3D voxel data a value of 0;

[0042] S472: Create a marker matrix of the same size as the 3D voxel data, with all initial values ​​set to 0; set an empty list list to store the voxels to be processed;

[0043] S473: For each unlabeled voxel in the preliminary abnormal body, assign it a new unique label and add the voxel to the list;

[0044] S474: When the list is not empty, take a voxel from the list, check its adjacent voxels, and assign the same label to each adjacent voxel that is located in the preliminary abnormal volume and is not labeled, and add it to the list;

[0045] S475: Repeat step S474 until the list is empty;

[0046] S476: Repeat steps S473-S475 until all voxels in the preliminary abnormal volume are labeled, and output a labeling matrix in which the value of each voxel is the label of the preliminary abnormal volume;

[0047] S477: Voxels with the same label are grouped into a 3D anomaly volume to complete connectivity analysis.

[0048] Based on the above technical solution, preferably, step S5 includes:

[0049] S51: Analyze the spatial distribution of 3D anomalies and soluble rock 3D geological bodies, and select representative locations, including areas with obvious resistivity anomalies, soluble rock distribution areas, and model boundary areas;

[0050] S52: Use rotary drilling method to conduct core drilling at representative locations, record core characteristics, and analyze the distribution characteristics of soluble rocks and the degree of karst development in the core;

[0051] S53: Based on the core drilling results, adjust the boundaries of the karst development area, update the spatial distribution of the soluble rock 3D geological body, and mark the karst development degree to update the 3D geological model. The karst development degree is divided into strong development, moderate development, and weak development.

[0052] Based on the above technical solution, preferably, step S6 includes:

[0053] For karst development areas with a medium or higher degree of karst development, avoid them on the plane;

[0054] If the tunnel needs to pass through a karst development area with a moderate or higher karst development degree, the route should be set at a large angle and short distance to quickly pass through the karst development area;

[0055] According to the spatial distribution of the three-dimensional soluble rock geological body, the longitudinal section of the line is optimized, and the tunnel is selected to pass through the karst vertical seepage zone;

[0056] For areas where karst development is not observed in limestone interlayers, detours are still carried out.

[0057] The method of the present invention has the following beneficial effects compared to the prior art:

[0058] (1) It solves the shortcoming of poor accuracy in analyzing the spatial distribution and development degree of karst by combining conventional drilling with two-dimensional geophysical profiles, and realizes the three-dimensional exploration of the karst development in limestone interlayers in tunnel engineering, providing an accurate basis for geological line selection and engineering design and construction of tunnel engineering;

[0059] (2) During the implementation process, the karst distribution of limestone interlayers can be determined without a large amount of on-site drilling, which shortens the survey period and reduces engineering investment. It is economical and applicable, and has important social and economic significance and practical value.

[0060] (3) The proposed comprehensive survey method is simple to construct and easy to operate. It can also scientifically and rationally solve the problem of high difficulty in surveying limestone interlayers in tunnel engineering, and has a good prospect of being widely used in engineering surveys. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0062] Figure 1 is a flow chart of the method of the present invention;

[0063] Figure 2 It is a technical implementation diagram of the present invention;

[0064] Figure 3 This is a schematic diagram of the survey line arrangement scheme of the present invention;

[0065] Figure 4 is a two-dimensional cross-sectional view of the present invention;

[0066] Figure 5 It is a schematic diagram of the three-dimensional cross-section of the limestone interlayer of the present invention. DETAILED DESCRIPTION

[0067] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0068] See also Figure 1 The present invention provides a comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity, comprising the following steps:

[0069] S1: Collect regional geological data of the tunnel site, analyze and preliminarily determine the regional geological data, and determine the karst development area;

[0070] S2: Select typical areas in the tunnel site to conduct high-density resistivity testing and establish a database of resistivity values ​​for different rock and soil bodies;

[0071] S3: Design a survey line layout plan based on the tunnel direction of the tunnel site, measure the two-dimensional profile resistivity of the tunnel site based on the survey line layout plan, and obtain two-dimensional profile resistivity data;

[0072] S4: Preprocess the 2D profile resistivity data based on the resistivity database of different rock and soil bodies, and construct a 3D geological model. Extract 3D anomalies and soluble rock 3D geological bodies based on the 3D geological model and karst development areas;

[0073] S5: Drill and verify the 3D anomaly bodies and soluble rock 3D geological bodies, revise the karst development area based on the verification results, and update the 3D geological model;

[0074] S6: Combine the updated 3D geological model to guide engineering geological line selection.

[0075] See also Figure 2 The technical solution adopted by the present invention is: first, the karst development location is preliminarily determined based on regional geological data analysis, geological surveys, etc.; then, geophysical detection lines are arranged in the survey area, and based on the two-dimensional geophysical profile data, the principle of strip division is used to form a three-dimensional resistivity contour line body, and according to the resistivity difference between the limestone and other rock cores, the position of the limestone interlayer and the degree of karst development of the limestone interlayer are determined. Then, a small amount of drilling is performed to correct and improve the geophysical results. After comprehensive analysis, the spatial position of the limestone interlayer and the degree of karst development can be revealed, forming a comprehensive survey result to guide engineering line selection and formulate appropriate engineering disposal plans.

[0076] Specifically, in one embodiment of the present invention, step S1 includes:

[0077] S11: Collect regional geological data of the tunnel site, including regional geological maps, geological reports, drilling data and auxiliary information;

[0078] Regional geological maps: including geological maps at a scale of 1:50,000 or larger. Geological reports: including regional geological survey reports and hydrogeological reports. Borehole data: collection of existing borehole data in and around the tunnel site. Supporting materials: such as satellite images, aerial photographs, and historical geological disaster records.

[0079] S12: Develop a survey plan based on regional geological data and conduct on-site geological surveys based on the survey plan, including field investigations, topographic and geomorphological surveys, and establishing an image database;

[0080] Based on available data, determine key survey areas and routes. Observe and record surface rock outcrops, paying particular attention to the distribution of limestone layers. Measure the rock formations and draw simple cross-sections. Collect representative rock samples. Record topographic features, such as mountain orientation and river distribution. Observe geomorphic features potentially associated with karst development, such as funnels and caves. Detailed photography of important geological phenomena will be conducted to establish an image database.

[0081] S13: Organize and standardize the drilling data and draw a drilling histogram. Based on the drilling histogram, analyze the characteristics and distribution of the limestone layer and preliminarily determine the degree of karst development.

[0082] Enter the borehole data into a database. Standardize the coordinate system and depth units. Use specialized software (such as LogPlot) to create standardized borehole histograms. Calculate the depth and thickness of the limestone layers in each borehole. Analyze the continuity and changing trends of the limestone layers. Based on information about caves, fissures, and other features found in the borehole data, make a preliminary assessment of the extent of karst development.

[0083] S14: Establish a GIS project based on the regional geological map and drill hole data, outline the distribution range of the limestone layer in the GIS layer, and use the interpolation method to estimate the thickness distribution of the limestone layer;

[0084] Create a new project using software such as ArcGIS or QGIS. Import the regional geological map and drillhole data into the GIS. Create a polygon layer in the GIS to outline the planar distribution of the limestone layers. Estimate the thickness of the limestone layers using kriging or inverse distance weighted methods.

[0085] S15: Based on topographic features, drilling data and field geological survey results, mark potential karst development areas and use different symbols or colors to indicate the degree of karst development;

[0086] Identify topographic features potentially associated with karst development, such as funnels and karst gullies. Annotate information such as karst caves and fissures discovered during drilling on a GIS layer. Add karst phenomena discovered during field surveys to the GIS. Use different symbols or colors to indicate areas of weak, moderate, and strong karst development.

[0087] S16: According to step S14 and step S15, a limestone distribution range map is generated, and a karst development area distribution map is produced.

[0088] Based on the results of S14, generate a map containing the distribution range and thickness contours of the limestone. Integrate the karst development information annotated in S15 to generate a karst development regional distribution map. Save this map as a high-resolution image in an editable GIS file format.

[0089] Specifically, in one embodiment of the present invention, step S2 includes:

[0090] S21: Select a typical area within the tunnel site based on regional geological data;

[0091] Analyze regional geological data to determine the main rock and soil types and distribution within the tunnel site. Select areas that represent different terrain characteristics. Select a representative area that covers all major rock and soil types and geological characteristics.

[0092] S22: Sampling of rock and soil in typical areas, where the sampling points cover all rock and soil types within the tunnel site;

[0093] Design sampling plan: Develop a detailed sampling plan based on the geological characteristics of typical areas.

[0094] Determine sampling points: Select sampling points that can represent various types of rock and soil in typical areas.

[0095] Selection of sampling method: Select appropriate sampling method (such as core drilling, surface sampling, etc.) according to the type of rock and soil.

[0096] Sample collection:

[0097] For rocks: Collect fresh, unweathered rock samples.

[0098] For soil: collect undisturbed soil samples.

[0099] Sample numbering and recording: Assign a unique number to each sample and record information such as sampling location, depth, and lithology description.

[0100] S23: Conduct high-density resistivity testing on the collected samples, calculate the resistivity values ​​of various rock and soil bodies, and draw box plots to show the resistivity distribution characteristics;

[0101] Use a high-precision resistivity tester. Test method: Use the four-electrode method to test and ensure good contact between the electrodes and the sample.

[0102] S24: Obtain reference standards for resistivity ranges of different rock and soil bodies, and determine resistivity ranges of various rock and soil bodies based on the test results and the reference standards, so as to establish a resistivity value database for different rock and soil bodies.

[0103] Specifically, the reference standard adopted is Appendix A of the "Railway Engineering Object Exploration Specifications":

[0104]

[0105] Compare test results to reference standards. Analyze differences between the resistivity characteristics of local geotechnical materials and standard values. Based on test results and reference standards, determine a reasonable resistivity range for each geotechnical material type. Create a database table containing fields such as geotechnical material type, resistivity range, and typical value.

[0106] Specifically, in one embodiment of the present invention, step S3 includes:

[0107] S31: Design a survey line layout plan, where the geophysical detection lines are parallel to the tunnel axis, the spacing between the geophysical detection lines is D, and n geophysical detection lines are arranged according to the spacing D. The survey line layout plan covers the karst development area;

[0108] First, determine the precise tunnel axis direction based on the engineering design drawings. The distribution of karst areas determines the required survey coverage. Survey lines should extend a certain distance on both sides of the tunnel axis to obtain more comprehensive geological information. The distance D between survey lines should be set between 10 and 50 meters. A total of n survey lines should be set, ensuring that the number of survey lines is sufficient to cover the entire karst area.

[0109] In this embodiment, the schematic diagram of the survey line arrangement scheme is as follows Figure 3 According to the designed survey line layout plan, draw a survey line layout diagram and mark the coordinates of the starting and ending points of each survey line.

[0110] S32: Use GPS or total station to locate the position of the geophysical detection line, arrange electrodes at the designed geophysical detection line position according to the spacing D, and record the coordinates of each electrode;

[0111] For high-precision requirements, use RTK-GPS or a total station. In areas with poor GPS signals, a total station is preferred.

[0112] Survey the survey line on-site and arrange the electrodes along the survey line at the predetermined spacing D. Use GPS or a total station to accurately measure the three-dimensional coordinates (X, Y, Z) of each electrode. Record the electrode number and corresponding coordinate information. Connect the electrode cables according to the requirements of the high-density resistivity meter.

[0113] S33: Use a high-density resistivity meter to perform multiple rounds of resistivity testing, and take the average value of the multiple test results as the resistivity data;

[0114] The instrument parameters are configured based on the measurement line length and electrode spacing. In this embodiment, the Wenner method or the dipole-dipole method is selected as the measurement method, and at least three independent tests are performed on each measurement line.

[0115] During the resistivity test, if the resistivity value of the current round is abnormal, the abnormal area is located, and additional geophysical detection lines are added. The spacing D between the geophysical detection lines in this area is set to 3-5m.

[0116] S34: forming cross-sectional data of the tunnel site region based on the coordinates of the electrodes, and generating two-dimensional cross-sectional resistivity data by combining the resistivity data and the cross-sectional data.

[0117] Import the coordinate information of all electrodes into the data processing software. Based on the electrode coordinates, construct a two-dimensional profile frame for each survey line. Use an appropriate interpolation method (such as Kriging) to convert the point resistivity data into continuous profile data. Map the interpolated resistivity data onto the profile frame. Use a color scale to represent different resistivity values. Generate a high-quality two-dimensional resistivity profile. The two-dimensional profile in this embodiment is as follows: Figure 4 shown.

[0118] Specifically, in one embodiment of the present invention, step S4 includes:

[0119] S41: importing the survey line terrain data, correcting the two-dimensional profile resistivity data according to the terrain correction algorithm, and removing noise using the filtering algorithm to obtain the two-dimensional profile resistivity data after preliminary processing;

[0120] Specifically, terrain elevation data for each survey line is imported into processing software. A terrain correction algorithm (such as the parallel electric field method) is used to correct the raw data. The resistivity distribution is adjusted to account for the impact of terrain undulation on current paths. An appropriate filtering algorithm (such as a median filter or wavelet transform) is applied to remove noise from the data. This yields preliminarily processed two-dimensional resistivity profile data.

[0121] S42: performing inversion calculation on the two-dimensional profile resistivity data after preliminary processing using an inversion algorithm, and outputting the inverted two-dimensional profile resistivity data;

[0122] Select inversion algorithm: Select an appropriate inversion algorithm based on the data characteristics, such as least squares method or smoothing constrained inversion.

[0123] Set inversion parameters: define parameters such as grid size, number of iterations, damping factor, etc. Consider prior information, such as the resistivity range of known geological units.

[0124] Perform inversion calculations: Use RES2DINV to perform inversion calculations. Monitor the iteration process to ensure convergence.

[0125] Evaluate inversion results: Check fitting errors and model resolution. Generate a 2D cross-sectional resistivity map after inversion.

[0126] S43: Compare the inverted 2D profile resistivity data with a database of resistivity values ​​for different rock and soil bodies, quantify the differences using statistical methods, and mark areas with significant differences;

[0127] Recall the geotechnical resistivity database established in step S24. Compare the inversion results with the resistivity range in the database. Quantify the differences using statistical methods (such as the Z-score or chi-square test). Define a significance threshold for the differences. Mark areas on the profile that exceed the threshold.

[0128] S44: Correct the areas with significant differences and return to step S42 until the inversion result is consistent with the resistivity database, and output the processed two-dimensional profile resistivity data;

[0129] Evaluate whether the discrepancy is caused by geological anomalies or data processing errors. Adjust inversion parameters (such as constraints and weights) for areas with significant discrepancies. Return to step S42 and re-perform the inversion calculation using the adjusted parameters. Repeat step S43 to check the new inversion results. Repeat the above process until the inversion results are consistent with the resistivity database. Set a maximum number of iterations to avoid infinite loops.

[0130] S45: Input the processed two-dimensional profile resistivity data into MATLAB and create a structure array, where each structure corresponds to a two-dimensional profile;

[0131] Convert the processed 2D cross-section resistivity data into a MATLAB-compatible format. Create a structure for each 2D cross section. The structure fields include: cross section ID, X coordinate, Y coordinate, depth, resistivity value, etc.

[0132] S46: Determine the grid range in the X, Y, and Z directions according to the range of the structure, and set the grid resolution. Merge all processed two-dimensional profile resistivity data into a three-dimensional point set. Each point in the three-dimensional point set includes its X, Y, and Z coordinates and the corresponding resistivity value. Use the interpolation algorithm to perform linear interpolation on the three-dimensional point set to obtain three-dimensional voxel data.

[0133] Analyze the spatial extent of all sections and determine the maximum and minimum values ​​in the X, Y, and Z directions. Select an appropriate grid resolution based on survey accuracy requirements and computing resources. Combine all section data into a single 3D point set. Each point contains X, Y, and Z coordinates and a resistivity value. Interpolate the point set using a linear interpolation algorithm (such as trilinear interpolation or kriging). This generates regular 3D voxel data.

[0134] S47: setting a resistivity threshold of the anomaly, extracting an isosurface from the three-dimensional voxel data according to the resistivity threshold of the anomaly, and taking the area surrounded by the isosurface as a preliminary anomaly, wherein the isosurface is a collection of points with the same resistivity value; performing connectivity analysis on the preliminary anomaly, classifying the connected areas as one anomaly, and obtaining a three-dimensional anomaly;

[0135] Determine the resistivity threshold for the anomaly based on the geological background and exploration objectives. Use the MATLAB isosurface function to extract the isosurface. Mark the area enclosed by the isosurface as a preliminary anomaly.

[0136] Connectivity analysis of preliminary anomalies:

[0137] S471: Convert the 3D voxel data into binary data, assign the data points on the preliminary abnormal volume a value of 1, and assign the other data points in the 3D voxel data a value of 0;

[0138] Use a triple loop to iterate over all voxels in the X, Y, and Z directions. For each voxel, check whether its resistivity value is within the defined anomaly range. If the voxel belongs to the preliminary anomaly, assign it a value of 1. Otherwise, assign it a value of 0. Generate a binary matrix of the same size as the original 3D voxel data.

[0139] S472: Create a marker matrix of the same size as the 3D voxel data, with all initial values ​​set to 0; set an empty list list to store the voxels to be processed;

[0140] Use the zeros function to create a 3D matrix of the same size as the 3D voxel data. All elements are initialized to 0. Use the MATLAB cell array or queue data structure to create an empty list.

[0141] S473: For each unlabeled voxel in the preliminary abnormal volume, assign it a new unique label and add the voxel to the list;

[0142] Traverse the binary data: Use a triple loop to find voxels with a value of 1 that are not labeled. Assign unique labels: Assign a new unique label (starting from 1 and increasing) to each newly discovered unlabeled anomalous voxel. Update the label matrix: Record the label at the corresponding position in the label matrix. Add to the processing list: Add the coordinates of the voxel to the list. Record label information: Maintain a label counter to ensure that each new anomalous voxel receives a unique label.

[0143] S474: When the list is not empty, take a voxel from the list, check its adjacent voxels, and assign the same label to each adjacent voxel that is located in the preliminary abnormal volume and is not labeled, and add it to the list;

[0144] Use a pop or dequeue operation to remove a voxel from the list. Check the current voxel's neighbors based on 26-connectivity. For each neighboring voxel with a value of 1 and not labeled: a. Assign the same label as the current voxel in the labeling matrix. b. Add the neighboring voxel to the list. Ensure that already processed voxels are not added to the list repeatedly.

[0145] S475: Repeat step S474 until the list is empty;

[0146] S476: Repeat steps S473-S475 until all voxels in the preliminary abnormal volume are labeled, and output a labeling matrix in which the value of each voxel is the label of the preliminary abnormal volume;

[0147] S477: Group voxels with the same label into a 3D anomaly to complete the connectivity analysis. Use the unique function to obtain all unique labels in the label matrix. For each unique label, find all voxels with that label in the label matrix. Create a data structure (such as a structure) for each label, containing: a. Anomaly ID (i.e., label value) b. Voxel coordinate list c. Number of voxels (volume) d. Center of gravity coordinates e. Boundary coordinates (minimum and maximum X, Y, Z values). Generate a data structure containing all anomaly information

[0148] S48: Set the resistivity range of the limestone body, preliminarily extract the limestone body according to the isosurface method, optimize the extraction results using the 3D morphology method, and perform continuity analysis on the optimized extraction results to identify the limestone interlayer, i.e., the 3D soluble rock geological body.

[0149] According to the database of step S24, the resistivity range of the limestone body is determined. The area that meets the resistivity range is extracted using the isosurface method. The extraction results are optimized using dilation and erosion operations. Small isolated areas are removed and internal cavities are filled. The spatial distribution and continuity of the limestone body are analyzed. The possible location of the limestone interlayer is identified. In combination with the results of the abnormal body analysis, the possible karst development area is determined. A three-dimensional geological body model of soluble rock is generated. In this embodiment, the limestone interlayer of the three-dimensional section is as follows Figure 5 shown.

[0150] Specifically, in one embodiment of the present invention, step S5 includes:

[0151] S51: Analyze the spatial distribution of 3D anomalies and soluble rock 3D geological bodies, and select representative locations, including areas with obvious resistivity anomalies, soluble rock distribution areas, and model boundary areas;

[0152] Use VTK to load the 3D anomaly model, observe the spatial distribution, size, and shape of the anomaly, and identify the area with the most significant resistivity anomaly.

[0153] Superimpose the 3D geological model of soluble rocks to identify the areas with the most concentrated distribution of soluble rocks.

[0154] Identify boundary areas of the model, especially areas where anomalous bodies or soluble rocks extend beyond the boundary.

[0155] Select representative locations that cover: a. The area with the most obvious resistivity anomaly b. The area with the most concentrated distribution of soluble rocks c. The key area of ​​the model boundary d. Representative areas at different depths.

[0156] S52: Use rotary drilling method to conduct core drilling at representative locations, record core characteristics, and analyze the distribution characteristics of soluble rocks and the degree of karst development in the core;

[0157] Drilling is carried out at representative locations, using coring equipment to place cores into core boxes in sequence, marking the depth and direction.

[0158] Record detailed information about each borehole, including: a. Borehole coordinates and depth b. Drilling rate variations c. Groundwater levels encountered d. Special geological phenomena (such as cavities, fissures, etc.).

[0159] Provide a detailed description of each core section, including: a. lithology b. color c. texture and structure d. degree of weathering e. fracture characteristics f. type and distribution of soluble rocks.

[0160] The degree of karst development is assessed based on the characteristics of the soluble rocks in the core: a. The size and density of the dissolution pores b. The nature of the fracture fillings c. The integrity and recovery rate of the core.

[0161] Summarize all drilling data and core analysis results.

[0162] S53: Based on the core drilling results, adjust the boundaries of the karst development area, update the spatial distribution of the soluble rock 3D geological body, and mark the karst development degree to update the 3D geological model. The karst development degree is divided into strong development, moderate development, and weak development.

[0163] Compare drilling results with the original 3D geological model. Mark the drill hole locations and actual geological information.

[0164] Based on the core analysis results, the boundaries of the karst development area are redefined.

[0165] Based on the actual soluble rock distribution verified by drilling, the soluble rock 3D geological model is modified to adjust the spatial extent and morphology of the soluble rock.

[0166] Based on the results of core analysis, the degree of karst development is divided into three levels: a. Strong development: a large number of dissolution pores, highly fragmented or fissure-developed b. Moderate development: obvious dissolution features, but on a small scale c. Weak development: slight signs of dissolution, and relatively intact structure.

[0167] Use appropriate interpolation methods (such as Kriging) to expand the point-like drilling data to the entire model space, generating a continuous distribution model of karst development degree.

[0168] Calculate the consistency of the updated model with the drilling data.

[0169] Ultimately, the engineering geological route selection is guided by the degree of karst development and the spatial relationship between the limestone interlayer and the tunnel body. Specifically: ① For areas with moderate or higher karst development, avoid them in plan view whenever possible. ② When a tunnel must pass through an area with moderate or higher karst development, the route should be traversed at a high angle and over a short distance, allowing for rapid passage through the karst zone. ③ Based on the spatial distribution of the karst interlayer, optimize the route's longitudinal profile and choose to pass through the karst vertical seepage zone. ④ Even if the limestone interlayer has no karst development but runs parallel to the tunnel for a long distance, avoid it whenever possible.

[0170] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A comprehensive survey method for tunnel engineering in limestone interlayers based on three-dimensional resistivity, characterized by: The steps include: S1: Collect regional geological data of the tunnel site, analyze and preliminarily determine the regional geological data, and determine the karst development area; S2: Select typical areas in the tunnel site to conduct high-density resistivity testing and establish a database of resistivity values ​​for different rock and soil bodies; S3: Design a survey line layout plan based on the tunnel direction of the tunnel site, measure the two-dimensional profile resistivity of the tunnel site based on the survey line layout plan, and obtain two-dimensional profile resistivity data; S4: Preprocess the 2D profile resistivity data based on the resistivity database of different rock and soil bodies, and construct a 3D geological model. Extract 3D anomalies and soluble rock 3D geological bodies based on the 3D geological model and karst development areas; Step S4 includes: S41: importing the survey line terrain data, correcting the two-dimensional profile resistivity data according to the terrain correction algorithm, and removing noise using the filtering algorithm to obtain the two-dimensional profile resistivity data after preliminary processing; S42: performing inversion calculation on the two-dimensional profile resistivity data after preliminary processing using an inversion algorithm, and outputting the inverted two-dimensional profile resistivity data; S43: Compare the inverted 2D profile resistivity data with a database of resistivity values ​​for different rock and soil bodies, quantify the differences using statistical methods, and mark areas with significant differences; S44: Correct the areas with significant differences and return to step S42 until the inversion result is consistent with the resistivity database, and output the processed two-dimensional profile resistivity data; S45: Input the processed two-dimensional profile resistivity data into MATLAB and create a structure array, where each structure corresponds to a two-dimensional profile; S46: Determine the grid range in the X, Y, and Z directions according to the range of the structure, and set the grid resolution. Merge all processed two-dimensional profile resistivity data into a three-dimensional point set. Each point in the three-dimensional point set includes its X, Y, and Z coordinates and the corresponding resistivity value. Use the interpolation algorithm to perform linear interpolation on the three-dimensional point set to obtain three-dimensional voxel data. S47: setting a resistivity threshold of the anomaly, extracting an isosurface from the three-dimensional voxel data according to the resistivity threshold of the anomaly, and taking the area surrounded by the isosurface as a preliminary anomaly, wherein the isosurface is a collection of points with the same resistivity value; performing connectivity analysis on the preliminary anomaly, classifying the connected areas as one anomaly, and obtaining a three-dimensional anomaly; S48: Set the resistivity range of the limestone body, preliminarily extract the limestone body using the isosurface method, optimize the extraction results using the 3D morphology method, and perform continuity analysis on the optimized extraction results to identify the limestone interlayer, i.e., the 3D soluble rock geological body; S5: Drill and verify the 3D anomaly bodies and soluble rock 3D geological bodies, revise the karst development area based on the verification results, and update the 3D geological model; S6: Combine the updated 3D geological model to guide engineering geological line selection.

2. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 1, characterized in that: Step S1 includes: S11: Collect regional geological data of the tunnel site, including regional geological maps, geological reports, drilling data and auxiliary information; S12: Develop a survey plan based on regional geological data and conduct on-site geological surveys based on the survey plan, including field investigations, topographic and geomorphological surveys, and establishing an image database; S13: Organize and standardize the drilling data and draw a drilling histogram. Based on the drilling histogram, analyze the characteristics and distribution of the limestone layer and preliminarily determine the degree of karst development. S14: Establish a GIS project based on the regional geological map and drill hole data, outline the distribution range of the limestone layer in the GIS layer, and use the interpolation method to estimate the thickness distribution of the limestone layer; S15: Based on topographic features, drilling data and field geological survey results, mark potential karst development areas and use different symbols or colors to indicate the degree of karst development; S16: According to step S14 and step S15, a limestone distribution range map is generated, and a karst development area distribution map is produced.

3. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 1, characterized in that: Step S2 includes: S21: Select a typical area within the tunnel site based on regional geological data; S22: Sampling of rock and soil in typical areas, where the sampling points cover all rock and soil types within the tunnel site; S23: Conduct high-density resistivity testing on the collected samples, calculate the resistivity values ​​of various rock and soil bodies, and draw box plots to show the resistivity distribution characteristics; S24: Obtain reference standards for resistivity ranges of different rock and soil bodies, and determine resistivity ranges of various rock and soil bodies based on the test results and the reference standards, so as to establish a resistivity value database for different rock and soil bodies.

4. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 1, characterized in that: Step S3 includes: S31: Design a survey line layout plan, where the geophysical detection lines are parallel to the tunnel axis, the spacing between the geophysical detection lines is D, and n geophysical detection lines are arranged according to the spacing D. The survey line layout plan covers the karst development area; S32: Use GPS or total station to locate the position of the geophysical detection line, arrange electrodes at the designed geophysical detection line position according to the spacing D, and record the coordinates of each electrode; S33: Use a high-density resistivity meter to perform multiple rounds of resistivity testing, and take the average value of the multiple test results as the resistivity data; S34: forming cross-sectional data of the tunnel site region based on the coordinates of the electrodes, and generating two-dimensional cross-sectional resistivity data by combining the resistivity data and the cross-sectional data.

5. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 4, characterized in that: The spacing D is 10-20m.

6. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 4, characterized in that: During the resistivity test, if the resistivity value of the current round is abnormal, the abnormal area is located, and additional geophysical detection lines are added. The spacing D between the geophysical detection lines in this area is set to 3-5m.

7. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 1, characterized in that: In step S47, the process of connectivity analysis is as follows: S471: Convert the 3D voxel data into binary data, assign the data points on the preliminary abnormal volume a value of 1, and assign the other data points in the 3D voxel data a value of 0; S472: Create a marker matrix of the same size as the 3D voxel data, with all initial values ​​set to 0; set an empty list list to store the voxels to be processed; S473: For each unlabeled voxel in the preliminary abnormal body, assign it a new unique label and add the voxel to the list; S474: When the list is not empty, take a voxel from the list, check its adjacent voxels, and assign the same label to each adjacent voxel that is located in the preliminary abnormal volume and is not labeled, and add it to the list; S475: Repeat step S474 until the list is empty; S476: Repeat steps S473-S475 until all voxels in the preliminary abnormal volume are labeled, and output a labeling matrix in which the value of each voxel is the label of the preliminary abnormal volume; S477: Voxels with the same label are grouped into a 3D anomaly volume to complete connectivity analysis.

8. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 1, characterized in that: Step S5 includes: S51: Analyze the spatial distribution of 3D anomalies and soluble rock 3D geological bodies, and select representative locations, including areas with obvious resistivity anomalies, soluble rock distribution areas, and model boundary areas; S52: Use rotary drilling method to conduct core drilling at representative locations, record core characteristics, and analyze the distribution characteristics of soluble rocks and the degree of karst development in the core; S53: Based on the core drilling results, adjust the boundaries of the karst development area, update the spatial distribution of the soluble rock 3D geological body, and mark the karst development degree to update the 3D geological model. The karst development degree is divided into strong development, moderate development, and weak development.

9. The comprehensive survey method for limestone interlayer tunnel engineering based on three-dimensional resistivity according to claim 8, characterized in that: Step S6 includes: For karst development areas with a medium or higher degree of karst development, avoid them on the plane; If the tunnel needs to pass through a karst development area with a moderate or higher karst development degree, the route should be set at a large angle and short distance to quickly pass through the karst development area; According to the spatial distribution of the three-dimensional soluble rock geological body, the longitudinal section of the line is optimized, and the tunnel is selected to pass through the karst vertical seepage zone; For areas where karst development is not observed in limestone interlayers, detours are still carried out.

Images