Method and system for optimally selecting carbon storage site based on multi-frequency band seismic data and equipment

Abstract

The invention belongs to the field of environmental monitoring, and in particular relates to a method for optimally selecting a carbon storage site based on multi-frequency band seismic data. The method comprises the steps of: performing seismic wavelet spread spectrum simulation based on three-dimensional post-stack seismic data to obtain spread spectrum simulated wavelets; building an isochronous stratigraphic framework model of a target horizon, and calculating the geometric structure and spatial distribution of a fault-karst; then performing waveform-indicated inversion to obtain a wave impedance inversion data volume, and obtaining a stable stratum wave impedance data volume through a virtual well cross-well wave impedance interpolation; calculating the difference between the stable stratum wave impedance data volume and the wave impedance inversion data volume to obtain an abnormal wave impedance data volume, then obtaining a fault-karst reservoir bed interpretation model, and determining the position of a carbon storage box.

Claims

1. A method for optimally selecting a carbon storage site based on multi-frequency band seismic data, comprising: using well logging instruments to generate original seismic data and well logging data at a known well site; configuring an electronic equipment, comprising at least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to realize the steps of: obtaining the original seismic data and well logging data of the known well site, and then obtaining near-well geological interpretation results and three-dimensional post-stack seismic data; based on the three-dimensional post-stack seismic data, obtaining the depth data of a maker layer of a target horizon and the isochronous three-dimensional distribution of the marker layer; based on the well logging data, eliminating the abnormal values and performing standardized pretreatment to obtain standardized well logging data; based on the three-dimensional post-stack seismic data, performing seismic wavelet spread spectrum simulation to obtain spread spectrum simulated wavelets; based on the spread spectrum simulated wavelets, building an isochronous stratigraphic framework model of a target horizon, and calculating the geometric structure and spatial distribution of the fault-karst; based on the three-dimensional post-stack seismic data and the spread spectrum simulated wavelets, performing multi-scale decomposition to obtain multi-scale Bimf components of the seismic data; based on the spread spectrum simulated wavelets and standardized well logging data, performing well-to-seismic calibration and ignoring small-scale Bimf components in the multi-scale Bimf components of the seismic data, thereby obtaining the time-depth conversion relationship; based on the well logging data, the frequency-boosted seismic data, and the time-depth conversion relationship, and on the basis of stratigraphic trend information provided by the isochronous stratigraphic framework model, performing waveform-indicated seismic wave impedance inversion to obtain a wave impedance inversion data volume; based on the three-dimensional post-stack seismic data, calculating a three-dimensional variance attribute volume, delimiting a sedimentary stable area, selecting a virtual-well site based on the sedimentary stable area, and obtaining a stable stratum wave impedance data volume through a virtual well cross-well wave impedance interpolation; calculating the difference between the stable stratum wave impedance data volume and the wave impedance inversion data volume to obtain the abnormal wave impedance data volume, and highlighting the influence of the fault zone on the wave impedance of the stable sedimentary strata; by removing areas lower than the average value in the three-dimensional variance attribute volume, retaining the abnormal wave impedance data in the spatial geometric contour of the fault-karst to obtain a fault-karst wave impedance data volume including the geometric structure and internal wave impedance characteristics of the fault-karst; comparing the near-well geological interpretation result with the fault-karst wave impedance data volume, delineating a characteristic value interval of a cave reservoir, a characteristic value interval of a crack fractured zone, and a characteristic value interval of surrounding rock, and obtaining a fault-karst reservoir bed interpretation model; and according to a fault-karst reservoir bed interpretation model, calculating migration and reservoir space after carbon dioxide injection, and finding a trap structure with a preset size threshold value, so that the migration and reservoir space where the injected carbon dioxide is stored in a high-density free phase under the conditions of being blocked by the tight carbonate rocks at the side of the fault zone and sealed by a top marl cover coat is selected as the position of the carbon storage box.

2. The method for optimally selecting a carbon storage site based on multi-frequency band seismic data of claim 1, wherein a method of obtaining spread spectrum simulated wavelet comprises: expressing a seismic record convolution model of the original seismic data in a frequency domain as:
s()=()() where, s() represents the seismic record frequency spectrum after Fourier transformation, () represents the seismic wavelet frequency spectrum after Fourier transformation, () represents the reflection coefficient frequency spectrum after Fourier transformation, and represents the angular frequency; transforming the expression of the seismic record convolution model in the frequency domain into the expression in the frequency domain of a seismic record linear system:
lns()=ln()+ln() where, lns() is the expression in the frequency domain of the seismic record linear system, ln() is the expression in the frequency domain of a seismic wavelet linear system, and ln() is the expression in the frequency domain of a reflection coefficient linear system; Fourier transformation is:
{tilde over (s)}(t)={tilde over ()}(t)+{tilde over ()}(t) where, {tilde over (s)}(t) represents a complex cepstrum sequence of a seismic record frequency spectrum, {tilde over ()}(t) represents a complex cepstrum sequence of a seismic wavelet frequency spectrum, {tilde over ()}(t) represents a complex cepstrum sequence of a reflection coefficient spectrum, and t represents the arrival time of seismic waves; separating the wavelet complex cepstrum sequence and the reflection coefficient complex cepstrum sequence in a complex cepstrum by a low-pass filter to extract wavelet amplitude; based on the wavelet amplitude, simulating a seismic wavelet amplitude spectrum by the least square method:
(f)=A(f)f.sup.e.sup.H(f) where, represents a constant greater than or equal to 0, (f) represents the seismic wavelet frequency spectrum, which is obtained from the Fourier transform of (t), and H(f) and A(f) are polynomials about the frequency f to be solved, which are used to fit the seismic wavelet amplitude spectrum; obtaining the maximum phase component and the minimum phase component of the wavelet based on the simulated seismic wavelet amplitude spectrum; letting the maximum phase component of the wavelet (t) be u(t) and the minimum phase component be v(t), then the wavelet (t) is:
(t)=u(t).Math.v(t) the complex cepstrum of the amplitude spectrum is expressed as:
2(t)=(t)+{tilde over (v)}(t)+(t)+{tilde over (v)}(t) where, the complex cepstrum (t) of the amplitude spectrum is symmetrically displayed on positive and negative axes of the complex cepstrum, {tilde over (v)}(t) is the complex cepstrum of the maximum phase function corresponding to the minimum phase component v(t) of the seismic wavelet, and (t) is the complex cepstrum of the minimum phase function corresponding to the maximum phase component of the seismic wavelet; and determining a set of mixed-phase wavelet sets with the same amplitude spectrum based on the complex cepstrum of the amplitude spectrum, adjusting Yu wavelet parameters, raising the main frequency by keeping the low frequency and expanding the high frequency on the premise of ensuring the integrity of the main frequency of the seismic wavelet, and raising the effective bandwidth to the preset bandwidth threshold to obtain the spread spectrum simulated wavelets.

