High-resolution wave-impedance inversion method and apparatus based on spectral fusion

Abstract

A high-resolution wave-impedance inversion method and apparatus based on spectral fusion. The method comprises: under the constraints of horizon data and a spatial variation function, performing Kriging spatial interpolation processing on a low frequency band of well logging data to obtain a low-frequency-band wave-impedance data body zL; according to seismic wavelets and a mid-frequency model of Cauchy norm-constrained sparse spike inversion, performing Cauchy norm-constrained sparse spike inversion on seismic data to obtain a reflection coefficient, and performing trace integration to obtain a mid-frequency-band relative wave-impedance data body zM; under the constraints of the horizon data and the spatial variation function, performing sequential Gaussian stochastic simulation on a high frequency band of the well logging data to obtain a high-frequency-band wave-impedance data body zH; performing Fourier transform on zL, zM and zH to obtain frequency domain data bodies ZL, ZM and ZH; performing spectral fusion to obtain a full-frequency-band wave-impedance data body ZMGI; and performing inverse Fourier transform to obtain a high-resolution wave-impedance inversion result zMGI. High-resolution wave-impedance inversion is therefore effectively performed.

Description

High-resolution wave impedance inversion method and device based on spectrum fusion

[0001] Related applications

[0002] This application claims priority to Chinese patent application No. 202310677711.1 filed on June 8, 2023, and cites the contents disclosed in the above patent application as part of this application. Technical Field

[0003] The present application relates to the field of geophysical exploration technology, and in particular to a high-resolution wave impedance inversion method and device based on spectrum fusion. Background Art

[0004] This section is intended to provide a background or context to the embodiments of the present application that are recited in the claims. No admission is made that the description herein is prior art by virtue of its inclusion in this section.

[0005] Seismic inversion is the process of combining seismic, well logging, and geological data to reveal the spatial structure and physical properties of subsurface rock formations. It is currently one of the most effective methods for finding oil and gas reservoirs. Wave impedance (the product of velocity and density) is a key formation property that links geology and geophysics. Determining wave impedance through seismic inversion eliminates the influence of seismic wavelets and converts interface information from seismic data into stratigraphic information, thereby characterizing the spatial distribution of reservoirs and evaluating their properties.

[0006] Seismic impedance inversion can generally be divided into two categories: deterministic inversion and stochastic inversion. Deterministic inversion involves obtaining a single optimal solution estimate under specific solution constraints. This method is computationally efficient, but due to the limitations of the seismic frequency band, the inversion results cannot identify thin reservoir information outside the seismic frequency band and cannot evaluate the uncertainty of the inversion results. Stochastic inversion, on the other hand, uses a sampling strategy to fully explore the posterior solution space, capable of outputting high-resolution inversion results. It has great application potential in thin reservoir prediction and has attracted considerable attention in recent years.

[0007] Depending on whether the inverse problem is linear or nonlinear, stochastic inversion can be divided into two categories: iterative optimization nonlinear stochastic inversion and Bayesian linearized stochastic inversion. In the field of iterative optimization nonlinear stochastic inversion, Haas and Dubrule (1994) first proposed geostatistical inversion. This method uses seismic data as a constraint to control the selection of random simulation results of logging data in the full frequency band. The computational efficiency is low, and the track-by-track calculation method makes the inversion results have poor lateral continuity and lack of low-frequency information in the vertical direction, which is inconsistent with the actual geological conditions. Sambridge and Mosegaard (2002) proposed a seismic inversion method based on Monte Carlo random sampling. This method also requires a large number of "rejection-acceptance" processes and has low computational efficiency. Azevedo et al. (2020) combined co-simulation with direct sequential Gaussian random simulation to perform random realization of the entire volume, improving the lateral continuity of the inversion results, but faced the problem of difficulty in convergence. In Bayesian linearized stochastic inversion, it is often assumed that model parameters follow a multidimensional Gaussian distribution and that a linear forward relationship exists between the model parameters and the seismic data. This allows for the analytical expression of the posterior probability density function of the model parameters under Bayesian theory. This posterior probability density function is then sampled to obtain the stochastic inversion results (e.g., Buland and Omre, 2003; Hansen et al., 2006; Grana et al., 2021). This method has the advantage of high computational efficiency. However, when performing calculations on the entire data volume, if lateral spatial correlation is considered, the memory overhead of the model covariance matrix becomes extremely large, making it difficult to implement on conventional computers. Therefore, this method generally only considers the longitudinal variogram during inversion, resulting in poor lateral continuity in the stochastic implementation results. Furthermore, the assumptions of the Gaussian prior model and the linear forward operator are often inconsistent with actual conditions.

[0008] Judging from the published literature, in the field of seismic exploration, there is currently no effective method to integrate seismic inversion and geostatistics in the frequency domain to perform high-resolution wave impedance inversion.

[0009] Summary of the Invention

[0010] The present application provides a high-resolution impedance inversion method based on spectral fusion to improve the accuracy of traditional deterministic inversion, solve the problem that traditional random inversion requires a large number of iterations and has low computational efficiency, and effectively perform high-resolution impedance inversion. The method includes:

[0011] Acquire seismic data, seismic wavelets, well logging data, and horizon data;

[0012] Under the constraints of the layer data and spatial variation function, the low-frequency band of the logging data is processed by Kriging spatial interpolation to obtain the low-frequency band wave impedance data volume z L ;

