Joint inversion method of microseismic source and velocity based on simulated annealing method

Abstract

The present invention discloses a joint inversion method for microseismic sources and velocities based on the simulated annealing method, which can simultaneously calculate the source location and the average velocity of the underground medium, avoid the positioning error caused by inaccurate velocity parameters, and at the same time update the average velocity of the underground medium to reduce the impact on the positioning accuracy caused by the change of the underground medium. Before source location, this method does not need to give the accurate average velocity of the underground medium, only provides a rough numerical range of the average velocity. During the iterative calculation of the location inversion, it simultaneously searches for the true average velocity, jointly inverses the source location and velocity parameters, eliminates the influence of the velocity parameters on the accuracy of the positioning result, and updates the velocity model to improve the positioning accuracy.

Description

Technical Field

[0001] The present invention relates to the technical field of geophysical exploration, and particularly relates to a joint inversion method for microseismic sources and velocities based on the simulated annealing method. Background Art

[0002] The microseismic monitoring technology is used to monitor the signals of tiny fractures in underground media, and can provide technical support for the early warning and prediction of large-scale dynamic geological disasters. An important task of microseismic monitoring is to describe the three-dimensional spatial distribution of underground fractures. Therefore, an accurate microseismic positioning method is an important link in the microseismic data processing process.

[0003] The positioning method based on the seismic wave travel time establishes an objective function by using the relationship between the seismic wave travel time and the spatial coordinate positions between potential microseismic sources and several geophones, and calculates the global minimum value of the objective function through an optimization algorithm to obtain the true source position. The positioning method based on the seismic wave travel time has high calculation efficiency and accurate results, and is one of the classic methods in the field of microseismic positioning. However, such methods are easily affected by the seismic wave velocity parameters of underground media. If the velocity measurement error of underground media is large, or the velocity of underground media changes during the microseismic monitoring process, it will lead to an increase in the positioning result error and positioning failure. Summary of the Invention

[0004] In view of the above problems existing in the prior art, the present invention provides a joint inversion method for microseismic sources and velocities based on the simulated annealing method, which can simultaneously calculate the source position and the average velocity of underground media, avoid the positioning error caused by inaccurate velocity parameters, and at the same time can update the average velocity of underground media to reduce the impact on the positioning accuracy caused by the change of underground media. This method does not need to give an accurate average velocity of underground media before source positioning, only provides a rough numerical range of the average velocity. During the positioning inversion iterative calculation process, it searches for the true average velocity at the same time, jointly inverses the source position and velocity parameters, eliminates the influence of velocity parameters on the positioning result accuracy, and updates the velocity model to improve the positioning accuracy.

[0005] The present invention discloses a joint inversion method for microseismic sources and velocities based on the simulated annealing method, including:

[0006] Step 1: Set the estimated range of the equivalent average velocity in the microseismic monitoring area, where the maximum value of the equivalent average velocity is v max , the minimum value of the equivalent average velocity is v min , and uniformly discretize the estimated range of the equivalent average velocity into a total of N points to obtain N discrete estimated values of the equivalent average velocity v n , v n , where the subscript n of v is the discrete serial number, n = 1, 2,... N;

[0007] Step 2: Set the initial temperature T0 of the simulated annealing iteration, the termination temperature T M , the total number of thermal equilibrium iteration times K, and the optimal iteration vector where the subscript best represents the current optimal iteration vector (value), is the optimal iteration vector of the source location, (x best , y best , z best ) is the three-dimensional coordinate of the source location, and v best is the optimal iteration value of the equivalent average velocity;

[0008] Step 3: Establish the source iteration vector m at the instantaneous temperature T n , the discrete equivalent average velocity estimate is v in the k-th thermal equilibrium iteration and the m-th temperature drop, and the energy function and let the temperature drop times m = 0; the instantaneous temperature T m = T0;

[0009] Step 4: Let the iteration times k = 0, randomly assign values to the source location iteration vector within the monitoring area, and combine it with the N discrete equivalent average velocity estimates v n obtained in Step 1 to form N iteration vectors and substitute them into the energy function respectively to obtain N energy function values Take the smallest energy function value and the corresponding iteration vector, and record them as