3. The method for optimally selecting a carbon storage site based on multi-frequency band seismic data of claim 1, wherein a method of obtaining the Bimf components of the seismic data comprises: setting an initial pre-stack gather matrix;
r.sub.0(x,y)=U(x,y) where, x is the matrix row coordinate, y is the matrix column coordinate, U is the value of the post-stack three-dimensional seismic data, and r.sub.0 is the matrix element value; setting a variable l as the number of decomposed Bimf layers, and the initial value of l to 1, defining an initialization matrix variable as h, and obtaining a variable matrix h(x, y) as follows:
h(x,y)=r.sub.l-1(x,y) calculating the local maximum matrix and local minimum matrix about h(x, y); observing with an observation matrix of a preset size in the variable matrix, extracting the element value in the observation matrix, and obtaining the maximum value of the observation matrix, wherein the small matrix with the preset size of 44 can be used for observation; sliding the observation matrix until the center of the observation matrix traverses all matrix variables h, and assigning the measured maximum value of the observation matrix to the central element position of the corresponding observation matrix to obtain a local maximum matrix Max(x, y); obtaining the minimum value in the observation matrix by setting the observation matrix, and then obtaining a local minimum matrix Min(x, y); obtaining the maximum envelope surface and the minimum envelope surface by a spline interpolation method; the interpolation formula is:
N(i+m.sub.model,j+n.sub.model)=ABC where, A, B, and C represent an interpolation transition matrix: $\begin{array}{c}A=\left[\begin{array}{cccc}\left(1+m_{model}\right)& \left(m_{model}\right)& \left(1-m_{model}\right)& \left(2-m_{model}\right)\end{array}\right]\\ B=\left[\begin{array}{cccc}\mathrm{Max}\left(i-1,j-2\right)& \mathrm{Max}\left(i,j-2\right)& \mathrm{Max}\left(i+1,j-2\right)& \mathrm{Max}\left(i+2,j-2\right)\\ \mathrm{Max}\left(i-1,j-1\right)& \mathrm{Max}\left(i,j-1\right)& \mathrm{Max}\left(i+1,j-1\right)& \mathrm{Max}\left(i+2,j-1\right)\\ \mathrm{Max}\left(i-1,j\right)& \mathrm{Max}\left(i,j\right)& \mathrm{Max}\left(i+1,j\right)& \mathrm{Max}\left(i+2,j\right)\\ \mathrm{Max}\left(i-1,j+1\right)& \mathrm{Max}\left(i,j+1\right)& \mathrm{Max}\left(i+1,j+1\right)& \mathrm{Max}\left(i+2,j+1\right)\end{array}\right]\\ C=\left[\begin{array}{cccc}\left(1+n_{\mathrm{model}}\right)& \left(n_{model}\right)& \left(1-n_{model}\right)& \left(2-n_{model}\right)\end{array}\right]\end{array}$ where, m.sub.model and n.sub.model are preset values greater than 0 and less than 1, and 8 is the spread spectrum simulated wavelet; obtaining N.sub.Max as the maximum envelope surface, replacing the local maximum matrix Max(x, y) in a transition formula B with the local minimum matrix Min(x, y), and obtaining Min(x, y) as the minimum envelope surface; calculating the average value of the maximum envelope surface and the minimum envelope surface to obtain the intermediate transition quantity;
m.sub.l-1(x,y)=[N.sub.max(x,y)+Min(x,y)]/2 subtracting the intermediate transition quantity from the variable matrix h(x, y) to obtain the Bimf components of the J-th layer corresponding to the seismic data;
Bimf.sub.l(x,y)=h(x,y)m.sub.l-1(x,y) obtaining the margin r.sub.l-1 as an input variable for calculating the Bimf components of the next layer;
r.sub.l(x,y)=r.sub.l-1(x,y)Bimf.sub.l(x,y) where, r.sub.l represents the margin of the l-th layer, l1 represents the margin of the l1-th layer, and the margin of l1-th layer serves as an input variable for calculating the Bimf components of the l-th layer; calculating the Bimf components iteratively until the termination condition is met: ${.Math.}_{x=1}^M{.Math.}_{y=1}^N\left[\frac{{.Math.\mathrm{Bim}{f}_l\left(x,y\right)-{{\mathrm{Bimf}}_l}_{-1}\left(x,y\right).Math.}^2}{Bim{f}_{l-1}^2\left(x,y\right)}\right]r;$ and obtaining the multi-scale Bimf components of the seismic data, where, r is the preset termination threshold.

4. The method for optimally selecting a carbon storage site based on multi-frequency band seismic data of claim 1, wherein a method of obtaining the time-depth conversion relationship comprises: performing product operation based on a sonic time difference curve and a density curve in the well logging data of each known well site to obtain a wave impedance curve, and further calculating a reflection coefficient curve; constructing a Ricker wavelet on the basis of the main seismic frequency of a target interval, and performing convolution calculation of the Ricker wavelet and the reflection coefficient curve to obtain the synthetic seismic record; making the depth data of the maker layer at the wellbore of each drilling well position correspond to the isochronous three-dimensional distribution of the maker layer, calculating the correlation between the synthetic seismic record and the spread spectrum simulated wavelets of a seismic trace near the well, when the waveform correlation is higher than the first correlation threshold, and the preliminary well-to-seismic calibration is completed to obtain the preliminary time-depth conversion relationship between the well logging depth and the two-way travel time of seismic reflection waves; $T_d=T_{H_0}+2\underset{i_{{=1}^{}}}{\overset{z_{}}{.Math.}}{T}_{},.Math.H$ where, T.sub.H.sub.0 represents the two-way travel time of the seismic data corresponding to the depth of the sonic well logging marker layer; T.sub.i represents the sonic time difference; H represents the well logging curve data sampling interval; T.sub.d represents the two-way travel time of the seismic wave; adding the Bimf components of the seismic data step by step in descending order of scale to obtain updated effective seismic data; adding a Bimf component every time to obtain the updated effective seismic data, and calculating the second correlation between the synthetic seismic records after one time updating and the effective seismic data; as the Bimf components are gradually added to the effective seismic data, the second correlation shows an upward trend at first, and when the second correlation shows a downward trend, the updated synthetic seismic record at the peak of the second correlation and the updated earthquake are taken to calculate the second time-depth conversion relationship; and taking the second time-depth conversion relationship as the finally measured time-depth conversion relationship.

5. The method for optimally selecting a carbon storage site based on multi-frequency band seismic data of claim 1, wherein a method of obtaining a wave impedance inversion data volume comprises: based on the spread spectrum simulated wavelets, calculating the waveform correlation between the waveform of the seismic trace to be discriminated and the synthetic seismic record of the known well, and establishing an initial model according to the wave impedance curve corresponding to the well with the highest waveform correlation; using white noise to satisfy the law of Gaussian distribution, and expressing the wave impedance data in well logging data as:
Y.sub.sim=X.sub.sim+ where, Y.sub.sim represents the well logging wave impedance curve, X.sub.sim represents the actual wave impedance value of the underground stratum to be solved, and represents random noise; according to the central limit theorem, X.sub.simX.sub.p.sup.2 also satisfies the Gaussian distribution, determining the initial objective function as: $J_1=\frac{{.Math.X_{sim}-X_p.Math.}^2}{2{}^2}$ where J.sub.1 represents the function related to posterior information, X.sub.p represents the selection of a wave impedance curve of a sample well to a sample well based on the optimal number of samples, after matched filtering, posterior probability statistical distribution density is obtained, then the expected wave impedance value is calculated, and .sup.2 represents the covariance of white noise; based on the initial objective function, introducing prior information into the objective function through maximum posterior estimation to obtain a stable objective function as follows:
J(Z.sub.sim)=J.sub.1(Z.sub.sim)+.sup.2J.sub.2(Z.sub.sim) where, Z.sub.sim represents the characteristic parameter to be simulated, J.sub.2 represents the function related to the prior information of geological and well logging data, and represents the smoothing parameter used to coordinate the interaction between J.sub.1 and J.sub.2; and taking the stable objective function as the input of the initial model, sampling the posterior probability distribution by a Markov chain Monte Carlo method (MCMC) and Metropolis-Hastings sampling criterion, continuously optimizing the parameters of the initial model, selecting the solution when the objective function takes the maximum value as a value randomly realized, taking the average value randomly realized many times as the expected value output, and taking the expected value output as the wave impedance inversion data volume.

6. The method for optimally selecting a carbon storage site based on multi-frequency band seismic data of claim 1, wherein a method of delineating the sedimentary stable area specifically comprises: letting the data of each sampling point in the three-dimensional post-stack seismic data be S.sub.pqk.sub. where p represents the seismic gird wire size, q represents the seismic grid trace number, and k represents the sampling point serial number of seismic records sampled at 1 ms; calculating the mean square error of sampling point data in a preset sampling area: $Q_{{\mathrm{pqk}}_{}}=\underset{p-1}{\overset{p+1}{.Math.}}\underset{q-1}{\overset{q+1}{.Math.}}\underset{k_{}1}{\overset{k+1}{.Math.}}{\left(S_{\mathrm{pqk}}-1/9\underset{p-1}{\overset{p+1}{.Math.}}\underset{q-1}{\overset{q+1}{.Math.}}\underset{k_{}-1}{\overset{k+1}{.Math.}}{S}_{pq{k}_{}}\right)}^2$ translating the sampling areas vertically and horizontally, transversely calculating the data mean square errors of all the sampling areas to obtain a three-dimensional variance attribute volume; and slicing the three-dimensional variance attribute volume, obtaining the distribution characteristics of variance attribute data on the plane, and taking the area where the variance attribute value is lower than the total energy average value as the sedimentary stable area.

