A structural facies constrained high resolution inversion method for river channel sand

Abstract

The application discloses a river sand facies controlled high-resolution inversion method of structural class attribute constraint, and the inversion method comprises the following steps: step S1, river facies facies belt division based on structural class attribute; step S2, well logging data preprocessing, and calculating a matching factor at a well site; step S3, calculating an inversion wavelet according to the matching factor, and performing space variation processing on the matching factor; step S4, establishing an initial impedance model based on the characteristics of a work area and structural class attribute; and step S5, performing facies controlled high-resolution inversion according to the initial impedance model. For the river sand target, the identification capability of a seismic trace beside a well is improved through well-to-seismic correlation and matching, and then high-resolution wave impedance data volume is obtained through the constraint high-resolution inversion, and a horizon slice exhibits good river sand delineation capability.

Description

Technical Field

[0001] This invention relates to inversion methods, seismic attribute extraction, well-seismic calibration, and belongs to the field of seismic inversion, particularly to a high-resolution inversion method for channel sand facies controlled by structural attribute constraints. Background Technology

[0002] Rivers can be classified according to different erosion stages in erosion cycles, sediment transport methods, river types based on deposition rates, and river channel morphology. Each classification scheme has its own characteristics but also its limitations, making direct comparisons between different disciplines impossible. Based on river morphology and sediment characteristics, rivers are classified into braided rivers, meandering rivers, distributary rivers, network rivers, and direct-flow rivers.

[0003] Currently, seismic data remains the primary means of directly obtaining information about channel sand bodies. Luo Hongmei et al. (2015) and Zhang Feifei et al. (2016) have achieved good results in describing and predicting thin reservoirs such as channel sand bodies. In 2018, Zhu Chao et al. used a frequency-division nonlinear inversion method with low-frequency model constraints to predict the location of sand bodies. In 2019, Ni Xianglong et al. applied AVF inversion technology to the Qaidam Basin, reducing the ambiguity of the inversion. While these methods have achieved good results to a certain extent, they have all failed to deeply extract the inherent channel information contained in seismic data, and the inverted data volumes all suffer from limitations such as unclear channel sand boundary identification.

[0004] Chinese patent application CN104502969A discloses a method for identifying channel sandstone reservoirs. This method includes: obtaining the seismic response characteristics of channel sand bodies through forward modeling of the reservoir, and determining seismic identification markers for channel sand bodies in different geological bodies; extracting corresponding seismic attributes based on the seismic identification markers to determine channel boundaries; determining the depositional stages of the channel sand using stratigraphic slicing technology within an isochronous stratigraphic framework; and determining the thickness of the channel sand body based on the statistical fitting relationship between the thickness of the sand body in actual wells and seismic attributes. Based on this method, in-depth seismic attribute extraction and analysis are conducted. Utilizing known well reservoir and oil-gas layer calibration results, well statistical analysis, and combined with model forward modeling results, a channel sand body identification and description technology is established. This method can accurately identify and describe ultra-deep, thin-layered, and narrow-channel channel sandstone reservoirs, and is less affected by seismic resolution.

[0005] Chinese patent application CN109541685A discloses a method for identifying channel sand bodies, belonging to the field of oilfield reservoir prediction. The method includes: seismic acquisition of a target area containing channel sand bodies to obtain seismic data; well logging within the target area to obtain well logging data; interpretation of seismic horizons using the seismic data; and interpretation of sand body thickness, sedimentary facies, and various development horizons using the well logging data; well-seismic calibration using the seismic horizons, sand body thickness, sedimentary facies, and various development horizons to obtain well-seismic calibration results; establishment of spatial horizons for sandstone groups, sub-layers, and sedimentary units in the time domain based on the well-seismic calibration results; acquisition of sedimentary facies types; and vertical combination of the spatial horizons of the above three types based on the sedimentary facies types to obtain a framework model; and identification of the distribution characteristics of channel sand bodies within the target area using the framework model. This method can identify the distribution characteristics of channel sand bodies of different thicknesses, improving the accuracy of channel sand body distribution characteristic identification.

[0006] Chinese patent application CN107656312A discloses a method and apparatus for predicting channel sand bodies based on azimuth-based stacking. The method includes: determining the main source direction of a target layer within the work area, and the fracture system characteristics of the target layer and its upper layers; obtaining a regularized migration gather based on the OVT domain for the target layer; dividing the regularized migration gather into multiple azimuth-stacked gathers according to the main source direction and the fracture system characteristics, and determining the azimuth-stacked gathers perpendicular to the main source direction; performing wave impedance inversion on the azimuth-stacked gathers perpendicular to the main source direction to obtain an inversion profile; and extracting the planar map from the inversion profile to obtain the planar distribution of channel sand bodies within the target layer. The embodiments of this application can improve the identification accuracy of channel sand body prediction. Summary of the Invention

[0007] In view of the above problems, the present invention is proposed to provide a high-resolution inversion method for channel sand facies control with structural class attribute constraints to overcome or at least partially solve the above problems.

[0008] According to one aspect of the present invention, a high-resolution inversion method for channel sand facies controlled by structural class attribute constraints is provided, the inversion method comprising:

[0009] Step S1: River facies zone division based on structural class attributes;

[0010] Step S2: Preprocess the logging data and calculate the matching factor at the well location;

[0011] Step S3: Calculate the inversion wavelet based on the matching factor, and perform space-variant processing on the matching factor;

[0012] Step S4: Establish an initial impedance model based on the characteristics of the work area and structural attributes;

[0013] Step S5: Perform phase-controlled high-resolution inversion based on the initial impedance model.

[0014] Optionally, step S1: river facies zone division based on structural class attributes specifically includes:

[0015] By using the spatial distribution of the phase axis of seismic data with wave crests and troughs as the center, the extreme points and zero-crossing points of amplitude within the time window are extracted;

[0016] Different mathematical calculation methods are used to obtain the waveform structure class attributes contained in three-dimensional data;

[0017] Based on mathematical calculation methods, waveform structure attributes are classified into integral, statistical, and difference categories.

[0018] Based on the characteristics of fluvial facies depositional models, structural data is used to divide fluvial facies zones and obtain information on facies-controlled constraints.

[0019] Optionally, the integral class includes waveform area and waveform length; the statistical class includes coefficient of variation, kurtosis, and skewness; and the difference class includes composite envelope difference, half-time curvature difference, and peak-valley kurtosis difference.

[0020] Optionally, step S2: well logging data preprocessing and calculating the matching factor at the well location specifically includes: preprocessing the well logging data, including depth correction, environmental correction and standardization of the well logging data;

[0021] Extract seismic traces near the well based on well logging data;

[0022] The matching factor at the well location is calculated using the step-by-step iterative method.

[0023] Optionally, step S3: calculating the inversion wavelet based on the matching factor and performing space-variant processing on the matching factor specifically includes:

[0024] The inversion wavelet is obtained by iteratively combining the matching factor obtained at the well location with seismic data near the well.

[0025] The spatial variation processing of matching factors is performed by the correlation coefficient constrained inverse distance interpolation method to obtain the three-dimensional matching factor data volume.

[0026] Optionally, the step of using correlation coefficient-constrained inverse distance interpolation to perform spatial variation processing on the matching factors to obtain the three-dimensional matching factor data volume specifically includes:

[0027] Using the extracted well location matching factors, interpolation is performed in the entire spatial domain to obtain the matching factor data volume of the entire three-dimensional work area.

[0028] The matching factor is spatially variable using inverse distance interpolation, and the calculation formula is as follows:

[0029]

