CN108645994B

CN108645994B - Geological random inversion method and device based on multipoint geostatistics

Info

Publication number: CN108645994B
Application number: CN201810379669.4A
Authority: CN
Inventors: 李景叶; 陈小宏; 刘兴业; 刘国昌; 马继涛; 王守东
Original assignee: China University of Petroleum Beijing
Current assignee: China University of Petroleum Beijing
Priority date: 2018-04-25
Filing date: 2018-04-25
Publication date: 2020-01-17
Anticipated expiration: 2038-04-25
Also published as: CN108645994A

Abstract

The invention provides a geological random inversion method and device based on multipoint geostatistics, and relates to the technical field of geological and petroleum exploration. The method comprises the steps of determining an initial lithofacies model through a multipoint geostatistical random simulation algorithm; carrying out probability disturbance updating to generate a lithofacies model after the probability disturbance updating; counting the distribution functions of the reservoir physical property parameters in different rock phases, and converting the rock phase model after the probability disturbance updating into a pseudo reservoir physical property parameter model; converting the pseudo reservoir physical property parameter model into an elastic parameter model; performing time-depth conversion to obtain a time domain elastic parameter model, determining a reflection coefficient, and forward modeling by adopting a convolution model to synthesize the seismic record; comparing the seismic record with a predetermined actual seismic record, and determining a lithofacies model result and a pseudo reservoir physical property parameter model result according to the maximum disturbance times and the maximum simulation times; and carrying out quantum annealing optimization on the pseudo reservoir physical property parameter model result to generate a reservoir physical property parameter inversion result.

Description

Geological random inversion method and device based on multipoint geostatistics

Technical Field

The invention relates to the technical field of geological and petroleum exploration, in particular to a geological random inversion method and device based on multi-point geostatistics.

Background

At present, reservoir physical parameters (including porosity, shale content, water saturation and the like) are important bases for reservoir prediction and evaluation, oil and gas reserve estimation and well position development determination. The seismic data contain abundant reservoir physical property information, have the characteristic of high transverse resolution, and are important data used in reservoir characterization and prediction, but because the relationship between the reservoir physical property parameters and the seismic data is abnormal and complex, great difficulty is brought to the acquisition of the reservoir physical property parameters by using the seismic data. Seismic inversion techniques can play an important role in reservoir characterization and modeling. In general, seismic inversion can be divided into deterministic inversion and stochastic inversion. The deterministic inversion influences the further improvement of the resolution ratio due to the band-limited property of the seismic data; compared with deterministic inversion, stochastic inversion is based on geological statistics, and well logging data are basic condition data, so that the stochastic inversion has certain advantages in the aspects of improving resolution and reducing uncertainty.

Geostatistical is an important principle of stochastic inversion. Two-point geostatistics represent the spatial correlation of regional variables through a variation function, and are widely applied to a plurality of fields of scientific research. In the field of geophysical, sequential Gaussian simulation and sequential indication simulation are two most commonly used methods in two-point geostatistics, and are mainly used for processing problems such as modeling of continuous variables (such as porosity, shale content and the like) and discrete variables (such as lithofacies). The two-point geostatistics method has high calculation efficiency and can simulate and invert discrete variables and continuous variables, but the two-point geostatistics method utilizes a variation function to describe the correlation and the variation of the geological variables, describes the change characteristic of space only through the variation relation of the geological variables between two points in a certain direction and has poor reproducibility to the form of a target body. And the multi-point geostatistical simulation is generated as the supplement and development of two-point geostatistical. The multipoint geostatistics method takes the pixel as a basic simulation unit, combines the advantages of a target-based algorithm, utilizes a training image to describe the spatial correlation of regional variables instead of a variation function, and can better depict the geometric form of a target body. At present, due to the fact that training images of continuous variables are difficult to obtain, multi-point geostatistics are mainly used for obtaining space models of discrete variables (such as lithofacies). The general approach of multipoint geostatistical stochastic inversion is to generate a large number of equally likely reservoir lithofacies models in a multipoint geostatistical simulation stage. Sampling from the cumulative probability distribution of the corresponding reservoir physical property parameters according to the lithofacies models to obtain a plurality of pseudo-physical parameter distributions; then converting the elastic information into elastic information by utilizing a rock physics theory; and then, the forward synthetic seismic record is matched with the actual seismic record, if the requirements cannot be met, simulation is carried out again, generally hundreds of times of simulation are needed to obtain a satisfactory inversion result, the calculated amount is very large, meanwhile, a large amount of forward simulation is needed, the calculation efficiency is very low due to the huge calculated amount, and the application of the multi-point geostatistics in reservoir characterization and modeling is hindered.

It can be seen that the current random inversion method based on multi-point geostatistics has the problems of limited application range, low calculation efficiency, capability of obtaining only lithofacies inversion results and the like.

Disclosure of Invention

The embodiment of the invention provides a geological random inversion method and device based on multipoint geostatistics, and aims to solve the problems that the application range is limited, the calculation efficiency is low and only a lithofacies inversion result can be obtained in the conventional random inversion method based on the multipoint geostatistics.

In order to achieve the purpose, the invention adopts the following technical scheme:

a geological stochastic inversion method based on multipoint geostatistics comprises the following steps:

step 1, determining an initial lithofacies model through a multipoint geostatistical random simulation algorithm, and acquiring preset maximum simulation times;

step 2, obtaining preset maximum disturbance times and weight factors, and performing probability disturbance updating on the initial lithofacies model according to a probability disturbance method to generate a lithofacies model after the probability disturbance updating;

step 3, counting distribution functions of reservoir physical property parameters in different rock phases, and converting the rock phase model after probability disturbance updating into a pseudo reservoir physical property parameter model according to a random sampling method;

step 4, converting the pseudo reservoir physical property parameter model into an elastic parameter model according to the preset relationship between the reservoir physical property parameters and the elastic parameters of the corresponding research area in the statistical rock physical model;

step 5, performing time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model, determining a reflection coefficient, and forward modeling and synthesizing the seismic record by adopting a convolution model according to the reflection coefficient;

step 6, comparing the seismic record with a predetermined actual seismic record, and determining a lithofacies model result and a pseudo reservoir physical property parameter model result according to the maximum disturbance times and the maximum simulation times;

