Stratum structure and ground pressure distribution detection method based on wave velocity tomography

Abstract

The invention provides a stratigraphic structure and ground pressure distribution detection method based on wave velocity tomography, and relates to the technical field of wave velocity tomography. The method comprises the following steps: firstly, arranging dense linear detector arrays in a monitoring area at an interval of 5m, and collecting the distribution condition of an aquifer, a lithology detection result and a drill core diagram of the monitoring area; recording seismic wave data, carrying out noise reduction processing on the seismic wave data, picking up arrival time, and associating the arrival time with a waveform event; performing path tracking on an observation result by adopting a wave velocity inversion algorithm, generating a wave velocity tomography image and substituting the wave velocity tomography image into the wave velocity-ground pressure evaluation model; and finally generating a ground pressure distribution diagram.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of wave velocity tomography technology, and particularly relates to a stratum structure and ground pressure distribution detection method based on wave velocity tomography. BACKGROUND

[0002] With the increase of mining depth, complex geological conditions such as high ground stress bring great challenges to construction safety. Traditional geophysical prospecting methods such as ground penetrating radar and borehole sampling have been difficult to meet the requirements of fine construction. Geological radar, direct current method and transient electromagnetic method can locate a small area, have weak anti-interference ability, and can only detect a depth of about 150 meters underground, which is far from meeting the needs of deep engineering. Geological anomaly areas such as goaf, fault fracture and aquifer are prone to cause rock burst, water inrush and other dangerous accidents, and if precise detection technology is not used, only relying on experience design may cause engineering delay and even major accidents.

[0003] Wave velocity tomography technology is a technology for three-dimensional high-precision inversion and dynamic monitoring of stratum internal structure by analyzing seismic wave propagation data. This technology can accurately identify geological anomaly areas such as goaf, fracture zone and aquifer, so as to optimize the design of shaft and roadway development construction and reduce the risk of accidents by using the detection results.

[0004] At the same time, since there is a power function relationship between wave velocity and ground pressure, dynamic monitoring of the position of stress concentration and rock burst in shafts, roadways and other places using wave velocity tomography technology can further reveal the ground pressure distribution in the surrounding rock, so as to provide early warning and prevention and control before rock burst, roof fall and other disasters occur, and improve the safety of mine production. Especially for deep mining, by monitoring the wave velocity for a long time in the whole development, preparation and mining area, the wave velocity and ground pressure inversion image covering the whole area and the whole life cycle of mining can be obtained, and then the whole area ground pressure evolution characteristics of the whole life cycle of deep mining can be obtained. SUMMARY

[0005] In view of the deficiencies of the prior art, the present application provides a stratum structure and ground pressure distribution detection method based on wave velocity tomography.

[0006] A stratum structure and ground pressure distribution detection method based on wave velocity tomography, comprising the following steps:

[0007] Step 1: densely arranging linear geophone arrays at intervals of five meters in the monitoring area;

[0008] Step 2: collecting the aquifer distribution, lithology detection results and borehole core map of the monitoring area;

[0009] Step 2.1: establishing an initial wave velocity model according to the collected lithology detection results of the monitoring area;

[0010] The initial wave velocity model is a geological profile velocity model at the exploration line generated according to the experimental wave velocity of seismic waves in different media;

[0011] Step 2.2: According to the collected aquifer distribution, the initial wave velocity model generated in step 2.1 is first corrected, and the wave velocity of the aquifer area is lowered to 1500 m / s;

[0012] Step 2.3: According to the collected drilling core map, the RQD value is calculated, the occurrence of stratigraphic structure and faults in each exploration line is inferred, and the second correction is made to the velocity model after the first correction in step 2.2, and the wave velocity of the area containing structure and fault is lowered to 2800 m / s;

[0013] Step 3: Record the seismic wave data, record the period greater than or equal to 2 days, set the acquisition frequency to 10 Hz, and record the waveform data for 2 Hz-10 Hz;

[0014] Step 4: Denoising processing and time picking of seismic wave data, and correlation with waveform events;

[0015] Step 4.1: After denoising processing of the detected seismic wave data, the P-wave arrival time and S-wave arrival time are picked up at the same time;

[0016] Step 4.2: Automatically match the seismic wave data with the waveforms in the PEER waveform template library to generate waveform events;

[0017] Step 4.3: According to the waveform arrival time sequence picked up in step 4.1, correlate with the waveform events generated in step 4.2.

