Multi-domain seismic inversion method and device for fractured porous media

Abstract

The present invention provides a multi-domain seismic inversion method and device for fractured porous media, and relates to the technical field of geological exploration. The present invention adopts a reflection coefficient calculation method based on an accurate Graebner equation, and an objective function constructed by combining a time domain and a frequency domain with a Cauchy prior constraint, thereby realizing simultaneous inversion of three parameter profiles of porosity, fracture weakness and water saturation. The accurate Graebner equation can accurately describe the reflection coefficient characteristics of fractured media, time domain inversion can enhance noise resistance, and frequency domain inversion can improve resolution, thereby taking into account the accuracy, noise resistance and resolution of the inversion result, and realizing high-precision and high-resolution estimation of porosity, fracture weakness and water saturation.

Description

Technical Field

[0001] The invention relates to the technical field of geological exploration, and in particular to a fracture porous medium multi-domain seismic inversion method and device. Background Art

[0002] In the process of oil and gas resource exploration and development, the identification and characterization of fractured reservoirs are of great significance. The pore structure of fractured reservoirs is complex, including multiple media types such as pores, fractures and matrix, and the physical properties of each medium have a significant impact on oil and gas storage and mobility.

[0003] When dealing with fractured porous media, traditional prestack inversion methods usually have difficulty in simultaneously obtaining parameter profiles such as porosity, fracture weakness, and water saturation with high precision, which limits the detailed description of fractured reservoirs. In addition, most existing inversion methods rely on the linearized inversion of the Ruger equation (also known as the Rüger equation). Since the Ruger equation itself has many assumptions, such as weak contrast and small angle incidence, there are certain deviations when calculating the reflection coefficient in the fracture medium model (the fracture medium model is a mathematical and physical model used to describe the fluid flow and elastic wave propagation behavior in a geological body containing fractures or fissures), resulting in limited accuracy of the inversion results. In addition, the inversion results obtained by the existing inversion methods are insufficient in terms of resolution and noise resistance. Summary of the invention

[0004] The purpose of the present invention is to provide a multi-domain seismic inversion method and device for fractured porous media, so as to take into account the accuracy, noise resistance and resolution of the inversion results and realize high-precision and high-resolution estimation of porosity, fracture weakness and water saturation.

[0005] In a first aspect, the present invention provides a multi-domain seismic inversion method for fractured porous media, comprising:

[0006] Obtain well logging data, seismic observation data and geological formation information in the target work area;

[0007] According to the logging data and geological formation information, an initial model of the target work area is constructed; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation. The initial model parameter data of the initial model at the logging position is obtained by filtering and interpolating the real model parameter data corresponding to the logging data;

[0008] According to the preset fracture-porous medium model, the reflection coefficient data at multiple positions of the initial model are calculated using the accurate Graebner equation;

[0009] According to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, the initial model parameter data at the logging position of the initial model are iteratively optimized to obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein the objective function is constructed under the Bayesian prestack inversion framework by combining the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain with the Cauchy prior constraint;

[0010] Determine the optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data;

[0011] According to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function, the initial model parameter data at each position of the initial model are iteratively optimized to obtain the target model parameter profile inversion result of the target work area.

[0012] Furthermore, based on the well logging data and geological formation information, an initial model of the target work area is constructed, including:

[0013] Perform logging interpretation conversion on the logging data to obtain real model parameter data, which includes porosity values, fracture weakness values ​​and water saturation values ​​corresponding to multiple sampling points. The sampling points are used to characterize the depth.

[0014] Perform filtering and interpolation processing on the real model parameter data to obtain the initial model parameter data at the logging location;

[0015] According to the geological layer information and the initial model parameter data at the logging position, the initial model parameter data at other positions in the target work area except the logging position are extrapolated to obtain the initial model.

[0016] Furthermore, according to the preset fracture-pore medium model, the reflection coefficient data at multiple positions of the initial model are calculated using the precise Graebner equation, including:

[0017] According to the fracture pore medium model and the set medium layer parameter data corresponding to the target work area, the stiffness matrix data and equivalent density data at each position of the initial model are determined; wherein the stiffness matrix data includes the stiffness matrix corresponding to each sampling point, and the equivalent density data includes the equivalent density value of the background material layer corresponding to each sampling point;

[0018] According to the stiffness matrix data and equivalent density data at each position of the initial model, the reflection coefficient data at each position of the initial model are calculated using the precise Graebner equation; wherein the reflection coefficient data includes the reflection coefficient values ​​at different incident angles corresponding to each sampling point.

[0019] Further, the correction parameters include a damping factor and a regularization parameter; according to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, the initial model parameter data of the initial model at the logging position is iteratively optimized to obtain the inversion model parameter data under different correction parameter combinations at the logging position, including:

[0020] Determine multiple correction parameter combinations according to the correction parameter value range; wherein the damping factors and / or regularization parameters in different correction parameter combinations are different;

[0021] Traverse each correction parameter combination in turn;

[0022] When traversing to the current correction parameter combination, the reflection coefficient data, sub-seismic data and objective function at the logging position are used to calculate the update amount of the model parameters under the current correction parameter combination at each iteration through the Levenberg-Marquardt algorithm, and the initial model parameter data of the initial model at the logging position are iteratively optimized to obtain the inversion model parameter data under the current correction parameter combination at the logging position;

[0023] Until all correction parameter combinations are traversed, the inversion model parameter data under each correction parameter combination at the logging position is obtained.

[0024] Furthermore, the objective function is:

[0025] ;

[0026] in, J (m) represents the objective function of the model parameter m, s e 2 represents the time domain variance of the noise, s eF 2 represents the frequency domain variance of the noise, s m 2 represents the variance of the model parameter m, d represents the observed seismic data in the time domain, and d F represents the observed seismic data in the frequency domain, S(θ,m) represents the synthetic seismic data in the time domain under the incident angle θ and the model parameter m, and S F'(θ,m) represents the synthetic seismic data in the target form in the frequency domain under the incident angle θ and the model parameter m. The target form is the form of the expansion of the real part and the imaginary part. S(θ,m) is obtained by forward modeling using the convolution of the seismic wavelet and the longitudinal wave reflection coefficient. m represents the regularization parameter, m i Indicates i model parameters, N is the number of sampling points in the time domain.

[0027] Furthermore, the update amount of the model parameters under the current correction parameter combination is calculated by the following formula:

[0028] ;