[0030] Where j represents the target seismic trace, e[j] is the matching factor of the target seismic trace after spatial variation processing, oper[i][j] represents the matching factor between the target seismic trace j and each point of the current seismic trace, and r 2 It is the distance weight between seismic trace j and the current seismic trace, and np represents the maximum number of seismic traces.

[0031] Optionally, step S4: establishing an initial impedance model based on the characteristics of the work area and structural attributes specifically includes:

[0032] Based on the characteristics of the work area, well logging and seismic weights are set. Based on the distribution range of thick sand bodies and fluvial facies zones, combined with well logging curve interpolation, an initial impedance model with geological understanding is obtained as an inversion constraint.

[0033] Based on the characteristics of the work area, corresponding weights are assigned to match well seismic data.

[0034] Optionally, step S4: establishing an initial impedance model based on the characteristics of the work area and structural attributes specifically includes:

[0035] Based on the characteristics of fluvial facies depositional models, facies zones are divided using structural attributes. The divided facies zones are used as constraints to establish thick and thin channel regions according to the target location.

[0036] Based on the inversion target, structural attributes that can reflect the characteristics of the target are selected as the initial impedance constraint information, and interpolation is performed in combination with well logging data according to the phase zone division rules;

[0037] mod(t0) = ∑well(ti+α)·[weight(well)+weight(structural attribute)]

[0038] In the formula, mod represents the impedance value of the current calculation channel, well represents the logging impedance value, weight represents the weight of the well and the weight of the structural attribute, ti represents the longitudinal time range of the well, and α represents the perturbation variable of the structural attribute in the facies region on the time position of the logging data.

[0039] The initial weight is calculated using the following formula:

[0040]

[0041] weight(well) = 2

[0042] Where b represents the relationship between structural attributes and wellpoint channel sand thickness, and the formula for calculating b is as follows:

[0043]

[0044] Wherein, the x vector represents the wellpoint channel sand thickness, arranged in ascending order, the y vector represents the structural attribute value, and n represents the statistical number of wellpoint sand body thicknesses; when b < 0.5, it indicates that the structural attribute and the wellpoint channel sand thickness have a nearly linear relationship, and the structural attribute accounts for a large proportion; when 0.5 ≤ b ≤ 1, it indicates that the structural attribute and the wellpoint channel sand thickness have a non-linear relationship, and the proportion of the structural attribute decreases; when b > 1, there is no obvious relationship between the structural attribute and the wellpoint channel sand thickness, and the structural attribute accounts for the smallest proportion.

[0045] Optionally, step S5: performing phase-controlled high-resolution inversion based on the initial impedance model specifically includes:

[0046] A reflection coefficient model is constructed based on the initial wave impedance model, and a high-resolution impedance inversion result is obtained by using broadband constrained inversion technology in combination with the inversion wavelet obtained in step S3.

[0047] Optionally, the process of obtaining high-resolution impedance inversion results by combining the inversion wavelet obtained in step S3 with broadband constrained inversion technology specifically includes:

[0048] In the inversion process, it is assumed that the result to be obtained is The constraints are set as follows: And require the inversion results Gradually moving towards constraints during the iteration process deflection;

[0049] Based on the bias requirement, the following relation is obtained:

[0050]

[0051] In the formula, G is the wavelet matrix, G T Let I be the transpose of the matrix, and let I be the identity matrix. The data is structured attribute data, where α and β are constraint factors;

[0052] By modifying the model, broadband, high-resolution information from well logging data at the wellbore location is transferred to the seismic data;

[0053] The final result of the inversion is the reflection coefficient profile, using a high-resolution seismic wavelet G. * The two are convolved to synthesize a high-resolution seismic data volume:

[0054]

[0055] Geological models constrained by structural properties are used as constraints for broadband constrained inversion to characterize channel sands.

[0056] This invention provides a high-resolution facies-controlled inversion method for channel sands constrained by structural attributes. The inversion method includes: Step S1: dividing the river facies zones based on structural attributes; Step S2: preprocessing well logging data and calculating the matching factor at the well location; Step S3: calculating the inversion wavelet based on the matching factor and performing spatial variation processing on the matching factor; Step S4: establishing an initial impedance model based on the characteristics of the work area and structural attributes; Step S5: performing facies-controlled high-resolution inversion based on the initial impedance model. For channel sand targets, well-seismic comparison matching improves the identification capability of seismic traces near the well, and then, through constrained high-resolution inversion, a high-resolution wave impedance data volume is obtained, whose stratigraphic slices demonstrate excellent channel sand characterization capabilities.

[0057] The above description is merely an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention and to implement it in accordance with the contents of the specification, and in order to make the above and other objects, features and advantages of the present invention more apparent and understandable, specific embodiments of the present invention are described below. Attached Figure Description

[0058] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0059] Figure 1 A flowchart illustrating a high-resolution inversion method for channel sand facies control with structural attribute constraints provided in an embodiment of the present invention;

[0060] Figure 2 Seismic profile data provided for embodiments of the present invention;

[0061] Figure 3 This is a profile of seismic wave peak properties provided in an embodiment of the present invention;

[0062] Figure 4 This is a well logging seismic matching result diagram provided in an embodiment of the present invention; Figure 4 (a) shows the original velocity, density, and reflection coefficient curves. Figure 4 (b) shows the velocity, density, and reflection coefficient curves after frequency reduction and matching with the earthquake;

[0063] Figure 5 This is a high-resolution inversion plan view of the channel sand facies controlled by the structural class attribute constraints of the Ngx oil layer in this embodiment of the invention;

[0064] Figure 6 This is a comparison and analysis diagram of the coefficient of variation attribute and amplitude attribute in an embodiment of the present invention;

[0065] Figure 7 This is a phase-controlled high-resolution inversion profile superimposed SP logging curve diagram with structural class attribute constraints in an embodiment of the present invention;

[0066] Figure 8 This is a comparative analysis diagram of the low-frequency model with structural attribute constraints and the conventional low-frequency model in the embodiments of the present invention. Detailed Implementation

[0067] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0068] The terms "comprising" and "having," and any variations thereof, in the specification, embodiments, claims, and drawings of this invention are intended to cover non-exclusive inclusion, such as including a series of steps or units.

[0069] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0070] Example 1

[0071] Practical application tests were conducted using data from a certain work area 1, and well logging and seismic data of the Guantao Formation riverbed sand strata were tested and processed.

[0072] like Figure 1 As shown, a high-resolution inversion method for channel sand facies controlled by structural attribute constraints includes:

[0073] Step S1: Extraction of seismic waveform structural attributes, including:

[0074] Waveform structure attributes mainly refer to a class of layer attributes extracted from seismic data. Based on different mathematical calculation methods, they can be categorized into three types: The first type is integral structure attributes that characterize waveform strength variations, such as waveform area and waveform length; the second type is statistical structure attributes that characterize waveform stability, sharpness, and symmetry, corresponding to coefficient of variation, kurtosis, and skewness, respectively. Figure 2 This refers to the seismic profile data in this embodiment. Figure 3This example shows a profile of the wave crest attribute in the seismic structural attributes category. The third category consists of differential structural attributes characterizing waveform structure changes, such as composite envelope difference, half-time curvature difference, and peak-valley kurtosis difference. Combined with fluvial sedimentary models and sand body distribution patterns, structural attributes can be used to delineate fluvial facies zones and identify fluvial sand body regions. The specific steps are as follows:

[0075] (1) Classify river sedimentary facies according to different waveforms and peak properties;

[0076] (2) Locate the distribution area of ​​river sand bodies by measuring waveform length, frequency, and phase change properties;

[0077] (3) The distribution area of ​​thicker fluvial sand bodies can be found by using waveform area, length and average curvature properties;