[0013] A Cauchy norm constrained sparse pulse inversion intermediate frequency model is established. Based on the selected seismic wavelet and the Cauchy norm constrained sparse pulse inversion intermediate frequency model, the seismic data is subjected to Cauchy norm constrained sparse pulse inversion to obtain the reflection coefficient. The reflection coefficient is integrated to obtain the intermediate frequency relative wave impedance data volume z. M ;

[0014] Under the constraints of the layer data and spatial variation function, the high frequency band of the logging data is simulated sequentially by Gaussian random simulation to obtain the high frequency band wave impedance data volume z H ;

[0015] The low frequency band wave impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、Z H ;

[0016] According to the selected frequency parameters and scale parameters, the frequency domain data volume Z L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ;

[0017] For the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0018] The present application also provides a high-resolution impedance inversion device based on spectrum fusion, which is used to improve the accuracy of traditional deterministic inversion, solve the problem that traditional random inversion requires a large number of iterations and has low computational efficiency, and effectively perform high-resolution impedance inversion. The device includes:

[0019] Data acquisition module, used to acquire seismic data, seismic wavelets, well logging data and horizon data;

[0020] The Kriging spatial interpolation module is used to perform Kriging spatial interpolation processing on the low-frequency band of the logging data under the constraints of the layer data and the spatial variation function to obtain the low-frequency band wave impedance data volume z L ;

[0021] The Cauchy norm constrained sparse pulse inversion module is used to establish a Cauchy norm constrained sparse pulse inversion intermediate frequency model. According to the selected seismic wavelet and the Cauchy norm constrained sparse pulse inversion intermediate frequency model, the seismic data is subjected to Cauchy norm constrained sparse pulse inversion to obtain the reflection coefficient; the reflection coefficient is integrated to obtain the intermediate frequency relative wave impedance data volume z M ;

[0022] The sequential Gaussian random simulation module is used to perform sequential Gaussian random simulation on the high frequency band of the logging data under the constraints of the layer data and the spatial variation function to obtain the high frequency band wave impedance data volume z H ;

[0023] Fourier transform module is used to transform the low frequency band wave impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、Z H ;

[0024] Spectral fusion module is used to fusion the frequency domain data volume Z according to the selected frequency parameters and scale parameters. L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ;

[0025] The inverse Fourier transform module is used to transform the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0026] An embodiment of the present application also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the above-mentioned high-resolution wave impedance inversion method based on spectrum fusion is implemented.

[0027] An embodiment of the present application further provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-mentioned high-resolution wave impedance inversion method based on spectrum fusion.

[0028] An embodiment of the present application further provides a computer program product, which includes a computer program. When the computer program is executed by a processor, it implements the above-mentioned high-resolution wave impedance inversion method based on spectrum fusion.

[0029] In the embodiment of the present application, seismic data, seismic wavelets, well logging data and horizon data are obtained; under the constraints of the horizon data and the spatial variation function, the low-frequency band of the well logging data is subjected to Kriging spatial interpolation processing to obtain the low-frequency band wave impedance data volume z L ; Establish a Cauchy norm constrained sparse pulse inversion intermediate frequency model, perform Cauchy norm constrained sparse pulse inversion on seismic data according to the selected seismic wavelet and Cauchy norm constrained sparse pulse inversion intermediate frequency model, and obtain the reflection coefficient; perform trace integration on the reflection coefficient to obtain the intermediate frequency relative wave impedance data volume z M Under the constraints of the layer data and spatial variation function, the high frequency band of the logging data is simulated sequentially by Gaussian random simulation to obtain the high frequency band wave impedance data volume z H ; respectively for the low frequency band wave impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、Z H ; According to the selected frequency parameters and scale parameters, the frequency domain data volume Z L 、Z M 、 Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ; For the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0030] By combining seismic inversion with geostatistical theory and adopting the idea of ​​"divide and conquer", that is, processing separately, decoupling the broadband inversion results, obtaining low-frequency band, medium-frequency band and high-frequency band results respectively, and finally integrating them to obtain high-precision wave impedance inversion results, the accuracy of traditional deterministic inversion is improved, and the problem of low computational efficiency caused by traditional random inversion, which requires a large number of iterations to achieve, is solved. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the following briefly introduces the drawings required for the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present application. For those skilled in the art, other drawings can be obtained based on these drawings without creative work. In the drawings:

[0032] FIG1 is a processing flow chart of a high-resolution wave impedance inversion method based on spectrum fusion in an embodiment of the present application;

[0033] FIG2 is a flow chart of a method for establishing a Cauchy norm constrained sparse pulse inversion intermediate frequency model in an embodiment of the present application;

[0034] FIG3 is a schematic diagram of the fusion of wave impedance data volumes and spectra in different frequency bands according to an embodiment of the present application;

[0035] FIG4 is a logic flow chart of a high-resolution wave impedance inversion method based on spectrum fusion in an embodiment of the present application;

[0036] FIG5 is a schematic diagram of an actual seismic profile of an oil field in an embodiment of the present application;

[0037] FIG6 (a) is a conventional sparse pulse inversion result in an embodiment of the present application; FIG6 (b) is a high-resolution impedance inversion result in an embodiment of the present application;

[0038] FIG7( a ) is a schematic diagram of a root mean square amplitude slice of the inversion result with a high cutoff frequency of 120 Hz in an embodiment of the present application, and FIG7( b ) is a schematic diagram of a root mean square amplitude slice of the inversion result with a high cutoff frequency of 150 Hz in an embodiment of the present application;

[0039] FIG8 is a schematic structural diagram of a high-resolution wave impedance inversion device based on spectrum fusion in an embodiment of the present application;

