Reservoir velocity frequency dispersion curve inversion method and device, equipment and storage medium

Abstract

The invention provides a reservoir velocity dispersion curve inversion method and device, and the method comprises the steps: calculating a single-interface reflection coefficient through employing a Zoeppritz equation, and carrying out the recursion of a layered medium reflection coefficient through combining with a reflectivity method; establishing an AVAF forward modeling equation based on the frequency domain convolution model; establishing an inversion objective function by using the difference between the frequency spectrum real part and the frequency spectrum imaginary part of the actual seismic record and the synthetic seismic record, and determining the upper limit and the lower limit of a to-be-inverted parameter; the method comprises the following steps: constructing a Metropolis-hot bath inversion framework by fusing an acceptance-rejection criterion of a Metropolis algorithm and a one-by-one parameter optimization thought of a hot bath algorithm, and carrying out multiple times of global optimization inversion based on an improved method to obtain an optimal reservoir thickness and speed frequency dispersion curve. According to the method provided by the invention, the reservoir thickness and the velocity values under different frequencies in the seismic frequency band can be inversely estimated, and a powerful tool is provided for reservoir prediction and fluid identification.

Description

Technical Field

[0001] The invention relates to the field of exploration geophysics, and in particular to a reservoir velocity dispersion curve inversion method, device, equipment and storage medium. Background Art

[0002] Existing rock physics theory shows that the velocity dispersion caused by seismic waves is a complex function of the properties of pore fluids, saturation, fluid mobility and their spatial distribution. When rocks are saturated with different fluids, their velocity dispersion characteristics will be significantly different. In fact, a large number of rock physics experimental studies have shown that the magnitude of the seismic wave reflection amplitude is not only related to the incident angle, but also closely related to the frequency. Especially in hydrocarbon-bearing reservoir rocks, the presence of pore fluids will cause seismic waves to disperse to varying degrees, which provides an important basis for using seismic dispersion attributes for reservoir prediction and fluid identification.

[0003] In recent years, the great potential of velocity dispersion attributes in oil and gas exploration has attracted widespread attention from scholars at home and abroad, and the dispersion AVO reservoir prediction technology with dispersion gradient as the hydrocarbon indicator factor has developed rapidly. This technology is based on the classical AVO theory, which introduces the frequency term to describe the influence of velocity dispersion, and combines the spectral decomposition technology to achieve dispersion gradient inversion, thereby effectively making up for the shortcomings of the classical AVO technology. The dispersion AVO technology was originally proposed based on the Zoeppritz approximation formula. It introduced the frequency term to add the influence of velocity dispersion to the classical AVO theory, established the dispersion gradient inversion equation, and laid an important foundation for the dispersion AVO inversion theory. Some people have combined the Wigner-Ville distribution spectral decomposition technology with the dispersion AVO technology to perform dispersion gradient inversion on actual seismic data, or apply the dispersion AVO technology to the quantitative prediction of reservoir gas saturation and porosity based on the Bayesian method. Furthermore, the velocity term is included in the dispersion gradient as a parameter to be inverted, and the dispersion AVO equation is re-derived using the Zoeppritz approximation. Compared with the dispersion AVO technical method, this method can still obtain accurate dispersion gradients under the condition of unknown velocity field, which is convenient for practical application. In addition, frequency-variant AVO inversion based on the Russell reflection coefficient approximation can construct a new frequency-variant fluid factor. The existing technology proposes a high-resolution sparse constrained inversion spectral decomposition method, and uses this method to extract reliable time-frequency phase information and amplitude information.

[0004] At present, the effectiveness and practicality of the dispersion AVO technology have been confirmed, and it has also been initially applied in actual production. However, the existing dispersion AVO inversion methods still have a series of difficult-to-overcome problems: first, it is difficult to explain physically. The dispersion AVO inversion only introduces the influence of frequency from a mathematical perspective and does not have a clear physical meaning. Second, the inversion method and the obtained results have poor mobility. The dispersion AVO inversion can only obtain the dispersion gradient at the reference frequency, and cannot quantitatively estimate the velocity values ​​at different frequencies. Therefore, its inversion results are difficult to be linked to reservoir physical parameters such as porosity and fluid saturation, and directly applied to the quantitative prediction of reservoir physical parameters. In addition, the dispersion AVO inversion is also affected by the thin layer tuning effect, and has problems such as low inversion accuracy and difficulty in distinguishing thin reservoirs. Therefore, it is very necessary to explore and develop velocity dispersion curve inversion algorithms. Summary of the invention

[0005] In view of this, in order to solve the problems of poor technical portability, low precision and difficulty in distinguishing thin reservoirs in the inversion methods in the prior art, the present invention provides a reservoir velocity dispersion curve inversion method, device, equipment and storage medium to achieve the inversion and quantitative estimation of thin reservoir velocity dispersion curve and reservoir thickness.

[0006] The first embodiment of the present invention provides a reservoir velocity dispersion curve inversion method, characterized in that the method comprises the following steps: Calculate single interface reflection coefficient; Recursive calculation of reflection coefficient of layered media; Establish AVAF forward equation; Determine the prior information of the inversion objective function and model parameters; The Metropolis-heat bath algorithm inversion framework is constructed and reservoir thickness and velocity dispersion curves are estimated.