7. The method for optimally selecting a carbon storage site based on multi-frequency band seismic data of claim 1, wherein a method of obtaining the stable stratum wave impedance data volume comprises: delineating grids with preset sizes based on the sedimentary stable area, and taking each grid node as a virtual well site; letting the unknown underground wave impedance model parameter m be the g-dimensional space vector m={m.sub.1, . . . , m.sub.g}; through an earthquake acquisition process, obtaining the observation data D as a K-dimensional data vector d={d.sub.1, . . . , d.sub.k}; through a nonlinear function kernel G, establishing an unknown underground wave impedance model parameter, and establishing the relationship between the parameter and the K-dimensional data vector to obtain a forward model:
d=G(m)+n where, n={n.sub.1, . . . , n.sub.k} represents random noise independent of the underground wave impedance model parameter m, which obeys Gaussian distribution; constructing the inversion objective function based on the forward model:
minF[d,G(m)]=mind,G(m).sub.2.sup.2 where, F is the mean square error between the observation data d and the predicted data G(m); performing linearized solution on the inversion objective function; performing Taylor expansion on the forward model and omitting higher-order terms of more than quadratic to obtain a brief expression of prediction data:
G(m)=G(m)G(m.sup.0)=Am=A(mm.sup.0), where, m.sup.0 represents the initial model established according to prior information, A is a Jacobian matrix, and the element of A is first-order partial differential G.sub.i/m.sub.j; letting d.sub.0=G(m.sub.0)+.sub.2, d=dd.sub.0, .sub.2 represent random noise, then the iterative equation of the forward model is:
m.sup.t+1=m.sup.t+m, where, m.sup.t represents the forward model after iteration t times, m.sup.t+1 represents the forward model after iteration t+1 times, and d.sub.0 is the noise-adding predicted data; obtaining a pre-stack depth and an offset seismic profile d through a seismic acquisition process, assuming that the inversion depth domain model parameters have prior probability distribution P(M=m), according to Bayesian formula, the probability distribution is as follows: $P\left(M=m.Math.D=d^{}\right)=\frac{P\left(D=d^{}.Math.M=m\right)P\left(M=m\right)}{P\left(D=d^{}\right)}$ the probability distribution of adding random noise is follows: $P\left(D=d^{}.Math.M=m\right)=P\left(n\right)=\frac{1}{{\left(2\right)}^{\frac{n}{2}}{.Math.{\det C}_n.Math.}^{\frac{1}{2}}}\exp \left\{-{\left[d^{}-G\left(m\right)\right]}^T{C}_n^{-1}\left[d^{}-G\left(m\right)\right]\right\}$ where, C.sub.n is the noise covariance matrix, C.sub.n.sup.1 is the inversion operation of a noise covariance matrix, and det is the determinant of the corresponding matrix; n is noise and obeys Gaussian distribution, the average value is 0, C.sub.n is an equivalent diagonal matrix, and diagonal elements are the variance a containing noise of data center, then the probability distribution of adding the random noise is deformed as follows: $P\left(D=d^{}.Math.M=m\right)=.Math.\frac{1}{{\left(2\right)}^{\frac{n}{2}}{}_n}\exp \left\{-{}_n^2.Math.{\left[d^{}-G\left(m\right)\right]}^2\right\}$ for the prior probability distribution P(M=m) of the inversion depth domain model parameter, setting the initial model m.sup.0, there is m=m.sup.0+m, and the equivalent depth domain model probability distribution is: P(M=m)=P(m=m.sup.0+m)=P(m=mm.sup.0); letting the probability distribution of the equivalent depth domain model obey Gaussian distribution, and the probability distribution of the depth domain model with the random noise added is: $P\left(M=m\right)=\frac{1}{{\left(2\right)}^{\frac{n}{2}}{.Math.{\det C}_m.Math.}^{1/2}}\exp \left\{-{\left[m-m^0\right]}^T{C}_m^{-1}\left[m-m^0\right]\right\}$ assuming that the average value of the depth domain model parameters is equal to the initial model m.sup.0, and they do not affect each other, there is a simplified probability distribution of the depth domain model, which is as follows: $P\left(M=m\right)=\frac{1}{{\left(2\right)}^{\frac{n}{2}}{}_m}\exp \left\{-{}_m^2\left[m-m^0\right]\right\}$ where, .sub.m.sup.2 represents the variance value of the real model relative to the given initial model; a denominator P(D=d) in the probability distribution deformation added with the random noise added depends on data acquisition and processing, and is set as a constant, and maximizing posterior probability is equivalent to maximizing a numerator, which maximizes posterior probability distribution: $P\left(M=m.Math.D=d^{}\right)=.Math.\frac{1}{{\left(2\right)}^n{}_n{}_m}\exp \left\{-{}_n^{-2}.Math.{\left[d^{}-G\left(m\right)\right]}^2-{}_m^2.Math.{\left[m-m^0\right]}^2\right\}$ making the maximum equivalent formula of posterior probability distribution as:
F(m)=.sub.n.sup.2[dG(m)].sup.2+.sub.m.sup.2[mm.sup.0].sup.2.fwdarw.0 letting the partial derivative of the model parameter m equal to zero, and G(m)=G(m.sup.0)+Am, m=mm.sup.0, thereby obtaining the matrix equation: $\begin{array}{c}\frac{F\left(m\right)}{m}={}_n^{-2}\left[d^{}-G\left(m^0\right)-Am\right]\mathrm{IA}-{}_m^{2}Im=0\\ d^{}I{}_n^{-2}A-{}_m^{2}Im-G\left(m^0\right){}_n^{-2}\mathrm{IA}-A\mathrm{mI}{}_n^{-2}A=0\end{array}$ where, $A_{\mathrm{ij}}=\frac{G_i}{m_j},$ I is an identity matrix, m is a wave impedance iterative matrix, arranging to obtain the basic formula of random inversion as follows:
[dG[m.sup.0]]A=[AA.sup.T+.sub.n.sup.2.sub.m.sup.2I]m
m=[AA.sup.T+.sub.n.sup.2.sub.m.sup.2I].sup.1A.sup.T[dG(m.sup.0)] the iterative formula of the basic formula of random inversion is: m.sup.t+1=m.sup.t+m; where, G[m.sup.0] is the depth domain synthetic seismic record formed by the parametric model; updating the parameters of the unknown underground wave impedance model through continuous iteration until positive and negative shocks occur in m, where, m t+l is low frequency wave impedance inversion data; based on the correlation between the plane coordinates of the virtual well site and the plane coordinates of the seismic data, further determining the one-to-one correspondence between the virtual well and the low-frequency wave impedance inversion data, and then assigning the low frequency wave impedance inversion data beside the virtual well to the virtual well to obtain the wave impedance data of the virtual well site; and based on the wave impedance inversion data volume, extracting the wave impedance data of the virtual well site, setting the calculation area as the whole range of the work area in the geological framework according to the interpolation calculation of all virtual well wave impedance data in the time window range defined by the stratigraphic framework to obtain the stable stratigraphic wave impedance data.

8. A non-transitory computer-readable storage medium stores computer instructions which are used to be executed by a computer to realize the method for optimally selecting the carbon storage site based on multi-frequency band seismic data of claim 1.