[0010] Step 5: Let the iteration times k = k + 1, randomly select a new source location iteration vector within the monitoring area and update the objective function value;

[0011] Step 6: Calculate the difference

[0012] between the minimum energy functions calculated in the k-th and (k - 1)-th iterations Update the current optimal iteration vector

[0013] Step 7: According to the relationship between the iteration times k and the iteration thermal equilibrium times K, update the iteration times k and determine whether to update the instantaneous temperature T m ;

[0014] Step 8: Update the temperature drop times m and reduce the instantaneous temperature T m, determine whether the stop iteration condition is satisfied and output the current optimal iteration vector as the joint inversion result.

[0015] Furthermore, in the first step, set the equivalent average velocity estimation range within the microseismic monitoring area and perform N-point uniform discretization, and the specific form is:

[0016]

[0017] Furthermore, in the second step, the start process of the joint inversion method iteration is controlled by the initial iteration temperature T0, and the end process is controlled by the iteration termination temperature T M control. At the same temperature, the number of iterations is controlled by the total number of thermal equilibrium iterations K, and the optimal iteration vector is continuously updated during the iteration process where is the optimal iteration vector of the source location, and v best is the optimal iteration value of the equivalent average velocity.

[0018] Furthermore, in the third step, the joint inversion method obtains the minimum value of the energy function through iterative calculation to realize the joint inversion of the source location and the average propagation velocity of the source. The specific definition of the energy function is:

[0019]

[0020] where j is the detector number, J is the number of detectors, a is the amplification coefficient of the objective function, v n is the nth discrete equivalent average velocity estimation value, t j is the arrival time of the microseismic wave received by the jth detector, is the estimated value of the microseismic source occurrence time, is the distance between two vectors in three-dimensional space; at the beginning of the inversion iteration, let the number of temperature drops m = 0; the instantaneous temperature T m = T0.

[0021] Furthermore, in the fourth step, let the number of iterations k = 0, and randomly select a source location iteration vector in the monitoring area and N discrete equivalent average velocity estimation values v n to form N iteration vectors Substitute them into the energy function to obtain N energy function values Take the minimum energy function value and the corresponding iteration vector, and denote them as The specific definitions are:

[0022]

[0023]

[0024] Further, in the sixth step, calculate the difference between the minimum energy function values of the k-th and (k - 1)-th iterations. The specific method is as follows:

[0025]

[0026] Further, in the seventh step, according to the difference between the minimum energy function values of the k-th and (k - 1)-th iterations judgment factor Select a random number R within the closed interval of 0 to 1, and update the current optimal iteration vector The specific method is as follows:

[0027]

[0028] where the judgment factor When the iteration vector and the current optimal iteration vector After the update is completed, enter the eighth step.

[0029] Further, in the eighth step, judge the relationship between the iteration number k and the iteration thermal equilibrium number K. If k < K, enter the fifth step; if k = K, enter the ninth step.

[0030] Further, in the ninth step, cool down the temperature parameter that controls the iteration, let m = m + 1, T m = T0 × 0.99 m , if T m < T M , let k = 0, enter the fourth step; if T m ≥ T M , then stop the iteration and output the current optimal iteration vector

[0031] The present invention has at least the following beneficial effects:

[0032] The present invention first gives a rough velocity range of the seismic wave, uniformly discretizes the velocity range into N points, and forms N iteration vectors with the potential microseismic source spatial position vectors. When performing iterative calculations, while using the simulated annealing method to randomly search for the best position of the source in space, it also traverses and finds the velocity value that can make the objective function obtain the current minimum value. This method uses the initial temperature, termination temperature, and thermal equilibrium iteration number of the iteration to control the iterative process of the joint inversion algorithm. When the iterative calculation stops, the accurate spatial position of the source and the average propagation velocity of the seismic wave in the underground medium can be obtained. This method can effectively avoid the influence of inaccurate velocity parameters on the accuracy of the microseismic positioning algorithm based on seismic wave travel time. While obtaining an accurate positioning result, the velocity result obtained by the joint inversion can update the velocity model in real time.