[0078] (4) Identify the thin sand body superposition region in the thick sand body region by using kurtosis, coefficient of variation and skewness attributes;

[0079] (5) The river channel superposition area is identified by verifying the composite envelope difference, half-time curvature difference, and peak-valley kurtosis difference attributes.

[0080] Step 2: Preprocessing of well seismic data and determination of matching factors at well locations, including:

[0081] Step 2.1: Well Seismic Data Preprocessing

[0082] Well logging and seismic data differ in many aspects, such as measurement scale and frequency band range, making direct comparison impossible. Therefore, a series of preprocessing steps are required for well logging data to match with seismic traces, thus better serving subsequent seismic interpretation. Well logging data preprocessing is a crucial step before well logging interpretation, mainly including depth correction, standardization, and environmental correction of well logging curves. Based on well location information, seismic data is extracted from the well path. Since well logging data is generally depth-domain while seismic data is generally time-domain, accurate well-seismic calibration, depth-time conversion, and resampling are necessary for the well logging data.

[0083] Step 2.2: Frequency reduction of well logging data using a step-by-step iterative method.

[0084] Based on the Zoeppritz equation, the well logging reflection coefficient sequence is calculated using sonic logging curves and density logging curves. This calculated initial well logging reflection coefficient sequence has high resolution and exhibits broadband characteristics, with most high-frequency components exceeding the frequency range of surface seismic data. Furthermore, well logging data is affected by numerous factors during excitation and acquisition, resulting in significant high-frequency noise in the original data, inevitably leading to discrepancies with the actual subsurface conditions. Directly using the initial well logging reflection coefficient sequence for frequency extension processing may result in unreliable extension results. Therefore, frequency reduction processing of the well logging data is necessary to match it with surface seismic data. From the perspective of the well logging reflection coefficient sequence, using the wellside trace as a basis, under constraints, the well logging reflection coefficient sequence and the estimated seismic wavelet are continuously modified through step-by-step iteration, gradually reducing the frequency of the well logging reflection sequence during the matching process with the seismic trace. Finally, a reflection coefficient sequence that simultaneously contains well logging and seismic information and matches the seismic trace is obtained.

[0085] Figure 4 This diagram shows the well logging reflection coefficient and the matched reflection coefficient.

[0086] Figure 4 (a) shows the original velocity, density, and reflection coefficient curves. Figure 4 (b) shows the velocity, density, and reflection coefficient curves after frequency reduction and matching with the earthquake. Figure 4 It can be seen that the vertical resolution of the matched velocity, density, and reflection coefficient curves is significantly lower than that of the original logging data.

[0087] From the sonic logging curves and density logging curves after deep-time conversion resampling, the wave propagation velocity v at each sampling point can be obtained. i and formation density ρ i The wave impedance is expressed as ρ i v i According to the Zoeppritz equation, the reflection coefficient can be expressed as:

[0088]

[0089] Within the time-depth range of the seismic trace near the well, the reflection coefficient at each sampling point is calculated using the above formula, resulting in a reflection coefficient sequence r. o (i), r o (i) is called the initial reflection coefficient sequence of the well logging. Assuming the seismic wavelet is w(t) and the actual wellside seismic trace at the well location is S(i), the synthesized seismic trace is expressed as:

[0090] S(t)=r(t)*w(t)=∑r(τ)w(t-τ)

[0091] There must be a certain error between the synthesized seismic trace S(t) and the well-side seismic trace S(i). The difference of the squared errors between the two can be expressed as:

[0092]

[0093] In the formula: N is the wavelet length, and M is the length of the ground seismic data. The smaller the error E value, the higher the matching degree between the well-side seismic trace and the synthetic seismic trace. To minimize the sum of squared errors, i.e., to take the extreme value of the above formula, we have:

[0094]

[0095] Converting the above equation into matrix form G1r=d1:

[0096]

[0097] When extracting the reflection coefficient sequence from the seismic traces near the well using the above formula, the initial reflection coefficient sequence of the well logging is used as a constraint condition:

[0098]

[0099] In the formula, I is the identity matrix. Let be the transpose matrix, and β be the constraint factor for the initial reflection coefficient sequence of the well logging. Similarly, by taking the partial derivative with respect to w(i), we obtain:

[0100]

[0101] In the formula, α is the constraint factor of the wavelet, and W(t) is the seismic wavelet extracted from the well-side trace using the autocorrelation method. By iteratively modifying the reflection coefficient sequence and the seismic wavelet multiple times until the sum of squared errors E between the convolutional seismic trace and the well-side seismic trace reaches its minimum value, the matching factor at the well location is obtained.

[0102] Step 2.3: Calculation of well-seismic matching factor at well location

[0103] Matching factors are extracted based on inter-well and well-side seismic traces. First, a pseudo-inter-well seismic pattern is constructed using the iteratively derived reflection coefficient sequence and the inter-well seismic wavelet. Second, the well-side seismic data undergoes SFMOD processing; the sparse pulse-processed seismic data can be used as reflection coefficients to match the pseudo-inter-well seismic records, thus determining the matching factor at the well location. The specific implementation process is as follows.

[0104] Using well logging data and well-side seismic traces, an optimization operator is designed. Given that the pseudo-well-to-well seismic record synthesized from the iterative reflection coefficient sequence is y(t) = {y(0), y(1), ..., y(np)}, and the corresponding well-side seismic trace data, after SFMOD processing, is x(t) = {x(0), x(1), ..., x(np)}, an optimization operator ω(t) = {ω(-l), ω(-l+1), ..., ω(0), ω(1), ..., ω(l)} is designed to continuously approach y(t) with the convolution x(t)*ω(t) while minimizing the sum of squared errors. Assuming z(t) = x(t)*ω(t), then:

[0105]

[0106] To minimize the sum of squared errors Q, ω(t) should satisfy:

[0107]

[0108] Substituting and simplifying, we get:

[0109]

[0110] The above formula can also be expressed as:

[0111]

[0112] Represented in matrix form as follows:

[0113]

[0114] In the formula, γ xx (i) represents the autocorrelation of the well-side seismic trace data after SFMOD processing, γ xy (i) represents the cross-correlation between well-side seismic traces and inter-well seismic data after SFMOD processing, and np represents the length of the matching factor.

[0115] When using a convolution model to extract matching factors, it is actually an optimization operator.

[0116] Let a(t) = {a(-l), a(-l+1), ..., a(0), a(1), ..., a(l)}, such that the convolution y(t)*a(t) continuously approaches x(t) while minimizing the sum of squared errors. Then the above equation can be expressed as:

[0117]

[0118] The matrix equation is called the Toeplitz matrix equation, and the matching factor at the well location is obtained by recursion.

[0119] Step 3: Obtaining the matching factor spatial variation and inversion wavelet

[0120] Step 3.1: Perform spatial variation of the matching factor using inverse distance interpolation.

[0121] The spatial variation technique for matching factors mainly utilizes the extracted matching factors at the well locations, interpolating them across the entire spatial domain to obtain the matching factor data volume for the entire three-dimensional work area. The inverse distance interpolation method is used for spatial variation processing of the matching factors, and the calculation formula is as follows:

[0122]

[0123] Where j represents the target seismic trace, e[j] is the matching factor of the target seismic trace after spatial variation processing, oper[i][j] represents the matching factor between the target seismic trace j and each point of the current seismic trace, and r 2 It is the distance weight between seismic trace j and the current seismic trace, and np represents the maximum number of seismic traces.

[0124] Step 3.2: Iteratively calculate the inversion wavelet

[0125] Using the reduced-frequency logging reflection coefficient and wellside seismic trace obtained in step 2.2 as input, the convolution model is used with the least squares method to take the 30Hz zero-phase Ricker wavelet as the initial wavelet. By continuously adjusting the Ricker wavelet parameters, the optimal wavelet is finally obtained for subsequent inversion calculations.