[0029] in, Indicates t The update amount of the model parameters at the iteration, A represents the Hessian matrix of the objective function, the calculation formula of A includes the damping factor and the regularization parameter, and C represents the gradient of the objective function.

[0030] Further, according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data, the optimal correction parameter combination is determined, including:

[0031] According to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data, the optimal correction parameter combination is selected from each correction parameter combination by the two-dimensional L curve method.

[0032] In a second aspect, the present invention further provides a multi-domain seismic inversion device for fractured porous media, comprising:

[0033] An acquisition module is used to obtain well logging data, observation seismic data and geological formation information of the target work area;

[0034] A construction module is used to construct an initial model of the target work area based on the logging data and geological layer information; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation. The initial model parameter data of the initial model at the logging position is obtained by filtering and interpolating the real model parameter data corresponding to the logging data;

[0035] A calculation module, used to calculate the reflection coefficient data at multiple positions of the initial model using the precise Graebner equation according to the preset fracture porous medium model;

[0036] An optimization module is used to iteratively optimize the initial model parameter data at the logging position according to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, so as to obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein the objective function is constructed under the Bayesian prestack inversion framework by combining the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain with the Cauchy prior constraint;

[0037] A determination module, used to determine the optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data;

[0038] The inversion module is used to iteratively optimize the initial model parameter data at each position of the initial model according to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function, and obtain the target model parameter profile inversion result of the target work area.

[0039] In a third aspect, the present invention further provides an electronic device, including a memory and a processor, wherein the memory stores a computer program that can be run on the processor, and when the processor executes the computer program, the multi-domain seismic inversion method for fractured porous media of the first aspect is implemented.

[0040] In a fourth aspect, the present invention further provides a computer-readable storage medium having a computer program stored thereon, and when the computer program is executed by a processor, the multi-domain seismic inversion method for fractured porous media of the first aspect is executed.

[0041] The multi-domain seismic inversion method and device for fractured porous media provided by the embodiment of the present invention can obtain well logging data, observed seismic data and geological layer information of a target work area; construct an initial model of the target work area according to the well logging data and geological layer information; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation, and the initial model parameter data of the initial model at the well logging position is obtained by filtering and interpolating the real model parameter data corresponding to the well logging data; according to the preset fractured porous medium model, the reflection coefficient data at multiple positions of the initial model are calculated using the precise Graebner equation; according to the reflection coefficient data at the well logging position, the sub-seismic data corresponding to the well logging position in the observed seismic data, and the pre-constructed The objective function and the preset correction parameter value range are used to iteratively optimize the initial model parameter data at the logging position of the initial model, and obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein, the objective function is constructed under the Bayesian prestack inversion framework with the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain combined with the Cauchy prior constraint; the optimal correction parameter combination is determined according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data; the initial model parameter data at each position of the initial model are iteratively optimized according to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function, and the target model parameter profile inversion result of the target work area is obtained.

[0042] In this way, the reflection coefficient calculation method based on the precise Graebner equation and the objective function constructed by combining the time domain and the frequency domain with the help of the Cauchy prior constraint are adopted to realize the simultaneous inversion of the three parameter profiles of porosity, fracture weakness and water saturation. The precise Graebner equation can accurately describe the reflection coefficient characteristics of the fractured medium, the time domain inversion can enhance the noise resistance, and the frequency domain inversion can improve the resolution. Therefore, the accuracy, noise resistance and resolution of the inversion results are taken into account, and high-precision and high-resolution estimation of porosity, fracture weakness and water saturation is realized. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0044] Figure 1 A schematic flow chart of a multi-domain seismic inversion method for fractured porous media provided by an embodiment of the present invention;

[0045] Figure 2 A schematic flow chart of another multi-domain seismic inversion method for fractured porous media provided by an embodiment of the present invention;

[0046] Figure 3 A schematic structural diagram of a multi-domain seismic inversion device for fractured porous media provided by an embodiment of the present invention;

[0047] Figure 4 A schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0048] The technical solution of the present invention will be clearly and completely described below in conjunction with the embodiments. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0049] At present, the traditional prestack inversion method is difficult to simultaneously obtain parameter profiles such as porosity, fracture weakness and water saturation with high precision, and relies on the linear inversion of the Ruger equation, resulting in the accuracy of the inversion results being limited; in addition, the inversion process is usually only carried out in the time domain or frequency domain, ignoring the potential of the combination of the two, making the inversion results insufficient in terms of resolution or noise resistance. Based on this, a multi-domain seismic inversion method and device for fractured porous media provided in an embodiment of the present invention adopts an accurate Graebner equation, which can more accurately describe the reflection coefficient characteristics of fractured media, and uses the sparsity prior constraint of the Cauchy distribution to jointly construct the objective function in the time domain and the frequency domain, taking into account the accuracy, noise resistance and resolution of the inversion results, and can achieve high-precision and high-resolution estimation of porosity, fracture weakness and water saturation, providing a new technical means for the effective characterization of complex fractured reservoirs.

[0050] To facilitate understanding of this embodiment, a multi-domain seismic inversion method for fractured porous media disclosed in an embodiment of the present invention is first introduced in detail.

[0051] The embodiment of the present invention provides a multi-domain seismic inversion method for fractured porous media, which can be executed by an electronic device with data processing capability. Figure 1 The flowchart of a multi-domain seismic inversion method for fractured porous media is shown in FIG. 1 , and the method mainly includes the following steps S110 to S160:

[0052] Step S110, acquiring well logging data, observed seismic data and geological layer information of the target work area.

[0053] The target work area is the work area that needs to perform parameter profile inversion. The parameter profiles include porosity profile, fracture weakness profile and water saturation profile. Among them, porosity refers to the ratio of the pore volume in the rock to the total rock volume; fracture weakness refers to the impact of naturally existing cracks in the rock or other geological structures on the stability of the formation; water saturation refers to the proportion of pore space occupied by water in the rock.

[0054] The above-mentioned logging data may include logging parameter data of a logging position selected in the target work area at multiple sampling points (i.e., the time after the logging depth is converted to the time domain), and the logging parameter data of each sampling point includes multiple logging parameter values, and the logging parameters may include one or more of natural gamma logging (GR), resistivity logging (Rt), sonic logging (DT), density logging (RHOB), neutron porosity logging (NPHI) and natural potential logging (SP). The horizontal axis of the logging data is the logging parameter value, and the vertical axis is the time (representing the depth).

