A two-dimensional numerical simulation method for fault detection in geological prediction

Abstract

The present invention belongs to the technical field of tunnel geophysical detection, and provides a two-dimensional numerical simulation method for faults in geological prediction detection. The method uses an unequally spaced grid difference method to perform numerical calculations on a model, adopts a PML absorbing boundary, continuously iterates velocity and stress for numerical simulation, and outputs seismic wave fields and seismic records. The method can effectively improve the accuracy of numerical simulations for advanced geological prediction detection of tunnels.

Description

Technical Field

[0001] The invention belongs to the technical field of tunnel geophysical detection, and in particular relates to a two-dimensional numerical simulation method for geological forecast detection faults. Background Art

[0002] It is of great significance to conduct advanced geological prediction before tunnel excavation, and it is necessary to study the advanced geological prediction of tunnels. During the construction process, the shape of the surrounding rock around the tunnel is undulating, and the traditional finite difference numerical simulation of tunnel advanced geological prediction detection treats the tunnel shape as a regular shape, such as a cylinder, rectangle, etc., which is inconsistent with the actual conditions and directly affects the results of numerical simulation. At present, there are limited numerical simulation methods for tunnel undulation conditions. Therefore, it is urgent to develop numerical simulation technology methods for tunnel undulation earthquake advanced geological prediction detection to further improve the accuracy of numerical simulation. Summary of the invention

[0003] In order to overcome the above-mentioned shortcomings of the prior art, the purpose of the present invention is to provide a two-dimensional numerical simulation method for geological prediction and fault detection, which can effectively improve the accuracy of numerical simulation of advanced geological prediction and detection of tunnels.

[0004] The technical solution adopted by the present invention to solve the technical problem is:

[0005] A two-dimensional numerical simulation method for geological prediction and fault detection comprises the following steps:

[0006] Step 1: Measure the ups and downs of the tunnel in the opposite direction of the tunnel face and build an ups and downs tunnel model;

[0007] Step 2: Establish a three-dimensional model of uniform medium velocity, with velocity parameter v p1 ,v s1 ,ρ1, where v p1 ,v s1 ,ρ1 are the longitudinal wave velocity, shear wave velocity and density of the three-dimensional model respectively;

[0008] Step 3: Insert the tunnel elevation parameter information into the velocity model to establish the velocity model of the undulating tunnel. The velocity parameter of the tunnel model is v p2 ,v s2 ,ρ2, where v p2 ,v s2 ,ρ2 are the longitudinal wave velocity, shear wave velocity and density of the tunnel space respectively;

[0009] Step 4: Set a fault with a strike of 90 degrees in the X direction, and the velocity parameter of the fault part is v p3 ,v s3 ,ρ3, where v p3 ,v s3,ρ3 are the P-wave velocity, S-wave velocity and density of the fault part respectively;

[0010] Step 5: Establish a two-dimensional coordinate system along the tunnel, with the X direction pointing to the tunnel face and the Y direction perpendicular to the tunnel sidewall. Establish a two-dimensional coordinate system with the center of the tunnel face on the model boundary as the origin, and establish an advanced geological prediction observation system for tunnel excavation;

[0011] Step 6: Given parameters such as grid spacing and time step, sampling law and stability conditions need to be satisfied;

[0012] Step 7: Calculate the distance between the real terrain point and the grid point on the tunnel surface through the real elevation;

[0013] Step 8: Use unequal distance finite difference scheme at the tunnel and conventional finite difference scheme at the rest;

[0014] Step 9: Arrange the detectors on the side wall of the undulating tunnel. Evenly distribute n three-component detectors and m seismic sources in sequence from the face to the tunnel. The number of the three-component detector is set to R n , the source number is S m , build a basement earthquake advanced geological prediction observation system;

[0015] Step 10: Use Ricker wavelet, source frequency is Z Hz, sampling interval is l, sampling time is t;

[0016] Step 11: Use conventional PML absorbing boundary processing;

[0017] Step 12: Iterate velocity and stress, and output wavefield snapshots and seismic records.

[0018] Preferably, in step 1:

[0019] (1) According to the tunnel conditions, measure the elevation parameters of the entire tunnel within a range of 90m from the tunnel face to the tunnel direction;

[0020] (2) Convert the elevation parameters into relative elevations to construct the basic form of the entire tunnel;

[0021] (3) Construct an undulating tunnel model based on the basic form of the tunnel.