[0007] Specifically, the process of calculating the single interface reflection coefficient includes: According to the elastic wave theory, the Zoeppritz equation under the condition of downward P wave incidence is used to solve the single interface reflection and transmission coefficients under the condition of downward P wave incidence: According to the elastic wave theory, the Zoeppritz equation under the condition of downward P wave incidence is used to solve the single interface reflection and transmission coefficients under the condition of downward P wave incidence:

[0008] in, denote the single interface reflection coefficient and transmission coefficient, respectively. The upper index It indicates that the incident wave is a downgoing wave, the first subscript indicates the wave type of the incident wave, and the second subscript indicates the wave type of the converted wave; Given by:

[0009] and are the vertical slowness of P-wave and S-wave, respectively. , "0" and "1" represent medium 0 above the incident interface and medium 1 below the incident interface, respectively. is the medium density, is the shear modulus, are the media where the incident wave is located P-wave velocity and shear-wave velocity; is the horizontal slowness, where is the P-wave incident angle; The Zoeppritz equation for the case of downlink S-wave incidence is used to solve the single-interface reflection and transmission coefficients for the case of downlink S-wave incidence:

[0010] The Zoeppritz equation for the case of upgoing P-wave incidence is used to solve the single-interface reflection and transmission coefficients for the case of upgoing P-wave incidence:

[0011] in, The upper index It indicates that the incident wave is an upgoing wave. The first subscript indicates the wave type of the incident wave, and the second subscript indicates the wave type of the converted wave. Given by:

[0012] Using the Zoeppritz equation for the case of upward S-wave incidence, the single-interface reflection and transmission coefficients for the case of upward S-wave incidence are solved as follows:

[0013] Specifically, the process of recursively calculating the reflection coefficient of the layered medium includes: According to the reflectivity method, the reflection coefficient matrix of the top interface of the model is The reflection coefficient matrix of the bottom layer medium Step by step recursion upwards we can obtain:

[0014] in, Indicates the total number of layers of layered media; For the The phase shift matrix of the layer is given by: , is the angular frequency, For the The thickness of the layer, Respectively P-wave and S-wave vertical slowness of the layer; matrix elements Interface The single interface reflection coefficient and transmission coefficient at Respectively represent the interface The reflection coefficient matrix and transmission coefficient matrix when the upgoing wave is incident at Respectively represent the interface The reflection coefficient matrix and transmission coefficient matrix for the downlink wave incident at are given by: , ; Since the velocity of a fluid-containing reservoir is usually frequency-dependent, the reflection coefficient matrix of the top interface of the model is It’s also about angles and angular frequency Function matrix of: .

[0015] Specifically, the process of establishing the AVAF forward equation includes: According to the reflection coefficient matrix of the top interface of the model , get the PP wave reflection coefficient ; According to the frequency domain convolution model, the synthetic data of the time domain PP wave angle gather is obtained by inverse Fourier transform: , in, is the frequency domain seismic wavelet, For time; By adjusting the phase shift matrix, the reflection time of the non-zero angle gather and the zero angle gather are unified, and the phase axis of the angle gather is flattened, and the dynamic correction processing of the synthetic data is realized to match the pre-stack migration or dynamic correction that has been completed for the actual angle gather data; the phase shift matrix is ​​adjusted to .

[0016] Specifically, the process of determining the prior information of the inversion objective function and model parameters includes: The objective function is defined by the difference between the real and imaginary parts of the spectrum of the synthetic data and the actual data: , in, is the actual observed data vector; is the forward synthetic data vector; Spectra recorded for earthquakes; For vector Length; Using the well logging data and core test data of the target area as inversion constraints, determine the upper limit of the model parameters. and lower limit .

[0017] Specifically, the process of constructing the Metropolis-heat bath algorithm inversion framework and estimating the velocity dispersion curve and reservoir thickness includes: Analyze the inversion framework of the Metropolis algorithm, including: Analyze the method of obtaining a new model parameter vector in the classic Metropolis algorithm, which is based on the initial model Apply some disturbance To obtain the new model parameter vector: ; Analysis of the acceptance-rejection rule of the Metropolis algorithm: According to The size of the model determines whether to accept the new model. , accept the new model, otherwise, the probability of accepting the new model is ,in, Based on the model parameters The error between the synthetic seismic record obtained by forward modeling and the actual seismic record, is temperature; The inversion framework for the analytical heat bath algorithm includes: Given an initial temperature And the range of model parameters and step length ; Let the model parameter vector be , and the model parameters The possible values ​​of , any model parameter value in the model space is expressed as ,in and , randomly select model parameters from the model space and build the initial model ; Perform parameter-by-parameter optimization by visiting the model each time All possible values ​​of a model parameter in , determine its best value, while keeping other parameters that have not been visited unchanged, and determine the best value of each model parameter in the same way, and finally get a new model The process includes: fixing Other model parameters , select , and calculate the corresponding marginal probability: , And the cumulative probability distribution: ; In uniform distribution Pick a random number from , select Closest No. Model parameter values Replace the original Take the value, let ; At the same temperature, using the same method, determine the model parameters in turn The latest value of ; Use the acceptance-rejection rule of the Metropolis algorithm to determine whether to accept the new model If accepted, then the Assign to the ; The temperature is gradually lowered according to the cooling criteria of the heat bath algorithm. The above process is repeated until the accuracy requirements are met or the maximum number of iterations is reached, and the final algorithm model is output. ; By integrating the acceptance-rejection criterion of the Metropolis algorithm and the idea of ​​parameter-by-parameter optimization of the heat bath algorithm, an inversion framework of the Metropolis-heat bath algorithm is constructed to form an algorithm flow, and the algorithm flow is used to estimate the velocity dispersion curve and reservoir thickness.