and 7, performing quantum annealing optimization on the pseudo reservoir physical property parameter model result to generate a reservoir physical property parameter inversion result.

Specifically, the determining an initial lithofacies model by a multipoint geostatistical stochastic simulation algorithm and obtaining a preset maximum simulation number includes:

step 101, establishing a simulation grid of the geology of a research area, obtaining a training image, distributing the condition data of the geology of the research area to grid nodes of the training image, and setting a random path;

step 102, scanning the training image according to a preset data template to obtain a data event; the data template is a spatial combination relation among grid nodes;

step 103, retrieving in the training image according to the data event to obtain a conditional probability distribution of a point to be estimated; the conditional probability distribution of the point to be estimated is:

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)

wherein, c (d)_n) A number of repetitions in a valid training image of the training images for the data event; c. C_k(d_n) F (u) ═ F in the number of all repetitions of the data event_kThe number of (2); p { F (u) ═ F_k|d_nThe conditional probability distribution of the point u to be estimated is obtained;

step 104, sampling is carried out according to the conditional probability distribution of the point to be estimated, a simulated lithofacies of a simulated node is obtained and used as a primary simulation result, and the simulation result is added into the conditional data;

105, randomly selecting another simulation node in the points to be estimated, and repeating the steps 103 to 105 until all the points to be estimated are simulated, so as to obtain an initial lithofacies model of the research area;

and 106, acquiring the preset maximum simulation times.

Specifically, the obtaining of the preset maximum disturbance times and the preset weight factor, and performing probability disturbance updating on the initial lithofacies model according to a probability disturbance method to generate a lithofacies model after the probability disturbance updating includes:

step 201, acquiring preset maximum disturbance times and weight factors;

step 202, carrying out probability disturbance updating on the initial lithofacies model, and determining the joint probability after the probability disturbance updating;

the updated joint probability of the probability perturbation is as follows: p (F)_k|d_obs)＝(1-r)*i₀(u,f_k)+r*P(F_k) (ii) a Wherein the content of the first and second substances,

is an indicator variable; p (F)_k) Representing the prior probability of lithofacies distribution at an unknown point u in space; p (F)_k|d_obs) Perturb the updated joint probability for the probability of the unknown point u; r represents a weight factor, wherein the weight factor is a weight factor adopting nonlinear transformation;

step 203, updating the joint probability P (F) after the probability perturbation_k|d_obs) And a previously obtained multipoint probability P (F)_k| G) fusion, determining the conditional probability distribution after the disturbance fusion;

the conditional probability distribution after disturbance fusion is

Wherein tau is₁And τ₂As a model parameter, a ═ 1-P (F)_k))/P(F_k)，b＝(1-P(F_k|G))/P(F_k|G)，c＝(1-P(F_k|d_obs))/P(F_k|d_obs)；

And 204, converting the conditional probability distribution after the disturbance fusion into cumulative probability distribution, and sampling the simulation nodes in sequence by adopting a random sampling method according to the cumulative probability distribution to generate a lithofacies model after the probability disturbance updating.

Specifically, the counting of the distribution functions of the reservoir physical parameters in different rock phases and the conversion of the probability disturbance updated rock phase model into the pseudo reservoir physical parameter model according to a random sampling method includes:

301, according to the pre-obtained logging data, counting the distribution characteristics of the reservoir physical property parameters in different rock phases to obtain the distribution functions of the reservoir physical property parameters in different rock phases;

and 302, randomly sampling each grid node in the lithofacies model after the probability disturbance updating according to a random sampling method to obtain a physical parameter state value at each corresponding grid node, and forming the pseudo reservoir physical parameter model by the physical parameter state value at each corresponding grid node.

Specifically, the converting the pseudo reservoir physical property parameter model into an elastic parameter model according to a preset relationship between reservoir physical property parameters and elastic parameters of a corresponding research area in a statistical petrophysical model includes:

according to the preset relationship between the reservoir physical property parameter and the elastic parameter of the corresponding research area in the statistical rock physical model: e ═ f_RPM(r) + epsilon, converting the pseudo reservoir physical property parameter model into an elastic parameter model;

wherein f is_RPM(r) is a deterministic petrophysical model obtained by theoretical petrophysical models or empirical petrophysical relationships; e is an elasticity parameter; r is a reservoir physical property parameter; ε is the error term of the deterministic petrophysical model.

Specifically, the performing time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model, determining a reflection coefficient, and forward modeling and synthesizing the seismic record by using a convolution model according to the reflection coefficient includes:

step 501, performing time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model;

502, calculating a reflection coefficient according to the time domain elastic parameter model;

step 503, forward modeling and synthesizing the seismic record by using the convolution model

Wherein, w is the seismic record synthesized by forward modeling by adopting a convolution model; b is seismic wavelets extracted from the actual seismic record for synthesizing the seismic record; r is the reflection coefficient.

Specifically, the comparing the seismic record with a predetermined actual seismic record, and determining a lithofacies model result and a pseudo reservoir physical property parameter model result according to the maximum disturbance times and the maximum simulation times includes:

step 601, comparing the seismic record with a predetermined actual seismic record, and determining the error E | | | d between the seismic record and the predetermined actual seismic record_obs-d_syn||²(ii) a Wherein d is_obsIs the actual seismic record; d_synA synthetic seismic record for the forward actor;

step 602, judging whether the error E is less than or equal to a preset error tolerance;

603, when the error E is less than or equal to the error tolerance, taking the lithofacies model after the probability disturbance updating as the lithofacies model result, and taking the pseudo reservoir physical property parameter model as the pseudo reservoir physical property parameter model result;

step 604, when the error E is greater than the error tolerance, judging whether the frequency of probability disturbance updating reaches the maximum disturbance frequency;

step 605, when the frequency of probability disturbance updating does not reach the maximum disturbance frequency, returning to the process of probability disturbance updating in the step 2 until the error E in the step 602 is less than or equal to the preset error tolerance;

step 606, when the number of times of probability disturbance updating reaches the maximum disturbance number, judging whether the number of times of current simulation reaches the maximum simulation number;

