A CT inversion method for the integrated fusion of dual-source vibration waves

Abstract

The present invention discloses a CT inversion method with integrated fusion of dual-source shock waves. The method adopts the mutual fusion of active source and passive source. A small number of artificial seismic sources are deployed to obtain the active source data. Then, the passive source data collected by the microseismic monitoring system is obtained. The two are fused for joint inversion, thereby obtaining an inversion result. Due to the combination of the passive source data, only a small number of artificial seismic sources need to be deployed during inversion, so the scope of downtime required for inversion is relatively small. At the same time, due to the addition of the active seismic source, the required passive source data is greatly reduced. Under the premise of ensuring the inversion accuracy, the inversion cycle of the passive source CT can be shortened at the same time, thereby eliminating the need to wait for a long time for inversion each time. Therefore, under the premise of ensuring the inversion accuracy, the present invention can reduce the number of artificial seismic sources deployed and the scope of downtime required for inversion, and can greatly shorten the inversion cycle required for each inversion of the passive source, thereby facilitating long-term monitoring of dangerous areas.

Description

Technical Field

[0001] The invention relates to a CT inversion method for integrated fusion of dual-source shock waves, and belongs to the technical field of safe mining of coal mines. Background Art

[0002] Currently, shockwave CT inversion technology, which uses microseismic location information and the arrival times of channel markers to determine the stress distribution in the monitored area, is widely used for coal mine rock burst monitoring and early warning. This technology requires only a large number of shockwave rays to pass through the monitored area to accurately determine its distribution, thus potentially identifying hidden rock burst risk areas. However, it requires sufficient mine earthquakes to generate sufficient shockwave ray coverage within the monitored area. Because mine earthquakes occur randomly during production and their source locations are unknown, current passive-source shockwave CT inversion technology must calculate the arrival times of shockwaves at each station, resulting in a certain deviation between the calculated source location and the actual location. Furthermore, based on previous engineering experience, in order to obtain a large number of shockwaves generated by mine earthquakes (typically, each inversion requires approximately 300 microseismic events, generating approximately 1000 rays and covering a range greater than 500 meters), the interval between mine shockwave CT inversions is typically two weeks. This long warning interval significantly impacts the timeliness of early warnings. To address the shortcomings of the aforementioned passive seismic sources, artificial blasting is used to actively stimulate shock wave signals. This method can directly obtain the source location and the warning time is controllable. By setting up a large number of artificial seismic sources, radiation coverage of the monitoring area can be quickly achieved. However, the disadvantage of actively stimulating shock waves is that in order to fully cover the monitoring area with radiation generated by the artificial seismic sources, a large number of artificial seismic sources need to be carried down the mine and deployed according to the scope of the monitoring area. This results in a large amount of labor for workers and time-consuming deployment. Once the monitoring area changes, the position of each artificial seismic source needs to be adjusted. In addition, since the coal mine in the entire monitoring area needs to be shut down to reduce interference during inversion detection, it is not conducive to long-term research in dangerous areas.

[0003] Therefore, how to provide a new method that can greatly reduce the number of artificial seismic sources deployed while ensuring the timeliness of early warning, thereby reducing the workload of workers, reducing the time required for deployment, and facilitating long-term monitoring of dangerous areas, is one of the research directions of this industry. Summary of the Invention

[0004] In response to the problems existing in the above-mentioned prior art, the present invention provides a CT inversion method with integrated fusion of dual-source vibration waves. Under the premise of ensuring the timeliness of early warning, it can not only reduce the number of artificial seismic sources deployed and the scope of work suspension required for inversion, but also greatly shorten the inversion cycle required for each inversion of the passive source, thereby facilitating long-term monitoring of dangerous areas.

[0005] To achieve the above objectives, the present invention adopts a technical solution: a CT inversion method with integrated fusion of dual-source shock waves, the specific steps of which are as follows:

[0006] (1) Using the mine earthquakes induced during mining production recorded by the installed microseismic monitoring system underground as passive sources, the vibration wave signals of artificial blasting are collected as active sources; according to the range of wave velocities to be solved, a certain number of active and passive sources are selected to form the initial data for dual-source integrated inversion, and the model size and number of grid divisions are determined according to the distribution characteristics of the dual sources;

[0007] (2) performing active tomography on the active source data obtained in step (1), and selecting a ray tracing algorithm to calculate and obtain an initialization velocity model;

[0008] (3) Solve the passive source parameters θ(x0, y0, z0, t0) obtained in step (1) using the initial wave velocity model obtained in step (2) according to the microseismic positioning algorithm, and solve the covariance matrix C of θ θ , by the covariance matrix C θ The eigenvalue of Calculate the volume v of the error ellipsoid, where

[0009] (4) using the ray tracing algorithm determined in step (2), determine the distance matrix G of the active source and the passive source's respective rays passing through the model unit grid;

[0010] (5) According to the source type and the error ellipsoid volume v, each ray is assigned a corresponding weight factor w i , and calculate the dual-source cumulative correction factor Δp for each pixel j on the ray j ;

[0011] (6) According to the inversion algorithm (G T G+ε 2 I) S=G T (dT 0 ) realizes the adjustment of the velocity model under the dual-source integration fusion, and repeats steps (3) to (5) after obtaining the new velocity model until the inversion iteration termination condition is met, and finally outputs the inversion result, where T 0 is the earthquake occurrence time vector with length m. For the active source T 0 = 0, for passive source T 0 The value of is determined by the microseismic location algorithm.

[0012] Furthermore, in step (1), the method for determining the size of the model and the number of grid divisions according to the distribution characteristics of the dual sources includes an equal-interval division method and an adaptive division method.

[0013] Furthermore, the ray tracing algorithm selected in step (2) is specifically selected from a straight ray algorithm, a ray bending algorithm, a shortest path algorithm and a hybrid algorithm.

[0014] Furthermore, the specific process of obtaining the initialization velocity model in step (2) is: using active tomography, the shock wave theory to time t i Discretization, where V(x,y,z) is the velocity distribution, t i is the theoretical calculation of the arrival time of the shock wave, Γ i is the propagation ray, and S(x,y,z) is equal to 1 / V(x,y,z), which is the slowness value at the point (x,y,z); the linear equation GS=d is obtained, where d represents the calculated arrival time vector with a length of m, and S represents the model parameter vector with a length of n, i.e., the slowness; and the equation form of the least squares solution is established (G T G+ε 2 I) S=G T (d) A stable initialization velocity model is obtained under the selected ray tracing algorithm, where ε is a small positive scalar.

[0015] Furthermore, the specific process of step (5) is as follows: wherein the method for determining the weight of each ray is: if the ray belongs to an active source, the weight factor w i =1; if the ray belongs to a passive source, the weight factor w is determined by the ellipsoid volume i Value, that is The larger the volume v of the ellipsoid, the larger the w i The smaller, the worse. i The smaller it is, the weight factor w is given to each ray. i ;

[0016] Calculate the dual-source cumulative correction factor Δp for each pixel j on the ray j Specifically:

[0017] First calculate the slowness correction factor Δp applied to grid pixel j on the path of ray i ij , by adding the slowness correction factor, we can get the cumulative correction factor Δp on each pixel. j , the specific formula is: where Δt i is the residual between the theoretical arrival time and the marked arrival time of the shock wave on each ray, G ij is the path segment length, N p is the number of rays passing through pixel j, w i Weight factors determined for different ray types.

[0018] Furthermore, the objective function of the microseismic location algorithm in step (6) is:

[0019]

[0020] Where: w i is the weight function of the observation values of each station in the microseismic monitoring system, n is the number of stations in the microseismic monitoring system marked by P waves, x0, y0, z0 are the coordinates of the passive source, x i ,y i , z i is the station coordinate of the microseismic monitoring system; the positional relationship between the coordinates of the passive source and the coordinates of each station is calculated according to the above formula to determine the passive source T 0 value.

[0021] Since active source CT inversion has a clear excitation position and time, the input parameter error of its wave velocity inversion model is small and the inversion result is accurate. However, its signal quality is contradictory to production activities. In order to ensure clear and stable signals, production activities within the inversion range need to be stopped, and it requires the deployment of more artificial seismic sources. Passive source CT inversion has the advantages of a wide detection range and no need for shutdown, but each inversion takes a long time (generally at least 10 days as a cycle), and its detection accuracy is lower than that of active source CT inversion. Compared with the prior art, the present invention adopts a fusion of active and passive sources, by deploying a small number of artificial seismic sources to excite active seismic waves and obtain the monitoring data, and then obtains the passive source data monitored and collected by the microseismic monitoring system, and fuses the two for joint inversion to obtain the inversion result. Since the present invention combines the data of passive sources, only a small number of artificial seismic sources need to be deployed to excite active seismic waves during inversion, so the scope of shutdown required is small, the impact on the entire coal mine is small, and the deployment work is convenient. At the same time, due to the addition of active seismic sources, the required passive source data is greatly reduced. Under the premise of ensuring the inversion accuracy, the inversion cycle of passive source CT can be shortened at the same time. The inversion cycle can be shortened to 3 days or even shorter, so there is no need to wait for a long time for inversion each time. In summary, the present invention can reduce the number of artificial seismic sources deployed and the scope of shutdown required for inversion while ensuring inversion accuracy, and can greatly shorten the inversion cycle required for each passive source inversion, thereby facilitating long-term monitoring of dangerous areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is a schematic diagram of the overall process of the present invention. DETAILED DESCRIPTION