[0033] Other beneficial effects of the present invention will be described in detail in the specific implementation part. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0035] Figure 1 is a flowchart of the joint inversion method of microseismic source and velocity based on the simulated annealing method disclosed in the preferred embodiment of the present invention. SPECIFIC IMPLEMENTATION MANNER

[0036] To make the objectives, technical solutions, and advantages of the present invention clearer, the following will describe the technical solutions of the present invention in detail. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other implementation manners obtained by those of ordinary skill in the art without creative efforts fall within the scope of protection of the present invention.

[0037] Figure 1 As shown, the present invention discloses a joint inversion method of microseismic source and velocity based on the simulated annealing method, including:

[0038] Step 1: Set the equivalent average velocity estimation range within the microseismic monitoring area, where the maximum value of the equivalent average velocity is v max , and the minimum value of the equivalent average velocity is v min . And uniformly discretize the equivalent average velocity estimation range into N points to obtain N discrete equivalent average velocity estimation values v n , v n . The subscript n of v is the discrete serial number, n = 1, 2,... N;

[0039] Step 2: Set the initial temperature T0 of the simulated annealing iteration, the termination temperature T M , the total number of thermal equilibrium iteration times K, and the optimal iteration vector . Where the subscript best represents the current optimal iteration vector (value), is the optimal iteration vector of the source position, (x best , y best , z best ) is the three-dimensional coordinate of the source position, and v best is the optimal iteration value of the equivalent average velocity;

[0040] Step 3: Establish at an instantaneous temperature of T m , with the discrete equivalent average velocity estimate being v n , the k-th heat balance iteration, the source iteration vector at the m-th temperature drop and the iterative energy function and let the number of temperature drops m = 0; the instantaneous temperature T m = T0;

[0041] Step 4: Let the iteration number k = 0, randomly assign values to the source position iteration vector within the monitoring area, and combine it with the N discrete equivalent average velocity estimates v n obtained in Step 1 to form N iteration vectors and substitute them into the energy function respectively to obtain N energy function values Take the smallest energy function value and the corresponding iteration vector, and denote them as

[0042] Step 5: Let the iteration number k = k + 1, randomly select a new source position iteration vector within the monitoring area and update the objective function value;

[0043] Step 6: Calculate the difference between the minimum energy functions calculated in the k-th and (k - 1)-th iterations

[0044] Step 7: Update the current optimal iteration vector according to the difference between the minimum energy functions

[0045] Step 8: According to the relationship between the iteration number k and the iteration heat balance number K, update the iteration number k and judge whether to update the instantaneous temperature T m ;

[0046] Step 9: Update the number of temperature drops m, reduce the instantaneous temperature T m

[0047] and judge whether the stop iteration condition is satisfied and output the current optimal iteration vector as the joint inversion result.Preferably, in Step 1, set the equivalent average velocity estimation range within the microseismic monitoring area and perform N-point uniform discretization, in the specific form of:

[0048]

[0049] Preferably, in Step 2, the start process of the joint inversion method iteration is controlled by the initial iteration temperature T0, and the end process is controlled by the iteration termination temperature T MControl, at the same temperature, the number of iterations is controlled by the total number of thermal equilibrium iterations K, and the optimal iteration vector is continuously updated during the iteration process. Among them, is the optimal iteration vector of the source position, and v best is the optimal iteration value of the equivalent average velocity.

[0050] Preferably, in the third step, the joint inversion method obtains the minimum value of the energy function through iterative calculation to realize the joint inversion of the source position and the average propagation velocity of the source. The specific definition of the energy function is:

[0051]