step 607, when the current simulation times do not reach the maximum simulation times, modifying the random path, and returning to the process of determining the initial lithofacies model through the multipoint geostatistical random simulation algorithm in the step 1 until the error E in the step 602 is less than or equal to the preset error tolerance;

and 608, when the current simulation times reach the maximum simulation times, acquiring a lithofacies model updated by probability disturbance corresponding to the error E closest to the preset error tolerance as the lithofacies model result, and acquiring a pseudo reservoir physical property parameter model corresponding to the error E closest to the preset error tolerance as the pseudo reservoir physical property parameter model result.

Specifically, the performing quantum annealing optimization on the pseudo reservoir physical property parameter model result to generate a reservoir physical property parameter inversion result includes:

701, acquiring a physical property parameter model result of a pseudo reservoir to be processed, and setting an initial system kinetic energy C gamma (t) and an annealing termination condition;

step 702, according to the formula

Computing the objective function E (m)^(k)) (ii) a Wherein d is_jIs observed seismic data; m is^(k)Model parameters for the kth iteration;

step 703, according to the formula m^(k+1)＝m^(k)+ xi alpha updates the model parameter to determine the model parameter m of the (k + 1) th iteration^(k+1)(ii) a Xi is a random number between 0 and 1, and alpha is a step length;

step 704, according to the formula Δ E ═ E (m)^(k+1))-E(m^(k)) Calculating system potential energy delta E;

step 705, at Δ E<At 0, m is adopted^(k+1)Replacement of m^(k)；

Step 706, when Δ E is larger than or equal to 0, according to the probability acceptance criterion ρ (H) ═ E when the system temperature is T^-H/PT)^PAt m^(k ⁺¹⁾And m^(k)Determining the model parameter corresponding to the larger value of the probability rho (H); wherein H ═ Δ E + C Γ (t); p is the total number of model parameters；

Step 707, when the model parameter of the current iteration meets the annealing end condition, dividing m^(k+1)As a result of inversion of reservoir physical property parameters;

step 708, when the current iterative model parameter does not satisfy the annealing termination condition, returning to the model parameter updating process of step 703 until the iterative model parameter satisfies the annealing termination condition.

A geosynchronous inversion apparatus based on multi-point geostatistics, comprising:

the initial lithofacies model determining unit is used for determining an initial lithofacies model through a multi-point geostatistical random simulation algorithm and acquiring preset maximum simulation times;

the probability disturbance updating unit is used for acquiring preset maximum disturbance times and weight factors, and performing probability disturbance updating on the initial lithofacies model according to a probability disturbance method to generate a lithofacies model after the probability disturbance updating;

the pseudo reservoir physical property parameter model conversion unit is used for counting the distribution functions of reservoir physical property parameters in different rock phases and converting the rock phase model after probability disturbance updating into a pseudo reservoir physical property parameter model according to a random sampling method;

the elastic parameter model conversion unit is used for converting the pseudo reservoir physical property parameter model into an elastic parameter model according to the preset relationship between the reservoir physical property parameters and the elastic parameters of the corresponding research area in the statistical rock physical model;

the forward modeling synthesis unit is used for carrying out time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model, determining a reflection coefficient and forward modeling and synthesizing the seismic record by adopting a convolution model according to the reflection coefficient;

the result determining unit is used for comparing the seismic record with a predetermined actual seismic record and determining a lithofacies model result and a pseudo reservoir physical property parameter model result according to the maximum disturbance times and the maximum simulation times;

and the annealing optimization unit is used for carrying out quantum annealing optimization on the pseudo reservoir physical property parameter model result to generate a reservoir physical property parameter inversion result.

Furthermore, the initial lithofacies model determination unit is specifically configured to perform:

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)

and 106, acquiring the preset maximum simulation times.

Furthermore, the probability perturbation updating unit is specifically configured to perform:

step 201, acquiring preset maximum disturbance times and weight factors;

the conditional probability distribution after disturbance fusion is

In addition, the pseudo reservoir physical property parameter model conversion unit is specifically configured to perform:

Furthermore, the elastic parametric model conversion unit is specifically configured to:

Furthermore, the forward synthesis unit is specifically configured to perform:

Furthermore, the result determination unit is specifically configured to perform:

605, returning to the process of probability disturbance updating of the probability disturbance updating unit when the frequency of probability disturbance updating does not reach the maximum disturbance frequency until the error E in the step 602 is less than or equal to the preset error tolerance;

step 607, when the current simulation times do not reach the maximum simulation times, modifying the random path, and returning to the process of determining the initial lithofacies model by the multi-point geostatistical random simulation algorithm of the initial lithofacies model determining unit until the error E in the step 602 is less than or equal to the preset error tolerance;

Furthermore, the annealing optimization unit is specifically configured to perform:

step 702, according to the formula

step 705, at Δ E<At 0, m is adopted^(k+1)Replacement of m^(k)；