9. A system for optimally selecting a carbon storage site based on multi-frequency band seismic data, comprising: well logging instruments configured to generate original seismic data and well logging data of a known well site; a data obtaining module configured for obtaining the original seismic data and well logging data of the known well site, then obtaining near-well geological interpretation results and three-dimensional post-stack seismic data, and obtaining the depth data of a maker layer of a target horizon and the isochronous three-dimensional distribution of the marker layer based on the three-dimensional post-stack seismic data; a pretreatment module configured for eliminating the abnormal values and performing standardized pretreatment to obtain standardized well logging data based on the well logging data; a spread spectrum simulation module configured for performing seismic wavelet spread spectrum simulation to obtain spread spectrum simulated wavelets based on the three-dimensional post-stack seismic data; an isochronous stratigraphic framework model obtaining module configured for building an isochronous stratigraphic framework model of a target horizon, and calculating the geometric structure and spatial distribution of the fault-karst based on the spread spectrum simulated wavelets; a Bimf component obtaining module configured for performing multi-scale decomposition to obtain multi-scale Bimf components of the seismic data based on the three-dimensional post-stack seismic data and the spread spectrum simulated wavelets; a time-depth conversion relationship obtaining module configured for performing well-to-seismic calibration and ignoring small-scale Bimf components in the multi-scale Bimf components of the seismic data to obtain the time-depth conversion relationship based on the spread spectrum simulated wavelets and standardized well logging data; a wave impedance inversion data volume obtaining module configured for performing waveform-indicated seismic wave impedance inversion to obtain a wave impedance inversion data volume based on the well logging data, the frequency-boosted seismic data and the time-depth conversion relationship, and on the basis of stratigraphic trend information provided by the isochronous stratigraphic framework model; a stable stratum wave impedance data volume obtaining module configured for based on the three-dimensional post-stack seismic data, calculating a three-dimensional variance attribute volume, delineating the a sedimentary stable area, selecting the virtual-well site based on the sedimentary stable area, and obtaining the stable stratum wave impedance data volume through the virtual well cross-well wave impedance interpolation; an abnormal wave impedance data volume obtaining module configured for calculating the difference between the stable stratum wave impedance data volume and the wave impedance inversion data volume to obtain the abnormal wave impedance data volume; a fault-karst wave impedance data volume obtaining module configured for retaining the abnormal wave impedance data in the spatial geometric contour of the fault-karst by removing areas lower than the average value in the three-dimensional variance attribute volume to obtain the fault-karst wave impedance data volume including the geometric structure and internal wave impedance characteristics of the fault-karst; a model interpretation module configured for comparing the near-well geological interpretation result with the fault-karst wave impedance data volume, delineating the characteristic value interval of a cave reservoir, the characteristic value interval of a fractured zone and the characteristic value interval of a surrounding rock, and obtaining a fault-karst reservoir bed interpretation model; and according to a fault-karst reservoir bed interpretation model, the migration and reservoir space is calculated after carbon dioxide injection, and a trap structure with the preset size threshold value is found, so that the migration and reservoir space where the injected carbon dioxide is stored in a high-density free phase under the conditions of being blocked by the tight carbonate rocks at the side of the fault zone and sealed by a top marl cover coat is selected as the position of the carbon storage box.

Description

BRIEF DESCRIPTION OF DRAWINGS

(1) FIG. 1A and FIG. 1B is a flow schematic diagram of a method for optimally selecting a carbon storage site based on multi-frequency band seismic data according to an embodiment of the present invention;

(2) FIG. 2 is a schematic diagram of original seismic data according to an embodiment of the present invention;

(3) FIG. 3 is a schematic diagram of an isochronous stratigraphic framework model according to an embodiment of the present invention;

(4) FIG. 4 is a schematic diagram of multi-scale Bimf components of seismic data when r is set to 0.2 according to an embodiment of the present invention;

(5) FIG. 5 is a schematic diagram of a virtual well site grid according to an embodiment of the present invention;

(6) FIG. 6 is an effect schematic diagram of a stable stratum wave impedance data volume according to an embodiment of the present invention;

(7) FIG. 7 is an effect schematic diagram of an abnormal wave impedance data volume according to an embodiment of the present invention; and

(8) FIG. 8 is a schematic diagram of the site selection result of carbon dioxide storage and the internal geometric structure of a storage box according to an embodiment of the present invention.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(9) The present invention provides a method for optimally selecting a carbon storage site based on multi-frequency band seismic data, which comprises: obtaining original seismic data and well well logging data of a known well site, and then obtaining near-well geological interpretation results and three-dimensional post-stack seismic data; based on the three-dimensional post-stack seismic data, obtaining the depth data of a maker layer of a target horizon and the isochronous three-dimensional distribution of the marker layer; based on the well logging data, eliminating abnormal values and performing standardized pretreatment to obtain standardized well logging data; based on the three-dimensional post-stack seismic data, performing seismic wavelet spread spectrum simulation to obtain spread spectrum simulated wavelets; based on the spread spectrum simulated wavelets, building an isochronous stratigraphic framework model of the target horizon, and calculating the geometric structure and spatial distribution of the fault-karst; based on the three-dimensional post-stack seismic data and the spread spectrum simulated wavelets, performing multi-scale decomposition to obtain multi-scale Bimf components of the seismic data; based on the spread spectrum simulated wavelets and the standardized well logging data, performing well-to-seismic calibration and ignoring small-scale Bimf components in the multi-scale Bimf components of the seismic data to obtain a time-depth conversion relationship; based on the well logging data, the frequency-boosted seismic data and the time-depth conversion relationship, and on the basis of stratigraphic trend information provided by the isochronous stratigraphic framework model, performing waveform-indicated seismic wave impedance inversion to obtain a wave impedance inversion data volume; based on the three-dimensional post-stack seismic data, calculating a three-dimensional variance attribute volume, delimiting a sedimentary stable area, selecting a virtual well site based on the sedimentary stable area, and obtaining a stable stratum wave impedance data volume through a virtual well cross-well wave impedance interpolation; calculating the difference between the stable stratum wave impedance data volume and the wave impedance inversion data volume to obtain the abnormal wave impedance data volume, and highlighting the influence of the fault zone on the wave impedance of the stable sedimentary strata; by removing areas lower than the average value in the three-dimensional variance attribute volume, retaining the abnormal wave impedance data in the spatial geometric contour of the fault-karst to obtain a fault-karst wave impedance data volume including the geometric structure and internal wave impedance characteristics of the fault-karst; comparing the near-well geological interpretation result with the fault-karst wave impedance data volume, delineating a characteristic value interval of a cave reservoir, a characteristic value interval of a fractured zone and a characteristic value interval of surrounding rock, and obtaining a fault-karst reservoir bed interpretation model; and according to the fault-karst reservoir bed interpretation model, calculating migration and reservoir space after carbon dioxide injection, and finding a trap structure with a preset size threshold value, so that the migration and reservoir space where the injected carbon dioxide is stored in a high-density free phase under the conditions of being blocked by tight carbonate rocks at the side of the fault zone and sealed by a top marl cover coat is selected as the position of the carbon storage box.

(10) The method for optimally selecting the carbon storage site based on multi-frequency band seismic data solves the problems by performing difference operation between a background stratum wave impedance model and the wave impedance inversion data volume that as the carbon storage boxes developed near the fault zones are complex in distribution and have various geometric structures, and it is difficult to accurately locate and identify the spatial geometric structures of the carbon storage boxes to assist carbon dioxide storage.

(11) In order to more clearly explain the method for optimally selecting the carbon storage site based on multi-frequency band seismic data, the steps in the embodiment of the present invention will be described in detail below with reference to FIG. 1.

(12) Each step of the method for optimally selecting the carbon storage site based on multi-frequency band seismic data in the first embodiment of the present invention is described in detail as follows: obtaining the original seismic data and well logging data of the known well site, and then obtaining the near-well geological interpretation results and the three-dimensional post-stack seismic data; based on the three-dimensional post-stack seismic data, obtaining the depth data of the maker layer of the target horizon and the isochronous three-dimensional distribution of the marker layer, wherein the original seismic data is as shown in FIG. 2; determining a top boundary marker layer of the Yingshan Formation by observing seismic waveform data, compared with the mudstone layer of the overlying Bachu Formation, the waveform reflection in the inner part of the carbonate rock is irregular and has no certain direction, the amplitude can be strong or weak, and the event can be long or short and has poor continuity; and there are non-systematic event reflection termination and bifurcation phenomena. According to the major axis of the interface between them, the top surface of Yingshan Formation is determined as the marker layer.

(13) The well logging data are nine conventional well logging curve data obtained by means of conventional well logging instruments within the detection depth range of 5500-5750 m at the wellbore of a work area; the well logging data comprises data of lithology characterization curves, data of well diameter curves, natural gamma curves, spontaneous potential curves and resistivity curves, data of deep lateral well logging curves, shallow lateral logging curves, micro lateral well logging curves, and physical property characterization curves, and data of compensated neutron curves, compensated acoustic wave curves and density curves. The sampling interval of the well logging curve data is 0.01 m. With the help of imaging well logging, well drilling, mud logging, cores, and other information, determinate lithology and physical property information of individual depth sections can be obtained, then various parameter threshold ranges of surrounding rocks, cave type reservoir beds, fracture type reservoir beds, and hole type reservoir beds are demarcated, and the near-well geological interpretation results can be obtained according to various parameter threshold ranges.