[0126] Step 4: Establishing the initial impedance model based on the characteristics of the work area and structural attributes

[0127] Based on the characteristics of fluvial sedimentary models, facies zones are divided using structural attributes. These facies zones serve as constraints, and thick and thin channel regions are established according to the target location. Structural attributes reflecting the characteristics of the inverted target are selected as initial impedance constraint information, and interpolation is performed based on well logging data according to the facies zone division rules.

[0128] mod(t0) = ∑well(ti+α)·[weight(well)+weight(structural attribute)]

[0129] In the formula, mod represents the impedance value of the current calculation channel, well represents the logging impedance value, weight represents the weight of the corresponding well and the weight of the structural attribute, ti represents the longitudinal time range of the well, and α represents the perturbation variable of the structural attribute in the phasor region on the time position of the logging data.

[0130] The initial weight is calculated using the following formula:

[0131]

[0132] weight(well) = 2

[0133] Where b represents the relationship between structural attributes and wellpoint channel sand thickness, and the formula for calculating b is as follows:

[0134]

[0135] In the above formula, the x vector represents the wellpoint channel sand thickness (arranged from smallest to largest), the y vector represents the structural attribute value, and n represents the statistical number of wellpoint sand body thicknesses. When b < 0.5, it indicates that the structural attribute and the wellpoint channel sand thickness have an approximately linear relationship, and the structural attribute accounts for a large proportion; when 0.5 ≤ b ≤ 1, it indicates that the structural attribute and the wellpoint channel sand thickness have a non-linear relationship, and the proportion of the structural attribute decreases; when b > 1, there is no significant relationship between the structural attribute and the wellpoint channel sand thickness, and the structural attribute accounts for the smallest proportion.

[0136] Step 5: Phased-array high-resolution inversion

[0137] A reflection coefficient model is constructed using an initial wave impedance model. The reflection coefficient model is then compared with the structural attribute profile extracted from seismic data. The error between the synthesized seismic structural attribute profile and the well track is then compared. The minimum error required is used as a constraint to iteratively modify the model. By modifying the model, the weights of different attributes in the structural attributes are adjusted so that the high-frequency information in the constraint conditions (well logging data) can guide and enhance the effective information in the seismic data. Various frequency components in the data move closer to the standard of the constraint conditions, ultimately achieving the goal of improving the resolution of seismic data while reducing the difficulty of channel sand characterization.

[0138] In the inversion process, it is assumed that the result to be obtained is The constraints are set as follows: And require the inversion results Gradually moving towards constraints during the iteration process Deflection. Based on the bias requirements of both, the following relationship is obtained:

[0139]

[0140] In the formula, G is the wavelet matrix, G T Let I be the transpose of the matrix, and let I be the identity matrix. The data represents structural attributes, with α and β being constraint factors.

[0141] By modifying the model, broadband, high-resolution information from well logging data at the wellbore location is gradually transferred to the seismic data. The final result of the inversion is a reflection coefficient profile, using a high-resolution seismic wavelet G. * The two are convolved to synthesize a high-resolution seismic data volume:

[0142]

[0143] Using geological models constrained by structural properties as constraints for broadband constrained inversion can enhance and guide the effective information in seismic data, thereby improving the resolution of seismic data and achieving the goal of channel sand characterization.

[0144] This invention effectively improves the identification capability of channel sands by combining seismic and well logging data, extracting and analyzing seismic structural attributes, dividing river facies zones, and constraining high-resolution inversion techniques. Well logging currently offers the highest vertical resolution data, while seismic data offers the highest horizontal resolution. Although geophysicists have been working to improve the vertical resolution of seismic data, the complexity of subsurface conditions, the ambiguity of geophysics, and noise in seismic data have made it difficult to more precisely characterize small-scale geological bodies, such as channel sands, using seismic data. Combining the ability of well logging curves to distinguish strata and identify lithology vertically with the ability to control the distribution of geological bodies in three-dimensional seismic space to better characterize the distribution of channels and channel sands in three-dimensional subsurface space remains a pressing problem for oil and gas resource exploration.

[0145] The phase-controlled high-resolution inversion module based on channel sand structure class attribute constraints combines well logging and seismic analysis. It utilizes seismic structure class attribute facies zone division and well-seismic matching information as constraints, and employs constrained high-resolution inversion to solve the aforementioned problems. For channel sand targets, well-seismic comparison and matching improves the identification capability of well-side seismic traces. Furthermore, constrained high-resolution inversion yields a high-resolution acoustic impedance data volume, whose stratigraphic slices demonstrate excellent channel sand characterization capabilities.

[0146] Figure 5 This embodiment shows a phase-controlled high-resolution inversion planar image of the Ngx oil layer, constrained by structural class attributes. Figure 5 The planar distribution of the channel sandstone in the figure is shown for two oil layers, Ngx1 oil layer 2 and Ngx1 oil layer 3. The bright color in the figure represents the channel sandstone development area.

[0147] Example 2

[0148] Practical application tests were conducted using data from a certain work area 2, and well logging and seismic data of the Guantao Formation river sand strata were tested and processed.

[0149] The specific implementation steps are as follows:

[0150] Step 1: Extraction of seismic waveform structural attributes

[0151] Waveform structure attributes mainly refer to a class of layer attributes extracted from seismic data. Based on different mathematical calculation methods, they can be categorized into three types: The first type is integral-type structure attributes that characterize waveform strength variations, such as waveform area and waveform length; the second type is statistical-type structure attributes that characterize waveform stability, sharpness, and symmetry, corresponding to coefficient of variation, kurtosis, and skewness, respectively. Figure 6 This is a comparative analysis chart of the coefficient of variation and amplitude attributes in this embodiment. Figure 6 The comparison shows that the coefficient of variation attribute more clearly characterizes fault morphology and sand body distribution; the third category is the differential structural attribute characterizing waveform structure changes, such as composite envelope difference, half-time tortuosity difference, and peak-valley kurtosis difference. Combined with fluvial sedimentary models and sand body distribution patterns, structural attributes can delineate fluvial facies zones and identify fluvial sand body regions. The specific steps are as follows:

[0152] (1) Classify river sedimentary facies according to different waveforms and peak properties;

[0153] (2) Locate the distribution area of ​​river sand bodies by measuring waveform length, frequency, and phase change properties;

[0154] (3) The distribution area of ​​thicker fluvial sand bodies can be found by using waveform area, length and average curvature properties;

[0155] (4) Identify the thin sand body superposition region in the thick sand body region by using kurtosis, coefficient of variation and skewness attributes;

[0156] (5) The river channel superposition area is identified by verifying the composite envelope difference, half-time curvature difference, and peak-valley kurtosis difference attributes.

[0157] Step 2: Preprocessing of seismic data and determination of matching factors at well locations

[0158] Step 2.1: Well Seismic Data Preprocessing

[0159] Due to differences in measurement scale and frequency band range between well logging and seismic data, direct comparison is not possible. Therefore, a series of preprocessing steps are required for well logging data to match the seismic traces, thus better serving subsequent seismic interpretation. Well logging data preprocessing is a crucial step before well logging interpretation, mainly including depth correction, standardization, and environmental correction of well logging curves. Based on well location information, seismic data is extracted from the well path. Since well logging data is generally depth-domain while seismic data is generally time-domain, accurate well-seismic calibration, depth-time conversion, and resampling are necessary.

[0160] Step 2.2: Frequency reduction of well logging data using a step-by-step iterative method.