Step 706, when Δ E is larger than or equal to 0, according to the probability acceptance criterion ρ (H) ═ E when the system temperature is T^-H/PT)^PAt m^(k ⁺¹⁾And m^(k)Determining the model parameter corresponding to the larger value of the probability rho (H); wherein H ═ Δ E + C Γ (t); p is the total number of model parameters;

A computer-readable storage medium, on which a computer program is stored which, when executed by a processor, carries out the steps of:

A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the steps when executing the program of:

According to the geological random inversion method and device based on the multi-point geostatistics, the probability disturbance inversion strategy is introduced into the multi-point geological random inversion, and a geological random inversion method based on the multi-point geostatistics is formed by utilizing the statistical rock physics and earthquake forward theory and combining the quantum annealing optimization algorithm; the invention can utilize multi-point geostatistics to randomly simulate lithofacies, considers the correlation among multiple points in the space and can better depict the facies state of the target body; by utilizing the probability disturbance inversion strategy, the simulation times of multipoint geological statistics can be reduced, the calculated amount is reduced, random inversion results are obtained in fewer iteration times, and the practicability is high. In addition, the generated reservoir physical property parameter model with poor space continuity is quickly optimized through a quantum annealing optimization algorithm, so that a final reservoir physical property parameter inversion result is obtained. The invention can simultaneously obtain the inversion result of the lithofacies and the physical property parameters with higher resolution and precision in fewer iteration times. The random inversion result can provide relatively reliable reference information for reservoir characterization and reservoir description. Therefore, the method can solve the problems that the existing random inversion method based on the multi-point geostatistics is limited in application range and low in calculation efficiency, and only a lithofacies inversion result can be obtained.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to these drawings without creative efforts.

Fig. 1 is a first flowchart of a geological stochastic inversion method based on multi-point geostatistics according to an embodiment of the present invention;

FIG. 2 is a flowchart illustrating an embodiment of step 1;

FIG. 3 is a diagram illustrating an event of scanning a training image to obtain data according to an embodiment of the present invention;

FIG. 4(a) is a schematic diagram of a lithofacies section of test model data in an embodiment of the present invention;

FIG. 4(b) is a schematic porosity profile of test model data in an embodiment of the present invention;

FIG. 5 is a flowchart illustrating an embodiment of step 2;

FIG. 6(a) is a schematic diagram of a lithofacies model obtained according to a probability perturbation method;

FIG. 6(b) is a schematic diagram of a porosity inversion result of a lithofacies model obtained according to a probability perturbation method;

FIG. 6(c) is a schematic diagram of a porosity inversion result error of a lithofacies model obtained according to a probability perturbation method;

FIG. 7 is a flowchart illustrating an embodiment of step 3;

FIG. 8 is a schematic of an A-well data for multi-point geostatistical stochastic inversion in an embodiment of the invention;

FIG. 9 is a schematic of a B-well data for multi-point geostatistical stochastic inversion in an embodiment of the invention;

FIG. 10 is a flowchart illustrating an embodiment of step 5;

FIG. 11 is a schematic illustration of an actual seismic recording in an embodiment of the invention;

FIG. 12(a) is a schematic diagram of the results of a lithofacies inversion of an actual seismic record in an embodiment of the present invention;

FIG. 12(b) is a schematic diagram of the results of porosity inversion for an actual seismic record in an embodiment of the invention;

FIG. 13 is a flowchart illustrating an embodiment of step 6;

FIG. 14 is a flowchart illustrating an embodiment of step 7;

fig. 15 is a schematic structural diagram of a geosynchronous inversion apparatus based on multi-point geostatistical in an embodiment of the present invention.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

As shown in fig. 1, an embodiment of the present invention provides a method for geologic stochastic inversion based on multi-point geostatistics, including:

step 1, determining an initial lithofacies model through a multipoint geostatistical random simulation algorithm, and acquiring preset maximum simulation times.

And 2, acquiring preset maximum disturbance times and weight factors, and performing probability disturbance updating on the initial lithofacies model according to a probability disturbance method to generate a lithofacies model after the probability disturbance updating.

And 3, counting distribution functions of the reservoir physical property parameters in different rock phases, and converting the rock phase model after the probability disturbance updating into a pseudo reservoir physical property parameter model according to a random sampling method.

And 4, converting the pseudo reservoir physical property parameter model into an elastic parameter model according to the preset relationship between the reservoir physical property parameters and the elastic parameters of the corresponding research area in the statistical rock physical model.

And 5, performing time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model, determining a reflection coefficient, and forward modeling and synthesizing the seismic record by adopting a convolution model according to the reflection coefficient.

And 6, comparing the seismic record with a predetermined actual seismic record, and determining a lithofacies model result and a pseudo reservoir physical property parameter model result according to the maximum disturbance times and the maximum simulation times.

The embodiment of the invention introduces a probability disturbance inversion strategy into multipoint geostatistical stochastic inversion, and forms a geostatistical inversion method based on multipoint geostatistical by using statistical petrophysics and seismic forward theory and combining a quantum annealing optimization algorithm; the invention can utilize multi-point geostatistics to randomly simulate lithofacies, considers the correlation among multiple points in the space and can better depict the facies state of the target body; by utilizing the probability disturbance inversion strategy, the simulation times of multipoint geological statistics can be reduced, the calculated amount is reduced, random inversion results are obtained in fewer iteration times, and the practicability is high. In addition, the generated reservoir physical property parameter model with poor space continuity is quickly optimized through a quantum annealing optimization algorithm, so that a final reservoir physical property parameter inversion result is obtained. The invention can simultaneously obtain the inversion result of the lithofacies and the physical property parameters with higher resolution and precision in fewer iteration times. The random inversion result can provide relatively reliable reference information for reservoir characterization and reservoir description. Therefore, the method can solve the problems that the existing random inversion method based on the multi-point geostatistics is limited in application range and low in calculation efficiency, and only a lithofacies inversion result can be obtained.

In order to make the present invention better understood by those skilled in the art, the following steps are described in detail with reference to the above step 1 to step 7 in fig. 1:

firstly, the method is established on the basis of multi-point geostatistics simulation, wherein the multi-point geostatistics takes a pixel as a basic simulation unit and utilizes hard data as conditional data, so that a simulation result is easy to realize well logging data; meanwhile, the method combines the advantages of the target-based algorithm and can better reproduce the structural form of the simulated target body. At present, multi-point geostatistics are mainly used for researching discrete variables. The SNESIM algorithm is a multipoint geostatistical random simulation method with strong applicability, is high in calculation speed and efficiency, can effectively fuse multi-source information and perform multi-grid simulation, and is popular with the earth scientists. The key to the SNESIM algorithm is to obtain a conditional probability distribution. The tool for representing spatial correlation in multi-point geostatistics is a training image, the training image is a set of geologic patterns, the training image is a digitized image which can express the structure, the geometric form and the distribution pattern of an actual reservoir, the digitized image is generally regarded as a priori geologic concept, the prior geologic concept mainly comprises a main characteristic pattern of an area to be simulated, and the prior geologic concept does not need to be consistent with accurate logging information or obey certain distribution. The training image can reflect general characteristics of geologic body space distribution, is a quantitative representation of characteristic patterns, and can embody the spatial relationship of different characteristic patterns.

Therefore, as shown in fig. 2, the initial lithofacies model is determined by the multi-point geostatistical stochastic simulation algorithm in step 1, and the preset maximum simulation times is obtained, which may be specifically implemented by the following method:

step 101, establishing a simulation grid of the geology of the research area, obtaining a training image, distributing the condition data of the geology of the research area to grid nodes of the training image, and setting a random path.

And step 102, scanning the training image according to a preset data template to obtain a data event.

Wherein the data template is a spatial combination relationship between grid nodes.

Scanning the training image to obtain data events may be as shown in fig. 3, where u is the point to be estimated in the training image. As shown in fig. 4 a and 4B, the facies (shown in fig. 4 a) and porosity (shown in fig. 4B) sections of the test model data are obtained, the condition data are extracted therefrom, the distribution positions thereof are shown by straight lines corresponding to a and B in the figure, Depth in fig. 4 a and 4B, CDP is Common Depth Point, F is facies, and Φ is porosity.

And 103, retrieving in the training image according to the data event to obtain the conditional probability distribution of the point to be estimated.

The conditional probability distribution of the point to be estimated is as follows:

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)

wherein, c (d)_n) A number of repetitions in a valid training image of the training images for the data event; c. C_k(d_n) F (u) ═ F in the number of all repetitions of the data event_kThe number of (2); p { F (u) ═ F_k|d_nThe conditional probability distribution of the point u to be estimated.

A conditional probability distribution P { F (u) ═ F of the point to be estimated_k|d_n}≈c_k(d_n)/c(d_n) Is from the following process:

in the multi-point geostatistical stochastic simulation, the conditional probability distribution function at the point to be estimated can be expressed as:

wherein the denominator represents a data event

The probability of occurrence of α ═ 1, …, n, corresponding to the number of repetitions of the data event in the valid training image, c (d)_n) And the size N of the effective training image (e.g., the gray area in FIG. 3)_nRatio c (d) of_n)/N_n，u_αRepresenting other points in the data template; the numerator of the formula represents a data event and a point u to be estimated is F_kProbability of simultaneous occurrence of states, i.e. F (u) ═ F of all repetitions of data events_kNumber of (c)_k(d_n) And effective training image size N_nCan be expressed as c_k(d_n)/N_nSo that the conditional probability distribution of the point to be estimated is:

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)。

and 104, sampling according to the conditional probability distribution of the point to be estimated to obtain a simulated lithofacies of a simulated node as a primary simulation result, and adding the simulation result into the conditional data.

And 105, randomly selecting another simulation node in the points to be estimated, and repeating the steps 103 to 105 until all the points to be estimated are simulated, so as to obtain an initial lithofacies model of the research area.

In addition, in the embodiment of the invention, other initial lithofacies models can be obtained by changing the random path.

And 106, acquiring the preset maximum simulation times.

In addition, probability perturbation is a stochastic inversion strategy in an iterative manner. By adopting the Tau model, constraint information from different sources can be fused into a joint conditional probability in a conditional probability mode, the operation speed is high, and the calculated amount is small. Usually, about dozens of times of simulation can obtain satisfactory inversion results.

Therefore, as shown in fig. 5, in step 2, the preset maximum perturbation frequency and the preset weight factor are obtained, and the probability perturbation update is performed on the initial lithofacies model according to the probability perturbation method to generate a lithofacies model after the probability perturbation update, which may be specifically implemented by the following method:

step 201, obtaining the preset maximum disturbance times and the weight factor.

Step 202, probability disturbance updating is carried out on the initial lithofacies model, and the joint probability after probability disturbance updating is determined.

The updated joint probability of the probability perturbation is as follows: p (F)_k|d_obs)＝(1-r)*i₀(u,f_k)+r*P(F_k) (ii) a Wherein the content of the first and second substances,is an indicator variable; p (F)_k) Representing the prior probability of lithofacies distribution at an unknown point u in space; p (F)_k|d_obs) Updated joint probability for probability perturbation of the unknown point u(ii) a And r represents a weight factor, wherein the weight factor is a weight factor adopting nonlinear transformation, and the probability disturbance can be better carried out on the initial lithofacies model.

Step 203, updating the joint probability P (F) after the probability perturbation_k|d_obs) And a previously obtained multipoint probability P (F)_k| G) fusion, and determining the conditional probability distribution after the disturbance fusion.

The conditional probability distribution after disturbance fusion is

It uses a Tau model, where Tau₁And τ₂As a Tau model parameter, a ═ 1-P (F)_k))/P(F_k)，b＝(1-P(F_k|G))/P(F_k|G)，c＝(1-P(F_k|d_obs))/P(F_k|d_obs)。

For example, the lithofacies model of FIG. 6(a), the porosity inversion result of FIG. 6(b), and the porosity inversion result error of FIG. 6(c) may be obtained according to a probability perturbation method. Where Δ φ is the porosity inversion result error.

Specifically, as shown in fig. 7, in step 3, the distribution functions of the reservoir physical parameters in different rock phases are counted, and the rock phase model after the probability perturbation update is converted into the pseudo reservoir physical parameter model according to a random sampling method, which can be implemented in the following manner:

step 301, according to the pre-obtained logging data, counting the distribution characteristics of the reservoir physical property parameters in different rock phases to obtain the distribution functions of the reservoir physical property parameters in different rock phases.

Since the spatial correlation of the physical property parameter change is not considered in the conversion process, a physical property parameter distribution model is obtained only by random sampling, and the model has poor spatial correlation, so that the model is called a pseudo reservoir physical property parameter model. And performing sampling on all grid nodes once to obtain a pseudo reservoir physical property parameter model, and performing random sampling for multiple times to obtain multiple pseudo reservoir physical property parameter models.

In addition, after the pseudo reservoir physical property parameter model is obtained, the pseudo reservoir physical property parameter model can be converted into elastic parameters (such as speed, density, wave impedance and the like) according to the rock physics theory, and forward modeling is facilitated. In actual seismic exploration, underground reservoir conditions are very complex, and factors such as different pore structures, roundness of mineral particles, formation temperature, pressure condition changes, shale content and the like can cause deviation of the physical relationship of the elastic parameters and the reservoir physical property parameters rock to a certain extent. Deterministic petrophysical models, do not take into account the effects of such deviations. To reduce such effects, a statistical petrophysical model is constructed by adding random error terms to the deterministic petrophysical relationships.

Therefore, in the step 4, the pseudo reservoir physical property parameter model is converted into the elastic parameter model according to the preset relationship between the reservoir physical property parameter and the elastic parameter of the corresponding research area in the statistical rock physical model, and the method can be implemented as follows:

according to the preset relationship between the reservoir physical property parameter and the elastic parameter of the corresponding research area in the statistical rock physical model: e ═ f_RPM(r) + epsilon, and converting the pseudo reservoir physical property parameter model into an elastic parameter model.

Wherein f is_RPM(r) is a deterministic petrophysical model obtained by theoretical petrophysical models or empirical petrophysical relationships; e is an elasticity parameter; r is a reservoir physical property parameter; epsilon is an error term of a deterministic petrophysical model, can be estimated through the relative difference between actual logging information and a deterministic petrophysical relation, and a truncated Gaussian distribution with a mean value of zero is often selected to weaken the underground complex structureThe influence on the relationship between the two.

For the characteristics of different research areas, a more suitable petrophysical model or statistical petrophysical model should be selected to reduce the uncertainty in the conversion process from physical parameters to elastic parameters, and therefore the petrophysical model should be carefully selected.

In addition, because the well data is acquired in the depth domain (for example, the well data A in FIG. 8 and the well data B in FIG. 9, where φ is porosity, Vsh is shale relative volume, Sw is water saturation, Vp is pore volume, ρ is density, Facies is seismic Facies), the generated lithofacies model, the pseudo-reservoir physical property parameter model and the elastic parameter model are conditional data, and are thus data in the depth domain, and the seismic records are data in the time domain, therefore, between the synthetic seismic data, the depth domain data is converted into time domain data, and after the elastic parameter model in the time domain is obtained, forward modeling can be carried out to obtain the synthetic seismic record d_synTo facilitate recording with actual earthquake d_obsAnd comparing to select the lithofacies distribution and the physical parameter distribution which can be matched with the actual seismic record most, thereby realizing the random inversion of the lithofacies and the physical parameters.

Therefore, as shown in fig. 10, in the step 5, time-depth conversion is performed on the elastic parameter model to obtain a time domain elastic parameter model, reflection coefficients are determined, and forward modeling is performed on the synthetic seismic record by using a convolution model according to the reflection coefficients, which can be implemented as follows:

and 501, performing time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model.

And 502, calculating a reflection coefficient according to the time domain elastic parameter model.

The actual seismic record may be as shown in FIG. 11, the lithofacies inversion result of the actual seismic record is as shown in FIG. 12(a), and the porosity inversion result of the actual seismic record is as shown in FIG. 12 (b).

Specifically, as shown in fig. 13, in step 6, comparing the seismic record with a predetermined actual seismic record, and determining a lithofacies model result and a pseudo reservoir property parameter model result according to the maximum disturbance times and the maximum simulation times, may be implemented as follows:

step 601, comparing the seismic record with a predetermined actual seismic record, and determining the error E ═ d of the seismic record and the predetermined actual seismic record_obs-d_syn‖²。

Wherein d is_obsIs the actual seismic record; d_synThe synthetic seismic record is forward-acting.

Step 602, determining whether the error E is less than or equal to a preset error tolerance.

And 603, when the error E is less than or equal to the error tolerance, taking the lithofacies model after the probability disturbance updating as the lithofacies model result, and taking the pseudo reservoir physical property parameter model as the pseudo reservoir physical property parameter model result.

And step 604, judging whether the frequency of probability disturbance updating reaches the maximum disturbance frequency when the error E is greater than the error tolerance.

And 605, when the frequency of the probability disturbance updating does not reach the maximum disturbance frequency, returning to the process of the probability disturbance updating in the step 2 until the error E in the step 602 is less than or equal to the preset error tolerance.

And 606, judging whether the current simulation times reach the maximum simulation times or not when the probability disturbance updating times reach the maximum disturbance times.

And 607, when the current simulation times do not reach the maximum simulation times, modifying the random path, and returning to the process of determining the initial lithofacies model through the multipoint geostatistical random simulation algorithm in the step 1 until the error E in the step 602 is less than or equal to the preset error tolerance.

The pseudo reservoir physical property parameter model output in the previous stage can meet the preset error tolerance, so that the synthesized seismic record is closer to the actual seismic record and closer to the underground actual physical property parameter distribution condition, but the pseudo reservoir physical property parameter model is obtained according to random sampling without considering spatial correlation, so that the spatial continuity is poor, which is not consistent with the actual condition, and the seismic inversion generally has multi-solution property. Therefore, the output pseudo reservoir physical property parameter model result is optimized, and the spatial continuity of the pseudo reservoir physical property parameter model result is improved. The quantum annealing has certain advantages in the aspects of annealing convergence speed and local minimum avoidance, and has important value for improving the accuracy and efficiency of random inversion.

Therefore, as shown in fig. 14, quantum annealing optimization is performed on the pseudo reservoir property parameter model result in step 7 to generate a reservoir property parameter inversion result, which can be implemented as follows:

step 701, obtaining a physical property parameter model result of a pseudo reservoir to be processed, and setting an initial system kinetic energy C gamma (t) and an annealing termination condition.

Step 702, according to the formula

Computing the objective function E (m)^(k)) (ii) a Wherein d is_jIs observed seismic data; m is^(k)Is the model parameter for the kth iteration.

Step 703, according to the formula m^(k+1)＝m^(k)+ xi α advances the model parameterLine updating and determining the model parameter m of the (k + 1) th iteration^(k+1)(ii) a ξ is a random number between 0 and 1 and α is the step size.

Step 704, according to the formula Δ E ═ E (m)^(k+1))-E(m^(k)) And calculating the system potential energy delta E.

Step 705, at Δ E<At 0, m is adopted^(k+1)Replacement of m^(k)。