[0022] Preferably, in the step 2: in the three-dimensional model, the lengths in the X direction and the Y direction are 300m and 100m respectively;

[0023] In step three: p2 =340m / s,v s2 =0m / s,ρ2=1.29kg / m 2 .

[0024] Preferably, in step six:

[0025] (1) Using the wave equation to grid the grid, the size is h;

[0026] (2) Set the numerical simulation parameters: frequency is f, unit is Hz, sampling interval is l, unit is ms;

[0027] (3) To avoid numerical dispersion, the sampling theorem and stability conditions must be met, specifically the grid size and time step:

[0028] Ⅰ. The grid spacing size meets the following requirements:

[0029] Among them, v min / (2f c n) represents the minimum wavelength of the simulation, v min is the minimum velocity parameter among the numerical simulation parameters, and n is the correlation factor of the difference order taken;

[0030] II. The time step size satisfies:

[0031] Where m is a factor related to the order of difference taken, v max is the minimum speed parameter among the numerical simulation parameters.

[0032] Preferably, in step eight:

[0033] (1) After meshing the model, the conventional equidistant finite difference method is used to process the periphery of the undulating tunnel sidewall;

[0034] (2) Unequally spaced finite differences are performed near the undulating tunnel, which are mainly divided into unequally spaced finite differences of stress and velocity, as follows:

[0035] Ⅰ. Through U'=R β U rotates the undulating tunnel coordinates to horizontal conditions, where U = (v x ,v y ,σ xx ,σ xy ,σ yy ) T ,U'=(v' x ,v' y ,σ' xx ,σ' xy ,σ' yy ) T , v x ,v y Respectively represent the velocity components, σ xx ,σ xy ,σ yy is the stress component, v'x ,v' y is the converted velocity component, σ' xx ,σ' xy ,σ' yy is the converted stress component, R β is the coordinate rotation matrix, β is the rotation angle, specifically:

[0036]

[0037] Ⅱ. When calculating stress, there are 8 cases of unequally spaced finite differences:

[0038] ① Under undulating conditions σ xx , σ xy , σ yy When there are four grid points around the tunnel undulation position, there is no unequally spaced finite difference, then σ xx , σ xy , σ yy The expression is as follows:

[0039]

[0040] Among them, Δl and Δh are the sampling interval and grid size respectively, i and j represent the number of grid points in the X and Y directions respectively, λ and μ are Lamé parameters, and the sampling interval is l;

[0041] ②Under undulating conditions σ xx , σ yy There are four grid points around the tunnel undulation position, but σ xy When there are only three grid points around the tunnel undulation position, calculate σ xx , σ yy The conventional equidistant finite difference scheme is used for calculation; while the calculation of σ xy When there is an unequally spaced finite difference calculation, then σ xx , σ xy , σ yy The expression is as follows:

[0042]

[0043] Among them, y0(i) is the coordinate parameter of the tunnel undulation, and the following y1(i) is the tunnel coordinate rounded down;

[0044] ③ Under undulating conditions σ xx , σ yy There are three grid points around the tunnel undulation position, σ xy When there are only two grid points around the tunnel undulation position, the difference format is an unequally spaced finite difference format, then σ xx , σ xy , σyy The expression is as follows:

[0045]

[0046] ④ Under undulating conditions σ xx , σ yy There are three grid points around the tunnel undulation position. xy When there is only one grid point around the tunnel undulation position, the difference format is changed to the unequally spaced finite difference format, and only σ xx , σ yy , the expression is as follows:

[0047]

[0048] ⑤Under undulating conditionsσ xx , σ yy There are three grid points around the tunnel undulation position. xy When there are only four grid points around the tunnel undulation position, calculate σ xx , σ yy The calculation is performed using the unequally spaced finite difference format; while the calculation of σ xy Using the conventional difference method, σ xx , σ xy , σ yy The expression is as follows:

[0049]

[0050] ⑥Under undulating conditionsσ xx , σ yy There are two grid points around the tunnel undulation position. xy When there are only three grid points around the tunnel undulation position, calculate σ xx , σ xy , σ yy Using the unequal-distance finite difference method, σ xx , σ xy , σ yy The expression is as follows:

[0051]

[0052] ⑦Under undulating conditionsσ xx , σ yy There are two grid points around the tunnel undulation position. xy When there are only two grid points around the tunnel undulation position, calculate σ xx , σ xy , σ yy Using the unequal-distance finite difference method, σ xx, σ xy , σ yy The expression is as follows:

[0053]

[0054] ⑧Under fluctuating conditions σ xx , σ yy There are two grid points around the tunnel undulation position. xy When there is only one grid point around the tunnel undulation position, only σ xx , σ yy Using the unequal-distance finite difference method, σ xx , σ yy The expression is as follows:

[0055]

[0056] III. When calculating speed, there are 5 different expressions in the horizontal or vertical direction, specifically:

[0057] ① In calculating v x and v y When , both horizontal and vertical directions can be calculated normally, and the conventional equidistant finite difference format is used for calculation, as follows:

[0058]

[0059] ② When calculating v y The conventional equidistant finite difference scheme is used when calculating v x When there is unequal spacing in the vertical direction, the unequal spacing finite difference format is used instead, as follows:

[0060]

[0061] ③ In calculating v y The conventional equidistant finite difference scheme is used when calculating v x When there are unequal spacings in both the horizontal and vertical directions, the unequal spacing finite difference format is used instead, as follows:

[0062]

[0063] ④ In calculating v y The conventional equidistant finite difference scheme is used when v x No calculation is required when outside the surface, as follows:

[0064]

[0065] ⑤Calculate v yWhen there is unequal spacing in the horizontal direction, the unequal spacing finite difference format is used instead; and v x No calculation is required when outside the surface, as follows:

[0066]

[0067] IV. The undulation of the other side wall of the tunnel is handled according to the above steps I-III;

[0068] Ⅴ. The undulating conditions at the tunnel face shall be handled according to the above steps Ⅰ-Ⅲ.

[0069] Preferably, in step nine: the interval between adjacent three-component geophones is three meters, the tunnel face is provided with a first hole position for setting the three-component geophone, the hole depth of the first hole position is two meters, and the three-component geophone is arranged at the center position of the side wall of the first hole position.

[0070] Preferably, in step nine: the number of the three-component detectors is no more than 31.

[0071] Preferably, in step nine: the interval between adjacent seismic sources is five meters, the tunnel face is provided with a second hole position for setting the seismic source, the hole depth of the second hole position is three meters, and the seismic source is arranged at the center of the side wall of the second hole position.

[0072] Preferably, in step nine: the number of seismic sources is not greater than 17.

[0073] Preferably, in step twelve:

[0074] (1) This paper adopts the first-order velocity stress two-dimensional elastic wave equation of isotropic medium:

[0075]

[0076] In the above formula, σ xx , σ yy are the normal stresses along the X and Y directions, λ and μ are the Lame parameters, ρ is the density, and v x ,v y are the vibration speeds of the particle along the X and Y directions respectively, and t is the time;

[0077] (2) Solve using high-order differences under staggered grids. The specific difference method is shown in step 8.

[0078] (3) Continuously iterate velocity and stress, and output wave field snapshots and seismic record files.

[0079] Compared with the prior art, the beneficial effects of the present invention include:

[0080] This scheme is based on the advanced geological prediction detection of tunnel fluctuation earthquakes using unequally spaced finite differences. Modeling is carried out through the actual measured tunnel parameters, and velocity values ​​are assigned to different grid point positions. According to the different undulating boundaries of the tunnel, the model is numerically calculated using the unequally spaced grid difference method. The PML absorbing boundary is used, and the velocity and stress are continuously iterated for numerical simulation. The seismic wave field and seismic records are output. This method can effectively improve the accuracy of numerical simulation of tunnel advanced geological prediction detection. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0082] Figure 1 This is the scenario in which the solution of the present invention is applied.

[0083] Figure 2 It is the structure of the mesh;

[0084] Figure 3 It is a schematic diagram of unequal-distance finite-difference stress expression;

[0085] Figure 4 Schematic diagram of unequally spaced finite difference velocity expression.

[0086] in:

[0087] 1-undulating tunnel, 2-fault structure, 3-difference operator, 4-difference operator of undulating tunnel, 5-unequally spaced finite difference center points, 6-equidistant difference center points, 7-undulating surface of tunnel, 8-grid spacing in X direction, 9-grid spacing in Z direction, 10-elevation Z0. DETAILED DESCRIPTION