[0040] FIG9 is a schematic diagram of the computer device structure according to an embodiment of the present application. DETAILED DESCRIPTION

[0041] In order to make the purpose, technical solutions and advantages of the embodiments of the present application more clear, the embodiments of the present application are further described in detail below with reference to the accompanying drawings. Here, the illustrative embodiments of the present application and their descriptions are used to explain the present application, but are not intended to limit the present application.

[0042] The applicant found that, regardless of whether it is deterministic inversion or random inversion, under the action of seismic wavelet convolution (filtering), the inversion results in the mid-frequency band controlled by the seismic data are consistent, that is, the high-frequency band of the inversion results will not affect the degree of matching between the synthetic record and the actual seismic record. This is also the embodiment of the multi-solution nature of geophysical inverse problems.

[0043] Based on this, in order to solve the problems of low computational efficiency of traditional geostatistical inversion, poor lateral continuity of inversion results, missing low-frequency information in the vertical direction, and inconsistency with actual geological conditions, the applicant combined deterministic seismic inversion with geostatistical Kriging estimation and random simulation technology, and comprehensively utilized logging, seismic and geological spatial structure information to provide a high-resolution wave impedance inversion method and device based on spectral fusion, aiming to improve the resolution of wave impedance inversion results, thereby improving the prediction accuracy of thin reservoirs in oil field development.

[0044] The specific technical solution of this application is introduced below.

[0045] FIG1 is a processing flow chart of a high-resolution wave impedance inversion method based on spectrum fusion in an embodiment of the present application. As shown in FIG1 , the high-resolution wave impedance inversion method based on spectrum fusion in an embodiment of the present application may include:

[0046] Step 101: Acquire seismic data, seismic wavelets, well logging data, and horizon data;

[0047] Step 102: Under the constraints of the layer data and the spatial variation function, perform Kriging spatial interpolation processing on the low-frequency band of the logging data to obtain the low-frequency band wave impedance data volume z L ;

[0048] Step 103: Establish a Cauchy norm constrained sparse pulse inversion intermediate frequency model. Perform Cauchy norm constrained sparse pulse inversion on the seismic data based on the selected seismic wavelet and the Cauchy norm constrained sparse pulse inversion intermediate frequency model to obtain the reflection coefficient. Perform trace integration on the reflection coefficient to obtain the intermediate frequency relative wave impedance data volume z. M ;

[0049] Step 104: Under the constraints of the layer data and the spatial variation function, perform sequential Gaussian random simulation on the high frequency band of the logging data to obtain the high frequency band wave impedance data volume z H ;

[0050] Step 105: respectively calculate the low frequency impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、Z H ;

[0051] Step 106: Based on the selected frequency parameters and scale parameters, the frequency domain data volume Z is L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ;

[0052] Step 107: Full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0053] After obtaining the seismic data, seismic wavelets, well logging data and stratigraphic data of the specified underground rock formation, the Kriging low-frequency background model can be established, the low-frequency soft-constrained sparse pulse inversion intermediate-frequency volume calculation can be performed, the sequential Gaussian random simulation high-frequency model can be established, and frequency domain integration can be performed.

[0054] Kriging spatial interpolation technology is a generalized regression method that obtains estimation results with minimum error variance according to the least squares criterion. It not only considers the distance between known points and points to be estimated, but also considers the geological correlation between points. It is widely used in the field of earth science.

[0055] Under different data assumptions, there are different types of Kriging interpolation methods, such as simple Kriging interpolation, ordinary Kriging interpolation, trend Kriging interpolation, etc. Simple Kriging interpolation is suitable for estimating stationary variables whose regional mean is known and constant, while ordinary Kriging interpolation is suitable for estimating non-stationary variables whose regional mean is unknown. During the calculation process, the local mean is automatically calculated from the neighborhood. Since wave impedance is non-stationary in space, this application adopts the ordinary Kriging interpolation method to perform spatial interpolation on the wave impedance logging curve under the constraint of the layer to obtain the low-frequency wave impedance data volume z L .

[0056] For step 102, in one embodiment, the low frequency band of the logging data is subjected to Kriging spatial interpolation processing to obtain the low frequency band wave impedance data volume z L , which may include: using the ordinary Kriging interpolation method to perform Kriging spatial interpolation processing on the low frequency band of the logging data to obtain the low frequency band wave impedance data volume z L .

[0057] Regarding step 103, the process of establishing the intermediate frequency model of the Cauchy norm constrained sparse pulse inversion is described below:

[0058] FIG2 is a flow chart of a method for establishing a Cauchy norm constrained sparse pulse inversion intermediate frequency model in an embodiment of the present application. As shown in FIG2 , in one embodiment, the process of establishing a Cauchy norm constrained sparse pulse inversion intermediate frequency model includes:

[0059] Step 201: construct an objective function with the goal of minimizing the error between seismic data and synthetic seismic records, wherein the seismic data is the actual observed seismic data and the synthetic seismic records are obtained by convolution of seismic wavelets and reflection coefficients;

[0060] Step 202: Based on the objective function, add the Cauchy norm constraint of the reflection coefficient to obtain a new objective function;

[0061] Step 203: Obtain the derivative of the new objective function with respect to the reflection coefficient, and set the new objective function to 0, thereby obtaining a Cauchy norm constrained sparse pulse inversion intermediate frequency model.

[0062] In specific implementation, assuming that the post-stack seismic record can be described by the Robinson convolution model, then d obs =Wr+n (1)

[0063] Among them, d obs represents the observed seismic data, W represents the wavelet convolution matrix, r represents the reflection coefficient, and n is the noise term.