[0018] Preferably, the process of estimating the velocity dispersion curve and the reservoir thickness using the algorithm flow comprises: Set the maximum number of iterations for the inversion calculation to ; Extract the target reservoir segment information from the seismic profile data and extract the seismic wavelet ; The information of the target reservoir segment is Fourier transformed to obtain the actual observed data vector in the frequency domain ; Obtain the model parameter vector consisting of the parameters to be inverted: the target reservoir segment is composed of The medium is composed of layers, and the parameter describing the medium model is the P wave velocity of the top elastic layer medium , S wave speed and density , the frequency-dependent P-wave velocity of the intermediate dispersive layer medium , S wave speed ,density and thickness , and the P wave velocity of the bottom elastic layer medium , S wave speed and density , then the model parameter vector composed of these parameters to be inverted is ,in, A series of discrete points , , …, Interpolation is obtained, and the parameters Expands to: , , …, ; Given an initial temperature And the search lower limit of the model parameters and upper limit , and take the initial number of iterations ; Randomly select model parameters within the upper and lower limits of the model parameters to generate an initial model vector ; Using AVAF to synthesize seismic data , calculate its difference with the actual data The error function between ; At temperature Next, we search for the best parameters one by one according to the idea of ​​the heat bath algorithm, and decide whether to accept the new model according to the Metropolis criterion. The process includes: ,Change The model parameters, namely ,in To satisfy the uniform distribution of random numbers, the updated ;calculate ,like ,or and , then accept the new model and take and ,in is another random number that satisfies the uniform distribution; otherwise, reject the new model; if ,but , and repeat the above steps until all model parameters are visited; According to the cooling guidelines , gradually lowering the temperature, is the cooling rate; like , , repeat the above process until the maximum number of iterations is reached; Different initial models are randomly selected and the above process is repeated to obtain a series of inversion results. Finally, the inversion results are averaged to obtain the final model parameter values, and the inverted reservoir thickness and velocity dispersion curves are further obtained.

[0019] The present invention also protects a reservoir velocity dispersion curve inversion device, which uses the steps of the above method to perform reservoir velocity dispersion curve inversion, and the device includes the following modules: The first module: used to calculate the single interface reflection coefficient; The second module: used to recursively calculate the reflection coefficient of layered media; The third module: used to establish the AVAF forward equation; Module 4: Prior information used to determine the inversion objective function and model parameters; The fifth module is used to build the Metropolis-heat bath algorithm inversion framework and invert the reservoir thickness and velocity dispersion curves.

[0020] The present invention also protects a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps of the above-mentioned reservoir velocity dispersion curve inversion method when executing the computer program.

[0021] On the other hand, the present invention protects a storage medium having a computer program stored thereon, which implements the steps of the aforementioned reservoir velocity dispersion curve inversion method when executed by a processor.

[0022] In summary, the present invention proposes a reservoir velocity dispersion curve inversion method. Compared with the prior art, the method of the present invention combines two classic simulated annealing methods, the Metropolis algorithm and the heat bath algorithm. By making full use of the high efficiency of the acceptance-rejection criterion of the Metropolis algorithm and the advantage of the heat bath algorithm in parameter-by-parameter optimization, the stable inversion and calculation of the velocity dispersion curve and thickness of the thin reservoir are achieved. While improving the technical portability of the inversion method, the accuracy of the reservoir velocity dispersion curve inversion is mainly improved, and thin reservoirs can be more easily distinguished. The method has a good application and promotion prospect. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0024] Figure 1 This is a flow chart of the reservoir velocity dispersion curve inversion method in the first embodiment of this specification; Figure 2 Schematic diagram of the structure of the reservoir medium model under the P wave incident condition and the S wave incident condition in the first embodiment of the present invention, wherein: Indicates the total number of layers in the layered medium. The model parameters of each layer include: P wave velocity , S wave speed ,density and the thickness of the layered dielectric , ; Figure 3 Schematic diagram of single interface reflection and transmission in the first embodiment of the present invention, wherein ① represents the incidence and refraction of a downstream incident P wave from medium 0 into medium 1, ② represents the incidence and refraction of a downstream incident S wave from medium 0 into medium 1, ③ represents the incidence and refraction of an upstream incident P wave from medium 1 into medium 0, and ④ represents the incidence and refraction of an upstream incident S wave from medium 1 into medium 0; Figure 4 is an inversion flow chart of the Metropolis-heat bath algorithm in the first embodiment of this specification; Figure 5 Schematic diagram of a three-layer dielectric model structure and a P-wave velocity curve in the second embodiment of this specification, wherein (a) is a schematic diagram of the three-layer dielectric model structure, and (b) is a schematic diagram of the P-wave velocity curve corresponding to the second dielectric layer of the three-layer dielectric model structure; Figure 6 The angle gather data of the three-layer medium model in the second embodiment of this specification; Figure 7 This is a curve diagram showing the variation of the error function of the Metropolis-heat bath algorithm with the number of iterations in the second embodiment of this specification; Figure 8 : is a comparison diagram of the dispersion layer parameters obtained by inversion based on the Metropolis-heat bath algorithm in the second embodiment of this specification and the actual model parameters, wherein (i) is a comparison diagram of the P-wave velocity of the second layer, (ii) is a comparison diagram of the thickness of the layered medium layer, (iii) is a comparison diagram of the density of the second layer, and (iv) is a comparison diagram of the S-wave velocity of the second layer; Fig. 9 This is a comparison chart of the observed data and the simulated data calculated based on the inversion results of the optimal model parameters in the second embodiment of this specification. DETAILED DESCRIPTION

[0025] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0026] The present invention provides a reservoir velocity dispersion curve inversion method, device, equipment and storage medium, which better solve the problems of poor technical mobility, low precision and difficulty in distinguishing thin reservoirs in the inversion method in the prior art.

[0027] In the first embodiment, referring to Figure 1 The present invention proposes a reservoir velocity dispersion curve inversion method, comprising the following steps: Calculate single interface reflection coefficient; Recursive calculation of reflection coefficient of layered media; Establish AVAF forward equation; Determine the prior information of the inversion objective function and model parameters; The Metropolis-heat bath algorithm inversion framework is constructed and reservoir thickness and velocity dispersion curves are estimated.

[0028] The general reservoir medium model contains multiple layered medium layers, such as Figure 2 As shown, use Indicates the total number of layers of layered media. Different layers of layered media are connected by interfaces. The model parameters include: P-wave (longitudinal wave) velocity , S wave (transverse wave) speed ,density and the thickness of the layered dielectric , . Figure 2 Two cases of S-wave and P-wave incident on layered media are shown. In order to explore the physical principles behind and use their laws, we can first discuss the reflection and transmission of different wave types on a single interface (hereinafter referred to as a single interface), obtain the reflection coefficient of the corresponding single interface, and then generalize it to multiple interfaces and multiple layered media layers, and recursively calculate the reflection coefficient of the layered media.

[0029] Specifically, the calculation of the single interface reflection coefficient may include the following steps: For a single interface (single interface) in a model, such as Figure 3 As shown, this can be done by solving the Zoeppritz equation for the case of downlink P-wave incidence: , (1) Calculate the reflection coefficient and transmission coefficient : (2) in, denote the single interface reflection coefficient and transmission coefficient, respectively, and their superscripts indicates that the incident wave is a downlink wave, the first subscript indicates the wave type of the incident wave, and the second subscript indicates the wave type of the converted wave; further, Given by: (3) In the above formula, and are the vertical slowness of P-wave and S-wave, respectively. , “0” and “1” represent the upper medium (medium 0) and lower medium (medium 1) of the incident interface, respectively. is the medium density, is the shear modulus, where are the media where the incident wave is located P-wave velocity and shear-wave velocity; is the horizontal slowness, where is the P-wave incident angle; By solving the Zoeppritz equation for the case of downlink S-wave incidence: , (4) Calculate the reflection coefficient and transmission coefficient : (5) Similarly, by solving the Zoeppritz equation for the case of upgoing P-wave incidence: , (6) Calculate the reflection coefficient and transmission coefficient : (7) Furthermore, by solving the Zoeppritz equation for the case of upgoing S-wave incidence: (8) Calculate the reflection coefficient and transmission coefficient : (9) in: (10) Specifically, the recursive calculation of the layered medium reflection coefficient includes the following steps: For Figure 2 In the layered medium model shown in the figure, when the P wave When the incident angle is downward into the layered medium, according to the reflectivity method, by stepping up recursion, the reflection coefficient matrix of the top interface of the model is The reflection coefficient matrix of the bottom layer medium can be Calculation results: (11) in, For the The phase shift matrix of the layer is , (12) is the angular frequency; For the Thickness of the layer; Respectively P-wave and S-wave vertical slowness of the layer; matrix elements Interface The single interface reflection coefficient and transmission coefficient at Respectively represent the interface The reflection coefficient matrix and transmission coefficient matrix when the upgoing wave is incident at ; Respectively represent the interface The reflection coefficient matrix and transmission coefficient matrix when the downlink wave is incident; that is , (12) , (13) Among them, the matrix elements are interfaces The single interface inversion coefficient and transmission coefficient at are calculated by equations (2), (5), (7) and (9), where the superscript “ "and" " respectively represent the upgoing wave incident situation and the downgoing wave incident situation, the first letter of the subscript represents the incident wave type, and the second letter represents the converted wave type.

[0030] In addition, since the velocity of a fluid-bearing reservoir is usually frequency-dependent, the reflection coefficient matrix of the top interface of the model is It’s also about angles and angular frequency The function matrix of . (14) Specifically, the establishment of the AVAF forward equation includes the following steps: First, according to the calculated reflection coefficient matrix , get the PP wave reflection coefficient .

[0031] Then, the frequency domain convolution model is used to obtain the time domain PP wave angle trace composite data: , (15) In the above formula, is the frequency domain seismic wavelet; Considering that the actual angle gather data usually has completed pre-stack migration or dynamic correction, in order to match it, the synthetic data also needs to be processed by dynamic correction, so the phase shift matrix needs to be adjusted to , (16) This unifies the reflection times of the non-zero angle gathers and the zero angle gathers, and flattens the event axes of the angle gathers.

