CN106772577A - Source inversion method based on microseism data and SPSA optimized algorithms - Google Patents

Abstract

The invention relates to a source inversion method based on microseismic data and SPSA optimization algorithm. The method includes the following steps: establishing a horizontal layered velocity model according to well logging data and underground rock property data analysis; performing forward modeling on microseismic events Simulation, adjust the ray angle by dichotomy to optimize the ray tracing path, calculate the first arrival travel time, and then calculate the travel time difference between two adjacent geophones; add the travel time difference to random disturbance as the observation value; apply the SPSA algorithm to iteratively update the obtained The microseismic source position is calculated according to the aforementioned method to obtain the difference between the first arrival travel time, and the source position is continuously optimized by fitting the calculated data and the observed value, and the difference between the calculated data and the observed value is checked. If the accuracy is satisfied or the maximum number of iterations is reached, the update is stopped, and the output As a result, positioning is done, otherwise proceed to the next iteration. The present invention improves the accuracy of seismic source location by quickly and efficiently solving microseismic event sources.

Description

Source inversion method based on microseismic data and SPSA optimization algorithm

technical field

The invention belongs to the field of oil and gas seismic exploration, and in particular relates to a seismic source inversion method based on microseismic data and an SPSA optimization algorithm.

Background technique

In the production process of low-permeability oil and gas fields, fracturing technology is an important measure to increase production and injection, and the fractures generated by fracturing (the fracture orientation is closely related to in-situ stress) and fracturing scale are important factors for well pattern deployment. Therefore, in-depth study of fracture azimuth and shape, and timely adjustment of well patterns have always been urgent issues to be solved in oilfields.

Artificial microseismic real-time monitoring and evaluation technology for fracturing wells is a dynamic monitoring technology for oilfield production based on microseismic monitoring technology. Microseismic monitoring is currently the most accurate, timely and informative monitoring method in reservoir fracturing. Real-time microseismic imaging can guide fracturing projects in a timely manner and adjust fracturing parameters in a timely manner; track and locate the range of fracturing, the direction and size of fracture development, objectively evaluate the effect of fracturing projects, and provide effective information for the next step of production and development. guidance, reduce development costs, and provide a better basis for subsequent oilfield management.

The core issue of data processing in microseismic monitoring is source location, and the most commonly used inversion method based on travel time, such as longitudinal and transverse wave time difference method or homotype wave time difference method, is a conventional method in microseismic source inversion. However, this method needs to solve linear or nonlinear equations, and there may be overdetermined, underdetermined, or ill-conditioned problems in the process of solving the equations. In addition, this method relies heavily on the event time detected by the geophone, and the positioning accuracy cannot meet the requirements. .

There are also various inversion methods commonly used at present, including some gradient algorithms, such as Newton’s method and conjugate gradient method. The selection is highly dependent; there are also some stochastic algorithms, such as particle swarm algorithm, genetic algorithm, simulated annealing method, etc., which can solve the problem of local optimal solution well, but the solution speed of this type of algorithm is relatively slow.

Contents of the invention

In order to overcome the defects in the prior art, the present invention provides a seismic source inversion method based on microseismic data and SPSA optimization algorithm, which can quickly and efficiently realize seismic source positioning.

To achieve the above object, the present invention adopts the following scheme:

The microseismic source location method based on the SPSA algorithm includes the following steps:

Step 1: Establish a horizontal layered velocity model based on well logging data and analysis of underground rock properties;

Step 2: Perform forward modeling on microseismic events, calculate the first arrival travel time, and then calculate the travel time difference between two adjacent geophones;

Step 3: Add the travel time difference calculated in Step 2 to a random disturbance composed of a group of random numbers subject to Gaussian distribution within (0, 1) as the observed value;

Step 4: Apply the SPSA algorithm to optimize and solve the optimal source location by fitting the calculated data and the observed values.

In terms of ray tracing positioning, aiming at the shortcomings of the local ray tracing method that relies heavily on the ray density factor to cause poor accuracy, the present invention proposes to use the dichotomy method to improve the ray tracing shooting method based on Snell's law. First, the initial ray angle is given, and then According to the continuous adjustment of the dichotomy method, it has been tested that this method not only improves the convergence speed, but also greatly improves the accuracy, which improves the speed and accuracy of positioning.

In the selection of the objective function, the present invention uses the difference in travel time of the same seismic source monitored by adjacent geophones as a fitting parameter, which not only avoids the cumbersome calculation of the occurrence time of the seismic source, but also reduces the seismic wave received during the stratum propagation process. The degree of interference improves the noise immunity.