Step 706, when Δ E is larger than or equal to 0, according to the probability acceptance criterion ρ (H) ═ E when the system temperature is T^-H/PT)^PAt m^(k ⁺¹⁾And m^(k)Determining the model parameter corresponding to the larger value of the probability rho (H); wherein H ═ Δ E + C Γ (t); p is the total number of model parameters; c is a constant and M is the number of seismic traces.

Step 707, when the model parameter of the current iteration meets the annealing end condition, dividing m^(k+1)As a result of inversion of reservoir physical parameters.

In summary, the geological stochastic inversion method based on multi-point geostatistics provided by the embodiment of the present invention has the following advantages: 1. applying multipoint geostatistics random simulation, considering the correlation among multiple points in the space, and better representing the characteristics of the target such as the shape and the like; 2. a probability disturbance inversion strategy is introduced, constraint information of multiple sources is fused by utilizing a Tau model, the simulation times are reduced, the inversion process is accelerated, the inversion efficiency can be improved, and the calculation time is shortened; 3. the method combines the quantum annealing algorithm, optimizes the obtained pseudo-physical parameter distribution, reduces the search space and obtains the inversion result of the physical parameters more quickly; 4. the invention can obtain the inversion result of the lithofacies and reservoir physical property parameters with higher resolution in fewer iteration times, thereby providing a more reliable reservoir model for reservoir evaluation and fine description of an oil reservoir; 5. the invention can also output the inversion result of the elastic parameter according to the requirement, which is convenient for the follow-up research; 6. on the premise of not changing inversion process, the random inversion can be easily popularized from post-stack to pre-stack by utilizing pre-stack seismic data and forward operators.

In addition, as shown in fig. 15, an embodiment of the present invention further provides a geological stochastic inversion apparatus based on multi-point geostatistics, including:

and the initial lithofacies model determining unit 11 is configured to determine an initial lithofacies model through a multi-point geostatistical random simulation algorithm, and acquire a preset maximum simulation number.

And the probability disturbance updating unit 12 is configured to obtain preset maximum disturbance times and weight factors, perform probability disturbance updating on the initial lithofacies model according to a probability disturbance method, and generate a lithofacies model after the probability disturbance updating.

And the pseudo reservoir physical property parameter model conversion unit 13 is used for counting the distribution functions of the reservoir physical property parameters in different rock phases and converting the rock phase model after the probability disturbance updating into the pseudo reservoir physical property parameter model according to a random sampling method.

And the elastic parameter model conversion unit 14 is used for converting the pseudo reservoir physical property parameter model into an elastic parameter model according to the preset relationship between the reservoir physical property parameters and the elastic parameters of the corresponding research area in the statistical rock physical model.

And the forward modeling synthesis unit 15 is configured to perform time-depth conversion on the elastic parameter model to obtain a time domain elastic parameter model, determine a reflection coefficient, and forward model synthesize a seismic record by using a convolution model according to the reflection coefficient.

And the result determining unit 16 is used for comparing the seismic record with a predetermined actual seismic record and determining a lithofacies model result and a pseudo reservoir physical property parameter model result according to the maximum disturbance times and the maximum simulation times.

And the annealing optimization unit 17 is configured to perform quantum annealing optimization on the pseudo reservoir property parameter model result to generate a reservoir property parameter inversion result.

Furthermore, the initial lithofacies model determining unit 11 is specifically configured to perform:

Step 102, scanning the training image according to a preset data template to obtain a data event; the data template is a spatial combination relation among grid nodes.

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)

And 106, acquiring the preset maximum simulation times.

Furthermore, the probability perturbation updating unit 12 is specifically configured to perform:

step 201, obtaining the preset maximum disturbance times and the weight factor.

is an indicator variable; p (F)_k) Representing the prior probability of lithofacies distribution at an unknown point u in space; p (F)_k|d_obs) Perturb the updated joint probability for the probability of the unknown point u; r represents a weight factor, which is a weight factor using a nonlinear transformation.

The conditional probability distribution after disturbance fusion is

Wherein tau is₁And τ₂As a model parameter, a ═ 1-P (F)_k))/P(F_k)，b＝(1-P(F_k|G))/P(F_k|G)，c＝(1-P(F_k|d_obs))/P(F_k|d_obs)。

In addition, the pseudo reservoir property parameter model conversion unit 13 is specifically configured to perform:

Furthermore, the elastic parameter model transformation unit 14 is specifically configured to:

Furthermore, the forward synthesis unit 15 is specifically configured to perform:

Furthermore, the result determining unit 16 is specifically configured to perform:

step 601, comparing the seismic record with a predetermined actual seismic record, and determining the error E | | | d between the seismic record and the predetermined actual seismic record_obs-d_syn||²(ii) a Wherein d is_obsIs the actual seismic record; d_synThe synthetic seismic record is forward-acting.

And 605, when the frequency of the probability disturbance updating does not reach the maximum disturbance frequency, returning to the process of the probability disturbance updating unit 12 for probability disturbance updating until the error E in the step 602 is less than or equal to the preset error tolerance.

And 607, when the current simulation times do not reach the maximum simulation times, modifying the random path, and returning to the process of determining the initial lithofacies model by the multi-point geostatistical random simulation algorithm of the initial lithofacies model determining unit 11 until the error E in the step 602 is less than or equal to the preset error tolerance.

Furthermore, the annealing optimization unit 17 is specifically configured to perform:

Step 702, according to the formulaComputing the objective function E (m)^(k)) (ii) a Wherein d is_jIs observed seismic data; m (m)^k)Is the model parameter for the kth iteration.

Step 703, according to the formula m^(k+1)＝m^(k)+ xi alpha updates the model parameter to determine the model parameter m of the (k + 1) th iteration^(k+1)(ii) a ξ is a random number between 0 and 1 and α is the step size.

Step 705, at Δ E<At 0, m is adopted^(k+1)Replacement of m^(k)。