[0032] Specifically, the process of determining the prior information of the inversion objective function and the model parameters includes the following steps: The ultimate goal of this method is to find a model parameter vector , so that the synthetic data forward modeled by it Compared with the actual observation data The error between In order to make full use of the information of different frequency components of seismic data, the difference between the real and imaginary parts of the spectrum of synthetic data and actual data is used to define the objective function: , (17) In the above formula, is the actual observed data vector; is the forward synthetic data vector; Spectra recorded for earthquakes; For vector In addition, in order to reduce the multi-solution of inversion and improve the inversion accuracy, the upper limit of the model parameters can be determined by using the well logging data and core test data of the target area as inversion constraints. and lower limit .

[0033] Specifically, the process of constructing the Metropolis-heat bath algorithm inversion framework and estimating the velocity dispersion curve and reservoir thickness includes the following steps: First, the inversion framework of the Metropolis algorithm is analyzed. Suppose the vector consisting of the model parameters to be inverted is In the classic Metropolis algorithm, the initial model is usually Apply some disturbance To obtain the new model parameter vector: . (18) Based on the model parameters The error between the synthetic seismic record obtained by forward modeling and the actual seismic record is , then the key to the Metropolis algorithm is: according to The size of determines whether to accept the new model, that is, if , then the new model is accepted; otherwise, the probability of accepting the new model is , (19) In the above formula, is the temperature. When the temperature is high, the probability of the Metropolis method accepting a new model is high; as the temperature decreases, the probability of accepting a new model gradually decreases until the temperature approaches zero, and the Metropolis algorithm degenerates into a greedy algorithm, at which point new models are almost no longer accepted. Therefore, in the Metropolis algorithm, many models will be rejected, especially when the temperature is low, the probability of rejecting a new model is very high.

[0034] Then, the inversion framework of the heat bath algorithm is analyzed. The heat bath algorithm avoids extremely high rejection probabilities by calculating the relative acceptance probabilities of all possible values ​​of each model parameter. Its characteristics are: at each temperature, the algorithm will sequentially visit all possible values ​​of each model parameter to determine the best new model. The specific search process is as follows: 1) Given the initial temperature And the range of model parameters and step length ; 2) Let the model parameter vector be , and the model parameters The possible values ​​of Therefore, any model parameter value in the model space can be expressed as ,in and , randomly select model parameters from the model space and build the initial model ; 3) While keeping other model parameters unchanged, visit all possible values ​​of each model parameter in turn and determine its optimal value. For example: First, fix Other model parameters ; Then, select , and calculate the corresponding marginal probability: , (20) And the cumulative probability distribution: .(twenty one) In uniform distribution Pick a random number from , select Closest No. Model parameter values Replace the original Take the value, that is At the same temperature, the model parameters are determined in the same way. The latest value of .

[0035] Finally, the acceptance-rejection criterion of the Metropolis algorithm and the idea of ​​the heat bath algorithm to optimize each parameter at each temperature are fully integrated to construct Figure 4 The Metropolis-heat bath algorithm inversion process shown in Figure 4 improves the computational efficiency and stability of the inversion.

[0036] Inverting reservoir thickness and velocity dispersion curves based on the improved algorithm can include the following steps: (a) Set the maximum number of iterations for the inversion calculation to ; (b) Extract the target reservoir segment from the seismic profile data in the actual seismic record and extract the seismic wavelet ; (c) Perform Fourier transform on the seismic records of the target reservoir segment to obtain the actual observed data vector in the frequency domain ; (d) The target reservoir segment is composed of The medium is composed of layers, and the parameter describing the medium model is the P wave velocity of the top elastic layer medium , S wave speed and density , the frequency-dependent P-wave velocity of the intermediate dispersive layer medium , S wave speed ,density and thickness , and the P wave velocity of the bottom elastic layer medium , S wave speed and density , then the model parameter vector composed of these parameters to be inverted is , among which due to A series of discrete points , , …, Interpolation is obtained, and the parameters can be Expand to , , …, ; (e) Given the initial temperature And the search lower limit of the model parameters and upper limit , and take the current number of iterations as ; (f) Randomly select model parameters within the upper and lower limits of the model to generate the initial model vector , and synthesized seismic data using AVAF forward modeling , calculate its difference with the actual data The error function between ; (g) at temperature Under this condition, we search for the best parameters one by one according to the idea of ​​the heat bath algorithm, and decide whether to accept the new model according to the Metropolis criterion. The specific steps are as follows: ① Let ,Change The ( ) model parameters, that is, ,in To satisfy the uniform distribution of random numbers, the updated ; ② Calculation ,like ,or and , then accept the new model and take and ,in is another random number that satisfies uniform distribution; otherwise, reject the new model; ③If ,but , and repeat the above steps until all model parameters are visited; (h) According to the cooling criteria , gradually lowering the temperature, is the cooling rate; (i) If , , repeat the above process until the maximum number of iterations is reached; (j) Randomly select different initial models and repeat the above process to obtain a series of inversion results. Finally, average the inversion results to obtain the optimal model parameter values, and further obtain the inverted reservoir thickness and velocity dispersion curves.