[0018] Step 5: Use the wave velocity inversion algorithm to track the path of the observation results;

[0019] Step 5.1: Assume that the seismic wave propagates along the path with the shortest propagation time, then the seismic wave satisfies Fermat's principle:

[0020] (1);

[0021] In the formula, t is time, l is ray position, and v(r) is wave velocity;

[0022] According to the extended form of Snell's law, the process of adjusting the ray direction with the velocity gradient satisfies the differential equation:

[0023] (2);

[0024] In the formula, v is the medium velocity, l is the ray position, s is the ray arc length, and ▽ is the gradient operator.

[0025] The formula (2) is substituted into the wave velocity model after the second correction in step 2.3 to perform iterative correction;

[0026] Step 5.2: Discretize the monitoring area into a grid to obtain a shortest path tree;

[0027] From the source, first calculate the travel time of each node adjacent to the source point, take the secondary node as a sub-wave source point, calculate the travel time of each node adjacent to the sub-wave source point at the next level of the sub-wave source point, and select the ray path with the minimum travel time to form a shortest path tree; the ray path and travel time from the source point to all receiving points are obtained from the position of the receiving point; wherein the wave velocity is regarded as uniform in the grid, and the travel time is equal to the Euclidean distance between two adjacent nodes divided by the wave velocity in the grid;

[0028] Step 5.3: Perform path smoothing processing on the shortest path tree obtained in step 5.2 to obtain a simulated ray path;

[0029] Step 5.4: After obtaining the simulated ray path, a theoretical best matching model is defined using a full waveform inversion algorithm :

[0030] (3);

[0031] wherein, is the theoretical best matching model, i is the receiver number, n src represents the number of sources, F i (m) represents the forward synthetic data, d i is the observation data; is the error function of the observation data and the synthetic data;

[0032] The wave field residual of the receiver point is assumed to be an accompanying source, and it is considered that this accompanying source releases an accompanying wave field V; then the accompanying wave field exists in the time domain and the frequency domain;

[0033] The equation of the time domain wave motion is as follows:

[0034] (4);

[0035] wherein, m(x) is the time domain of the wave, x is a particle in the wave propagation path, s(x,t) and u(x,t) represent the source term and the wave field in the time domain respectively, and Δ represents the Laplace operator;

[0036] The equation of the frequency domain wave motion is as follows:

[0037] (5);

[0038] In the formula, B(s, ω) is the frequency domain of the wave, ω is the angular frequency of the wave, s(x, ω) and u(x, ω) represent the source term and the wave field in the frequency domain;

[0039] Step 5.5: Solve the minimum value of the residual, that is, the minimum value of the adjoint wave field v;

[0040] Calculate the gradient vector g of the adjoint wave field:

[0041] (6);

[0042] The equation of the adjoint wave field v is as follows:

[0043] (7);

[0044] In the formula, B is the frequency domain of the wave, s is the source term in the frequency domain, d is the observation data, P represents the projection matrix of the microseismic sensor position, is the operator for taking the real part of a complex number, and then the adjoint wave field v is iteratively calculated until the iterative calculation converges;

[0045] Step 5.6: When the number of iterations exceeds 500 or the amplitude order of the adjoint wave field v is less than 10 -3 i (m).

[0046] Step 6: Generate a wave velocity tomography image and substitute it into a wave velocity-ground pressure evaluation model;

[0047] The wave velocity tomography image is the final wave velocity image after the initial wave velocity model is iteratively calculated in step 5.5 until the adjoint wave field V converges;

[0048] Step 6.1: According to the geology, lithology, and depth of the monitoring area, select a wave velocity-ground pressure empirical formula;

[0049] The selection basis of the wave velocity-ground pressure empirical formula is:

[0050] When the rock RQD value is > 75% and the pressure range is 10 MPa - 30 MPa, select the Marcus formula;

[0051] When the rock RQD value is < 50% and the pressure range is < 10 MPa, select the Ericsson formula;

[0052] When the rock RQD value is < 50% and the pressure range is > 30 MPa, select the effective stress formula;

[0053] When it does not belong to the first three cases, select the Ericsson formula;

[0054] ​Step 6.2: according to the selection result of step 6.1, one of formulas (8), (9), (10), (11) is selected for calculation to generate the ground pressure result;

[0055] Marcus formula: V p =aσ 0.5 (8);

[0056] Ericsson formula: V p = V p0 +aσ (9);

[0057] Ericsson formula: V p = V p0 +aσ+ bσ 2 (10);

[0058] Effective stress formula: V p = V0+DF-B0exp(-kP) (11);

[0059] Wherein V p is the longitudinal wave velocity, V p0 is the propagation speed of longitudinal wave in the medium without pressure loading, exp() is the natural exponential function, a and b are constants related to the properties of rock, sigma is the ground pressure, V0 is the wave velocity without pressure, pore and fracture, D is the internal derivative of velocity to pressure in the linear elastic range, B0 is the initial velocity drop value caused by pore and fracture without pressure, k is the constant controlling the attenuation degree of velocity drop, and F is the effective pressure applied.

[0060] Step 7: generate a ground pressure distribution map;

[0061] The ground pressure distribution map is an image obtained by calculating the wave velocity tomography image by the formula selected in step 6.2, which converts the wave velocity distribution into the ground pressure distribution.

[0062] The beneficial effects produced by the above technical solution are:

[0063] The present application provides a kind of based on wave velocity tomography stratum structure and ground pressure distribution detection method, for the defects that prior art needs to be realized by excavation or drilling hard rock stratum, the present application can detect the internal structure and ground pressure distribution of stratum without damaging original hard rock stratum, realize the fine detection of stratum to meet the construction needs of current deep mining engineering. BRIEF DESCRIPTION OF DRAWINGS

[0064] Figure 1 It is the overall flow chart of the present application deep mining stratum structure and ground pressure distribution detection method;

[0065] Figure 2 It is the wave velocity model schematic diagram of geological profile;

[0066] Figure 3 Wave velocity uniform ray path node distribution diagram for the present application;

[0067] Figure 4 Ray tracing path profile diagram for the present application;

[0068] Figure 5 Stratum wave velocity tomographic image and technical diagram for the present application. DETAILED DESCRIPTION

[0069] The specific embodiments of the present application are described in further detail below in conjunction with the accompanying drawings and examples. The following examples are used to illustrate the present application, but are not used to limit the scope of the present application.

[0070] A stratum structure and ground pressure distribution detection method based on wave velocity tomography, as shown in Figure 1 includes the following steps:

[0071] Step 1: A dense linear geophone array is arranged at five-meter intervals in the monitoring area;

[0072] Step 2: Collect the aquifer distribution, lithology detection results, and borehole core map of the monitoring area;

[0073] Step 2.1: Establish an initial wave velocity model according to the collected lithology detection results of the monitoring area;

[0074] The initial wave velocity model is a geological profile velocity model at the exploration line, which is generated according to the experimental wave velocity of seismic waves in different media. Since geological exploration in mining engineering is generally in the form of exploration lines, the obtained aquifer distribution and lithology detection results are generally geological profile data of the exploration line. Therefore, the initial wave velocity model is a geological profile velocity model at the exploration line, with a basic form Figure 2 As shown in, the color is redder, the wave speed is higher, and the color is bluer, the wave speed is lower.

[0075] Step 2.2: According to the collected aquifer distribution, the initial wave velocity model generated in step 2.1 is corrected for the first time, and the wave speed in the aquifer area is lowered to 1500 m / s;

[0076] Step 2.3: Calculate the RQD value according to the collected borehole core map, and speculate the occurrence of stratum structure surface and fault of each exploration line. The velocity model corrected for the first time in step 2.2 is corrected for the second time, and the wave speed in the area containing structure surface and fault is lowered to 2800 m / s;

[0077] Step 3: Record the seismic wave data, record the period greater than or equal to 2 days, set the acquisition frequency to 10 Hz, and record the waveform data of 2 Hz-10 Hz.

[0078] Step 4: denoising the seismic wave data, picking up the arrival time, and associating with waveform events;

[0079] Step 4.1: denoising the detected seismic wave data, and picking up the arrival time of P wave and S wave;

[0080] Step 4.2: automatically matching the seismic wave data with the waveforms in the PEER waveform template library to generate waveform events;

[0081] Step 4.3: associating the waveform events generated in step 4.2 with the waveform arrival time picked up in step 4.1 in order.

[0082] Step 5: using the wave velocity inversion algorithm to track the path of the observation results;

[0083] Step 5.1: assuming that the seismic wave propagates along the path with the shortest propagation time, the seismic wave satisfies Fermat's principle:

