Ground and cross-hole earthquake combined tomography method for improving karst detection precision

Abstract

The invention discloses a ground and cross-hole earthquake combined tomography method for improving karst detection precision. The method comprises the following steps: arranging a seismic source and a detector on the ground surface to realize ground exploration, respectively arranging a seismic source and a detector in different holes to realize cross-hole exploration, and combining seismic data acquired by the two exploration modes; picking up first-motion wave travel time in each seismic record to prepare for inversion calculation; and establishing an initial velocity model and carrying out inversion iteration, and explaining an inversion result in combination with karst features and geological backgrounds. According to the method, karst distribution under complex geological conditions can be accurately positioned and recognized, and the problems that inversion imaging is uncomfortable and prominent in nature and geological interpretation multiplicity is high due to the fact that cross-hole tomography used in traditional karst detection is uneven in seismic ray distribution, insufficient in cross degree and the like are solved; the method has remarkable advantages in fine exploration of karst in a complex geological environment, and has wide practical value.

Description

Technical Field

[0001] The invention relates to the field of exploration geophysical technology, and in particular to a ground and cross-hole seismic combined tomography method for improving karst detection accuracy. Background Art

[0002] As urban construction in my country accelerates, the demand for transparency of underground space has increased significantly. However, in the process of underground space engineering construction, adverse geological conditions, especially karst problems, seriously affect the quality and safety of the project. The formation and development process of karst is complex, and its diverse forms underground pose challenges to detection. If the karst distribution is not accurately identified during engineering construction in karst areas, major safety accidents such as sudden water inrush and surface collapse may occur. Therefore, it is of great significance to carry out refined detection of karst areas, clarify the scope and distribution characteristics of karst, and realize the visualization of underground geology for the development of urban underground space and the smooth implementation of engineering projects.

[0003] At present, karst exploration usually uses cross-hole data collection. This method is close to the detection target and the seismic wave energy propagation distance is short, so the seismic data obtained has a relatively high signal-to-noise ratio and resolution, which is conducive to clearer imaging of underground karst. However, traditional cross-hole tomography also has obvious disadvantages. Its observation angle is limited, the ray distribution is often uneven, and the cross-degree is insufficient. It cannot meet the accuracy requirements of seismic wave tomography inversion, resulting in prominent instability problems in inversion imaging and strong multi-solution of geological interpretation, making it difficult to achieve refined exploration of karst in complex geological environments. Summary of the invention

[0004] The object of the present invention is to provide a ground and cross-hole seismic combined tomography method for improving the accuracy of karst detection in view of the deficiencies of the prior art.

[0005] The object of the present invention is achieved through the following technical scheme: A ground and cross-hole seismic combined tomography method for improving the accuracy of karst detection, comprising the following steps:

[0006] 1. Arrange seismic sources and geophones on the surface to achieve ground exploration, and arrange seismic sources and geophones in different holes to achieve cross-hole exploration, and combine the seismic data collected by the two exploration methods;

[0007] 2. Pick up the first arrival wave travel time b in each earthquake record i , the travel time of the first arrival wave is b i As observed values;

[0008] 3. Establish the initial velocity model and perform inversion iteration, and use ray tracing to calculate the travel time of the first arrival wave of the model According to the With the bi The difference The velocity model is modified and iterated repeatedly until the travel time error of the first arrival wave reaches the accuracy requirement or the iteration reaches the maximum number of iterations, and the final inversion result is obtained, and the inversion result is interpreted in combination with the karst characteristics and geological background.

[0009] Furthermore, step 3 is implemented by the following sub-steps:

[0010] 3.1. Establish an initial velocity model based on the existing geological information, select a rectangular grid for discretization, arrange the discretized rectangular grid in the order of rows first and columns, and obtain N pixels. Set the slowness in each pixel as a constant. The slowness is expressed as s. j (j=1,2,…,N) means that an image basis function g is defined for each pixel j (j=1,2,…,N) j (x,z):

[0011]

[0012] Among them, j represents the index of the pixel, x represents the horizontal distance, and z represents the depth;

[0013] Then the seismic wave slowness image s(x,z) of the underground medium can be obtained by s j and g j It can be expressed as a linear combination of, as shown in the following formula:

[0014]

[0015] Among them, s(x,z) represents the seismic wave slowness image, g j (x,z) represents the basis function on each pixel, s j represents the constant value within the jth pixel, and N represents the total number of pixels;

[0016] 3.2. Using the parabolic interpolation ray tracing algorithm to find the minimum travel time of each ray and the length a of each ray in each grid ij (i=1,2,…,I,j=1,2,…,N), when the mth ray does not pass through the nth grid, a mn = 0; wherein the ray represents the propagation path of the seismic wave in seismic tomography, and the ray is bent during the propagation process due to the velocity change of the underground medium;

[0017] 3.3. Calculate the observed value b of the first arrival wave travel time i and Difference: Will Substitute into the tomographic equations and solve using the inversion algorithm Get the slowness correction in the tomography during the kth iteration The corrected slowness vector is calculated

[0018] 3.4. Use the corrected slowness vector as the new initial value and repeat steps 3.2 to 3.4 until the error of the first arrival wave travel time reaches the accuracy requirement or the maximum number of iterations is reached, then the iteration is stopped;

[0019] 3.5. Interpret the inversion results based on karst characteristics and geological background.

[0020] Furthermore, in step 3.3, the process of establishing the tomographic equation group is as follows:

[0021] If the first arrival wave ray path from each shot point to each detection point records a total of I rays, each ray is represented by i (i = 1, 2, ..., I), the corresponding first arrival wave travel time b i (i=1,2,…,I) is the curve integral of the first arrival wave slowness s(x,z) along the i-th ray, so b i (i=1,2,…,I) is the generalized Radon transform of s(x,z), as shown in the following formula:

[0022]

[0023] in, is the length of the i-th first arrival ray in the j-th grid, expressed as a ij (i=1,2,…,I,j=1,2,…,N) represents, the following system of equations is obtained:

[0024]

[0025] The above equations are recorded as As = b, where A = (a ij ), A represents the matrix composed of the length of each ray in each grid, s=(s1,s2,...,s N ) T represents the vector composed of the first arrival slowness in each grid; b = (b1, b2, ..., b I ) T represents the vector consisting of the travel time of the first arrival wave of each ray.

[0026] Furthermore, in step 3.3, solve The iterative algorithm used is the damped least squares QR factorization method.

[0027] Furthermore, in the step 3.5, the karst feature is mainly that when karst is developed, the wave velocity value of the karst is significantly lower than that of the surrounding rock mass, resulting in a wave velocity value difference; the wave velocity value difference provides a favorable geophysical condition prerequisite for karst exploration, and whether the rock mass has karst development can be judged based on the change in wave velocity.

[0028] The beneficial effects of the present invention are as follows: the combined ground and cross-hole seismic tomography method adopted by the present invention fully utilizes the complementarity of ground and cross-hole data, effectively overcomes the problems of uneven ray distribution and insufficient cross-degree in traditional single cross-hole tomography, and significantly improves the stability and accuracy of seismic tomography. The method can meet the urgent technical needs of fine detection of underground space, effectively solve the problem of fine detection of small targets, and provide strong support for fine karst exploration in complex geological environments. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 It is a flow chart of the combined surface and cross-hole seismic tomography method;

[0030] Figure 2 It is a schematic diagram of the scope of the work area and the layout of different observation systems, among which (a) is a schematic diagram of the ground exploration observation system, and (b) is a schematic diagram of the cross-hole seismic exploration observation system;

[0031] Figure 3 is the initial velocity model and gridding schematic diagram;

[0032] Figure 4 This is the result of the combined surface and cross-hole seismic tomography method;

[0033] Figure 5 are the tomography results and ray distribution diagrams of cross-hole seismic, where (a) is the tomography diagram and (b) is the ray distribution diagram;

[0034] Figure 6 These are the combined tomography results of ground and cross-hole seismic and the ray distribution map, where (a) is the tomography map and (b) is the ray distribution map. DETAILED DESCRIPTION

[0035] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0036] This embodiment discloses a combined ground and cross-hole seismic tomography method for improving karst detection accuracy, such as Figure 1 As shown, the specific steps include:

[0037] 1. Arrange seismic sources and geophones on the surface to realize ground exploration, and arrange seismic sources and geophones in different holes to realize cross-hole exploration, and combine the seismic data collected by the two exploration methods;

[0038] In this embodiment, the work area is located in Puyan Street, Binjiang District, Hangzhou City, Zhejiang Province. A survey line is arranged in the work area based on the drilling data. The survey line length is X=124m and the depth is Z=41m. The sources and detectors of ground and cross-hole earthquakes are arranged for data collection, respectively: ① Ground earthquake: The layout of the sources and detectors is as follows Figure 2 As shown in (a), the source (red "*" symbol) is excited on the surface, and the detectors (blue "." symbol) are also arranged on the surface with a spacing of 2m; ② Cross-hole earthquake: The layout of the source and the detector is as follows Figure 2 As shown in (b), the source (red "*" symbol) is excited in the borehole (42m of the survey line), and the detector (blue "." symbol) is deployed in the borehole (28m of the survey line) to receive it.

[0039] 2. Pick up the travel time of the first arrival wave in the earthquake record i , as the observed value.

[0040] 3. Establish an initial velocity model and perform inversion iterations. Rays represent the propagation path of seismic waves in seismic tomography. During propagation, they bend due to the change in velocity of the underground medium. Ray tracing is used to calculate the travel time of the first arrival wave of the model. according to The observed first arrival wave travel time b i The difference r i k Modify the velocity model and repeat the iteration until the first arrival wave travel time error reaches the accuracy requirement or the iteration reaches the maximum number of iterations, and then obtain the final inversion result, which is interpreted in combination with the karst characteristics and geological background. It includes the following sub-steps:

[0041] (3.1) The initial velocity model is established based on the existing geological information. A rectangular grid is selected for discretization. The discretized rectangular grid is arranged in the order of rows first and columns to obtain N pixels. The slowness in each pixel is set as a constant and expressed as s. j (j=1,2,…,N) indicates that the slowness is the inverse of the speed. For each pixel j (j=1,2,…,N), an image basis function g is defined j (x,z):

[0042]

[0043] Among them, j represents the index of the pixel, x represents the horizontal distance, and z represents the depth.

[0044] Then the seismic wave slowness image s(x,z) of the underground medium can be obtained by s j and g j It can be expressed as a linear combination of:

[0045]

[0046] Among them, s(x,z) represents the seismic wave slowness image, g j (x,z) represents the basis function on each pixel, s j represents the constant value (speed value) within the j-th pixel, and N is the total number of pixels.

[0047] In this embodiment, the initial velocity model of ground exploration is as follows: Figure 3 As shown, the grid size (dx, dz) is (0.5m, 0.5m), the velocity variation range is 400m / s-4000m / s, and the velocity is increased every 0.5m along the depth Z direction to construct an initial velocity model with continuously changing velocity.

[0048] (3.2) Use the parabolic interpolation ray tracing algorithm to find the minimum travel time of each ray and the length a of each ray in each grid ij (i=1,2,…,I,j=1,2,…,N), this information is used to invert the velocity distribution of the underground medium, so as to understand the characteristics of the geological structure. When the mth ray does not pass through the nth grid, a mn =0.

[0049] (3.3) Calculate the observed value b of the first arrival wave travel time i With calculated value difference Substitute into the tomographic equations and solve using the least squares QR factorization method The slowness correction in tomography during the kth iteration can be obtained as The new slowness correction vector is calculated

[0050] The process of establishing the tomographic equations is as follows:

[0051] If the first arrival wave ray path from each shot point to each detection point records a total of I rays, each ray is represented by i (i = 1, 2, ..., I), the corresponding first arrival wave travel time b i (i=1,2,…,I) is the curve integral of the first arrival wave slowness s(x,z) along the i-th ray, so b i (i=1,2,…,I) is the generalized Radon transform of s(x,z), as shown in the following formula:

[0052]

[0053] in, is the length of the i-th first-arrival ray in the j-th grid, expressed as a ij (i=1,2,…,I,j=1,2,…,N) represents, the following system of equations is obtained:

[0054]

[0055] The above equations are recorded as As = b, where A = (a ij ), A represents the matrix composed of the length of each ray in each grid, s=(s1,s2,...,s N ) T represents the vector composed of the first arrival wave slowness in each grid; b = (b1, b2, ..., b I ) T represents the vector consisting of the travel time of the first arrival wave of each ray.

[0056] (3.4) The corrected slowness vector is used as the new initial value, and steps (3.2) to (3.4) are repeated until the travel time error of the first arrival wave reaches the accuracy requirement or the iteration reaches the maximum number of iterations, then the iteration is stopped; wherein the accuracy requirement and the maximum number of iterations can be set by oneself.

[0057] In this embodiment, the maximum number of iterations is set to 12, and the accuracy requirement of the root mean square error is 10 -6 .

[0058] (3.5) Interpret the inversion results based on karst characteristics and geological background.

[0059] The karst characteristics are mainly that when karst is developed, the wave velocity value of the karst is significantly lower than that of the surrounding rock mass. This difference in velocity provides a favorable geophysical condition prerequisite for karst exploration, and whether the rock mass has karst development can be judged based on the change in wave velocity.

[0060] In this embodiment, the combined tomographic imaging results of ground and cross-hole seismic (28m-42m area of ​​the survey line) are as follows: Figure 4 As shown, Figure 4 The velocity of the black frame marked in the figure is significantly lower than that of the surrounding areas at the same depth, so it can be inferred that it is a karst development area, and the drilling verification is consistent with the inversion result. The inversion result provides a strong basis for the evaluation of the karst development status in the region, which helps to guide risk assessment and decision-making in engineering construction.

[0061] 4. In this embodiment, to further illustrate the advantages of the present invention over the prior art, the cross-hole seismic data of step 1 is taken, and steps 2 and 3 are repeated. The initial velocity model is taken as Figure 3 In the area between 28m and 42m of the survey line, the other designs are the same, and conventional cross-hole tomography is performed. The imaging results are similar to the ray distribution. Figure 5 As shown, Figure 6 The combined surface and cross-hole seismic tomography results are compared with the ray distribution map as follows:

[0062] ①Comparison of imaging results: Figure 5 (a) Conventional single-well tomography and Figure 6 (a) Both ground and cross-hole seismic tomography can show the presence of karst near the depth of 30m-40m. Figure 5 (a) The low-velocity anomaly area shown is larger than the actual karst distribution area, and shows obvious anomaly symmetry near the depth of 30 to 35 meters, which is inconsistent with the actual geological conditions. Figure 6 (a) The range of the low-velocity anomaly area is consistent with the actual karst distribution range. Secondly, combined with the drilling data, there is a harder limestone area at a depth of 25 to 30 meters. Figure 5 (a) Figure 6 In (a), the area shows obvious high-speed anomalies, which can be better identified and distinguished from karst. In addition, the combined tomography results clearly show the undulating interface between the shallow surface (0-10m) topsoil and the gravel-containing soil layer, which more accurately depicts the structure of the shallow surface medium, while cross-hole tomography alone fails to accurately depict the interface morphology.

[0063] ② Comparison of ray distribution diagrams: Figure 5 (b) The ray coverage in the shallow surface is obviously insufficient. The nonlinear inversion algorithm involved in tomography is very sensitive to the quality of the input data. The insufficient shallow ray distribution means that during the data inversion process, the shallow layer constraints are insufficient, resulting in a decrease in the inversion accuracy of the entire model. Figure 6 (b) The ray distribution is relatively uniform, which is why the combined ground and cross-hole seismic tomography method has more advantages in karst fine-grained detection than traditional single cross-hole imaging.

[0064] Those skilled in the art can understand that the above are only preferred examples of the invention and are not intended to limit the invention. Although the invention is described in detail with reference to the above examples, those skilled in the art can still modify the technical solutions recorded in the above examples or replace some of the technical features therein with equivalents. Any modifications, equivalent replacements, etc. made within the spirit and principles of the invention should be included in the protection scope of the invention.

Claims