[0037] In this embodiment, the present invention proposes a reservoir velocity dispersion curve inversion method that integrates the Metropolis algorithm and the heat bath algorithm based on the AVAF forward equation constructed by the reflectivity method. The method improves the Metropolis algorithm by combining the idea of ​​the heat bath algorithm. On the premise of maintaining the acceptance-rejection criterion of the Metropolis algorithm, the advantage of the heat bath algorithm in parameter-by-parameter optimization is integrated to construct a simulated annealing inversion framework based on the Metropolis-heat bath algorithm, thereby improving the computational efficiency and stability of the reservoir velocity dispersion curve inversion.

[0038] In the second embodiment of the present invention, a three-layer medium model is taken as an example to further illustrate the steps of the reservoir velocity dispersion curve inversion method based on the Metropolis-heat bath algorithm of the present invention, and verify the effect of the method of the first embodiment. , the top and bottom layers are elastic layers, such as Figure 5 As shown in (i), the P wave velocity at the top layer is , the S wave speed is , the density is ; The P-wave velocity of the bottom layer is , the S wave speed is , the density is ; The middle layer is the dispersion layer, with a thickness of 20 m and a P-wave velocity of As frequency changes Figure 5 As shown in (ii), the S wave speed is , the density is The angle gather data corresponding to the three-layer model is as follows: Figure 6 To test the reliability of the Metropolis-heat bath algorithm, the P-wave velocity, S-wave velocity, density and thickness of the dispersion layer of the model are assumed to be unknown parameters to be inverted, and its synthetic angle gathers are used as actual observation data to carry out global optimization inversion.

[0039] Specifically, in the process of inverting the reservoir thickness and velocity dispersion curve based on the improved simulated annealing algorithm, the step (c) includes: assuming that the target reservoir segment is composed of three layers of media, and the parameter describing the media model is the P wave velocity of the top elastic layer medium , S wave speed and density , the frequency-dependent P-wave velocity of the intermediate dispersive layer medium , S wave speed ,density and thickness , and the P wave velocity of the bottom elastic layer medium , S wave speed and density ,in A series of discrete points , , …, The model parameter vector composed of these parameters to be inverted is .

[0040] For the three-layer model described in this embodiment, during the inversion process, the seismic wavelet is taken as the Ricker wavelet with a main frequency of 30 Hz, and the initial temperature of the Metropolis-heat bath algorithm is set to The cooling rate is , the maximum number of iterations is , and set the search range of elastic layer longitudinal wave velocity to [5.5, 6.5] km / s, the search range of shear wave velocity to [3.6, 4.4] km / s, and the search range of density to [2.2, 2.8] g / cm 3 ; The search range for the dispersion layer P-wave velocity is [3.5, 4.5] km / s, the S-wave velocity is [2.2, 3.0] km / s, and the density is [2.2, 2.8] g / cm 3 , thickness is [10, 30] m. Figure 7 The error function of the Metropolis-heat bath algorithm changes with the number of iterations. Figure 7 It can be seen that when the number of iterations is greater than 400, the Metropolis-heat bath algorithm has basically reached a stable state; thereafter, with the continued increase in the number of iterations or the further decrease in temperature, the error function hardly changes significantly. Figure 8 The inversion results of the dispersion layer model parameters obtained by 100 random inversions based on the Metropolis-heat bath algorithm are shown in Figure 2. Figure 8 The two parallel dashed lines in (i)-(iv) represent the search range of each corresponding model parameter. The dashed line parallel to the true solution curve of the model parameter in the figure represents the average value of the model parameter for 100 inversion results, that is, the optimal inversion result. Figure 8 It can be seen that the average values ​​of the P-wave velocity dispersion curve and thickness of the dispersion layer inverted by the Metropolis-heat bath algorithm are in the best agreement with the true values, while the accuracy of the S-wave velocity and density is slightly inferior, with inversion errors of 2.7% and 3.2% respectively. Fig. 9 The real and imaginary part curves of the spectrum of the simulated data calculated based on the optimal inversion result and the actual observed data are compared. It can be seen from the figure that the simulated data and the actual data are in very good agreement, with an error of only 2.3×10 -5 , which proves the reliability of the Metropolis-heat bath inversion algorithm.

[0041] On the other hand, in a third embodiment, the present invention protects a reservoir velocity dispersion curve inversion device, the device comprising the following modules: The first module: used to calculate the single interface reflection coefficient; The second module: used to recursively calculate the reflection coefficient of layered media; The third module: used to establish the AVAF forward equation; Module 4: Prior information used to determine the inversion objective function and model parameters; Module 5: Used to build the Metropolis-heat bath algorithm inversion framework and estimate reservoir thickness and velocity dispersion curves.

[0042] For the convenience of description, the above device is described in various units according to their functions. Of course, when implementing the technical solution recorded in the present invention, the functions of each unit can be implemented in the same or multiple software and / or hardware.

[0043] Compared with the prior art, the present invention provides a reservoir velocity dispersion curve inversion method, which combines two classic simulated annealing methods, the Metropolis algorithm and the heat bath algorithm. By making full use of the high efficiency of the acceptance-rejection criterion of the Metropolis algorithm and the advantage of the heat bath algorithm in parameter-by-parameter optimization, the stable inversion and calculation of the velocity dispersion curve and thickness of the thin reservoir are achieved, solving the problems of poor technical transferability, low precision and difficulty in distinguishing thin reservoirs in the inversion method in the prior art.