(14) According to the present embodiment, the three-dimensional post-stack seismic data with a work area of about 27 km.sup.2 is obtained by utilizing a three-dimensional seismic exploration method with the help of a seismic wave excitation source and a seismic signal detector, the double-way travel time of the signal record is 4 s, the time interval of the sampling point is 1 ms, and the detection depth exceeds 6000 m.

(15) At present, a 25 m25 m trace spacing three-dimensional seismic grid technology is widely used in the field of petroleum exploration. Generally, seismic wave reflection signals are received vertically at a sampling interval of 1 ms, and the total sampling time is within 6 s, so that the geological characteristics of different depth intervals are detected. The present embodiment is used for the detection of target bodies above 5 m, and has a high requirement for the main frequency of the seismic data, which should be within the range of 50-60 Hz. When the main frequency of the amplitude data volume of the seismic data in a fractured-cave body development interval is lower than 50 HZ, it is necessary to use a one-dimensional empirical mode decomposition algorithm based on seismic wavelet simulation to expand the frequency and reduce the noise of the seismic data for obtaining spread spectrum seismic data with a high resolution and a high signal-to-noise ratio.

(16) Based on the well logging data, eliminating the abnormal values and performing standardized pretreatment to obtain the standardized well logging data; the abnormal values are outliers, based on the obtained nine conventional well logging curve data, a statistical histogram is drawn for all wellbore data of each kind curves in the nine conventional logging curve data, the interval threshold is adjusted reasonably, the top 5% of data with the largest deviation from a median is removed, and data points in a box-selected interval, that is, the well logging data with abnormal values removed are retained.

(17) In the standardization process, based on the well logging data with abnormal values removed, single kind of data of all known well sites in the work area with the abnormal values removed are superimposed to draw a curve histogram, and standardized well logging data is obtained by integrating thresholds.

(18) The steps of abnormal value removal and standardization can eliminate the influence of the cross-well instruments.

(19) Based on the three-dimensional post-stack seismic data, seismic wavelet spread spectrum simulation is performed to obtain the spread spectrum simulated wavelets; and seismic wavelet simulation is a data processing method that improves the resolution of seismic signals by widening the effective frequency band on the premise of ensuring the relatively high fidelity of the processed seismic data.

(20) According to the present invention, a method of obtaining the spread spectrum simulated wavelets comprises: expressing a seismic record convolution model of the original seismic data in a frequency domain as:
s()=6()() where, s() represents the seismic record frequency spectrum after Fourier transformation, () represents the seismic wavelet frequency spectrum after Fourier transformation, () represents the reflection coefficient frequency spectrum after Fourier transformation, and represents the angular frequency; transforming the expression of the seismic record convolution model in the frequency domain into the expression in the frequency domain of a seismic record linear system:
lns()=ln()+ln() where, lns() is the expression in the frequency domain of the seismic record linear system, ln() is the expression in the frequency domain of a seismic wavelet linear system, and ln() is the expression in the frequency domain of a reflection coefficient linear system; Fourier transformation is:
{tilde over (s)}(t)={tilde over ()}(t)+{tilde over ()}(t) where, {tilde over (s)}(t) represents a complex cepstrum sequence of a seismic record frequency spectrum, {tilde over ()}(t) represents a complex cepstrum sequence of a seismic wavelet frequency spectrum, {tilde over ()}(t) represents a complex cepstrum sequence of a reflection coefficient spectrum, and t represents the arrival time of seismic waves; separating the wavelet complex cepstrum sequence from the reflection coefficient complex cepstrum sequence in a complex cepstrum by a low-pass filter to extract wavelet amplitude, wherein the wavelet amplitude is extracted by using the characteristic that the difference of smoothness between the wavelets and the reflection coefficient sequence is easy to distinguish in the complex cepstrum, the wavelet energy appears near the origin, while the reflection coefficient sequence is far from the origin; based on the wavelet amplitude, simulating a seismic wavelet amplitude spectrum by the least square method:
(f)=A(f)f.sup.e.sup.H(f) where, represents a constant greater than or equal to 0, (f) represents the seismic wavelet frequency spectrum, which is obtained from the Fourier transform of (t), and H(f) and A(f) are polynomials about the frequency f to be solved, which are used to fit the seismic wavelet amplitude spectrum; obtaining the maximum phase component and the minimum phase component of the wavelet based on the simulated seismic wavelet amplitude spectrum; letting the maximum phase component of the wavelet (t) be u(t) and the minimum phase component be (t), then the wavelet v(t) is:
(t)=u(t).Math.v(t) the complex cepstrum of the amplitude spectrum is expressed as:
2(t)=(t)+{tilde over (v)}(t)+(t)+{tilde over (v)}(t) where, the complex cepstrum (t) of the amplitude spectrum is symmetrically displayed on positive and negative axes of the complex cepstrum, {tilde over (v)}(t) is the complex cepstrum of the minimum phase function corresponding to the maximum phase component u(t) of the seismic wavelet, and {tilde over (v)}(t) is the complex cepstrum of the maximum phase function corresponding to the minimum phase component u(t) of the seismic wavelet; determining a set of mixed-phase wavelet sets with the same amplitude spectrum based on the complex cepstrum of the amplitude spectrum, adjusting Yu wavelet parameters, raising the main frequency by keeping the low frequency and expanding the high frequency on the premise of ensuring the integrity of the main frequency of the seismic wavelet, and raising the effective bandwidth to the preset bandwidth threshold to obtain the spread spectrum simulated wavelets. Usually, the effective bandwidth is set to 0-60 Hz.

(21) Based on the spread spectrum simulated wavelets, building the isochronous stratigraphic framework model of the target horizon, and calculating the geometric structure and spatial distribution of the carbon storage box, wherein the isochronous stratigraphic framework model is as shown in FIG. 3;

(22) based on the three-dimensional post-stack seismic data and the spread spectrum simulated wavelets, performing multi-scale decomposition to obtain multi-scale Bimf components of the seismic data; according to the present embodiment, the multi-scale decomposition of the seismic signal is realized by the step-by-step screening of the spread spectrum simulated wavelets, characteristic signals of different scales can be quickly extracted, noise are separated from useful signals, and the extraction of the internal trend information of the seismic signal is realized. Further, the ability of the seismic data to represent the details of the underground geological body is improved. The purpose of the present invention is to obtain a high-resolution result; high-resolution means that the seismic data with the high-frequency band is needed, but now the distribution range of the effective frequency band is 20-35 Hz, and there is a lot of noise in the high-frequency band, which affects the observation of the fault-karst; and by ignoring the small-scale Bimf components, the influence of noise in high resolution can be avoided, and the accuracy rate is improved. In the present invention, a method of obtaining the Bimf components of the seismic data comprises: slicing the three-dimensional post-stack seismic data along the direction of gather to obtain several pre-stack gather matrices; U(x, y), x=1, . . . , M; y=1, . . . , N; obtaining Bimf components by empirical mode decomposition; setting an initial pre-stack gather matrix:
r.sub.0(x,y)=U(x,y), where, x is the matrix row coordinate, y is the matrix column coordinate, U is the value of the post-stack three-dimensional seismic data, and r.sub.0 is the matrix element value; setting a variable l as the number of decomposed Bimf layers, and the initial value of l to 1, defining an initialization matrix variable as h, and obtaining a variable matrix h(x, y) as follows:
h(x,y)=r.sub.l-1(x,y); calculating the local maximum matrix and local minimum matrix about h(x, y); observing with an observation matrix of a preset size in the variable matrix h(x, y), extracting the element value in the observation matrix, and obtaining the maximum value of the observation matrix, wherein the small matrix with the preset size of 44 can be used for observation; sliding the observation matrix until the center of the observation matrix traverses all matrix variables h, and assigning the measured maximum value of the observation matrix to the central element position of the corresponding observation matrix to obtain a local maximum matrix Max(x, y); obtaining the minimum value in the observation matrix by setting the observation matrix, and then obtaining a local minimum matrix Min(x, y); obtaining the maximum envelope surface and the minimum envelope surface by a spline interpolation method; the interpolation formula is:
N(i+m,j+n)=ABC where, A, B, and C represent an interpolation transition matrix:

(23) $\begin{array}{c}A=\left[\begin{array}{cccc}\left(1+m\right)& \left(m\right)& \left(1-m\right)& \left(2-m\right)\end{array}\right]\\ B=\left[\begin{array}{cccc}\mathrm{Max}\left(i-1,j-2\right)& \mathrm{Max}\left(i,j-2\right)& \mathrm{Max}\left(i+1,j-2\right)& \mathrm{Max}\left(i+2,j-2\right)\\ \mathrm{Max}\left(i-1,j-1\right)& \mathrm{Max}\left(i,j-1\right)& \mathrm{Max}\left(i+1,j-1\right)& \mathrm{Max}\left(i+2,j-1\right)\\ \mathrm{Max}\left(i-1,j\right)& \mathrm{Max}\left(i,j\right)& \mathrm{Max}\left(i+1,j\right)& \mathrm{Max}\left(i+2,j\right)\\ \mathrm{Max}\left(i-1,j+1\right)& \mathrm{Max}\left(i,j+1\right)& \mathrm{Max}\left(i+1,j+1\right)& \mathrm{Max}\left(i+2,j+1\right)\end{array}\right]\\ C=[\begin{array}{cccc}\left(1+n\right)& \left(n\right)& \left(1-n\right)& \left(2-n\right)],\end{array}\end{array}$ where, u and v are preset values greater than 0 and less than 1, and 6 is the spread spectrum simulated wavelet; obtaining N.sub.Max as the maximum envelope surface, replacing the local maximum matrix Max(x, y) in a transition formula B with the local minimum matrix Min(x, y), and obtaining Min(x, y) as the minimum envelope surface; calculating the average value of the maximum envelope surface and the minimum envelope surface to obtain the intermediate transition quantity;
m.sub.l-1(x,y)=[N.sub.Max(x,y)+Min(x,y)]/2; subtracting the intermediate transition quantity from the variable matrix h(x, y) to obtain the Bimf components of the J-th layer corresponding to the seismic data:
Bimf.sub.l(x,y)=h(x,y)m.sub.l-1(x,y); obtaining the margin r.sub.l-1 as an input variable for calculating the Bimf components of the next layer;
r.sub.l(x,y)=r.sub.l-1(x,y)Bimf.sub.l(x,y); where, r.sub.l represents the margin of the l-th layer, l1 represents the margin of the l1-th layer, and the margin of the l1-th layer serves as an input variable of the Bimf component of the l-th layer; calculating the Bimf components iteratively until the termination condition is met:

(24) ${.Math.}_{x=1}^M{.Math.}_{y=1}^N\left[\frac{{.Math.\mathrm{Bim}{f}_l\left(x,y\right)-Bim{f}_{l-1}\left(x,y\right).Math.}^2}{Bim{f}_{l-1}^2\left(x,y\right)}\right]r;$
and obtaining the multi-scale Bimf components of the seismic data, where, r is the preset termination threshold; usually, the value of r is 0.2-0.3; the present embodiment takes 0.2 as an example, as shown in F 4, which can ensure the number and quality of Bimf and ensures that the bimf better reflects the details of waveform.

(25) The small-scale Bimf components contain a large amount of noise and detail information and edge information of the signal, while medium and large-scale Bimf components characterize the internal structural characteristics and trend characteristics of the signal.

(26) By ignoring the small-scale Bimf components of the seismic data, it is found that the correlation coefficients of synthetic seismic records and near-well earthquakes first increase and then decrease, and the corresponding reduced seismic data with the highest correlation coefficient is selected for subsequent calculation.

(27) After wavelet spread spectrum simulation, the effective frequency band of the seismic data is expanded, and the high frequency part is reasonably strengthened. The number of the events increases on the seismic waveform, which makes it easier to reflect the detailed changes of seismic reflection information, and improves the consistency of waveforms of the same reflection wave group in terms of amplitude, phase and frequency. The beaded reflection feature is particularly obvious on the seismic response of fractured-cave bodies, and the details of the internal shape of beads can be clearly displayed, which represents the seismic reflection of complex fractured-cave reservoir beds with different structural features and fillings, and is helpful for later fine geological interpretation.

(28) based on the spread spectrum simulated wavelets and the standardized well logging data, performing well-to-seismic calibration and ignoring small-scale Bimf components in the multi-scale Bimf components of the seismic data to obtain the time-depth conversion relationship; and based on the well logging data, the frequency-boosted seismic data and the time-depth conversion relationship, and on the basis of stratigraphic trend information provided by the isochronous stratigraphic framework model, performing waveform-indicated seismic wave impedance inversion to obtain a wave impedance inversion data volume.

(29) According to the present embodiment, a method of obtaining the time-depth conversion relationship comprises: performing product operation based on a sonic time difference curve and a density curve in the well logging data of each known well site to obtain a wave impedance curve, and further calculating a reflection coefficient curve; constructing a Ricker wavelet on the basis of the main seismic frequency of a target interval, preferably using the Ricker wavelet of 25 Hz, and performing convolution calculation of the Ricker wavelet and the reflection coefficient curve to obtain the synthetic seismic record; making the depth data of the maker layer at the wellbore of each drilling well position correspond to the isochronous three-dimensional distribution of the maker layer, calculating the correlation between the synthetic seismic record and the spread spectrum simulated wavelets of a seismic trace near the well, when the waveform correlation is higher than the first correlation threshold, the first correlation threshold is preferably greater than 85%, and the preliminary well-to-seismic calibration is completed to obtain the preliminary time-depth conversion relationship between the well logging depth and the two-way travel time of seismic reflection waves;
T.sub.d=T.sub.H.sub.0+2.sub.i=1.sup.nT.sub.i.Math.H, where, T.sub.H.sub.0 represents the two-way travel time of the seismic data corresponding to the depth of the sonic well logging marker layer, T, represents the sonic time difference, H represents the well logging curve data sampling interval; T.sub.d represents the two-way travel time of the seismic wave; adding the Bimf components of the seismic data step by step in descending order of scale to obtain updated effective seismic data; adding a Bimf component every time to obtain the updated effective seismic data, and calculating the second correlation between the synthetic seismic record after one time updating and the effective seismic data;

(30) As the Bimf components are gradually added to the effective seismic data, the second correlation shows an upward trend at first, and when the second correlation shows a downward trend, the updated synthetic seismic record at the peak of the second correlation and the updated earthquake are taken to calculate the second time-depth conversion relationship; and taking the second time-depth conversion relationship as the finally measured time-depth conversion relationship.

(31) According to the present embodiment, a method of obtaining a wave impedance inversion data volume comprises: based on the spread spectrum simulated wavelets, calculating the waveform correlation between the waveform of the seismic trace to be discriminated and the synthetic seismic record of the known well, and establishing an initial model according to the wave impedance curve corresponding to the well with the highest waveform correlation; using white noise to satisfy the law of Gaussian distribution, and expressing the wave impedance data in well logging data as:
Y.sub.sim=X.sub.sim+, where, Y.sub.sim represents the well logging wave impedance curve, X.sub.sim represents the actual wave impedance value of the underground stratum to be solved, and represents random noise;

(32) according to the central limit theorem, X.sub.simX.sub.p.sup.2 also satisfies Gaussian distribution, determining the initial objective function as:

(33) $J_1=\frac{{.Math.x_{sim}-x_p.Math.}^2}{2{}^2},$ where J.sub.1 represents the function related to posterior information, X.sub.p represents the selection of a wave impedance curve of a sample well to a sample well based on the optimal number of samples, after matched filtering, posterior probability statistical distribution density is obtained, then the expected wave impedance value is calculated, and represents the covariance of white noise; based on the initial objective function, introducing prior information into the objective function through maximum posterior estimation to obtain a stable objective function as follows:
J(Z.sub.sim)=J.sub.1(Z.sub.sim)+.sup.2J.sub.2(Z.sub.sim), where, Z.sub.sim represents the characteristic parameter to be simulated, J.sub.2 represents the function related to the prior information of geological and well logging data, and represents the smoothing parameter used to coordinate the interaction between J.sub.1 and J.sub.2; taking the stable objective function as the input of the initial model, sampling the posterior probability distribution by a Markov chain Monte Carlo method (MCMC) and Metropolis-Hastings sampling criterion, continuously optimizing the parameters of the initial model, selecting the solution when the objective function takes the maximum value as a value randomly realized, taking the average value randomly realized many times as the expected value output, and taking the expected value output as the wave impedance inversion data volume. based on the three-dimensional post-stack seismic data, calculating a three-dimensional variance attribute volume, delineating a sedimentary stable area, selecting a virtual well site based on the sedimentary stable area, and obtaining a stable stratum wave impedance data volume through a virtual well cross-well wave impedance interpolation; and the virtual well site division is shown in FIG. 5;

(34) According to the present embodiment, a method of delineating the sedimentary stable area specifically comprises: letting the data of each sampling point in the three-dimensional post-stack seismic data be S.sub.ijk, where p represents the seismic gird wire size, q represents the seismic grid trace number, and k represents the sampling point serial number of seismic records sampled at 1 ms; calculating the mean square error of sampling point data in a preset sampling area:
Q.sub.pqk=.sub.p1.sup.p+1.sub.k1.sup.k+1(S.sub.pqk 1/9.sub.p1.sup.p+1.sub.q1.sup.q+1S.sub.pqk).sup.2, translating the sampling areas vertically and horizontally, transversely calculating the data mean square errors of all the sampling areas to obtain a three-dimensional variance attribute volume; slicing the three-dimensional variance attribute volume, obtaining the distribution characteristics of variance attribute data on the plane, and taking the area where the variance attribute value is lower than the total energy average value as the sedimentary stable area.

(35) According to the present invention, a method of obtaining the stable stratum wave impedance data volume specifically comprises: delineating grids with preset sizes based on the sedimentary stable area, and taking each grid node as a virtual well site; according to the present embodiment, establishing 5050 m grids for the structural working drawing of the sedimentary stable area, and taking each grid node as the virtual well site.

(36) Letting the unknown underground wave impedance model parameter m be the g-dimensional space vector m={m.sub.1, . . . , m.sub.g}; through an earthquake acquisition process, obtaining the observation data D as a K-dimensional data vector d={d.sub.1, . . . , d.sub.k}; through a nonlinear function kernel G, establishing an unknown underground wave impedance model parameter, and establishing the relationship between the parameter and the K-dimensional data vector to obtain a forward model:
d=G(m)+n, where, n={n.sub.1, . . . , n.sub.r} represents random noise independent of the underground wave impedance model parameter m, which obeys Gaussian distribution; constructing the inversion objective function based on the forward model:
minF[d,G(m)]=mind,G(m).sub.2.sup.2, where, F is the mean square error between the observation data d and the predicted data G(m); there is a serious nonlinearity between the observation data d and the model parameter M to be inverted, so it is necessary to perform linearized solution on the inversion objective function, and the accuracy of the solution depends on the reliability of prior information. performing linearized solution on the inversion objective function; performing Taylor expansion on the forward model and omitting higher-order terms of more than quadratic to obtain a brief expression of prediction data:
G(m)=G(m)G(m.sup.0)=Am=A(mm.sup.0), where, m.sup.0 represents the initial model established according to prior information, A is a Jacobian matrix, and the element of A is first-order partial differential G.sub.i/m.sub.j; letting d.sub.0=G(m.sub.0)+,d=dd.sub.0, then the iterative equation of the forward model is:
m.sup.t+1=m.sup.t+m where, m.sup.k represents the forward model after iteration t times, m.sup.k+1 represents the forward model after iteration k+1 times, and d.sub.0 is the noise-adding predicted data; obviously, the linearized inversion method can obtain accurate results only when the initial model is relatively close to the real model, and from the point of mathematical calculation, the linear equation system tends to be ill-conditioned and has poor stability. When the prior information is insufficient, quantitative constraints can be performed on the model parameter to be inverted to control the size of the solution space so as to obtain stable and accurate inversion results; and this is the basic idea of a broadband constrained inversion method, and its essence is to obtain the optimal solution by using maximum likelihood distribution. obtaining a pre-stack depth and an offset seismic profile d through a seismic acquisition process, assuming that the inversion depth domain model parameters have prior probability distribution P(M=m), according to Bayesian formula, the probability distribution is:

(37) $P\left(M=m.Math.D=d\right)=\frac{P\left(D=d.Math.M=m\right)P\left(M=m\right)}{P\left(D=d\right)},$

(38) From the perspective of probability theory, the purpose of inversion is to obtain the maximum posterior probability density P(M=m|D=d), from the perspective of geophysics, when the depth migration is known, there are many depth domain model parameters that can form the depth domain migration profile through the forward modeling process, the solution space is very large, and the purpose of inversion is to find the solution closest to the real earth model; the probability distribution of adding random noise is follows:

(39) $P\left(D=d.Math.M=m\right)=P\left(n\right)=\frac{1}{{\left(2\right)}^{\frac{n}{2}}{.Math.{\det C}_n.Math.}^{\frac{1}{2}}}\exp \left\{-{\left[d-G\left(m\right)\right]}^T{C}_n^{-1}\left[d-G\left(m\right)\right]\right\},$ where, C.sub.n is the noise covariance matrix, C.sub.n.sup.1 is the inversion operation of a noise covariance matrix, and det is the determinant of the corresponding matrix; n is noise and obeys Gaussian distribution, the average value is 0, C.sub.n is an equivalent diagonal matrix, and diagonal elements are the variance a containing noise of data center, then the probability distribution of adding the random noise is deformed as follows:

(40) $P\left(D=d.Math.M=m\right)=.Math.\frac{1}{{\left(2\right)}^{\frac{n}{2}}{}_n}\exp \left\{-{}_n^2.Math.{\left[d-G\left(m\right)\right]}^2\right\},$

(41) For the prior probability distribution P(M=m) of the inversion depth domain model parameter, setting the initial model m.sup.0, there is m=m.sup.0+m, and the equivalent depth domain model probability distribution is: P(M=m)=P(m=m.sup.0+m)=P(m=mm.sup.0); letting the probability distribution of the equivalent depth domain model obey Gaussian distribution, and the probability distribution of the depth domain model with the random noise added is:

(42) $P\left(M-m\right)=\frac{1}{{\left(2\right)}^{\frac{n}{2}}{.Math.{\det C}_m.Math.}^{1/2}}\exp \left\{-{\left[m-m^0\right]}^T{C}_m^{-1}\left[m-m^0\right]\right\},$ assuming that the average value of the depth domain model parameters is equal to the initial model m.sup.0, and they do not affect each other, there is simplified probability distribution of the depth domain model, which is as follows:

(43) 0$P\left(M-m\right)=\frac{1}{{\left(2\right)}^{\frac{n}{2}}{}_m}\exp \left\{-{}_m^2\left[m-m^0\right]\right\},$ where, .sub.m represents the variance value of the real model relative to the given initial model; a denominator P(D=d) in the probability distribution deformation with the random noise added depends on data acquisition and processing, and is set as a constant, and maximizing posterior probability is equivalent to maximizing a numerator, which maximizes posterior probability distribution:

(44) $P\left(M=m.Math.D=d\right)=.Math.\frac{1}{{\left(2\right)}^n{}_n{}_m}\exp \left\{-{}_n^{-2}.Math.{\left[d-G\left(m\right)\right]}^2-{}_m^2.Math.{\left[m-m^0\right]}^2\right\},$ making the maximum equivalent formula of posterior probability distribution as:
F(m)=.sub.n.sup.2[dG(m)].sup.2+.sub.m.sup.2[mm.sup.0].sup.2.fwdarw.0, letting the partial derivative of the model parameter m equal to zero, and G(m)=G(m.sup.0)+Am, m=mm.sup.0, thereby obtaining the matrix equation:

(45) $$\begin{gathered} \frac{F\left(m\right)}{m}={}_n^{-2}\left[d-G\left(m^0\right)-Am\right]\mathrm{IA}-{}_m^{2}Im=0, \\ \mathrm{dI}{}_n^{-2}A-{}_m^{2}Im-G\left(m^0\right){}_n^{-2}\mathrm{IA}-A\mathrm{mI}{}_n^{-2}A=0, \end{gathered}$$ where,

(46) $A_{\mathrm{ij}}=\frac{G_i}{m_j},$ I is an identity matrix, m is a wave impedance iterative matrix, arranging to obtain the basic formula of random inversion as follows:
[dG[m.sup.0]]A=[AA.sup.T+.sub.n.sup.2.sub.m.sup.2I]m
m=[AA.sup.T+.sub.n.sup.2m.sup.2I].sup.1A.sup.T[dG(m.sup.0)]; the iterative formula of the basic formula of random inversion is: m.sup.k+1=m.sup.k+m; where, G[m ] is the depth domain synthetic seismic record formed by the parametric model, and d is observation data; updating the parameters of the unknown underground wave impedance model through continuous iteration until positive and negative shocks occur in m, where, m is low frequency wave impedance inversion data; based on the correlation between the plane coordinates of the virtual well site and the plane coordinates of the seismic data, further determining the one-to-one correspondence between the virtual well and the low-frequency wave impedance inversion data, and then assigning the low frequency wave impedance inversion data beside the virtual well to the virtual well to obtain the wave impedance data of the virtual well site; based on the wave impedance inversion data volume, extracting the wave impedance data of the virtual well site, setting the calculation area as the whole range of the work area in the geological framework according to the interpolation calculation of all virtual well wave impedance data in the time window range defined by the stratigraphic framework to obtain the stable stratigraphic wave impedance data, and the effect of the stable stratum wave impedance data is as shown in FIG. 6; calculating the difference between the stable stratum wave impedance data volume and the wave impedance inversion data volume to obtain the abnormal wave impedance data volume, and the effect of the abnormal wave impedance data volume is shown in FIG. 7; by removing areas lower than the average value in the three-dimensional variance attribute volume, retaining the abnormal wave impedance data in the spatial geometric contour of the fault-karst to obtain the fault-karst wave impedance data volume including the geometric structure and internal wave impedance characteristics of the fault-karst; comparing the near-well geological interpretation result with the fault-karst wave impedance data volume, delineating the characteristic value interval of a cave reservoir bed, the characteristic value interval of the crack fractured zone, and the characteristic value interval of the surrounding rock, and obtaining the fault-karst reservoir bed interpretation model; and according to the fault-karst bed interpretation model, evaluating the lithology and thickness of the caprock and selecting the position of the carbon storage box.

(47) Well logging interpretation results, such as lithologic interpretation in the depth section where drilling is blocked, including fracture reservoir beds, cave reservoir beds, transition zones, and surrounding rocks are input; and according to the time-depth conversion relationship among the well site, the well logging data and the seismic data, the drilling trajectory and the well logging interpretation results are projected on the profile of the fault-karst wave impedance data.

(48) The well logging interpretation results are compared with the eigenvalue energy data of a structure-eigenvalue model of the fault-karst, an area with the eigenvalue greater than 0.82 is defined as a crack reservoir bed, an area with the eigenvalue between 0.63 and 0.82 is defined as the cave reservoir bed, an area with the eigenvalue between 0.31 and 0.63 is defined as a transition zone, and an area with the eigenvalue less than 0.31 is defined as the surrounding rock, which is used as the final fault-karst reservoir bed interpretation model; and

(49) according to the fault-karst reservoir bed interpretation model, evaluating the lithology and thickness of the caprock and selecting the position of the carbon storage box. The site selection result of carbon dioxide storage and the internal geometric structure of the storage box are as shown in FIG. 8.

(50) Although the steps are described in the above-mentioned sequence in the above-mentioned embodiment, those skilled in the art can understand that in order to achieve the effect of the present embodiment, different steps do not need to be performed in such an order, but can be performed simultaneously (in parallel) or in reverse order, and these simple changes are within the scope of protection of the present invention.

(51) A second embodiment of the present invention provides a system for optimally selecting the carbon storage site based on multi-frequency band seismic data, which comprises: a data obtaining module configured for obtaining the original seismic data and well logging data of the known well site, then obtaining near-well geological interpretation results and three-dimensional post-stack seismic data, and obtaining the depth data of the maker layer of the target horizon and the isochronous three-dimensional distribution of the marker layer based on the three-dimensional post-stack seismic data; a pretreatment module configured for eliminating the abnormal values and performing standardized pretreatment to obtain standardized well logging data based on the well logging data; a spread spectrum simulation module configured for performing seismic wavelet spread spectrum simulation to obtain spread spectrum simulated wavelets based on the three-dimensional post-stack seismic data; an isochronous stratigraphic framework model obtaining module configured for building the isochronous stratigraphic framework model of the target horizon, and calculating the geometric structure and spatial distribution of the fault-karst based on the spread spectrum simulated wavelets; a Bimf component obtaining module configured for performing multi-scale decomposition to obtain multi-scale Bimf components of the seismic data based on the three-dimensional post-stack seismic data and the spread spectrum simulated wavelets; a time-depth conversion relationship obtaining module configured for performing well-to-seismic calibration and ignoring small-scale Bimf components in the multi-scale Bimf components of the seismic data to obtain the time-depth conversion relationship based on the spread spectrum simulated wavelets and standardized well logging data; a wave impedance inversion data volume obtaining module configured for performing waveform-indicated seismic wave impedance inversion to obtain a wave impedance inversion data volume based on the well logging data, the frequency-boosted seismic data and the time-depth conversion relationship, and on the basis of stratigraphic trend information provided by the isochronous stratigraphic framework model; a stable stratum wave impedance data volume obtaining module configured for based on the three-dimensional post-stack seismic data, calculating a three-dimensional variance attribute volume, delineating the sedimentary stable area, selecting the virtual well site based on the sedimentary stable area, and obtaining the stable stratum wave impedance data volume through the virtual well cross-well wave impedance interpolation; an abnormal wave impedance data volume obtaining module configured for calculating the difference between the stable stratum wave impedance data volume and the wave impedance inversion data volume to obtain the abnormal wave impedance data volume; a fault-karst wave impedance data volume obtaining module configured for retaining the abnormal wave impedance data in the spatial geometric contour of the fault-karst by removing areas lower than the average value in the three-dimensional variance attribute volume to obtain the fault-karst wave impedance data volume including the geometric structure and internal wave impedance characteristics of the fault-karst; and a model interpretation module configured for comparing the near-well geological interpretation result with the fault-karst wave impedance data volume, delineating the characteristic value interval of the cave reservoir, the characteristic value interval of the fractured zone and the characteristic value interval of the surrounding rock, and obtaining the fault-karst reservoir bed interpretation model.

(52) According to the fault-karst reservoir bed interpretation model, the migration and reservoir space is calculated after carbon dioxide injection, and a trap structure with the preset size threshold is found, so that the migration and reservoir space where the injected carbon dioxide is stored in a high-density free phase under the conditions of being blocked by tight carbonate rocks at the side of the fault zone and sealed by the top marl cover coat is selected as the position of a carbon storage box. A person skilled in the art can clearly understand that for convenience and conciseness of description, specific working processes and related explanations of the foregoing described system can refer to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

(53) It should be noted that the system for optimally selecting the carbon storage site based on multi-frequency band seismic data provided by the above embodiment is only exemplified by the division of the above-mentioned functional modules, in practical application, above functions may be allocated to different functional modules according to needs, that is, the modules or steps in the embodiment of the present invention can be decomposed or combined again; and for example, the modules in the above embodiments can be combined into one module or further split into multiple sub-modules to complete all or part of the above-mentioned functions. The names of the modules and steps involved in the embodiments of the present invention are only for distinguishing each module or step, and are not regarded as improper restrictions on the present invention.

(54) Electronic equipment of the third embodiment of the present invention comprises at least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to realize the method for optionally selecting the carbon storage site based on multi-frequency band seismic data.

(55) A fourth embodiment of the present invention provides a computer-readable storage medium; and the computer-readable storage medium stores the computer instructions which are used to be executed by a computer to realize the method for optionally selecting the carbon storage site based on multi-frequency band seismic data.