Step 706, when Δ E is larger than or equal to 0, according to the probability acceptance criterion ρ (H) ═ E when the system temperature is T^-H/PT)^PAt m^(k ⁺¹⁾And m^(k)Determining the model parameter corresponding to the larger value of the probability rho (H); wherein H ═ Δ E + C Γ (t); p is the total number of model parameters.

It should be noted that, for a specific implementation manner of the geological stochastic inversion apparatus based on multi-point geostatistics provided in the embodiment of the present invention, reference may be made to the above method embodiment, and details are not described here again.

The embodiment of the invention introduces a probability disturbance inversion strategy into multipoint geostatistical stochastic inversion, and forms a geostatistical inversion device based on multipoint geostatistical by using statistical petrophysics and seismic forward theory and combining a quantum annealing optimization algorithm; the invention can utilize multi-point geostatistics to randomly simulate lithofacies, considers the correlation among multiple points in the space and can better depict the facies state of the target body; by utilizing the probability disturbance inversion strategy, the simulation times of multipoint geological statistics can be reduced, the calculated amount is reduced, random inversion results are obtained in fewer iteration times, and the practicability is high. In addition, the generated reservoir physical property parameter model with poor space continuity is quickly optimized through a quantum annealing optimization algorithm, so that a final reservoir physical property parameter inversion result is obtained. The invention can simultaneously obtain the inversion result of the lithofacies and the physical property parameters with higher resolution and precision in fewer iteration times. The random inversion result can provide relatively reliable reference information for reservoir characterization and reservoir description. Therefore, the method can solve the problems that the existing random inversion method based on the multi-point geostatistics is limited in application range and low in calculation efficiency, and only a lithofacies inversion result can be obtained.

In addition, an embodiment of the present invention further provides a computer-readable storage medium, on which a computer program is stored, where the computer program, when executed by a processor, implements the following steps:

In addition, an embodiment of the present invention further provides a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, where the processor implements the following steps when executing the program:

As will be appreciated by one skilled in the art, embodiments of the present invention may be provided as a method, system, or computer program product. Accordingly, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.

The present invention is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each flow and/or block of the flow diagrams and/or block diagrams, and combinations of flows and/or blocks in the flow diagrams and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.

These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.

These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.

The principle and the implementation mode of the invention are explained by applying specific embodiments in the invention, and the description of the embodiments is only used for helping to understand the method and the core idea of the invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, there may be variations in the specific embodiments and the application scope, and in summary, the content of the present specification should not be construed as a limitation to the present invention.

Claims

1. A geological stochastic inversion method based on multi-point geostatistics is characterized by comprising the following steps:

step 7, performing quantum annealing optimization on the pseudo reservoir physical property parameter model result to generate a reservoir physical property parameter inversion result;

wherein, the step 5 specifically comprises:

2. The method of claim 1, wherein the determining an initial lithofacies model by a multi-point geostatistical stochastic simulation algorithm and obtaining a preset maximum number of simulations comprises:

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)

and 106, acquiring the preset maximum simulation times.

3. The method according to claim 1, wherein the obtaining of the preset maximum perturbation times and the preset weighting factor, and the probability perturbation updating of the initial lithofacies model according to the probability perturbation method to generate the lithofacies model after the probability perturbation updating comprises:

step 201, acquiring preset maximum disturbance times and weight factors;

the updated joint probability of the probability perturbation is as follows: p (F)_k|d_obs)＝(1-r)*i₀(u,f_k)+r*P(F_k) (ii) a Wherein the content of the first and second substances,is an indicator variable; p (F)_k) Representing the prior probability of lithofacies distribution at an unknown point u in space; p (F)_k|d_obs) Perturb the updated joint probability for the probability of the unknown point u; r represents a weight factor, wherein the weight factor is a weight factor adopting nonlinear transformation;

the conditional probability distribution after disturbance fusion is

4. The method of claim 1, wherein the step of counting distribution functions of the reservoir physical parameters in different rock phases and converting the probability disturbance updated rock phase model into a pseudo reservoir physical parameter model according to a random sampling method comprises the following steps:

5. The method according to claim 1, wherein the converting the pseudo reservoir physical property parameter model into an elastic parameter model according to a preset relationship between reservoir physical property parameters and elastic parameters of a corresponding research region in the statistical petrophysical model comprises:

6. The method of claim 1, wherein comparing the seismic record to a predetermined actual seismic record and determining lithofacies model results and pseudo-reservoir property parameter model results based on the maximum number of perturbations and the maximum number of simulations comprises:

7. The method of claim 1, wherein the performing quantum annealing optimization on the pseudo reservoir property parameter model result to generate a reservoir property parameter inversion result comprises:

step 702, according to the formula

step 705, when delta E is less than 0, m is adopted^(k+1)Replacement of m^(k)；

Step 706, when Δ E is larger than or equal to 0, according to the probability acceptance criterion ρ (H) ═ E when the system temperature is T^-H/PT)^PAt m^(k+1)And m^(k)Determining the model parameter corresponding to the larger value of the probability rho (H); wherein H ═ Δ E + C Γ (t); p is the total number of model parameters;

8. A geosynchronous inversion apparatus based on multi-point geostatistics, comprising:

the annealing optimization unit is used for carrying out quantum annealing optimization on the pseudo reservoir physical property parameter model result to generate a reservoir physical property parameter inversion result;

wherein the forward synthesis unit is specifically configured to perform:

9. The apparatus according to claim 8, wherein the initial lithofacies model determination unit is specifically configured to perform:

P{F(u)＝F_k|d_n}≈c_k(d_n)/c(d_n)

and 106, acquiring the preset maximum simulation times.

10. The apparatus according to claim 8, wherein the probability perturbation updating unit is specifically configured to perform:

step 201, acquiring preset maximum disturbance times and weight factors;

the conditional probability distribution after disturbance fusion is

11. The apparatus according to claim 8, wherein the pseudo reservoir physical property parameter model transformation unit is specifically configured to perform:

12. The apparatus according to claim 8, wherein the elastic parametric model transformation unit is specifically configured to:

13. The apparatus according to claim 8, wherein the result determination unit is specifically configured to perform:

step 601, comparing the seismic record with a predetermined actual seismic record, and determining the error of the seismic record and the predetermined actual seismic recordDifference E | | | d_obs-d_syn||²(ii) a Wherein d is_obsIs the actual seismic record; d_synA synthetic seismic record for the forward actor;

14. The apparatus according to claim 8, wherein the annealing optimization unit is specifically configured to perform:

step 702, according to the formulaComputing the objective function E (m)^(k)) (ii) a Wherein d is_jIs observed seismic data; m is^(k)Model parameters for the kth iteration;

15. A computer-readable storage medium, on which a computer program is stored, which program, when executed by a processor, carries out the steps of:

wherein, the step 5 specifically comprises:

16. A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein the processor when executing the program performs the steps of:

wherein, the step 5 specifically comprises:

step 503, forward modeling and synthesizing the seismic record by using the convolution modelWherein, w is the seismic record synthesized by forward modeling by adopting a convolution model; b is seismic wavelets extracted from the actual seismic record for synthesizing the seismic record; r is the reflection coefficient.