阅读说明：本技术 井震联合评价深层古岩溶储层充填特征的方法与系统 (Method and system for evaluating filling characteristics of deep paleo-karst reservoir stratum by well-seismic combination ) 是由 田飞 张江云 底青云 郑文浩 王中兴 杨永友 张文秀 于 2021-09-13 设计创作，主要内容包括：本发明属于数据识别、记录载体的处理领域,具体涉及了一种井震联合评价深层古岩溶储层充填特征的方法与系统,旨在解决现有的石油勘探技术无法预测横向变化快的储层、无法识别大范围内复杂盆地碳酸盐岩洞穴型储层的发育特征的问题。本发明包括：获取标准化测井曲线数据；通过混合相位子波反褶积和扩散滤波,获得高精度的三维地震振幅数据体；通过PCA方法分析对所述对储层敏感的特征参数进行降维,获得第一组PCA数据；计算缝洞型储层结构特征；通过解释结论门槛值进行交会分析获得古岩溶洞穴空间展布和不同填充物类型的洞穴内部充填的发育特征。本发明达到识别大范围内复杂盆地碳酸盐岩岩溶洞穴型储层的发育特征的效果提高了刻画的精度。(The invention belongs to the field of data identification and record carrier processing, and particularly relates to a method and a system for evaluating filling characteristics of deep ancient karst reservoirs by well-seismic combination, aiming at solving the problems that the existing petroleum exploration technology cannot predict reservoirs with rapid transverse change and cannot identify development characteristics of complex basin carbonate rock cave type reservoirs in a large range. The invention comprises the following steps: acquiring standardized logging curve data; obtaining a high-precision three-dimensional seismic amplitude data volume through deconvolution and diffusion filtering of mixed phase wavelets; performing dimensionality reduction on the characteristic parameters sensitive to the reservoir through analysis of a PCA method to obtain a first group of PCA data; calculating structural characteristics of the fracture-cavity reservoir; and performing intersection analysis by explaining the conclusion threshold value to obtain the space distribution of the ancient karst cave and the development characteristics of the filling inside the cave with different filler types. The method achieves the effect of identifying the development characteristics of the complex basin carbonate karst cave type reservoir in a large range, and improves the precision of the carving.)

1. A method for evaluating deep paleo-karst reservoir filling characteristics through well-seismic combination is characterized by comprising the following steps:

step S100, obtaining original geophysical logging data: obtaining raw logging data for each sample well by a logging device, comprising: measuring the natural potential SP of each sample well through a measuring electrode, measuring the natural gamma GR of each sample well through a natural gamma underground device and a natural gamma ground instrument, and obtaining the well diameter CAL of each sample well through a well diameter arm; obtaining resistivity curve data by conventional logging equipment: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and physical property characterization curve data: a compensated neutron CNL, a compensated acoustic curve AC and a density curve DEN; (ii) a

Obtaining determined lithology information and physical property information of individual depth sections based on imaging logging information, drilling information, logging information and core information, and further determining depth data of a target horizon marker layer;

step S200, acquiring seismic data, acquiring original seismic wave reflection signal data through a seismic wave excitation device and a receiving device, and acquiring isochronous three-dimensional spread of a target layer position mark layer according to the waveform of the original seismic wave reflection signal data;

step S300, preprocessing original geophysical logging data: drawing logging curve data based on all original logging data of the sample well, and performing abnormal value processing and standardization processing to obtain standardized logging curve data;

step S400, seismic data preprocessing: based on the seismic wave reflection signal data, obtaining a high-precision three-dimensional seismic amplitude data volume through mixed phase wavelet deconvolution and diffusion filtering;

step S500, well seismic calibration and characteristic parameter selection: obtaining a wave impedance curve of a sample well based on a compensation acoustic curve and a density curve DEN in the standardized logging curve data, further calculating a reflection coefficient curve, obtaining the preferred frequency of a Rake wavelet to keep the same with the main frequency of a high-precision three-dimensional seismic amplitude data body, performing convolution operation on the Rake wavelet and the reflection coefficient curve to obtain a synthetic seismic record, comparing the depth data of a target layer position mark layer with the isochronous three-dimensional spread of the target layer position mark layer for well seismic calibration, calculating the correlation between the synthetic seismic record and a well side seismic channel waveform, judging that the well seismic calibration result is qualified when the correlation is greater than or equal to a preset first threshold, and obtaining the time-depth conversion relation between the logging curve data and the seismic record and characteristic parameter sensitive to a reservoir layer

The characteristic parameter sensitive to the reservoir is obtained by the method comprising the following steps:

drawing a histogram by responding to logging parameters generated by different geologic bodies beside a well, and selecting the standardized logging curve data as characteristic parameters sensitive to a reservoir when the numerical value of the standardized logging curve data can distinguish data points above a second threshold value preset by different logging interpretation conclusions; the characteristic parameters sensitive to the reservoir at least comprise wave impedance IMP, well diameter CAL, natural gamma GR, natural potential SP, resistivity curve data: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and the data of the physical property characterization curve: one or more of a compensated neutron CNL, a compensated acoustic curve AC, a density curve DEN;

step S600, obtaining a first group of PCA data: optimizing a preset number of pieces of standardized well logging curve data from the characteristic parameters sensitive to the reservoir, and performing dimensionality reduction on the characteristic parameters sensitive to the reservoir through a PCA method to obtain a first set of PCA data sensitive to the reservoir;

step S700, constructing an isochronous trellis model: constructing an isochronous trellis model based on sedimentary stratum rules reflected by the high-precision three-dimensional seismic amplitude data volume and the time-depth conversion relation;

step S800, interwell reservoir parameter simulation: determining an optimal sample number parameter capable of reflecting the overall geological condition based on a first group of PCA data of each sample well, selecting the first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample number parameter to construct an initial model, continuously correcting parameters of the initial model, and outputting a first group of high-precision characteristic value simulation result data volumes which are data volumes corresponding to the first group of PCA data one by one;

step S900, depicting the boundary of the karst cave system between wells: performing intersection analysis based on a first group of PCA data of the sample well to obtain a first group of geologic body discrete data point distribution maps;

obtaining an interpretation conclusion of a first group of geologic body discrete data point distribution diagrams based on the lithological information and the physical property information, dividing a karst cave and surrounding rocks, further constructing a first group of PCA data intersection layouts, and obtaining a first group of PCA data threshold values required by the division of the karst cave;

performing intersection analysis on the first group of high-precision characteristic value simulation result data volumes based on the first group of PCA data threshold values to obtain three-dimensional space morphological characteristics of the karst cave;

step S1000, depicting a lithologic property boundary filled in the cave: analyzing characteristic parameters sensitive to a reservoir, which correspond to the karst cave data points, by a PCA method based on the karst cave data points of the first group of PCA data intersection layout to obtain a second group of PCA data sensitive to the response of the filling;

generating a second group of high-precision characteristic value simulation result data volume by the method of S800 based on the second group of PCA data;

performing intersection analysis based on the second group of PCA data sensitive to the filler response to obtain a second group of geologic body discrete data point distribution diagrams;

obtaining an interpretation conclusion of a second group of geologic body discrete data point distribution diagrams based on the lithological information and the physical property information, distinguishing filler types, further constructing a second group of PCA data intersection layout, and obtaining a second group of PCA data threshold values required by the type division of the filler in the ancient karst cave;

performing intersection analysis on the second group of high-precision characteristic value simulation result data volumes based on the second group of PCA data threshold values to obtain three-dimensional space morphological characteristics of different filler types in the karst cave;

step S1100, archaeological cavern structure and filling description: and engraving the development characteristics of the ancient karst cave space distribution and the internal filling of different filler types by adopting a lithologic shielding technology and a three-dimensional engraving technology based on the three-dimensional space morphological characteristics of the karst cave and the three-dimensional space morphological characteristics of different filler types in the karst cave.

2. The method for jointly evaluating the filling characteristics of deep paleo-karst reservoirs through well-seismic according to claim 1, wherein the measuring of the natural potential SP of each sample well through the measuring electrode and the measuring of the natural gamma GR of each sample well through the natural gamma downhole device and the natural gamma surface instrument specifically comprise:

measuring the natural potential SP of each sample well through the measuring electrode:

arranging a measuring electrode N on the ground, and arranging a measuring motor M underground through a cable;

lifting the measuring electrode M along the well axis to measure the change of the natural potential along with the well depth;

the calculation method of the natural potential value comprises the following steps:

wherein the content of the first and second substances,the total natural potential is the total natural potential,in order to be a diffusion potential coefficient,in order to obtain a diffusion adsorption potential coefficient,the resistivity value of the slurry filtrate is shown as,is the formation water resistivity value;

the natural gamma GR of each sample well is measured through the natural gamma underground device and the natural gamma ground instrument:

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

acquiring natural gamma rays through a detector, converting the natural gamma rays into electric pulse signals, and amplifying the electric pulse signals through an amplifier;

the ground instrument converts the electric pulse count formed every minute into a potential difference for recording.

3. The method for jointly evaluating the filling characteristics of deep paleo-karst reservoirs through well-seismic analysis according to claim 1, wherein the PCA analysis specifically comprises:

step B100, calculating a sample mean value:

wherein the content of the first and second substances,represents the total number of samples of the curve, i represents the sample count parameter,represents the characteristic parameter data of the ith sample,represents the sample mean;

step B200, carrying out normalization processing on each parameter curve to obtain normalized sample characteristic parameter data：

Wherein x represents sample feature parameter sample point data;

step B300, based on the normalized sample characteristic parameter dataCalculating the covariance of the characteristic parameters of the sample：

Wherein the content of the first and second substances,represents the c-th X sample characteristic parameter data,representing the c sample characteristic parameter data;

step B400, constructing a sample characteristic parameter covariance matrix based on the sample characteristic parameter covariance:

z represents sample characteristic parameter data;

Includedpthe covariance matrix of the seed sample characteristic parameter data is:

wherein the content of the first and second substances,representing a plurality of sample characteristic parameter data,representing the characteristic parameter data of the 1 st sample and thepCovariance of the sample feature parameter data;

step B500, based on the inclusionpThe covariance matrix of the sample characteristic parameter data is planted, and the eigenvalue of the covariance matrix is obtainedAndpa feature vector，，；

Eigenvalue and eigenvector computations based on the covariance matrix includepPCA data for species eigenvalues:

wherein the content of the first and second substances,represents the jth PCA data vector and,representing the jth sample characteristic parameter data vector,a feature vector representing the jth sample feature parameter data;j 1,2, pis a counting parameter;

step B600 based onpAnd (4) sequencing the corresponding eigenvalues from large to small according to the PCA data of the variety eigenvalues to obtain the eigenvalues sequenced from large to small，…，Calculating the variance cumulative contribution:

wherein the content of the first and second substances,representing the jth eigenvalue sorted from large to small,to accumulate the first k covariance accumulated contributions,；

when k is increased from 1, the k is selectedAnd when the first time is greater than a preset third threshold value, corresponding PCA data vectors corresponding to the first k covariances are PCA data obtained by PCA analysis in the current round.

4. The method for jointly evaluating the filling characteristics of the deep paleo-karst reservoir by well-seismic according to claim 1, wherein the step S400 specifically includes:

step S410, based on the seismic wave reflection signal data, expressing a frequency domain seismic record convolution model as:

wherein the content of the first and second substances,representing the fourier transformed seismic records,representing the wavelet after the fourier transform,a frequency spectrum representing the fourier transformed reflection coefficient,represents angular frequency;

step S420, logarithm taking two sides of the equation of the frequency domain seismic record convolution model is converted into a linear system, and a linear seismic record convolution model is obtained:

step S430, performing inverse Fourier transform on the linear seismic record convolution model to obtain a cepstrum sequence:

wherein the content of the first and second substances,representing a repeating spectral sequence of seismic waveform recordings,a cepstrum sequence representing seismic wavelets,a repeating spectral sequence representing the reflection coefficients of the formation,representing a seismic waveform recording time;

step S440, based on the cepstrum sequence, performing wavelet and reflection coefficient separation through a low-pass filter, and extracting a wavelet amplitude spectrum;

step S450, obtaining the amplitude spectrum of the simulated seismic wavelet by a least square method:

wherein, least square method is used for fitting parametersIs a constant number of times, and is,a representation of the amplitude spectrum of the wavelet is obtained,anda polynomial expression representing f, which represents the frequency of the seismic wave;

step S460, obtaining a wavelet maximum phase component and a wavelet minimum phase component based on the simulated seismic wavelet amplitude spectrum;

wavelet setting deviceHas a maximum phase component ofThe minimum phase component isWavelet of fundamental waveComprises the following steps:

the magnitude spectrum is represented in the cepstrum as:

wherein the repetition spectrum of the amplitude spectrumSymmetrically displayed on the positive and negative axes of the match score,for maximum phase component of seismic waveletThe cepstrum of the corresponding minimum phase function,for minimum phase component of seismic waveletsThe corresponding cepstrum of the maximum phase function;

step S470, determining a group of mixed phase wavelet sets with the same amplitude spectrum based on the cepstrum in the amplitude spectrum, continuously adjusting the parameters of Shu' S wavelets, maintaining low frequency, expanding high frequency and properly improving dominant frequency to construct an expected output wavelet form, searching for an optimal balance point between resolution and fidelity by taking a signal-to-noise ratio spectrum as a reference under the control of a well curve, and obtaining waveform data after shaping;

step S480, constructing a tensor diffusion model based on the shaped waveform data:

wherein the content of the first and second substances,it is shown that the time of diffusion,denotes a divergence operator, D denotes a tensor-type diffusion coefficient of the diffusion filter, U denotes a diffusion filtering result,to representThe result of diffusion filtering when =0,to representThe waveform data after the shaping at that time is used as an initial condition of the tensor diffusion model,a gradient representing a result of the diffusion filtering;

constructing a gradient structure tensor based on the tensor diffusion model:

where, U represents the result of the diffusion filtering,representing a gradientA vector tensor product;

the expression scale isGaussian function of (d):

wherein r represents the calculated radius;

the eigenvectors of the structure tensor are:

wherein the content of the first and second substances,、andthe 3 eigenvectors, expressed as gradient structure tensors, can be considered as local orthogonal coordinate systems,pointing in the direction of the gradient of the seismic signal,andthe planes of composition are parallel to the local structural features of the seismic signal,、andare respectively connected with、Andcorresponding three characteristic values;

step S490, calculating a linear structure confidence metric, a planar structure confidence metric, and a diffusion tensor, respectively, based on the feature vectors of the structure tensor;

the presence structure confidence measureComprises the following steps:

the planar structure confidence measureComprises the following steps:

the diffusion tensor D is:

wherein the content of the first and second substances,、andthree non-negative eigenvalues representing the diffusion tensor, each representing a diffusion filter edge、Andthe filter strengths of the three characteristic directions;

and S4100, repeating the steps of S480-S490 until a preset iteration number is reached, and obtaining a diffusion filtering result, namely the high-precision three-dimensional seismic amplitude data volume.

5. The well-seismic joint evaluation method for deep paleo-karst reservoir filling characteristics according to claim 4, wherein the calculation method of the simulated seismic wavelet amplitude spectrum is as follows:

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

obtaining parameters by fitting the maximum value of the seismic signal amplitude spectrum and the simulated seismic wavelet amplitude spectrum in a least square method modeAndobtaining the corresponding frequency amplitude value of the fitted maximum value by the coefficients of the polynomial;

dividing the maximum value of the seismic signal amplitude spectrum by the fitted amplitude value of the corresponding frequency, and further using a quotient fitting polynomialThe coefficient of (a).

6. The method for evaluating the filling characteristics of the deep paleo-karst reservoir by well-seismic combination according to claim 1, wherein the step S800 specifically includes:

step S810, selecting a sample well as a reference target well, and setting an initial sample number parameter to be 1;

step S820, selecting a first group of PCA data of the sample wells with the quantity being the sample number parameter and a first group of PCA data of the reference target well mark according to the waveform similarity principle to perform correlation analysis to obtain a first group of PCA data correlation values of the sample number parameter-reference target well;

step S830, increasing the sample number parameters 1 by 1, repeating the method of the step S720 to obtain a first group of PCA data correlation values of the sample number parameter-reference target well corresponding to each sample number parameter, connecting the first group of PCA data correlation values of all the sample number parameters-reference target wells, and obtaining a correlation curve of the first group of PCA data correlation of the reference wells along with the change of the sample number parameters;

step S840, selecting another sample well as a reference target well, repeating the steps S810-S830, obtaining a correlation curve of the first group of PCA data correlations of the plurality of reference wells along with the change of the sample number parameters, fitting the correlation curve of the first group of PCA data waveform correlations of all the reference wells along with the change of the sample number parameters into an overall correlation curve, selecting an inflection point at which the correlation in the overall correlation curve rises along with the increase of the sample number parameters and finally keeps stable, and determining the optimal sample number parameters;

step S850, calculating the waveform correlation of the point to be detected and the sample well position based on the high-precision three-dimensional seismic amplitude data volume and the isochronous grid model, and sequencing the waveform correlation from large to small; constructing an initial model based on the sample well corresponding to the seismic waveform characteristic data of the sample well with the highest correlation through an inter-well characteristic parameter interpolation mode;

step S860, based on the initial model, selecting a first group of PCA data of the sample well with the optimal sample number parameter bar and the highest seismic waveform correlation degree as prior information;

step S870, performing matched filtering on the initial model and the prior information to obtain a maximum likelihood function;

step S880, based on the maximum likelihood function and the prior information, obtaining the posterior probability statistical distribution density under a Bayes framework, and sampling the posterior probability statistical distribution density to obtain a target function;

step 890, taking the target function as the input of the initial model, sampling posterior probability distribution by a Markov chain Monte Carlo method MCMC and Metropolis-Hastings sampling criterion, continuously optimizing parameters of the initial model, selecting a solution of the target function when the maximum value is taken as random realization, taking the average value of multiple random realization as expected value output, and taking the expected value output as a high-precision characteristic value simulation result data body; and parameters in the high-precision characteristic value simulation result data volume correspond to characteristic parameters corresponding to the first group of PCA data one to one.

7. The method for jointly evaluating the filling characteristics of the deep paleo-karst reservoir by well-seismic according to claim 6, wherein the step S880 specifically comprises:

step S881, using white noise to satisfy the rule of gaussian distribution, representing the parameters of the high-precision eigenvalue simulation result data volume as:

y represents parameters of a logging curve high-precision characteristic value simulation result data volume, X represents actual characteristic parameter values of an underground stratum to be solved, and N represents random noise;

step S882, becauseAlso satisfying a gaussian distribution, the initial objective function can be determined as:

wherein the content of the first and second substances,a function relating to a posteriori information is represented,showing that the characteristic curve of the sample well is matched and filtered by selecting the sample well based on the optimal sample number, the posterior probability statistical distribution density is obtained, and the expected value of the characteristic parameter is further calculated,a covariance representing white noise;

step S883, based on the initial objective function, introducing prior information into the objective function through maximum posterior estimation, and obtaining a stable objective function as:

wherein the content of the first and second substances,representing the characteristic parameter to be simulated,representing functions related to prior information such as geological and well log data,representation for coordinationAndthe smoothing parameters of the mutual influence between them.

8. The well-seismic joint evaluation method for deep paleo-karst reservoir filling characteristics according to claim 6, wherein the step S890 includes the following specific steps:

step S891, setting M as a target space, n as the total number of samples, and M as the number of samples when the Markov chain tends to be stable;

step S892, presetting a Markov chain to make the Markov chain converge to a stable distribution;

step S893, starting from a certain point in MStarting from, sampling simulation is performed through a Markov chain to generate a point sequence:；

step S894, functionThe expected estimate of (c) is:

wherein n represents the total number of generated samples, m represents the number of samples when the Markov chain reaches a plateau, and k represents an accumulation parameter;

step S895, selecting a transfer functionAnd an initial valueIf the parameter value at the beginning of the ith iteration isThen, the ith iteration process is:

fromExtract an alternative valueCalculating alternative valuesProbability of acceptance of：

Step S896, in order toIs arranged atBy probabilityIs arranged at；

Step S897, continuously disturbing the parameters of the initial model, repeating the method of the steps C792-C796 to reach a preset iteration number n, and obtaining a posterior sampleAnd further calculating each order matrix of the delay distribution to obtain an expected output value, and outputting the expected value as a first group of high-precision characteristic value simulation result data volume.

9. The method for jointly evaluating the filling characteristics of deep paleo-karst reservoirs through well-seismic according to claim 1, wherein the step S300 includes:

step S310, drawing original logging curve data based on the original logging data;

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

and S330, superposing single logging curve histogram data of all sample well positions in the work area based on the logging curve data without outliers, and obtaining standardized logging curve data by integrating threshold values.

10. A system for jointly evaluating filling characteristics of deep paleo-karst reservoirs through well-seismic, the system comprising: 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 seismic calibration and characteristic parameter selection module, a first group of PCA data acquisition modules, an isochronous grid model construction module, an interwell reservoir parameter simulation module, an interwell karst cave system boundary description module, a lithologic property boundary description module for describing filling inside a cave and an ancient karst cave structure and filling description module;

the original geophysical logging data acquisition module is configured to acquire original logging data of each sample well through logging equipment, and comprises: measuring the natural potential SP of each sample well through a measuring electrode, measuring the natural gamma GR of each sample well through a natural gamma underground device and a natural gamma ground instrument, and obtaining the well diameter CAL of each sample well through a well diameter arm; obtaining resistivity curve data by conventional logging equipment: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and physical property characterization curve data: a compensated neutron CNL, a compensated acoustic curve AC and a density curve DEN;

obtaining determined lithology information and physical property information of individual depth sections based on imaging logging information, drilling information, logging information and core information, and further determining depth data of a target horizon marker layer;

the seismic data acquisition module is configured to acquire original seismic wave reflection signal data through a seismic wave excitation device and a receiving device, and acquire isochronous three-dimensional spread of a target layer position mark layer according to the waveform of the original seismic wave reflection signal data;

the original geophysical logging data preprocessing module is configured to draw logging curve data based on original logging data of the sample well, perform abnormal value processing and standardization processing, and obtain standardized logging curve data;

the seismic data preprocessing module is configured to obtain a high-precision three-dimensional seismic amplitude data volume through mixed phase wavelet deconvolution and diffusion filtering based on the seismic wave reflection signal data;

the well-seismic calibration and characteristic parameter selection module is configured to perform well-seismic calibration and characteristic parameter selection: acquiring a wave impedance curve of a sample well based on a compensation acoustic curve and a density curve DEN in the standardized logging curve data, further calculating a reflection coefficient curve, acquiring the preferred frequency of a Rake wavelet to keep the same as the main frequency of a high-precision three-dimensional seismic amplitude data body, performing convolution operation on the Rake wavelet and the reflection coefficient curve to obtain a synthetic seismic record, comparing the depth data of a target layer position mark layer with the isochronous three-dimensional spread of the target layer position mark layer for well seismic calibration, calculating the correlation between the synthetic seismic record and a well-side seismic channel waveform, and judging that a well seismic calibration result is qualified when the correlation is greater than or equal to a preset first threshold value, thereby acquiring the time-depth conversion relation between the logging curve data and the seismic record and characteristic parameters sensitive to a reservoir;

the characteristic parameter sensitive to the reservoir is obtained by the method comprising the following steps:

drawing a histogram by responding to logging parameters generated by different geologic bodies beside a well, and selecting the standardized logging curve data as characteristic parameters sensitive to a reservoir when the numerical value of the standardized logging curve data can distinguish data points above a second threshold value preset by different logging interpretation conclusions; the characteristic parameters sensitive to the reservoir at least comprise wave impedance IMP, well diameter CAL, natural gamma GR, natural potential SP, resistivity curve data: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and the data of the physical property characterization curve: one or more of a compensated neutron CNL, a compensated acoustic curve AC, a density curve DEN;

the first group of PCA data acquisition modules are configured to optimize a preset number of pieces of standardized well logging curve data from the characteristic parameters sensitive to the reservoir, and perform dimensionality reduction on the characteristic parameters sensitive to the reservoir through PCA analysis to obtain a first group of PCA data;

the isochronous trellis model building module is configured to build an isochronous trellis model based on sedimentary stratum rules reflected by the high-precision three-dimensional seismic amplitude data volume and the time-depth conversion relation;

the interwell reservoir parameter simulation module is configured to determine an optimal sample number parameter capable of reflecting the overall geological condition based on a first group of PCA data of each sample well, select a first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample number parameter to construct an initial model, continuously correct parameters of the initial model, and output 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 are data bodies corresponding to the first group of PCA data one by one;

the boundary delineation module of the inter-well karst cave system is configured to perform intersection analysis based on a first group of PCA data of the sample well to obtain a first group of geologic body-like discrete data point distribution maps;

obtaining an interpretation conclusion of a first group of geologic body discrete data point distribution diagrams based on the lithological information and the physical property information, dividing a karst cave and surrounding rocks, further constructing a first group of PCA data intersection layouts, and obtaining a first group of PCA data threshold values required by the division of the karst cave;

performing intersection analysis on the first group of high-precision characteristic value simulation result data volumes based on the first group of PCA data threshold values to obtain three-dimensional space morphological characteristics of the karst cave;

the module for depicting the rock property boundary filled in the cave is configured to depict the rock property boundary filled in the cave: analyzing the characteristic parameters sensitive to the reservoir, which correspond to the karst cave data points, by a PCA method to obtain a second group of PCA data sensitive to the response of the filler;

generating a second set of high-precision characteristic value simulation result data volume through the function as the interwell reservoir parameter simulation module based on the second set of PCA data;

performing intersection analysis based on the second set of PCA data sensitive to the filler response to obtain a second set of geologic body-like discrete data point distribution diagram;

obtaining an interpretation conclusion of a second group of geologic body-like discrete data point distribution diagrams based on the lithological information and the physical property information, distinguishing filler types, further constructing a second group of PCA data intersection layout, and obtaining a second group of PCA data threshold values required by the division of the filler types in the karst cave;

performing intersection analysis on the second group of high-precision characteristic value simulation result data volumes based on the second group of PCA data threshold values to obtain three-dimensional space morphological characteristics of different filler types in the karst cave;

the paleo-karst cavern structure and fill description module is configured to: and engraving the development characteristics of the ancient karst cave space distribution and the internal filling of different filler types by adopting a lithologic shielding technology and a three-dimensional engraving technology based on the three-dimensional space morphological characteristics of the karst cave and the three-dimensional space morphological characteristics of different filler types in the karst cave.

Technical Field

The invention belongs to the field of data identification and record carrier processing, and particularly relates to a method and a system for evaluating deep paleo-karst reservoir filling characteristics by well-seismic combination.

Background

The deep layer of the oil-gas-containing basin (the buried depth is more than 4500 m) contains rich oil-gas resources, and the deep layer of the oil-gas-containing basin becomes an important field for exploration in the oil industry. Carbonate reservoirs account for 52% of the worldwide proven oil reserves and 60% of the production, which are the main targets of deep-layer oil and gas exploration.

The paleo-karst reservoir has stronger heterogeneity due to the complex lithology, the variable diagenesis and random fractures of the carbonate rock. The sparse logging curve has high longitudinal resolution, can only reflect stratum information (generally less than 10 m) of limited depth near a shaft, and the paleo-karst reservoir far away from the shaft still has great change. Therefore, the detection techniques and experience obtained for those shallow clastic and porous carbonate reservoirs cannot be directly used for the detection of deep paleo-karst reservoirs. However, this is often overlooked by those skilled in the art, resulting in poor evaluation of deep paleo-karst reservoirs.

The current evaluation of the filling characteristics of the paleo-karst reservoir has the following problems: as the ancient karst cave type reservoir is buried deeply, the seismic wave energy absorption attenuation effect is strong, so that the seismic data dominant frequency is obviously reduced (less than 25 Hz), the resolution is greatly reduced (generally 10 m), and the signal-to-noise ratio in a deep time window is too low. And drilling, well logging and core samples show that most of caverns are not more than 5 meters in height, and the resolution of the existing three-dimensional seismic data is difficult to achieve the accurate identification effect. In addition, compared with a logging curve, the lithological physical property information content contained in the seismic attribute is low, and the paleo-karst reservoir filling characteristics cannot be effectively identified.

The key problem of the research of the paleo-karst reservoir stratum is how to accurately know the structure of the deep-buried paleo-karst system, how to effectively identify the position of the paleo-karst system by utilizing geophysical data and how to accurately describe the filling characteristics in the paleo-karst cave system.

Disclosure of Invention

The method aims to solve the technical problem of evaluation of filling characteristics of the deep paleo-karst reservoir, namely the existing method for analyzing the filling characteristics of the cave-type reservoir for petroleum exploration only takes geological data as constraint on a random simulation conclusion, does not consider contribution of seismic information to well logging sample optimization, and adopts all well logging characteristic parameters to describe the reservoir too much, so that an inversion result is directly controlled by inter-well interpolation in a high-frequency section, the transverse prediction randomness is strong, and the prediction accuracy of the deep paleo-karst reservoir with rapidly changing transverse filling characteristics is poor. The invention provides a well-seismic joint evaluation method for filling characteristics of deep paleo-karst reservoirs, which comprises the following steps:

step S100, obtaining original geophysical logging data: obtaining raw logging data for each sample well by a logging device, comprising: measuring the natural potential SP of each sample well through a measuring electrode, measuring the natural gamma GR of each sample well through a natural gamma underground device and a natural gamma ground instrument, and obtaining the well diameter CAL of each sample well through a well diameter arm; obtaining resistivity curve data by conventional logging equipment: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and physical property characterization curve data: a compensated neutron CNL, a compensated acoustic curve AC and a density curve DEN;

obtaining determined lithology information and physical property information of individual depth sections based on imaging logging information, drilling information, logging information and core information, and further determining depth data of a target horizon marker layer;

step S200, acquiring seismic data, acquiring original seismic wave reflection signal data through a seismic wave excitation device and a receiving device, and acquiring isochronous three-dimensional spread of a target layer position mark layer according to the waveform of the original seismic wave reflection signal data;

step S300, preprocessing original geophysical logging data: drawing logging curve data based on all original logging data of the sample well, and performing abnormal value processing and standardization processing to obtain standardized logging curve data;

step S400, seismic data preprocessing: based on the seismic wave reflection signal data, obtaining a high-precision three-dimensional seismic amplitude data volume through mixed phase wavelet deconvolution and diffusion filtering;

step S500, well seismic calibration and characteristic parameter selection: acquiring a wave impedance curve of a sample well based on a compensation acoustic curve AC and a density curve DEN in the standardized logging curve data, further calculating a reflection coefficient curve, acquiring the preferred frequency of a Rake wavelet to keep the same as the main frequency of a high-precision three-dimensional seismic amplitude data body, performing convolution operation on the Rake wavelet and the reflection coefficient curve to obtain a synthetic seismic record, comparing the depth data of a target layer position mark layer with the isochronous three-dimensional spread of the target layer position mark layer for well seismic calibration, calculating the correlation between the synthetic seismic record and a well-side seismic channel waveform, judging that a well seismic result is qualified when the correlation is greater than or equal to a preset first threshold, and acquiring the time-depth conversion relation between the logging curve data and the seismic record and characteristic parameters sensitive to a reservoir;

the characteristic parameter sensitive to the reservoir is obtained by the method comprising the following steps:

drawing a histogram by responding to logging parameters generated by different geologic bodies beside a well, and selecting the standardized logging curve data as characteristic parameters sensitive to a reservoir when the numerical value of the standardized logging curve data can distinguish data points above a second threshold value preset by different logging interpretation conclusions; the characteristic parameters sensitive to the reservoir at least comprise wave impedance IMP, well diameter CAL, natural gamma GR, natural potential SP, resistivity curve data: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and the data of the physical property characterization curve: one or more of a compensated neutron CNL, a compensated acoustic curve AC, a density curve DEN;

step S600, obtaining a first group of PCA data: preferably selecting a preset number of pieces of standardized well logging curve data from the characteristic parameters sensitive to the reservoir, and analyzing the characteristic parameters sensitive to the reservoir by a PCA method to obtain a first group of PCA data sensitive to the reservoir;

step S700, constructing an isochronous trellis model: constructing an isochronous trellis model based on sedimentary stratum rules reflected by the high-precision three-dimensional seismic amplitude data volume and the time-depth conversion relation;

step S800, interwell reservoir parameter simulation: determining an optimal sample number parameter capable of reflecting the overall geological condition based on a first group of PCA data of each sample well, selecting the first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample number parameter to construct an initial model, continuously correcting parameters of the initial model, and outputting a first group of high-precision characteristic value simulation result data volumes which are data volumes corresponding to the first group of PCA data one by one;

step S900, depicting the boundary of the karst cave system between wells: performing intersection analysis based on a first group of PCA data of the sample well to obtain a first group of geologic body discrete data point distribution maps;

obtaining an interpretation conclusion of a first group of geologic body discrete data point distribution diagrams based on the lithological information and the physical property information, dividing a karst cave and surrounding rocks, further constructing a first group of PCA data intersection layouts, and obtaining a first group of PCA data threshold values required by the division of the karst cave;

performing intersection analysis on the first group of high-precision characteristic value simulation result data volumes based on the first group of PCA data threshold values to obtain three-dimensional space morphological characteristics of the karst cave;

step S1000, depicting a lithologic property boundary filled in the cave: analyzing the characteristic parameters sensitive to the reservoir corresponding to the karst cave data points by a PCA method to obtain a second group of PCA data sensitive to the response of the filling;

generating a second group of high-precision characteristic value simulation result data volume by the method of S800 based on the second group of PCA data;

performing intersection analysis based on the second group of PCA data sensitive to the filler response to obtain a second group of geologic body discrete data point distribution diagrams;

obtaining an interpretation conclusion of a second group of geologic body discrete data point distribution diagrams of various types based on the lithological information and the physical property information, distinguishing filler types, further constructing a second group of PCA data intersection layout, and obtaining a second group of PCA data threshold values required by the division of the filler types in the karst cave;

performing intersection analysis on the second group of high-precision characteristic value simulation result data volumes based on the second group of PCA data threshold values to obtain three-dimensional space morphological characteristics of different filler types in the karst cave;

step S1100, archaeological cavern structure and filling description: and engraving the development characteristics of the ancient karst cave space distribution and the internal filling of different filler types by adopting a lithologic shielding technology and a three-dimensional engraving technology based on the three-dimensional space morphological characteristics of the karst cave and the three-dimensional space morphological characteristics of different filler types in the karst cave.

In some preferred embodiments, the measuring the natural potential SP of each sample well by the measuring electrode and the measuring the natural gamma GR of each sample well by the natural gamma downhole device and the natural gamma surface instrument specifically include:

measuring the natural potential SP of each sample well through the measuring electrode:

arranging a measuring electrode N on the ground, and arranging a measuring motor M underground through a cable;

lifting the measuring electrode M along the well axis to measure the change of the natural potential along with the well depth;

the calculation method of the natural potential value comprises the following steps:

wherein the content of the first and second substances,the total natural potential is the total natural potential,in order to be a diffusion potential coefficient,in order to obtain a diffusion adsorption potential coefficient,the resistivity value of the slurry filtrate is shown as,is the formation water resistivity value;

the natural gamma GR of each sample well is measured through the natural gamma underground device and the natural gamma ground instrument:

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

acquiring natural gamma rays through a detector, converting the natural gamma rays into electric pulse signals, and amplifying the electric pulse signals through an amplifier;

the ground instrument converts the electric pulse count formed every minute into a potential difference for recording.

In some preferred embodiments, the PCA assay specifically comprises:

step B100, calculating a sample mean value:

wherein the content of the first and second substances,the total number of samples of the curve is shown,ia sample count parameter is represented that is representative of,is shown asiThe characteristic parameter data of each sample is obtained,represents the sample mean;

step B200, for each stripThe characteristic parameter curve is normalized to obtain normalized sample characteristic parameter data：

Wherein the content of the first and second substances,xrepresenting sample feature parameter sample point data;

step B300, based on the normalized sample characteristic parameter dataCalculating the covariance of the characteristic parameters of the sample：

Wherein the content of the first and second substances,represents the c-th X sample characteristic parameter data,representing the c sample characteristic parameter data;

step B400, constructing a sample characteristic parameter covariance matrix based on the sample characteristic parameter covariance:

z represents sample characteristic parameter data;

Includedpthe covariance matrix of the seed sample characteristic parameter data is:

wherein the content of the first and second substances,representing a plurality of sample characteristic parameter data,representing the characteristic parameter data of the 1 st sample and thepCovariance of seed sample characteristic parameter data;

step B500, based on the inclusionpThe covariance matrix of the sample characteristic parameter data is planted, and the eigenvalue of the covariance matrix is obtainedAndpa feature vector，，；

Eigenvalue and eigenvector computations based on the covariance matrix includepPCA data for species eigenvalues:

wherein the content of the first and second substances,represents the jth PCA data vector and,representing the jth sample characteristic parameter data vector,represents the jth sample bitA feature vector of the feature parameter data,j 1,2, pis a counting parameter;

step B600, based on includingpAnd (4) sequencing the corresponding eigenvalues from large to small according to the PCA data of the variety eigenvalues to obtain the eigenvalues sequenced from large to small，…，Calculating the variance cumulative contribution:

wherein the content of the first and second substances,representing the jth characteristic value sorted from large to small,to accumulate the first k covariance accumulated contributions,；

when k is increased from 1, the k is selectedAnd when the first time is greater than a preset third threshold value, corresponding PCA data vectors corresponding to the first k covariances are PCA data obtained by PCA analysis in the current round.

In some preferred embodiments, the step S400 specifically includes:

step S410, based on the seismic wave reflection signal data, expressing a frequency domain seismic record convolution model as:

wherein the content of the first and second substances,representing the fourier transformed seismic records,representing the wavelet after the fourier transform,a frequency spectrum representing the fourier transformed reflection coefficient,represents angular frequency;

step S420, logarithm taking two sides of the equation of the frequency domain seismic record convolution model is converted into a linear system, and a linear seismic record convolution model is obtained:

step S430, performing inverse Fourier transform on the linear seismic record convolution model to obtain a cepstrum sequence:

wherein the content of the first and second substances,representing a repeating spectral sequence of seismic waveform recordings,a cepstrum sequence representing seismic wavelets,a repeating spectral sequence representing the reflection coefficients of the formation,representing a seismic waveform recording time;

step S440, based on the cepstrum sequence, performing wavelet and reflection coefficient separation through a low-pass filter, and extracting a wavelet amplitude spectrum;

step S450, obtaining the amplitude spectrum of the simulated seismic wavelet by a least square method:

wherein, least square method is used for fitting parametersIs a constant number of times, and is,a representation of the amplitude spectrum of the wavelet is obtained,anda polynomial expression representing f, which represents the frequency of the seismic wave;

step S460, obtaining a wavelet maximum phase component and a wavelet minimum phase component based on the simulated seismic wavelet amplitude spectrum;

wavelet setting deviceHas a maximum phase component ofThe minimum phase component isWavelet of fundamental waveComprises the following steps:

the magnitude spectrum is represented in the cepstrum as:

wherein the repetition spectrum of the amplitude spectrumSymmetrically displayed on the positive and negative axes of the match score,for maximum phase component of seismic waveletThe cepstrum of the corresponding minimum phase function,for minimum phase component of seismic waveletsThe corresponding cepstrum of the maximum phase function;

step S470, determining a group of mixed phase wavelet sets with the same amplitude spectrum based on the cepstrum in the amplitude spectrum, continuously adjusting the parameters of Shu' S wavelets, maintaining low frequency, expanding high frequency and properly improving dominant frequency to construct an expected output wavelet form, searching for an optimal balance point between resolution and fidelity by taking a signal-to-noise ratio spectrum as a reference under the control of a well curve, and obtaining waveform data after shaping;

step S480, constructing a tensor diffusion model based on the shaped waveform data:

wherein the content of the first and second substances,it is shown that the time of diffusion,denotes a divergence operator, D denotes a tensor-type diffusion coefficient of the diffusion filter, U denotes a diffusion filtering result,to representThe result of diffusion filtering when =0,to representThe waveform data after the shaping at that time is used as an initial condition of the tensor diffusion model,a gradient representing a result of the diffusion filtering;

constructing a gradient structure tensor based on the tensor diffusion model:

where, U represents the result of the diffusion filtering,representing a gradient vector tensor product;

to representHas a dimension ofThe gaussian function of (d) is:

wherein, r represents the calculated radius,

the eigenvectors of the structure tensor are:

wherein the content of the first and second substances,、andthe 3 eigenvectors, expressed as gradient structure tensors, can be considered as local orthogonal coordinate systems,pointing in the direction of the gradient of the seismic signal,andthe planes of composition are parallel to the local structural features of the seismic signal,、andare respectively connected with、Andcorresponding three characteristic values;

step S490, calculating a linear structure confidence metric, a planar structure confidence metric, and a diffusion tensor, respectively, based on the feature vectors of the structure tensor;

the presence structure confidence measureComprises the following steps:

the planar structure confidence measureComprises the following steps:

the diffusion tensor D is:

wherein the content of the first and second substances,、andthree non-negative eigenvalues representing the diffusion tensor, each representing a diffusion filter edge、Andthe filter strengths of the three characteristic directions;

and S4100, repeating the steps of S480-S490 until a preset iteration number is reached, and obtaining a diffusion filtering result, namely the high-precision three-dimensional seismic amplitude data volume.

In some preferred embodiments, the method for calculating the amplitude spectrum of the simulated seismic wavelet comprises:

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

obtaining parameters by fitting the maximum value of the seismic signal amplitude spectrum and the simulated seismic wavelet amplitude spectrum in a least square method modeAndobtaining the corresponding frequency amplitude value of the fitted maximum value by the coefficients of the polynomial;

dividing the maximum value of the seismic signal amplitude spectrum by the fitted amplitude value of the corresponding frequency, and further using a quotient fitting polynomialThe coefficient of (a).

In some preferred embodiments, the time-depth conversion relationship of the well log data to the time-depth conversion relationship of the seismic record is:

wherein the content of the first and second substances,represents the seismic signal time corresponding to the initial depth of acoustic logging,the time difference of the sound wave is represented,indicating the sequence number representing the time operation between each sample point,representing the interval of sampling of the log data,representing the seismic wave two-way travel time.

In some preferred embodiments, the step S800 specifically includes, step S810-step S890:

step S810, selecting a sample well as a reference target well, and setting an initial sample number parameter to be 1;

step S820, selecting a first group of PCA data of the sample wells with the quantity being the sample number parameter and a first group of PCA data of the reference target well mark according to the waveform similarity principle to perform correlation analysis to obtain a first group of PCA data correlation values of the sample number parameter-reference target well;

step S830, increasing the sample number parameters 1 by 1, repeating the method of the step S720 to obtain a first group of PCA data correlation values of the sample number parameter-reference target well corresponding to each sample number parameter, connecting the first group of PCA data correlation values of all the sample number parameters-reference target wells, and obtaining a correlation curve of the first group of PCA data correlation of the reference wells along with the change of the sample number parameters;

step S840, selecting another sample well as a reference target well, repeating the steps S810-S830, obtaining a correlation curve of the first group of PCA data correlations of the plurality of reference wells along with the change of the sample number parameters, fitting the correlation curve of the first group of PCA data waveform correlations of all the reference wells along with the change of the sample number parameters into an overall correlation curve, selecting an inflection point at which the correlation in the overall correlation curve rises along with the increase of the sample number parameters and finally keeps stable, and determining the optimal sample number parameters;

step S850, based on the high-precision three-dimensional seismic amplitude data volume and the isochronous lattice model, calculating the waveform correlation between the point to be detected and the sample well position, sorting the waveform correlation from large to small, and selecting a first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample quantity parameter strip; constructing an initial model based on the sample well corresponding to the seismic waveform characteristic data of the sample well with the highest correlation through an inter-well characteristic parameter interpolation mode;

step S860, based on the initial model, selecting a first group of PCA data of the sample well with the optimal sample number parameter bar and the highest seismic waveform correlation degree as prior information;

step S870, performing matched filtering on the initial model and the prior information to obtain a maximum likelihood function;

step S880, based on the maximum likelihood function and the prior information, obtaining the posterior probability statistical distribution density under a Bayes framework, and sampling the posterior probability statistical distribution density to obtain a target function;

step S890, using the target 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 of the target function when the maximum value is taken as the random realization, taking the average value of multiple random realizations as the expected value output, and taking the expected value output as the high-precision characteristic value simulation result data volume; and parameters in the high-precision characteristic value simulation result data body correspond to the first group of PCA data one by one.

In some preferred embodiments, the step S880 specifically includes:

step S881, using white noise to satisfy the rule of gaussian distribution, representing the parameters of the high-precision eigenvalue simulation result data volume as:

y represents parameters of a logging curve high-precision characteristic value simulation result data volume, X represents actual characteristic parameter values of an underground stratum to be solved, and N represents random noise;

step S882, becauseAlso satisfying a gaussian distribution, the initial objective function can be determined as:

wherein the content of the first and second substances,a function relating to a posteriori information is represented,showing that the characteristic curve of the sample well is matched and filtered by selecting the sample well based on the optimal sample number, the posterior probability statistical distribution density is obtained, and the expected value of the characteristic parameter is further calculated,a covariance representing white noise;

step S883, based on the initial objective function, introducing prior information into the objective function through maximum posterior estimation, and obtaining a stable objective function as:

wherein the content of the first and second substances,representing the characteristic parameter to be simulated,representing functions related to prior information such as geological and well log data,representation for coordinationAndthe smoothing parameters of the mutual influence between them.

In some preferred embodiments, the step S890 includes the specific steps of:

step S891, setting M as a target space, n as the total number of samples, and M as the number of samples when the Markov chain tends to be stable;

step S892, presetting a Markov chain to make the Markov chain converge to a stable distribution;

step S893, starting from a certain point in MStarting from, sampling simulation is performed through a Markov chain to generate a point sequence:；

step S894, functionThe expected estimate of (c) is:

wherein n represents the total number of generated samples, m represents the number of samples when the Markov chain reaches a plateau, and k represents an accumulation parameter;

step S895, selecting a transfer functionAnd an initial valueIf the parameter value at the beginning of the ith iteration isThen, the ith iteration process is:

fromExtract an alternative valueCalculating alternative valuesProbability of acceptance of：

Step S896, in order toIs arranged atBy probabilityIs arranged at；

Step S897, continuously disturbing the parameters of the initial model,repeating the steps S892 to S896 until the preset iteration number n is reached, and obtaining posterior samplesAnd further calculating each order matrix of posterior distribution to obtain expected output values, and outputting the expected values as a first group of high-precision characteristic value simulation result data volume.

In some preferred embodiments, the step S300 includes:

step S310, drawing original logging curve data based on the original logging data;

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

and S330, superposing single logging curve histogram data of all sample well positions in the work area based on the logging curve data without outliers, and obtaining standardized logging curve data by integrating threshold values.

In another aspect of the present invention, a system for evaluating deep paleo-karst reservoir filling characteristics by well-seismic combination is provided, the system comprising: 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 seismic calibration and characteristic parameter selection module, a first group of PCA data acquisition modules, an isochronous grid model construction module, an interwell reservoir parameter simulation module, an interwell karst cave system boundary description module, a lithologic property boundary description module for describing filling inside a cave and an ancient karst cave structure and filling description module;

the original geophysical logging data acquisition module is configured to acquire original logging data of each sample well through logging equipment, and comprises: measuring the natural potential SP of each sample well through a measuring electrode, measuring the natural gamma GR of each sample well through a natural gamma underground device and a natural gamma ground instrument, and obtaining the well diameter CAL of each sample well through a well diameter arm; obtaining resistivity curve data by conventional logging equipment: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and physical property characterization curve data: a compensated neutron CNL, a compensated acoustic curve AC and a density curve DEN;

obtaining determined lithology information and physical property information of individual depth sections based on imaging logging information, drilling information, logging information and core information, and further determining depth data of a target horizon marker layer;

the seismic data acquisition module is configured to acquire original seismic wave reflection signal data through a seismic wave excitation device and a receiving device, and acquire isochronous three-dimensional spread of a target layer position mark layer according to the waveform of the original seismic wave reflection signal data;

the original geophysical logging data preprocessing module is configured to draw logging curve data based on original logging data of the sample well, perform abnormal value processing and standardization processing, and obtain standardized logging curve data;

the seismic data preprocessing module is configured to obtain a high-precision three-dimensional seismic amplitude data volume through mixed phase wavelet deconvolution and diffusion filtering based on the seismic wave reflection signal data;

the well-seismic calibration and characteristic parameter selection module is configured to perform well-seismic calibration and characteristic parameter selection: acquiring a wave impedance curve of a sample well based on a compensation acoustic curve and a density curve DEN in the standardized logging curve data, further calculating a reflection coefficient curve, acquiring the preferred frequency of a Rake wavelet to keep the same as the main frequency of a high-precision three-dimensional seismic amplitude data body, performing convolution operation on the Rake wavelet and the reflection coefficient curve to obtain a synthetic seismic record, comparing the depth data of a target layer position mark layer with the isochronous three-dimensional spread of the target layer position mark layer for well seismic calibration, calculating the correlation between the synthetic seismic record and a well-side seismic channel waveform, and judging that a well seismic calibration result is qualified when the correlation is greater than or equal to a preset first threshold value, thereby acquiring the time-depth conversion relation between the logging curve data and the seismic record and characteristic parameters sensitive to a reservoir;

the characteristic parameter sensitive to the reservoir is obtained by the method comprising the following steps:

drawing a histogram by responding to logging parameters generated by different geologic bodies beside a well, and selecting a certain standardized logging curve data as a characteristic parameter sensitive to a reservoir when the data point above a second threshold preset by different logging interpretation conclusions can be distinguished by the standardized logging curve data; the characteristic parameters sensitive to the reservoir at least comprise wave impedance IMP, well diameter CAL, natural gamma GR, natural potential SP, resistivity curve data: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and the data of the physical property characterization curve: one or more of a compensated neutron CNL, a compensated acoustic curve AC, a density curve DEN;

the first group of PCA data acquisition modules are configured to preferably select a preset number of pieces of standardized well logging curve data from the characteristic parameters sensitive to the reservoir, and perform dimensionality reduction on the characteristic parameters sensitive to the reservoir through PCA analysis to obtain a first group of PCA data;

the isochronous trellis model building module is configured to build an isochronous trellis model based on sedimentary stratum rules reflected by the high-precision three-dimensional seismic amplitude data volume and the time-depth conversion relation;

the interwell reservoir parameter simulation module is configured to determine an optimal sample number parameter capable of reflecting the overall geological condition based on a first group of PCA data of each sample well, select a first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample number parameter to construct an initial model, continuously correct parameters of the initial model, and output 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 are data bodies corresponding to the first group of PCA data one by one;

the boundary delineation module of the inter-well karst cave system is configured to perform intersection analysis based on a first group of PCA data of the sample well to obtain a first group of discrete data point distribution maps of various geologic bodies;

obtaining an interpretation conclusion of a first group of geologic body discrete data point distribution diagrams based on the lithological information and the physical property information, dividing a karst cave and surrounding rocks, further constructing a first group of PCA data intersection layouts, and obtaining a first group of PCA data threshold values required by the division of the karst cave;

performing intersection analysis on the first group of high-precision characteristic value simulation result data volumes based on the first group of PCA data threshold values to obtain three-dimensional space morphological characteristics of the karst cave;

the module for depicting the rock property boundary filled in the cave is configured to depict the rock property boundary filled in the cave: analyzing the characteristic parameters sensitive to the reservoir, which correspond to the karst cave data points, by a PCA method to obtain a second group of PCA data sensitive to the response of the filler;

generating a second set of high-precision characteristic value simulation result data volume through the function as the interwell reservoir parameter simulation module based on the second set of PCA data;

performing intersection analysis based on the second group of PCA data sensitive to the filler response to obtain a second group of geologic body discrete data point distribution diagrams;

obtaining an interpretation conclusion of a second group of geologic body discrete data point distribution diagrams of various types based on the lithological information and the physical property information, distinguishing filler types, further constructing a second group of PCA data intersection layout, and obtaining a second group of PCA data threshold values required by the division of the filler types in the karst cave;

performing intersection analysis on the second group of high-precision characteristic value simulation result data volumes based on the second group of PCA data threshold values to obtain three-dimensional space morphological characteristics of different filler types in the karst cave;

the paleo-karst cavern structure and fill description module is configured to: and engraving the development characteristics of the ancient karst cave space distribution and the internal filling of different filler types by adopting a lithologic shielding technology and a three-dimensional engraving technology based on the three-dimensional space morphological characteristics of the karst cave and the three-dimensional space morphological characteristics of different filler types in the karst cave.

In a third aspect of the present invention, an electronic device is provided, including:

at least one processor;

and at least one memory communicatively coupled to the processor; wherein the content of the first and second substances,

the memory stores instructions executable by the processor for performing the method for joint well-seismic evaluation of deep paleo-reservoir packing characteristics described above.

In a fourth aspect of the present invention, a computer-readable storage medium is provided, which stores computer instructions for execution by the computer to implement the method for joint well-seismic evaluation of deep paleo-karst reservoir filling characteristics as described above.

The invention has the beneficial effects that:

the method fully utilizes logging and seismic data, adopts a machine learning algorithm, takes a waveform indication simulation method as a core, and utilizes a principal component analysis algorithm (PCA algorithm) to establish sensitive parameters so as to realize the accurate identification of the logging data on the reservoir structure and filling; and then, establishing a three-dimensional seismic simulation data volume by using the sensitive parameters and the three-dimensional seismic data, extrapolating a high-resolution logging interpretation result to a three-dimensional large-range area between wells, and finally evaluating the filling characteristics of the deep paleo-karst reservoir.

The invention carries out three-dimensional parameter simulation on well-seismic combination on the basis of an initial isochronous geological framework model by utilizing the relevance among high-precision three-dimensional seismic waveform data, transforms a coordinate system of a plurality of preprocessed logging curves by applying two-time PCA analysis (PCA), and carries out dimensionality reduction analysis on lithological and physical property information reflected by a plurality of conventional logging curve data, so that a small amount of PCA data has the characterization capability of the plurality of logging curves on a reservoir.

The karst cave longitudinal and transverse structure and the internal filling feature description conclusion of the karst cave under the PCA intersection analysis constraint are obtained by well distinguishing the karst cave and the internal filling data points of the karst cave through a two-dimensional intersection chart constructed based on the PCA data so as to obtain a karst reservoir identification intersection chart, and applying reservoir division threshold value information embodied by the chart to three-dimensional parameter simulation data. The method improves the precision of the depiction, and can predict the reservoir with quick transverse change. And a reliable theoretical basis is provided for large-scale evaluation of the storage performance of the carbonate rock ancient karst cave and oil and gas migration.

According to the invention, a mapping relation between the geophysical detection method is established by a Markov chain sampling criterion and a Monte Carlo estimation method under the guidance of a Bayesian theory, so that the effect of identifying the development characteristics of the carbonate rock ancient karst cave type reservoir in the complicated basin in a large range is achieved.

Drawings

Other features, objects and advantages of the present application will become more apparent upon reading of the following detailed description of non-limiting embodiments thereof, made with reference to the accompanying drawings in which:

FIG. 1 is a schematic flow diagram of an embodiment of a method and system for joint well-seismic evaluation of deep paleo-karst reservoir filling characteristics;

FIG. 2 is a schematic diagram of a log with outliers removed according to the present invention;

FIG. 3 is a schematic diagram of the present invention showing the superposition of all curves during the normalization process;

FIG. 4 is a plot of a single well seismic calibration in an embodiment of the present invention;

FIG. 5 is a graph of an isochronous trellis model below the T74 marker level in an embodiment of the present invention;

FIG. 6 is a schematic diagram of a PCA for fitting a correlation curve of waveform correlations for all reference wells as a function of a sample number parameter to an overall correlation curve in an embodiment of the present invention;

FIG. 7 is a schematic diagram illustrating a first set of PCA data intersection analysis capturing karst cave data points and gating threshold values in accordance with an embodiment of the present invention;

FIG. 8 is a cross-sectional view of the three-dimensional distribution and structure of an archaeological cavern in an embodiment of the present invention;

FIG. 9 is a graphical illustration of a second set of PCA data intersection analysis for dividing the interior fill of a karst cavity and for gating threshold values in accordance with an embodiment of the present invention;

fig. 10 is a cross-section depicting the development characteristics of the filler in an archaeological cavern in accordance with an embodiment of the present invention.

Detailed Description

The present application will be described in further detail with reference to the following drawings and examples. It is to be understood that the specific embodiments described herein are merely illustrative of the relevant invention and not restrictive of the invention. It should be noted that, for convenience of description, only the portions related to the related invention are shown in the drawings.

It should be noted that the embodiments and features of the embodiments in the present application may be combined with each other without conflict. The present application will be described in detail below with reference to the embodiments with reference to the attached drawings.

The method fully utilizes logging and seismic data, adopts a machine learning algorithm, takes a waveform indication simulation method as a core, and utilizes a principal component analysis algorithm (PCA algorithm) to establish sensitive parameters so as to realize the accurate identification of the logging data on the reservoir structure and filling; and then, establishing a three-dimensional seismic simulation data volume by using the sensitive parameters and the three-dimensional seismic data, extrapolating a high-resolution logging interpretation result to a three-dimensional large-range area between wells, and finally evaluating the filling characteristics of the deep paleo-karst reservoir.

In order to more clearly describe the method for intelligently dividing the filling characteristics of the paleo-karst cavern reservoir, the following describes the steps in the embodiment of the invention in detail with reference to fig. 1.

The method and the system for evaluating the filling characteristics of the deep paleo-karst reservoir by combining well and seismic in the first embodiment of the invention comprise the steps S100-S1100, and the steps are described in detail as follows:

the practical seismic data statistical analysis shows that the seismic waveform characteristics generated by the ancient karst cave generally consist of a plurality of groups of wave troughs and wave crests, and research shows that the characteristic reflection is caused by seismic wave interference. And the amplitude of the reflected wave is related to the cavern filling combination. Xu et al (2016) can also find that ancient karst cave type objects with different shapes, thicknesses and combination relations can generate various seismic waveform reflection characteristics through physical model simulation, namely the reflection form characteristics are related to the cave diameter and the cave form and the distribution rule. The scholars further discovered that within the same facies of the isochronous stratigraphic framework, similarity of waveform features may represent lithology combinations with certain correlations. Therefore, the method carries out the partition inversion of the isochronous interface and the simulation of the characteristic parameters based on the idea of waveform indication, utilizes the transverse change information of the seismic waveform, better embodies the constraint of the sedimentary environment and better accords with the sedimentary geological rule.

Step S100, obtaining original geophysical logging data: obtaining raw logging data for each sample well by a logging device, comprising: measuring the natural potential SP of each sample well through a measuring electrode, measuring the natural gamma GR of each sample well through a natural gamma underground device and a natural gamma ground instrument, and obtaining the well diameter CAL of each sample well through a well diameter arm; obtaining resistivity curve data by conventional logging equipment: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and physical property characterization curve data: a compensated neutron CNL, a compensated acoustic curve AC and a density curve DEN;

obtaining determined lithology information and physical property information of individual depth sections based on imaging logging information, drilling information, logging information and core information, and further determining depth data of a target horizon marker layer;

the natural potential SP of each sample well is measured through the measuring electrode, and the natural gamma GR of each sample well is measured through the natural gamma underground device and the natural gamma ground instrument, and the method specifically comprises the following steps:

measuring the natural potential SP of each sample well through the measuring electrode:

arranging a measuring electrode N on the ground, and arranging a measuring motor M underground through a cable;

lifting the measuring electrode M along the well axis to measure the change of the natural potential along with the well depth;

the calculation method of the natural potential value comprises the following steps:

wherein the content of the first and second substances,the total natural potential is the total natural potential,in order to be a diffusion potential coefficient,in order to obtain a diffusion adsorption potential coefficient,the resistivity value of the slurry filtrate is shown as,is the formation water resistivity value;

the natural gamma GR of each sample well is measured through the natural gamma underground device and the natural gamma ground instrument:

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

acquiring natural gamma rays through a detector, converting the natural gamma rays into electric pulse signals, and amplifying the electric pulse signals through an amplifier;

the ground instrument converts the electric pulse count formed every minute into a potential difference for recording.

The dual lateral electrodes are composed of a main electrode, a monitor electrode, a ring-shaped shield electrode and a column electrode. Main electrodeCentrally and vertically symmetrically distributed monitoring electrodes、And、and a ring-shaped shield electrode、In a、Two columnar electrodes are added at the symmetrical positions of the outer side of the electrode. Two columnar electrodes in the deep lateral electrode system are shielding electrodes、(ii) a Two columnar electrodes in the shallow lateral electrode system are loop electrodesAnd. A contrast electrode N and a return electrode B of a deep lateral electrode system are arranged at the far position of the electrode system; main electrode of micro-spherical focusing electrode systemIs a rectangular sheet electrode, and the rectangular frame electrodes which are sequentially outward are measuring electrodesAuxiliary electrodeMonitor electrodeAndand the loop electrode B is arranged on the instrument shell or the polar plate support frame. Measuring the potential difference change between any monitoring electrode and the contrast electrode, namely reflecting the change of the medium resistivity, wherein the apparent resistivity expression is as follows:

in the formula (I), the compound is shown in the specification,the electrode coefficient is different from deep lateral logging, shallow lateral logging and micro-spherical focusing logging;is the potential of the monitor electrode M;is the main current.

The transmitting transducer crystals of the downhole tool vibrate causing particles of the surrounding medium to vibrate, producing acoustic waves that propagate into the mud and rock formation in the well. The receiving transducer R, R2 may be used to receive the slip waves sequentially in the well and record the time differenceAnd thereby measuring the acoustic velocity of the formation. Arrive atAndrespectively at the times ofAndthen the time difference of arrival at the two receiving transducersComprises the following steps:

in the formula (I), the compound is shown in the specification,is the mud sound velocity;at the speed of sound of the earth formation

The density logging instrument comprises a gamma source, two detectors for receiving gamma rays, namely a long source distance detector and a short source distance detector. They are mounted on a skid plate and are pushed against the borehole wall during logging. An auxiliary electronic circuit is arranged above the downhole instrument. Generally 137Cs is used as the gamma source, which emits gamma rays of moderate energy (0.661MeV) and which only produce compton scattering and photoelectric effects when irradiated on a substance. The density of the formation is different, and the scattering and absorption capabilities of the gamma photons are different, so that the counting rates of the gamma photons received by the detector are different. The count rate N of a gamma photon with a known distance L is known as:

if only Compton scattering is present, thenNamely, the Compton scattering absorption coefficient is obtained by transformation:

where C is a constant associated with the formation and L is the distance of the receiving source from the gamma source.

After the source distance is selected, the instrument is calibrated in scale to findAnd N, recording the count rate N of the scattered gamma photons to obtain the density of the stratum

Compensated neutron logging by a neutron emitting sourceAnd two receiving detectors. By counting two detectors, two counting rates are obtainedAndtaking the ratio may reflect the hydrogen content of the formation:

in the formula (I), the compound is shown in the specification,andthe distance from the two detectors to the neutron source;to slow down the length, the formation hydrogen content is reflected.

During the hole diameter measurement, the instrument is lowered to the expected depth of a human well, then the hole diameter arms are opened in a traditional mode, so that four hole diameter legs which are mutually at 90 degrees are outwards stretched under the action of an elastic force , and the tail ends of the four hole diameter legs are tightly attached to the wall of the well. Along with the upward lifting of the instrument, the well diameter arm can expand and contract due to the change of the well diameter and drive the connecting rod to move up and down. The connecting rod is connected with the sliding end of a potentiometer, so that the change of the well diameter can be converted into the change of the resistance. When a current of a certain intensity is applied to the movable resistor, the potential difference between a fixed end and a sliding end of the movable resistor changes along with the change of the resistance value between the fixed end and the sliding end. Thus, measuring this potential difference reflects the borehole diameter:

in the formula (I), the compound is shown in the specification,is the potential difference of a borehole diameter measuring instrument;is an initial value; and c is an instrument constant.

Step S200, acquiring seismic data, acquiring original seismic wave reflection signal data through a seismic wave excitation device and a receiving device, and acquiring isochronous three-dimensional spread of a target layer position mark layer according to the waveform of the original seismic wave reflection signal data;

in the embodiment, nine pieces of conventional logging curve data with the depth range of 5500-5750 m are detected at a work area well bore, and the sampling interval is set to be 0.01 m; the area of a work area is about 27km by using a three-dimensional seismic exploration method and using a seismic wave excitation source and a seismic signal detector²The signal recording time of the three-dimensional seismic data is 4s in a double-travel mode, the time interval of sampling points is 1ms, and the detection depth exceeds 6000 m.

Step S300, preprocessing original geophysical logging data: drawing logging curve data based on all original logging data of the sample well, and performing abnormal value processing and standardization processing to obtain standardized logging curve data;

in this embodiment, the step S300 includes:

step S310, drawing original logging curve data based on the original logging data;

step S320, based on the original logging curve data, removing outliers to obtain logging curve data with the outliers removed; outliers, that is, unreasonable values exist in the data set, a single curve parameter histogram of all the well data is counted, the threshold value of the retention interval is reasonably adjusted to remove the outliers, in this embodiment, the data with the maximum deviation from the median of the first 5 percent is preferably removed, and the logging curve histogram with the outliers removed is shown in fig. 2;

step S330, superposing single logging curve histogram data of all sample well positions in the work area based on the logging curve data without outliers, and obtaining standardized logging curve data by integrating a threshold value; taking the AC normalization process as an example, the normalized well log data is obtained. Due to instrument differences or other factors, conventional logging data among different wells are larger overall or smaller overall, and the curves need to be standardized, and all curves are overlapped as shown in fig. 3.

Step S300, seismic data preprocessing: based on the seismic wave reflection signal data, obtaining a high-precision three-dimensional seismic amplitude data volume through mixed phase wavelet deconvolution and diffusion filtering;

at present, a three-dimensional seismic network measuring technology with a channel interval of 25m multiplied by 25m is widely applied to the field of petroleum exploration, seismic wave reflection signals are received at a sampling interval of 2ms in the vertical direction, and the total sampling time is within 6s so as to detect geological features of intervals with different depths. The method is usually used for target body detection of more than 2m, has higher requirements on seismic data dominant frequency, and should be within the range of 50-60 Hz. When the dominant frequency of the seismic data amplitude data body in the development interval of the ancient karst cave is lower than 50HZ, mixed-phase wavelet deconvolution and maximum stereo deconvolution and diffusion filtering are adopted to carry out frequency extension and noise reduction processing on the three-dimensional seismic data, and a three-dimensional seismic data body with high resolution and high signal-to-noise ratio is obtained.

The mixed-phase wavelet deconvolution is a data processing method for increasing the resolution of seismic signals by widening the effective frequency band on the premise of ensuring that the processed seismic data has higher fidelity, and is equivalent to S410-step S1100 in this embodiment.

In this embodiment, step S400 specifically includes:

step S410, based on the seismic wave reflection signal data, expressing a frequency domain seismic record convolution model as:

wherein the content of the first and second substances,representing the fourier transformed seismic records,representing the wavelet after the fourier transform,a frequency spectrum representing the fourier transformed reflection coefficient,represents angular frequency;

step S420, logarithm taking two sides of the equation of the frequency domain seismic record convolution model is converted into a linear system, and a linear seismic record convolution model is obtained:

step S430, performing inverse Fourier transform on the linear seismic record convolution model to obtain a cepstrum sequence:

wherein the content of the first and second substances,representing a repeating spectral sequence of seismic waveform recordings,a cepstrum sequence representing seismic wavelets,a repeating spectral sequence representing the reflection coefficients of the formation,representing a seismic waveform recording time;

step S440, based on the cepstrum sequence, performing wavelet and reflection coefficient separation through a low-pass filter, and extracting a wavelet amplitude spectrum; the difference in the smoothness of the wavelet and the sequence of reflection coefficients is easily distinguished in the cepstrum: the wavelet energy appears near the origin and the sequence of reflection coefficients is far from the origin. The wavelet in the cepstrum can be separated from the reflection coefficient by using a low-pass filter, so that the purpose of extracting the wavelet amplitude spectrum is achieved.

Step S450, obtaining the amplitude spectrum of the simulated seismic wavelet by a least square method:

wherein, least square method is used for fitting parametersIs a constant number of times, and is,a representation of the amplitude spectrum of the wavelet is obtained,anda polynomial expression representing f, which represents the frequency of the seismic wave;

in this embodiment, the method for calculating the amplitude spectrum of the simulated seismic wavelet includes:

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

obtaining parameters by fitting the maximum value of the seismic signal amplitude spectrum and the simulated seismic wavelet amplitude spectrum in a least square method modeAndthe coefficients of the polynomial to obtain the corresponding frequency of the fitted maximumAn amplitude value;

dividing the maximum value of the seismic signal amplitude spectrum by the fitted amplitude value of the corresponding frequency, and further using a quotient fitting polynomialThe coefficient of (a);

step S460, obtaining a wavelet maximum phase component and a wavelet minimum phase component based on the simulated seismic wavelet amplitude spectrum;

wavelet setting deviceHas a maximum phase component ofThe minimum phase component isWavelet of fundamental waveComprises the following steps:

the magnitude spectrum is represented in the cepstrum as:

wherein the repetition spectrum of the amplitude spectrumSymmetrically displayed on the positive and negative axes of the match score,for maximum phase component of seismic waveletThe cepstrum of the corresponding minimum phase function,for minimum phase component of seismic waveletsThe corresponding cepstrum of the maximum phase function;

step S470, determining a group of mixed phase wavelet sets with the same amplitude spectrum based on the cepstrum in the amplitude spectrum, continuously adjusting the parameters of Shu' S wavelets, maintaining low frequency, expanding high frequency and properly improving dominant frequency to construct an expected output wavelet form, searching for an optimal balance point between resolution and fidelity by taking a signal-to-noise ratio spectrum as a reference under the control of a well curve, and obtaining waveform data after shaping; the wavelet amplitude spectrum is shown in FIG. 4;

after deconvolution of the mixed phase wavelets, the effective frequency band of the seismic data is expanded, and the high-frequency part is reasonably strengthened. The number of the same-phase axes represented on the seismic waveform is increased, the detail change of seismic wave reflection information is reflected more easily, and the consistency of the same reflection wave group waveform is improved in the aspects of amplitude, phase and frequency. On the ancient karst cave seismic response, the 'beaded' reflection characteristic is particularly obvious, the details of the internal form of beads can be clearly displayed, the complex ancient karst cave type reservoir seismic reflection of different structural characteristics and filler combination is represented, and the later fine geological interpretation is facilitated.

Fhemers and Hocker first applied diffusion filtering techniques in seismic data processing interpretation in 2003. The technology not only can effectively suppress noise, but also can keep details in the seismic data as much as possible: such as geologic body edge, fault, unconformity surface, pinch-out and the like, provides reliable seismic data for subsequent seismic interpretation and reservoir prediction work, and greatly improves the success rate of oil and gas exploration and development.

In order to attenuate seismic noise and enhance the diffusion filtering effect of the geological structure characteristics, seeking diffusion tensor is the most key step of the method:

step S480, constructing a tensor diffusion model based on the shaped waveform data:

wherein the content of the first and second substances,it is shown that the time of diffusion,denotes a divergence operator, D denotes a tensor-type diffusion coefficient of the diffusion filter, U denotes a diffusion filtering result,to representThe result of diffusion filtering when =0,to representThe waveform data after the shaping at that time is used as an initial condition of the tensor diffusion model,a gradient representing a result of the diffusion filtering;

constructing a gradient structure tensor based on the tensor diffusion model:

where, U represents the result of the diffusion filtering,representing a gradient vector tensor product;

the expression scale isThe gaussian function of (d) is:

wherein r represents the calculated radius;

the eigenvectors of the structure tensor are:

wherein the content of the first and second substances,、andthe 3 eigenvectors, expressed as gradient structure tensors, can be considered as local orthogonal coordinate systems,pointing in the direction of the gradient of the seismic signal,andthe planes of composition are parallel to the local structural features of the seismic signal,、andare respectively connected with、Andcorresponding three characteristic values;

step S490, calculating a linear structure confidence metric, a planar structure confidence metric, and a diffusion tensor, respectively, based on the feature vectors of the structure tensor;

the presence structure confidence measureComprises the following steps:

the planar structure confidence measureComprises the following steps:

the diffusion tensor D is:

wherein the content of the first and second substances,、andthree non-negative eigenvalues representing the diffusion tensor, each representing a diffusion filter edge、Andthe filter strengths of the three characteristic directions;

and S4100, repeating the steps of S480-S490 until a preset iteration number is reached, and obtaining a diffusion filtering result, namely the high-precision three-dimensional seismic amplitude data volume. The diffusion filtering algorithm reserves the geological characteristics of a beaded reflective ancient karst cave-type reservoir and enhances the imaging capability of seismic data on a target geologic body. Meanwhile, the effects of suppressing noise and improving the transverse continuity of the in-phase axis and the signal-to-noise ratio of the seismic signal are achieved.

Step S500, well seismic calibration and characteristic parameter selection: acquiring a wave impedance curve of a sample well based on a compensation acoustic curve and a density curve DEN in the standardized logging curve data, further calculating a reflection coefficient curve, acquiring the preferred frequency of a Rake wavelet to keep the same as the main frequency of a high-precision three-dimensional seismic amplitude data body, performing convolution operation on the Rake wavelet and the reflection coefficient curve to obtain a synthetic seismic record, comparing the depth data of a target layer position mark layer with the isochronous three-dimensional spread of the target layer position mark layer for well seismic calibration, calculating the correlation between the synthetic seismic record and a well-side seismic channel waveform, and judging that a well seismic calibration result is qualified when the correlation is greater than or equal to a preset first threshold value, thereby acquiring the time-depth conversion relation between the logging curve data and the seismic record and characteristic parameters sensitive to a reservoir; the first threshold is preferably selected to be 85%;

the characteristic parameter sensitive to the reservoir is obtained by the method comprising the following steps:

drawing a histogram by responding to logging parameters generated by different geologic bodies beside a well, and selecting the standardized logging curve data as characteristic parameters sensitive to a reservoir when the numerical value of the standardized logging curve data can distinguish data points above a second threshold value preset by different logging interpretation conclusions; the characteristic parameters sensitive to the reservoir at least comprise wave impedance IMP, well diameter CAL, natural gamma GR, natural potential SP, resistivity curve data: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and the data of the physical property characterization curve: one or more of compensated neutrons CNL, compensated acoustic curve AC, density curve DEN. The second threshold is preferably 70%.

When a carbonate rock ancient karst cave reservoir stratum is searched, a small karst cave reservoir body with the height of 0.5-5.0 m has the characteristics of wide distribution and large quantity and is used as a key object for reservoir prediction. The seismic reflection and inversion wave impedance characteristics only have certain response to reservoirs above 5m and cannot be identified deterministically; meanwhile, although conventional well logging is characterized by low resistivity, low density, increased neutron and acoustic wave time difference, etc., its detection range is limited. By the well seismic calibration technology, the well logging curve is used as hard data, the three-dimensional seismic waveform is used as soft data, the well logging constraint is established for the interpretation of the seismic data of the reservoir in a large range, and the identification precision of the reservoir can be greatly improved.

Firstly, calibrating and integrating well logging and seismic data, calculating a wave impedance curve according to the compensated acoustic curve AC and the density curve DEN, and then calculating a reflection coefficient curve. Constructing Rake wavelets based on the seismic dominant frequency of the target interval, synthesizing seismic records, comparing the artificially synthesized records with well-side seismic channels from the two aspects of the position of the marker layer and the seismic reflection event axis, and finally obtaining the time-depth conversion relation through the quality control of the correlation coefficients of the artificially synthesized records and the well-side seismic channels.

In this embodiment, a synthetic seismic record is obtained from the wave impedance information calculated from the compensated acoustic curve and the density value curve using a 25HZ (dominant frequency of the seismic data in the target interval).

Determining a eagle mountain group top boundary marker layer by observing seismic waveform data: compared with an overlying Bachu group mudstone layer, the waveform reflection of the inner curtain area of the carbonate rock is irregular and has no certain direction, the amplitude can be strong or weak, and the length of the in-phase axis can be long or short, and the continuity is poor; and has non-systematic in-phase axis reflection termination and bifurcation phenomena.

In this embodiment, the time-depth transformation relationship between the well log data and the seismic record is as follows:

wherein the content of the first and second substances,represents the seismic signal time corresponding to the initial depth of acoustic logging,the time difference of the sound wave is represented,representing the interval of sampling of the log data and T representing the seismic wave two-way travel time.

The single well seismic calibration chart is shown in FIG. 5;

and determining the paleo-karst cave position of the target interval according to the drilling, logging and core data, and dividing the reservoir types of sedimentary filling, collapse filling, chemical filling and mixed filling according to the lithology and combination relation of the fillers. And (4) optimizing logging curve parameters sensitive to reservoir type division by using a statistical mode of a logging curve histogram and an intersection graph. For example, the wave impedance curve reflects the difference of rock physical properties and can be used for distinguishing an ancient karst cave from surrounding rocks; the natural gamma curve utilizes the radioactive characteristics of the geologic body to effectively identify the fillers of the sedimentary sandstone and mudstone reservoir.

Step S600, obtaining a first group of PCA data: preferably selecting a preset number of pieces of standardized well logging curve data from the characteristic parameters sensitive to the reservoir, and analyzing the characteristic parameters sensitive to the reservoir by a PCA method to obtain a first group of PCA data sensitive to the reservoir;

in this embodiment, the PCA analysis specifically includes:

in this embodiment, 5 logging characteristic parameters, which are sensitive to reservoir classification, such as a wave impedance curve, a compensated neutron curve, a deep lateral resistivity curve, a shallow lateral resistivity curve, and a natural gamma curve are preferably used for PCA analysis.

Step B100, calculating a sample mean value:

wherein the content of the first and second substances,the total number of samples of the curve is shown,ia sample count parameter is represented that is representative of,is shown asiThe characteristic parameter data of each sample is obtained,represents the sample mean;

different logging curves reflect different physical parameters of the geologic body, and the influence of the dimensions on data analysis cannot be ignored due to different detection methods of various logging curves.

Step B200, carrying out normalization processing on each parameter curve to obtain normalized sample characteristic parameter data：

Wherein the content of the first and second substances,xrepresenting sample feature parameter sample point data;

step B300, based on the normalized sample characteristic parameter dataCalculating the covariance of the characteristic parameters of the sample：

Wherein the content of the first and second substances,represents the c-th X sample characteristic parameter data,representing the c sample characteristic parameter data;

step B400, based on the sample characteristic parameter covariance, a covariance matrix is constructed:

z represents sample characteristic parameter data;

Includedpthe covariance matrix of the seed sample characteristic parameter data is:

wherein the content of the first and second substances,representing a plurality of sample characteristic parameter data,representing the characteristic parameter data of the 1 st sample and thepCovariance of the sample feature parameter data;Kthere is no specific numerical meaning;

step B500, based on the inclusionpThe covariance matrix of the sample characteristic parameter data is planted, and the eigenvalue of the covariance matrix is obtainedAndpa feature vector，，；

Eigenvalue and eigenvector computations based on the covariance matrix includepPCA data for species eigenvalues:

wherein the content of the first and second substances,is shown asjA vector of the PCA data, and,is shown asjA vector of characteristic parameters of the individual samples,a feature vector representing the jth sample feature parameter data;j 1,2, pis a counting parameter;

step B600, based on includingpAnd (4) sequencing the corresponding eigenvalues from large to small according to the PCA data of the variety eigenvalues to obtain the eigenvalues sequenced from large to small，…，Calculating the variance cumulative contribution:

wherein the content of the first and second substances,represents the j-th characteristic value after sorting from big to small,to accumulate the first k covariance accumulated contributions,；

when k is increased from 1, the k is selectedAnd when the first time is greater than a preset third threshold value, corresponding PCA data vectors corresponding to the first k covariances are PCA data obtained by PCA analysis in the current round. In this embodiment, the third threshold is preferably 85%. The PCA data corresponding to the first k characteristic parameters reflect various parameters for characterizing the reservoirPart of the information, preferably the PCA data, can be used for further analyzing the lithology and physical properties of the reservoir by utilizing the PCA data.

In this embodiment, preferably, 5 well logging sample characteristic parameters are analyzed by the PCA method to obtain 5 PCA data. The variance contribution rate of the PC1 and the PC2 reaches more than 85%, the PC1 accounts for 48.72%, and the PC2 accounts for 36.96%, and can represent various logging response characteristics of the reservoir.

In this embodiment, the PC1 and PC2 obtained by the first round of PCA method are:

wherein PC1 represents the first data in the first set of PCA data for the first PCA analysis of the characteristic parametric curve for all data points; PC2 represents the second data in the first set of PCA data for the first PCA analysis of the characteristic parametric curve for all data points;natural gamma values for all data points;deep lateral resistivity values for all data points;shallow lateral resistivity values for all data points;a compensated neutron value for all data points;wave impedance values for all data points;

step S700, constructing an isochronous trellis model: and constructing an isochronous trellis model based on the sedimentary stratum rule reflected by the high-precision three-dimensional seismic amplitude data volume and the time-depth conversion relation. And after the position of the target layer section is determined, selecting an interface with continuous in-phase axis and stable deposition environment or a marker layer as a top-bottom interface of the range to be predicted according to the seismic data profile. And constructing an isochronous trellis model based on the actual geological structure background and the time-depth transformation relation.

Step S800, interwell reservoir parameter simulation: determining an optimal sample number parameter capable of reflecting the overall geological condition based on a first group of PCA data of each sample well, selecting the first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample number parameter to construct an initial model, continuously correcting parameters of the initial model, and outputting a first group of high-precision characteristic value simulation result data volumes which are data volumes corresponding to the first group of PCA data one by one;

in this embodiment, the step S800 specifically includes steps S810 to S890:

step S810, selecting a sample well as a reference target well, and setting an initial sample number parameter to be 1;

step S820, selecting a first group of PCA data of the sample wells with the quantity being the sample number parameter and a first group of PCA data of the reference target well mark according to the waveform similarity principle to perform correlation analysis to obtain a first group of PCA data correlation values of the sample number parameter-reference target well;

step S830, increasing the sample number parameters 1 by 1, repeating the method of the step S820 to obtain a first group of PCA data correlation values of the sample number parameter-reference target well corresponding to each sample number parameter, connecting the first group of PCA data correlation values of all the sample number parameters-reference target wells, and obtaining a correlation curve of the first group of PCA data correlation of the reference wells along with the change of the sample number parameters;

step S840, selecting another sample well as a reference target well, repeating the steps S810-S830, obtaining a correlation curve of the first group of PCA data correlations of the plurality of reference wells along with the change of the sample number parameters, fitting the correlation curve of the first group of PCA data waveform correlations of all the reference wells along with the change of the sample number parameters into an overall correlation curve, selecting an inflection point at which the correlation in the overall correlation curve rises along with the increase of the sample number parameters and finally keeps stable, and determining the optimal sample number parameters; the effect of fitting a correlation curve of the first set of PCA data waveform correlations versus sample number parameters for all reference wells to an overall correlation curve is shown in fig. 6;

in the embodiment, the well with similar low-frequency structure is preferably selected as the spatial estimation sample by using the waveform similarity and spatial distance bivariate in the sample well, so that the low-frequency change of the spatial structure can be well reflected.

Two wells with similar seismic waveforms indicate that the deposition environments are similar, the low-frequency components of the wells have commonality, the certainty of the low-frequency section of the inversion result can be enhanced, the value range of high frequency is restricted, and the reliability of the inversion and simulation result is improved.

Step S850, based on the high-precision three-dimensional seismic amplitude data volume and the isochronous lattice model, calculating the waveform correlation between the point to be detected and the sample well positions, sorting the waveform correlation from large to small, and selecting a first group of PCA data of sample wells with the optimal sample quantity and the highest seismic waveform correlation; constructing an initial model based on the sample well corresponding to the seismic waveform characteristic data of the sample well with the highest correlation through an inter-well characteristic parameter interpolation mode;

step S860, based on the initial model, selecting a first group of PCA data of the sample well with the optimal sample number parameter bar and the highest seismic waveform correlation degree as prior information;

step S870, performing matched filtering on the initial model and the prior information to obtain a maximum likelihood function;

step S880, based on the maximum likelihood function and the prior information, obtaining the posterior probability statistical distribution density under a Bayes framework, and sampling the posterior probability statistical distribution density to obtain a target function;

in this embodiment, the step S880 specifically includes:

step S881, using white noise to satisfy the rule of gaussian distribution, representing the parameters of the high-precision eigenvalue simulation result data volume as:

y represents parameters of a logging curve high-precision characteristic value simulation result data volume, X represents actual characteristic parameter values of an underground stratum to be solved, and N represents random noise;

step S882, becauseAlso satisfies the GaussDistribution, the initial objective function can be determined as:

wherein the content of the first and second substances,a function relating to a posteriori information is represented,showing that the characteristic curve of the sample well is matched and filtered by selecting the sample well based on the optimal sample number, the posterior probability statistical distribution density is obtained, and the expected value of the characteristic parameter is further calculated,a covariance representing white noise;

step S883, based on the initial objective function, introducing prior information into the objective function through maximum posterior estimation, and obtaining a stable objective function as:

wherein the content of the first and second substances,representing the characteristic parameter to be simulated,representing functions related to prior information such as geological and well log data,representation for coordinationAndthe smoothing parameters of the mutual influence between them.

Step 890, taking the target function as the input of the initial model, sampling posterior probability distribution by a Markov chain Monte Carlo method MCMC and Metropolis-Hastings sampling criterion, continuously optimizing parameters of the initial model, selecting a solution of the target function when the maximum value is taken as random realization, taking the average value of multiple random realization as expected value output, and taking the expected value output as a high-precision characteristic value simulation result data body; and parameters in the high-precision characteristic value simulation result data volume correspond to characteristic parameters corresponding to the first group of PCA data one to one.

Step S900, depicting the boundary of the karst cave system between wells: performing intersection analysis based on a first group of PCA data of the sample well to obtain a first group of geologic body discrete data point distribution maps;

obtaining an interpretation conclusion of a first group of geologic body-like discrete data point distribution diagrams based on the lithological information and the physical property information, dividing a karst cave and surrounding rocks, further constructing a first group of PCA data intersection layouts, and obtaining a first group of PCA data threshold values required by the division of the karst cave;

in this embodiment, the PCA data PC1 and PC2 are subjected to intersection analysis to obtain the discrete data point distribution of each type of geologic body. And obtaining a preferred PCA data intersection graph interpretation conclusion based on the determined lithology and physical property information of the individual depth section obtained by the imaging logging, the drilling logging and the core data, distinguishing the karst cave from the surrounding rock, constructing a preferred PCA data intersection graph plate, and obtaining the PC1 and PC2 data threshold values required by the karst cave division. A schematic diagram of the first set of PCA data intersection analysis capturing karst cavern data points and gating threshold values is shown in fig. 7;

and performing intersection analysis on the first group of high-precision characteristic value simulation result data volumes based on the first group of PCA data threshold values to obtain three-dimensional space morphological characteristics of the karst cave. The structural characteristics and the spatial distribution rule of the ancient karst cave type reservoir stratum with the size of more than 2m can be obtained. The cross-sectional effect of three-dimensional morphology and structure is shown in fig. 8.

Step S1000, depicting a lithologic property boundary filled in the cave: analyzing the characteristic parameters sensitive to the reservoir, which correspond to the karst cave data points, by a PCA method to obtain a second group of PCA data sensitive to the response of the filler;

generating a second group of high-precision characteristic value simulation result data volume by the method of S700 based on the second group of PCA data;

performing intersection analysis based on the second group of PCA data sensitive to the filler response to obtain a second group of geologic body discrete data point distribution diagrams;

obtaining an interpretation conclusion of a second group of geologic body-like discrete data point distribution diagrams based on the lithological information and the physical property information, distinguishing filler types, further constructing a second group of PCA data intersection layout, and obtaining a second group of PCA data threshold values required by the division of the filler types in the karst cave;

performing intersection analysis on the second group of high-precision characteristic value simulation result data volumes based on the second group of PCA data threshold values to obtain three-dimensional space morphological characteristics of different filler types in the karst cave; a schematic diagram of the intersection analysis of the second set of PCA data to partition the interior filling of the karst cavity and to select the threshold values is shown in fig. 9;

in this embodiment, the karst cave interpretation data points obtained by intersection of the PCA data PC1 and the PC2 are extracted, and the corresponding five characteristic parameter curves are subjected to PCA analysis again to obtain a second set of PCA data F1 and F2 sensitive to the response of the filling. And (4) carrying out intersection analysis on the PCA data F1 and F2 to obtain the discrete data point distribution of various geologic bodies. And obtaining a second group of PCA data intersection graph interpretation conclusion based on the determined lithology and physical property information of the individual depth section obtained by the imaging logging, the drilling logging and the core data, distinguishing the sedimentary filling and the collapse filling, constructing a second group of PCA data intersection graph, and obtaining F1 and F2 data threshold values required by the division of the sedimentary filling and the collapse filling in the karst cave.

In this example, the F1 and F2 data obtained by the second round of PCA method are:

wherein F1 represents a first PCA data of the second set of PCA data for the PCA analysis of the cavern characterizing data point characteristic parameter curve; f2 is second PCA data in the second group of PCA data for carrying out PCA analysis on the characteristic parameter curve of the cavern characterization data point for the second time;characterizing a natural gamma value of a data point for the cavern;characterizing the deep lateral resistivity value of the data point for the cavern;shallow lateral resistivity values characterizing data points for caverns;a compensated neutron value for a cavern characterization data point;wave impedance values for cavern characterization data points

Performing intersection analysis on the second group of PCA data simulation result data bodies through an interpretation conclusion threshold value for dividing the sedimentary filling and the collapsed filling in the karst cave to obtain three-dimensional space morphological characteristics of the sedimentary filling and the collapsed filling in the karst cave;

step S1100, archaeological cavern structure and filling description: and engraving the development characteristics of the ancient karst cave space distribution and the internal filling of different filler types by adopting a lithologic shielding technology and a three-dimensional engraving technology based on the three-dimensional space morphological characteristics of the karst cave and the three-dimensional space morphological characteristics of different filler types in the karst cave. The section of the development characteristic of the filler in the ancient karst cave is shown in figure 10;

the system for evaluating the filling characteristics of the deep paleo-karst reservoir by well-seismic combination according to the second embodiment of the invention comprises: 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 seismic calibration and characteristic parameter selection module, a first group of PCA data acquisition modules, an isochronous grid model construction module, an interwell reservoir parameter simulation module, an interwell karst cave system boundary description module, a lithologic property boundary description module for describing filling inside a cave and an ancient karst cave structure and filling description module;

the original geophysical logging data acquisition module is configured to acquire original logging data of each sample well through logging equipment, and comprises: measuring the natural potential SP of each sample well through a measuring electrode, measuring the natural gamma GR of each sample well through a natural gamma underground device and a natural gamma ground instrument, and obtaining the well diameter CAL of each sample well through a well diameter arm; obtaining resistivity curve data by conventional logging equipment: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and physical property characterization curve data: a compensated neutron CNL, a compensated acoustic curve AC and a density curve DEN;

obtaining determined lithology information and physical property information of individual depth sections based on imaging logging information, drilling information, logging information and core information, and further determining depth data of a target horizon marker layer;

the original geophysical logging data preprocessing module is configured to draw logging curve data based on original logging data of the sample well, perform abnormal value processing and standardization processing, and obtain standardized logging curve data;

the seismic data preprocessing module is configured to obtain a high-precision three-dimensional seismic amplitude data volume through mixed phase wavelet deconvolution and diffusion filtering based on the seismic wave reflection signal data;

in this embodiment, based on the correlation between the borehole coordinates and the coordinates of the seismic signal detectors, the one-to-one correspondence between the borehole and the borehole-side seismic waveform is determined, and then the plane positions of the 9 logging curves and the borehole-side seismic waveforms corresponding to the plane positions are determined.

The well-seismic calibration and characteristic parameter selection module is configured to perform well-seismic calibration and characteristic parameter selection: acquiring a wave impedance curve of a sample well based on a compensation acoustic curve and a density curve DEN in the standardized logging curve data, further calculating a reflection coefficient curve, acquiring the preferred frequency of a Rake wavelet to keep the same as the main frequency of a high-precision three-dimensional seismic amplitude data body, performing convolution operation on the Rake wavelet and the reflection coefficient curve to obtain a synthetic seismic record, comparing the depth data of a target layer position mark layer with the isochronous three-dimensional spread of the target layer position mark layer for well seismic calibration, calculating the correlation between the synthetic seismic record and a well-side seismic channel waveform, and judging that a well seismic calibration result is qualified when the correlation is greater than or equal to a preset first threshold value, thereby acquiring the time-depth conversion relation between the logging curve data and the seismic record and characteristic parameters sensitive to a reservoir;

the characteristic parameter sensitive to the reservoir is obtained by the method comprising the following steps:

drawing a histogram by responding to logging parameters generated by different geologic bodies beside a well, and selecting the standardized logging curve data as characteristic parameters sensitive to a reservoir when the numerical value of the standardized logging curve data can distinguish data points above a second threshold value preset by different logging interpretation conclusions; the characteristic parameters sensitive to the reservoir at least comprise wave impedance IMP, well diameter CAL, natural gamma GR, natural potential SP, resistivity curve data: deep lateral logging RLLD, shallow lateral logging RLLS and micro lateral logging RLLM, and the data of the physical property characterization curve: one or more of a compensated neutron CNL, a compensated acoustic curve AC, a density curve DEN;

the first group of PCA data acquisition modules are configured to optimize a preset number of pieces of standardized well logging curve data from the characteristic parameters sensitive to the reservoir, and perform dimensionality reduction on the characteristic parameters sensitive to the reservoir through PCA analysis to obtain a first group of PCA data;

the isochronous trellis model building module is configured to build an isochronous trellis model based on sedimentary stratum rules reflected by the high-precision three-dimensional seismic amplitude data volume and the time-depth conversion relation;

the interwell reservoir parameter simulation module is configured to determine an optimal sample number parameter capable of reflecting the overall geological condition based on a first group of PCA data of each sample well, select a first group of PCA data of the sample well with the highest seismic waveform correlation of the optimal sample number parameter to construct an initial model, continuously correct parameters of the initial model, and output 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 are data bodies corresponding to the first group of PCA data one by one;

the boundary delineation module of the inter-well karst cave system is configured to perform intersection analysis based on a first group of PCA data of the sample well to obtain a first group of geologic body-like discrete data point distribution maps;

obtaining an interpretation conclusion of a first group of geologic body discrete data point distribution diagrams based on the lithological information and the physical property information, dividing a karst cave and surrounding rocks, further constructing a first group of PCA data intersection layouts, and obtaining a first group of PCA data threshold values required by the division of the karst cave;

performing intersection analysis on the first group of high-precision characteristic value simulation result data volumes based on the first group of PCA data threshold values to obtain three-dimensional space morphological characteristics of the karst cave;

the module for depicting the rock property boundary filled in the cave is configured to depict the rock property boundary filled in the cave: analyzing the characteristic parameters sensitive to the reservoir, which correspond to the karst cave data points, by a PCA method to obtain a second group of PCA data sensitive to the response of the filler;

generating a second set of high-precision characteristic value simulation result data volume through the function as the interwell reservoir parameter simulation module based on the second set of PCA data;

performing intersection analysis based on the second set of PCA data sensitive to the filler response to obtain a second set of geologic body-like discrete data point distribution diagram;

obtaining an interpretation conclusion of a second group of geologic body-like discrete data point distribution diagrams based on the lithological information and the physical property information, distinguishing filler types, further constructing a second group of PCA data intersection layout, and obtaining a second group of PCA data threshold values required by the division of the filler types in the karst cave;

performing intersection analysis on the second group of high-precision characteristic value simulation result data volumes based on the second group of PCA data threshold values to obtain three-dimensional space morphological characteristics of different filler types in the karst cave;

the paleo-karst cavern structure and fill description module is configured to: and engraving the development characteristics of the ancient karst cave space distribution and the internal filling of different filler types by adopting a lithologic shielding technology and a three-dimensional engraving technology based on the three-dimensional space morphological characteristics of the karst cave and the three-dimensional space morphological characteristics of different filler types in the karst cave.

It can be clearly understood by those skilled in the art that, for convenience and brevity of description, the specific working process and related description of the system described above may refer to the corresponding process in the foregoing method embodiments, and will not be described herein again.

It should be noted that, the system for intelligently dividing filling characteristics of an ancient karst cave reservoir provided in the foregoing embodiment is only illustrated by the division of the above functional modules, and in practical applications, the functions may be allocated to different functional modules according to needs, that is, the modules or steps in the embodiment of the present invention are further decomposed or combined, for example, the modules in the foregoing embodiment may be combined into one module, or may be further split into multiple sub-modules, so as to complete all or part of the functions described above. The names of the modules and steps involved in the embodiments of the present invention are only for distinguishing the modules or steps, and are not to be construed as unduly limiting the present invention.

An electronic apparatus according to a third embodiment of the present invention includes: at least one processor; and a memory communicatively coupled to at least one of the processors; wherein the memory stores instructions executable by the processor for performing the method for joint well-seismic evaluation of deep paleo-reservoir packing characteristics as described above.

A computer-readable storage medium according to a fourth embodiment of the present invention stores computer instructions for execution by the computer to implement the method for joint well-seismic evaluation of deep paleo-karst reservoir packing characteristics as described above.

It can be clearly understood by those skilled in the art that, for convenience and brevity of description, the specific working processes and related descriptions of the storage device and the processing device described above may refer to the corresponding processes in the foregoing method embodiments, and are not described herein again.

Computer program code for carrying out operations for aspects of the present application may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C + + or the like and conventional procedural programming languages, such as the "C" programming language or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the case of a remote computer, the remote computer may be connected to the user's computer through any type of network, including a Local Area Network (LAN) or a Wide Area Network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet service provider).

The flowchart and block diagrams in the figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present application. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems which perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.

The terms "first," "second," and the like are used for distinguishing between similar elements and not necessarily for describing or implying a particular order or sequence.

The terms "comprises," "comprising," or any other similar term are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus.

So far, the technical solutions of the present invention have been described in connection with the preferred embodiments shown in the drawings, but it is easily understood by those skilled in the art that the scope of the present invention is obviously not limited to these specific embodiments. Equivalent changes or substitutions of related technical features can be made by those skilled in the art without departing from the principle of the invention, and the technical scheme after the changes or substitutions can fall into the protection scope of the invention.