[0064] The purpose of seismic inversion is to obtain the observed seismic data from obs The reflection coefficient r is estimated so that the error between the synthetic record and the actual observation data is minimized. At this time, the objective function can be written as: min J(r)=(Wr-d obs ) T (Wr-d obs ) (2)

[0065] Among them, minJ(r) represents the minimum error between seismic data and synthetic seismic records, that is, the minimum value of the objective function, d obs represents seismic data, W represents the seismic wavelet convolution matrix, and r represents the reflection coefficient.

[0066] However, due to the influence of noise in the observation data, directly solving the above objective function will cause strong instability in the inversion result. It is necessary to further add a constraint term on the reflection coefficient in the objective function. Compared with the L2 norm constraint term, the Cauchy norm constraint can produce sparse reflection coefficients with full frequency band and better noise resistance. This application adds the Cauchy norm constraint on the reflection coefficient on the basis of the traditional inversion objective function. The new objective function is

[0067] Where minJ(r) represents the minimum error between seismic data and synthetic seismic records, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, J s (r) is the seismic data constraint term, which is used to control the similarity between the synthetic seismic record and the seismic data. Cauchy (r) is the Cauchy norm constraint term of the reflection coefficient, which is used to control the sparsity of the reflection coefficient. μ is the weight coefficient of the sparsity constraint. n is the number of sampling points of the seismic trace. i represents the sampling point sequence number of the seismic trace. σ represents the variance of the reflection coefficient.

[0068] Obtain the derivative of the new objective function (Formula 3) with respect to the reflection coefficient r and set it to 0. Then we can get (W T W+μQ)r=W T d obs (4)

[0069] Among them, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, μ is the weight coefficient of the sparse constraint, Q is a diagonal weighted matrix, and the diagonal elements i represents the sampling point number of the seismic trace, and σ represents the variance of the reflection coefficient. Since the matrix Q contains the reflection coefficient r to be determined, Equation 4 is an implicit nonlinear system of equations and cannot be solved directly. This application uses the iterative reweighted least squares method (IRLS) to solve it.

[0070] Finally, the relative wave impedance can be obtained by performing trace integral calculation on the reflection coefficient inversion result. The generated mid-frequency relative impedance data volume is denoted as z M Where z(t) represents the relative wave impedance, r represents the reflection coefficient, t represents the time corresponding to a certain sampling point of the seismic trace, and τ represents the seismic trace sampling point before time point t.

[0071] Regarding step 104, the specific process of establishing the sequential Gaussian random simulation high-frequency model is described below:

[0072] Geostatistical stochastic simulation can generate multiple high-resolution geological models that conform to known data and underground spatial structural characteristics. Each model has equal probability and can be used to establish static and dynamic reservoir models, as well as to perform uncertainty assessment. It is one of the key core technologies for detailed reservoir description. For continuous variables (such as wave impedance), sequential Gaussian simulation is the most commonly used geostatistical stochastic simulation algorithm, which perfectly combines the Gaussian random function model with the sequential simulation algorithm. The advantage of the Gaussian random function is that it can be fully described using only the mean vector and covariance matrix (calculated by the variogram); the sequential simulation algorithm can simplify the problem of obtaining a complex high-dimensional space joint distribution to the problem of obtaining a low-dimensional space conditional distribution, thereby achieving random sampling of high-dimensional data.

[0073] This application uses sequential Gaussian simulation to build a high-frequency model for well logging data under the control of horizons. The specific implementation steps are:

[0074] ①Perform normal transformation on logging data;

[0075] Since the basic assumption of the Sequential Gaussian Simulation algorithm is that the data follows a Gaussian distribution, it is necessary to first perform a normal (Gaussian) transformation on the logging data to make it meet the basic assumption of Sequential Gaussian Simulation. This step also obtains logging data, but the data distribution has changed and will be transformed back in the end. This step is used to prepare the logging data for Sequential Gaussian Simulation, thereby obtaining a high-frequency body. This is the step of establishing a high-frequency model. In the specific embodiment of this application, only the high-frequency body of the Sequential Gaussian model is reflected, and the step of establishing the high-frequency model is not further explored.

[0076] ②Determine a random path and visit all nodes;

[0077] Each sampling point in a three-dimensional data volume can be called a "node", and the path is the order in which all sampling points are calculated. This order is random; for example, the 100th sampling point may be calculated first, and then the 1000th sampling point may be calculated next. This random path is the calculation order of all sampling points, which is the same as the Gaussian distribution and is also required by the Sequential Gaussian Simulation Algorithm.

[0078] ③ At each node, conditional data consisting of the wave impedance of adjacent hard data points (logging data) and simulated points is obtained, and the conditional cumulative probability distribution function (CCDF) of the current node is calculated using the Kriging interpolation method and the variogram.

[0079] The hard data points here refer to the points measured by logging data, which are equivalent to known points.

[0080] ④ Randomly obtain a value from CCDF as the value of the node and add the simulated value to the dataset;

[0081] The dataset here refers to the collection of hard data points and previously calculated points. These points are treated as known points when the next point is calculated. The dataset is the collection of known points. When all points (wave impedance) are calculated, the high-frequency wave impedance data volume is obtained.

[0082] ⑤ Return to step ③ until all nodes are calculated;

[0083] ⑥ The simulation results are transformed into the original space through inverse normal transformation, and finally the high-frequency wave impedance data volume is obtained, which is recorded as z H .

[0084] For step 105, after obtaining the low-frequency wave impedance data volume z according to the above steps, L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Afterwards, the low-frequency impedance data volume z can be L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、Z H Then in step 106, the frequency domain data volume Z is selected according to the frequency parameter and scale parameter. L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI .

[0085] FIG3 is a schematic diagram of the spectrum fusion of wave impedance data volumes of different frequency bands in an embodiment of the present application. Accordingly, in one embodiment, the frequency domain data volume Z can be fused according to the following formula based on the selected frequency parameters and scale parameters: L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI : Z MGI =Z L H(0,0,f1,f2)+αZ M H(f1,f2,f3,f4)+βZ H H(f1,f2,f3,f4)

[0086] Among them, Z MGI It represents the full-band wave impedance data volume, α represents the scale factor of the mid-frequency data volume, β represents the scale factor of the high-frequency data volume, which is used to adjust the amplitude of the components in each frequency band, H represents a filter with an amplitude of 1, and f1, f2, f3, and f4 represent the frequency parameters of the frequency domain filter.

[0087] Finally, in step 107, the full-band impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0088] To simply and clearly demonstrate the logic of this application, Figure 4 shows a logic flow chart of the high-resolution wave impedance inversion method based on spectrum fusion in an embodiment of this application. Figure 4 helps quickly understand the specific processing logic of the high-resolution wave impedance inversion method based on spectrum fusion in an embodiment of this application.

[0089] This application can obtain high-resolution wave impedance inversion results, which has a good application effect in predicting thin reservoirs in oil field development. Figure 5 is a schematic diagram of an actual seismic profile of an oil field in an embodiment of this application. As shown in Figure 5, the curve inserted in the profile is a wave impedance logging curve. The main frequency of the seismic data is approximately 30Hz. The target layer is a thin sandstone reservoir. Due to the limited range of the seismic frequency band, these thin reservoirs cannot be directly identified in the seismic profile.

[0090] Figure 6 (a) shows the results of conventional sparse pulse inversion in an embodiment of the present application, while Figure 6 (b) shows the results of high-resolution impedance inversion in an embodiment of the present application. As indicated by the arrows, it can be seen that the resolution of the impedance inversion results obtained in this application is significantly higher than that of conventional sparse pulse inversion results, clearly depicting the spatial distribution characteristics of thin reservoirs.

[0091] (a) in FIG7 is a schematic diagram of the root mean square amplitude slice of the inversion result with a high cutoff frequency of 120 Hz in an embodiment of the present application, and (b) in FIG7 is a schematic diagram of the root mean square amplitude slice of the inversion result with a high cutoff frequency of 150 Hz in an embodiment of the present application. As shown in (a) in FIG7 and (b) in FIG7 , the root mean square amplitude slice comparison of the inversion results with different high cutoff frequencies is shown. From the comparison between the two, it can be seen that although the high-frequency components of the inversion results are different, the overall consistency on the slice is good, and the lateral distribution of the inversion results is clear and natural, without the "bull's eye" shaped modeling phenomenon, which also reflects that the lateral information of the inversion results of the present application mainly depends on seismic data and has good predictability.

[0092] The present application also provides a high-resolution wave impedance inversion device based on spectrum fusion, as described in the following embodiments. Since the principle of solving the problem of this device is similar to that of the high-resolution wave impedance inversion method based on spectrum fusion, the implementation of this device can refer to the implementation of the high-resolution wave impedance inversion method based on spectrum fusion, and the repeated parts will not be repeated here.

[0093] FIG8 is a schematic diagram of the structure of a high-resolution wave impedance inversion device based on spectrum fusion in an embodiment of the present application. As shown in FIG8 , the electric vehicle use risk detection device in an embodiment of the present application may specifically include:

[0094] Data acquisition module 801, used to acquire seismic data, seismic wavelets, well logging data and horizon data;

[0095] The Kriging spatial interpolation module 802 is used to perform Kriging spatial interpolation processing on the low frequency band of the logging data under the constraints of the layer data and the spatial variation function to obtain the low frequency band wave impedance data volume z L ;

[0096] The Cauchy norm constrained sparse pulse inversion module 803 is used to establish a Cauchy norm constrained sparse pulse inversion intermediate frequency model, perform Cauchy norm constrained sparse pulse inversion on the seismic data according to the selected seismic wavelet and the Cauchy norm constrained sparse pulse inversion intermediate frequency model, and obtain the reflection coefficient; perform trace integration on the reflection coefficient to obtain the intermediate frequency band relative wave impedance data volume z M ;

[0097] The sequential Gaussian random simulation module 804 is used to perform sequential Gaussian random simulation on the high frequency band of the logging data under the constraints of the layer data and the spatial variation function to obtain the high frequency band wave impedance data volume z H ;

[0098] Fourier transform module 805 is used to transform the low frequency band wave impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume zH Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、Z H ;

[0099] Spectrum fusion module 806 is used to fusion the frequency domain data volume Z according to the selected frequency parameters and scale parameters. L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ;

[0100] The inverse Fourier transform module 807 is used to transform the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0101] In one embodiment, the Kriging spatial interpolation module 802 is specifically configured to:

[0102] The ordinary Kriging method is used to perform Kriging spatial interpolation processing on the low-frequency band of the logging data to obtain the low-frequency band wave impedance data volume z L .

[0103] In one embodiment, the Cauchy norm constrained sparse pulse inversion module 803 is specifically configured to:

[0104] An objective function is constructed with the goal of minimizing the error between seismic data and synthetic seismic records, where seismic data refers to actual observed seismic data and synthetic seismic records are obtained by convolution of seismic wavelets and reflection coefficients.

[0105] Based on the objective function, the Cauchy norm constraint of the reflection coefficient is added to obtain a new objective function;

[0106] The derivative of the new objective function with respect to the reflection coefficient is obtained and the new objective function is set to 0 to obtain the Cauchy norm constrained sparse pulse inversion intermediate frequency model.

[0107] In one embodiment, the objective function can be expressed as follows: minJ(r)=(Wr-d obs ) T (Wr-d obs )

[0108] Among them, minJ(r) represents the minimum error between seismic data and synthetic seismic records, that is, the minimum value of the objective function, d obs represents seismic data, W represents the seismic wavelet convolution matrix, and r represents the reflection coefficient.

[0109] In one embodiment, the new objective function can be expressed as follows:

[0110] Where minJ(r) represents the minimum error between seismic data and synthetic seismic records, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, J s (r) is the seismic data constraint term, which controls the similarity between the synthetic seismic record and the seismic data. Cauchy (r) is the Cauchy norm constraint term of the reflection coefficient, which controls the sparsity of the reflection coefficient, μ is the weight coefficient of the sparsity constraint, n is the number of sampling points of the seismic trace, i represents the sampling point sequence number of the seismic trace, and σ represents the variance of the reflection coefficient.

[0111] In one embodiment, the Cauchy norm constrained sparse pulse inversion intermediate frequency model is expressed as follows: (W T W+μQ)r=W T d obs

[0112] Among them, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, μ is the weight coefficient of the sparse constraint, Q is a diagonal weighted matrix, and the diagonal elements i represents the sampling point number of the seismic trace, and σ represents the variance of the reflection coefficient.

[0113] In one embodiment, the Cauchy norm constrained sparse pulse inversion module 803 is specifically configured to:

[0114] Perform trace integration on the reflection coefficient to obtain the relative wave impedance;

[0115] According to the relative wave impedance, the mid-frequency band relative wave impedance data volume z is obtained M .

[0116] In one embodiment, the Cauchy norm constrained sparse pulse inversion module 803 is specifically configured to perform trace integration on the reflection coefficient according to the following formula to obtain the relative wave impedance:

[0117] Where z(t) represents the relative wave impedance, r represents the reflection coefficient, t represents the time corresponding to a certain sampling point of the seismic trace, and τ represents the seismic trace sampling point before time point t.

[0118] In one embodiment, the spectrum fusion module 806 is specifically configured to: according to the following formula, based on the selected frequency parameters and scale parameters, the frequency domain data volume Z is: L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI : Z MGI =ZL H(0,0,f1,f2)+αZ M H(f1,f2,f3,f4)+βZ H H(f1,f2,f3,f4)

[0119] Among them, Z MGI It represents the full-band wave impedance data volume, α represents the scale factor of the mid-frequency data volume, β represents the scale factor of the high-frequency data volume, which is used to adjust the amplitude of the components in each frequency band, H represents a filter with an amplitude of 1, and f1, f2, f3, and f4 represent the frequency parameters of the frequency domain filter.

[0120] Based on the aforementioned inventive concept, as shown in FIG9 , the present application further proposes a computer device 900, comprising a memory 910, a processor 920, and a computer program 930 stored in the memory 910 and executable on the processor 920. When the processor 920 executes the computer program 930, the aforementioned high-resolution wave impedance inversion method based on spectral fusion is implemented.

[0121] An embodiment of the present application further provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-mentioned high-resolution wave impedance inversion method based on spectrum fusion.

[0122] An embodiment of the present application further provides a computer program product, which includes a computer program. When the computer program is executed by a processor, it implements the above-mentioned high-resolution wave impedance inversion method based on spectrum fusion.

[0123] In summary, in the embodiment of the present application, seismic data, seismic wavelets, well logging data and horizon data are obtained; under the constraints of the horizon data and the spatial variation function, the low-frequency band of the well logging data is subjected to Kriging spatial interpolation processing to obtain the low-frequency band wave impedance data volume z L ; Establish a Cauchy norm constrained sparse pulse inversion intermediate frequency model, perform Cauchy norm constrained sparse pulse inversion on seismic data according to the selected seismic wavelet and Cauchy norm constrained sparse pulse inversion intermediate frequency model, and obtain the reflection coefficient; perform trace integration on the reflection coefficient to obtain the intermediate frequency relative wave impedance data volume z M Under the constraints of the layer data and spatial variation function, the high frequency band of the logging data is simulated sequentially by Gaussian random simulation to obtain the high frequency band wave impedance data volume z H ; respectively for the low frequency band wave impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L 、Z M 、ZH ; According to the selected frequency parameters and scale parameters, the frequency domain data volume Z L 、Z M 、Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ; For the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

[0124] By combining seismic inversion with geostatistical theory and adopting the idea of ​​"divide and conquer", the broadband inversion results are decoupled to obtain low-frequency, medium-frequency and high-frequency band results respectively, and finally integrated to obtain high-precision wave impedance inversion results, which improves the accuracy of traditional deterministic inversion and solves the problem that traditional random inversion requires a large number of iterations and has low computational efficiency.

[0125] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0126] The present application is described with reference to the flow chart and / or block diagram of the method, device (system), and computer program product according to the embodiment of the present application. It should be understood that each flow process and / or box in the flow chart and / or block diagram and the combination of the flow process and / or box in the flow chart and / or block diagram can be realized by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processing machine or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device produce a device for realizing the function specified in one flow chart flow or multiple flows and / or one box or multiple boxes of the block diagram.

[0127] These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce a product including an instruction device that implements the functions specified in one or more processes in the flowchart and / or one or more boxes in the block diagram.

