Reservoir parameter multichannel direct inversion method

Abstract

The invention provides a multi-channel direct inversion method for reservoir parameters. The inversion method comprises the following steps: inputting angle superposition seismic data, angle wavelets, a reservoir parameter inversion initial model and logging data; inputting a seismic interpretation horizon, analyzing and setting a deposition mode between layers, and calculating a sequence micro-layer grid; constructing a reservoir parameter inversion objective function constrained by the sequence stratigraphic grid based on the logarithmic linear rock physical model and the Bayesian theory; preferably selecting a space constraint weight based on a sequence stratigraphic grid; solving an objective function based on block coordinate descent and iterative reweighting; and outputting inversion results of the sand content, the total porosity and the oil and gas saturation. The method has the advantages of high inversion stability, strong anti-noise capability and good transverse continuity, and can meet the requirements of geological research.

Description

Technical Field

[0001] The invention relates to the field of reservoir inversion of oil and gas exploration seismic data, and in particular to a reservoir parameter multi-channel direct inversion method. Background Art

[0002] Compared with elastic parameters, reservoir parameters can directly characterize the lithology, physical properties and oil and gas content of underground media, and have clearer geological significance. Traditional reservoir parameter inversion methods are generally single-channel processing, which is easily affected by noise factors, resulting in poor stability of inversion results, low signal-to-noise ratio and poor lateral continuity, making it difficult to directly apply.

[0003] Considering that underground reservoir parameters show spatial correlation, that is, the properties of reservoir parameters with similar spatial positions are relatively close, spatial constraints can be constructed to improve the stability and lateral continuity of the inversion results.

[0004] At present, there are two common ways to introduce structural and sedimentary characteristic information in the spatial constraint term: (1) guiding the spatial constraint of the inversion result based on the local dip information of the post-stack seismic data waveform; (2) guiding the spatial constraint of the inversion result based on the local dip information of the interpreted horizon. In fact, both methods characterize the spatial correlation between any two points in the neighborhood by estimating the local dip information of the stratum and combining it with the coordinate rotation technology. However, this type of technology has certain defects. The accuracy of the local dip information estimated based on the waveform of post-stack seismic data is often affected by the signal-to-noise ratio of the seismic data and the size of the cross-correlation time window, and additional calculations are required to calculate the local dip. By using the local dip information of the layer instead of the calculation of the dip attribute, the amount of calculation can be reduced to a certain extent and the robustness of the dip estimation can be improved. However, both methods cannot avoid the strict condition that the inversion time window size of each CDP point is consistent, otherwise the spatial constraint term based on coordinate rotation cannot be constructed. This is not a problem for processing model data inversion or actual data inversion with a fixed time window. In actual situations, the time difference between two sets of layers is often space-varying. When dealing with such a variable time window determined by two layers, it becomes very difficult to impose such spatial constraints. These two types of spatial constraints often fail to achieve the purpose of reasonable spatial constraints when dealing with the inversion problem of sedimentary reservoirs with fast spatial variation, which will reduce the resolution and accuracy of the inversion results to a certain extent and make it difficult to meet the needs of actual production. Summary of the invention

[0005] In view of the above problems, the present invention is proposed to provide a reservoir parameter multi-channel direct inversion method that overcomes the above problems or at least partially solves the above problems.

[0006] According to one aspect of the present invention, a reservoir parameter multi-channel direct inversion method is provided, the inversion method comprising:

[0007] Step 1: Input angle stack seismic data, angle wavelet, reservoir parameter inversion initial model and logging data;

[0008] Step 2: Input the seismic interpretation horizon, analyze and set the sedimentation pattern between layers, and calculate the sequence micro-layer grid;

[0009] Step 3: Construct a reservoir parameter inversion objective function based on a log-linear rock physics model and sequence stratigraphic grid constraints based on Bayesian theory;

[0010] Step 4: Optimizing spatial constraint weights based on sequence stratigraphic grids;

[0011] Step 5: Solve the objective function based on block coordinate descent and iterative reweighting;

[0012] Step 6: Output the inversion results of sand content, total porosity and oil and gas saturation.

[0013] Optionally, the step 1: inputting angle stack seismic data, angle wavelet, reservoir parameter inversion initial model and logging data specifically includes:

[0014] Input seismic data volumes stacked at different angles;

[0015] Angle wavelets extracted from the superimposed data at each angle;

[0016] Initial models of sand content, total porosity, and hydrocarbon saturation along with well log data.

[0017] Optionally, the step 2: inputting seismic interpretation horizons, analyzing and setting the sedimentation patterns between the layers, and calculating the sequence microlayer grids specifically includes: inputting the number of seismic interpretation horizons M≥2,

[0018] Based on the given high-precision seismic interpretation horizon and sedimentary pattern, the interlayers are isochronously divided to obtain the sedimentary microlayer grid.

[0019] Optionally, the deposition mode between the layers includes three optional situations: parallel to the top, parallel to the bottom, and parallel to the top and bottom.

[0020] Optionally, the step 3: constructing a reservoir parameter inversion objective function based on a log-linear rock physics model and a sequence stratigraphic grid constraint of the Bayesian theory specifically comprises:

[0021] Establishment of angular reflection coefficient R based on log-linear rock physics model and FattiAVO formula pp (θ) has an approximately linear relationship with reservoir parameters;

[0022] According to the angle reflection coefficient R pp (θ) and the angle wavelet W(θ) are convolved to establish the forward modeling formula from reservoir parameters to angle gathers;

[0023] Spatial constraints of inversion are introduced in combination with sequence stratigraphic grid.

[0024] Optionally, the angle reflection coefficient R is established based on the logarithmic linear rock physics model and the FattiAVO formula pp The approximate linear relationship between (θ) and reservoir parameters specifically includes:

[0025] Establishment of angular reflection coefficient R based on log-linear rock physics model and FattiAVO formula pp The approximate linear formula of (θ) and reservoir parameters is:

[0026]

[0027] Where: C s , φ and Sg are sand content, total porosity and oil and gas saturation respectively. A(θ), B(θ) and C(θ) are functions of the incident angle θ, the background P-wave velocity ratio and rock physical parameters. The calculation formula is:

[0028]

[0029] The rock physical parameter ω in formula (2) I ,ω J ,ω ρ , I , J , ρ , I , J and ψ ρ It is obtained by linear regression of the following log-linear rock physics model and logging data,

[0030]

[0031] In formula (3), I, J and ρ represent longitudinal wave impedance, transverse wave impedance and density respectively.

[0032] Optionally, the angle reflection coefficient R pp The forward modeling formula from reservoir parameters to angle gathers is established by convolution of the angle wavelet W(θ) and W(θ), which specifically includes:

[0033] According to the angle reflection coefficient R pp (θ) and the angle wavelet W(θ) are convolved to establish the forward modeling formula from reservoir parameters to angle gathers:

[0034] d=Gm (4)

[0035] in,

[0036] K is the number of angle-stacked seismic data volumes.

[0037] Optionally, the spatial constraint item introduced into the inversion in combination with the sequence stratigraphic grid specifically includes:

[0038] Combined with the sequence stratigraphic grid, the spatial constraint term can be expressed as:

[0039]

[0040] in, represents the inversion result of reservoir parameters of the jth channel resampled to the sequence stratigraphic grid, where represents all seismic traces within the inversion area, δ j It represents the peripheral tracks of the jth track. The range of the peripheral tracks is reduced to the four tracks directly adjacent to the jth track in the longitudinal and lateral directions and the contact line directions by combining the properties of the first-order Markov field;

[0041] The inversion results of reservoir parameters under the seismic grid are expressed in the spatial constraint terms:

[0042]

[0043] In formula (6), m j represents the inversion result of reservoir parameters of the jth channel under the seismic grid, R j Represents the resampling matrix of the jth channel. Through the resampling matrix, the m under the seismic grid are j Resample to

[0044] The resampling matrix R represents the linear interpolation operator in the form of a matrix, realizes linear interpolation by matrix multiplication, and realizes resampling the inversion results under the seismic grid to the sequence stratigraphic grid;

[0045] Assume that the sequence stratigraphic grid nodes to be interpolated are Where, subscript j represents the jth track, subscript k represents the kth point in the vertical direction of the sequence stratigraphic grid, and the corresponding time is and The vertically adjacent two seismic grid nodes are m j,n and m j,n+1 , the inversion result of the nth sampling point of the jth channel and the inversion result of the n+1th sampling point, the corresponding time is T j,n and T j,n+1 , then according to the definition of linear interpolation, we can get:

[0046]

[0047] The resampling matrix R can be expressed as the following double diagonal matrix form:

[0048]

[0049] Combining the spatial constraint term of formula (8) with formula (4), the single-channel inversion objective function can be written as the following reservoir parameter multi-channel direct inversion objective function based on sparse constraints, low-frequency constraints and spatial constraints:

[0050]

[0051] Among them, L is the low-pass filter matrix, and R is the resampling matrix, which resamples the inversion results to the sequence stratigraphic grid to facilitate spatial constraints;

[0052] The value of the sparse constraint coefficient μ is generally between 0.5 and 2, and the value range of the low-frequency constraint coefficient α is generally between 0.1 and 10. The default values ​​of the coefficients μ and α are generally set to 1. The quality of the inversion result is not highly dependent on the values ​​of these two parameters.

[0053] In formula (9), β is the spatial constraint coefficient, which is a 1×3 vector representing the spatial constraint coefficients of the three quantities of sand content, total porosity and oil and gas saturation in the inversion.

[0054] Optionally, the step 4: preferably based on the spatial constraint weight of the sequence stratigraphic grid specifically includes:

[0055] A survey line passing through the well in the seismic data is selected for inversion, and other parameters remain unchanged. The sand content, total porosity and oil and gas saturation spatial constraint term coefficients are set to default parameters, and the reservoir parameter inversion is performed using formula (9);

[0056] The signal-to-noise ratio, spatial continuity and consistency of the inversion results of sand content and hydrocarbon saturation are evaluated based on the geological background and well information to determine whether the spatial constraint capacity is over-constrained or under-constrained.

[0057] When adjusting the coefficients, it is necessary to keep the spatial constraint coefficient values ​​of sand content and hydrocarbon saturation equal, and keep the spatial constraint coefficient of total porosity about 10 times the sand content or hydrocarbon saturation, and use the newly set parameters to perform reservoir parameter inversion using formula (9);

[0058] Repeat the steps until the signal-to-noise ratio and spatial continuity of the inversion results of sand content and oil and gas saturation meet the requirements of reservoir research. On the basis of keeping the sand content and oil and gas saturation parameters unchanged, further fine-tune the value of the spatial constraint coefficient of total porosity. Since the sensitivity of total porosity is much greater than that of sand content and oil and gas saturation, the adjustment step length should be reduced as much as possible when adjusting the spatial constraint coefficient of this item.

[0059] The spatial constraint coefficients of sand content, total porosity and oil and gas saturation are further fine-tuned until the signal-to-noise ratio, spatial continuity and other aspects of the inversion results of the three are basically consistent and meet the needs of geological research.

[0060] Optionally, the step 5: solving the objective function based on block coordinate descent and iterative reweighting specifically includes:

[0061] According to the theoretical algorithm of block coordinate descent, when solving the multi-channel inversion objective function of formula (9), the reservoir parameter vectors of all channels except the j-th channel are fixed, and only the reservoir parameter vector of the j-th channel is updated. Formula (9) is transformed into the following formula for solution:

[0062]

[0063] Formula (10) is the quasi-single-channel inversion objective function to be solved. To solve this formula, since the three-variable Cauchy constraint is a nonlinear source, it is necessary to use the iterative reweighted least squares optimization algorithm. j Taking partial derivatives, we get:

[0064]

[0065] Q in formula (11) j is the inconsistency weighting matrix, and its expression is:

[0066]

[0067] In formula (12) [] x,y It means taking the element of the xth row and yth column of the matrix;

[0068] The iterative reweighted least squares algorithm is used to solve the inversion parameter m j include:

[0069] Given the inversion parameter m j The initial solution is

[0070] According to the inversion parameter m j And formula (12), construct the non-consistent weight matrix Q j ;

[0071] According to the non-consistency weighting matrix Qj , use Cholesky decomposition or conjugate gradient algorithms to solve the linear equations (11) and obtain the new inversion parameter m j ;

[0072] Repeat the steps 2 to 5 times in an iterative manner to complete the inversion parameter m j The solution.

[0073] The invention provides a multi-channel direct inversion method for reservoir parameters, the inversion method comprising: inputting angle stacking seismic data, angle wavelet, reservoir parameter inversion initial model and logging data; inputting seismic interpretation horizon, analyzing and setting the sedimentation mode between each layer and calculating the sequence micro-layer grid; constructing a reservoir parameter inversion objective function based on the logarithmic linear rock physics model and the sequence stratigraphic grid constraint of the Bayesian theory; optimizing the spatial constraint weight based on the sequence stratigraphic grid; solving the objective function based on block coordinate descent and iterative reweighting; and outputting the inversion results of sand content, total porosity and oil and gas saturation. The inversion method has the advantages of high inversion stability, strong anti-noise ability and good lateral continuity, and can meet the needs of geological research.

[0074] The above description is only an overview of the technical solution of the present invention. In order to more clearly understand the technical means of the present invention, it can be implemented according to the contents of the specification. In order to make the above and other purposes, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are listed below. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0076] Figure 1 A flow chart of a reservoir parameter multi-channel direct inversion method provided by an embodiment of the present invention;

[0077] Figure 2 To obtain a comparison diagram between sedimentary microlayer grid and seismic grid by isochronous subdivision of interlayers according to seismic interpretation horizon and sedimentary pattern;

[0078] Figure 3 The solution process of multi-channel inversion method for reservoir parameters based on block coordinate descent and iterative reweighting is presented;

[0079] Figure 4 The synthetic angle gathers are generated by reservoir parameters based on the exact Zoeppritz equation and using Ricker wavelet with a main frequency of 35 Hz.

[0080] Figure 5 In order to utilize the method of the present invention Figure 4 The results of sand content, total porosity and oil and gas saturation obtained by inverting the synthetic angle gathers. The black curve in the figure is the measured curve, and the gray curve is the inversion curve;

[0081] Figure 6 The figure is a comparison diagram of the sand content inversion results of the actual data of a certain area using the method of the present invention. The left figure is a single-channel inversion result with the spatial constraint coefficient in the inversion objective function set to zero, and the right figure is a multi-channel inversion result using the sequence stratigraphic grid constraint of the present invention;

[0082] Figure 7 A comparison diagram of the well sections of sand content inverted by the method of the present invention for the actual data of a certain area;

[0083] Figure 8 A well cross-section comparison diagram of total porosity inverted by the method of the present invention for actual data of a certain area;

[0084] Fig. 9 The invention is used to compare the well sections with the oil and gas saturation inverted by the method of the present invention for the actual data of a certain area. DETAILED DESCRIPTION

[0085] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0086] The terms "comprises" and "having" and any variations thereof in the description embodiments, claims and drawings of the present invention are intended to cover non-exclusive inclusions, for example, including a series of steps or units.

[0087] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments.

[0088] like Figure 1 As shown, the present invention provides a multi-channel direct inversion method for reservoir parameters based on sequence stratigraphic grid constraints, which specifically includes:

[0089] Step 1: Input angle stack seismic data, angle wavelet, reservoir parameter inversion initial model and logging data;

[0090] Step 2: Input the seismic interpretation horizon, analyze and set the sedimentation pattern between layers, and calculate the sequence micro-layer grid;

[0091] On the basis of completing step 1, input the high-precision seismic interpretation horizons, analyze and set the sedimentary patterns between the horizons, and calculate the stratigraphic microlayer grids, where the number of input seismic interpretation horizons M ≥ 2, and the sedimentary patterns between the layers include three optional situations: parallel to the top, parallel to the bottom, and parallel to the top and bottom. According to the given high-precision seismic interpretation horizons and sedimentary patterns, the layers are isochronously divided to obtain the sedimentary microlayer grids.

[0092] Step 3: construct the objective function of reservoir parameter inversion constrained by sequence stratigraphic grid based on log-linear rock physics model and Bayesian theory;

[0093] After completing step 2, linear regression is performed based on the following log-linear rock physics model and logging data to obtain the rock physics parameter ω I ,ω J ,ω ρ , I , J , ρ , I , J and ψ ρ .

[0094]

[0095] Based on the calculated rock physical parameters, the following functions A(θ), B(θ) and C(θ) about the incident angle θ, the background P-wave velocity ratio and the rock physical parameters are calculated:

[0096]

[0097] After calculating A(θ), B(θ) and C(θ), combined with the initial model, the reflection coefficient R is substituted into the following formula: pp R is calculated by the approximate linear formula of (θ) and reservoir parameters pp (θ),

[0098]

[0099] Get the angular reflection coefficient R pp (θ) is then convolved with the input angle wavelet W(θ), and the two are convolved to obtain the forward modeling of reservoir parameters to angle gathers.

[0100] d=Gm

[0101] Combine the input low-pass filter parameters and the calculated sedimentary microlayer grid, as well as the default values ​​of the sparse constraint term coefficient, low-frequency constraint term coefficient, and spatial constraint term coefficient, and bring them into the inversion objective function

[0102]

[0103] Step 4: Optimize the spatial constraint weight coefficient.

[0104] After completing step 3, perform the optimization process of the spatial constraint weight coefficient as follows:

[0105] (1) Select a survey line passing through the well in the seismic data for inversion, keep other parameters unchanged, set the sand content, total porosity and oil and gas saturation spatial constraint term coefficients to default parameters, and perform reservoir parameter inversion using the method of step 5;

[0106] (2) Based on the geological background and well information, the signal-to-noise ratio, spatial continuity and consistency of the inversion results of sand content and oil and gas saturation are evaluated, and the consistency of geological knowledge is determined to determine whether the spatial constraint capacity is over-constrained or under-constrained, so as to appropriately reduce or increase the spatial constraint coefficient. When adjusting the coefficient, it is generally necessary to keep the spatial constraint coefficient values ​​of sand content and oil and gas saturation equal, and keep the spatial constraint coefficient of total porosity about 10 times that of sand content or oil and gas saturation, and use the newly set parameters to perform reservoir parameter inversion;

[0107] (3) Repeat step (2) until the signal-to-noise ratio and spatial continuity of the inversion results of sand content and hydrocarbon saturation meet the requirements of reservoir research. On the basis of keeping the sand content and hydrocarbon saturation parameters unchanged, further fine-tune the value of the spatial constraint coefficient of total porosity. Since the sensitivity of total porosity is much greater than that of sand content and hydrocarbon saturation, the adjustment step size should be reduced as much as possible when adjusting the spatial constraint coefficient of this item.

[0108] (4) Further fine-tune the spatial constraint coefficients of sand content, total porosity, and hydrocarbon saturation until the signal-to-noise ratio, spatial continuity, etc. of the inversion results of the three are basically consistent and meet the needs of geological research.

[0109] Step 5: Solve the objective function using block coordinate descent and iterative reweighting;

[0110] After completing step 4, solve the objective function as follows:

[0111] (1) Given the inversion parameter m j The initial solution is

[0112] (2) According to the inversion parameter m j And the following formula to construct the inconsistency weighted matrix Q j ;

[0113]

[0114] (3) According to the inconsistency weighting matrix Q j, use Cholesky decomposition or conjugate gradient algorithms to solve the following linear equation and obtain the new inversion parameter m j ;

[0115]

[0116] (4) Repeat steps (2) and (3) 2 to 5 times in an iterative manner to complete the solution of the inversion parameter mj.

[0117] Complete three cycles, the first cycle is the number of iterations of the block coordinate descent algorithm, the second cycle is the traversal of the three-dimensional inversion work area, and the third cycle is the iterative reweighted least squares solution process in each solution process. Finally, the final output results of the inverted sand content, total porosity and hydrocarbon saturation are obtained.

[0118] The following specifically describes a reservoir parameter multi-channel direct inversion method based on sequence stratigraphic grid constraints provided by the present invention in combination with model data and actual data of a certain area. Figure 2-Figure 3 As shown, Figure 2 In order to obtain the sedimentary micro-layer grid by isochronous division of the layers according to the input actual seismic interpretation horizon and the sedimentary pattern obtained by analysis, the spatial constraint in the inversion objective function is performed on the sedimentary micro-layer grid. Figure 3 It is an iterative solution process based on block coordinate descent and iterative reweighting strategy for the inversion objective function. When using the block coordinate descent method to calculate and update the inversion result of the current track, it is necessary to input the information of the surrounding four tracks for spatial constraints. The iterative reweighted least squares algorithm is used to solve the inversion update result of the track. The number of iterations in the data processing of this example is set to 5. The block coordinate descent algorithm is used to solve each track on the angle stacked seismic data body. After all tracks are updated, the next round of data is updated until the preset number of updates is reached. Completing a round of updates can be regarded as completing an iteration. In this example, the number of iterations of the block coordinate descent algorithm is set to 5.

[0119] Figure 4 The angle gathers obtained by forward modeling are given. The reservoir parameters are inverted using this method. The inversion results are shown in Figure 5 shown.

[0120] from Figure 5 It can be seen that the trends of the inversion results and theoretical values ​​are well consistent, thus verifying the correctness of the method. The method of the present invention is used to invert reservoir parameters on actual data, and the results are as follows: Figure 6-9 shown.

[0121] Figure 6This is a comparison of the results of sand content inversion of actual data in a certain area using single-channel direct inversion of reservoir parameters and multi-channel direct inversion of reservoir parameters based on sequence stratigraphic grid constraints. It can be seen from the figure that the single-channel inversion result has a low signal-to-noise ratio and poor spatial continuity, while the multi-channel inversion result has a high signal-to-noise ratio and good spatial continuity. Further comparison with the fault development area in the elliptical frame in the figure, three sets of small faults can be clearly identified in the multi-channel inversion results, while the small faults in the single-channel inversion results are poorly presented and difficult to identify.

[0122] Figure 7 The well sections are compared with the results of sand content inversion of actual data using the multi-channel direct inversion method of reservoir parameters based on sequence stratigraphic grid constraints. It can be seen that the spatial continuity of the multi-channel inversion results of sand content is good, the resolution is high, and the consistency between the well bypass and the well is good.

[0123] Figure 8 This is a comparison of the well sections after the porosity inversion results of the actual data using the multi-channel direct inversion method of reservoir parameters based on sequence stratigraphic grid constraints. It can be seen from the figure that the spatial continuity of the porosity inversion results is good, and the wellside channel inversion results are highly consistent with the well curve.

[0124] Fig. 9 The results of the well section comparison after the oil and gas saturation inversion were performed using the multi-channel direct inversion method of reservoir parameters based on sequence stratigraphic grid constraints for the actual data of a certain area. From the comparison in the figure, it can be seen that the inversion results are generally consistent with the well saturation information. The reservoir parameters inverted by this method have high resolution and good spatial continuity, which can directly predict the favorable oil and gas development zones and reservoir research.

[0125] Beneficial effects: The present invention provides a multi-channel direct inversion method for reservoir parameters based on sequence stratigraphic grid constraints, which can improve the stability, continuity and noise resistance of reservoir parameter inversion. The method fully considers geological information such as the structural morphology and sedimentary characteristics of the underground medium, calculates the sequence micro-layer grid, and imposes spatial constraints on the inversion process on the sequence micro-layer grid based on the spatial similarity characteristics of the underground reservoir parameters, and adaptively integrates the stratigraphic structure and sedimentary pattern information, thereby improving the inversion quality of the reservoir parameters. The direct inversion of reservoir parameters using this method has the advantages of high stability, good continuity and strong noise resistance.

[0126] The above specific implementation methods further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above are only specific implementation methods of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A multi-channel direct inversion method for reservoir parameters, characterized in that: The inversion method comprises: Step 1: Input angle stack seismic data, angle wavelet, reservoir parameter inversion initial model and logging data; Step 2: Input the seismic interpretation horizon, analyze and set the sedimentation pattern between layers, and calculate the sequence micro-layer grid; Step 3: Construct a reservoir parameter inversion objective function based on a log-linear rock physics model and sequence stratigraphic grid constraints based on Bayesian theory; Step 4: Optimizing spatial constraint weights based on sequence stratigraphic grids; Step 5: Solve the objective function based on block coordinate descent and iterative reweighting; Step 6: Output the inversion results of sand content, total porosity and oil and gas saturation.

2. A reservoir parameter multi-channel direct inversion method according to claim 1, characterized in that: The step 1: inputting angle stacking seismic data, angle wavelet, reservoir parameter inversion initial model and logging data specifically includes: Input seismic data volumes stacked at different angles; Angle wavelets extracted from the superimposed data at each angle; Initial models of sand content, total porosity, and hydrocarbon saturation along with well log data.

3. A reservoir parameter multi-channel direct inversion method according to claim 1, characterized in that: The step 2: inputting the seismic interpretation horizon, analyzing and setting the sedimentation mode between each layer, and calculating the stratigraphic micro-layer grid specifically includes: inputting the number of seismic interpretation horizons M≥2, Based on the given high-precision seismic interpretation horizon and sedimentary pattern, the interlayers are isochronously divided to obtain the sedimentary microlayer grid.

4. A reservoir parameter multi-channel direct inversion method according to claim 1, characterized in that: The deposition modes between the layers include three optional situations: parallel to the top, parallel to the bottom, and parallel to the top and bottom.

5. A reservoir parameter multi-channel direct inversion method according to claim 1, characterized in that: The step 3: constructing a reservoir parameter inversion objective function based on a logarithmic linear rock physics model and a sequence stratigraphic grid constraint of the Bayesian theory specifically includes: Establishment of angular reflection coefficient R based on log-linear rock physics model and FattiAVO formula pp (θ) has an approximately linear relationship with reservoir parameters; According to the angle reflection coefficient R pp (θ) and the angle wavelet W(θ) are convolved to establish the forward modeling formula from reservoir parameters to angle gathers; Spatial constraints of inversion are introduced in combination with sequence stratigraphic grid.

6. A reservoir parameter multi-channel direct inversion method according to claim 5, characterized in that: The angle reflection coefficient R is established based on the logarithmic linear rock physics model and the FattiAVO formula pp The approximate linear relationship between (θ) and reservoir parameters specifically includes: Establishment of angular reflection coefficient R based on log-linear rock physics model and FattiAVO formula pp The approximate linear formula of (θ) and reservoir parameters is: Where: C s , φ and Sg are sand content, total porosity and oil and gas saturation respectively. A(θ), B(θ) and C(θ) are functions of the incident angle θ, the background P-wave velocity ratio and rock physical parameters. The calculation formula is: The rock physical parameter ω in formula (2) I ,ω J ,ω ρ , I , J , ρ , I , J and ψ ρ It is obtained by linear regression of the following log-linear rock physics model and logging data, In formula (3), I, J and ρ represent longitudinal wave impedance, transverse wave impedance and density respectively.

7. A reservoir parameter multi-channel direct inversion method according to claim 5, characterized in that: The angle reflection coefficient R pp The forward modeling formula from reservoir parameters to angle gathers is established by convolution of the angle wavelet W(θ) and W(θ), which specifically includes: According to the angle reflection coefficient R pp (θ) and the angle wavelet W(θ) are convolved to establish the forward modeling formula from reservoir parameters to angle gathers: d=Gm (4) in, K is the number of angle-stacked seismic data volumes.

8. A reservoir parameter multi-channel direct inversion method according to claim 5, characterized in that: The spatial constraints introduced into the inversion by combining the sequence stratigraphic grid specifically include: Combined with the sequence stratigraphic grid, the spatial constraint term can be expressed as: in, represents the inversion result of reservoir parameters of the jth channel resampled to the sequence stratigraphic grid, where represents all seismic traces within the inversion area, δj represents the peripheral traces of the jth trace, and the range of the peripheral traces is reduced to the four traces directly adjacent to the jth trace in the longitudinal and lateral directions and the contact line directions by combining the properties of the first-order Markov field; The inversion results of reservoir parameters under the seismic grid are expressed in the spatial constraint terms: In formula (6), m j represents the inversion result of reservoir parameters of the jth channel under the seismic grid, R j Represents the resampling matrix of the jth channel. Through the resampling matrix, the m under the seismic grid are j Resample to The resampling matrix R represents the linear interpolation operator in the form of a matrix, realizes linear interpolation by matrix multiplication, and realizes resampling the inversion results under the seismic grid to the sequence stratigraphic grid; Assume that the sequence stratigraphic grid nodes to be interpolated are Where, subscript j represents the jth track, subscript k represents the kth point in the vertical direction of the sequence stratigraphic grid, and the corresponding time is and The vertically adjacent two seismic grid nodes are m j,n and m j,n +1, the inversion result of the nth sampling point of the jth channel and the inversion result of the n+1th sampling point, the corresponding time is T j,n and T j,n+1 , then according to the definition of linear interpolation, we can get: The resampling matrix R can be expressed as the following double diagonal matrix form: Combining the spatial constraint term of formula (8) with formula (4), the single-channel inversion objective function can be written as the following reservoir parameter multi-channel direct inversion objective function based on sparse constraints, low-frequency constraints and spatial constraints: Among them, L is the low-pass filter matrix, and R is the resampling matrix, which resamples the inversion results to the sequence stratigraphic grid to facilitate spatial constraints; The value of the sparse constraint coefficient μ is generally between 0.5 and 2, and the value range of the low-frequency constraint coefficient α is generally between 0.1 and 10. The default values ​​of the coefficients μ and α are generally set to 1. The quality of the inversion result is not highly dependent on the values ​​of these two parameters. In formula (9), β is the spatial constraint coefficient, which is a 1×3 vector representing the spatial constraint coefficients of the three quantities of sand content, total porosity and oil and gas saturation in the inversion.

9. A reservoir parameter multi-channel direct inversion method according to claim 5, characterized in that: The step 4: preferably based on the spatial constraint weight of the sequence stratigraphic grid specifically comprises: A survey line passing through the well in the seismic data is selected for inversion, and other parameters remain unchanged. The sand content, total porosity and oil and gas saturation spatial constraint term coefficients are set to default parameters, and the reservoir parameter inversion is performed using formula (9); The signal-to-noise ratio, spatial continuity and consistency of the inversion results of sand content and hydrocarbon saturation are evaluated based on the geological background and well information to determine whether the spatial constraint capacity is over-constrained or under-constrained. When adjusting the coefficients, it is necessary to keep the spatial constraint coefficient values ​​of sand content and hydrocarbon saturation equal, and keep the spatial constraint coefficient of total porosity about 10 times the sand content or hydrocarbon saturation, and use the newly set parameters to perform reservoir parameter inversion using formula (9); Repeat the steps until the signal-to-noise ratio and spatial continuity of the inversion results of sand content and oil and gas saturation meet the requirements of reservoir research. On the basis of keeping the sand content and oil and gas saturation parameters unchanged, further fine-tune the value of the spatial constraint coefficient of total porosity. Since the sensitivity of total porosity is much greater than that of sand content and oil and gas saturation, the adjustment step length should be reduced as much as possible when adjusting the spatial constraint coefficient of this item. The spatial constraint coefficients of sand content, total porosity and oil and gas saturation are further fine-tuned until the signal-to-noise ratio, spatial continuity and other aspects of the inversion results of the three are basically consistent and meet the needs of geological research.

10. A reservoir parameter multi-channel direct inversion method according to claim 5, characterized in that: The step 5: solving the objective function based on block coordinate descent and iterative reweighting specifically includes: According to the theoretical algorithm of block coordinate descent, when solving the multi-channel inversion objective function of formula (9), the reservoir parameter vectors of all channels except the j-th channel are fixed, and only the reservoir parameter vector of the j-th channel is updated. Formula (9) is transformed into the following formula for solution: Formula (10) is the quasi-single-channel inversion objective function to be solved. To solve this formula, since the three-variable Cauchy constraint is a nonlinear source, it is necessary to use the iterative reweighted least squares optimization algorithm. j Taking partial derivatives, we get: Q in formula (11) j is the inconsistency weighting matrix, and its expression is: In formula (12) [] x,y It means taking the element of the xth row and yth column of the matrix; The iterative reweighted least squares algorithm is used to solve the inversion parameter m j include: Given the inversion parameter m j The initial solution is According to the inversion parameter m j And formula (12), construct the non-consistent weight matrix Q j ; According to the non-consistency weighting matrix Q j , use Cholesky decomposition or conjugate gradient algorithms to solve the linear equations (11) and obtain the new inversion parameter m j ; Repeat the steps 2 to 5 times in an iterative manner to complete the inversion parameter m j The solution.

Images