[0088] In order to more clearly understand the above-mentioned objects, features and advantages of the present invention, the present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that, in the absence of conflict, the embodiments of the present application and the features in the embodiments can be combined with each other. In the following description, many specific details are set forth in order to fully understand the present invention, and the embodiments described are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in the art without making creative work are within the scope of protection of the present invention.

[0089] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those commonly understood by those skilled in the art of the present invention. The terms used in the specification of the present invention herein are only for the purpose of describing specific embodiments and are not intended to limit the present invention.

[0090] Example:

[0091] like Figure 1-4 As shown, this embodiment provides a two-dimensional numerical simulation method for geological prediction fault detection, including the following steps:

[0092] Step 1: Measure the ups and downs of the tunnel in the opposite direction of the tunnel face and build an ups and downs tunnel model;

[0093] Step 2: Establish a three-dimensional model of uniform medium velocity, with velocity parameter v p1 ,v s1 ,ρ1, where v p1 ,v s1 ,ρ1 are the longitudinal wave velocity, shear wave velocity and density of the three-dimensional model respectively;

[0094] Step 3: Insert the tunnel elevation parameter information into the velocity model to establish the velocity model of the undulating tunnel. The velocity parameter of the tunnel model is v p2 ,v s2 ,ρ2, where v p2 ,v s2 ,ρ2 are the longitudinal wave velocity, shear wave velocity and density of the tunnel space respectively;

[0095] Step 4: Set a fault with a strike of 90 degrees in the X direction, and the velocity parameter of the fault part is v p3 ,v s3 ,ρ3, where v p3 ,v s3 ,ρ3 are the P-wave velocity, S-wave velocity and density of the fault part respectively;

[0096] Step 5: Establish a two-dimensional coordinate system along the tunnel, with the X direction pointing to the tunnel face and the Y direction perpendicular to the tunnel sidewall. Establish a two-dimensional coordinate system with the center of the tunnel face on the model boundary as the origin, and establish an advanced geological prediction observation system for tunnel excavation;

[0097] Step 6: Given parameters such as grid spacing and time step, sampling law and stability conditions need to be satisfied;

[0098] Step 7: Calculate the distance between the real terrain point and the grid point on the tunnel surface through the real elevation;

[0099] Step 8: Use unequal distance finite difference scheme at the tunnel and conventional finite difference scheme at the rest;

[0100] Step 9: Arrange the detectors on the side wall of the undulating tunnel. Evenly distribute n three-component detectors and m seismic sources in sequence from the face to the tunnel. The number of the three-component detector is set to R n , the source number is S m , build a basement earthquake advanced geological prediction observation system;

[0101] Step 10: Use Ricker wavelet, source frequency is Z Hz, sampling interval is l, sampling time is t;

[0102] Step 11: Use conventional PML absorbing boundary processing;

[0103] Step 12: Iterate velocity and stress, and output wavefield snapshots and seismic records.

[0104] Specifically, in step one:

[0105] (1) According to the tunnel conditions, measure the elevation parameters of the entire tunnel within a range of 90m from the tunnel face to the tunnel direction;

[0106] (2) Convert the elevation parameters into relative elevations to construct the basic form of the entire tunnel;

[0107] (3) Construct an undulating tunnel model based on the basic form of the tunnel.

[0108] Specifically, in step 2: in the three-dimensional model, the lengths in the X and Y directions are 300m and 100m respectively;

[0109] In step three: p2 =340m / s,v s2 =0m / s,ρ2=1.29kg / m 2 .

[0110] Specifically, in step six:

[0111] (1) Using the wave equation to grid the grid, the size is h;

[0112] (2) Set the numerical simulation parameters: frequency is f, unit is Hz, sampling interval is l, unit is ms;

[0113] (3) To avoid numerical dispersion, the sampling theorem and stability conditions must be met, specifically the grid size and time step:

[0114] Ⅰ. The grid spacing size meets the following requirements:

[0115] Among them, v min / (2f c n) represents the minimum wavelength of the simulation, v minis the minimum velocity parameter among the numerical simulation parameters, and n is the correlation factor of the difference order taken;

[0116] II. The time step size satisfies:

[0117] Where m is a factor related to the order of difference taken, v max is the minimum speed parameter among the numerical simulation parameters.

[0118] Specifically, in step eight:

[0119] (1) After meshing the model, the conventional equidistant finite difference method is used to process the periphery of the undulating tunnel sidewall;