[0128] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, so that the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one or more processes in the flowchart and / or one or more boxes in the block diagram.

[0129] The specific embodiments described above further illustrate the purpose, technical solutions and beneficial effects of the present application. It should be understood that the above description is only a specific embodiment of the present application and is not intended to limit the scope of protection of the present application. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application should be included in the scope of protection of the present application.

Claims

1. A high-resolution wave impedance inversion method based on spectrum fusion, It is characterized in that include: Acquire seismic data, seismic wavelets, well logging data and horizon data; Under the constraints of the layer data and the spatial variation function, the low-frequency band of the logging data is processed by Kriging spatial interpolation to obtain the low-frequency band wave impedance data volume z L ; A Cauchy norm constrained sparse pulse inversion intermediate frequency model is established. According to the selected seismic wavelet and the Cauchy norm constrained sparse pulse inversion intermediate frequency model, the seismic data is subjected to Cauchy norm constrained sparse pulse inversion to obtain the reflection coefficient. The reflection coefficient is integrated to obtain the intermediate frequency band relative wave impedance data volume z. M ; Under the constraints of the layer data and the spatial variation function, a sequential Gaussian random simulation is performed on the high-frequency band of the logging data to obtain the high-frequency band wave impedance data volume z H ; The low frequency band wave impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L , Z M , Z H ; According to the selected frequency parameters and scale parameters, the frequency domain data volume Z L , Z M , Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ; For the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

2. The method according to claim 1, It is characterized in that The low-frequency band of the logging data is processed by Kriging spatial interpolation to obtain the low-frequency band wave impedance data volume z L ,include: The ordinary Kriging method is used to perform Kriging spatial interpolation processing on the low-frequency band of the logging data to obtain the low-frequency band wave impedance data volume z L .

3. The method according to claim 1, It is characterized in that Establish a Cauchy norm constrained sparse pulse inversion intermediate frequency model, including: The objective function is constructed with the goal of minimizing the error between seismic data and synthetic seismic records, where seismic data refers to the actual observed seismic data and synthetic seismic records are obtained by convolution of seismic wavelets and reflection coefficients. Based on the objective function, the Cauchy norm constraint of the reflection coefficient is added to obtain a new objective function; The derivative of the new objective function with respect to the reflection coefficient is obtained, and the new objective function is set to 0, and the Cauchy norm constrained sparse pulse inversion intermediate frequency model is obtained.

4. The method according to claim 3, It is characterized in that The objective function is expressed as follows: minJ(r)=(Wr-d obs ) T (Wr-d obs ) Where minJ(r) represents the minimum error between seismic data and synthetic seismic records, that is, the minimum value of the objective function, d obs represents seismic data, W represents the seismic wavelet convolution matrix, and r represents the reflection coefficient.

5. The method according to claim 3, It is characterized in that The new objective function is expressed as follows: Where minJ(r) represents the minimum error between seismic data and synthetic seismic records, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, J s (r) is the seismic data constraint term, which controls the similarity between the synthetic seismic record and the seismic data. Cauchy (r) is the Cauchy norm constraint term of the reflection coefficient, which controls the sparsity of the reflection coefficient, μ is the weight coefficient of the sparse constraint, n is the number of sampling points of the seismic trace, i represents the sampling point sequence number of the seismic trace, and σ represents the variance of the reflection coefficient.

6. The method according to claim 3, It is characterized in that The Cauchy norm constrained sparse pulse inversion intermediate frequency model can be expressed as follows: (W T W+μQ)r=W T d obs Among them, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, μ is the weight coefficient of the sparse constraint, Q is a diagonal weighted matrix, and the diagonal elements i represents the sampling point number of the seismic trace, and σ represents the variance of the reflection coefficient.

7. The method according to claim 1, It is characterized in that Perform channel integration on the reflection coefficient to obtain the mid-frequency relative impedance data volume z M ,include: Perform trace integration on the reflection coefficient to obtain the relative wave impedance; According to the relative wave impedance, the mid-frequency band relative wave impedance data volume z is obtained. M .

8. The method according to claim 7, It is characterized in that Including integrating the reflection coefficient according to the following formula to obtain the relative wave impedance: z(t)=exp[2∫ 0 t r(τ)dτ] Among them, z(t) represents the relative wave impedance, r represents the reflection coefficient, t represents the time corresponding to a certain sampling point of the seismic trace, and τ represents the seismic trace sampling point before time point t.

9. The method according to claim 1, It is characterized in that Including according to the following formula, based on the selected frequency parameters and scale parameters, the frequency domain data volume Z L , Z M , Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI : Z MGI =Z L ·H(0,0,f 1 ,f 2 )+αZ M ·H(f 1 ,f 2 ,f 3 ,f 4 )+βZ H ·H(f 1 ,f 2 ,f 3 ,f 4 ) Among them, Z MGI represents the full-band wave impedance data volume, α represents the scale factor of the mid-frequency data volume, β represents the scale factor of the high-frequency data volume, which is used to adjust the amplitude of each frequency band component, H represents a filter with an amplitude of 1, and f 1 ,f 2 ,f 3 ,f 4 Represents the frequency parameters of the frequency domain filter.