In addition, the present invention also proposes to use the SPSA algorithm to invert the source position. This method combines the gradient algorithm and the random algorithm together, uses the random algorithm to give the search direction and the random gradient, and then uses the gradient optimization to solve the problem. It can consider the global optimum, and can realize the fast solution of the objective function, has a fast convergence speed, can realize fast and accurate seismic source positioning, and has the advantages of fast calculation speed and high efficiency.

Description of drawings

Figure 1 is a schematic diagram of the source inversion method based on microseismic data and SPSA optimization algorithm;

Fig. 2 is a schematic diagram of ray tracing positioning based on Snell's law;

Figure 3 is a flow chart of SPSA algorithm optimization inversion;

Figure 4 is a ray tracing path diagram obtained by forward modeling;

Figure 5 is the diagram of the travel time of the first arrival obtained from the forward modeling;

Figure 6 is a diagram of the first arrival travel time difference obtained by the forward modeling;

Figure 7 is an illustration of the observation data obtained after adding random disturbance;

Figure 8 is a diagram of the real source location and the source location obtained by inversion.

detailed description

As shown in Figure 1, the source inversion method based on microseismic data and SPSA optimization algorithm includes the following steps:

Step 1: Establish a horizontal layered velocity model based on the analysis of logging data and underground rock property data;

The specific method is as follows:

① According to the logging data and analysis of underground rock properties, strata are classified and divided into layers. The stratum is divided into M layers, which are numbered from deep to shallow. The medium interfaces are all horizontal planes, and the depth of each layer from the ground is Layer_h= [Layer_h ₁ Layer_h ₂ ... Layer_h _M ];

② Assuming that the oil layer and overlying rock are homogeneous and the porous medium is filled with fluid, the formation velocity field can be established according to the Raymer model:

v＝(1-φ) ² ·v _solid +φ·v _fluid

in,

K _flud =K _w ×S _w +K _o ×S _o +K _g ×S _g ;

ρ _fliud = ρ _w × S _w + ρ _o × S _o + ρ _g × S _g ;

In the formula, v is velocity (m/s), φ is porosity, K is elastic modulus (Pa), μ is shear modulus (Pa), ρ is density (kg/cm ³ ), f _i is solid The volume fraction of the i-th composition in , N _solid is the number of components contained in the solid, S _m is the saturation of phase m; the subscripts solid and fluid represent solid and fluid respectively, and w, o and g represent water respectively , oil and gas three-phase.

It can be seen from the above formula that the Raymer model comprehensively considers the basic characteristics of underground solids and fluids, and the above parameter values can be obtained through the analysis of well logging data and underground rock properties, so that the velocity of seismic wave propagation in each layer can be calculated. Since the P wave (P wave) has the advantages of fast propagation speed and easy identification, the P wave at the layer interface is usually selected for source location, so that the velocity field from shallow to deep is Layer_v=[Layer_v ₁ Layer_v ₂ ...Layer_v _M ] .

Step 2: Perform forward modeling on microseismic events, calculate the first arrival travel time, and then calculate the travel time difference between two adjacent geophones;

The specific method is as follows:

① Ray tracing based on Snell's law is defined as follows:

The strata in the test area are divided into M layers, the propagation velocity of P wave (longitudinal wave) in each layer from shallow to deep is Layer_v=[Layer_v ₁ Layer_v ₂ ... Layer_v _M ], the number of receivers is N, and the seismic wave starts from point S The calculation formulas for the propagation angle, propagation path and travel time of the ray reaching the i-th stage geophone in the formation are as follows:

Among them, s _i is the path of the seismic wave starting from point S to the i-th level geophone and the ray in the stratum, n _i is the geophone number (the number in the corresponding stratum from top to bottom), h _ij is the longitudinal direction of the ray Projection height, _θij is the angle (the angle with the vertical direction) when the seismic wave passes through the jth layer when it propagates to the i-level geophone, and Layer_v _j is the velocity of the seismic wave propagating in the jth layer, is the travel time of the ray in the formation when the seismic wave departs from point S and arrives at the i-th level geophone.

② When the improved shooting method is used for ray tracing, the launch angle θ of the initial position of the ray is given first:

a=0,b=90

θ=(a+b)/2(θ∈[a,b])

Among them, θ is the launch angle of the initial position of the ray, and a and b are the upper and lower limits of the launch angle of the initial position.

Starting from the source point, the propagation of seismic waves between layers follows Snell’s law, and it is judged whether the end position of the ray is the position of the receiving point. If the end position of the ray is not the position of the receiving point, if it is higher than the position of the receiving point, then:

b=(a+b)/2

θ=(a+b)/2(θ∈[a,b])

Conversely, there are:

a=(a+b)/2

θ=(a+b)/2(θ∈[a,b])

Constantly adjust the ray angle according to the dichotomy, and recalculate until it meets the requirements.

③Because the duration of the fracturing process is relatively long, the occurrence time of the seismic source is not easy to determine. Therefore, the time difference between the first arrivals of the geophones can be calculated to avoid the occurrence time of the seismic source, and it can also effectively reduce the degree of interference in the formation. The calculation is as follows :

in, is the travel time difference between two adjacent detectors.

Step 3: Add the difference in travel time in step 2 to the random disturbance composed of a group of random numbers subject to Gaussian distribution within (0,1) as the observed value;

The specific method is as follows:

Generate random disturbances composed of a set of random numbers subject to Gaussian distribution within (0,1)

Δt＝[Δt ₁ Δt ₂ ... Δt _N-1 ]

Then the observed value is

Among them, Δt is the random disturbance vector, and ΔT _obs is the observed data vector.

Step 4: Apply the SPSA algorithm to optimize and solve the optimal source location by fitting the calculated data and the observed values;

The specific method is as follows:

① The initial position of the microseismic event X ₀ = (x ₀ , y ₀ , z ₀ ) is randomly given, and the travel time of the first arrival is calculated as described in step (2) Then get the difference of travel time of first arrival

Get the difference between calculated data and observed data

Given k=0 (k is the current number of searches), maxIter=100 (the maximum number of searches).

Among them, ΔT _cal0 is the difference vector of the first arrival travel time at the initial position, and ||ΔT ₀ || is the difference between the calculated data and the observed data at this time.

② Add 1 to the number of searches k, randomly given an iteration direction α (between [-1,1] and not 0) and direction step size ΔX, then

k=k+1

X ₁ =X ₀ +αΔX

Also calculate the travel time of the first arrival as described in step 2 Then get the difference of travel time of first arrival

The difference between the calculated data and the observed data can also be obtained

If ||ΔT ₀ ||＞||ΔT ₁ ||, the selected direction is otherwise

Among them, ΔT _cal1 is the difference vector of the first arrival travel time obtained by calculation and update, ||ΔT ₁ || is the difference between the calculated data and the observed data at this time, α is the iteration direction, and ΔX is the step size of the iteration direction.

③ Calculate and update the location of microseismic events

X _update = X ₀ +αgStep

Among them, X _update is the updated source position, and Step is the iteration step.

④ Calculate the travel time of the first arrival according to step 2 Then get the difference of travel time of first arrival

Get the difference between calculated data and observed data

Among them, g(X _update ) is the objective function of iterative update, ΔT _update is the updated first arrival travel time difference vector, and ΔT _obs is the observation data vector.

⑤ Determine whether g(X _update ) satisfies the precision, if so, output X _best ＝X _update ; otherwise, X ₀ ＝X _update , return to continue execution ③Iteratively update the microseismic event position until g(X _update ) no longer drops, then return to execute ②, judge whether k is greater than max Iter, if so, output X _best = X _update , otherwise reselect the iteration direction α, and add 1 to the search times k. In addition, if g(X _update ) meets the precision requirement during the iterative update process, output X _best =X _update and stop the update.

Example

Based on step 1, a horizontal layered velocity model is established based on well logging data and analysis of underground rock properties: In this calculation example, an area with a length of 400 meters, a width of 400 meters, and a height of 400 meters is selected. They are numbered sequentially from shallow to shallow, and the medium interface is a horizontal plane. The depth of each layer from the ground is Layer_h=[2100 2150 2300 2400]m, and the P wave (longitudinal wave) velocity of the layer interface is Layer_v=[1600 2000 2400] from shallow to deep. 2800] m/s.

Based on the second step, forward modeling is carried out for a given microseismic event, and the first arrival travel time is calculated, and then the travel time difference between two adjacent geophones is calculated:

Set the detector parameters as shown in Table 1 below:

Table 1 Geophone parameters

series 1 2 3 4 5 6 7 8 9 10 11 12 x 0 0 0 0 0 0 0 0 0 0 0 0 Y 0 0 0 0 0 0 600 600 600 600 600 600 Z 2017 2032 2047 2062 2077 2092 2107 2122 2137 2152 2167 2182 series 13 14 15 16 17 18 19 20 twenty one twenty two twenty three twenty four x 600 600 600 600 600 600 600 600 600 600 600 600 Y 0 0 0 0 0 0 600 600 600 600 600 600 Z 2197 2212 2227 2242 2257 2272 2287 2302 2317 2332 2347 2362

The Z coordinate is the depth relative to the ground. The depth of the first stage geophone is -2017 meters, the depth of the 24th stage geophone is -2362 meters, and the height difference of each stage geophone is 15 meters.

Set the real source event location as shown in Table 2 below:

Table 2 Parameters of real source events

According to the improved ray tracing shooting method and adjusting the ray angle based on the dichotomy method, the ray tracing path obtained by forward modeling is shown in Figure 4.