[0161] Based on the Zoeppritz equation, the well logging reflection coefficient sequence is calculated using sonic logging curves and density logging curves. This calculated initial well logging reflection coefficient sequence has high resolution and exhibits broadband characteristics, with most high-frequency components exceeding the frequency range of surface seismic data. Furthermore, well logging data is affected by numerous factors during excitation and acquisition, resulting in significant high-frequency noise in the original data, inevitably leading to discrepancies with the actual subsurface conditions. Directly using the initial well logging reflection coefficient sequence for frequency extension processing may result in unreliable extension results. Therefore, frequency reduction processing of the well logging data is necessary to match it with surface seismic data. From the perspective of the well logging reflection coefficient sequence, using the wellside trace as a basis, under constraints, the well logging reflection coefficient sequence and the estimated seismic wavelet are continuously modified through step-by-step iteration, gradually reducing the frequency of the well logging reflection sequence during the matching process with the seismic trace. Finally, a reflection coefficient sequence that simultaneously contains well logging and seismic information and matches the seismic trace is obtained.

[0162] From the sonic logging curves and density logging curves after deep-time conversion resampling, the wave propagation velocity v at each sampling point can be obtained. i and formation density ρ i The wave impedance is expressed as ρ i v i According to the Zoeppritz equation, the reflection coefficient can be expressed as:

[0163]

[0164] Within the time-depth range of the seismic trace near the well, the reflection coefficient at each sampling point is calculated using the above formula, resulting in a reflection coefficient sequence r. o (i), r o (i) is called the initial reflection coefficient sequence of the well logging. Assuming the seismic wavelet is w(t) and the actual wellside seismic trace at the well location is S(i), the synthesized seismic trace is expressed as:

[0165] S(t)=r(t)*w(t)=∑r(τ)w(t-τ)

[0166] There must be a certain error between the synthesized seismic trace S(t) and the well-side seismic trace S(i). The difference of the squared errors between the two can be expressed as:

[0167]

[0168] In the formula: N is the wavelet length, and M is the length of the ground seismic data. The smaller the error E value, the higher the matching degree between the well-side seismic trace and the synthetic seismic trace. To minimize the sum of squared errors, i.e., to take the extreme value of the above formula, we have:

[0169]

[0170] Converting the above equation into matrix form G1r=d1:

[0171]

[0172] When extracting the reflection coefficient sequence from the seismic traces near the well using the above formula, the initial reflection coefficient sequence of the well logging is used as a constraint condition:

[0173]

[0174] In the formula, I is the identity matrix. Let be the transpose matrix, and β be the constraint factor for the initial reflection coefficient sequence of the well logging.

[0175] Similarly, by taking the partial derivative with respect to w(i), we can derive:

[0176]

[0177] In the formula, α is the constraint factor of the wavelet, and W(t) is the seismic wavelet extracted from the well bypass using the autocorrelation method.

[0178] By iteratively modifying the reflection coefficient sequence and seismic wavelet multiple times until the sum of squared errors E between the two convolved seismic traces and the well-side seismic traces reaches its minimum value, the matching factor at the well location is obtained.

[0179] Step 2.3: Calculation of the well-seismic matching factor at the well location, including:

[0180] Matching factors are extracted based on inter-well seismic and well-side seismic traces. First, pseudo-inter-well seismic data is constructed using the iterative reflection coefficient sequence and inter-well seismic wavelet. Second, the well-side seismic data is processed by SFMOD, and the seismic data after sparse pulse processing is used as the reflection coefficient to match the pseudo-inter-well seismic records and obtain the matching factor at the well location.

[0181] The specific implementation process is as follows:

[0182] Using well logging data and wellside seismic trace design optimization operators, it is known that the pseudo-well seismic record synthesized using the iterated reflection coefficient sequence is y(t)={y(0),y(1),…,y(np)}, and the corresponding wellside seismic trace data after SFMOD processing is x(t)={x(0),x(1),…,x(np)};

[0183] Design an optimal operator ω(t) = {ω(-l), ω(-l+1), ..., ω(0), ω(1), ..., ω(l)} such that the convolution x(t)*ω(t) continuously approaches y(t) while minimizing the sum of squared errors.

[0184] Assume z(t) = x(t) * ω(t), then:

[0185]

[0186] To minimize the sum of squared errors Q, ω(t) should satisfy:

[0187]

[0188] Substituting and simplifying, we get:

[0189]

[0190] The above formula can also be expressed as:

[0191]

[0192] Represented in matrix form as follows:

[0193]

[0194] In the formula, γ xx (i) represents the autocorrelation of the well-side seismic trace data after SFMOD processing, γ xy (i) represents the cross-correlation between well-side seismic traces and inter-well seismic data after SFMOD processing, and np is the length of the matching factor. When using the convolution model to extract the matching factor, an optimization operator a(t) = {a(-l), a(-l+1), ..., a(0), a(1), ..., a(l)} is actually used to make the convolution y(t)*a(t) continuously approach x(t) while minimizing the sum of squared errors.

[0195] The above expression can be represented as:

[0196]

[0197] The matrix equation is called the Toeplitz matrix equation, and the matching factor at the well location is obtained by recursion.

[0198] Step 3: Calculation of the matching factor spatial variation and inversion wavelet, including:

[0199] Step 3.1: Perform spatial variation of the matching factor using inverse distance interpolation.

[0200] The spatial variation technique for matching factors mainly utilizes the extracted matching factors at the well locations, interpolating them across the entire spatial domain to obtain the matching factor data volume for the entire three-dimensional work area. The inverse distance interpolation method is used for spatial variation processing of the matching factors, and the calculation formula is as follows:

[0201]

[0202] Where j represents the target seismic trace, e[j] is the matching factor of the target seismic trace after spatial variation processing, oper[i][j] represents the matching factor between the target seismic trace j and each point of the current seismic trace, and r 2 It is the distance weight between seismic trace j and the current seismic trace, and np represents the maximum number of seismic traces.

[0203] Step 3.2: Iteratively calculate the inversion wavelet

[0204] Using the reduced-frequency logging reflection coefficient and wellside seismic trace obtained in step 2.2 as input, the convolution model is used with the least squares method to take the 30Hz zero-phase Ricker wavelet as the initial wavelet. By continuously adjusting the Ricker wavelet parameters, the optimal wavelet is finally obtained for subsequent inversion calculations.

[0205] Step 4: Establishing the initial impedance model based on the characteristics of the work area and structural attributes

[0206] Based on the characteristics of fluvial sedimentary models, facies zones are divided using structural attributes. These facies zones serve as constraints, and thick and thin channel regions are established according to the target location. Structural attributes reflecting the characteristics of the inverted target are selected as initial impedance constraint information, and interpolation is performed based on well logging data according to the facies zone division rules.

[0207] mod(t0) = ∑well(ti+α)·[weight(well)+weight(structural attribute)]

[0208] In the formula, mod represents the impedance value of the current calculation channel, well represents the logging impedance value, weight represents the weight of the well and the weight of the structural attribute, ti represents the longitudinal time range of the well, and α represents the perturbation variable of the structural attribute in the facies zone on the time position of the logging data.

[0209] The initial weight is calculated using the following formula:

[0210]

[0211] weight(well) = 2

[0212] Where b represents the relationship between structural attributes and wellpoint channel sand thickness, and the formula for calculating b is as follows:

[0213]

[0214] In the above formula, the x vector represents the wellpoint channel sand thickness (arranged from smallest to largest), the y vector represents the structural attribute value, and n represents the statistical number of wellpoint sand body thicknesses. When b < 0.5, it indicates that the structural attribute and the wellpoint channel sand thickness have an approximately linear relationship, and the structural attribute accounts for a large proportion; when 0.5 ≤ b ≤ 1, it indicates that the structural attribute and the wellpoint channel sand thickness have a non-linear relationship, and the proportion of the structural attribute decreases; when b > 1, there is no significant relationship between the structural attribute and the wellpoint channel sand thickness, and the structural attribute accounts for the smallest proportion.