[0023] The present invention will be further described below.

[0024] like Figure 1 As shown, the specific steps of the present invention are:

[0025] (1) Using the mine earthquakes induced during mining production recorded by the installed microseismic monitoring system underground as the passive source, the vibration wave signals of artificial blasting are collected as the active source; according to the range of wave velocities to be solved, a certain number of active and passive sources are selected to form the initial data for dual-source integrated inversion, and the model size and number of grid divisions are determined by using the equal-interval division method or the adaptive division method according to the distribution characteristics of the dual sources;

[0026] (2) Perform active tomography on the active source data obtained in step (1), and select one of the algorithms from the straight ray algorithm, ray bending algorithm, shortest path algorithm and hybrid algorithm as the ray tracing algorithm for this inversion according to the actual situation of the active source data, and obtain the initialization velocity model after calculation; the specific process is: using active tomography, the shock wave theory to time t i Discretization, where V(x,y,z) is the velocity distribution, t i is the theoretical calculation of the arrival time of the shock wave, Γ i is the propagation ray, and S(x,y,z) is equal to 1 / V(x,y,z), which is the slowness value at the point (x,y,z); the linear equation GS=d is obtained, where d represents the calculated arrival time vector with a length of m, and S represents the model parameter vector with a length of n, i.e., the slowness; and the equation form of the least squares solution is established (G T G+ε 2 I) S=G T (d) A stable initialization velocity model is obtained under the selected ray tracing algorithm, where ε is a small positive scalar.

[0027] (3) Solve the passive source parameters θ(x0, y0, z0, t0) obtained in step (1) using the initial wave velocity model obtained in step (2) according to the microseismic positioning algorithm, and solve the covariance matrix C of θ θ , by the covariance matrix C θ The eigenvalue of Calculate the volume v of the error ellipsoid, where

[0028] (4) using the ray tracing algorithm determined in step (2), determine the distance matrix G of the active source and the passive source's respective rays passing through the model unit grid;

[0029] (5) According to the source type and the error ellipsoid volume v, each ray is assigned a corresponding weight factor w i , and calculate the dual-source cumulative correction factor Δp for each pixel j on the ray j ; The method for determining the weight of each ray is: if the ray belongs to an active source, the weight factor w i =1; if the ray belongs to a passive source, the weight factor w is determined by the ellipsoid volumei Value, that is The larger the volume v of the ellipsoid, the larger the w i The smaller, the worse. i The smaller it is, the weight factor w is given to each ray. i ;

[0030] Calculate the dual-source cumulative correction factor Δp for each pixel j on the ray j Specifically:

[0031] First calculate the slowness correction factor Δp applied to grid pixel j on the path of ray i ij , by adding the slowness correction factor, we can get the cumulative correction factor Δp on each pixel. j , the specific formula is: where Δt i is the residual between the theoretical arrival time and the marked arrival time of the shock wave on each ray, G ij is the path segment length, N p is the number of rays passing through pixel j, w i Weight factors determined for different ray types.

[0032] (6) According to the inversion algorithm (G T G+ε 2 I) S=G T (dT 0 ) realizes the adjustment of the velocity model under the dual-source integration fusion, and repeats steps (3) to (5) after obtaining the new velocity model until the inversion iteration termination condition is met, and finally outputs the inversion result, where T 0 is the earthquake occurrence time vector with length m. For the active source T 0 = 0, for passive source T 0 The value of is determined by the microseismic positioning algorithm, and the objective function of the microseismic positioning algorithm is:

[0033]

[0034] Where: w i is the weight function of the observation values of each station in the microseismic monitoring system, n is the number of stations in the microseismic monitoring system marked by P waves, x0, y0, z0 are the coordinates of the passive source, x i ,y i , z i is the station coordinate of the microseismic monitoring system; the positional relationship between the coordinates of the passive source and the coordinates of each station is calculated according to the above formula to determine the passive source T 0 value.

[0035] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A dual-source seismic wave integrated fusion CT inversion method, characterized in that: The specific steps are: (1) Using the mine earthquakes induced during mining production recorded by the installed microseismic monitoring system underground as passive sources, the vibration wave signals of artificial blasting are collected as active sources; according to the range of wave velocities to be solved, a certain number of active and passive sources are selected to form the initial data for dual-source integrated inversion, and the model size and number of grid divisions are determined according to the distribution characteristics of the dual sources; (2) performing active tomography on the active source data obtained in step (1), and selecting a ray tracing algorithm to calculate and obtain an initialization velocity model; (3) Solve the passive source parameters θ(x0, y0, z0, t0) obtained in step (1) using the initial wave velocity model obtained in step (2) according to the microseismic positioning algorithm, and solve the covariance matrix C of θ θ , by the covariance matrix C θ The eigenvalue of Calculate the volume v of the error ellipsoid, where (4) using the ray tracing algorithm determined in step (2), determine the distance matrix G of the active source and the passive source's respective rays passing through the model unit grid; (5) According to the source type and the error ellipsoid volume v, each ray is assigned a corresponding weight factor w i , and calculate the dual-source cumulative correction factor Δp for each pixel j on the ray j ; (6) According to the inversion algorithm (G T G+ε 2 I) S=G T (dT 0 ) realizes the adjustment of the velocity model under the dual-source integrated fusion. After obtaining the new velocity model, steps (3) to (5) are repeated until the inversion iteration termination condition is met, and finally the inversion result is output, where ε is a small positive scalar; d represents the calculated arrival time vector with a length of m; S represents the model parameter vector with a length of n, that is, slowness; T 0 is the earthquake occurrence time vector with length m. For the active source T 0 = 0, for passive source T 0 The value of is determined by the microseismic location algorithm.

2. The CT inversion method of dual-source seismic wave integration fusion according to claim 1 is characterized in that: In the step (1), the method for determining the size of the model and the number of grid divisions according to the distribution characteristics of the dual sources includes an equal-interval division method and an adaptive division method.

3. The CT inversion method with integrated dual-source seismic wave fusion according to claim 1 is characterized in that: The ray tracing algorithm selected in step (2) specifically includes: selecting one of the straight ray algorithm, ray bending algorithm, shortest path algorithm and hybrid algorithm.

4. The CT inversion method of dual-source seismic wave integration fusion according to claim 1 is characterized in that: The specific process of obtaining the initialization velocity model in step (2) is: using active tomography to convert the shock wave theory to time t i Discretization, where V(x,y,z) is the velocity distribution, t i is the theoretical calculation of the arrival time of the shock wave, Γ i is the propagation ray, and S(x,y,z) is equal to 1 / V(x,y,z), which is the slowness value at the point (x,y,z); the linear equation GS=d is obtained, where d represents the calculated arrival time vector with a length of m, and S represents the model parameter vector with a length of n, i.e., the slowness; and the equation form of the least squares solution is established (G T G+ε 2 I) S=G T (d) A stable initialization velocity model is obtained under the selected ray tracing algorithm, where ε is a small positive scalar.

5. The CT inversion method of dual-source seismic wave integration fusion according to claim 1 is characterized in that: The specific process of step (5) is as follows: the method for determining the weight of each ray is: if the ray belongs to an active source, the weight factor w i =1; if the ray belongs to a passive source, the weight factor w is determined by the ellipsoid volume i Value, that is The larger the volume v of the ellipsoid, the larger the w i The smaller, the worse. i The smaller it is, the weight factor w is given to each ray. i ; Calculate the dual-source cumulative correction factor Δp for each pixel j on the ray j Specifically: First calculate the slowness correction factor Δp applied to grid pixel j on the path of ray i ij , by adding the slowness correction factor, we can get the cumulative correction factor Δp on each pixel. j , the specific formula is: where Δt i is the residual between the theoretical arrival time and the marked arrival time of the shock wave on each ray, G ij is the path segment length, N p is the number of rays passing through pixel j, w i Weight factors determined for different ray types.

6. The CT inversion method of dual-source seismic wave integration fusion according to claim 1 is characterized in that: The objective function of the microseismic location algorithm in step (6) is: Where: w i is the weight function of the observation values of each station in the microseismic monitoring system, n is the number of stations in the microseismic monitoring system marked by P waves, x0, y0, z0 are the coordinates of the passive source, x i ,y i , z i is the station coordinate of the microseismic monitoring system; the positional relationship between the coordinates of the passive source and the coordinates of each station is calculated according to the above formula to determine the passive source T 0 value.

Images