[0052] Among them, j is the detector number, J is the number of detectors, a is the amplification coefficient of the objective function, and v n is the nth discrete equivalent average velocity estimate value, t j is the arrival time of the microseismic wave received by the jth detector, is the estimated value of the microseismic source occurrence time, is the distance between two vectors in three-dimensional space; at the beginning of the inversion iteration, let the number of temperature drops m = 0; the instantaneous temperature T m = T0.

[0053] Preferably, in the fourth step, let the number of iterations k = 0, and randomly select an iteration vector of the source position in the monitoring area and N discrete equivalent average velocity estimate values v n to form N iteration vectors Substitute them into the energy function to obtain N energy function values Take the minimum energy function value and the corresponding iteration vector, and record them as Specifically defined as:

[0054]

[0055]

[0056] Preferably, in the sixth step, calculate the difference between the minimum energy function values of the kth and k-1th iterations, and the specific method is:

[0057] Preferably, in the seventh step, according to the difference between the minimum energy function values of the kth and k-1th judgment factor select a random number R in the closed interval range of 0 to 1, and update the current optimal iteration vector The specific method is:

[0058]

[0059] Among them, the judgment factor When the iteration vector and the current optimal iteration vector After the update is completed, enter Step Eight described above.

[0060] Preferably, in Step Eight described above, judge the relationship between the iteration times k and the iteration heat balance times K. If k < K, enter Step Five described above; if k = K, enter Step Nine described above.

[0061] Preferably, in Step Nine described above, perform a temperature reduction process on the temperature parameter that controls the iteration, let m = m + 1, T m = T0 × 0.99 m , if T m < T M , let k = 0, enter Step Four described above; if T m ≥ T M , then stop the iteration and output the current optimal iteration vector

[0062] The embodiments of the present invention disclose the verification of a microseismic source location method based on multi-level grid division and similarity coefficient.

[0063] I. Forward model

[0064] The forward model is a monitoring area with a range of 1000m × 1000m × 500m, and the distribution coordinates of the geophones are shown in Table 1. The average P-wave velocity of the microseismic wave is 5m / ms. Six microseismic events are set in the monitoring area, and the distribution coordinates are shown in Table 2. The seismogenic time of the microseismic is 5ms, and the arrival times of the microseismic events recorded by each geophone are shown in Table 3.

[0065] Table 1 Spatial distribution coordinates of geophones

[0066]

[0067] Table 2 Spatial distribution coordinates of microseismic events

[0068]

[0069] Table 3 Arrival times of the first arrival events picked up by each geophone

[0070]

[0071] II. Inversion results

[0072] Based on the data in Table 1 and Table 3 of the forward model, jointly invert the six source parameters and the equivalent average velocity, and the inversion parameters are shown in Table 4.

[0073] Table 4 Inversion parameters

[0074]

[0075] Through joint inversion calculation, the joint inversion results of the 6 source positions and velocities are shown in Table 5.

[0076] Table 5 Joint inversion results and errors

[0077]

[0078] From the positioning and velocity joint inversion results in Table 5, it can be seen that the joint inversion results all converge to the vicinity of the true source position and the equivalent average velocity value, indicating that the joint inversion method of microseismic source and velocity based on the simulated annealing method is accurate and effective.

[0079] As described above, it is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered within the protection scope of the present invention.

Claims

