CROSS-REFERENCE TO RELATED APPLICATIONS

The application claims priority to Chinese patent application No. 2021110663557, filed on Sep. 13, 2021, the entire contents of which are incorporated herein by reference.

TECHNICAL FIELD

The present invention belongs to the field of treatment for data identification and recording carriers, and specifically relates to a method and system for evaluating the filling characteristics of a deep paleokarst reservoir through well-to-seismic integration.

BACKGROUND

Deep strata (the buried depth over 4,500 m) of oil and gas bearing basins are rich in oil and gas resources and have become the important field of industrial exploration of petroleum. Carbonate reservoirs account for 52% of proven oil reserves and 60% of production in the world, which are the main target of deep oil and gas exploration.

Due to the complex lithology, changeable diagenesis and random fractures of carbonate rocks, a paleokarst reservoir is more highly heterogeneous. Sparse logging curves have high longitudinal resolution ratio and only can reflect information of the strata with limited depth (less than 10 m generally) near a wellbore, and the paleokarst reservoir that is far away from the wellbore still has great change. Therefore, detection technologies and experience that are gained from shallow clastic rocks and porous carbonate reservoirs cannot be directly used in the detection for a deep paleokarst reservoir. However, the above descriptions are often neglected by those skilled in the art, which leads to the poorer evaluation effect of the deep paleokarst reservoir.

The current evaluation for the filling characteristics of the paleokarst reservoir has the following problems: the paleokarst cave type reservoir is buried deeply, and the absorption and attenuation of seismic wave energy are strong, which lead to the significant decrease of the dominant frequency of seismic data (less than 25 Hz), greatly reduced resolution ratio (10 m generally) and too low signal-to-noise ratio in a deep window. However, drilling, logging and rock-core samples indicate that the height of most caves is less than 5 m, so the existing resolution ratio of 3D (Three-Dimensional) seismic data is difficult to achieve the accurate identification effect. In addition, compared with the logging curves, seismic attributes contain lower content of lithological and physical information, so that the filling characteristics of the paleokarst reservoir cannot be effectively identified.

The key problems of the research on the paleokarst reservoir are: how to accurately understand the structure of a deeply buried paleokarst system, how to effectively identify the position of the paleokarst reservoir by utilizing geophysical data, and how to characterize the internal filling characteristics of a paleokarst cave system with high precision.

SUMMARY

In order to solve the technical problems of evaluation for the filling characteristics of a deep paleokarst reservoir at present, namely, the problems that the inversion result is directly controlled by inter-well interpolation in the high-frequency band, the lateral prediction is highly random, and the prediction for the deep paleokarst reservoir with fast change of lateral filling characteristics is poor in precision, which are caused by that in the existing methods for analyzing the filling characteristics of a cave type reservoir in petroleum exploration, only geological data is regarded as the constraint for a random simulation conclusion, the contribution of seismic information to the optimization for logging samples is not considered, and data is too complex as the reservoir is described by adopting all logging characteristic parameters, the present invention provides a method for evaluating the filling characteristics of the deep paleokarst reservoir through well-to-seismic integration, which comprises the following steps:

Step S100, acquiring original geophysical logging data: acquiring original logging data of each sample well by logging equipment, comprising: measuring the self-potential SP of each sample well by a measuring electrode, measuring the natural gamma ray GR of each sample well by a natural gamma ray downhole device and a natural gamma ray ground instrument, and acquiring the caliper CAL of each sample well by a caliper arm; and acquiring data of a resistivity curve: deep lateral logging RLLD, shallow lateral logging RLLS and micro-lateral logging RLLM, and data of a physical property characterization curve: compensated neutrons CNL, a compensated acoustic curve AC and a density curve DEN by traditional logging equipment; and

acquiring definite lithological information and physical information of an individual depth segment based on imaging logging information, drilling information, mud logging information and rock-core information, so as to determine depth data of a marker bed at a target horizon;

Step S200, acquiring seismic data: acquiring data of an original seismic wave reflected signal by a seismic wave excitation device and a receiving device and acquiring isochronous 3D (Three-Dimensional) distribution of the marker bed at the target horizon according to the waveform of the data of the original seismic wave reflected signal;

Step S300, preprocessing original geophysical logging data: drawing data of logging curves based on all original logging data of the sample wells and carrying out abnormal value processing and standardization processing, so as to obtain data of standardized logging curves;

Step S400, preprocessing the seismic data: obtaining a high-precision 3D seismic amplitude data body by mixed-phase wavelet estimation and maximum posteriori deconvolution and enhancing diffusion filtering based on the data of the seismic wave reflected signal;

Step S500, carrying out well-to-seismic calibration and selecting the characteristic parameters: acquiring a wave impedance curve of each sample well based on the compensated acoustic curve AC and the density curve DEN in the data of the standardized logging curve, so as to calculate a reflection coefficient curve; acquiring the preferred frequency of the ricker wavelet, so that the preferred frequency is consistent with the dominant frequency of the high-precision 3D seismic amplitude data body; carrying out convolution operation on the ricker wavelet and the reflection coefficient curve, so as to obtain a synthetic seismic record; comparing the depth data of the marker bed at the target horizon and the isochronous 3D distribution of the marker bed at the target horizon for well-to-seismic calibration; calculating the correlation between the synthetic seismic record and the waveform of a near-well seismic channel; and when the correlation is more than or equal to a preset first threshold, judging that the result of well-to-seismic calibration is qualified, so as to obtain the time-depth transforming relationship between the data of the logging curve and the seismic record and the characteristic parameters sensitive to the reservoir,

wherein a method for obtaining the characteristic parameters sensitive to the reservoir comprises:

drawing a columnar statistical graph of responses of the logging parameters, which are generated to different geological bodies near the wells; and when the numerical value of data of a certain standardized logging curve can be used for distinguishing data points that are more than a second threshold preset by different logging interpretation conclusions, selecting the data of the standardized logging curve as the characteristic parameters sensitive to the reservoir, wherein the characteristic parameters sensitive to the reservoir at least comprise the wave impedance IMP and can also comprise one or more of the caliper CAL, the natural gamma ray GR, the self-potential SP, the data of the resistivity curve: the deep lateral logging RLLD, the shallow lateral logging RLLS and the micro-lateral logging RLLM, and the data of the physical property characterization curve: the compensated neutrons CNL, the compensated acoustic curve AC and the density curve DEN;

Step S600, acquiring a first group of PCA (Principal Component Analysis) data: optimizing the data of the standardized logging curves from the characteristic parameters sensitive to the reservoir, wherein the number of pieces of the data is preset; and analyzing the characteristic parameters sensitive to the reservoir by a PCA method, so as to obtain the first group of PCA data sensitive to the reservoir;

Step S700, constructing an isochronous framework model: constructing the isochronous framework model based on the law of sedimentary strata reflected by the high-precision 3D seismic amplitude data body and the time-depth transforming relationship;

Step S800, simulating the parameters of an inter-well reservoir: based on the first group of PCA data of each sample well, determining the parameter of the optimal sample number, which can reflect the overall geological condition; selecting the first group of PCA data of the sample wells to construct an initial model, wherein the first group of PCA data has the highest correlation of the seismic waveforms, and the number of pieces of first group of PCA data is equal to the parameter of the optimal sample number; constantly correcting the parameters of the initial model; and outputting a first group of high-precision characteristic value simulation result data bodies, wherein the first group of high-precision characteristic value simulation result data bodies is data bodies corresponding to the first group of PCA data one to one;

Step S900, characterizing the boundary of an inter-well karst cave system: carrying out cross analysis based on the first group of PCA data of the sample wells, so as to obtain a distribution diagram of a first group of discrete data points of various geological bodies;

based on the lithological information and the physical information, acquiring the interpretation conclusion of the distribution diagram of the first group of discrete data points of various geological bodies, dividing a karst cave and a surrounding rock to construct a cross plot of the first group of PCA data, and acquiring the threshold of the first group of PCA data, which is required for division of the karst cave; and

based on the threshold of the first group of PCA data, carrying out cross analysis on the first group of high-precision characteristic value simulation result data bodies, so as to obtain the characteristics of the 3D spatial morphologies of the karst cave;

Step S1000, characterizing the lithological and physical boundaries of internal filling of the cave: analyzing the characteristic parameters sensitive to the reservoir, which correspond to data points of the karst cave, by the PCA method, so as to obtain a second group of PCA data sensitive to the response of fillings;

based on the second group of PCA data, generating a second group of high-precision characteristic value simulation result data bodies by the method described in Step S800;

based on the second group of PCA data sensitive to the response of the fillings, carrying out cross analysis to obtain a distribution diagram of a second group of discrete data points of various geological bodies;

based on the lithological information and the physical information, acquiring the interpretation conclusion of the second group of discrete data points of various geological bodies, and distinguishing the types of the fillings to construct a cross plot of the second group of PCA data, so as to obtain the threshold of the second group of PCA data, which is required for division of the types of the internal fillings of the karst cave; and

based on the threshold of the second group of PCA data, carrying out cross analysis on the second group of high-precision characteristic value simulation result data bodies, so as to obtain the characteristics of the 3D spatial morphologies of different types of internal fillings of the karst cave; and

Step S1100, describing the structure and the filling of the paleokarst cave: based on the characteristics of the 3D spatial morphologies of the karst cave and the characteristics of the 3D spatial morphologies of different types of internal fillings of the karst cave, engraving the spatial distribution of the paleokarst cave and the development characteristics of the internal filling adopting different types of internal fillings by adopting a lithology shielding technology and a 3D engraving technology.

In some preferred implementation manners, the step of measuring the self-potential SP of each sample well by the measuring electrode and measuring the natural gamma ray GR of each sample well by the natural gamma ray downhole device and the natural gamma ray ground instrument specifically comprises:

measuring the self-potential SP of each sample well by the measuring electrode:

arranging a measuring electrode Qelectrodeon the ground and arranging a measuring electrode Zelectrodeunder the well by a cable;

lifting the measuring electrode Zelectrodealong the well axis to measure the change of the self-potential with the well depth, wherein

a calculation method of the value of the self-potential is shown as follows:

EA=Ed-Eda=(Kd-Kda)⁢lg⁢RmfRw=Kdiffusion⁢lg⁢RmfRw

wherein EArepresents the total self-potential; Kdrepresents the coefficient of diffusion potential; Kdarepresents the coefficient of diffusion and adsorption potential; Rmfrepresents the value of resistivity of mud filtrate; and RWrepresents the value of resistivity of stratum water; and

measuring the natural gamma ray GR of each sample well by the natural gamma ray downhole device and the natural gamma ray ground instrument:

wherein the natural gamma ray downhole device comprises a detector, an amplifier and a high-voltage power supply;

acquiring the natural gamma ray by the detector, transforming the natural gamma ray into an electric pulse signal and amplifying the electric pulse signal by the amplifier; and

transforming the number of electric pulses formed in each minute into the potential difference and recording the potential difference by the ground instrument.

In some preferred implementation manners, the PCA method specifically comprises:

Step B100, calculating the average value of the samples:

x_=1N⁢∑i=1Nxi

wherein N represents the total sample number of the curve; i represents the parameter of the sample number, xirepresents the ithdata of the characteristic parameters of the sample, andxrepresents the average value of the samples;

Step B200, normalizing each characteristic parameter curve to obtain the normalized data xSof the characteristic parameters of the sample:

xS=x-x_N

wherein x represents the data of sampling points of the characteristic parameters of the sample;

Step B300, based on the normalized data xSof the characteristic parameters of the samples, calculating the covariance Cov(X, Y) of the characteristic parameters of the samples:

Cov(X,Y)=E[(X-E⁡(X))⁢(Y-E⁡(Y))]=1N-1⁢∑c=1N(xc-x¯)⁢(yc-y¯)

wherein xcrepresents the cthdata of the characteristic parameters of an X sample, and ycrepresents the cthdata of the characteristic parameters of a Y sample;

Step B400, based on the covariance of the characteristic parameters of the samples, constructing a covariance matrix of the characteristic parameters of the samples:

Cov(X,Y,Z)=Cov(X,X)Cov(X,Y)Cov(X,Z)[Cov(Y,X)Cov(Y,Y)Cov(Y,Z)]Cov(Z,X)Cov(Z,Y)Cov(Z,Z)

wherein Z represents the data of the characteristic parameters of the sample, and

the covariance matrix comprising the data of the characteristic parameters of p types of samples is:

Cov(K)=[Cov(K1,1)…Cov⁢(K1,p)⋮⋱⋮Cov⁢(Kp,1)…Cov⁢(Kp,p)]

wherein K represents the multiple data of the characteristic parameters of the samples, and Cov(K1,p) represents the covariance of first data of the characteristic parameters of the sample and pthdata of the characteristic parameters of the sample;

Step B500, based on the covariance matrix comprising p types of data of the characteristic parameters of the samples, acquiring the characteristic values,, . . . ,and p characteristic vectors,, . . .of the covariance matrix; and

based on the characteristic values and the characteristic vectors, calculating the PCA data comprising p types of characteristic values:

mj=∑j=1pêJ⁢nj,

wherein mjrepresents a jthvector of the PCA data; njrepresents a jthvector of the data of the characteristic parameters of the sample;represents a jthcharacteristic vector of the data of the characteristic parameters of the sample; and j meets the condition: j∈1, 2, . . . p and represents the counting parameter; and

Step B600, based on the PCA data comprising p types of characteristic values, sorting the corresponding characteristic values from large to small, so as to obtain the characteristic values,, . . . ,sorted from large to small, and calculating the accumulative contribution of the covariance:

Ck=∑j=1kλjΣjp⁢λj

wherein λjrepresents the jthcharacteristic value sorted from large to small, Ckrepresents the accumulative contribution of the front k covariance, and k meets the condition: k∈1, 2, . . . , p; and

when k is constantly increased from 1, selecting the vector of the PCA data, which corresponds to the corresponding front k covariance when Ckis more than the third threshold for the first time, i.e. the PCA data obtained by the PCA method.

In some preferred implementation manners, Step S400specifically comprises:

Step S410, based on the data of the seismic wave reflected signal, expressing a frequency domain seismic record convolution model as:
s(ω)=σ(ω)ε(ω)

wherein s(ω) represents the seismic record after Fourier transformation; σ(ω) represents the wavelet after Fourier transformation; ε(ω) represents the frequency spectrum of the reflection coefficient after Fourier transformation; and co represents the angular frequency;

Step S420, transforming logarithms on the two sides of the equation of the frequency domain seismic record convolution model into a linear system to obtain a linear seismic record convolution model:
lns(ω)=lnσ(co)+lnε(ω)

Step S430, carrying out inverse Fourier transformation on the linear seismic record convolution model to obtain a complex cepstrum sequence:
{tilde over (s)}(t)={tilde over (σ)}(t)+{tilde over (ε)}(t)

wherein {tilde over (s)}(t) represents a complex cepstrum sequence of a seismic waveform record; {tilde over (σ)}(t) represents a complex cepstrum sequence of the seismic wavelet; {tilde over (ε)}(t) represents a complex cepstrum sequence of the stratum reflection coefficient; and t represents the time of the seismic waveform record;

Step S440, based on the complex cepstrum sequence, separating the wavelet and the reflection coefficient by a low-pass filter and extracting a wavelet amplitude spectrum;

Step S450, acquiring a simulated seismic wavelet amplitude spectrum by a least square method:
σ(f)=A(f)fαeG(f)

wherein the fitting parameter α≥0 of the least square method is a constant; σ(f) represents the wavelet amplitude spectrum; H(f) and A(f) represent the polynomials off; and f represents the frequency of the seismic wavelet;

Step S460, based on the simulated seismic wavelet amplitude spectrum, obtaining the maximum phase component and the minimum phase component of the wavelet,

wherein it is assumed that the maximum phase component of the wavelet σ(t) is u(t), and the minimum phase component of the wavelet σ(t) is v(t), the wavelet σ(t) is:
σ(t)=u(t)·v(t);

and the complex cepstrum of the amplitude spectrum is expressed as follows:
2(t)=ũ(t)+{tilde over (v)}(t)+ũ(−t)+{tilde over (v)}(−t)

wherein the complex cepstrum(t) of the amplitude spectrum is displayed symmetrically on the positive axis and the negative axis; {tilde over (v)}(−t) represents the complex cepstrum of the minimum phase function corresponding to the maximum phase component v(t) of the seismic wavelet; and ũ(−t) represents the complex cepstrum of the maximum phase function corresponding to the minimum phase component u(t) of the seismic wavelet;

Step S470, based on the complex cepstrum in the amplitude spectrum, determining a group of mixed-phase wavelet sets with the same amplitude spectrum; constantly adjusting the parameters of the Yu wavelet; keeping the low frequency, broadening the high frequency and properly increasing the dominant frequency to construct the expected output wavelet form; and looking for the best balance point for improving the resolution ratio and the fidelity with reference to a signal-to-noise ratio spectrum under the control of the logging curves, so as to obtain data of the waveform after shaping;

Step S480, based on the data of the waveform after shaping, constructing a tensor diffusion model:

{∂∪∂τ=div(D·∇U)U|τ=0=U0

wherein T represents the diffusion time; div represents the divergence operator; D represents the diffusion tensor; U represents the diffusion filtering result; U|τ=0represents the diffusion filtering result when τ=0; U0represents the data of the waveform after shaping when τ=0; and as the initial condition of the tensor diffusion model, ∇U represents the gradient of the diffusion filtering result;

based on the tensor diffusion model, constructing the gradient structure tensor:
Jβ(∇U)=Gρ*J0=Gρ*(∇U(∇U)T)

wherein U represents the diffusion filtering result, and J0represents the tensor product of the gradient vector;

the Gaussian function with the scale of ρ is expressed by Gρ:
Gρ(r)=(2πρ2)−3/2exp(−|r|2/(2ρ2))

wherein r represents the calculation radius;

the characteristic vector of the structure tensor is:

Jρ(∇U)=(v1⁢v2⁢v3)⁢(ξ1000ξ2000ξ3)⁢(v1⁢v2⁢v3)T

wherein v1, v2and v3represent three characteristic vectors of the gradient structure tensor and can be regarded as a local orthogonal coordinate system; v1points to the gradient direction of the seismic signal; a plane formed by v2and v3is in parallel to the local structure characteristics of the seismic signal; and ξ1, ξ2and ξ3represent three characteristic values corresponding to v1, v2and v3respectively;

Step S490, based on the characteristic vector of the structure tensor, respectively calculating the confidence degree of a linear structure, the confidence degree of a face structure and the diffusion tensor, wherein

the confidence degree Clineof the linear structure is:

Cl⁢i⁢n⁢e=λ2-λ3λ2+λ3;

the confidence degree Cplaneof the face structure is:

Cp⁢l⁢a⁢n⁢e=λ1-λ2λ1+λ2;

the diffusion tensor D is:

D=(v1⁢v2⁢v3)⁢(μ1000μ2000μ3)⁢(v1⁢v2⁢v3)T=(D1⁢1D1⁢2D1⁢3D2⁢1D2⁢2D2⁢3D3⁢1D3⁢2D3⁢3)

wherein μ1, μ2and μ3represent three non-negative characteristic values of the diffusion tensor and respectively represent the intensity of filtering of the diffusion filter in the three characteristic directions v1, v2and v3; and

Step S4100, repeating the steps S480-S490until the preset number of iterations is achieved, so as to obtain the diffusion filtering result, i.e. the high-precision 3D seismic amplitude data body.

In some preferred implementation manners, a calculation method of the simulated seismic wavelet amplitude spectrum comprises:

positioning the maximum value of an amplitude spectrum in the data of the seismic wave reflected signal and the frequency corresponding to the maximum value;

fitting the maximum value of the amplitude spectrum of the seismic signal and the simulated seismic wavelet amplitude spectrum by the least square method to obtain the coefficients of the polynomials of the parameters a and H (f), so as to obtain the fit amplitude value of the corresponding frequency of the maximum value; and

dividing the fit amplitude value of the corresponding frequency by the maximum value of the amplitude spectrum of the seismic signal, so as to fit the coefficient of the polynomial A(f) by the quotient.

In some preferred implementation manners, the time-depth transforming relationship between the data of the logging curve and the seismic record is:

Ts=TH0+2⁢∑i=1nTi·ΔH

wherein TH0represents the time of the seismic signal corresponding to the initial depth of acoustic logging; Tirepresents the time difference of the acoustic wave; i represents the sequence number of time operation between every two sampling points; ΔH represents the sampling interval of the data of the logging curve; and Tsrepresents the two-way travel time of the seismic wave.

In some preferred implementation manners, Step S800specifically comprises Step S810-Step S890:

Step S810, optionally selecting a sample well as a reference target well and setting the parameter of the initial sample number as 1;

Step S820, according to the principle of similarity of the waveforms, selecting the first group of PCA data of the sample wells and a first group of PCA data of the reference target well for correlation analysis, so as to obtain the correlation value of the parameter of the sample number and the first group of PCA data of the reference target well, wherein the number of pieces of first group of PCA data of the sample wells is equal to the parameter of the sample number;

Step S830, increasing the parameter of the sample number1by1; repeating the method of Step S820to obtain the correlation value of the parameter of the sample number and the first group of PCA data of the reference target well, which corresponds to each parameter of the sample number; and connecting all the correlation values of the parameter of the sample number and the first group of PCA data of the reference target well, so as to obtain a correlation curve with the correlation of the first group of PCA data of the reference well changing with the parameter of the sample number;

Step S840, selecting another sample well as a reference target well; repeating the methods of Step S810-Step S830to obtain multiple correlation curves with the correlation of the first group of PCA data of the reference well changing with the parameter of the sample number; fitting all the correlation curves with the correlation of the waveforms of the first group of PCA data of the reference well changing with the parameter of the sample number into an overall correlation curve; selecting an inflection point capable of enabling the correlation in the overall correlation curve to be increased with the increase of the parameter of the sample number and to be finally kept at a stable position; and determining the parameter of the optimal sample number;

Step S850, based on the high-precision 3D seismic amplitude data body and the isochronous framework model, calculating the correlations of the waveforms of point positions to be detected and the positions of the sample wells; sorting the correlations of the waveforms from large to small; selecting the first group of PCA data of the sample wells having the highest correlation of the seismic waveforms and the number of pieces of first group of PCA data being equal to the parameter of the optimal sample number; and based on the sample wells corresponding to the characteristic data of the seismic waveforms of the sample wells having the highest correlation, constructing the initial model by an inter-well characteristic parameter interpolation manner;

Step S860, based on the initial model, selecting the first group of PCA data of the sample wells as prior information, wherein the first group of PCA data has the highest correlation of the seismic waveforms, and the number of pieces of first group of PCA data is equal to the parameter of the optimal sample number;

Step S870, carrying out matched filtering on the initial model and the prior information, so as to obtain the maximum likelihood function;

Step S880, based on the maximum likelihood function and the prior information, solving the statistical distribution density of the posterior probabilities under a Bayesian framework, and sampling the statistical distribution density of the posterior probabilities to obtain the objective function; and

Step S890, sampling the distribution of the posterior probabilities by an MCMC (Markov Chain Monte Carlo) method and the criterion for Metropolis-Hastings sampling and by taking the objective function as the input of the initial model; constantly optimizing the parameters of the initial model; selecting the solution when the maximum value is selected for the objective function as random implementation; selecting the average value of multiple random implementations as the expected output value; and taking the expected output value as the high-precision characteristic value simulation result data bodies, wherein the parameters in the high-precision characteristic value simulation result data bodies correspond to the first group of PCA data one to one.

In some preferred implementation manners, Step S880specifically comprises:

Step S881, by utilizing the law that the white noise meets Gaussian distribution, expressing the parameter of the high-precision characteristic value simulation result data bodies as:
Ysim=Xsim+Nnoise

wherein Ysimrepresents the parameter of the high-precision characteristic value simulation result data bodies, Xsimrepresents the to-be-solved actual characteristic parameter value of the subsurface stratum, and Nnoiserepresents the random noise;

Step S882, as ∥Xsim−Xp∥2also meets Gaussian distribution, determining the initial objective function as:

J1=xs⁢i⁢m-xp22⁢σ2

wherein Jirepresents the function related to posterior information; Xprepresents the expected value of the characteristic parameter, which is obtained by calculation after the steps of selecting the sample wells based on the optimal sample number, carrying out matched filtering on characteristic curves of the sample wells and then solving the statistical distribution density of the posterior probabilities; and σ represents the covariance of the white noise; and

Step S883, based on the initial objective function, introducing the prior information into the objective function by maximum posterior estimation, so as to obtain the following stable objective function:
J(Zsim)=J1(Zsim)+λ2J2(Zsim)

wherein Zsimrepresents the characteristic parameter to be simulated; J2represents the function related to the prior information of the geological data and the logging data; and λ represents the smoothing parameter used for coordinating the mutual effect between J1and J2.

In some preferred implementation manners, Step S890specifically comprises:

Step S891, setting M as the target space, setting n as the total sample number and setting m as the sample number when the Markov chain is stable;

Step S892, presetting a Markov chain to enable the Markov chain to be converged to stable distribution;

Step S893, starting from a certain point P0in M, carrying out sampling simulation by the Markov chain and generating a point sequence: P1, . . . , Pn

Step S894, expressing the expected estimation of the function J2(Z) as:

E[J2(P)]=1n-m⁢∑L=m+1nf⁡(P(1))

wherein n represents the total sample number, m represents the sample number when the Markov chain is stable, and L represents the accumulated parameter;

Step S895, selecting the transition function q(P; P(h-1)) and the initial value P(0), wherein if the parameter value when the hthiteration is started is P(h-1), the hthiteration process comprises:

extracting the alternative value P′ from q(P; P(h-1)) and calculating the probability of acceptance α(P(h-1), P′) of the alternative value P′:

α⁡(P(h-1),P′)=min⁢{π⁡(P(h-1);P′)π⁡(P(h-1)⁢q⁡(P′;P(h-1)))};

Step S896, enabling P(h)=P with α(P(h-1), P′) and enabling P(h)=P(h-1)with the probability 1−α(P(h-1), P′); and

Step S897, constantly disturbing the parameters of the initial model, repeating the methods of Step S892-Step S896to achieve the preset number of iterations hs, and obtaining posterior samples P(1), P(2), P(3). . . P(hs)so as to calculate the matrix at each order of posterior distribution to obtain the expected output value, wherein the expected output value is taken as the first group of high-precision characteristic value simulation result data bodies.

In some preferred implementation manners, Step S300comprises:

Step S310, drawing data of an original logging curve based on the original logging data;

Step S320, removing outliers to obtain data of the logging curve with the outliers removed based on the data of the original logging curve; and

Step S330, superposing all data of a histogram of the single logging curve at the positions of the sample wells in a work area based on the data of the logging curve with the outliers removed, and integrating the thresholds to obtain the data of the standardized logging curve.

Another aspect of the present invention provides a system for evaluating the filling characteristics of the deep paleokarst reservoir through well-to-seismic integration. The system comprises: an original geophysical logging data acquisition module, a seismic data acquisition module, an original geophysical logging data preprocessing module, a seismic data preprocessing module, a well-to-seismic calibration and characteristic parameter selection module, a module for acquiring the first group of PCA data, an isochronous framework model constructing module, an inter-well reservoir parameter simulation module, a module for characterizing the boundary of the inter-well karst cave system, a module for characterizing the lithological and physical boundaries of internal filling of the cave and a module for describing the structure and the filling of the paleokarst cave;

the original geophysical logging data acquisition module is configured to acquire the original logging data of each sample well by the logging equipment, comprising: measuring the self-potential SP of each sample well by the measuring electrode, measuring the natural gamma ray GR of each sample well by the natural gamma ray downhole device and the natural gamma ray ground instrument, and acquiring the caliper CAL of each sample well by the caliper arm; and acquiring the data of the resistivity curve: the deep lateral logging RLLD, the shallow lateral logging RLLS and the micro-lateral logging RLLM, and the data of the physical property characterization curve: the compensated neutrons CNL, the compensated acoustic curve AC and the density curve DEN by the traditional logging equipment; and

acquiring the definite lithological information and physical information of the individual depth segment based on the imaging logging information, the drilling information, the mud logging information and the rock-core information, so as to determine the depth data of the marker bed at the target horizon;

the seismic data acquisition module is configured to acquire the data of the original seismic wave reflected signal by the seismic wave excitation device and the receiving device and acquire isochronous 3D distribution of the marker bed at the target horizon according to the waveform of the data of the original seismic wave reflected signal;

the original geophysical logging data preprocessing module is configured to draw the data of the logging curves based on the original logging data of the sample wells and carry out abnormal value processing and standardization processing, so as to obtain the data of the standardized logging curves;

the seismic data preprocessing module is configured to obtain the high-precision 3D seismic amplitude data body by mixed-phase wavelet estimation and maximum posteriori deconvolution and enhancing diffusion filtering based on the data of the seismic wave reflected signal;

the well-to-seismic calibration and characteristic parameter selection module is configured to carry out well-to-seismic calibration and characteristic parameter selection: acquiring the wave impedance curve of each sample well based on the compensated acoustic curve and the density curve DEN in the data of the standardized logging curve, so as to calculate the reflection coefficient curve; acquiring the preferred frequency of the ricker wavelet, so that the preferred frequency is consistent with the dominant frequency of the high-precision 3D seismic amplitude data body; carrying out convolution operation on the ricker wavelet and the reflection coefficient curve, so as to obtain the synthetic seismic record; comparing the depth data of the marker bed at the target horizon and the isochronous 3D distribution of the marker bed at the target horizon for well-to-seismic calibration; calculating the correlation between the synthetic seismic record and the waveform of the near-well seismic channel; and when the correlation is more than or equal to the preset first threshold, judging that the result of well-to-seismic calibration is qualified, so as to obtain the time-depth transforming relationship between the data of the logging curve and the seismic record and the characteristic parameters sensitive to the reservoir,

wherein the method for obtaining the characteristic parameters sensitive to the reservoir comprises:

drawing the columnar statistical graph of the responses of logging parameters, which are generated to the different geological bodies near the wells; and when the data of the certain standardized logging curve can be used for distinguishing the data points that are more than the second threshold preset by the different logging interpretation conclusions, selecting the data of the standardized logging curve as the characteristic parameters sensitive to the reservoir, wherein the characteristic parameters sensitive to the reservoir at least comprise the wave impedance IMP and can also comprise one or more of the caliper CAL, the natural gamma ray GR, the self-potential SP, the data of the resistivity curve: the deep lateral logging RLLD, the shallow lateral logging RLLS and the micro-lateral logging RLLM, and the data of the physical property characterization curve: the compensated neutrons CNL, the compensated acoustic curve AC and the density curve DEN;

the module for acquiring the first group of PCA data is configured to optimize the data of the standardized logging curves from the characteristic parameters sensitive to the reservoir, wherein the number of pieces of the data is preset; and carry out dimensionality reduction on the characteristic parameters sensitive to the reservoir by the PCA method, so as to obtain the first group of PCA data;

the isochronous framework model constructing module is configured to construct the isochronous framework model based on the law of sedimentary strata reflected by the high-precision 3D seismic amplitude data body and the time-depth transforming relationship;

the inter-well reservoir parameter simulation module is configured to determine the parameter of the optimal sample number, which can reflect the overall geological condition, based on the first group of PCA data of each sample well; select the first group of PCA data of the sample wells to construct the initial model, wherein the first group of PCA data has the highest correlation of the seismic waveforms, and the number of pieces of first group of PCA data is equal to the parameter of the optimal sample number; constantly correct the parameters of the initial model; and output the first group of high-precision characteristic value simulation result data bodies, wherein the first group of high-precision characteristic value simulation result data bodies is the data bodies corresponding to the first group of PCA data one to one;

the module for characterizing the boundary of the inter-well karst cave system is configured to carry out cross analysis based on the first group of PCA data of the sample wells, so as to obtain the distribution diagram of the first group of discrete data points of various geological bodies;

based on the lithological information and the physical information, acquire the interpretation conclusion of the distribution diagram of the first group of discrete data points of various geological bodies, divide the karst cave and the surrounding rock to construct the cross plot of the first group of PCA data, and acquire the threshold of the first group of PCA data, which is required for division of the karst cave; and

based on the threshold of the first group of PCA data, carry out cross analysis on the first group of high-precision characteristic value simulation result data bodies, so as to obtain the characteristics of the 3D spatial morphologies of the karst cave;

the module for characterizing the lithological and physical boundaries of internal filling of the cave is configured to characterize the lithological and physical boundaries of internal filling of the cave: analyzing the characteristic parameters sensitive to the reservoir, which correspond to the data points of the karst cave, by the PCA method, so as to obtain the second group of PCA data sensitive to the response of the fillings;

based on the second group of PCA data, generating the second group of high-precision characteristic value simulation result data bodies by the functions described in the inter-well reservoir parameter simulation module;

based on the second group of PCA data sensitive to the response of the fillings, carrying out cross analysis to obtain the distribution diagram of the second group of discrete data points of various geological bodies;

based on the lithological information and the physical information, acquiring the interpretation conclusion of the second group of discrete data points of various geological bodies, and distinguishing the types of the fillings to construct the cross plot of the second group of PCA data, so as to obtain the threshold of the second group of PCA data, which is required for division of the types of the internal fillings of the karst cave; and

based on the threshold of the second group of PCA data, carrying out cross analysis on the second group of high-precision characteristic value simulation result data bodies, so as to obtain the characteristics of the 3D spatial morphologies of different types of internal fillings of the karst cave; and

the module for describing the structure and the filling of the paleokarst cave is configured to engrave the spatial distribution of the paleokarst cave and the development characteristics of the internal filling adopting different types of internal fillings by adopting the lithology shielding technology and the 3D engraving technology based on the characteristics of the 3D spatial morphologies of the karst cave and the characteristics of the 3D spatial morphologies of different types of internal fillings of the karst cave.

The third aspect of the present invention provides electronic equipment, which comprises:

at least one processor; and

a memory that is in communication connection with the at least one processor, wherein

the memory is configured to store instructions capable of being executed by the processor(s); and the instructions are used for being executed by the processor(s) for realizing the method for evaluating the filling characteristics of the deep paleokarst reservoir through well-to-seismic integration.

The fourth aspect of the present invention provides a computer readable storage medium; the computer readable storage medium is configured to store computer instructions; and the computer instructions are used for being executed by a computer for realizing the method for evaluating the filling characteristics of the deep paleokarst reservoir through well-to-seismic integration.

The method and the system of the present invention have the beneficial effects that:

in the present invention, under the conditions of fully utilizing the logging data and the seismic data, adopting a machine learning algorithm and taking a wave indication simulation method as a core, sensitive parameters are established by using a PCA algorithm, so that the structure and the filling of the reservoir are precisely identified by the logging data; and then, 3D seismic simulation data bodies are established by using the sensitive parameters and the 3D seismic data, so that the logging interpretation result with high resolution is extrapolated to the inter-well 3D large-scale area, and the filling characteristics of the deep paleokarst reservoir are evaluated finally.

In the present invention, the 3D parameters are simulated through well-to-seismic integration based on an initial isochronous geological framework model by utilizing the correlation between high-precision 3D seismic waveform data; coordinate systems are transformed for multiple logging curves after being preprocessed by utilizing the PCA method twice; and dimensionality reduction analysis is carried out on the lithological information and the physical information that are reflected by the data of multiple conventional logging curves, so that a small amount of PCA data has the capability of characterizing the reservoir by multiple logging curves.

In the present invention, the karst cave and the data points of internal filling can be well distinguished by a 2D (Two-Dimensional) cross plot constructed based on the PCA data, so as to obtain a plate of a cross plot of identification for the karst reservoir; and threshold information for division of the reservoir, which is reflected by the plate, is applied to simulated data of the 3D parameters, so as to obtain the characterization conclusion of longitudinal and lateral structures of the karst cave and the internal filling characteristics thereof under the constraint of PCA cross analysis. Therefore, the characterization precision is improved, the reservoir with fast lateral change can be predicted, and the reliable theoretical basis is provided for evaluating the reservoir property and the oil-gas migration of a carbonate paleokarst cave on a large scale.

According to the method and the system of the present invention, the mapping relationship between geophysical detection methods is established by adopting the Markov chain sampling criterion and a Monte Carlo estimation method under the guidance of the Bayesian theory, so as to achieve the effect of identifying the development characteristics of a carbonate paleokarst cave type reservoir in a large-scale complex basin.

