CN111077571A - Porosity inversion method for improving resolution - Google Patents
Porosity inversion method for improving resolution Download PDFInfo
- Publication number
- CN111077571A CN111077571A CN201911288767.8A CN201911288767A CN111077571A CN 111077571 A CN111077571 A CN 111077571A CN 201911288767 A CN201911288767 A CN 201911288767A CN 111077571 A CN111077571 A CN 111077571A
- Authority
- CN
- China
- Prior art keywords
- porosity
- inversion
- logging
- reflection coefficient
- longitudinal wave
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 239000011148 porous material Substances 0.000 claims abstract description 39
- 239000011435 rock Substances 0.000 claims abstract description 34
- 239000012530 fluid Substances 0.000 claims abstract description 22
- 238000002922 simulated annealing Methods 0.000 claims abstract description 18
- 239000013598 vector Substances 0.000 claims description 24
- OIGNJSKKLXVSLS-VWUMJDOOSA-N prednisolone Chemical compound O=C1C=C[C@]2(C)[C@H]3[C@@H](O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 OIGNJSKKLXVSLS-VWUMJDOOSA-N 0.000 claims description 15
- 238000000137 annealing Methods 0.000 claims description 13
- 239000011159 matrix material Substances 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 10
- 238000012986 modification Methods 0.000 claims description 8
- 230000004048 modification Effects 0.000 claims description 8
- 229920006395 saturated elastomer Polymers 0.000 claims description 7
- 230000001419 dependent effect Effects 0.000 claims description 6
- 238000006073 displacement reaction Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 230000002238 attenuated effect Effects 0.000 claims description 3
- 230000001902 propagating effect Effects 0.000 claims description 3
- 230000008901 benefit Effects 0.000 abstract description 3
- 239000002245 particle Substances 0.000 description 8
- 239000007787 solid Substances 0.000 description 7
- 230000015572 biosynthetic process Effects 0.000 description 2
- 239000013078 crystal Substances 0.000 description 2
- 230000005283 ground state Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000000750 progressive effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 239000013049 sediment Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention provides a porosity inversion method for improving resolution, which comprises the following steps: (1) deducing a fast longitudinal wave reflection coefficient approximate expression vertically incident to a fluid-containing pore medium interface, and establishing a direct quantitative relation between the porosity and seismic records based on a rock physics theory; (2) according to the relation between the porosity and the seismic record, the logging porosity and the logging reflection coefficient are used as prior constraints, and a nonlinear inversion target function is constructed; (3) calculating porosity adjustment parameters by using the logging data information, and giving prior information such as constraint parameters, maximum inversion iteration times, initial porosity solution and the like; (4) and substituting the known parameters into a nonlinear inversion target function, and realizing direct quantitative inversion of porosity by using a simulated annealing algorithm. The method has the advantages that the limitation of the traditional theory depending on single-phase medium and equivalent medium is broken through, and the direct quantitative relation between the porosity and the seismic record is fundamentally established, so that the resolution of the porosity inversion is improved.
Description
Technical Field
The invention belongs to the technical field of rock porosity seismic inversion, and particularly relates to a porosity inversion method for improving resolution.
Background
Porosity is one of the important factors controlling the storage and transport of pore fluids in rock or sediment, and is defined as the ratio of the pore volume in rock to the total volume of rock. The actual storage capacity of rock pore space for hydrocarbons in reservoir predictions is defined as the effective porosity. The porosity prediction plays an important role in lithology identification, fluid identification, reserve estimation and the like of the oil-gas-containing reservoir.
Most of the existing porosity prediction methods are based on rock physics theory, and take an equivalent medium thought as a bridge, so that a quantitative relation between the porosity and actual seismic data is established, and indirect quantitative inversion and prediction of the porosity are realized. However, the method has limited prediction capability on the porosity, and the porosity inversion result has poor lateral resolution, especially for complex oil-bearing high-pore sandstone reservoirs or carbonate reservoirs. Therefore, there is a need in the field of the present invention to develop a method for establishing a direct quantitative relationship between porosity and actual seismic data, implementing direct quantitative inversion of porosity, improving porosity inversion resolution, and providing theoretical and methodological support for accurate reservoir prediction and fluid identification.
Disclosure of Invention
In order to solve the technical problem of poor porosity inversion resolution in the prior art, the invention provides a porosity inversion method for improving resolution, and the aim of realizing the method is to provide a high-resolution porosity inversion method based on the seismic plane wave approximate analysis theory of a pore medium containing saturated fluid. And (3) assuming that seismic data to be inverted, a logging reflection coefficient and logging porosity exist, establishing a direct quantitative relation between the porosity and actual seismic data, and improving the porosity inversion resolution.
In order to realize the purpose, the invention adopts the technical proposal that the flow for realizing high-precision porosity inversion by utilizing a simulated annealing algorithm is as follows on the assumption that the seismic data to be inverted, the logging reflection coefficient and the logging porosity exist,
(1) deducing a fast longitudinal wave reflection coefficient approximate expression vertically incident to a fluid-containing pore medium interface, and establishing a direct quantitative relation between the porosity and seismic records based on a rock physics theory;
(2) according to the relation between the porosity and the seismic record, the logging porosity and the logging reflection coefficient are used as prior constraints, and a nonlinear inversion target function is constructed;
(3) calculating porosity adjustment parameters by using the logging data information, and giving prior information such as constraint parameters, maximum inversion iteration times, initial porosity solution and the like;
(4) and (3) substituting the known parameters into a nonlinear inversion target function, and combining a rapid simulated annealing algorithm to realize high-precision direct quantitative inversion of porosity.
Further, 1, the quantitative relationship between porosity and seismic data
Based on the Biot elastic pore medium theory, the reflection coefficient of the fast longitudinal wave vertically incident on the interface of the pore medium containing saturated fluid is approximate,
wherein r is the reflection coefficient of fast longitudinal wave, Kl,μ l1,2 denote the bulk and shear modulus of the rock skeleton, αl,M l1,2 is the Biot parameter; vpljWhere 1,2, j 1,2 is the longitudinal wave velocity, γ lj1,2, j 1,2 is the ratio of the displacement amplitude of the pore fluid relative to the rock skeleton and the displacement amplitude of the rock skeleton; the subscript l is 1,2 represents the upper and lower pore media, respectively, and j is 1,2 represents the fast and slow longitudinal waves, respectively.
There is a link between the modulus of the rock skeleton, the modulus of the rock matrix and the porosity, and the invention uses the formulas (8) and (9) to describe the quantitative relationship between them.
Kd＝(1.0-f)Ks/(1.0+ck·f), (8)
μd＝(1.0-f)μs/(1.0+cμ·f)， (9)
Wherein f is porosity, Kd、μdRespectively the bulk modulus and shear modulus of the rock skeleton, Ks、μsThe bulk modulus and shear modulus of the rock matrix, respectively. c. Ck、cμParameters are adjusted for the pore media.
The fast longitudinal wave and the slow longitudinal wave exist in the pore medium containing the saturated fluid, and the slow longitudinal wave is quickly attenuated and disappears when propagating in the pore medium, so that the seismic record received by the ground geophone is actually the reflection of the fast longitudinal wave. The porosity f is a function f (t) of time t, and the fast longitudinal reflection coefficient r is a function r [ f (t) ] of the porosity f (t) as shown in equations (1) to (9). According to the convolution formula, the fast compressional synthesis seismic record d [ f (t) ] directly related to porosity can be expressed as,
d[f(t)]＝w(t)*r[f(t)]+n(t), (10)
the above equation is discretized and written into the form of a vector,
dN×1＝GN×N·rN×1+nN×1. (11)
wherein w (t) is a time domain wavelet, N (t) is noise, subscript Nx1 represents an N-dimensional column vector, NxN represents an N-dimensional matrix, and N is the number of sampling points of a seismic data; the symbol' denotes a convolution operator.
Further, 2, non-linear inversion of the objective function
According to the forward modeling process in the formula (11), the priori information in the well logging is taken as constraint, the fluid-containing pore medium porosity inversion objective function F is,
wherein f is the predicted porosity, R is the predicted reflection coefficient, and d is the synthetic seismic record; f. of0And R0Porosity and reflection coefficient, respectively, from well log datadataWhich is the actual seismic data, β is the reflection coefficient constraint parameter, gamma is the porosity constraint parameter, all vectors in the objective function are nx1 dimensional column vectors,
representing the L2 norm of the vector.
In the above formula, the first term is the fitting degree of observed seismic data and synthetic seismic records, the second term represents the similarity between the reflection coefficient calculated by the predicted porosity and the logging reflection coefficient, and the third term represents the similarity between the predicted porosity and the logging porosity.
Further, 3, simulated annealing algorithm porosity inversion
The simulated annealing algorithm is based on the solid annealing principle and simulates the annealing process of crystals in physics. The solids are first warmed to a sufficient level and then allowed to cool slowly. When the temperature is increased, the internal particles of the solid become an unordered state along with the temperature rise, and the internal energy is increased; when the particles are cooled, the particles gradually get ordered, and the temperature of each particle reaches an equilibrium state; finally, the ground state is reached at normal temperature, and the internal energy is reduced to the minimum. The algorithm has progressive convergence and is a global optimization algorithm which converges with the probability 1.
(1) Metropolis convergence criterion
The inversion algorithm is obtained by continuously iteratively updating the porosityAnd obtaining an optimal solution. At the i-th iteration, the porosity f is first calculatediCorresponding objective function FiAnd then the porosity is updated to be fi+1And calculating the corresponding objective function value Fi+1. The difference between the two function values is as follows,
ΔF＝Fi+1-Fi,i≥1, (13)
if the delta F is less than 0, the porosity updating direction enables the objective function value to be reduced, and then the porosity is modified; if the delta F is larger than or equal to 0, further judging whether the delta F meets the formula (14), if so, still accepting the modification of the porosity, and if not, rejecting the modification.
0＜p＜K, (14)
Wherein, p is probability density, K is a random number, and the value range is 0 < K < 1; k is a radical ofbThe method is a Boltzmann constant, and 1, exp {. cndot } is taken as an exponential function in practical application. T is the current annealing temperature, often expressed exponentially,
Ti＝T0·(0.95)i-1,i≥1, (16)
wherein, T0Is the initial temperature, TiThe ith annealing temperature.
(2) Iteration and update of solutions
The solution of the simulated annealing algorithm is independent of the initial state, but research experience shows that: if a proper initial solution is given, the iterative solving times of the algorithm can be reduced, and the operation speed of the algorithm can be improved. Herein porosity f will be logged0Smoothed to give an initial porosity f1. Given an initial solution, an update vector Δ f may be generated using a temperature-dependent Cauchy-like distributioniI.e. by
fi+1＝fi+ξ·Δfi,i≥1, (17)
Δfi＝T·sign(q-0.5)·[(1+1/T)|2·q-1|-1]. (18)
In the above equation, ξ is the update coefficient, fiFor the current prediction of porosity, fi+1To be moreNew porosity of value fi+1∈[0,0.3](ii) a q is obedience [0,1]And N-dimensional random vectors of the distribution, wherein sign is a sign function.
The advantage of using a temperature dependent Cauchy distribution to generate the update vector is: the porosity is searched in a large range under the high temperature condition, and the searching is performed near the current predicted value as far as possible under the low temperature condition, and a flat tail is distributed like Cauchy, so that the local minimum value is easy to jump out.
(3) Nonlinear porosity inversion process
Assuming the existence of seismic data d to be inverteddataWell logging reflection coefficient R0Logging porosity f0Initial porosity f1The maximum iteration number L and the logging constraint parameters β and gamma, the process for realizing high-precision porosity inversion based on the simulated annealing algorithm is as follows,
1) porosity f for the ith iterationiCalculating a new porosity value f using equation (17)i+1Then, based on the formula (12), calculating the inverse objective function value F before and after the porosity updatingiAnd Fi+1；
2) Calculating the difference delta F between the two objective function values by combining the formula (13), judging whether the condition delta F is less than 0, if so, accepting the updating of the porosity value, otherwise, continuously judging whether the formula (14) is met, if so, still accepting the updating of the porosity value, otherwise, refusing the updating;
3) judging whether the maximum iteration frequency L is reached, if so, exiting the loop to obtain an optimal porosity value; otherwise, iteratively updating the current annealing temperature by using the formula (16), and returning to 1) to continuously and iteratively update the porosity value f.
By adopting the technical scheme, the invention has the following beneficial effects: the invention realizes high-precision porosity inversion by taking the fluid-containing pore medium vertical incidence seismic plane wave approximate analytic expression as a theoretical basis and combining a rapid simulated annealing algorithm. Compared with the existing porosity inversion technology, the technology breaks through the traditional limitation depending on single-phase medium and equivalent medium theory, and fundamentally establishes the direct quantitative relation between the porosity and the seismic record, which is the theoretical advancement of the invention. By skillfully utilizing a simulated annealing algorithm, the prior information of logging data, rock physics, geology and the like is introduced into the inversion target function, so that the inversion is prevented from being trapped in a local minimum value, the inversion multi-solution is reduced, and the high-precision porosity inversion is finally realized.
Drawings
FIG. 1 is a flow chart of a porosity inversion method of the present invention;
FIG. 2 is a single-channel sandstone pore medium model, wherein (a), (b), (c) and (d) respectively show the porosity, reflection coefficient, synthetic seismic record and low-frequency porosity corresponding to the model;
FIG. 3 is an energy ratio of porosity for different initial solution cases;
FIG. 4 shows the results of porosity inversion for different initial solutions, where (a), (b), (c) and (d) are respectively the 100 th, 200 th, 300 th and 400 th iterative inversions;
FIG. 5 shows the reflection coefficient inversion results for different initial solutions, where (a), (b), (c), and (d) are respectively the 100 th, 200 th, 300 th, and 400 th iterative inversions;
FIG. 6 is a synthetic seismic record for different initial solution cases, and (a), (b), (c) and (d) are seismic records reconstructed from 100 th, 200 th, 300 th and 400 th iterative inversion results respectively;
FIG. 7 shows a parameter β of 0.3, where γ is the corresponding energy ratio of porosity at 0.005 (solid line), 0.01 (dashed line), 0.1 (dashed line), and 0.5, respectively;
FIG. 8 is a graph comparing the porosity inversion results (dashed line) with the original model (solid line) when the parameter β is 0.3 and γ is 0.005(a), 0.01(b), 0.1(c), and 0.5(d), respectively;
fig. 9 actual seismic data and porosity inversion results, (a) actual seismic data, (b) a low frequency porosity model, (c) porosity inverted by a conventional method, and (d) porosity inverted by the method of the present invention.
Detailed Description
The present invention is described in further detail below with reference to specific examples.
Example (b): referring to fig. 1, which is a flow chart of the method of the present invention, assuming that seismic data to be inverted, a logging reflection coefficient and a logging porosity exist, a flow for realizing high-precision porosity inversion by using a simulated annealing algorithm is as follows,
(1) deducing a fast longitudinal wave reflection coefficient approximate expression vertically incident to a fluid-containing pore medium interface, and establishing a direct quantitative relation between the porosity and seismic records based on a rock physics theory;
(2) according to the relation between the porosity and the seismic record, the logging porosity and the logging reflection coefficient are used as prior constraints, and a nonlinear inversion target function is constructed;
(3) calculating porosity adjustment parameters by using the logging data information, and giving prior information such as constraint parameters, maximum inversion iteration times, initial porosity solution and the like;
(4) and (3) substituting the known parameters into a nonlinear inversion target function, and combining a rapid simulated annealing algorithm to realize high-precision direct quantitative inversion of porosity.
Further, 1, the quantitative relationship between porosity and seismic data
Based on the Biot elastic pore medium theory, the reflection coefficient of the fast longitudinal wave vertically incident on the interface of the pore medium containing saturated fluid is approximate,
wherein r is the reflection coefficient of fast longitudinal wave, Kl,μ l1,2 denote the bulk and shear modulus of the rock skeleton, αl,M l1,2 is the Biot parameter; vpljWhere 1,2, j 1,2 is the longitudinal wave velocity, γ lj1,2, j 1,2 is the ratio of the displacement amplitude of the pore fluid relative to the rock skeleton and the displacement amplitude of the rock skeleton; the subscript l is 1,2 represents the upper and lower pore media, respectively, and j is 1,2 represents the fast and slow longitudinal waves, respectively.
There is a link between the modulus of the rock skeleton, the modulus of the rock matrix and the porosity, and the invention uses the formulas (8) and (9) to describe the quantitative relationship between them.
Kd＝(1.0-f)Ks/(1.0+ck·f), (8)
μd＝(1.0-f)μs/(1.0+cμ·f)， (9)
Wherein f is porosity, Kd、μdRespectively the bulk modulus and shear modulus of the rock skeleton, Ks、μsThe bulk modulus and shear modulus of the rock matrix, respectively. c. Ck、cμParameters are adjusted for the pore media.
The fast longitudinal wave and the slow longitudinal wave exist in the pore medium containing the saturated fluid, and the slow longitudinal wave is quickly attenuated and disappears when propagating in the pore medium, so that the seismic record received by the ground geophone is actually the reflection of the fast longitudinal wave. The porosity f is a function f (t) of time t, and the fast longitudinal reflection coefficient r is a function r [ f (t) ] of the porosity f (t) as shown in equations (1) to (9). According to the convolution formula, the fast compressional synthesis seismic record d [ f (t) ] directly related to porosity can be expressed as,
d[f(t)]＝w(t)*r[f(t)]+n(t), (10)
the above equation is discretized and written into the form of a vector,
dN×1＝GN×N·rN×1+nN×1. (11)
wherein w (t) is a time domain wavelet, N (t) is noise, subscript Nx1 represents an N-dimensional column vector, NxN represents an N-dimensional matrix, and N is the number of sampling points of a seismic data; the symbol' denotes a convolution operator.
Further, 2, non-linear inversion of the objective function
According to the forward modeling process in the formula (11), the priori information in the well logging is taken as constraint, the fluid-containing pore medium porosity inversion objective function F is,
wherein f is the predicted porosity, R is the predicted reflection coefficient, and d is the synthetic seismic record; f. of0And R0Porosity and reflection coefficient, respectively, from well log datadataWhich is the actual seismic data, β is the reflection coefficient constraint parameter, gamma is the porosity constraint parameter, all vectors in the objective function are nx1 dimensional column vectors,
representing the L2 norm of the vector.
In the above formula, the first term is the fitting degree of observed seismic data and synthetic seismic records, the second term represents the similarity between the reflection coefficient calculated by the predicted porosity and the logging reflection coefficient, and the third term represents the similarity between the predicted porosity and the logging porosity.
Further, 3, simulated annealing algorithm porosity inversion
The simulated annealing algorithm is based on the solid annealing principle and simulates the annealing process of crystals in physics. The solids are first warmed to a sufficient level and then allowed to cool slowly. When the temperature is increased, the internal particles of the solid become an unordered state along with the temperature rise, and the internal energy is increased; when the particles are cooled, the particles gradually get ordered, and the temperature of each particle reaches an equilibrium state; finally, the ground state is reached at normal temperature, and the internal energy is reduced to the minimum. The algorithm has progressive convergence and is a global optimization algorithm which converges with the probability 1.
(1) Metropolis convergence criterion
The inversion algorithm obtains an optimal solution by continually iteratively updating the porosity. At the i-th iteration, the porosity f is first calculatediCorresponding objective function FiAnd then the porosity is updated to be fi+1And calculating the corresponding objective function value Fi+1. The difference between the two function values is as follows,
ΔF＝Fi+1-Fi,i≥1, (13)
if the delta F is less than 0, the porosity updating direction enables the objective function value to be reduced, and then the porosity is modified; if the delta F is larger than or equal to 0, further judging whether the delta F meets the formula (14), if so, still accepting the modification of the porosity, and if not, rejecting the modification.
0＜p＜K, (14)
Wherein, p is probability density, K is a random number, and the value range is 0 < K < 1; k is a radical ofbThe method is a Boltzmann constant, and 1, exp {. cndot } is taken as an exponential function in practical application. T is the current annealing temperature, often expressed exponentially,
Ti＝T0·(0.95)i-1,i≥1, (16)
wherein, T0Is the initial temperature, TiThe ith annealing temperature.
(2) Iteration and update of solutions
The solution of the simulated annealing algorithm is independent of the initial state, but research experience shows that: if a proper initial solution is given, the iterative solving times of the algorithm can be reduced, and the operation speed of the algorithm can be improved. Herein porosity f will be logged0Smoothed to give an initial porosity f1. Given an initial solution, a temperature-dependent similarity may be employedCauchy distribution generates an update vector Δ fiI.e. by
fi+1＝fi+ξ·Δfi,i≥1, (17)
Δfi＝T·sign(q-0.5)·[(1+1/T)|2·q-1|-1]. (18)
In the above equation, ξ is the update coefficient, fiFor the current prediction of porosity, fi+1For updated porosity, the value range is fi+1∈[0,0.3](ii) a q is obedience [0,1]And N-dimensional random vectors of the distribution, wherein sign is a sign function.
The advantage of using a temperature dependent Cauchy distribution to generate the update vector is: the porosity is searched in a large range under the high temperature condition, and the searching is performed near the current predicted value as far as possible under the low temperature condition, and a flat tail is distributed like Cauchy, so that the local minimum value is easy to jump out.
(3) Nonlinear porosity inversion process
Assuming the existence of seismic data d to be inverteddataWell logging reflection coefficient R0Logging porosity f0Initial porosity f1The maximum iteration number L and the logging constraint parameters β and gamma, the process for realizing high-precision porosity inversion based on the simulated annealing algorithm is as follows,
1) porosity f for the ith iterationiCalculating a new porosity value f using equation (17)i+1Then, based on the formula (12), calculating the inverse objective function value F before and after the porosity updatingiAnd Fi+1；
2) Calculating the difference delta F between the two objective function values by combining the formula (13), judging whether the condition delta F is less than 0, if so, accepting the updating of the porosity value, otherwise, continuously judging whether the formula (14) is met, if so, still accepting the updating of the porosity value, otherwise, refusing the updating;
3) judging whether the maximum iteration frequency L is reached, if so, exiting the loop to obtain an optimal porosity value; otherwise, iteratively updating the current annealing temperature by using the formula (16), and returning to 1) to continuously and iteratively update the porosity value f.
Experiments are listed below to illustrate the technical effect of improving the resolution ratio by performing porosity inversion by using the method of the present invention.
1. Model testing
(1) Single-channel sandstone model
And establishing a sandstone pore medium model, and verifying the effectiveness and stability of the method by analyzing the value of the initial solution and different constraint parameters. FIG. 2 shows a single sandstone pore medium model, wherein (a), (b), (c) and (d) respectively show the porosity, reflection coefficient, synthetic seismic record and low-frequency porosity of the model. The sampling rate of the model is 2ms, the delay is 0.6s, and other parameters include: the bulk modulus, shear modulus and density of the sandstone matrix are respectively 37GPa, 44GPa and 2650kg/m3(ii) a The permeability and porosity factor of sandstone are 0.2darcy and 2.0 respectively; adjustment parameter c of pore mediumkAnd cμAre all constant 6.0.
(2) Inversion results corresponding to different initial solutions
The invention selects low frequency porosity (fig. 2d) and a list of random vectors as initial solutions to perform porosity inversion. Using the ratio of the energy of the predicted porosity and the error in the porosity log to the energy of the porosity log, i.e.
The convergence of the inversion results is described. FIG. 3 is an energy ratio of porosity for different initial solutions, where the dashed black line indicates the initial solution is low frequency porosity and the solid black line indicates the randomly generated initial solution. Fig. 4, 5, and 6 are respectively the inversion porosity, the inversion reflection coefficient, and the reconstructed seismic record under different initial solutions, where the graphs (a), (b), (c), and (d) respectively represent the 100 th, 200 th, 300 th, and 400 th iterative inversions, the black solid line represents the original model, and the black dotted line and the dotted line respectively represent the inversion or reconstruction results corresponding to the initial solution with low-frequency porosity and randomly generated solution.
As can be seen from the analysis of fig. 4 to 6, under the condition that the initial solution is not determined, the inversion method can still obtain a relatively ideal porosity inversion result, and the inversion reflection coefficient and the reconstructed seismic record are well matched with the original model, which indicates that the method is suitable for inversion without the initial solution, and the inversion result is effective and stable. When the inversion is carried out by adopting the initial solution generated randomly, although the algorithm adopts the Metropolis convergence criterion to avoid falling into the local minimum value solution, the inversion needs more iteration times and the inversion speed is relatively slow; when the low-frequency porosity is selected as the initial solution, the inversion porosity is high in resolution, the optimal solution can be quickly converged, and the inversion speed is improved, as shown in fig. 3.
(3) Inversion result test corresponding to different parameters
The invention analyzes the influence of the two constraint parameters on the porosity inversion result after 600 times of inversion iteration when the two constraint parameters are different values, fig. 8 shows the corresponding porosity inversion result (dotted line) when the parameter β is 0.3 and the gamma is respectively 0.005(a), 0.01(b), 0.1(c) and 0.5(d), the black solid line is the original porosity model, fig. 7 shows the porosity energy ratio corresponding to the parameter in fig. 8, wherein the black solid line, the dotted line, the dash-dot line, the short line and the asterisk continuous line respectively correspond to the values of the parameter gamma which are 0.005, 0.01, 0.1 and 0.5.
As can be seen from fig. 7, after 600 times of iterative inversion, the inversion results corresponding to different constraint parameters all converge to the porosity model. However, the optimal inversion result can be obtained only by properly selecting the constraint parameters of the logging reflection coefficient and the logging porosity. And gradually compensating low-frequency components in the porosity inversion result along with the gradual increase of the constraint parameter gamma, and finally obtaining the optimal inversion porosity.
2. Actual seismic data inversion results
And intercepting part of actual seismic data of the A well, and verifying the effectiveness of the invention. FIG. 9 shows the actual seismic data and porosity inversion results, (a) the actual seismic section through well A, time domain depth 3.2s-4.1s, sampling interval 4ms, CDP 400; (b) the method is characterized in that the method is a low-frequency porosity model and is constructed by logging porosity, horizons and seismic data; (c) is a porosity inversion result corresponding to the traditional method; (d) for the porosity inversion result after 2100 iterations, the constraint parameter of the logging reflection coefficient in the inversion is 0.3, and the constraint parameter of the logging porosity is 0.6.
When the actual seismic data are inverted, the contribution degree of the actual seismic data to the porosity inversion result can be increased by properly reducing the constraint parameters of logging, so that the inversion result is more dependent on the actual seismic data and is more real and reliable. As can be seen from fig. 9: after 2100 times of iterative inversion, each inversion converges to an optimal solution, and an ideal porosity inversion result is obtained; compared with the traditional porosity inversion method, the method has higher porosity resolution and better coincidence degree with the porosity of the well logging, and the effectiveness of the method is verified. The inversion of the porosity provides guiding significance for the prediction of lithology and pore fluid in the oil-gas reservoir, and provides theoretical and data support for the exploration and development of oil and gas fields.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (4)
1. A porosity inversion method for improving resolution is characterized by comprising the following steps: assuming that seismic data to be inverted, a logging reflection coefficient and logging porosity exist, utilizing a simulated annealing algorithm to realize porosity inversion:
(1) deducing a fast longitudinal wave reflection coefficient approximate expression vertically incident to a fluid-containing pore medium interface, and establishing a direct quantitative relation between the porosity and seismic records based on a rock physics theory;
(2) according to the relation between the porosity and the seismic record, the logging porosity and the logging reflection coefficient are used as prior constraints, and a nonlinear inversion target function is constructed;
(3) calculating porosity adjustment parameters by using the logging data information, and giving prior information such as constraint parameters, maximum inversion iteration times, initial porosity solution and the like;
(4) and (3) substituting the known parameters into a nonlinear inversion target function, and combining a rapid simulated annealing algorithm to realize high-precision direct quantitative inversion of porosity.
2. The improved resolution porosity inversion method of claim 1, wherein the step (1) comprises: quantitative relationship between porosity and seismic data:
based on the Biot elastic pore medium theory, the reflection coefficient approximation formula of the fast longitudinal wave vertically incident on the interface of the pore medium containing the saturated fluid is as follows:
wherein r is the reflection coefficient of fast longitudinal wave, Kl,μl1,2 denote the bulk and shear modulus of the rock skeleton, αl,Ml1,2 is the Biot parameter; vpljWhere 1,2, j 1,2 is the longitudinal wave velocity, γlj,l＝1,2J is 1,2 is the ratio of the displacement amplitude of the pore fluid relative to the rock skeleton and the displacement amplitude of the rock skeleton; subscript l is 1,2 represents upper and lower pore medium, j is 1,2 represents fast and slow longitudinal wave;
the rock skeleton modulus, the rock matrix modulus and the porosity are linked, and the quantitative relation between the rock skeleton modulus, the rock matrix modulus and the porosity is described by the formulas (8) and (9),
Kd＝(1.0-f)Ks/(1.0+ck·f), (8)
μd＝(1.0-f)μs/(1.0+cμ·f)， (9)
wherein f is porosity, Kd、μdRespectively the bulk modulus and shear modulus of the rock skeleton, Ks、μsBulk and shear modulus of the rock matrix, ck、cμAdjusting parameters for the pore media;
the fast longitudinal wave and the slow longitudinal wave exist in the saturated fluid pore medium, the slow longitudinal wave is rapidly attenuated and disappears when propagating in the pore medium, therefore, the seismic record received by the ground detector is actually the reflection of the fast longitudinal wave, the porosity f is a function f (t) of time t, the reflection coefficient r of the fast longitudinal wave is a function r [ f (t) ] of the porosity f (t) according to the formulas (1) to (9), and according to the convolution formula, the fast longitudinal wave synthetic seismic record d [ f (t) ] directly related to the porosity can be expressed as,
d[f(t)]＝w(t)*r[f(t)]+n(t), (10)
the above equation is discretized and written into the form of a vector,
dN×1＝GN×N·rN×1+nN×1. (11)
wherein w (t) is a time domain wavelet, N (t) is noise, subscript Nx1 represents an N-dimensional column vector, NxN represents an N-dimensional matrix, and N is the number of sampling points of a seismic data; the symbol' denotes a convolution operator.
3. The improved resolution porosity inversion method according to claim 1, wherein the step (2) specifically comprises:
according to the forward modeling process in the formula (11), the priori information in the well logging is taken as constraint, the fluid-containing pore medium porosity inversion objective function F is,
wherein f is the predicted porosity, R is the predicted reflection coefficient, and d is the synthetic seismic record; f. of0And R0Porosity and reflection coefficient, respectively, from well log datadataWhich is the actual seismic data, β is the reflection coefficient constraint parameter, gamma is the porosity constraint parameter, all vectors in the objective function are nx1 dimensional column vectors,
an L2 norm representing a vector;
in the above equation, the first term is the fitting degree of the observed seismic data and the synthetic seismic record, the second term represents the similarity between the reflection coefficient calculated from the predicted porosity and the logging reflection coefficient, the third term represents the similarity between the predicted porosity and the logging porosity, and the latter two terms are the logging constraints on the inversion result.
4. The improved resolution porosity inversion method of claim 1,
the step (4) specifically comprises:
(41) metropolis convergence criteria:
the inversion algorithm obtains the optimal solution by continuously iterating and updating the porosity, and in the ith iteration, the porosity f is firstly calculatediCorresponding objective function FiAnd then the porosity is updated to be fi+1And calculating the corresponding objective function value Fi+1The difference between the two function values is as follows,
ΔF＝Fi+1-Fi,i≥1, (13)
if the delta F is less than 0, the porosity updating direction enables the objective function value to be reduced, and then the porosity is modified; if the delta F is larger than or equal to 0, further judging whether the delta F meets the formula (14), if so, still accepting the modification of the porosity, and if not, refusing the modification;
0＜p＜K, (14)
wherein, p is probability density, K is a random number, and the value range is 0 < K < 1; k is a radical ofbIs a Boltzmann constant, 1 is taken as an exponential function in practical application, exp {. cndot.), T is the current annealing temperature and is often expressed in an exponential form,
Ti＝T0·(0.95)i-1,i≥1, (16)
wherein, T0Is the initial temperature, TiIs the ith annealing temperature;
(42) iteration and updating of the solution:
porosity f will be measured0Smoothed to give an initial porosity f1Given an initial solution, an update vector Δ f is generated using a temperature-dependent Cauchy distributioniI.e. by
fi+1＝fi+ξ·Δfi,i≥1, (17)
Δfi＝T·sign(q-0.5)·[(1+1/T)|2·q-1|-1]. (18)
In the above equation, ξ is the update coefficient, fiFor the current prediction of porosity, fi+1For updated porosity, the value range is fi+1∈[0,0.3](ii) a q is obedience [0,1]Distributed N-dimensional random vectors, sign is a sign function;
(43) nonlinear porosity inversion flow:
assuming the existence of seismic data d to be inverteddataWell logging reflection coefficient R0Logging porosity f0Initial porosity f1The maximum iteration number L and the logging constraint parameters β and gamma, and the process for realizing the porosity inversion based on the simulated annealing algorithm is as follows:
431) porosity f for the ith iterationiCalculating a new porosity value f using equation (17)i+1Then, based on the formula (12), calculating the inverse objective function value F before and after the porosity updatingiAnd Fi+1；
432) Calculating the difference delta F between the two objective function values by combining the formula (13), judging whether the condition delta F is less than 0, if so, accepting the updating of the porosity value, otherwise, continuously judging whether the formula (14) is met, if so, still accepting the updating of the porosity value, otherwise, refusing the updating;
433) judging whether the maximum iteration frequency L is reached, if so, exiting the loop to obtain an optimal porosity value; otherwise iteratively update the current annealing temperature using equation (16) and return to 431) continue to iteratively update the porosity value f.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911288767.8A CN111077571B (en) | 2019-12-12 | 2019-12-12 | Porosity inversion method for improving resolution |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911288767.8A CN111077571B (en) | 2019-12-12 | 2019-12-12 | Porosity inversion method for improving resolution |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111077571A true CN111077571A (en) | 2020-04-28 |
| CN111077571B CN111077571B (en) | 2020-11-06 |
Family
ID=70314527
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911288767.8A Active CN111077571B (en) | 2019-12-12 | 2019-12-12 | Porosity inversion method for improving resolution |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111077571B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112630835A (en) * | 2020-12-03 | 2021-04-09 | 重庆三峡学院 | High-resolution post-stack seismic wave impedance inversion method |
| CN113189645A (en) * | 2021-05-19 | 2021-07-30 | 中海石油(中国)有限公司深圳分公司 | Matrix mineral modulus determination method and device, electronic equipment and storage medium |
| CN114063149A (en) * | 2021-11-16 | 2022-02-18 | 成都理工大学 | Inverse Q filtering method and system for improving seismic resolution |
| CN114486680A (en) * | 2022-01-24 | 2022-05-13 | 北京市生态环境保护科学研究院 | Method and device for determining permeability coefficient of deep unsaturated zone rock soil and electronic equipment |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102636812A (en) * | 2012-04-18 | 2012-08-15 | 中国石油天然气股份有限公司 | Method for obtaining volume of reserving space of carbonate reservoir |
| CN103852785A (en) * | 2012-11-28 | 2014-06-11 | 中国石油集团长城钻探工程有限公司 | Evaluation method for stratum anisotropy |
| CN106405654A (en) * | 2016-10-26 | 2017-02-15 | 成都理工大学 | Seismic spectrum imaging method based on deconvolution generalized S transform |
| US9989661B2 (en) * | 2011-09-26 | 2018-06-05 | Saudi Arabian Oil Company | Methods for evaluating rock properties while drilling using drilling rig-mounted acoustic sensors |
| CN109839662A (en) * | 2019-01-24 | 2019-06-04 | 长江大学 | A kind of sandstone reservoir recognition methods |
-
2019
- 2019-12-12 CN CN201911288767.8A patent/CN111077571B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9989661B2 (en) * | 2011-09-26 | 2018-06-05 | Saudi Arabian Oil Company | Methods for evaluating rock properties while drilling using drilling rig-mounted acoustic sensors |
| CN102636812A (en) * | 2012-04-18 | 2012-08-15 | 中国石油天然气股份有限公司 | Method for obtaining volume of reserving space of carbonate reservoir |
| CN103852785A (en) * | 2012-11-28 | 2014-06-11 | 中国石油集团长城钻探工程有限公司 | Evaluation method for stratum anisotropy |
| CN106405654A (en) * | 2016-10-26 | 2017-02-15 | 成都理工大学 | Seismic spectrum imaging method based on deconvolution generalized S transform |
| CN109839662A (en) * | 2019-01-24 | 2019-06-04 | 长江大学 | A kind of sandstone reservoir recognition methods |
Non-Patent Citations (1)
| Title |
|---|
| 夏红敏，等: "区域特性约束下的油藏物性模拟", 《地球物理学进展》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112630835A (en) * | 2020-12-03 | 2021-04-09 | 重庆三峡学院 | High-resolution post-stack seismic wave impedance inversion method |
| CN112630835B (en) * | 2020-12-03 | 2023-10-17 | 重庆三峡学院 | High-resolution post-stack seismic wave impedance inversion method |
| CN113189645A (en) * | 2021-05-19 | 2021-07-30 | 中海石油(中国)有限公司深圳分公司 | Matrix mineral modulus determination method and device, electronic equipment and storage medium |
| CN114063149A (en) * | 2021-11-16 | 2022-02-18 | 成都理工大学 | Inverse Q filtering method and system for improving seismic resolution |
| CN114486680A (en) * | 2022-01-24 | 2022-05-13 | 北京市生态环境保护科学研究院 | Method and device for determining permeability coefficient of deep unsaturated zone rock soil and electronic equipment |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111077571B (en) | 2020-11-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111077571B (en) | Porosity inversion method for improving resolution | |
| Zhong et al. | Inversion of time‐lapse seismic reservoir monitoring data using CycleGAN: A deep learning‐based approach for estimating dynamic reservoir property changes | |
| US20210264262A1 (en) | Physics-constrained deep learning joint inversion | |
| Azevedo et al. | Geostatistical seismic Amplitude‐versus‐angle inversion | |
| Misra et al. | Global optimization with model-space preconditioning: Application to AVO inversion | |
| CN110618450B (en) | Intelligent gas-bearing property prediction method for tight reservoir based on rock physical modeling | |
| Stephen et al. | Improved normalization of time‐lapse seismic data using normalized root mean square repeatability data to improve automatic production and seismic history matching in the Nelson field | |
| Kaur et al. | Time-lapse seismic data inversion for estimating reservoir parameters using deep learning | |
| Aleardi | Using orthogonal Legendre polynomials to parameterize global geophysical optimizations: applications to seismic‐petrophysical inversion and 1D elastic full‐waveform inversion | |
| Azevedo et al. | Stochastic perturbation optimization for discrete-continuous inverse problems | |
| CN113874760A (en) | Generating velocity and density models of subsurface structures of a reservoir | |
| WO2022187685A1 (en) | Method and system for faster seismic imaging using machine learning | |
| US20210311223A1 (en) | System and method for seismic inversion | |
| Sun et al. | Intelligent ava inversion using a convolution neural network trained with pseudo-well datasets | |
| Mgimba et al. | Application of GMDH to Predict Pore Pressure from Well Logs Data: A Case Study from Southeast Sichuan Basin, China | |
| Bazulin et al. | Determination of the elastic parameters of a VTI medium from sonic logging data using deep learning | |
| Pradhan et al. | Approximate Bayesian inference of seismic velocity and pore-pressure uncertainty with basin modeling, rock physics, and imaging constraints | |
| Bao-Li et al. | Prestack seismic stochastic inversion based on statistical characteristic parameters | |
| Meng et al. | Seismic impedance inversion using a multi‐input neural network with a two‐step training strategy | |
| US20210215841A1 (en) | Bandwith Extension of Geophysical Data | |
| Li et al. | 2D high-resolution crosswell seismic traveltime tomography | |
| Li et al. | AVO inversion in orthotropic media based on SA-PSO | |
| Anthony et al. | An optimized staggered-grid finite-difference operator for seismic wave simulation in poroelastic media | |
| Stephen | Seismic history matching with saturation indicators combined with multiple objective function optimization | |
| Assari et al. | Shear wave velocity prediction from petrophysical logs using MLP-PSO algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |