A thin reservoir prediction method based on frequency division phasing

Abstract

The application discloses a thin reservoir prediction method based on frequency division phase control, relates to the field of seismic reservoir prediction, and is based on geophysical methods as a basic research means, and aims at thin reservoirs to establish variable gamma prestack deterministic inversion based on big data waveform driven low-frequency models, carry out phase control prestack geostatistics inversion combination, accurately depict internal details of thin sand layers, solve the problem of fine depiction and identification of thin layer reservoirs, and combine amplitude attributes, oil and gas detection and other methods to comprehensively verify each other, solve the geological problems of yield increase and speed increase, reduce exploration risks, and have universal applicability and broad application prospect.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic reservoir prediction, more particularly, it relates to a thin reservoir prediction method based on frequency division phase control. BACKGROUND

[0002] At present, domestic and foreign scholars have carried out a large amount of research on thin reservoir identification and prediction. Generally, it can be divided into two categories: qualitative or semi-quantitative thin reservoir prediction technology based on geophysical attributes and single karst transformation thin reservoir prediction technology based on comprehensive geological description.

[0003] The qualitative or semi-quantitative thin reservoir prediction technology based on geophysical attributes is to extract seismic attributes to analyze the reservoir lithology and physical property characteristics reflected by the seismic attributes according to the large amount of geological information implied in the original seismic data, and to predict the sedimentation and oil and gas bearing property of the reservoir. Since there are various types of seismic attributes, how to extract and optimize the characteristic attributes becomes a difficult problem for realizing thin reservoir identification and prediction.

[0004] The single karst transformation thin reservoir prediction technology is mainly based on high-resolution seismic data processing and inversion. Since thin reservoir research needs to be carried out on the basis of high resolution, the existing seismic data often has narrow frequency band and low main frequency, which cannot meet the accuracy of fine reservoir description. In order to solve this kind of problem, two aspects are generally adopted. One is to carry out frequency expansion processing on the original data, but this method has uncertainty in amplitude preservation and fidelity, which affects the reliability of reservoir prediction. The other is to improve the imaging accuracy of seismic data by high-density seismic data acquisition and processing from the stage of seismic data acquisition, but the cycle is long and cannot meet the fast-paced demand of oil reservoir development.

[0005] In summary, each of the above thin reservoir prediction methods has certain feasibility and limitations. Therefore, with the continuous deepening of oilfield exploration and development, it is urgent to form a set of thin reservoir and unconventional reservoir inversion technology suitable for fast-paced oil reservoir development. SUMMARY

[0006] In view of the deficiencies of the prior art, the purpose of the present application is to provide a thin reservoir prediction method based on frequency division phase control, which establishes variable gamma prestack deterministic inversion based on big data waveform driving low frequency model for thin reservoir, carries out phase-controlled prestack geostatistical inversion combination, accurately depicts the internal details of thin sand layer, and solves the problem of fine depiction and identification of thin layer reservoir.

[0007] To achieve the above purpose, the present application provides the following technical scheme: a thin reservoir prediction method based on frequency division phase control, comprising the following steps:

[0008] S1, quality control and consistency processing of well logging data. The prestack seismic data are quality controlled and optimized to improve the accuracy of the prestack seismic data.

[0009] S2, a low frequency model is established according to the data processed in S1.

[0010] S3, based on the low frequency model, combined with the data processed in S1, improved variable gamma prestack joint deterministic inversion is carried out to obtain a low frequency and medium frequency model, and the formula is as follows:

[0011]

[0012] In the formula, t is a time sampling point, t0 is an initial time, τ = τ0(s-1), s is a parameter of medium heterogeneity, V P is a P-wave velocity, V s is a S-wave velocity, V ρ is a density, R p is a P-wave velocity reflection coefficient, R s is a S-wave velocity reflection coefficient, R ρ is a density reflection coefficient.

[0013] R p , R s , R ρ , which is obtained by the following formula:

[0014] R pp (θ) = (1 + tan 2 θ) R p - 8γsin 2 θ R s + (1-4γsin 2 θ) R ρ

[0015] In the formula, θ is an incident angle, γ is a P / S-wave velocity ratio, R pp (θ) is a P-wave reflection coefficient varying with angle, θ, γ and R pp (θ) are all obtained according to the data processed in S1.

[0016] S4, sand thickness is depicted according to the inversion result, and a plane sand ratio map is made.

[0017] S5, according to the inversion result of S3, a plane sand ratio map is used as a constraint to carry out geostatistical inversion, and the inversion result is expanded to high frequency to obtain a full frequency band geological model.

[0018] The application is further provided with: S2, the specific steps are, first, performing big data analysis on the seismic waveform in the S1 result to obtain low-frequency data; then, according to the P-wave impedance, the P-wave to S-wave velocity ratio and the density at the well point in the S1 result, the physical properties and structural features of the underground medium are described, the low-frequency data are limited, and a low-frequency model is obtained.

[0019] The application is further provided with: after the S2 low-frequency model is completed, it is necessary to analyze whether the low-frequency model is consistent with the geological understanding, if yes, S3 is performed, and if not, S1 is returned; the judgment of whether the low-frequency model is consistent with the geological understanding includes verification according to the S1 result and checking of the consistency of each part of the model.

[0020] The application is further provided with: when the S2 low-frequency model is established, unknown seismic attributes need to be calculated according to known seismic attributes; the same type of seismic attributes of different wells are obtained through a space conversion formula, and the formula is as follows:

[0021]

[0022] In the formula, S n and S i are seismic attributes of different wells, S i is obtained according to the S1 result, W is a wellbore relationship function, W is obtained by deduction according to the S1 result, and n is the number of wells.

[0023] The application is further provided with: when the S2 low-frequency model is established, unknown seismic attributes need to be calculated according to known seismic attributes; the same type of seismic attributes of different wells are obtained through a space conversion formula, and the formula is as follows:

[0024]

[0025] In the formula, L j and L i are seismic attributes of different types, L i is obtained according to the S1 result, λ ij is an attribute correlation matching function, λ ij is obtained by deduction according to the S1 result, and n is the number of wells.

[0026] The application is further provided with: when the S2 low-frequency model is established, unknown seismic attributes need to be calculated according to known seismic attributes; the same type of seismic attributes of different wells are obtained through a space conversion formula, and the formula is as follows:

[0027]

[0028] In the formula, SL j and SL i are seismic attributes of different types of different wells, wij for SL j and SL i corresponding inter-shaft relationship function, w ij obtained by inferring the result of S1. λ ij for attribute correlation matching function, λ ij obtained by inferring the result of S1, and n is the number of wells.

[0029] The application further provides that the geostatistical inversion in S5 adopts a variogram, and the variogram formula is as follows:

[0030]

[0031] wherein, h is a given step length, N(h) is the number of points with a distance equal to h, Z(x i ) is a measured value of a variable at point x i , Z(x i +h) is a measured value of the variable deviating from point x i by h, and r(h) is an inversion result.

[0032] The application further provides that the processing of the well logging data in S1 includes well logging curve error correction and multi-well consistency processing.

[0033] The application further provides that the well logging curve error correction includes the following steps:

[0034] For the layer section with obvious deformation of the wellbore, the Castagna formula or the Smiths formula is used for correction.

[0035] The Castagna formula is as follows: wherein, p is the formation density, Vp is the P-wave velocity, a, b and c are all constants, and are obtained according to the well logging data analysis of the non-deformed layer section of the wellbore.

[0036] The Smiths formula is as follows: S(v) = dR e wherein, S(v) is the acoustic velocity, R is the resistivity, and d and e are both constants, and are obtained according to the well logging data analysis of the non-deformed layer section of the wellbore.

[0037] For the layer section with obvious collapse of the wellbore, the volume percentage model of the well logging response is used for correction, and the formula is as follows:

[0038] wherein, Ci is the well logging response value of the non-collapsed well section of the i-th rock component, Vi is the volume percentage of the i-th rock component, and C is the forward well logging response value.

[0039] The application further provides that the multi-well consistency processing formula is as follows: In the formula, Z(x, y) is a seismic attribute to be studied, (x, y) is the coordinate of a point on a plane, is a trend value, and g is a correction value.

[0040] Regression analysis is performed on the known Z(x, y) in the logging data to obtain a regression equation f(x, y), and the trend value is obtained according to f(x, y) Z(x, y) is subtracted from to obtain g.

[0041] The application is further provided as follows: the processing of the pre-stack seismic data in S1 includes pre-stack random noise suppression and time difference correction.

[0042] The application is further provided as follows: the pre-stack random noise suppression adopts 3D pre-stack random noise attenuation of the GeoEast system.

[0043] The application is further provided as follows: the time difference correction of the pre-stack seismic data in S1 is performed by using a non-rigid matching method.

[0044] The application is further provided as follows: the planar sand ratio map is drawn by using a gridding interpolation method. Firstly, the planar sand ratio map is divided into a series of small grids, and the sand content value of each sample point is distributed to the corresponding grid node. Then, according to the spatial position and the sand content value of the known data points, the sand content value of each grid node is calculated by using an interpolation method, and the contour map of the sand content distribution is drawn.

[0045] The application is further provided as follows: the gridding interpolation method is a contour mapping linear interpolation method based on a triangular grid.

[0046] In summary, compared with the prior art, the application has the following beneficial effects: the application establishes variable-g pre-stack deterministic inversion based on a large data waveform driven low-frequency model for a thin reservoir, performs a phase-controlled pre-stack geostatistical inversion combination, accurately depicts the internal details of a thin sand layer, solves the problem of fine depiction and identification of a thin layer reservoir, and solves the geological problems of yield increase and speed improvement by combining multiple methods such as amplitude attribute and oil and gas detection for mutual verification, reduces the exploration risk, and has universal applicability and broad application prospect. BRIEF DESCRIPTION OF DRAWINGS

[0047] Figure 1 It is a flowchart of an embodiment;

[0048] Figure 2 It is a comparison before and after the pre-stack seismic data processing;

[0049] Figure 3 It is a fine calibration of well logging and seismic data of TK235 well;

[0050] Figure 4Big data-driven geological framework models;

[0051] Figure 5 Distribution map of the CK1-1 sand body in the eastern part of the Santamu area of ​​the Tarim Oilfield;

[0052] Figure 6 Distribution map of the CK5-2 sand body in the eastern part of the Santamu area of ​​the Tarim Oilfield;

[0053] Figure 7 A schematic diagram of the characterization of the TKC1-2H well in the Carboniferous thin reservoir of the Tahe River;

[0054] Figure 8 A schematic diagram of the comprehensive characterization results of the thin Carboniferous reservoir in the Tahe River Basin. Detailed Implementation

[0055] The technical solution of the present invention will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are not all embodiments of the present invention. All other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the invention.

[0056] Example

[0057] This embodiment was applied in the Santamu working area, the most abundant block of clastic oil and gas reservoirs in the Tarim Oilfield. Several Ordovician wells showed good hydrocarbon shows in the Carboniferous system, but large-scale development has not yet materialized. The Carboniferous Karasai Formation is a transitional zone between marine and continental environments, with frequent interbedded sandstone and mudstone, exhibiting an overall mud-encased sand structure. The reservoir is characterized by rapid lateral variation and small differences in sand and mud impedance, showing no obvious response characteristics on seismic profiles. Describing the distribution characteristics of thin reservoirs is difficult, thus limiting the effectiveness of reservoir development.

[0058] like Figure 1 The diagram shown is a flowchart of a preferred embodiment of the present invention. A thin reservoir prediction method based on frequency division phase control includes the following steps:

[0059] S1. Perform quality control and consistency processing on well logging data, including well logging curve error correction and multi-well consistency processing. For example... Figure 2 As shown, pre-stack seismic data undergoes quality control and optimization to improve its accuracy, including pre-stack random noise suppression and time difference correction.

[0060] S2. Based on the data processed in S1, establish a low-frequency model. Based on the results of S1, analyze the seismic waveforms, such as... Figure 3 As shown, big data analysis is performed to obtain low-frequency data; then, based on the P-wave impedance, P-wave velocity ratio, and density at the well point in the S1 results, a detailed interpretation of the standard layer and cross-section is conducted to characterize the physical properties and structural features of the subsurface medium. The low-frequency data is then constrained to obtain a low-frequency model, such as...Figure 4 as shown.

[0061] S3, based on the low frequency model, combined with the data processed in S1, improved variable gamma pre-stack joint deterministic inversion was carried out to obtain a low frequency and medium frequency model, and the formula was as follows:

[0062]

[0063] In the formula, t is a time sampling point, t0 is an initial time, τ = τ0(s-1), s is a parameter of medium heterogeneity, V P is a P-wave velocity, V s is a S-wave velocity, V ρ is a density, R p is a P-wave velocity reflection coefficient, R s is a S-wave velocity reflection coefficient, R ρ is a density reflection coefficient.

[0064] R p , R s , R ρ , which was obtained by the following formula:

[0065] R pp (θ) = (1 + tan 2 θ) R p - 8γsin 2 θ R s + (1-4γsin 2 θ) R ρ

[0066] In the formula, θ is an incident angle, γ is a P-S wave velocity ratio, R pp (θ) is a P-wave reflection coefficient varying with angle, θ, γ and R pp (θ) are all obtained according to the data processed in S1.

[0067] S4, according to the inversion result, sand body thickness was depicted, and a plane sand ratio map was made, as shown in Figure 5 , 6 The plane sand ratio map was drawn by gridding interpolation method; first, the plane sand ratio map was divided into a series of small grids, and the sand content value of each sample point was distributed to the corresponding grid node; then, according to the spatial position and sand content value of the known data points, the sand content value of each grid node was calculated by interpolation method, and the contour map of sand content distribution was drawn.

[0068] S5, as shown in Figure 7 , 8As shown, according to the S3 inversion result, a geostatistical inversion is performed with a planar sand body ratio map as a constraint to expand the inversion result to high frequency and obtain a full-band geologic model. The geostatistical inversion uses a variogram, and the variogram formula is as follows:

[0069]

[0070] In the formula, h is a given step length, N(h) is the number of points with a distance equal to h, Z(x i ) is a measured value of a variable at point x i , Z(x i +h) is a measured value of the variable at a point deviating from point x i by h, and r(h) is an inversion result.

[0071] Specifically, for the S1 result, it is necessary to determine whether it is consistent with well-seismic calibration, AVO rules, and geologic understanding matching. If it is consistent, S2 is performed. If it is not consistent, S1 is performed again. For S2, it is necessary to analyze whether the low-frequency model is consistent with geologic understanding. If it is consistent, S3 is performed. If it is not consistent, S1 is returned. For the determination of whether the low-frequency model is consistent with geologic understanding, it includes verification according to the S1 result and checking consistency of each part of the model. For the S5 result, it is necessary to analyze whether it meets geologic understanding. If it is consistent, the result is output. If it is not consistent, S3 is returned.

[0072] Specifically, when the S2 low-frequency model is established, unknown seismic attributes are calculated according to known seismic attributes.

[0073] Different types of seismic attributes of different wells are obtained through a spatial conversion formula, and the formula is as follows:

[0074]

[0075] In the formula, S n and S i are seismic attributes of different wells, S i is obtained according to the S1 result; W is a well-to-well relationship function, W is obtained by deducing the S1 result, and n is the number of wells.

[0076] Different types of seismic physical properties of the same well are obtained through a physical property conversion formula, and the formula is as follows:

[0077]

[0078] In the formula, L j and L i are different types of seismic physical properties, L i is obtained according to the S1 result; λ ij is an attribute correlation matching function, λ ij is obtained by deducing the S1 result, and n is the number of wells.

[0079] Different types of seismic physical properties of different wells are obtained by full-area physical parameter conversion formula, and the formula is as follows:

[0080]

[0081] In the formula, SL j and SL v are different types of seismic physical properties of different wells, w ij is the wellbore interval relationship function corresponding to SL j and SL i , w ij is obtained by S1 result, and n is the number of wells.

[0082] Specifically, the well logging curve error correction in S1 is that for the layer section with obvious deformation of the wellbore, the gamma curve, resistivity curve and the like are less affected by the wellbore collapse alteration and have better quality, and the Castagna formula or Smiths formula is used for correction.

[0083] Castagna formula or Smiths formula is used for correction.

[0084] Castagna formula: In the formula, p is the formation density, Vp is the longitudinal wave velocity, a, b and c are constants, and are obtained according to the well logging data analysis of the wellbore un-deformed layer section.

[0085] Smiths formula: S(v)=dR e In the formula, S(v) is the acoustic velocity, R is the resistivity, and d and e are constants, and are obtained according to the well logging data analysis of the wellbore un-deformed layer section.

[0086] For the layer section with obvious collapse of the wellbore, a correction curve volume model is first established, a function response relationship between the target curve and the reference curve is established according to the component, fluid volume content and the corresponding logging response value, and the logging response volume percentage model is used for correction, and the formula is as follows:

[0087] In the formula, Ci is the logging response value of the i-th rock component in the un-collapsed well section, Vi is the volume percentage of the i-th rock component, and C is the forward logging response value.

[0088] The multi-well consistency processing in S1 selects the trend surface analysis method for multi-well correction, the mean response of the standard layer of several key wells is selected as the logging response standard at the well point, a “trend surface” is established, the mean response of the standard layer of the multi-well is made the same as the trend surface, and the abnormality between the wells caused by the instrument system error and the like is eliminated, and the formula is as follows:

[0089] where Z(x, y) is the seismic attribute to be studied, (x, y) is the coordinate of a point on the plane, is the trend value, and g is the correction value.

[0090] Regression analysis is performed on the known Z(x, y) in the logging data to obtain a regression equation f(x, y), and the trend value is calculated according to f(x, y) Z(x, y) is subtracted from to obtain g.

[0091] In S1, the prestack random noise suppression is performed by using the GeoEast system for 3D prestack random noise attenuation. The principle is that 3D prestack seismic data can be transformed into the frequency-space domain by Fourier transform. For each frequency component S(i, j, k), a prediction operator U(l, j, k) can be designed by using the least square principle, and the prediction error energy Q(f) is expressed as the following formula:

[0092]

[0093] is the conjugate of V(i, j, k), where L, M, and N are the lengths of the 3D prestack seismic data space window (longiudinal trace, transverse trace, and offset), respectively.

[0094] Partial derivatives are calculated for the above two formulas, and

[0095] where p = 1, 2, 3, …, LOPx; q = 1, 2, 3, …, LOPy; r = 1, 2, 3, …, LOPz. LOPx, LOPy, and LOPz are the lengths of the prediction operator in the three spatial directions, respectively, and they are usually selected to have the same length. Finally, a matrix equation is obtained:

[0096] H · U = R

[0097] where H and R are the covariance matrices from the 3D prestack seismic data; and U is the 3D prediction operator.

[0098] For each frequency component, the prediction operator is obtained by solving the matrix equation, and the prediction operator is used to perform prediction filtering on the seismic data. Finally, Fourier inverse transform is used to transform back to the time-space domain to obtain the denoised result.

[0099] In S1, the non-rigid matching method is used for time difference correction of the prestack seismic data, and the implementation process is as follows:

[0100] First, the pre-stack seismic data is sorted according to the line number, the far offset dynamic correction distortion zone is cut off, the multiple wave is removed after stacking, the stacked data is taken as the reference data of NRM, and the amplitude in the reference data is taken as a three-dimensional pixel, which is denoted as Sref(x, y, z, T) (the variables x, y and z represent the spatial coordinates, and T represents the time shift between the two collected data).

[0101] Then, the displacement field is edited to delete abnormal high-frequency changes. The multi-solution gradient technology is applied in the estimation of the displacement field, and the displacement field calculates the position difference of the pixels between the two data bodies, so that the pixel positions between the two data bodies are basically consistent.

[0102] Finally, all the sample points of the pre-stack gather data are corrected by using the optimized displacement field. The sample point values before and after the correction satisfy Sm(x, y, z, T) = S(x + dx(x, y, z, T) + dz(x, y, z, T), T). In the formula, S(x, y, z, T) and Sm(x, y, z, T) represent the data before and after the correction, respectively.

[0103] After the pre-stack gather optimization processing, the consistency of the seismic data is improved, the amplitude attribute abnormal boundary is clear, the amplitude attribute detail feature change is rich, and the overall quality is improved.

[0104] Through Figure 2 It can be seen that after the pre-stack gather optimization processing, the problems such as low signal-to-noise ratio of the original gather, uneven phase axis and inconsistent spectrum of near and far gathers are effectively solved, which lays a good foundation for providing high-quality data for pre-stack inversion and the like, and the thin reservoir prediction result based on the optimized gather is more consistent with the well information, which also verifies the rationality of the gather optimization processing.

[0105] Specifically, in S4, first, d = H i / H j is calculated, H i is the total thickness of the sandstone, H j is the formation thickness, and d is the sandstone-to-formation ratio. The sandstone thickness and the formation thickness of each well are extracted. Then, a contour mapping linear interpolation method based on a triangular mesh is used to draw a plane sandstone-to-formation ratio map through five steps of boundary search, triangular subdivision, linear interpolation, contour algorithm and Bezier curve smoothing.

[0106] Boundary search: the boundary of the data point can be manually read or identified by an algorithm. For large amount of data, the convex hull scanning algorithm invented by Professor Graham is used to improve the efficiency.

[0107] Triangular subdivision: the triangular mesh subdivision is based on the Delaunay criterion, and the data nodes are meshed under the control of the above boundary, and the triangular mesh vertices after the subdivision are the positions of the data nodes.

[0108] Linear interpolation: increase interpolation points in the original grid by linear interpolation method, and three new triangular units can be obtained by connecting the three vertices of the interpolation points with the original grid.

[0109] Contour algorithm: through the third step, the small triangular grid has been formed, for each contour, search all triangular grids to find the corresponding contour points, and connect all contour points to obtain a contour. Repeat the operation to obtain all contours. Further obtain the position relationship between the contour and the triangular unit.

[0110] Bezier curve smoothing: the contour obtained by the above steps is composed of broken lines, and the contour needs to be smoothed. The quadratic Bezier curve smoothing is defined as: B(t) = (1-t) 2 P0+2t(1-t)P1+t 2 P2, wherein t∈[0,1]. Through the above steps, the sand ratio contour plan of the area is obtained. Provide the trend for the spatial lateral distribution of the next step high-resolution reservoir inversion.

[0111] In summary, the embodiment establishes variable gamma prestack deterministic inversion based on large data waveform driving low-frequency model for thin reservoir, carries out phase-controlled prestack geostatistics inversion combination, accurately depicts the internal details of thin sand layer, solves the problem of fine depiction and identification of thin layer reservoir, and combines amplitude attribute, oil and gas detection and other methods to verify each other qualitatively and quantitatively, solves the geological problems of increasing production and speed, reduces the exploration risk, has universal applicability and broad application prospect.

[0112] The above only describes the preferred embodiments of the present application and is not used to limit the present application. For those skilled in the art, the present application can have various modifications and changes. Any modification, equivalent replacement, improvement, etc. within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method for thin reservoir prediction based on fractional phasing, characterized in that, It comprises the following steps: S1, quality control and consistency processing of logging data; quality control and optimization processing of pre-stack seismic data; S2, establishing a low-frequency model according to the data processed in step S1; S3, based on the low-frequency model, combined with the data processed in S1, carrying out improved variable gamma pre-stack joint deterministic inversion to obtain a low-frequency and medium-frequency model, the formula is as follows: where t is the time sample point, t0 is the initial time, τ = τ0(s-1), s is a parameter of the degree of inhomogeneity of the medium, V P is the longitudinal wave velocity, V s is the transverse wave velocity, V ρ is the density, R p is the longitudinal wave velocity reflection coefficient, R s is the transverse wave velocity reflection coefficient, R ρ is the density reflection coefficient; R p , R s , R ρ , by the following formula: R pp (θ) = (1 + tan 2 θ)R p - 8γsin 2 θR s + (1 - 4γsin 2 θ)R ρ where θ is the angle of incidence, γ is the ratio of P- to S-wave velocities, and R pp (θ) is the P-wave reflection coefficient as a function of angle, θ, γ, and R pp (θ) are obtained from the S1 processed data; S4, according to the inversion result, carrying out sand body thickness description, and making a plane sand ratio map; S5, according to the inversion result in S3, taking the plane sand ratio map as a constraint to carry out geostatistical inversion, and expanding the inversion result to high frequency to obtain a full-band geological model.

2. The thin reservoir prediction method based on frequency division phasing according to claim 1, characterized in that: The specific steps of S2 are as follows: first, carrying out big data analysis on the seismic waveform in the result of S1 to obtain low-frequency data; then, according to the P-wave impedance, P-S wave velocity ratio and density at the well point in the result of S1, describing the physical properties and structural features of the underground medium, limiting the low-frequency data to obtain a low-frequency model.

3. The thin reservoir prediction method based on frequency division phasing according to claim 2, characterized in that: After the low-frequency model in S2 is completed, it is analyzed whether the low-frequency model meets the geological understanding, if yes, S3 is carried out, if not, S1 is returned; the judgment of whether the low-frequency model meets the geological understanding includes verification according to the result of S1 and checking the consistency of each part of the model.

4. The thin reservoir prediction method based on frequency division phasing according to claim 2, characterized in that: When establishing the low-frequency model in S2, unknown seismic attributes are calculated according to known seismic attributes, the same type of seismic attributes of different wells are obtained through a spatial conversion formula, the formula is as follows: where S n and S i are seismic attributes of different wells, S i is obtained from the results of S1; W is a cross-well relationship function, W is inferred from the results of S1; and n is the number of wells.

5. The thin reservoir prediction method based on frequency-division phased control according to claim 2, characterized in that: When establishing the low-frequency model in S2, unknown seismic attributes are calculated according to known seismic attributes, different types of seismic physical properties of the same well are obtained through a physical property conversion formula, the formula is as follows: where L j and L i are different types of seismic petrophysical properties, L i is obtained from S1 results; λ ij is a property-dependent matching function, λ ij is inferred from S1 results; n is the number of wells.

6. The thin reservoir prediction method based on frequency-division phasing according to claim 2, characterized in that: When establishing the low-frequency model in S2, unknown seismic attributes are calculated according to known seismic attributes, different types of seismic physical properties of different wells are obtained through a full-area physical property parameter conversion formula, the formula is as follows: where SL j and SL i are different types of seismic physical properties for different wells, w ij is a well-to-well relationship function for SL j and SL i corresponding to w ij is obtained by inference from S1results; λ ij is a property correlation matching function for λ ij obtained by inference from S1results; n is the number of wells.

7. The thin reservoir prediction method based on frequency-division phased control according to claim 1, characterized in that: In S5, the geostatistical inversion adopts a variogram, the formula of the variogram is as follows: where h is a given step, N(h) is the number of pairs of points at a distance equal to h, Z(x i ) is the measured value of the variable at the point x i , Z(x i +h) is the measured value of the variable at a distance h from the point x i , and r(h) is the inversion result.

8. The thin reservoir prediction method based on frequency-division phased control according to claim 1, characterized in that: The processing of logging data in S1 includes logging curve error correction and multi-well consistency processing.

9. The thin reservoir prediction method based on frequency-division phased control according to claim 8, characterized in that: The steps of logging curve error correction are as follows: For the layer section with deformation of the wellbore, Castagna formula or Smiths formula is used for correction; Castagna formula: where p is the formation density; Vp is the P-wave velocity; a, b, and c are constants obtained from the analysis of wellbore logging data of the non-deformed section. Smith's formula: S(v) = dR e where S(v) is the acoustic wave velocity; R is the resistivity; d and e are constants obtained from the analysis of the wellbore logging data of the non-deformed section. For the wellbore with the collapsed section, the volume percentage model of logging response is used for correction, and the formula is: In the formula, Ci is the logging response value of the i th rock component in the non-collapsed section; Vi is the volume percentage of the i th rock component; and C is the forward logging response value.

10. The thin reservoir prediction method based on frequency-division phasing according to claim 8, characterized in that: The multi-well consistency processing formula is: In the formula, Z(x, y) is a seismic attribute under study, (x, y) is a coordinate of a point on a plane, is a trend value, and g is a correction value. Regression analysis is performed on the known Z(x, y) in the logging data to obtain a regression equation f(x, y), and a trend value is obtained according to f(x, y) Subtracting Z(x, y) from g is obtained.

11. The thin reservoir prediction method based on frequency-division phasing according to claim 1, characterized in that: The processing of pre-stack seismic data in S1 includes pre-stack random noise suppression and time difference correction.

12. The thin reservoir prediction method based on frequency-division phasing according to claim 11, characterized in that: The pre-stack random noise suppression adopts 3D pre-stack random noise attenuation of the GeoEast system.

13. The thin reservoir prediction method based on frequency-division phasing according to claim 11, characterized in that: In S1, the non-rigid matching method is used for time difference correction of pre-stack seismic data.

14. The thin reservoir prediction method based on frequency-division phasing according to claim 1, characterized in that: The plane sand ratio map is drawn by using a gridding interpolation method; first, the plane sand ratio map is divided into a series of small grids, and the sand content value of each sample point is distributed to the corresponding grid node; then, according to the spatial position and sand content value of the known data points, the sand content value of each grid node is calculated by using the interpolation method to draw the contour map of the sand content distribution.

15. The thin reservoir prediction method based on frequency-division phasing according to claim 14, characterized in that: The gridding interpolation method is a contour mapping linear interpolation method based on a triangular grid.

Images