After the path is obtained, divide it by the seismic wave propagation velocity of the layer to get the distribution of the first arrival travel time from the source to the geophone position, as shown in Figure 5, and then get the difference distribution of the first arrival travel time, as shown in Figure 6.

Based on step 3, add the difference in travel time calculated in step 2 to a random disturbance composed of a group of random numbers subject to Gaussian distribution within (0,1), and use this as the observed value;

Generate a random disturbance Δt=[Δt ₁ Δt ₂ ... Δt _N-1 ] composed of a group of random numbers subject to Gaussian distribution within (0,1), then the observed value is

After adding random disturbance, the observed data distribution is shown in Figure 7.

Based on step 4, the SPSA algorithm is applied to optimize and solve the optimal hypocenter position by fitting the calculated data and the observed values:

See Figure 3 for the flow chart of the inversion of the source position using the SPSA optimization method. Follow the steps shown in Figure 3 to perform the inversion. The source parameters obtained by the inversion are shown in Table 3:

Table 3 Inversion source location

The real source location and the source location obtained by inversion are shown in Fig. 8.

From Table 3 and Figure 8, it can be seen intuitively that after a certain disturbance is added to the real value, the maximum error between the real source position and the source position obtained by inversion is 21.9219m, and the minimum is 4.0864m. The total performance time is about 140s, and the inversion result is within a reasonable range.

Claims (5)

1. a source inversion method based on microseismic data and SPSA optimization algorithm, is characterized in that, comprises the following steps: Step 1: Establish a horizontal layered velocity model based on well logging data and analysis of underground rock properties; Step 2: Perform forward modeling on microseismic events, calculate the first arrival travel time, and then calculate the travel time difference between two adjacent geophones; Step 3: Add the travel time difference calculated in Step 2 to a random disturbance composed of a group of random numbers subject to Gaussian distribution within (0,1), and use this as the observed value; Step 4: Apply the SPSA algorithm to optimize and solve the optimal source location by fitting the calculated data and the observed values. 2. the seismic source inversion method based on microseismic data and SPSA optimization algorithm as claimed in claim 1, is characterized in that, step one comprises the following steps: ① According to the logging data and the analysis of underground rock properties, strata are classified and divided into layers. The strata are divided into M layers, numbered in sequence from deep to shallow, and the medium interfaces are all horizontal planes. The depth of each layer from the ground is given; (2) Assuming that the oil layer and overlying rock are homogeneous, and the porous medium is filled with fluid, the formation compressional wave (P wave) velocity field can be established according to the Raymer model. 3. the seismic source inversion method based on microseismic data and SPSA optimization algorithm as claimed in claim 1-2, is characterized in that, step 2 comprises the following steps: (1) Ray tracing positioning method based on Snell's law: Assuming that the formation in the test area is divided into M layers, according to Snell's law, the seismic wave starting from point S and arriving at the i-th level receiver can be calculated. The propagation angle of the ray in the formation, the ray propagation path and the travel time of the first arrival; (2) When improving the shooting method for ray tracing, firstly, the launch angle θ of the initial position of the ray is given: The propagation of seismic waves between layers follows Snell's law, and it is judged whether the end position of the ray is the position of the monitoring point. If the end position of the ray is not the position of the monitoring point, the angle of the ray is continuously adjusted according to the dichotomy method, and recalculated until it meets the requirements; (3) Since the duration of the fracturing process is relatively long, it is not easy to determine the occurrence time of the earthquake source. Therefore, by calculating the first-arrival travel time difference between geophones, the seismic source time can be avoided, and this method can effectively reduce the degree of interference in the formation and improve the robustness of the inversion. 4. the source inversion method based on microseismic data and SPSA optimization algorithm as claimed in claim 1-3, is characterized in that, step 3 comprises the following steps: Generate a random disturbance composed of a group of random numbers that obey the Gaussian distribution within (0,1), and add the calculated first arrival travel time difference to the random disturbance to obtain the observed value. 5. the seismic source inversion method based on microseismic data and SPSA optimization algorithm as claimed in claim 1-4, is characterized in that, step 4 comprises the following steps: (1), first, set the initial parameters, calculate the first arrival travel time according to step two, and then obtain the difference between the first arrival travel time, and then calculate the difference between the observed data and the calculated data; (2), update the iteration parameters, also calculate the first arrival travel time, the difference between the first arrival travel time and the difference between the calculated data and the observed data according to the calculation described in step 2, and select the search iteration direction; (3), calculate according to the selected iterative direction, obtain the updated microseismic event position; (4) Calculate the updated first arrival travel time according to step two, and then obtain the difference between the updated first arrival travel time and the difference between the calculated data and the observed data; (5) Judging whether the objective function meets the accuracy, if the convergence condition is met, stop updating; otherwise, select a new search direction to continue iterative calculation.