[0120] (2) Unequally spaced finite differences are performed near the undulating tunnel, which are mainly divided into unequally spaced finite differences of stress and velocity, as follows:

[0121] Ⅰ. Through U'=R β U rotates the undulating tunnel coordinates to horizontal conditions, where U = (v x ,v y ,σ xx ,σ xy ,σ yy ) T ,U'=(v' x ,v' y ,σ' xx ,σ' xy ,σ' yy ) T , v x ,v y Respectively represent the velocity components, σ xx ,σ xy ,σ yy is the stress component, v' x ,v' y is the converted velocity component, σ' xx ,σ' xy ,σ' yy is the converted stress component, R β is the coordinate rotation matrix, β is the rotation angle, specifically:

[0122]

[0123] Ⅱ. There are 8 cases of unequally spaced finite differences when calculating stress, see Figure 3 :

[0124] ① Under undulating conditions σ xx , σ xy , σ yyWhen there are four grid points around the tunnel undulation position, there is no unequally spaced finite difference, then σ xx , σ xy , σ yy The expression is as follows:

[0125]

[0126] Where Δl and Δh are the sampling interval and grid size, respectively, i and j represent the number of grid points in the X and Y directions, respectively, and λ and μ are the Lamé parameters;

[0127] ②Under undulating conditions σ xx , σ yy There are four grid points around the tunnel undulation position, but σ xy When there are only three grid points around the tunnel undulation position, calculate σ xx , σ yy The conventional equidistant finite difference scheme is used for calculation; while the calculation of σ xy When there is an unequally spaced finite difference calculation, then σ xx , σ xy , σ yy The expression is as follows:

[0128]

[0129] Among them, y0(i) is the coordinate parameter of the tunnel undulation, and the following y1(i) is the tunnel coordinate rounded down;

[0130] ③ Under undulating conditions σ xx , σ yy There are three grid points around the tunnel undulation position, σ xy When there are only two grid points around the tunnel undulation position, the difference format is an unequally spaced finite difference format, then σ xx , σ xy , σ yy The expression is as follows:

[0131]

[0132]

[0133] ④ Under undulating conditions σ xx , σ yy There are three grid points around the tunnel undulation position. xy When there is only one grid point around the tunnel undulation position, the difference format is changed to the unequally spaced finite difference format, and only σ xx , σ yy , the expression is as follows:

[0134]

[0135] ⑤Under undulating conditionsσ xx , σ yy There are three grid points around the tunnel undulation position. xy When there are only four grid points around the tunnel undulation position, calculate σ xx , σ yy The calculation is performed using the unequally spaced finite difference format; while the calculation of σ xy Using the conventional difference method, σ xx , σ xy , σ yy The expression is as follows:

[0136]

[0137] ⑥Under undulating conditionsσ xx , σ yy There are two grid points around the tunnel undulation position. xy When there are only three grid points around the tunnel undulation position, calculate σ xx , σ xy , σ yy Using the unequal-distance finite difference method, σ xx , σ xy , σ yy The expression is as follows:

[0138]

[0139] ⑦Under undulating conditionsσ xx , σ yy There are two grid points around the tunnel undulation position. xy When there are only two grid points around the tunnel undulation position, calculate σ xx , σ xy , σ yy Using the unequal-distance finite difference method, σ xx , σ xy , σ yy The expression is as follows:

[0140]

[0141] ⑧Under fluctuating conditions σ xx , σ yy There are two grid points around the tunnel undulation position. xy When there is only one grid point around the tunnel undulation position, only σ xx , σ yy Using the unequal-distance finite difference method, σ xx , σyy The expression is as follows:

[0142]

[0143] III. When calculating speed, there are 5 different expressions in the horizontal or vertical direction, see Figure 4 , specifically:

[0144] ① In calculating v x and v y When , both horizontal and vertical directions can be calculated normally, and the conventional equidistant finite difference format is used for calculation, as follows:

[0145]

[0146] (4-10);

[0147] ② When calculating v y The conventional equidistant finite difference scheme is used when calculating v x When there is unequal spacing in the vertical direction, the unequal spacing finite difference format is used instead, as follows:

[0148]

[0149] ③ In calculating v y The conventional equidistant finite difference scheme is used when calculating v x When there are unequal spacings in both the horizontal and vertical directions, the unequal spacing finite difference format is used instead, as follows:

[0150]

[0151] ④ When calculating v y The conventional equidistant finite difference scheme is used when v x No calculation is required when outside the surface, as follows:

[0152]

[0153] ⑤Calculate v y When there is unequal spacing in the horizontal direction, the unequal spacing finite difference format is used instead; and v x No calculation is required when outside the surface, as follows:

[0154]

[0155] IV. The undulation of the other side wall of the tunnel is handled according to the above steps I-III;

[0156] Ⅴ. The undulating conditions at the tunnel face shall be handled according to the above steps Ⅰ-Ⅲ.

[0157] Specifically, in step nine: the interval between adjacent three-component geophones is three meters, the tunnel face is provided with a first hole position for setting the three-component geophone, the hole depth of the first hole position is two meters, and the three-component geophone is arranged at the center of the side wall of the first hole position.

[0158] Specifically, in step nine: the number of three-component detectors is not greater than 31.

[0159] Specifically, in step nine: the interval between adjacent seismic sources is five meters, a second hole position for setting the seismic source is provided on the tunnel face, the hole depth of the second hole position is three meters, and the seismic source is arranged at the center of the side wall of the second hole position.

[0160] Specifically, in step nine: the number of earthquake sources is not greater than 17.

[0161] Specifically, in step 12:

[0162] (1) This paper adopts the first-order velocity stress two-dimensional elastic wave equation of isotropic medium:

[0163]

[0164] In the above formula, σ xx , σ yy are the normal stresses along the X and Y directions, λ and μ are the Lame parameters, ρ is the density, and v x ,v y are the vibration speeds of the particle along the X and Y directions respectively, and t is the time;

[0165] (2) Solve using high-order differences under staggered grids. The specific difference method is shown in step 8.

[0166] (3) Continuously iterate velocity and stress, and output wave field snapshots and seismic record files.

[0167] Figure 1 The scenario for application of this scheme includes undulating tunnel 1 and fault structure 2;

[0168] Figure 2 is the structure of the grid, specifically, the unequally spaced finite difference grid, which includes the difference operator 3 (v x ,v z ), the difference operator of the undulating tunnel 4(v x ,v z ), unequally spaced finite difference center points 5, equally spaced difference center points 6, tunnel undulating surface 7, X-direction grid spacing 8 and Z-direction grid spacing 9, the elevation Z0 of the Z-direction grid is labeled 10.

[0169] The damage to the surrounding rock during tunnel excavation is uneven. When conducting earthquake advance geological prediction detection, the detectors are arranged on the uneven tunnel walls. When the earthquake source is excited on the tunnel side wall, the complex tunnel conditions have an impact on the seismic wave field that cannot be ignored. Unlike the undulating surface conditions, ground three-dimensional seismic is a semi-space detection, and tunnel advance geological prediction detection is carried out under full-space conditions, with undulating tunnel space side walls. The simulation method proposed in the present invention models the tunnel parameters based on actual measurements, assigns velocities to different grid point positions, and numerically calculates the model using the unequally spaced grid difference method according to the different undulating boundaries of the tunnel. The PML absorbing boundary is used to continuously iterate the velocity and stress for numerical simulation, and output the seismic wave field and seismic records. This method can effectively improve the accuracy of the numerical simulation of tunnel advance geological prediction detection.

[0170] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Therefore, any modification, equivalent change and modification made to the above embodiment according to the technical essence of the present invention without departing from the technical solution of the present invention shall still fall within the scope of the technical solution of the present invention.

Claims

1. A two-dimensional numerical simulation method for geological prediction and fault detection, characterized in that: The following steps are involved: Step 1: Measure the ups and downs of the tunnel in the opposite direction of the tunnel face and build an ups and downs tunnel model; Step 2: Establish a three-dimensional model of uniform medium velocity, with velocity parameter v p1 ,v s1 ,ρ1, where v p1 ,v s1 ,ρ1 are the longitudinal wave velocity, shear wave velocity and density of the three-dimensional model respectively; Step 3: Insert the tunnel elevation parameter information into the velocity model to establish the velocity model of the undulating tunnel. The velocity parameter of the tunnel model is v p2 ,v s2 ,ρ2, where v p2 ,v s2 ,ρ2 are the longitudinal wave velocity, shear wave velocity and density of the tunnel space respectively; Step 4: Set a fault with a strike of 90 degrees in the X direction, and the velocity parameter of the fault part is v p3 ,v s3 ,ρ3, where v p3 ,v s3 ,ρ3 are the P-wave velocity, S-wave velocity and density of the fault part respectively; Step 5: Establish a two-dimensional coordinate system along the tunnel, with the X direction pointing to the tunnel face and the Y direction perpendicular to the tunnel sidewall. Establish a two-dimensional coordinate system with the center of the tunnel face on the model boundary as the origin, and establish an advanced geological prediction observation system for tunnel excavation; Step 6: Given the grid spacing and time step parameters, the sampling law and stability conditions need to be satisfied; Step 7: Calculate the distance between the real terrain point and the grid point on the tunnel surface through the real elevation; Step 8: Use unequal distance finite difference scheme at the tunnel and conventional equidistant finite difference scheme at the rest; (1) After meshing the model, the conventional equidistant finite difference method is used to process the periphery of the undulating tunnel sidewall; (2) Unequally spaced finite differences are performed near the undulating tunnel, which are divided into unequally spaced finite differences of stress and velocity, as follows: Ⅰ. Through U'=R β U rotates the undulating tunnel coordinates to horizontal conditions, where U = (v x ,v y ,σ xx ,σ xy ,σ yy ) T ,U'=(v' x ,v' y ,σ' xx ,σ' xy ,σ' yy ) T , v x ,v y Respectively represent the velocity components, σ xx ,σ xy ,σ yy is the stress component, v' x ,v' y is the converted velocity component, σ' xx ,σ' xy ,σ' yy is the converted stress component, R β is the coordinate rotation matrix, β is the rotation angle, specifically: Ⅱ. When calculating stress, there are 8 cases of unequally spaced finite differences: ① Under undulating conditions σ xx , σ xy , σ yy When there are four grid points around the tunnel undulation position, there is no unequally spaced finite difference, then σ xx , σ xy , σ yy The expression is as follows: Among them, Δl and Δh are the sampling interval and grid size respectively, i and j represent the number of grid points in the X and Y directions respectively, λ and μ are Lamé parameters, and the sampling interval is l; ②Under undulating conditions σ xx , σ yy There are four grid points around the tunnel undulation position, but σ xy When there are only three grid points around the tunnel undulation position, calculate σ xx , σ yy The conventional equidistant finite difference scheme is used for calculation; while the calculation of σ xy When there is an unequally spaced finite difference calculation, then σ xx , σ xy , σ yy The expression is as follows: Among them, y0(i) is the coordinate parameter of the tunnel undulation, and the following y1(i) is the tunnel coordinate rounded down; ③ Under undulating conditions σ xx , σ yy There are three grid points around the tunnel undulation position, σ xy When there are only two grid points around the tunnel undulation position, the difference format is the unequal distance finite difference format, h is the grid spacing, then σ xx , σ xy , σ yy The expression is as follows: ④ Under undulating conditions σ xx , σ yy There are three grid points around the tunnel undulation position. xy When there is only one grid point around the tunnel undulation position, the difference format is changed to the unequally spaced finite difference format, and only σ xx , σ yy , the expression is as follows: ⑤Under undulating conditionsσ xx , σ yy There are three grid points around the tunnel undulation position. xy When there are only four grid points around the tunnel undulation position, calculate σ xx , σ yy The calculation is performed using the unequally spaced finite difference format; while the calculation of σ xy Using the conventional equidistance finite difference method, σ xx , σ xy , σ yy The expression is as follows: ⑥Under undulating conditionsσ xx , σ yy There are two grid points around the tunnel undulation position. xy When there are only three grid points around the tunnel undulation position, calculate σ xx , σ xy , σ yy Using the unequal-distance finite difference method, σ xx , σ xy , σ yy The expression is as follows: ⑦Under undulating conditionsσ xx , σ yy There are two grid points around the tunnel undulation position. xy When there are only two grid points around the tunnel undulation position, calculate σ xx , σ xy , σ yy Using the unequal-distance finite difference method, σ xx , σ xy , σ yy The expression is as follows: ⑧Under fluctuating conditions σ xx , σ yy There are two grid points around the tunnel undulation position. xy When there is only one grid point around the tunnel undulation position, only σ xx , σ yy Using the unequal-distance finite difference method, σ xx , σ yy The expression is as follows: III. When calculating speed, there are 5 different expressions in the horizontal or vertical direction, specifically: ① In calculating v x and v y When , the conventional equidistant finite difference format is used for calculation in the horizontal and vertical directions, as follows: ρ in the above formula is the density; ② When calculating v y The conventional equidistant finite difference scheme is used when calculating v x When there is unequal spacing in the vertical direction, the unequal spacing finite difference format is used instead, as follows: ③ In calculating v y The conventional equidistant finite difference scheme is used when calculating v x When there are unequal spacings in both the horizontal and vertical directions, the unequal spacing finite difference format is used instead, as follows: ④ In calculating v y The conventional equidistant finite difference scheme is used when v x No calculation is required when outside the surface, as follows: ⑤Calculate v y When there is unequal spacing in the horizontal direction, the unequal spacing finite difference format is used instead; and v x No calculation is required when outside the surface, as follows: IV. The undulation of the other side wall of the tunnel is handled according to the above steps I-III; Ⅴ. The undulations at the tunnel face are handled according to the above steps Ⅰ-Ⅲ Step 9: Arrange the detectors on the side wall of the undulating tunnel. Evenly distribute n three-component detectors and m seismic sources in sequence from the face to the tunnel. The number of the three-component detector is set to R n , the source number is S m , build a basement earthquake advanced geological prediction observation system; Step 10: Use Ricker wavelet, source frequency is Z Hz, sampling interval is l, sampling time is t; Step 11: Use conventional PML absorbing boundary processing; Step 12: Iterate velocity and stress, and output wavefield snapshots and seismic records.

2. A two-dimensional numerical simulation method for fault detection in geological prediction according to claim 1, characterized in that: In step one: (1) According to the tunnel conditions, measure the elevation parameters of the entire tunnel within a range of 90m from the tunnel face to the tunnel direction; (2) Convert the elevation parameters into relative elevations to construct the basic form of the entire tunnel; (3) Construct an undulating tunnel model based on the basic form of the tunnel.

3. A two-dimensional numerical simulation method for geological prediction and fault detection according to claim 1, characterized in that: In the step 2: in the three-dimensional model, the lengths in the X direction and the Y direction are 300m and 100m respectively; In step three: p2 =340m / s,v s2 =0m / s,ρ2=1.29kg / m 2 .

4. A two-dimensional numerical simulation method for geological prediction and fault detection according to claim 1, characterized in that: In step six: (1) Using the wave equation to grid the grid, the size is h; (2) Set the numerical simulation parameters: frequency is f, unit is Hz, sampling interval is l, unit is ms; (3) To avoid numerical dispersion, the sampling theorem and stability conditions must be met, specifically the grid size and time step: Ⅰ. The grid spacing size meets the following requirements: Among them, v min / (2f c n) represents the minimum wavelength of the simulation, v min is the minimum velocity parameter among the numerical simulation parameters, and n is the correlation factor of the difference order taken; II. The time step size satisfies: Where m is a factor related to the order of difference taken, v max It is the minimum speed parameter among the numerical simulation parameters.

5. A two-dimensional numerical simulation method for fault detection in geological prediction according to claim 1, characterized in that: In step nine: the interval between adjacent three-component geophones is three meters, the tunnel face is provided with a first hole position for setting the three-component geophone, the hole depth of the first hole position is two meters, and the three-component geophone is arranged at the center of the side wall of the first hole position.

6. A two-dimensional numerical simulation method for fault detection in geological prediction according to claim 5, characterized in that: In the step nine: the number of the three-component detectors is no more than 31.

7. A two-dimensional numerical simulation method for geological prediction and fault detection according to claim 1, characterized in that: In step nine: the interval between adjacent seismic sources is five meters, the tunnel face is provided with a second hole position for setting the seismic source, the hole depth of the second hole position is three meters, and the seismic source is arranged at the center of the side wall of the second hole position.

8. A two-dimensional numerical simulation method for geological prediction and fault detection according to claim 1, characterized in that: In step nine: the number of the seismic sources is no more than 17.

9. A two-dimensional numerical simulation method for fault detection in geological prediction according to claim 1, characterized in that: In the step twelve: (1) This paper adopts the first-order velocity stress two-dimensional elastic wave equation of isotropic medium: In the above formula, σ xx , σ yy are the normal stresses along the X and Y directions, v x ,v y are the vibration speeds of the particle along the X and Y directions respectively; (2) Solve using high-order differences under staggered grids. The specific difference method is shown in step 8. (3) Continuously iterate velocity and stress, and output wave field snapshots and seismic record files.

Images