1. A joint inversion method for microseismic sources and velocities based on the simulated annealing method, characterized in that, Including: Step 1: Set the equivalent average velocity estimation range within the microseismic monitoring area, where the maximum value of the equivalent average velocity is v max , the minimum value of the equivalent average velocity is v min , and uniformly discretize the equivalent average velocity estimation range into a total of N points to obtain N discrete equivalent average velocity estimation values v n , v n , where the subscript n of v is the discrete serial number, n = 1, 2,... N; Step 2: Set the initial temperature T0 of the simulated annealing iteration, the termination temperature T M , the total number of thermal equilibrium iteration times K, and the optimal iteration vector where the subscript best represents the current optimal iteration vector, is the optimal iteration vector of the source location, (x best , y best , z best ) are the three-dimensional coordinates of the source location, and v best is the optimal iteration value of the equivalent average velocity; Step 3: Establish at the instantaneous temperature of T m , the discrete equivalent average velocity estimate is v n , the k-th thermal equilibrium iteration, the source iteration vector at the m-th temperature drop and the iteration energy function and let the number of temperature drops m = 0; the instantaneous temperature T m = T0; Step 4: Let the iteration number k = 0, and randomly assign values to the iteration vector of the source position within the monitoring area, and form N iteration vectors with the N discrete equivalent average velocity estimates v n obtained in Step 1 Substitute them into the energy function respectively, and obtain N energy function values Take the minimum energy function value and the corresponding iteration vector among them, and denote them as Step 5: Let the iteration number \(k = k + 1\), and randomly select a new seismic source position iteration vector within the monitoring area range and update the objective function value; Step Six: Calculate the difference in the minimum energy function between the k-th and (k - 1)-th iterative calculations Step Seven: Update the current optimal iteration vector according to the difference in the minimum energy function Update the current optimal iteration vector Step Eight: Update the iteration number k according to the relationship between the iteration number k and the iteration heat balance number K, and determine whether to update the instantaneous temperature T m ; Step Nine: Update the number of times of temperature drop \(m\) and reduce the instantaneous temperature \(T\). m , determine whether the stop iteration condition is satisfied and output the current optimal iteration vector as the joint inversion result.

2. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the first step described above, set the equivalent average velocity estimation range within the microseismic monitoring area and perform N-point uniform discretization, and the specific form is:

3. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the second step described above, the start process of the joint inversion method iteration is controlled by the initial iteration temperature T0, and the end process is controlled by the iteration termination temperature T M The number of iterations at the same temperature is controlled by the total number of thermal equilibrium iterations K, and the optimal iteration vector is continuously updated during the iteration process Among them, is the optimal iteration vector of the source position, and v best is the optimal iteration value of the equivalent average velocity.

4. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the third step described above, the joint inversion method obtains the minimum value of the energy function through iterative calculation to achieve the joint inversion of the source position and the average propagation velocity of the source. The specific definition of the energy function is: where j is the geophone serial number, J is the number of geophones, a is the amplification coefficient of the objective function, v n is the estimated value of the nth discrete equivalent average velocity, t j is the arrival time of the microseismic wave received by the jth geophone, is the estimated value of the occurrence time of the microseismic source, is the distance between two vectors in three-dimensional space; at the beginning of the inversion iteration, let the number of temperature drops m = 0; the instantaneous temperature T m = T0.

5. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the fourth step described above, let the iteration number k = 0, and randomly select a source position iteration vector in the monitoring area and N discrete equivalent average velocity estimates v n to form N iteration vectors Substitute them into the energy function to obtain N energy function values Take the minimum energy function value and the corresponding iteration vector, and denote them respectively as Specifically defined as:

6. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the fifth step described above, let k = k + 1, and randomly update the iterative vector of the seismic source position Recalculate the new iterative vector 7. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the sixth step described above, calculate the difference between the minimum energy functions of the k-th and (k-1)-th iterations, and the specific method is:

8. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the seventh step described above, according to the smallest difference in the energy function between the k-th and (k - 1)-th times judgment factor Select a random number R within the closed interval of 0 to 1, and update the current optimal iteration vector The specific method is as follows: Among them, the judgment factor When the iterative vector and the current optimal iterative vector After the update is completed, enter Step Eight described above.

9. The joint inversion method for microseismic sources and velocities based on the simulated annealing method according to claim 1, characterized in that, In the eighth step described above, judge the relationship between the iteration number k and the iteration thermal equilibrium number K. If k < K, enter the fifth step described above; if k = K, enter the ninth step described above.

10. The joint inversion method of microseismic source and velocity based on the simulated annealing method according to claim 1, characterized in that In the ninth step described above, the temperature parameter for controlling iteration is cooled down, and let m = m + 1, T m = T0 × 0.99 m , if T m < T M , let k = 0, and enter the fourth step described above; if T m ≥ T M , then stop the iteration and output the current optimal iteration vector

Images