[0084] (1);

[0085] In the formula, t is time, l is ray position, and v(r) is wave velocity.

[0086] In a continuously changing inhomogeneous medium, the wave velocity gradient will cause the ray path to bend. According to the extended form of Snell's law, the process of adjusting the ray direction with the velocity gradient satisfies the differential equation:

[0087] (2);

[0088] In the formula, v is the medium velocity, l is the ray position, s is the ray arc length, and ∇ is the gradient operator.

[0089] Substitute formula (2) into the wave velocity model after the second correction in step 2.3 to perform iterative correction.

[0090] Step 5.2: discretize the monitoring area into a grid to obtain the shortest path tree;

[0091] Starting from the source, first calculate the travel time of each node adjacent to the source point, then take the secondary node as a sub-wave source point, calculate the travel time of each node adjacent to the next level of the sub-wave source point from the sub-wave source point, and select the ray path with the minimum travel time to form a shortest path tree; the ray path and travel time from the source point to all receiving points are obtained from the position of the receiving point; the discretization grid node setting method is shown in Figure 3 , where the wave velocity is regarded as uniform within the grid, and the travel time is equal to the Euclidean distance between two adjacent nodes divided by the wave velocity within the grid.

[0092] Step 5.3: Path smoothing is performed on the shortest path tree obtained in step 5.2 to obtain the simulated ray path;

[0093] The simulated ray path result in this embodiment is as follows Figure 4 Different color lines represent different propagation paths, and the left starting point is the source of the vibration wave, and the right end point of each path is the detector;

[0094] Step 5.4: After obtaining the simulated ray path, a theoretical best matching model is defined using the full waveform inversion algorithm :

[0095] (3);

[0096] wherein, is the theoretical best matching model, i is the detector number, n src represents the number of sources, F i (m) represents the forward synthetic data, d i is the observation data; is the error function of the observation data and the synthetic data;

[0097] Full waveform inversion is essentially to find the minimum value of the residual, that is, an optimization problem. The wave field residual of the detector point is assumed to be an accompanied source, and it is considered that this accompanied source releases an accompanied wave field V; then the accompanied wave field exists in the time domain wave and the frequency domain wave;

[0098] The equation of the time domain wave is as follows:

[0099] (4);

[0100] wherein, m(x) is the time domain of the wave, x is a particle in the wave propagation path, s(x, t) and u(x, t) represent the source term and the wave field in the time domain respectively, and Δ represents the Laplace operator;

[0101] The equation of the frequency domain wave is as follows:

[0102] (5);

[0103] wherein, B(s, ω) is the frequency domain of the wave, ω is the angular frequency of the wave, s(x, ω) and u(x, ω) represent the source term and the wave field in the frequency domain respectively;

[0104] Step 5.5: The minimum value of the residual, that is, the minimum value of the accompanied wave field v, is solved;

[0105] The gradient vector g of the accompanied wave field is calculated:

[0106] (6);

[0107] The equation of the wave field v is as follows:

[0108] (7);

[0109] where B is the frequency domain of the wave, s is the source term in the frequency domain, d is the observation data, P represents the projection matrix of the microseismic sensor position, is the operator of taking the real part of a complex number, and then the wave field v is iteratively calculated until the iterative calculation converges;

[0110] Step 5.6: When the number of iterations exceeds 500 or the amplitude order of the wave field v is less than 10 -3 , the iteration stops and the forward synthetic data F i (m) is output.

[0111] Step 6: Generate a wave velocity tomography image and substitute it into a wave velocity-ground pressure evaluation model;

[0112] The wave velocity tomography image is the final wave velocity image after the initial wave velocity model is iteratively calculated in step 5.5 until the wave field V converges;

[0113] Step 6.1: According to the geology, lithology, and depth of the monitoring area, select a wave velocity-ground pressure empirical formula;

[0114] In this embodiment, the selection is as shown in Table 1;

[0115] Table 1. Selection of wave velocity-ground pressure empirical formula judgment basis:

[0116]

[0117] Step 6.2: According to the selection result of step 6.1, select one of formulas (8), (9), (10), and (11) to calculate and generate a ground pressure result;

[0118] Marcus formula: V p = aσ 0.5 (8);

[0119] Ericsson formula: V p = V p0 + aσ (9);