1. A combined ground and cross-hole seismic tomography method for improving karst detection accuracy, characterized in that: The following steps are involved: (1) Arrange seismic sources and geophones on the surface to achieve surface exploration, and arrange seismic sources and geophones in different holes to achieve cross-hole exploration, and combine the seismic data collected by the two exploration methods; (2) Pick the travel time b of the first arrival wave in each earthquake record i , the picked first arrival wave travel time b i As observed values; (3) Establish an initial velocity model and perform inversion iterations, and use ray tracing to calculate the travel time of the first arrival wave of the model According to the With the b i The difference The velocity model is modified and iterated repeatedly until the travel time error of the first arrival wave reaches the accuracy requirement or the iteration reaches the maximum number of iterations, and the final inversion result is obtained, and the inversion result is interpreted in combination with the karst characteristics and geological background.

2. A combined ground and cross-hole seismic tomography method for improving karst detection accuracy according to claim 1, characterized in that: The step (3) is implemented by the following sub-steps: (3.1) An initial velocity model is established based on the existing geological information. A rectangular grid is selected for discretization. The discretized rectangular grid is arranged in the order of rows first and columns to obtain N pixels. The slowness in each pixel is set to a constant. The slowness is expressed as s. j (j=1,2,…,N) means that an image basis function g is defined for each pixel j (j=1,2,…,N) j (x,z): Among them, j represents the index of the pixel, x represents the horizontal distance, and z represents the depth; Then the seismic wave slowness image s(x,z) of the underground medium can be obtained by s j and g j It can be expressed as a linear combination of, as shown in the following formula: Among them, s(x,z) represents the seismic wave slowness image, g j (x,z) represents the basis function on each pixel, s j represents the constant value within the jth pixel, and N represents the total number of pixels; (3.2) Use the parabolic interpolation ray tracing algorithm to find the minimum travel time of each ray and the length a of each ray in each grid ij (i=1,2,…,I,j=1,2,…,N), when the mth ray does not pass through the nth grid, a mn = 0; wherein the ray represents the propagation path of the seismic wave in seismic tomography, and the ray is bent during the propagation process due to the velocity change of the underground medium; (3.3) Calculate the observed value b of the first arrival wave travel time i and Difference: Will Substitute into the tomographic equations and solve using the inversion algorithm Get the slowness correction in the tomography during the kth iteration The corrected slowness vector is calculated (3.4) Using the corrected slowness vector as a new initial value, repeat steps (3.2) to (3.4) until the error in the travel time of the first arrival wave reaches the accuracy requirement or the maximum number of iterations is reached, then the iteration is stopped; (3.5) Interpret the inversion results based on karst characteristics and geological background.

3. A combined ground and cross-hole seismic tomography method for improving karst detection accuracy according to claim 2, characterized in that: In the step (3.3), the process of establishing the tomographic equations is as follows: If the first arrival wave ray path from each shot point to each detection point records a total of I rays, each ray is represented by i (i = 1, 2, ..., I), the corresponding first arrival wave travel time b i (i=1,2,…,I) is the curve integral of the first arrival wave slowness s(x,z) along the i-th ray, so b i (i=1,2,…,I) is the generalized Radon transform of s(x,z), as shown in the following formula: in, is the length of the i-th first-arrival ray in the j-th grid, expressed as a ij (i=1,2,…,I,j=1,2,…,N) represents, and the following system of equations is obtained: The above equations are recorded as As = b, where A = (a ij ), A represents the matrix composed of the length of each ray in each grid, s=(s1,s2,...,s N ) T represents the vector composed of the first arrival slowness in each grid; b = (b1, b2, ..., b I ) T represents the vector consisting of the travel time of the first arrival wave of each ray.

4. A combined ground and cross-hole seismic tomography method for improving karst detection accuracy according to claim 2, characterized in that: In step (3.3), solve The iterative algorithm used is the damped least squares QR factorization method.

5. The combined ground and cross-hole seismic tomography method for improving karst detection accuracy according to claim 2, characterized in that: In the step (3.5), the karst characteristic is that when karst is developed, the wave velocity value of the karst is significantly lower than that of the surrounding rock mass, forming a wave velocity value difference.

Images