[0215] Step 5: Phased-array high-resolution inversion

[0216] A reflection coefficient model is constructed using an initial wave impedance model. The reflection coefficient model is then compared with the structural attribute profile extracted from seismic data. The error between the synthesized seismic structural attribute profile and the well track is then compared. The minimum error required is used as a constraint to iteratively modify the model. By modifying the model, the weights of different attributes in the structural attributes are adjusted so that the high-frequency information in the constraint conditions (well logging data) can guide and enhance the effective information in the seismic data. Various frequency components in the data move closer to the standard of the constraint conditions, ultimately achieving the goal of improving the resolution of seismic data while reducing the difficulty of channel sand characterization.

[0217] In the inversion process, it is assumed that the result to be obtained is The constraints are set as follows: And require the inversion results Gradually moving towards constraints during the iteration process Deflection. Based on the bias requirements of both, the following relationship is obtained:

[0218]

[0219] In the formula, G is the wavelet matrix, G T Let I be the transpose of the matrix, and let I be the identity matrix. The data represents structural attributes, with α and β being constraint factors.

[0220] By modifying the model, broadband, high-resolution information from well logging data at the wellbore location is gradually transferred to the seismic data. The final result of the inversion is a reflection coefficient profile, using a high-resolution seismic wavelet G. * The two are convolved to synthesize a high-resolution seismic data volume:

[0221]

[0222] Using geological models constrained by structural properties as constraints for broadband constrained inversion can enhance and guide the effective information in seismic data, thereby improving the resolution of seismic data and achieving the goal of channel sand characterization.

[0223] This invention effectively improves the identification capability of channel sands by combining seismic and well logging data, extracting and analyzing seismic structural attributes, dividing river facies zones, and constraining high-resolution inversion techniques. Well logging currently offers the highest vertical resolution data, while seismic data offers the highest horizontal resolution. Although geophysicists have been working to improve the vertical resolution of seismic data, the complexity of subsurface conditions, the ambiguity of geophysics, and noise in seismic data have made it difficult to more precisely characterize small-scale geological bodies, such as channel sands, using seismic data. Combining the ability of well logging curves to distinguish strata and identify lithology vertically with the ability to control the distribution of geological bodies in three-dimensional seismic space to better characterize the distribution of channels and channel sands in three-dimensional subsurface space remains a pressing problem for oil and gas resource exploration.

[0224] The phase-controlled high-resolution inversion module based on channel sand structure class attribute constraints combines well logging and seismic analysis. It utilizes seismic structure class attribute facies zone division and well-seismic matching information as constraints, and employs constrained high-resolution inversion to solve the aforementioned problems. For channel sand targets, well-seismic comparison and matching improves the identification capability of well-side seismic traces. Furthermore, constrained high-resolution inversion yields a high-resolution acoustic impedance data volume, whose stratigraphic slices demonstrate excellent channel sand characterization capabilities. Figure 7 The image shown is a phase-controlled high-resolution inversion profile overlaid with SP logging curves, constrained by structural class attributes, in this embodiment. Figure 7 The high-resolution inversion results of the medium-structure attribute-constrained phase control clearly depict the spatial distribution morphology of the sand body and have a high degree of agreement with the SP logging curve.

[0225] Example 3

[0226] Practical application tests were conducted using data from a certain work area 3, and well logging and seismic data of the Guantao Formation river sand strata were tested and processed.

[0227] The specific steps include:

[0228] Step 1: Extraction of seismic waveform structural attributes

[0229] Waveform structure attributes mainly refer to a class of layer attributes extracted from seismic data. Based on different mathematical calculation methods, they can be categorized into three types: the first type is integral structure attributes characterizing waveform strength variations, such as waveform area and waveform length; the second type is statistical structure attributes characterizing waveform stability, sharpness, and symmetry, corresponding to coefficient of variation, kurtosis, and skewness, respectively. Combined with fluvial sedimentary models and sand body distribution patterns, structure attributes can be used to delineate fluvial facies zones and identify fluvial sand body regions. The specific steps are as follows:

[0230] (1) Classify river sedimentary facies according to different waveforms and peak properties;

[0231] (2) Locate the distribution area of ​​river sand bodies by measuring waveform length, frequency, and phase change properties;

[0232] (3) The distribution area of ​​thicker fluvial sand bodies can be found by using waveform area, length and average curvature properties;

[0233] (4) Identify the thin sand body superposition region in the thick sand body region by using kurtosis, coefficient of variation and skewness attributes;

[0234] (5) The river channel superposition area is identified by verifying the composite envelope difference, half-time curvature difference, and peak-valley kurtosis difference attributes.

[0235] Step 2: Preprocessing of well seismic data and determination of matching factors at well locations, including:

[0236] Step 2.1: Well Seismic Data Preprocessing

[0237] Due to differences in measurement scale and frequency band range between well logging and seismic data, direct comparison is not possible. Therefore, a series of preprocessing steps are required for well logging data to match the seismic traces, thus better serving subsequent seismic interpretation. Well logging data preprocessing is a crucial step before well logging interpretation, mainly including depth correction, standardization, and environmental correction of well logging curves. Based on well location information, seismic data is extracted from the well path. Since well logging data is generally depth-domain while seismic data is generally time-domain, accurate well-seismic calibration, depth-time conversion, and resampling are necessary.

[0238] Step 2.2: Frequency reduction of well logging data using a step-by-step iterative method.

[0239] Based on the Zoeppritz equation, the well logging reflection coefficient sequence is calculated using sonic logging curves and density logging curves. This calculated initial well logging reflection coefficient sequence has high resolution and exhibits broadband characteristics, with most high-frequency components exceeding the frequency range of surface seismic data. Furthermore, well logging data is affected by numerous factors during excitation and acquisition, resulting in significant high-frequency noise in the original data, inevitably leading to discrepancies with the actual subsurface conditions. Directly using the initial well logging reflection coefficient sequence for frequency extension processing may result in unreliable extension results. Therefore, frequency reduction processing of the well logging data is necessary to match it with surface seismic data. From the perspective of the well logging reflection coefficient sequence, using the wellside trace as a basis, under constraints, the well logging reflection coefficient sequence and the estimated seismic wavelet are continuously modified through step-by-step iteration, gradually reducing the frequency of the well logging reflection sequence during the matching process with the seismic trace. Finally, a reflection coefficient sequence that simultaneously contains well logging and seismic information and matches the seismic trace is obtained.

[0240] From the sonic logging curves and density logging curves after deep-time conversion resampling, the wave propagation velocity v at each sampling point can be determined. i and formation density ρ i The wave impedance is expressed as ρ i v i According to the Zoeppritz equation, the reflection coefficient can be expressed as:

[0241]

[0242] Within the time-depth range of the seismic trace near the well, the reflection coefficient at each sampling point is calculated using the above formula, resulting in a reflection coefficient sequence r. o (i), r o (i) is called the initial reflection coefficient sequence of the well logging. Assuming the seismic wavelet is w(t) and the actual wellside seismic trace at the well location is S(i), the synthesized seismic trace is expressed as:

[0243] S(t)=r(t)*w(t)=∑r(τ)w(t-τ)

[0244] There must be a certain error between the synthesized seismic trace S(t) and the well-side seismic trace S(i). The difference of the squared errors between the two can be expressed as:

[0245]