[0044] The present invention also provides a computer device in one embodiment, including a memory and a processor, the memory storing a computer program, and the processor implementing the steps of the reservoir velocity dispersion curve inversion method provided in any of the above embodiments when executing the computer program. The computer device may be a server. The computer device includes a processor, a memory, a network interface and a database connected via a system bus. Among them, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used to store sample data. The network interface of the computer device is used to communicate with an external terminal via a network connection.

[0045] Furthermore, in another embodiment, the present invention provides a computer-readable storage medium having a computer program stored thereon, and when the computer program is executed by a processor, the steps of the reservoir velocity dispersion curve inversion method provided in any of the above embodiments are implemented.

[0046] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0047] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0048] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0049] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0050] The memory may include non-permanent storage in a computer-readable medium, random access memory (RAM) and / or non-volatile memory in the form of read-only memory (ROM) or flash RAM. The memory is an example of a computer-readable medium.

[0051] Computer readable media include permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. Information can be computer readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, compact disk read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cassettes, disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device. As defined herein, computer readable media does not include temporary computer readable media (transitory media), such as modulated data signals and carrier waves.

[0052] It should also be noted that the terms "comprises", "includes" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, or device. In the absence of further restrictions, an element defined by the sentence "comprises a ..." does not exclude the presence of other identical elements in the process, method, or device including the element.

[0053] It should be understood by those skilled in the art that the embodiments of this specification may be provided as methods, systems or computer program products. Therefore, this specification may take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware. Moreover, this specification may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0054] This specification may be described in the general context of computer-executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform specific tasks or implement specific abstract data types. This specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices connected through a communication network. In a distributed computing environment, program modules may be located in local and remote computer storage media, including storage devices.

[0055] Matters not covered by the present invention are known technologies.

[0056] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0057] The above-mentioned embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be subject to the attached claims.

[0058] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A reservoir velocity dispersion curve inversion method, characterized in that: The steps of the method include: Calculate single interface reflection coefficient; Recursive calculation of reflection coefficient of layered media; Establish AVAF forward equation; Determine the prior information of the inversion objective function and model parameters; The Metropolis-heat bath algorithm inversion framework is constructed and reservoir thickness and velocity dispersion curves are estimated.

2. A reservoir velocity dispersion curve inversion method according to claim 1, characterized in that: The process of calculating the single interface reflection coefficient includes: According to the elastic wave theory, the Zoeppritz equation under the condition of downward P wave incidence is used to solve the single interface reflection and transmission coefficients under the condition of downward P wave incidence: in, denote the single interface reflection coefficient and transmission coefficient, respectively. The upper index It indicates that the incident wave is a downgoing wave, the first subscript indicates the wave type of the incident wave, and the second subscript indicates the wave type of the converted wave; Given by: and are the vertical slowness of P-wave and S-wave, respectively. , "0" and "1" represent medium 0 above the incident interface and medium 1 below the incident interface, respectively. is the medium density, is the shear modulus, are the media where the incident wave is located P-wave velocity and shear-wave velocity; is the horizontal slowness, where is the P-wave incident angle; The Zoeppritz equation for the case of downlink S-wave incidence is used to solve the single-interface reflection and transmission coefficients for the case of downlink S-wave incidence: The Zoeppritz equation for the case of upgoing P-wave incidence is used to solve the single-interface reflection and transmission coefficients for the case of upgoing P-wave incidence: in, The upper index It indicates that the incident wave is an upgoing wave. The first subscript indicates the wave type of the incident wave, and the second subscript indicates the wave type of the converted wave. Given by: Using the Zoeppritz equation for the case of upward S-wave incidence, the single-interface reflection and transmission coefficients for the case of upward S-wave incidence are solved as follows: 。 3. A reservoir velocity dispersion curve inversion method as claimed in claim 2, characterized in that: The process of recursively calculating the reflection coefficient of the layered medium includes: According to the reflectivity method, the reflection coefficient matrix of the top interface of the model is The reflection coefficient matrix of the bottom layer medium Step by step recursion upwards we can obtain: in, Indicates the total number of layers of layered media; For the The phase shift matrix of the layer is given by: ， is the angular frequency, For the The thickness of the layer, Respectively P-wave and S-wave vertical slowness of the layer; matrix elements Interface The single interface reflection coefficient and transmission coefficient at Respectively represent the interface The reflection coefficient matrix and transmission coefficient matrix when the upgoing wave is incident at Respectively represent the interface The reflection coefficient matrix and transmission coefficient matrix for the downlink wave incident at are given by: , ； Since the velocity of a fluid-containing reservoir is usually frequency-dependent, the reflection coefficient matrix of the top interface of the model is It’s also about angles and angular frequency Function matrix of: 。 4. A reservoir velocity dispersion curve inversion method as claimed in claim 3, characterized in that: The process of establishing the AVAF forward equation includes: According to the reflection coefficient matrix of the top interface of the model , get the PP wave reflection coefficient ; According to the frequency domain convolution model, the synthetic data of the time domain PP wave angle gather is obtained by inverse Fourier transform: ， in, is the frequency domain seismic wavelet, For time; By adjusting the phase shift matrix, the reflection time of the non-zero angle gather and the zero angle gather are unified, and the phase axis of the angle gather is flattened, and the dynamic correction processing of the synthetic data is realized to match the pre-stack migration or dynamic correction that has been completed for the actual angle gather data; the phase shift matrix is ​​adjusted to 。 5. A reservoir velocity dispersion curve inversion method as claimed in claim 4, characterized in that: The process of determining the prior information of the inversion objective function and the model parameters includes: The objective function is defined by the difference between the real and imaginary parts of the spectrum of the synthetic data and the actual data: ， in, is the actual observed data vector; is the forward synthetic data vector; Spectra recorded for earthquakes; For vector Length; Using the well logging data and core test data of the target area as inversion constraints, determine the upper limit of the model parameters. and lower limit .