[0055] The above-mentioned observed seismic data can be the angle track gather data of the target work area, the ordinate is time (that is, the sampling point, reflecting the time / depth from the emission to the reception of the reflected wave), and the abscissa is the position information of the common depth point (Common Depth Point, CDP for short) arranged along the survey line (that is, the position along the survey line direction); the angle track gather data of a stacked track (corresponding to a common depth point) includes the amplitude values ​​of different sampling points corresponding to multiple incident angles.

[0056] The above geological layer information is used to indicate the position of the target layer on the observed seismic data, and the target layer refers to the geological layer in the target work area that needs to be parameter profile inverted. In a possible implementation, the horizontal coordinate of the geological layer information can be the horizontal coordinate of the observed seismic data (i.e., the position information of the CDP, without the amplitude value), and the vertical coordinate can be the vertical coordinate of the observed seismic data (i.e., time).

[0057] Step S120, constructing an initial model of the target work area based on the logging data and geological layer information; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation. The initial model parameter data of the initial model at the logging position is obtained by filtering and interpolating the real model parameter data corresponding to the logging data.

[0058] In some possible embodiments, the above step S120 may include the following steps S121 to S123:

[0059] Step S121, performing logging interpretation conversion on the logging data to obtain true model parameter data, the true model parameter data including porosity values, fracture weakness values ​​and water saturation values ​​corresponding to a plurality of sampling points, the sampling points being used to characterize the depth.

[0060] The above logging interpretation conversion method can refer to the relevant existing technology and will not be described here. Through logging interpretation conversion, a set of parameter values ​​(including porosity value, fracture weakness value and water saturation value) can be calculated at each depth under the logging position.

[0061] Step S122, filtering and interpolation processing are performed on the real model parameter data to obtain initial model parameter data at the logging position.

[0062] In this step, the real model parameter data can be filtered and smoothed by a preset filtering algorithm to adjust the real model parameter data while retaining the overall trend of the data; the filtered real model parameter data is easier to use for interpolation and extrapolation when establishing the initial model, and can be used for comparison with the real model parameter data after inversion to determine the optimal correction parameter combination of the initial model. Based on the number of sampling points in the observed seismic data, new sampling point data can be constructed in the filtered real model parameter data by a preset interpolation algorithm, so that the number of sampling points of the initial model parameter data obtained after interpolation is consistent with that of the observed seismic data.

[0063] The above-mentioned filtering algorithm and interpolation algorithm can be selected according to actual needs, and this embodiment does not limit this. For example, the filtering algorithm can adopt moving average method or wavelet transform, and the interpolation algorithm can adopt linear interpolation, polynomial interpolation, spline interpolation or nearest neighbor interpolation.

[0064] Step S123, based on the geological layer information and the initial model parameter data at the logging position, the initial model parameter data at other positions in the target work area except the logging position are extrapolated to obtain an initial model.

[0065] The initial model parameter data at other locations in the target work area except the well logging location can be predicted by a preset extrapolation algorithm based on the initial model parameter data at the well logging location. The extrapolation algorithm can be selected according to actual needs, and this embodiment does not limit this. For example, extrapolation methods such as linear extrapolation, polynomial extrapolation, regression analysis or neural network are used. Each location in the target work area refers to the location of the target work area on the horizontal plane, such as the geographical location represented by the longitude and latitude coordinates.

[0066] Step S130, based on the preset fracture porous medium model, the reflection coefficient data at multiple positions of the initial model are calculated using the precise Graebner equation.

[0067] In this embodiment, a fractured porous medium model is established based on porosity, fracture weakness and water saturation. The fractured porous medium model includes multiple medium layer parameters, such as bulk modulus, shear modulus, density, porosity, fracture thickness and normalized frequency. In the fractured porous medium model, the fractured porous medium is divided into a relatively thick background material layer (with limited porosity) and a porous medium layer (with limited porosity). ) and a relatively thin layer of high-porosity material representing fractures ( →1) Two categories. The fracture-porous medium model characterizes the relationship between the stiffness matrix and porosity, fracture weakness and water saturation, as well as the calculation method of the equivalent density of the background material layer. The fracture-porous medium model will be described in detail later.

[0068] In some possible embodiments, the stiffness matrix data and equivalent density data at each position of the initial model can be determined based on the fracture pore medium model and the set medium layer parameter data corresponding to the target work area; wherein the stiffness matrix data includes the stiffness matrix corresponding to each sampling point, and the equivalent density data includes the equivalent density value of the background material layer corresponding to each sampling point; then, based on the stiffness matrix data and equivalent density data at each position of the initial model, the reflection coefficient data at each position of the initial model is calculated using the precise Graebner equation; wherein the reflection coefficient data includes the reflection coefficient values ​​at different incident angles corresponding to each sampling point. The specific calculation method of the reflection coefficient data will be described in detail later.

[0069] Step S140, iteratively optimize the initial model parameter data of the initial model at the logging position according to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, and obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein the objective function is constructed under the Bayesian prestack inversion framework with the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain combined with the Cauchy prior constraint.

[0070] In some possible embodiments, the correction parameters include a damping factor and a regularization parameter; the step S140 may be implemented by the following process:

[0071] Step S141, determining a plurality of correction parameter combinations according to the correction parameter value range; wherein the damping factors and / or regularization parameters in different correction parameter combinations are different;

[0072] Step S142, traversing each correction parameter combination in sequence;

[0073] Step S143, when traversing to the current correction parameter combination, using the reflection coefficient data, sub-seismic data and objective function at the logging position, the update amount of the model parameters under the current correction parameter combination at each iteration is calculated by the Levenberg-Marquardt algorithm, and the initial model parameter data of the initial model at the logging position is iteratively optimized to obtain the inversion model parameter data under the current correction parameter combination at the logging position;

[0074] Step S144, until all correction parameter combinations are traversed, the inversion model parameter data under each correction parameter combination at the logging position is obtained.

[0075] Alternatively, a possible objective function is:

[0076] ;

[0077] in, J (m) represents the objective function of the model parameter m, s e 2 represents the time domain variance of the noise, s eF 2 represents the frequency domain variance of the noise, s e 2 and s eF 2 All are preset values; s m 2 represents the variance of the model parameter m, d represents the observed seismic data in the time domain, and d F represents the observed seismic data in the frequency domain, S(θ,m) represents the synthetic seismic data in the time domain under the incident angle θ and the model parameter m, and S F '(θ,m) represents the synthetic seismic data in the target form of the frequency domain under the incident angle θ and the model parameter m. The target form is the form of the expansion of the real part and the imaginary part. S(θ,m) is obtained by forward modeling using the convolution of the seismic wavelet and the longitudinal wave reflection coefficient. m represents the regularization parameter, m i Indicates i model parameters, N is the number of sampling points in the time domain.