[0246] In the formula: N is the wavelet length, and M is the length of the ground seismic data. The smaller the error E value, the higher the matching degree between the well-side seismic trace and the synthetic seismic trace. To minimize the sum of squared errors, i.e., to take the extreme value of the above formula, we have:

[0247]

[0248] Converting the above equation into matrix form G1r=d1:

[0249]

[0250] When extracting the reflection coefficient sequence from the seismic traces near the well using the above formula, the initial reflection coefficient sequence of the well logging is used as a constraint condition:

[0251]

[0252] In the formula, I is the identity matrix. Let be the transpose matrix, and β be the constraint factor for the initial reflection coefficient sequence of the well logging. Similarly, by taking the partial derivative with respect to w(i), we obtain:

[0253]

[0254] In the formula, α is the constraint factor of the wavelet, and W(t) is the seismic wavelet extracted from the well-side trace using the autocorrelation method. By iteratively modifying the reflection coefficient sequence and the seismic wavelet multiple times until the sum of squared errors E between the convolutional seismic trace and the well-side seismic trace reaches its minimum value, the matching factor at the well location can be obtained.

[0255] Step 2.3: Calculation of well-seismic matching factor at well location

[0256] Matching factors are extracted based on inter-well and well-side seismic traces. First, a pseudo-inter-well seismic pattern is constructed using the iteratively derived reflection coefficient sequence and the inter-well seismic wavelet. Second, the well-side seismic data undergoes SFMOD processing; the sparse pulse-processed seismic data can be used as reflection coefficients to match the pseudo-inter-well seismic records, thus determining the matching factor at the well location. The specific implementation process is as follows.

[0257] Using well logging data and well-side seismic traces, an optimization operator is designed. It is known that the pseudo-well seismic record synthesized from the iterative reflection coefficient sequence is y(t)={y(0), y(1), …, y(np)}, and the corresponding well-side seismic trace data after SFMOD processing is x(t)={x(0), x(1), …, x(np)}. An optimal operator ω(t)={ω(-l), ω(-l+1), …, ω(0), ω(1), …, ω(l)} is designed so that the convolution x(t)*ω(t) continuously approaches y(t) while minimizing the sum of squared errors.

[0258] Assume z(t) = x(t) * ω(t), then:

[0259]

[0260] To minimize the sum of squared errors Q, ω(t) should satisfy:

[0261]

[0262] Substituting and simplifying, we get:

[0263]

[0264] The above formula can also be expressed as:

[0265]

[0266] Represented in matrix form as follows:

[0267]

[0268] In the formula, γ xx (i) represents the autocorrelation of the well-side seismic trace data after SFMOD processing, γ xy(i) represents the cross-correlation between well-side seismic traces and inter-well seismic data after SFMOD processing, and np represents the length of the matching factor.

[0269] When using a convolution model to extract matching factors, it is actually an optimization operator.

[0270] a(t) = {a(-l), a(-l+1), ..., a(0), a(1), ..., a(l)}, such that the convolution y(t)*a(t) continuously approaches x(t) while minimizing the sum of squared errors.

[0271] The above expression can be represented as:

[0272]

[0273] The matrix equation is called the Toeplitz matrix equation, and the matching factor at the well location is obtained by recursion.

[0274] Step 3: Obtaining the matching factor spatial variation and inversion wavelet

[0275] Step 3.1: Perform spatial variation of the matching factor using inverse distance interpolation.

[0276] The spatial variation technique for matching factors mainly utilizes the extracted matching factors at the well locations, interpolating them across the entire spatial domain to obtain the matching factor data volume for the entire three-dimensional work area. The inverse distance interpolation method is used for spatial variation processing of the matching factors, and the calculation formula is as follows:

[0277]

[0278] Where j represents the target seismic trace, e[j] is the matching factor of the target seismic trace after spatial variation processing, oper[i][j] represents the matching factor between the target seismic trace j and each point of the current seismic trace, and r 2 It is the distance weight between seismic trace j and the current seismic trace, and np represents the maximum number of seismic traces.

[0279] Step 3.2: Iteratively calculate the inversion wavelet

[0280] Using the reduced-frequency logging reflection coefficient and wellside seismic trace obtained in step 2.2 as input, the convolution model is used with the least squares method to take the 30Hz zero-phase Ricker wavelet as the initial wavelet. By continuously adjusting the Ricker wavelet parameters, the optimal wavelet is finally obtained for subsequent inversion calculations.

[0281] Step 4: Establishing the initial impedance model based on the characteristics of the work area and structural attributes

[0282] Based on the characteristics of fluvial sedimentary models, facies zones are divided using structural attributes. These facies zones serve as constraints, and thick and thin channel regions are established according to the target location. Structural attributes reflecting the characteristics of the inverted target are selected as initial impedance constraint information, and interpolation is performed based on well logging data according to the facies zone division rules.

[0283] mod(t0) = ∑Well(ti+α)·[weight(well)+weight(structural attribute)]

[0284] In the formula, mod represents the impedance value of the current calculation channel, well represents the logging impedance value, weight represents the weight of the well and the weight of the structural attribute, ti represents the longitudinal time range of the well, and α represents the perturbation variable of the structural attribute in the facies zone on the time position of the logging data.

[0285] The initial weight is calculated using the following formula:

[0286]

[0287] weight(well) = 2

[0288] Where b represents the relationship between structural attributes and wellpoint channel sand thickness, and the formula for calculating b is as follows:

[0289]

[0290] In the above formula, the x vector represents the wellpoint channel sand thickness (arranged from smallest to largest), the y vector represents the structural attribute value, and n represents the statistical number of wellpoint sand body thicknesses. When b < 0.5, it indicates that the structural attribute and the wellpoint channel sand thickness have an approximately linear relationship, and the structural attribute accounts for a large proportion; when 0.5 ≤ b ≤ 1, it indicates that the structural attribute and the wellpoint channel sand thickness have a non-linear relationship, and the proportion of the structural attribute decreases; when b > 1, there is no significant relationship between the structural attribute and the wellpoint channel sand thickness, and the structural attribute accounts for the smallest proportion.

[0291] This initial impedance model incorporates structural attribute constraints, and the inter-well interpolation better conforms to geological patterns. Figure 8 (a)); Based on the characteristics of the work area, different weights were assigned, and the resulting initial impedance model not only solved the problem of mismatch between well and seismic data, but also solved the problem that conventional low-frequency models of traditional inversion methods do not have channel morphology. Figure 8 (b)).

[0292] Step 5: Phased-array high-resolution inversion

[0293] A reflection coefficient model is constructed using an initial wave impedance model. The reflection coefficient model is then compared with the structural attribute profile extracted from seismic data. The error between the synthesized seismic structural attribute profile and the well track is then compared. The minimum error required is used as a constraint to iteratively modify the model. By modifying the model, the weights of different attributes in the structural attributes are adjusted so that the high-frequency information in the constraint conditions (well logging data) can guide and enhance the effective information in the seismic data. Various frequency components in the data move closer to the standard of the constraint conditions, ultimately achieving the goal of improving the resolution of seismic data while reducing the difficulty of channel sand characterization.

[0294] In the inversion process, it is assumed that the result to be obtained is The constraints are set as follows: And require the inversion results Gradually moving towards constraints during the iteration process Deflection. Based on the bias requirements of both, the following relationship is obtained:

[0295]

[0296] In the formula, G is the wavelet matrix, G T Let I be the transpose of the matrix, and let I be the identity matrix. The data represents structural attributes, with α and β being constraint factors.

[0297] By modifying the model, broadband, high-resolution information from well logging data at the wellbore location is gradually transferred to the seismic data. The final result of the inversion is a reflection coefficient profile, using a high-resolution seismic wavelet G. * The two are convolved to synthesize a high-resolution seismic data volume:

[0298]

[0299] By using geological models constrained by structural properties as constraints for broadband constrained inversion, the effective information of seismic data is enhanced and guided, thereby improving the resolution of seismic data and achieving the goal of channel sand characterization.

[0300] Beneficial Effects: This invention effectively improves the identification capability of channel sands by combining seismic and well logging data, extracting and analyzing seismic structural attributes, dividing river facies zones, and constraining high-resolution inversion techniques. Well logging currently offers the highest vertical resolution data, while seismic data offers the highest horizontal resolution. Although geophysicists have been committed to improving the vertical resolution of seismic data, the complexity of subsurface conditions, the ambiguity of geophysics, and noise in seismic data have made it difficult to more precisely characterize small-scale geological bodies, such as channel sands, using seismic data. Combining the ability of well logging curves to distinguish strata and identify lithology vertically with the ability to control the distribution of geological bodies in three-dimensional seismic space to better characterize the distribution of channels and channel sands in three-dimensional subsurface space remains a pressing problem for oil and gas resource exploration.

[0301] The phase-controlled high-resolution inversion module based on channel sand structure class attribute constraints combines well logging and seismic analysis. It utilizes seismic structure class attribute facies zone division and well-seismic matching information as constraints, and employs constrained high-resolution inversion to solve the aforementioned problems. For channel sand targets, well-seismic comparison and matching improves the identification capability of well-side seismic traces. Furthermore, constrained high-resolution inversion yields a high-resolution acoustic impedance data volume, whose stratigraphic slices demonstrate excellent channel sand characterization capabilities.

[0302] The above specific embodiments further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above are merely specific embodiments 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 within the scope of protection of the present invention.

Claims

1. A high-resolution inversion method for channel sand facies controlled by structural attribute constraints, characterized in that, The inversion method includes: Step S1: River facies zone division based on structural class attributes; Step S2: Preprocess the logging data and calculate the matching factor at the well location; Step S3: Calculate the inversion wavelet based on the matching factor, and perform space-variant processing on the matching factor; Step S4: Establish an initial impedance model based on the characteristics of the work area and structural attributes; Step S5: Perform phase-controlled high-resolution inversion based on the initial impedance model.

2. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, Step S1: River facies zone division based on structural class attributes specifically includes: By using the spatial distribution of the phase axis of seismic data with wave crests and troughs as the center, the extreme points and zero-crossing points of amplitude within the time window are extracted; Different mathematical calculation methods are used to obtain the waveform structure class attributes contained in three-dimensional data; Based on mathematical calculation methods, waveform structure attributes are classified into integral, statistical, and difference categories. Based on the characteristics of fluvial facies depositional models, structural data is used to divide fluvial facies zones and obtain information on facies-controlled constraints.

3. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, The integral class includes waveform area and waveform length; the statistical class includes coefficient of variation, kurtosis, and skewness; the difference class includes: composite envelope difference, half-time curvature difference, and peak-valley kurtosis difference.

4. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, Step S2: Well logging data preprocessing and calculation of the matching factor at the well location specifically includes: Preprocessing of well logging data includes depth correction, environmental correction, and standardization of well logging data; Extract seismic traces near the well based on well logging data; The matching factor at the well location is calculated using the step-by-step iterative method.

5. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, Step S3: Calculating the inversion wavelet based on the matching factor and performing space-variant processing on the matching factor specifically includes: The inversion wavelet is obtained by iteratively combining the matching factor obtained at the well location with seismic data near the well. The spatial variation processing of matching factors is performed by the correlation coefficient constrained inverse distance interpolation method to obtain the three-dimensional matching factor data volume.

6. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 5, characterized in that, The step of using correlation coefficient-constrained inverse distance interpolation to perform spatial variation processing on the matching factors to obtain the three-dimensional matching factor data volume specifically includes: Using the extracted well location matching factors, interpolation is performed in the entire spatial domain to obtain the matching factor data volume of the entire three-dimensional work area. The matching factor is spatially variable using inverse distance interpolation, and the calculation formula is as follows: where j represents a target seismic trace, e[j] is a matching factor of the target seismic trace after the deconvolution, oper[i][j] represents a matching factor of the target seismic trace j and each point on the current seismic trace, r 2 is a distance weight of the seismic trace j and the current seismic trace, and np represents a maximum number of seismic traces.

7. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, Step S4: Establishing an initial impedance model based on the characteristics of the work area and structural attributes specifically includes: Based on the characteristics of the work area, well logging and seismic weights are set. Based on the distribution range of thick sand bodies and fluvial facies zones, combined with well logging curve interpolation, an initial impedance model with geological understanding is obtained as an inversion constraint. Based on the characteristics of the work area, corresponding weights are assigned to match well seismic data.

8. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, Step S4: Establishing an initial impedance model based on the characteristics of the work area and structural attributes specifically includes: Based on the characteristics of fluvial facies depositional models, facies zones are divided using structural attributes. The divided facies zones are used as constraints to establish thick and thin channel regions according to the target location. Based on the inversion target, structural attributes that can reflect the characteristics of the target are selected as the initial impedance constraint information, and interpolation is performed in combination with well logging data according to the phase zone division rules; mod(t0) = ∑well(ti+α)·[weight(well) + weight (structural class attribute)] In the formula, mod represents the impedance value of the current calculation channel, well represents the logging impedance value, weight represents the weight of the well and the weight of the structural attribute, ti represents the longitudinal time range of the well, and α represents the perturbation variable of the structural attribute in the facies region on the time position of the logging data. The initial weight is calculated using the following formula: weight(well) = 2 Where b represents the relationship between structural attributes and wellpoint channel sand thickness, and the formula for calculating b is as follows: Wherein, the x vector represents the wellpoint channel sand thickness, arranged in ascending order, the y vector represents the structural attribute value, and n represents the statistical number of wellpoint sand body thicknesses; when b < 0.5, it indicates that the structural attribute and the wellpoint channel sand thickness have a nearly linear relationship, and the structural attribute accounts for a large proportion; when 0.5 ≤ b ≤ 1, it indicates that the structural attribute and the wellpoint channel sand thickness have a non-linear relationship, and the proportion of the structural attribute decreases; when b > 1, there is no obvious relationship between the structural attribute and the wellpoint channel sand thickness, and the structural attribute accounts for the smallest proportion.

9. The high-resolution inversion method for channel sand facies controlled by structural attribute constraints according to claim 1, characterized in that, Step S5: Performing phase-controlled high-resolution inversion based on the initial impedance model specifically includes: A reflection coefficient model is constructed based on the initial wave impedance model, and a high-resolution impedance inversion result is obtained by using broadband constrained inversion technology in combination with the inversion wavelet obtained in step S3.

10. The high-resolution inversion method for channel sand facies control with structural attribute constraints according to claim 8, characterized in that, The process of obtaining high-resolution impedance inversion results by combining the inversion wavelet obtained in step S3 with broadband constrained inversion technology specifically includes: In the inversion process, it is assumed that the result to be obtained is The constraints are set as follows: And require the inversion results Gradually moving towards constraints during the iteration process deflection; Based on the bias requirement, the following relation is obtained: In the formula, G is the wavelet matrix, G T Let I be the transpose of the matrix, and let I be the identity matrix. The data is structured attribute data, where α and β are constraint factors; By modifying the model, broadband, high-resolution information from well logging data at the wellbore location is transferred to the seismic data; The final result of the inversion is the reflection coefficient profile, using a high-resolution seismic wavelet G. * The two are convolved to synthesize a high-resolution seismic data volume: Geological models constrained by structural properties are used as constraints for broadband constrained inversion to characterize channel sands.

Images