6. A reservoir velocity dispersion curve inversion method as claimed in claim 5, characterized in that: The process of constructing the Metropolis-heat bath algorithm inversion framework and estimating the velocity dispersion curve and reservoir thickness includes: Analyze the inversion framework of the Metropolis algorithm, including: Analyze the method of obtaining a new model parameter vector in the classic Metropolis algorithm, which is based on the initial model Apply some disturbance To obtain the new model parameter vector: ; Analysis of the acceptance-rejection rule of the Metropolis algorithm: According to The size of the model determines whether to accept the new model. , accept the new model, otherwise, the probability of accepting the new model is ,in, Based on the model parameters The error between the synthetic seismic record obtained by forward modeling and the actual seismic record, is temperature; The inversion framework for the analytical heat bath algorithm includes: Given an initial temperature And the range of model parameters and step length ; Let the model parameter vector be , and the model parameters The possible values ​​of , any model parameter value in the model space is expressed as ,in and , randomly select model parameters from the model space and build the initial model ; Each time the model is accessed All possible values ​​of a model parameter in , determine its best value, while keeping other parameters that have not been visited unchanged, and determine the best value of each model parameter in the same way, and finally get a new model The process includes: fixing Other model parameters , select , and calculate the corresponding marginal probability: ， And the cumulative probability distribution: ； In uniform distribution Pick a random number from , select Closest No. Model parameter values Replace the original Take the value, let ; At the same temperature, using the same method, determine the model parameters in turn The latest value of ; Use the acceptance-rejection rule of the Metropolis algorithm to determine whether to accept the new model If accepted, then the Assign to the ; According to the cooling criteria, the temperature is gradually lowered, and the above process is repeated until the accuracy requirements are met or the maximum number of iterations is reached, and the final algorithm model is output. ; By integrating the acceptance-rejection criterion of the Metropolis algorithm and the idea of ​​parameter-by-parameter optimization of the heat bath algorithm, an inversion framework of the Metropolis-heat bath algorithm is constructed to form an algorithm flow, and the algorithm flow is used to estimate the velocity dispersion curve and reservoir thickness.

7. A reservoir velocity dispersion curve inversion method as claimed in claim 6, characterized in that: The process of estimating the velocity dispersion curve and reservoir thickness using the algorithm flow includes: Set the maximum number of iterations for the inversion calculation to ; Extract the target reservoir segment information from the seismic profile data and extract the seismic wavelet ; The information of the target reservoir segment is Fourier transformed to obtain the actual observed data vector in the frequency domain ; Obtain the model parameter vector consisting of the parameters to be inverted: the target reservoir segment is composed of The medium is composed of layers, and the parameter describing the medium model is the P wave velocity of the top elastic layer medium , S wave speed and density , the frequency-dependent P-wave velocity of the intermediate dispersive layer medium , S wave speed ,density and thickness , and the P wave velocity of the bottom elastic layer medium , S wave speed and density , then the model parameter vector composed of these parameters to be inverted is ,in, A series of discrete points , , …, Interpolation is obtained, and the parameters Expands to: 、 、…、 ; Given an initial temperature And the search lower limit of the model parameters and upper limit , and take the initial number of iterations ; Randomly select model parameters within the upper and lower limits of the model parameters to generate an initial model vector ; Using AVAF to synthesize seismic data , calculate its difference with the actual data The error function between ; At temperature Next, we search for the best parameters one by one according to the idea of ​​the heat bath algorithm, and decide whether to accept the new model according to the Metropolis criterion. The process includes: ,Change The model parameters, namely ,in To satisfy the uniform distribution of random numbers, the updated ;calculate ,like ,or and , then accept the new model and take and ,in is another random number that satisfies the uniform distribution; otherwise, reject the new model; if ,but , and repeat the above steps until all model parameters are visited; According to the cooling guidelines , gradually lowering the temperature, is the cooling rate; like , , repeat the above process until the maximum number of iterations is reached; Different initial models are randomly selected and the above process is repeated to obtain a series of inversion results. Finally, the inversion results are averaged to obtain the final model parameter values, and the inverted reservoir thickness and velocity dispersion curves are further obtained.

8. A reservoir velocity dispersion curve inversion device, characterized in that: The device uses the steps of the method of claim 1 to perform reservoir velocity dispersion curve inversion, and the device includes the following modules: The first module: used to calculate the single interface reflection coefficient; The second module: used to recursively calculate the reflection coefficient of layered media; The third module: used to establish the AVAF forward equation; Module 4: Prior information used to determine the inversion objective function and model parameters; The fifth module is used to build the Metropolis-heat bath algorithm inversion framework and invert the reservoir thickness and velocity dispersion curves.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the method according to any one of claims 1 to 7 are implemented.

10. A storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 7 are implemented.

Images