10. A high-resolution wave impedance inversion device based on spectrum fusion, It is characterized in that include: A data acquisition module, used to acquire seismic data, seismic wavelets, well logging data and layer data; The Kriging spatial interpolation module is used to perform Kriging spatial interpolation processing on the low-frequency band of the logging data under the constraints of the layer data and the spatial variation function to obtain the low-frequency band wave impedance data volume z L ; The Cauchy norm constrained sparse pulse inversion module is used to establish the Cauchy norm constrained sparse pulse inversion intermediate frequency model. According to the selected seismic wavelet and the Cauchy norm constrained sparse pulse inversion intermediate frequency model, the seismic data is subjected to Cauchy norm constrained sparse pulse inversion to obtain the reflection coefficient; the reflection coefficient is integrated to obtain the intermediate frequency band relative wave impedance data volume z M ; The sequential Gaussian random simulation module is used to perform sequential Gaussian random simulation on the high frequency band of the logging data under the constraints of the layer data and the spatial variation function to obtain the high frequency band wave impedance data volume z H ; Fourier transform module is used to transform the low frequency impedance data volume z L , mid-frequency relative impedance data volume z M , high frequency band wave impedance data volume z H Perform Fourier transform to obtain the corresponding frequency domain data volume Z L , Z M , Z H ; The spectrum fusion module is used to fusion the frequency domain data volume Z according to the selected frequency parameters and scale parameters. L , Z M , Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI ; The inverse Fourier transform module is used to transform the full-band wave impedance data volume Z MGI Perform inverse Fourier transform to obtain high-resolution wave impedance inversion result z MGI .

11. The device according to claim 10, It is characterized in that The Kriging spatial interpolation module is specifically used for: The ordinary Kriging method is used to perform Kriging spatial interpolation processing on the low-frequency band of the logging data to obtain the low-frequency band wave impedance data volume z L .

12. The device according to claim 10, It is characterized in that The Cauchy norm constrained sparse pulse inversion module is specifically used for: The objective function is constructed with the goal of minimizing the error between seismic data and synthetic seismic records, where seismic data refers to the actual observed seismic data and synthetic seismic records are obtained by convolution of seismic wavelets and reflection coefficients. Based on the objective function, the Cauchy norm constraint of the reflection coefficient is added to obtain a new objective function; The derivative of the new objective function with respect to the reflection coefficient is obtained, and the new objective function is set to 0, and the Cauchy norm constrained sparse pulse inversion intermediate frequency model is obtained.

13. The device according to claim 12, It is characterized in that The objective function is expressed as follows: minJ(r)=(Wr-d obs ) T (Wr-d obs ) Where minJ(r) represents the minimum error between seismic data and synthetic seismic records, that is, the minimum value of the objective function, d obs represents seismic data, W represents the seismic wavelet convolution matrix, and r represents the reflection coefficient.

14. The device according to claim 12, It is characterized in that The new objective function is expressed as follows: Where minJ(r) represents the minimum error between seismic data and synthetic seismic records, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, J s (r) is the seismic data constraint term, which controls the similarity between the synthetic seismic record and the seismic data. Cauchy (r) is the Cauchy norm constraint term of the reflection coefficient, which controls the sparsity of the reflection coefficient, μ is the weight coefficient of the sparse constraint, n is the number of sampling points of the seismic trace, i represents the sampling point sequence number of the seismic trace, and σ represents the variance of the reflection coefficient.

15. The device according to claim 12, It is characterized in that The Cauchy norm constrained sparse pulse inversion intermediate frequency model is expressed as follows: (W T W+μQ)r=W T d obs Among them, d obs represents seismic data, W represents the seismic wavelet convolution matrix, r represents the reflection coefficient, μ is the weight coefficient of the sparse constraint, Q is a diagonal weighted matrix, and the diagonal elements i represents the sampling point number of the seismic trace, and σ represents the variance of the reflection coefficient.

16. The device according to claim 10, It is characterized in that The Cauchy norm constrained sparse pulse inversion module is specifically used for: Perform trace integration on the reflection coefficient to obtain the relative wave impedance; According to the relative wave impedance, the mid-frequency band relative wave impedance data volume z is obtained. M .

17. The device according to claim 16, It is characterized in that The Cauchy norm constrained sparse pulse inversion module is specifically used to: perform trace integration on the reflection coefficient according to the following formula to obtain the relative wave impedance: z(t) = exp[2∫ 0 t r(τ)dτ] Among them, z(t) represents the relative wave impedance, r represents the reflection coefficient, t represents the time corresponding to a certain sampling point of the seismic trace, and τ represents the seismic trace sampling point before time point t.

18. The device according to claim 10, It is characterized in that The spectrum fusion module is specifically used to: According to the following formula, based on the selected frequency parameters and scale parameters, the frequency domain data volume Z L , Z M , Z H Perform spectrum fusion to obtain the full-band wave impedance data volume Z MGI : Z MGI =Z L ·H(0,0,f 1 ,f 2 )+αZ M ·H(f 1 ,f 2 ,f 3 ,f 4 )+βZ H ·H(f 1 ,f 2 ,f 3 ,f 4 ) Among them, Z MGI represents the full-band wave impedance data volume, α represents the scale factor of the mid-frequency data volume, β represents the scale factor of the high-frequency data volume, which is used to adjust the amplitude of each frequency band component, H represents a filter with an amplitude of 1, and f 1 ,f 2 ,f 3 ,f 4 Represents the frequency parameters of the frequency domain filter.

19. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, It is characterized in that When the processor executes the computer program, the method according to any one of claims 1 to 9 is implemented.

20. A computer-readable storage medium, It is characterized in that The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the method according to any one of claims 1 to 9 is implemented.

21. A computer program product, It is characterized in that The computer program product comprises a computer program, and when the computer program is executed by a processor, the method according to any one of claims 1 to 9 is implemented.

Images