[0078] Optionally, the update amount of the model parameters under the current correction parameter combination can be calculated by the following formula:

[0079] ;

[0080] in, Indicatest The update amount of the model parameters at the iteration, A represents the Hessian matrix of the objective function, the calculation formula of A includes the damping factor and the regularization parameter, and C represents the gradient of the objective function.

[0081] Step S150, determining the optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data.

[0082] The inversion model parameter data under different correction parameter combinations at the logging position can be compared with the real model parameter data corresponding to the logging data, and the correction parameter combination with the highest similarity and the smallest error can be selected as the correction parameter combination.

[0083] In some possible embodiments, the optimal correction parameter combination can be selected from various correction parameter combinations using a two-dimensional L curve method based on the inversion model parameter data under different correction parameter combinations at the logging location and the real model parameter data corresponding to the logging data.

[0084] Step S160, iteratively optimize the initial model parameter data at each position of the initial model according to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function, and obtain the target model parameter profile inversion result of the target work area.

[0085] The optimal correction parameter combination can be applied to all the superimposed traces in the target work area, and the model parameter inversion is performed using the angle trace data and the corresponding reflection coefficient data of each superimposed trace to obtain the inversion model parameter data at the corresponding position in the target work area; the inversion model parameter data at each position in the target work area is determined as the target model parameter profile inversion result of the target work area. The inverted porosity data at each position is the porosity profile at that position, the fracture weakness data at each position is the fracture weakness profile at that position, and the water saturation data at each position is the water saturation profile. It should be noted that each superimposed trace is inverted separately, and the specific inversion process can refer to the corresponding content in the aforementioned step S140, which will not be repeated here.

[0086] The multi-domain seismic inversion method for fractured porous media provided in an embodiment of the present invention adopts a reflection coefficient calculation method based on the precise Graebner equation and an objective function constructed by combining the time domain and the frequency domain with the help of the Cauchy prior constraint, thereby realizing the simultaneous inversion of three parameter profiles of porosity, fracture weakness and water saturation. The precise Graebner equation can accurately describe the reflection coefficient characteristics of the fractured medium, the time domain inversion can enhance the noise resistance, and the frequency domain inversion can improve the resolution. Therefore, the accuracy, noise resistance and resolution of the inversion result are taken into account, and high-precision and high-resolution estimation of porosity, fracture weakness and water saturation is realized.

[0087] For ease of understanding, the above-mentioned multi-domain seismic inversion method for fractured porous media is described in detail below. ), crack weakness ( d N ) and water saturation ( S w ) is a new method for prestack inversion with high accuracy and resolution. , d N , S w Based on the fracture porous medium model, the complete stiffness matrix is ​​obtained, and then the longitudinal wave reflection coefficient is calculated by the accurate Graebner equation to establish the , d N , S w The relationship between the seismic wavelet and the P-wave reflection coefficient is established by using the convolution of the seismic wavelet and the P-wave reflection coefficient. Under the Bayesian prestack inversion framework, the residuals of the simulated data in the time domain and frequency domain and the actual observed seismic data are combined with the Cauchy prior constraint to construct the objective function. The Levenberg-Marquardt (LM) method is used for least squares iteration to finally obtain , d N , S w Three parameters, time domain inversion can enhance noise immunity, and frequency domain inversion can improve resolution. Figure 2 Another flow chart of a multi-domain seismic inversion method for fractured porous media is shown in FIG. 1 , which includes the following seven steps:

[0088] Step 1: Get the data.

[0089] Obtain logging data, seismic data and geological formation information of the target area. Logging data include natural gamma logging (GR), resistivity logging (Rt), acoustic logging (DT), density logging (RHOB), neutron porosity logging (NPHI) and natural potential logging (SP) and other logging parameters. Through these basic logging parameters, logging interpretation can be performed to convert the required porosity ( ), crack weakness ( d N ) and water saturation ( S w). Seismic data is angular gather data with a high signal-to-noise ratio. It is used as actual observation data. The horizontal direction is the common depth point, including amplitude data of multiple incident angles, and the vertical direction is time (Time). The layer information is used to interpolate and extrapolate the initial model of the logging location.

[0090] Step 2: Establish a fracture-porous medium model.