[0120] Elkson formula: V p = V p0 + aσ + bσ 2 (10);

[0121] Effective stress formula: V p = V0 + DF - B0exp(-kP) (11);

[0122] where V0 is the initial wave velocity, D is the wave velocity change coefficient, F is the effective stress, B0 is the wave velocity change coefficient, k is the wave velocity change coefficient, and P is the ground pressure.p Vp is the P-wave velocity, p0 Vp is the P-wave velocity, exp() is the natural exponential function, a, b are constants related to the properties of the rock, σ is the ground stress, V0 is the wave velocity without pressure, pore and fracture, D is the intrinsic derivative of velocity to pressure in the linear elastic range, B0 is the initial velocity drop value caused by pore and fracture without pressure, k is the constant that controls the attenuation degree of velocity drop, F is the effective pressure applied.

[0123] Step 7: generate the ground stress distribution map;

[0124] The ground stress distribution map is the image obtained by converting the wave velocity distribution into the ground stress distribution after the wave velocity tomography image is calculated by the formula selected in step 6.2.

[0125] As Figure 5 Fig. 1 is a wave velocity tomography image and a technical schematic diagram of a certain gold mine in Shandong area, and the blue part is a low-velocity abnormal area, which is suspected to be an underground goaf. After subsequent drilling verification, it is confirmed that the remaining goaf caused by incomplete filling leads to subsidence on the ground.

[0126] The above description is only the preferred embodiments of the present disclosure and the explanation of the applied technical principles. Those skilled in the art should understand that the scope of the invention involved in the embodiments of the present disclosure is not limited to the technical solutions formed by the specific combinations of the above technical features, and should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the above inventive concept. For example, the technical solutions formed by mutually replacing the above features and the technical features disclosed in the embodiments of the present disclosure (but not limited to) having similar functions.

Claims

1. A method for detecting formation structure and pressure distribution based on wave velocity tomography, characterized in that, The method comprises the following steps: Step 1: arranging dense linear detector arrays at intervals of five meters in the monitoring area; Step 2: collecting the aquifer distribution of the monitoring area, the lithology detection results and the borehole core map; Step 3: recording the seismic wave data, the recording period is greater than or equal to 2 days, the acquisition frequency is set to 10 Hz, and the waveform data of 2 Hz-10 Hz is recorded; Step 4: performing noise reduction processing on the seismic wave data, picking up the arrival time, and associating with the waveform event; Step 5: performing path tracking on the observation results by using the wave velocity inversion algorithm; Step 6: generating a wave velocity tomography image and substituting it into a wave velocity-ground pressure evaluation model; The wave velocity tomography image is a final wave velocity image obtained after the initial wave velocity model is iteratively calculated until the accompanying wave field V converges; Step 7: generating a ground pressure distribution map.

2. The method according to claim 1, wherein, The step 2 comprises the following steps: Step 2.1: establishing an initial wave velocity model according to the collected lithology detection results of the monitoring area; The initial wave velocity model is a geological profile velocity model at the exploration line generated according to the experimental wave velocity of seismic waves in different media; Step 2.2: performing first correction on the initial wave velocity model generated in step 2.1 according to the collected aquifer distribution, and lowering the wave velocity in the aquifer area to 1500 m / s; Step 2.3: calculating the RQD value according to the collected borehole core map, inferring the occurrence of the stratigraphic structure surface and fault of each exploration line, and performing second correction on the velocity model after the first correction in step 2.2, and lowering the wave velocity in the area containing the structure surface and fault to 2800 m / s.

3. The method according to claim 1, wherein, The step 4 comprises the following steps: Step 4.1: picking up the P-wave arrival time and S-wave arrival time after noise reduction processing on the detected seismic wave data; Step 4.2: automatically matching the seismic wave data with the waveforms in the PEER waveform template library to generate a waveform event; Step 4.3: associating the waveform event generated in step 4.2 with the waveform arrival time picked up in step 4.1 in sequence.

4. The method according to claim 1, wherein, The step 5 comprises the following steps: Step 5.1: assuming that the seismic wave propagates along the path with the shortest propagation time, the seismic wave satisfies Fermat's principle: （1）； In the formula, t is time, l is ray position, and v(r) is wave velocity; According to the extended form of Snell's law, the process of adjusting the ray direction with the velocity gradient satisfies the differential equation: （2）； In the formula, v is the medium velocity, l is the ray position, s is the ray arc length, and ∇ is the gradient operator; Substitute formula (2) into the wave velocity model after the second correction in step 2.3 for iterative correction; Step 5.2: discretizing the grid of the monitoring area to obtain a shortest path tree; Starting from the source, first calculate the travel time of each node adjacent to the source point, then take the secondary node as a sub-wave source point, calculate the travel time of each node adjacent to the sub-wave source point at the next level of the sub-wave source point, and select the ray path with the minimum travel time to form a shortest path tree; the ray path and travel time from the source point to all receiving points are obtained from the position of the receiving point; wherein the wave velocity is regarded as uniform in the grid, and the travel time is equal to the Euclidean distance between two adjacent nodes divided by the wave velocity in the grid; Step 5.3: Path smoothing is performed on the shortest path tree obtained in step 5.2 to obtain the simulated ray path; Step 5.4: After obtaining the simulated ray paths, a full waveform inversion algorithm is used to define a theoretical best match model : （3）； where, is the best matching model, i is the receiver number, n src represents the number of sources, F i (m) represents the forward synthetic data, d i is the observed data; is the error function of the observed data and the synthetic data; The wave field residual of the detector point is assumed to be an accompanied source, and the accompanied source is considered to release the accompanied wave field V; the accompanied wave field exists in the time domain and the frequency domain; The equation of the time domain wave is as follows: （4）； In the formula, m(x) is the time domain of the wave, x is a particle in the wave propagation path, s(x, t) and u(x, t) represent the source term and the wave field in the time domain respectively, and Δ represents the Laplace operator; The equation of the frequency domain wave is as follows: （5）； In the formula, B(s, ω) is the frequency domain of the wave, ω is the angular frequency of the wave, s(x, ω) and u(x, ω) represent the source term and the wave field in the frequency domain respectively; Step 5.5: The minimum value of the residual, that is, the minimum value of the accompanied wave field v, is solved; The gradient vector g of the accompanied wave field v is calculated as follows: （6）； The equation of the accompanied wave field v is as follows: （7）； where B is the frequency domain of the wave, s is the source term in the frequency domain, d is the observed data, P represents the projection matrix of the microseismic sensor position, is the operator that takes the real part of a complex number, and then iteratively computes the adjoint wavefield v until the iterative computation converges. Step 5.6: When the number of iterations exceeds 500 or the amplitude order of the accompanying wavefield v is less than 10 -3 , the iteration stops and the forward synthetic data F i (m) is output.

5. The method according to claim 1, wherein, The step 6 comprises the following steps: Step 6.1: According to the geology, lithology and depth of the monitoring area, an empirical formula of wave velocity-pressure is selected; Step 6.2: According to the selection result of step 6.1, the empirical formula of wave velocity-pressure is selected for calculation to generate the pressure result.

6. The method according to claim 5, wherein, The selection basis of the empirical formula of wave velocity-pressure in step 6.1 is as follows: When the rock RQD value is greater than 75% and the pressure range is 10-30 MPa, the Marcus formula is selected; When the rock RQD value is less than 50% and the pressure range is less than 10 MPa, the Ericsson formula is selected; When the rock RQD value is less than 50% and the pressure range is greater than 30 MPa, the effective stress formula is selected; When the above three conditions are not met, the Ericsson formula is selected.

7. The method according to claim 6, wherein, The empirical formula of wave velocity-pressure is as follows: Marcus formula: V p = aσ 0.5 ; Ericsson formula: V p = V p0 + aσ; Elkson formula: V p = V p0 + aσ + bσ 2 ; Effective stress formula: V p = V0+DF-B0exp(-kP); where V p is the P-wave velocity, V p0 is the propagation velocity of the P-wave in the medium without pressure loading, exp() is the natural exponential function, a and b are constants related to the properties of the rock, σ is the ground pressure, V0 is the wave velocity without pressure, porosity and fractures, D is the intrinsic derivative of the velocity with respect to pressure in the linear elastic range, B0 is the initial velocity reduction value caused by porosity and fractures without pressure, k is a constant that controls the degree of attenuation of the velocity reduction, and F is the applied effective pressure.

8. The method according to claim 5, wherein, The pressure distribution map in step 7 is the image obtained by calculating the wave velocity tomography map by the empirical formula of wave velocity-pressure, which converts the wave velocity distribution into the pressure distribution.

Images