[0091] Fractured porous media include thicker layers of background material (with limited porosity) ) and a relatively thin layer of high-porosity material representing fractures ( →1), the pores and cracks are filled with different fluids. For sufficiently low wave frequencies (referring to the frequency of the compression waves transmitted from the seismic source to the medium), b and c Two ingredients ( b Refers to the composition of the background material layer and the fluid, that is, the background material layer part. c Refers to the compression modulus of a layered fluid-saturated fractured porous medium composed of a high-porosity material layer and a fluid, i.e., a high-porosity material layer. C 33 It can be written as follows:

[0092] (1);

[0093] in, h j Indications j The crack thickness, C j Indicates the composition of the fluid when it is saturated j The longitudinal wave modulus, α j Indications j The Biot effective stress coefficient, M j Indications j The pore space modulus, i represents the imaginary unit, oh represents the angular frequency of the oscillation frequency, or j Indications j The viscosity of L j Indications j The dry longitudinal wave modulus, k j Indications j The permeability of the fractured porous medium is assumed to be made of the same isotropic granular material. H represents the spatial period of the fractured porous medium, j The value of b orc , that is j = b , c , b and c The two components are saturated by different fluids 。

[0094] Composition of fluid saturated b and c The longitudinal wave modulus C j , can be calculated by the Gassmann equation:

[0095] (2);

[0096] in, K j Indications j The dry bulk modulus, K g represents the bulk modulus of the fractured porous medium, Indications j The porosity, K fj Indications j The bulk modulus of the fluid, m j Indications j The shear modulus, j = b , c.c. The composition represents the crack part, its thickness (volume fraction) h c and dry crack material modulus ( K c , m c ) are very small, and formula (1) can be approximated as:

[0097] (3);

[0098] Where Ω represents the normalized frequency, d N Indicates crack weakness 0-1, F is a dimensionless parameter.

[0099] (4);

[0100] in, Z N Represents the normal crack flexibility.

[0101] The fractured porous medium can be represented by a frequency-dependent function using a relaxed (low frequency) and non-relaxed (high frequency) model, which can itself be equivalent to a transversely isotropic medium (VTI). The complete elastic stiffness matrix is C ij As shown below:

[0102] (5);

[0103] Among them, the stiffness matrix C The elements in are elastic parameters; and denote the elastic parameter limits of the non-relaxed model and the relaxed model, respectively. R ( oh ) represents a C 33 , and The scalar relaxation function composed of is as follows:

[0104] (6).

[0105] In the low-frequency limit, low-frequency fluctuations provide sufficient time for the flow of fluids. The fluid pressures in the cracks and background pores are balanced within the seismic wave cycle, and they can flow to each other. Based on this, we can obtain:

[0106] (7);

[0107] (8);

[0108] (9);

[0109] (10).

[0110] In the high-frequency limit, the fluids do not have time to flow into each other during the seismic wave cycle, and the fluid pressures in the pores and fractures are no longer equal, showing acoustic and elastic characteristics that are different from those in the low-frequency limit. Based on this, we can obtain:

[0111] (11);

[0112] (12).

[0113] According to Krief formula, we can get b Dry bulk modulus and shear modulus of ingredients K b , m b as follows:

[0114] (13);

[0115] in, m g Represents the shear modulus of fractured porous media.

[0116] Assuming that the fluid in the medium is composed of water and air, the bulk modulus and density are expressed as K 1. K 2. r 1. r 2 means that according to the Wood equation and the density weighted average formula, it can be calculated b Equivalent fluid bulk modulus of the components K fb and equivalent density r b (The effect of cracks on density is negligible):

[0117] (14);

[0118] in, r fb represents the equivalent fluid density of a mixture of water and air, r g Represents the density of fractured porous media.

[0119] This step has completed the modeling of the fractured porous medium and obtained a complete stiffness matrix. The parameter data of the medium layer in the embodiment of the present invention are as follows: K c =0.01; K g =37; KFC =2; m g =44; m c =0.01; r g =2.65 g / cm³; =0.99; h c =0.001; K 1=2; K 2=0.0015;Ω=0.5Hz; r 1 = 1 g / cm³; r 2=0.0012 g / cm³, the unit of modulus is GPa. , d N , S w Combined ( ,S w For fractured porous media b The porosity and water saturation of the components can be calculated based on the fracture pore model established in this step. , d N , S w The complete elastic stiffness matrix corresponding to the combination C And the equivalent density r b , which is used for subsequent reflection coefficient calculation.

[0120] Step 3: Calculate the reflection coefficient using the exact Graebner equation.

[0121] A transversely isotropic medium (VTI) with a vertical symmetry axis can be described by five elastic stiffness parameters, among which C 11 , C 13 , C 33 and C 55 Determines the speed of quasi-longitudinal waves (qP) and quasi-transverse waves (qSV), C 66 Determines the speed of horizontal shear waves (SH). For qP waves incident on the interface separating two transversely isotropic homogeneous media (located in the xy plane), Graebner proposed an accurate analytical expression for the reflection coefficient of VTI media to describe its amplitude variation with the incident angle (Amplitude VersusAngle, AVA for short). This expression defines the reflection and transmission coefficient vectors corresponding to the incident P wave when a plane wave is incident, as follows:

[0122] (15);

[0123] ;

[0124] in, r PP represents the reflected longitudinal wave coefficient, r PS represents the reflected shear wave coefficient, t PP represents the transmitted longitudinal wave coefficient, t PS represents the transmission shear wave coefficient, the superscript U , D Indicates upper layer, lower layer, subscript P , S Represents P wave, S wave, and intermediate parameters in matrix S aj , b j , c j , d j ( j = U , D ) is expressed as follows:

[0125] (16);

[0126] Among them, direction cosines l P , l S , m P , m S , l and m They represent the components of the wave in the vertical and horizontal directions respectively, and the expressions are as follows:

[0127] (17);

[0128] in, and They represent the vertical slowness of qP wave and qSV wave respectively, and their expressions can be expressed by horizontal slowness p And elastic parameter C is expressed as:

[0129] (18);

[0130] ;

[0131] ;

[0132] ;

[0133] ;

[0134] ;

[0135] ;

[0136] in, i represents the angle of incidence. Equations (15)-(18) give the reflection coefficient r PP about C 11 , C 13 , C 33 , C44 and r b precise expression of .

[0137] Step 4: Construct the objective function.

[0138] In the Bayesian inversion framework, the objective function is constructed by combining the time domain and frequency domain information with the maximum a posteriori probability of the model parameters combined with the Cauchy distribution. Specifically:

[0139] The forward model of pre-stack angle gather data in the time domain is the wavelet matrix and the reflection coefficient rPP The convolution under different incident angles is called the Robinson convolution model, and the expression is as follows:

[0140] (19);

[0141] Among them, S(θ) represents the synthetic seismic data in the time domain, and W(θ) represents the wavelet matrix (calculated from the angle gather data). In the frequency domain, the Robinson convolution model can be Fourier transformed to obtain:

[0142] (20);

[0143] in, r PP ( t ,θ) represents the reflection coefficient in the time domain, α ( t ) represents the time sampling point, and the formula (20) is discretized. At a single incident angle, it can be written in the following matrix form:

[0144] (twenty one);

[0145] in, r PPk ( k =1, 2, … , N ) indicates that N sampling points are N reflection coefficients; oh i ( i =1, 2, …, M ) represents the ones involved in the inversion M The matrix S is F and W F is a complex number. To simplify the inversion process, it is expanded into expressions of real and imaginary parts:

[0146] (twenty two).

[0147] Therefore, under multiple incident angles, the frequency domain forward equation is as follows:

[0148] (twenty three).

[0149] The embodiment of the present invention applies the joint inversion method of time-frequency domain, transforms the seismic data in the time domain into the frequency domain through Fourier transformation, and combines the response characteristics of the two domains. In the Bayesian inversion framework, the construction of the likelihood function usually depends on the probability density distribution model of the background noise. Assuming that the seismic data in the time domain and the frequency domain are independent of each other, their joint likelihood function can be expressed by multiplying the probability distribution between the independent variables:

[0150] (twenty four);

[0151] Where d represents the time domain seismic data, d F represents the frequency domain seismic data after selecting the effective frequency, m=[ , d N , S w ]T represents the model parameters to be inverted, and the likelihood function p (d,d F |m) represents the matching degree between the synthetic data of the inversion result and the actual time-frequency domain seismic data. p (d|m) and p (d F |m) represent the matching degree between the synthetic data of the inversion results and the seismic data in the time domain and frequency domain, respectively.

[0152] Assuming random noise e , and the time domain seismic data d and the frequency domain seismic data d F All conform to Gaussian distribution, then the likelihood function can be expressed as:

[0153] (25);

[0154] in, s e 2 and s eF 2 Respectively represent the variance of noise in the time domain and frequency domain, N and M Respectively represent the number of sampling points in the time domain and the number of effective frequency components in the frequency domain, S(θ,m) and S F ′(θ,m) represents the forward modeling process of the model parameter m in the time domain and frequency domain respectively. Combining formulas (24) and (25), the likelihood function of the joint inversion in the time and frequency domain is p (d,d F |m) can be written as:

[0155] (26);

[0156] The Bayesian inversion framework builds the posterior distribution p (m | d,d F ) and the likelihood function and prior distribution p Cauchy (m), the specific formula is as follows:

[0157] (27).

[0158] Prior information is obtained by probabilistically analyzing the existing data. In Bayesian inversion theory, prior information is used to constrain and stabilize the inversion problem, so that the posterior probability distribution is closer to the actual situation. The embodiment of the present invention uses the Cauchy prior distribution as a constraint, and has better performance and stronger noise resistance.

[0159] (28);

[0160] in, s m 2 represents the variance of the model parameter m. Combining equations (26)-(28) and taking the logarithm of the posterior distribution, we get the objective function J (m):

[0161] (29);

[0162] in, m represents the regularization parameter of the Cauchy distribution.

[0163] Step 5: Optimize the objective function.

[0164] First, the initial model is established, and then the LM algorithm is used to optimize and derive the corresponding partial derivative matrix. Specifically:

[0165] First, an initial model of the entire target area needs to be established, which can be obtained by combining geological layer information and logging data through interpolation and extrapolation. The initial model needs to be consistent with the geological trend to ensure its rationality and provide a reliable starting point for inversion. The Levernberg-Marquardt (LM) algorithm is used to solve the objective function J (m), calculate the model update amount , t Represents the number of iterations, and the expression is as follows:

[0166] (30);

[0167] ;

[0168] ;

[0169] Among them, A represents the Hessian matrix of the objective function, C represents the gradient of the objective function, l and m They represent the damping factor and regularization parameter respectively. The two-dimensional L-curve method is used to select the appropriate correction parameter combination. I represents the unit matrix. and They represent the first-order derivative (gradient) and second-order derivative (Hessian matrix) of the Cauchy regularization parameter, respectively, and are expressed as follows:

[0170] (31);

[0171] Q=W×X,Q F =W F′ ×X, where X represents the Jacobian matrix, which is the reflection coefficient and the inversion model parameters ( , d N , S w ), the calculation accuracy of which directly affects the final result of the LM algorithm. In order to improve the calculation accuracy and expand the angle range in which the inversion can be applied, the embodiment of the present invention establishes the precise Jacobian matrix of the Graebner equation of the fracture void model:

[0172] (32);

[0173] in, L =3, indicating that there are three parameters to be inverted. H represents the number of incident angles, , , According to formula (15), we can get:

[0174] (33).

[0175] According to the first step obtained fracture pore model stiffness matrix C And the chain rule gives and The expression is as follows:

[0176] ;

[0177] ;

[0178] ;

[0179] ;

[0180] ;

[0181] .

[0182] Step 6: Quality control of inversion results.

[0183] The inversion results at the logging position are compared with the real logging curve. Within the value range of the correction parameters λ and μ, the two-dimensional L-curve method is used to select the most appropriate correction parameter combination so that the inversion results have the highest similarity with the real logging and the smallest error.

[0184] Step 7: Complete the inversion.

[0185] If the expectations are met, the correction parameter combination obtained in the sixth step is applied to all the trace sets (i.e., stacked traces) in the work area. The embodiment of the present invention does not require the approximate conditions of traditional linear inversion, nor does it require trace integration of the inversion results to cause error accumulation. Instead, it can directly obtain the inversion profiles of porosity, fracture weakness and water saturation with high accuracy and resolution.

[0186] The embodiment of the present invention constructs a fracture porous medium model and calculates the reflection coefficient based on the precise Graebner equation. Compared with the traditional Ruger equation, the result is more accurate. The Cauchy distribution is used to perform joint inversion in the time domain and the frequency domain under the Bayesian inversion framework, and the result is more accurate and has a higher resolution.

[0187] Corresponding to the above-mentioned multi-domain seismic inversion method for fractured porous media, an embodiment of the present invention further provides a multi-domain seismic inversion device for fractured porous media. Figure 3 The structure diagram of a multi-domain seismic inversion device for fractured porous media is shown, and the device comprises:

[0188] Acquisition module 301, used to acquire well logging data, observed seismic data and geological layer information of the target work area;

[0189] The construction module 302 is used to construct an initial model of the target work area according to the logging data and geological layer information; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation. The initial model parameter data of the initial model at the logging position is obtained by filtering and interpolating the real model parameter data corresponding to the logging data;

[0190] The calculation module 303 is used to calculate the reflection coefficient data at multiple positions of the initial model using the precise Graebner equation according to the preset fracture porous medium model;

[0191] The optimization module 304 is used to iteratively optimize the initial model parameter data at the logging position according to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, so as to obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein the objective function is constructed under the Bayesian prestack inversion framework by combining the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain with the Cauchy prior constraint;

[0192] A determination module 305 is used to determine an optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data;

[0193] The inversion module 306 is used to iteratively optimize the initial model parameter data at each position of the initial model according to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function, so as to obtain the target model parameter profile inversion result of the target work area.

[0194] The multi-domain seismic inversion device for fractured porous media provided in an embodiment of the present invention adopts a reflection coefficient calculation method based on the precise Graebner equation and an objective function constructed by combining the time domain and the frequency domain with the help of the Cauchy prior constraint, thereby realizing the simultaneous inversion of three parameter profiles of porosity, fracture weakness and water saturation. The precise Graebner equation can accurately describe the reflection coefficient characteristics of the fractured medium, the time domain inversion can enhance the noise resistance, and the frequency domain inversion can improve the resolution. Therefore, the accuracy, noise resistance and resolution of the inversion result are taken into account, and high-precision and high-resolution estimation of porosity, fracture weakness and water saturation is realized.

[0195] Furthermore, the above-mentioned construction module 302 is specifically used for: performing logging interpretation conversion on the logging data to obtain real model parameter data, the real model parameter data including porosity values, fracture weakness values ​​and water saturation values ​​corresponding to multiple sampling points, and the sampling points are used to characterize the depth; filtering and interpolation processing are performed on the real model parameter data to obtain initial model parameter data at the logging position; based on the geological layer information and the initial model parameter data at the logging position, the initial model parameter data at other positions in the target work area except the logging position are extrapolated to obtain the initial model.

[0196] Furthermore, the above-mentioned calculation module 303 is specifically used to: determine the stiffness matrix data and equivalent density data at each position of the initial model according to the set medium layer parameter data corresponding to the fracture pore medium model and the target work area; wherein the stiffness matrix data includes the stiffness matrix corresponding to each sampling point, and the equivalent density data includes the equivalent density value of the background material layer corresponding to each sampling point; according to the stiffness matrix data and the equivalent density data at each position of the initial model, the reflection coefficient data at each position of the initial model is calculated using the precise Graebner equation; wherein the reflection coefficient data includes the reflection coefficient values ​​under different incident angles corresponding to each sampling point.

[0197] Furthermore, the correction parameters include damping factors and regularization parameters; the optimization module 304 is specifically used to: determine multiple correction parameter combinations according to the correction parameter value range; wherein the damping factors and / or regularization parameters in different correction parameter combinations are different; traverse each correction parameter combination in turn; when traversing to the current correction parameter combination, use the reflection coefficient data, sub-seismic data and objective function at the logging position, calculate the update amount of the model parameters under the current correction parameter combination at each iteration through the Levenberg-Marquardt algorithm, iteratively optimize the initial model parameter data of the initial model at the logging position, and obtain the inversion model parameter data under the current correction parameter combination at the logging position; until all correction parameter combinations are traversed, the inversion model parameter data under each correction parameter combination at the logging position is obtained.

[0198] Furthermore, the above objective function is:

[0199] ;

[0200] in, J (m) represents the objective function of the model parameter m, s e 2 represents the time domain variance of the noise, s eF 2 represents the frequency domain variance of the noise, s m 2 represents the variance of the model parameter m, d represents the observed seismic data in the time domain, and d F represents the observed seismic data in the frequency domain, S(θ,m) represents the synthetic seismic data in the time domain under the incident angle θ and the model parameter m, and S F '(θ,m) represents the synthetic seismic data in the target form in the frequency domain under the incident angle θ and the model parameter m. The target form is the form of the expansion of the real part and the imaginary part. S(θ,m) is obtained by forward modeling using the convolution of the seismic wavelet and the longitudinal wave reflection coefficient. m represents the regularization parameter,m i Indicates i model parameters, N is the number of sampling points in the time domain.

[0201] Furthermore, the update amount of the model parameters under the above current correction parameter combination is calculated by the following formula:

[0202] ;

[0203] in, Indicates t The update amount of the model parameters at the iteration, A represents the Hessian matrix of the objective function, the calculation formula of A includes the damping factor and the regularization parameter, and C represents the gradient of the objective function.

[0204] Furthermore, the determination module 305 is specifically used to select the optimal correction parameter combination from various correction parameter combinations by a two-dimensional L curve method according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data.

[0205] The multi-domain seismic inversion device for fractured porous media provided in this embodiment has the same implementation principle and technical effects as those of the aforementioned multi-domain seismic inversion method for fractured porous media embodiment. For the sake of brief description, for matters not mentioned in the embodiment of the multi-domain seismic inversion device for fractured porous media, reference may be made to the corresponding contents in the aforementioned multi-domain seismic inversion method for fractured porous media embodiment.

[0206] like Figure 4 As shown, an electronic device 400 provided by an embodiment of the present invention includes: a processor 401, a memory 402 and a bus, the memory 402 stores a computer program that can be run on the processor 401, when the electronic device 400 is running, the processor 401 and the memory 402 communicate through the bus, and the processor 401 executes the computer program to implement the above-mentioned multi-domain seismic inversion method for fractured porous media.

[0207] Specifically, the memory 402 and the processor 401 can be general-purpose memories and processors, which are not specifically limited here.

[0208] The embodiment of the present invention further provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the multi-domain seismic inversion method for fractured porous media in the previous method embodiment is executed. The computer-readable storage medium includes: a USB flash drive, a mobile hard disk, a read-only memory (ROM), a RAM, a magnetic disk, or an optical disk, and other media that can store program codes.

[0209] The term "and / or" herein is only a description of the association relationship of the associated objects, indicating that there may be three relationships. For example, A and / or B can represent: A exists alone, A and B exist at the same time, and B exists alone. In addition, the term "at least one" herein represents any combination of at least two of any one or more of a plurality of. For example, including at least one of A, B, and C can represent including any one or more elements selected from the set consisting of A, B, and C.

[0210] In all examples shown and described herein, any specific values ​​should be interpreted as merely exemplary and not as limiting, and thus other examples of the exemplary embodiments may have different values.

[0211] The flowcharts and block diagrams in the accompanying drawings show the possible architecture, functions and operations of the devices, methods and computer program products according to multiple embodiments of the present invention. In this regard, each box in the flowchart or block diagram can represent a module, a program segment or a part of a code, and a part of a module, a program segment or a code contains one or more executable instructions for realizing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the box can also occur in a different order from the order marked in the accompanying drawings. For example, two consecutive boxes can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the block diagram and / or the flowchart, and the combination of the boxes in the block diagram and / or the flowchart can be implemented with a dedicated hardware-based system that performs a specified function or action, or can be implemented with a combination of dedicated hardware and computer instructions.

[0212] In the several embodiments provided in the present application, it should be understood that the disclosed devices and methods can be implemented in other ways. The device embodiments described above are only schematic. For example, the division of modules is only a logical function division. There may be other division methods in actual implementation. For example, multiple modules or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some communication interfaces, and the indirect coupling or communication connection of devices or modules can be electrical, mechanical or other forms.

[0213] The modules described as separate components may or may not be physically separated, and the components shown as modules may or may not be physical modules, that is, they may be located in one place or distributed on multiple network modules. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0214] In addition, each functional module in each embodiment of the present invention may be integrated into one processing module, or each module may exist physically separately, or two or more modules may be integrated into one module.

[0215] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A multi-domain seismic inversion method for fractured porous media, characterized in that: include: Obtain well logging data, seismic observation data and geological formation information in the target work area; An initial model of the target work area is constructed based on the well logging data and the geological layer information; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation, and the initial model parameter data of the initial model at the well logging position is obtained by filtering and interpolating the real model parameter data corresponding to the well logging data; According to the preset fracture porous medium model, the reflection coefficient data at multiple positions of the initial model are calculated using the accurate Graebner equation; According to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, the initial model parameter data of the initial model at the logging position is iteratively optimized to obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein the objective function is constructed under the Bayesian prestack inversion framework by combining the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain with the Cauchy prior constraint; Determining an optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data; Iteratively optimize the initial model parameter data at each position of the initial model according to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function to obtain the target model parameter profile inversion result of the target work area; The initial model of the target work area is constructed based on the well logging data and the geological layer information, including: Performing logging interpretation conversion on the logging data to obtain the real model parameter data, wherein the real model parameter data includes porosity values, fracture weakness values ​​and water saturation values ​​corresponding to a plurality of sampling points, wherein the sampling points are used to characterize the depth; Performing filtering and interpolation processing on the real model parameter data to obtain initial model parameter data at the logging position; According to the geological layer information and the initial model parameter data at the logging position, the initial model parameter data at other positions in the target work area except the logging position are extrapolated to obtain the initial model.

2. The method according to claim 1, characterized in that The reflection coefficient data at multiple positions of the initial model are calculated using the precise Graebner equation according to the preset fracture pore medium model, including: According to the fracture pore medium model and the set medium layer parameter data corresponding to the target work area, the stiffness matrix data and the equivalent density data at each position of the initial model are determined; wherein the stiffness matrix data includes the stiffness matrix corresponding to each sampling point, and the equivalent density data includes the equivalent density value of the background material layer corresponding to each sampling point; According to the stiffness matrix data and the equivalent density data at each position of the initial model, the reflection coefficient data at each position of the initial model are calculated using the precise Graebner equation; wherein the reflection coefficient data includes the reflection coefficient values ​​at different incident angles corresponding to each sampling point.

3. The method according to claim 1, characterized in that The correction parameters include a damping factor and a regularization parameter; the initial model parameter data of the initial model at the logging position is iteratively optimized according to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, the pre-constructed objective function and the preset correction parameter value range, to obtain the inversion model parameter data under different correction parameter combinations at the logging position, including: Determining a plurality of correction parameter combinations according to the correction parameter value range; wherein the damping factors and / or regularization parameters in different correction parameter combinations are different; Traversing each correction parameter combination in turn; When traversing to the current correction parameter combination, the reflection coefficient data, the sub-seismic data and the objective function at the logging position are used to calculate the update amount of the model parameters under the current correction parameter combination at each iteration through the Levenberg-Marquardt algorithm, and the initial model parameter data of the initial model at the logging position are iteratively optimized to obtain the inversion model parameter data under the current correction parameter combination at the logging position; Until all the correction parameter combinations are traversed, the inversion model parameter data under each correction parameter combination at the logging position is obtained.

4. The method according to claim 3, characterized in that The objective function is: ； in, J (m) represents the objective function of the model parameter m, σ e 2 represents the time domain variance of the noise, σ eF 2 represents the frequency domain variance of the noise, σ m 2 represents the variance of the model parameter m, d represents the observed seismic data in the time domain, and d F represents the observed seismic data in the frequency domain, S(θ,m) represents the synthetic seismic data in the time domain under the incident angle θ and the model parameter m, and S F '(θ,m) represents synthetic seismic data in the target form of the frequency domain under the incident angle θ and the model parameter m, wherein the target form is the form of the expansion of the real part and the imaginary part. S(θ,m) is obtained by forward modeling using the convolution of the seismic wavelet and the longitudinal wave reflection coefficient. μ represents the regularization parameter, m i Indicates i model parameters, N is the number of sampling points in the time domain.

5. The method according to claim 3, characterized in that: The update amount of the model parameters under the current correction parameter combination is calculated by the following formula: ； in, Indicates t The update amount of the model parameters at the iteration, A represents the Hessian matrix of the objective function, the calculation formula of A includes the damping factor and the regularization parameter, and C represents the gradient of the objective function.

6. The method according to claim 1, characterized in that The determining of the optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data comprises: According to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data, the optimal correction parameter combination is selected from each correction parameter combination by a two-dimensional L curve method.

7. A multi-domain seismic inversion device for fractured porous media, characterized in that: include: An acquisition module is used to obtain well logging data, observation seismic data and geological formation information of the target work area; A construction module is used to construct an initial model of the target work area according to the well logging data and the geological layer information; wherein the initial model is used to characterize the size of model parameters at different positions, and the model parameters include porosity, fracture weakness and water saturation, and the initial model parameter data of the initial model at the well logging position is obtained by filtering and interpolating the real model parameter data corresponding to the well logging data; A calculation module, used to calculate the reflection coefficient data at multiple positions of the initial model using the precise Graebner equation according to a preset fracture pore medium model; An optimization module, for iteratively optimizing the initial model parameter data of the initial model at the logging position according to the reflection coefficient data at the logging position, the sub-seismic data corresponding to the logging position in the observed seismic data, a pre-constructed objective function and a preset correction parameter value range, so as to obtain the inversion model parameter data under different correction parameter combinations at the logging position; wherein the objective function is constructed under the Bayesian prestack inversion framework by combining the residuals of the simulated seismic data and the observed seismic data corresponding to the time domain and the frequency domain with the Cauchy prior constraint; A determination module, used to determine an optimal correction parameter combination according to the inversion model parameter data under different correction parameter combinations at the logging position and the real model parameter data corresponding to the logging data; An inversion module, for iteratively optimizing the initial model parameter data at each position of the initial model according to the optimal correction parameter combination, the reflection coefficient data at each position of the initial model, the observed seismic data and the objective function, to obtain an inversion result of the target model parameter profile of the target work area; The construction module is specifically used for: performing logging interpretation conversion on the logging data to obtain the real model parameter data, the real model parameter data including porosity values, fracture weakness values ​​and water saturation values ​​corresponding to a plurality of sampling points, the sampling points being used to characterize the depth; performing filtering and interpolation processing on the real model parameter data to obtain the initial model parameter data at the logging position; and extrapolating the initial model parameter data at other positions in the target work area except the logging position according to the geological layer information and the initial model parameter data at the logging position to obtain the initial model.

8. An electronic device comprising a memory and a processor, wherein the memory stores a computer program that can be run on the processor, characterized in that: When the processor executes the computer program, the multi-domain seismic inversion method for fractured porous media according to any one of claims 1 to 6 is implemented.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the multi-domain seismic inversion method for fractured porous media according to any one of claims 1 to 6 is performed.

Images