CN112213775A - Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data - Google Patents

Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data Download PDF

Info

Publication number
CN112213775A
CN112213775A CN202010974716.7A CN202010974716A CN112213775A CN 112213775 A CN112213775 A CN 112213775A CN 202010974716 A CN202010974716 A CN 202010974716A CN 112213775 A CN112213775 A CN 112213775A
Authority
CN
China
Prior art keywords
frequency
domain
stack
gather
data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010974716.7A
Other languages
Chinese (zh)
Other versions
CN112213775B (en
Inventor
包燚
王成
吴佳忆
陈树民
裴江云
陈志德
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Petrochina Co Ltd
Daqing Oilfield Co Ltd
Original Assignee
Petrochina Co Ltd
Daqing Oilfield Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Petrochina Co Ltd, Daqing Oilfield Co Ltd filed Critical Petrochina Co Ltd
Priority to CN202010974716.7A priority Critical patent/CN112213775B/en
Publication of CN112213775A publication Critical patent/CN112213775A/en
Application granted granted Critical
Publication of CN112213775B publication Critical patent/CN112213775B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/20—Trace signal pre-filtering to select, remove or transform specific events or signal components, i.e. trace-in/trace-out
    • G01V2210/24—Multi-trace filtering
    • G01V2210/244—Radon transform
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G01V2210/322—Trace stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G01V2210/324—Filtering

Abstract

The invention relates to the technical field of seismic exploration, in particular to a fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data. The method comprises the following steps: carrying out applicability judgment on pre-stack gather data of seismic data to be processed, and if the pre-stack gather data meets the conditions, carrying out residual time difference elimination on a pre-stack CRP gather; carrying out random noise suppression treatment based on prediction errors on the pre-stack CRP gathers with the residual time difference removed; carrying out stability inverse Q filtering processing based on time-frequency analysis; carrying out weighted least square radon transform denoising treatment; performing time-varying pulse deconvolution processing; and denoising and coherent superposition processing. The fidelity frequency-extracting method for the pre-stack seismic data with high coverage times can effectively improve the signal-to-noise ratio of the high-frequency components of the pre-stack CRP gather. The resolution is greatly improved, the transverse resolution is not lost, the comprehensive quality of the CRP gather is improved, and the pre-stack inversion requirement can be met.

Description

Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data
Technical Field
The invention relates to the technical field of seismic exploration, in particular to a fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data.
Background
The higher the resolution of the earthquake processing result is, the stronger the thin layer identification capability and the micro geologic body depicting capability are, but the current common method for improving the resolution after stacking has a serious statistical effect, can reduce the identification capability of the transverse geologic body and destroy the original vertical amplitude relationship and wave resistance relationship. The resolution enhancement processing performed on the pre-stack seismic data can effectively avoid the problems, but the pre-stack CRP gather signal-to-noise ratio is generally low, the signal-to-noise ratio is more seriously reduced when the frequency is enhanced by only adopting one or two frequency enhancement methods, and the single frequency enhancement technology can greatly amplify the high-frequency components of the pre-stack seismic data, so that the resolution of the pre-stack seismic data is enhanced by adopting a plurality of technologies, and the signal-to-noise ratio is still within the acceptable range after stacking.
At present, a plurality of prestack data frequency boosting technologies exist, but the following problems generally exist: the signal to noise ratio is reduced when the resolution is improved, the AVO characteristics of the trace set are changed after amplitude compensation and equalization, and the visual resolution of the superposed profile is reduced after the frequency band is widened.
The three problems are that any pre-stack gather frequency extracting technology cannot avoid, and the reason for this is that the initial signal-to-noise ratio of the high-frequency component in the pre-stack CRP gather is very low, and the random noise mixed into the high-frequency component is greatly amplified by adopting any frequency extracting means, so that the signal-to-noise ratio is greatly reduced. The pre-stack CRP gather has an important role in predicting reservoir lithology and fluid by seismic pre-stack inversion except that a post-stack section is generated by stacking, the theoretical basis of the pre-stack CRP gather is a Zoeppritz equation set, and Qstrander finds the phenomena that the amplitude of gas-containing sandstone is increased along with the increase of offset and the amplitude of water-containing sandstone is reduced along with the increase of offset in 1984, marks the start of AVO analysis, and more analysis means are moved from the post-stack to the pre-stack. Shuey proposed a simplified formula of AVO analysis in 1985, and the reflection coefficients are expressed by far, middle and near incidence angles, and the relation between the amplitude and the incidence angle is visually reflected. Later Connolly proposed Elastic Impedance (EI) inversion and mathematical expressions in 1999, further enriching the connotation of AVO analysis. However, the application bases of the above prediction techniques are that the AVO characteristics of the prestack gather are to be maintained, the AVO characteristics of the prestack gather after amplitude compensation equalization are usually changed, and the original AVO characteristics of the prestack CRP gather can be maintained only by avoiding the use of the transverse amplitude energy equalization technique in the prestack gather frequency boosting technique. The problem that the visual resolution of the superposed section is reduced after the frequency band is widened is caused by the fact that the dominant frequency band energy of the original superposed section is prominent, the signal-to-noise ratio of the low-frequency component after the frequency band is widened is generally low, and the visual resolution of the superposed section is influenced, so that the frequency band is wide, and the capability of visually identifying the thin layer is reduced.
Disclosure of Invention
Technical problem to be solved
The invention provides a fidelity frequency-boosting method for prestack seismic data with high coverage times, which aims to overcome the defects that the signal-to-noise ratio of a prestack CRP gather is low, the AVO characteristic of the gather is changed after amplitude compensation and equalization, the visual resolution of a stacking section is reduced after a frequency band is widened and the like in the prior art.
(II) technical scheme
In order to solve the problems, the invention provides a fidelity frequency-extracting method of high-coverage-number pre-stack seismic data, which comprises the following steps:
s1, judging the applicability of the seismic data to be processed pre-stack gather data, if the applicability meets the condition, performing S2, and if the applicability does not meet the condition, terminating the operation;
step S2, residual time difference elimination is carried out on the CRP gather before folding;
step S3, carrying out random noise suppression treatment on the prestack CRP gather which is processed in the step S2 and eliminates the residual time difference based on the prediction error;
step S4, performing stability inverse Q filtering processing based on time-frequency analysis on the pre-stack CRP gather processed in the step S3;
s5, carrying out weighted least square Radon transformation denoising processing on the pre-stack CRP gather processed in the S4;
step S6, performing time-varying pulse deconvolution processing on the pre-stack CRP gathers processed in the step S5;
and S7, performing denoising and coherent superposition processing on the pre-stack CRP gather processed in the step S6.
Preferably, the step S1 of determining the applicability of the seismic data prestack gather data to be processed specifically includes:
when the number of times of coverage of the CMP gather of the seismic data to be processed is more than or equal to 200 times, or the number of times of coverage of the CRP gather is more than or equal to 60 times, executing the step S2;
the operation is terminated when the number of CMP gather overlays is less than 200 or the number of CRP gather overlays is less than 60.
Preferably, in step S2, the residual moveout elimination for the pre-stack CRP gather specifically includes:
when the residual time difference of the target layer in-phase axis of the pre-stack CRP gather is greater than or equal to 20ms when the reference offset distance is 3000m, the processing is terminated, and the velocity analysis is carried out again;
when the residual time difference is less than 20ms and more than or equal to 10ms, adjusting the speed by utilizing the speed analysis based on the spatial continuity, and eliminating the residual time difference;
when the residual time difference is less than 10ms, eliminating the residual time difference by using a spatial weighted median filtering method, wherein the output amplitude of each sampling point of the CRP gather in-phase axis after the spatial weighted median filtering is as follows:
Figure BDA0002685364760000031
wherein, VKIs the non-abnormal amplitude value in the aperture, N is the number of non-abnormal amplitude values in the aperture, WeightkIs the space weight;
after the treatment, the residual time difference of the gather is eliminated.
Preferably, the step S4 is performedThe stability inverse Q filtering process of time-frequency analysis specifically includes:
Figure BDA0002685364760000032
Figure BDA0002685364760000033
the formula is an inverse Q filtering equation based on Gabor transformation, tau is depth, j is an imaginary unit, omega is angular frequency, Q is stratum quality factor, for wave fields with different time depths, Gabor transformation is carried out on surface seismic records, then inverse Q filtering is carried out in a time-frequency domain by combining an amplitude compensation operator and a phase compensation operator, and the seismic records P (tau) processed by inverse Q filtering can be obtained by utilizing inverse Gabor transformation on the obtained wave fields P (tau, omega).
Preferably, the weighted least squares radon transform denoising processing performed in step S5 specifically includes:
step S51, least squares radon transform, comprising: the least square regularization method is adopted in a frequency domain, the least square discrete parabola Ladong positive transformation is established, and the mathematical expression is as follows:
M＝(LLH+λ2I)-1LD
in the formula: d is f-x domain data, M is f-q domain data, (LL)H)-1Is a deconvolution operator; wherein
Figure BDA0002685364760000041
Figure BDA0002685364760000042
Converting the pre-stack seismic data of the time-space domain into a tau-p domain by using the method;
step S52, tau-p domain weighted coherent noise suppression, comprising: introducing a nonlinear regularization diagonal matrix W into an operator of least square parabola Radon transform to realize, and attenuating tau-p domain coherent noise by time difference of multiples and primary waves: (LL)H+WHW)M＝LD
In the formula, the physical meaning of a nonlinear diagonal matrix W is a weight vector of a curvature parameter q, effective reflection is predicted in a tau-p domain, the energy of multiple coherent noise in a tau-p domain horizontal ray parameter nonzero region is inverted through a least square method, the energy of the nonzero region is suppressed in a data-driven automatic weighting mode, and the coherent noise is suppressed in the tau-p domain;
step S53, Stau-p domain residual noise suppression: in step S52, the weighted least-squares radon transform noise suppression is performed again on the full-frequency data other than the dominant frequency band by the same method.
And step S54, reconstructing the tau-p domain data processed in the steps back to a time-space domain through the control of a complex conjugate transpose matrix, and finishing the whole processing.
Preferably, the performing of the variable pulse deconvolution process in step S6 specifically includes:
calculating a deconvolution operator a (t) by using the basic equation of the pulse deconvolution,
calculating a deconvolution operator a (t) by using a fundamental equation of deconvolution,
Figure BDA0002685364760000051
in the formula, rxx(m) is the seismic reflection coefficient, a (t) is the pulse deconvolution operator; utilizing convolution of a '(t) and seismic record x (t) to obtain S (t) ═ a' (t) × (t), wherein S (t) is seismic trace output by pulse deconvolution, performing windowing calculation on input seismic record in the processing process to obtain pulse deconvolution operator a (t), wherein each 2 analysis calculation time windows have an overlapping transition region with a certain length, and the length of the window overlapping region is more than or equal to 4 seismic wavelet lengths; performing pulse deconvolution in a frequency domain by participating all frequency components of an input seismic record in time-varying pulse deconvolution calculation; the deconvolution prediction step size is equal to the sampling interval of the input data.
(III) advantageous effects
The fidelity frequency-extracting method for the pre-stack seismic data with high coverage times can effectively improve the signal-to-noise ratio of the high-frequency components of the pre-stack CRP gather. According to the seismic data processed by the method, after the superposition, the middle-shallow layer broadened frequency band is coherently superposed, the middle-shallow layer broadened frequency band is more than 20Hz, the deep layer broadened frequency band is more than 10Hz, the resolution is greatly improved, the transverse resolution is not lost, the comprehensive quality of the CRP gather is improved, and the pre-stack inversion requirement can be met.
Drawings
FIG. 1 is a flowchart of a fidelity frequency-boosting method for high-coverage pre-stack seismic data according to an embodiment of the invention.
Detailed Description
The present invention will be described in detail below with reference to the drawings and examples.
As shown in fig. 1, an embodiment of the present invention provides a method for fidelity frequency extraction of high coverage number pre-stack seismic data, including:
step S1, method applicability analysis: the method has the advantages that the maximum covering times of the CMP gather of the target layer and the maximum covering times of the CRP gather are determined by analyzing the type of the acquisition observation system applying the target data and combining an offset distance grouping scheme during offset processing, the method is high in applicability and low in requirement on data, and the data effect is better when the covering times of the CMP gather are more than or equal to 200 times and the covering times of the CRP gather are more than or equal to 60 times (the high covering times can ensure the signal-to-noise ratio of the post-stack data).
Step S2, residual time difference elimination is carried out on the CRP gather: residual moveout can be eliminated by adopting a method based on amplitude energy or a method based on residual velocity, and the CRP gather in-phase axis is leveled out. The processing methods of the energy-based gather leveling are numerous, and the core requirement is to eliminate residual time difference as far as possible, so that the application of subsequent processing technology is facilitated to obtain expected internal effects.
Classifying residual time difference existing in the CRP gather before stacking, for example, at a displacement distance of 3000m, if amplitude time of a same-phase axis mark of the point and self-excitation and self-collection time of the same-phase axis are greater than 20ms, judging that the CRP gather before stacking has large residual time difference, performing velocity analysis again, and performing migration processing again after a migration velocity is formed again; if the residual time difference is more than 10ms and less than 20ms, the pre-stack CRP gather is regarded as having moderate residual time difference, and the speed can be adjusted by utilizing the technologies such as Space Continuity Velocity Analysis (SCVA) and the like to eliminate the residual time difference; if the residual time difference is less than 10ms, the pre-stack CRP gather is regarded as having a tiny residual time difference, and the residual time difference can be eliminated by methods such as prior parameter median filtering or weighted median filtering.
Before median filtering processing, the number x of spatial median filtering channels of the pre-stack CRP channel set and an abnormal value rejection coefficient y are given, wherein x is an integer of the number of spatial channels, and y is the proportion of abnormal amplitude values at two ends of a sequence, and is a percentage.
Median filtered output amplitude aoutput＝(Amp1+Amp2+…+Ampn)/n (1)
Amp in formula (1)1，Amp2，…，AmpnThe amplitude values left after the abnormal values (including the abnormal large value and the abnormal small value) of the total number x of data are removed through the prior parameter y are not the full amplitude values of the x of data.
Take a set of input data with x equal to 5 as an example, a position is a3The seismic channel amplitude of (1) is 78, 2 data of each front and back channel are selected as input data, and the amplitude of 5 data is a1＝90，a2＝110，a3＝78，a4＝258，a5121, the prior outlier rejection parameter y is 40%, in the first step, 5 data are sorted according to size, and the result is b1＝78，b2＝90，b3＝110，b4＝121，b5258, the prior parameter of abnormal amplitude rejection is 40%, namely 20% abnormal large value and 20% abnormal small value are removed, and the result after removal is c1＝90，c2＝110，c3121, the amplitude of the channel after median filtering is obtained by calculation of the formula (1)output＝(90+110+121)/3＝107。
The residual time difference can be better suppressed by the weighted median filtering for the three-dimensional data, the pre-stack CRP gather can directly push the CRP gather two-dimensional prediction method to a three-dimensional space to not only count the continuity of the same-phase axis of one CRP gather, but also judge whether the corresponding same-phase axes of several or more than ten CRP gathers on the same-phase axis in the space still have the same spatial continuity, and the residual time difference is suppressed by the spatial 3D weighted median filtering.
The number of dominant line weighted median filter bins is:
Figure BDA0002685364760000071
the number of cross-line weighted median filter bins is:
Figure BDA0002685364760000072
in the formulas (2) and (3), INT represents the rounding of the calculation result in parentheses, x and y are the calculation ranges of the main measuring line and the contact measuring line respectively, the unit is m, and Δ x and Δ y are the bins in the directions of the main measuring line and the contact measuring line, and the unit is m.
The major survey line direction ellipse radius is:
Rinline＝[INT(No._cells_inline/2)+0.5]Δx (4)
the contact measurement line ellipse radius is:
Rcrossline＝[INT(No._cells_crossline/2)+0.5]Δy (5)
the weight is:
Figure BDA0002685364760000073
in the formula (6), x and y are absolute distances from any point to the calculation point in the main measuring line/cross measuring line direction, and the unit is m and RinlineAnd RcrosslineThe unit of the ellipse radius of the point in the direction of the main measuring line/cross measuring line is m, Decay is an attenuation coefficient, and the weight value of the weighted median filter is determined by the size of the ellipse radius.
The output amplitude of each weighted median filter calculation point is:
Figure BDA0002685364760000081
formula (7) VKIs the non-abnormal amplitude in the apertureThe value, N, is the number of non-anomalous amplitude values in the aperture, and Weight is the calculated Weight.
After the treatment, the residual time difference of the prestack CRP gather is basically eliminated.
Step S3, CRP gather random noise suppression: on the basis of eliminating residual time difference, a two-dimensional or three-dimensional space prediction error filtering technology can be adopted, the effect of a three-dimensional random noise suppression technology is better, if the step of random noise suppression is lacked, the frequency band is greatly widened by the subsequent processing technology, and the high-frequency random noise energy is amplified in an exponential order, so that the whole prestack fidelity frequency boosting method fails.
The random noise suppression technology can effectively suppress interference such as random noise and the like, improves the signal-to-noise ratio of the pre-stack CRP gather, can perform random noise suppression on the two-dimensional CRP gather according to three-dimensional data through way head conversion in actual calculation, and firstly performs Fourier forward transform on a data volume to obtain:
Figure BDA0002685364760000082
in the formula (8), D (x, y, ω) is a fourier forward transform of the three-dimensional time domain seismic data D (x, y, t), and ω is a frequency;
any frequency spectrum component omega in f-x domainiData D (x, y, ω) given by equation (8)i) Only a space position function is recorded as f (x, y), a rectangular frame (namely a two-dimensional predictor) is designed, and n groups of different-speed homophase axes R exist in the windowi(t) the time difference of each in-phase axis in the x and y directions is Deltatx(i) And Δ ty(i) Where i is 1,2, …, n, then the two-dimensional seismic data f (x, y) for a particular frequency component may be represented in the form of a two-dimensional Z-transform:
Figure BDA0002685364760000091
f (Z) in the formula (9)1,Z2) As a result of Z transformation of the signal f (x, y), r (f) is a spectral domain representation of the in-phase axis r (t), and the right end of equation (9) is expanded to obtain:
E(Z1,Z2)＝F(Z1,Z2)P(Z1,Z2) (10)
in the formula (10), E (Z)1,Z2) For the result of the Z transformation of the prediction error, P (Z)1,Z2) For the prediction error filter, F (Z)1,Z2) For the input signal, equation (10) divides the signal F (Z)1,Z2) Applied as a prediction error filter, and filtered to obtain E (Z)1,Z2) And (6) obtaining the result.
Prediction error filter P (Z)1,Z2) Actually, the prediction error output by the three-dimensional seismic data after being predicted by the rectangular predictor is the minimum, an objective function is constructed to enable the partial derivative of the objective function to the predictor to be 0, and a matrix equation is obtained: rH·P＝R (11)
In the formula (11), RHThe method is characterized in that the method is an f-x domain seismic data multi-channel autocorrelation Hermite matrix, R is an f-x domain seismic data multi-channel autocorrelation array, and P is each component array of a prediction operator. And (3) performing inverse Fourier transform on the data processed by the formula (11) to obtain a pre-stack CRP gather after t-x domain random noise suppression.
Step S4, gather bluing: the step adopts a processing means mainly based on inverse Q filtering, and aims to expand the dominant frequency band to a high-frequency end. Note that this step does not necessarily use deconvolution-like frequency broadening means, and does not broaden the frequency band towards both high and low frequencies. In the bluing process, the degree of broadening to the high band is properly controlled with reference to the number of covering times, and taking the time-varying inverse Q filtering process technique as an example, if the data spectrum frequency-down relationship is natural and there is a higher number of covering times, the Q value range can be properly reduced at the destination layer. The trace gather data processed by the step has the advantages that the effective reflection energy of the high-frequency end is improved, the difference between the level of the high-frequency reflection energy and the random noise (or background noise) of the high-frequency end is reduced, and the residual random noise in the CRP trace gather is amplified again after the step 3, so that the noise is continuously suppressed in the next step, and the signal-to-noise ratio of the high-frequency effective reflection is improved.
When a plane wave propagates in a viscoelastic medium, the single-pass wave analysis solution is as follows:
P(z+Δz,ω)＝P(z,ω)exp[-jk(ω)Δz] (12)
in formula (12): z is depth, j is an imaginary unit, and ω is angular frequency, where the effect of formation Q can be expressed by introducing a complex wave k (ω) as:
Figure BDA0002685364760000101
in formula (13): qrAnd upsilonrQuality factor and phase velocity, ω, respectively, of the reference frequencyhFor the tuning parameter, r ═ 1/π Q, so the inverse Q filtering process result can be expressed as:
Figure BDA0002685364760000102
when the Q value changes continuously along with the travel time, the surface wave field can be extended to a time depth tau, and the extended wave field is as follows:
Figure BDA0002685364760000103
further modification (15) gives,
Figure BDA0002685364760000104
Figure BDA0002685364760000105
Figure BDA0002685364760000106
σ in formula (17)2For the stability factor, it limits the gain of the actual datamin(db) the relationship is:
σ2＝exp[-(0.23Gmin+1.63)] (19)
in order to improve the operation efficiency of inverse Q filtering, when inverse Q filtering is performed in the time-frequency domain using Gabor transform, equation (15) is changed to:
Figure BDA0002685364760000107
wherein the content of the first and second substances,
Figure BDA0002685364760000108
and the equation (20) is an inverse Q filtering equation based on Gabor transformation, Gabor transformation is carried out on the surface seismic records for wave fields with different time depths, inverse Q filtering is carried out on the time-frequency domain by combining an amplitude compensation operator and a phase compensation operator, and the seismic records P (tau) processed by inverse Q filtering can be obtained by utilizing inverse Gabor transformation on the obtained wave fields P (tau, omega).
Step S5, tau-p domain gather vertical energy equalization processing, wherein the step and the following 2 steps are all carried out in tau-p domain, firstly, the prestack gather is converted from x-t domain (time-space domain) to tau-p domain (Ladong domain, intercept-gradient domain) through Ladong transformation, the prestack gather is carried out vertical energy equalization in tau-p domain, and methods such as large time window AGC gain and the like can be adopted, so as to facilitate the prediction of effective reflection in the following steps, suppress multiple waves and high frequency random interference and low frequency offset arc, facilitate the gather after tau-p domain vertical energy equalization to predict signals according to the principle that the horizontal coordinate p (x) is 0 (the horizontal ray parameter p (x) is zero, and the tau-p domain represents the same phase axis signal with infinite time-space domain speed, namely the same phase axis of reflected waves), and suppress other coherent noise and residual noise in non-zero regions.
Step S6, tau-p domain coherent noise suppression: and predicting effective reflection in a tau-p domain, inverting energy of multiple waves and other non-zero regions of tau-p domain horizontal ray parameters by a least square method, suppressing the energy of the non-zero regions by a data-driven automatic weighting mode, viewing through QC, integrally reconstructing the tau-p domain signals processed in the step back to a time-space domain, improving the signal-to-noise ratio of the reconstructed trace set, effectively suppressing multiple waves and short-period inclined linear noise, and reserving the high-frequency reflection energy developed in the step 4.
Step S7, tau-p domain residual noise suppression: the phenomenon of displacement arcing caused by the existence of covering times or abnormal amplitudes in the displacement process is obvious on the CRP gather, and the arcs of the time-space domain of low-frequency and ultralow frequencies are very easy to suppress in the tau-p domain. On the basis of the step 6, the result of suppressing the over-coherent noise in the tau-p domain is inverted into low-frequency residual noise by using the least square method again, and the step needs to be noted that the low frequency of the frequency band range processed by the least square inversion process is low enough, and can be used as the processing starting frequency from 1Hz or 2Hz, so that the low frequency shift arc drawing energy is not missed as much as possible. And (4) reconstructing the signals obtained in the step 5/6/7 to obtain the high-frequency signals subjected to bluing treatment in the trace concentration after the x-t domain, suppressing low-frequency-end noise and high-frequency-end random noise and multiples, and obviously improving the high-frequency energy signal-to-noise ratio of the trace concentration after the step 2-7.
Step S8, reconstructing a time-space domain channel set by tau-p domain signals: and 5-8, performing amplitude equalization inverse processing in the step 5 on the tau-p domain signal, namely removing amplitude equalization by using the same parameters in the step 5, and recovering the original amplitude and relative relation of the prestack gather. And then, performing stability control through a complex conjugate transpose matrix to perform time-space domain reconstruction. Since tau-p transform is not completely reversible, it is possible to reduce artifacts caused by tau-p domain-time space domain data reconstruction only when the offset is infinite, but obviously, the infinite offset is an unrealizable condition, so the artifacts are unavoidable when tau-p domain-time space domain signals are reconstructed. The method has the advantages of high stability and high efficiency of data reconstruction, can basically eliminate false images caused by tau-p conversion algorithm, and can be used for eliminating interference of imaging and prestack inversion influenced by multiple waves and the like, wherein high-frequency signals of the reconstructed CRP gather are abundant, noise is obviously reduced.
The theoretical basis of steps S4-S8 is discrete parabolic radon transform, and the forward and backward transform operators of the discrete parabolic radon transform performed in the frequency domain are non-orthogonal, not true reciprocal operators, and need to be solved by generalized inverse. The conventional parabola radon transformation adopts a least square regularization method to establish least square discrete parabola radon positive transformation, and the mathematical expression of the least square discrete parabola radon positive transformation is as follows:
M＝(LLH+λ2I)-1LD (22)
in formula (22): d is f-x domain data, and M is f-q domain data;
Figure BDA0002685364760000121
Figure BDA0002685364760000122
in least squares parabolic radon transform (LL)H)-1The method can be regarded as a deconvolution operator, and can improve the resolution of the Ladong domain data to a certain extent. The method for attenuating the multiples by using the time difference of the multiples and the primaries depends on the resolution and the focusing of seismic data in a Ladong domain, and is realized by introducing a nonlinear regularization diagonal matrix W in an operator of least square parabola Ladong transformation mathematically, namely:
(LLH+WHW)M＝LD (23)
the physical meaning of the nonlinear diagonal matrix W in the formula (23) is a weight vector of the curvature parameter q, and the focusing performance of the data in the Latin field can be improved. Because discrete Radon transformation has a Toeplitz structure, least square constraint inversion Radon transformation is a frequency domain space sparse constraint algorithm, and in the inversion iteration process, a weighting matrix is connected with the result of the previous iteration through a Bayes principle according to the result of the previous iteration to obtain a new weighting matrix; then, the weighting matrix equation is solved to obtain sparse solution of the frequency domain. The reconstructed CRP gather has abundant high-frequency signals, obviously reduced noise, and elimination of interference of multiple waves and the like which influence imaging and prestack inversion.
Step S9, CRP gather high-frequency signal-to-noise ratio evaluation: the step can be carried out by solving an FK spectrum analysis on the trace set, or analyzing a frequency spectrum after FFT (fast Fourier transform) of the trace set, and carrying out qualitative and main evaluation by combining the actual condition of the trace set. If the signal-to-noise ratio at the high-frequency end is higher, the processing of the step 10 is carried out to simultaneously widen the frequency band towards high and low frequencies, and if the CRP gather high-frequency potential is not completely mined through evaluation, the blueing processing stage noise suppression and the like of the steps 4-9 can be repeated.
Step S10, time-varying pulse deconvolution processing: the core of the step is pulse deconvolution, which is the most effective means for widening the frequency band, and the maximum characteristics of the deconvolution in the invention are 3: (1) the frequency band is widened to the maximum extent, and the processing upper limit is not set; (2) all frequency components including signals and noise participate in operation, and the potential of data is mined to the maximum extent; (3) a large amount of deconvolution parameter experiments are not needed, the sampling rate is directly adopted as the prediction step length, and the wavelet consistency of the result is ensured to the maximum extent. In order to ensure that the frequency bands of the CRP gathers of the shallow and middle deep layers are widened to the maximum extent, the invention also uses a time-sharing window time-varying control method in the step, deconvolution factors are obtained in different time windows, and the current time window is independently applied, so that the frequency bands are widened to the maximum extent no matter the shallow and middle deep layers. Note that when deconvolution word calculation and application are performed in the time-sharing window, the overlapping length of the upper and lower time windows is fully considered, and the overlapping region length of every two windows is greater than or equal to 4 wavelet lengths, so as to ensure that no significant distortion occurs in the vertical direction. Deconvolution is the most common technique in seismic data processing, and is also the easiest technique to implement, and the technical principle thereof is not described herein again.
Step S11, coherent addition processing: and superposing the prestack gather according to the coherence coefficient, improving the coherence coefficient of a given high-frequency end component, and reserving the high-frequency component developed through the process as much as possible, so that the purposes of widening a frequency band after coherent superposition, improving true resolution and enhancing thin layer identification capability are achieved.
Step S12, result evaluation step: and comparing and performing spectrum analysis on the pre-stack CRP gather processed in the steps to check whether the frequency bandwidth of the result gather and the frequency bandwidth of the original CRP gather and the signal-to-noise ratio of high-frequency energy are obviously improved. The original CRP gather and the CRP gather processed by the processing method of the invention are superimposed by adopting the same method parameters, the resolution of the evaluation profile and the result of the frequency spectrum analysis are compared, and objective evaluation on the improvement degree of the gather and the improvement degree of the result profile is given.
The following describes in detail a specific application method of the fidelity frequency-lifting method for high-coverage pre-stack seismic data:
and (3) adopting the initial seismic data prestack CRP gather to perform processing step S1 method applicability analysis: the method has the advantages that the maximum covering times of the CMP gather of the target layer and the maximum covering times of the CRP gather are determined by analyzing the type of the acquisition observation system applying the target data and combining an offset distance grouping scheme during offset processing, the method is high in applicability and low in requirement on data, and the data effect is better when the covering times of the CMP gather are more than 200 times and the covering times of the CRP gather are more than 60 times (the high covering times can ensure the signal-to-noise ratio of the post-stack data).
Then, the pre-stack CRP gather is processed by steps S2 and S3, namely the CRP gather is processed by residual time difference elimination and random noise suppression: and (4) eliminating residual time difference by adopting a method based on amplitude energy or a method based on residual velocity, and flattening the same-phase axis of the CRP gather. The processing methods of the energy-based gather leveling are numerous, and the core requirement is to eliminate residual time difference as far as possible, so that the application of subsequent processing technology is facilitated to obtain expected internal effects. Step S3, CRP gather random noise suppression: on the basis of eliminating the CRP gather with the residual time difference, a processing method can be adopted, and a two-dimensional or three-dimensional space prediction error filtering technology can be referred to. The core purpose of the step is to suppress random noise on the same-phase axis leveled gather, improve the signal-to-noise ratio of the prestack gather and provide a foundation for subsequent bluing treatment, if the step of suppressing random noise is lacked, the subsequent treatment technology will widen the frequency band greatly, the high-frequency random noise energy will be amplified by exponential order, and the whole prestack fidelity frequency-extracting method will fail.
The signal-to-noise ratio of the pre-stack CRP gather obtained by the two treatments is obviously improved, and a foundation is laid for the subsequent treatment steps.
The result of the gather bluing process of step S4: the step adopts a processing means mainly based on inverse Q filtering, and aims to expand the dominant frequency band to a high-frequency end. Note that this step does not necessarily use deconvolution-like frequency broadening means, and does not broaden the frequency band towards both high and low frequencies. In the bluing process, the degree of broadening to the high band is properly controlled with reference to the number of covering times, and taking the time-varying inverse Q filtering process technique as an example, if the data spectrum frequency-down relationship is natural and there is a higher number of covering times, the Q value range can be properly reduced at the destination layer. The trace gather data processed by the step has the advantages that the effective reflection energy of the high-frequency end is improved, the difference between the level of the high-frequency reflection energy and the random noise (or background noise) of the high-frequency end is reduced, and the residual random noise in the CRP trace gather is amplified again after the step S3, so that the noise is continuously suppressed in the next step, and the signal-to-noise ratio of the high-frequency effective reflection is improved.
In order to suppress noise such as interbed multiples and offset arcs, the gather subjected to time-varying inverse Q filtering processing in step 4 based on time-frequency analysis is denoised in 4 steps (step S5-S8: step S5 and tau-p domain gather vertical energy equalization processing: the steps and the following 2 steps are all carried out in the tau-p domain, firstly, the prestack gather is converted from an x-t domain (time-space domain) to a tau-p domain (Ladong domain, intercept-gradient domain) through Radon transformation, and vertical energy equalization is carried out on the prestack gather in the tau-p domain, and methods such as large time window AGC (automatic gain control) can be adopted to facilitate gain and the like, so that effective reflection is predicted in the next step, multiple times of suppressing wave, high-frequency random interference and low-frequency offset arcs, and the gather subjected to tau-p domain vertical energy equalization is more beneficial to representing the tau-p domain by the horizontal coordinate p (horizontal ray parameter p (x) of zero The time-space domain velocity of the inphase axis signal with infinite magnitude, namely the inphase axis of the reflected wave), and suppresses other coherent noise and residual noise in a non-zero region.
Step S6, tau-p domain coherent noise suppression: and predicting effective reflection in a tau-p domain, inverting energy of multiple waves and other non-zero regions of tau-p domain horizontal ray parameters by a least square method, suppressing the energy of the non-zero regions by a data-driven automatic weighting mode, viewing through QC, integrally reconstructing the tau-p domain signals processed in the step back to a time-space domain, improving the signal-to-noise ratio of the reconstructed trace set, effectively suppressing multiple waves and short-period inclined linear noise, and reserving the high-frequency reflection energy developed in the step 4.
Step S7, tau-p domain residual noise suppression: the phenomenon of displacement arcing caused by the existence of covering times or abnormal amplitudes in the displacement process is obvious on the CRP gather, and the arcs of the time-space domain of low-frequency and ultralow frequencies are very easy to suppress in the tau-p domain. On the basis of the step 6, the result of suppressing the over-coherent noise in the tau-p domain is inverted into low-frequency residual noise by using the least square method again, and the step needs to be noted that the low frequency of the frequency band range processed by the least square inversion process is low enough, and can be used as the processing starting frequency from 1Hz or 2Hz, so that the low frequency shift arc drawing energy is not missed as much as possible. And (4) reconstructing the signals obtained in the step 5/6/7 to obtain the high-frequency signals subjected to bluing treatment in the trace concentration after the x-t domain, suppressing low-frequency-end noise and high-frequency-end random noise and multiples, and obviously improving the high-frequency energy signal-to-noise ratio of the trace concentration after the step 2-7.
Step S8, reconstructing a time-space domain channel set by tau-p domain signals: and 5-7, performing amplitude equalization inverse processing in the step 5 on the tau-p domain signal, namely removing amplitude equalization by using the same parameters in the step 5, and recovering the original amplitude and relative relation of the prestack gather. And then, performing stability control through a complex conjugate transpose matrix to perform time-space domain reconstruction. Since tau-p transform is not completely reversible, it is possible to reduce artifacts caused by tau-p domain-time space domain data reconstruction only when the offset is infinite, but obviously, the infinite offset is an unrealizable condition, so the artifacts are unavoidable when tau-p domain-time space domain signals are reconstructed. The method has the advantages of high stability and high efficiency of data reconstruction, can basically eliminate false images caused by tau-p conversion algorithm, and can be used for eliminating interference of imaging and prestack inversion influenced by multiple waves and the like, wherein high-frequency signals of the reconstructed CRP gather are abundant, noise is obviously reduced.
The interbed multiples of the prestack CRP gather after the processing steps S5-S8 are suppressed, and the low-frequency and ultra-low-frequency noise is also effectively suppressed.
Step S9, CRP gather high-frequency signal-to-noise ratio evaluation: the step can be carried out by solving an FK spectrum analysis on the trace set, or analyzing a frequency spectrum after FFT (fast Fourier transform) of the trace set, and carrying out qualitative and main evaluation by combining the actual condition of the trace set. If the signal-to-noise ratio at the high-frequency end is higher, the processing of the step S10 is carried out to simultaneously widen the frequency band towards high and low frequencies, if the CRP gather high-frequency potential is not completely mined by evaluation, the blue processing level noise suppression and other processing of the step 4-the step S9 can be repeated,
step S10, time-varying pulse deconvolution processing: the core of the step is pulse deconvolution, which is the most effective means for widening the frequency band, and the maximum characteristics of the deconvolution in the invention are 3: (1) the frequency band is widened to the maximum extent, and the processing upper limit is not set; (2) all frequency components including signals and noise participate in operation, and the potential of data is mined to the maximum extent; (3) a large amount of deconvolution parameter experiments are not needed, the sampling rate is directly adopted as the prediction step length, and the wavelet consistency of the result is ensured to the maximum extent. In order to ensure that the frequency bands of the CRP gathers of the shallow and middle deep layers are widened to the maximum extent, the invention also uses a time-sharing window time-varying control method in the step, deconvolution factors are obtained in different time windows, and the current time window is independently applied, so that the frequency bands are widened to the maximum extent no matter the shallow and middle deep layers. Note that when deconvolution word calculation and application are performed in the time-sharing window, the overlapping length of the upper and lower time windows is fully considered, and the overlapping region length of every two windows is greater than or equal to 4 wavelet lengths, so as to ensure that no significant distortion occurs in the vertical direction.
Step S11, coherent addition processing: and superposing the prestack gather according to the coherence coefficient, improving the coherence coefficient of a given high-frequency end component, and reserving the high-frequency component developed through the process as much as possible, so that the purposes of widening a frequency band after coherent superposition, improving true resolution and enhancing thin layer identification capability are achieved.
Step S12, result evaluation step: and comparing and performing spectrum analysis on the pre-stack CRP gather processed in the steps to check whether the frequency bandwidth of the result gather and the frequency bandwidth of the original CRP gather and the signal-to-noise ratio of high-frequency energy are obviously improved. The original CRP gather and the CRP gather processed by the processing method of the invention are superimposed by adopting the same method parameters, the resolution of the evaluation profile and the result of the frequency spectrum analysis are compared, and objective evaluation on the improvement degree of the gather and the improvement degree of the result profile is given.
Furthermore, the inventors have studied the possible effects of certain elements of the invention on the final results when they are changed, and these conclusions should be patented, and have the following recognition and conclusion: a. before the frequency extraction of the method, technical applicability evaluation is required, and if the pre-stack CMP covering times are more than or equal to 200 times, or the CRP gather covering times are more than or equal to 60 times, the method is more applicable; the method of the invention is equally applicable if the pre-stack CRP gathers have a higher signal-to-noise ratio if the maximum coverage times are lower for older data.
b. Before the pre-stack CRP gather bluing treatment adopts Time-varying inverse Q filtering, Time-window spectrum analysis is required, Time-Q value parameters used in the step need to be tested repeatedly, the principle is to widen the frequency band of high-frequency components after the inverse Q filtering, but a double-peak-shaped spectrum can not appear in a target layer after the inverse Q filtering, and if the situation shows that the Q value in the Time window is too small, the Q value in the Time window applied to the analysis is required to be amplified properly. The method of performing spectrum analysis on the time-sharing window before the pre-stack inverse Q filtering and performing time-sharing window evaluation on the result after the inverse Q filtering is protected.
c. When multiple wave interference existing in a pre-stack trace set is inverted in a tau-p domain by adopting a least square inversion method, the iteration number of the least square inversion of data with the recorded trace length of 5s is not less than 30.
d. Taking data with a recording track length of 5s as an example, the first arrival intensity energy removal process is required before time-varying pulse deconvolution, the number of analysis time windows and application time windows and the overlapping time of adjacent time windows must be strictly consistent, the size of a time window is more than 4 times of the wavelength of an analysis wavelet, if the length of the analysis wavelet is 120ms, a single analysis time window needs to be more than 500ms, and in order to meet the influence caused by waveform distortion frequency reduction in the wavelet propagation process, the number of the analysis time windows should not be less than one fifth of the recording track length, namely, pre-stack CRP data with a 5s track length, the analysis application time windows should be more than 5.
The above embodiments are only for illustrating the invention and are not to be construed as limiting the invention, and those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention, therefore, all equivalent technical solutions also belong to the scope of the invention, and the scope of the invention is defined by the claims.

Claims (6)

1. A prestack fidelity frequency-boosting method for high-coverage seismic data is characterized by comprising the following steps:
s1, judging the applicability of the seismic data to be processed pre-stack gather data, if the applicability meets the condition, performing S2, and if the applicability does not meet the condition, terminating the operation;
step S2, residual time difference elimination is carried out on the CRP gather before folding;
step S3, carrying out random noise suppression treatment based on prediction error on the pre-stack CRP gather which is processed in the step S2 and eliminates the residual time difference;
step S4, performing stability inverse Q filtering processing based on time-frequency analysis on the pre-stack CRP gather processed in the step S3;
s5, carrying out weighted least square Radon transformation denoising processing on the pre-stack CRP gather processed in the S4;
step S6, performing time-varying pulse deconvolution processing on the pre-stack CRP gathers processed in the step S5;
and S7, performing denoising and coherent superposition processing on the pre-stack CRP gather processed in the step S6.
2. The method according to claim 1, wherein the step S1 of determining the applicability of the pre-stack gather data of the seismic data to be processed specifically comprises:
when the number of times of coverage of the CMP gather of the seismic data to be processed is more than or equal to 200 times, or the number of times of coverage of the CRP gather is more than or equal to 60 times, executing the step S2;
the operation is terminated when the number of CMP gather overlays is less than 200 or the number of CRP gather overlays is less than 60.
3. The high coverage seismic data prestack fidelity frequency extraction method of claim 1,
in step S2, the residual moveout elimination for the pre-stack CRP gather specifically includes:
when the residual time difference of the target layer in-phase axis of the pre-stack CRP gather is greater than or equal to 20ms when the reference offset distance is 3000m, the processing is terminated, and the velocity analysis is carried out again;
when the residual time difference is less than 20ms and more than or equal to 10ms, adjusting the speed by using a speed analysis technology based on spatial continuity, and eliminating the residual time difference;
when the residual time difference is less than 10ms, eliminating the residual time difference by using a spatial weighted median filtering method, wherein the output amplitude of each sampling point of the CRP gather in-phase axis after the spatial weighted median filtering is as follows:
Figure FDA0002685364750000021
wherein, VKIs the non-abnormal amplitude value in the aperture, N is the number of non-abnormal amplitude values in the aperture, WeightkIs the space weight;
after the treatment, the residual time difference of the CRP gather is eliminated.
4. The method according to claim 1, wherein the step S4 of performing the stability inverse Q filtering based on the time-frequency analysis specifically comprises:
Figure FDA0002685364750000022
Figure FDA0002685364750000023
the formula is an inverse Q filtering equation based on Gabor transformation, tau is depth, j is an imaginary unit, omega is angular frequency, Q is stratum quality factor, for wave fields with different time depths, Gabor transformation is carried out on surface seismic records, then inverse Q filtering is carried out in a time-frequency domain by combining an amplitude compensation operator and a phase compensation operator, and the seismic records P (tau) processed by inverse Q filtering can be obtained by utilizing inverse Gabor transformation on the obtained wave fields P (tau, omega).
5. The method for high coverage seismic data prestack fidelity frequency extraction as claimed in claim 1, wherein the step S5 of performing weighted least squares radon transform denoising process specifically includes:
step S51, least squares radon transform, comprising: the least square regularization method is adopted in a frequency domain, the least square discrete parabola Ladong positive transformation is established, and the mathematical expression is as follows:
M＝(LLH+λ2I)-1LD
in the formula: d is f-x domain data, M is f-q domain data, (LL)H)-1Is a deconvolution operator; wherein
Figure FDA0002685364750000031
Figure FDA0002685364750000032
Converting the pre-stack seismic data of the time-space domain into a tau-p domain by using the method;
step S52, tau-p domain weighted coherent noise suppression, comprising: introducing a nonlinear regularization diagonal matrix W into an operator of least square parabola Radon transform to realize, and attenuating tau-p domain coherent noise by time difference of multiples and primary waves: (LL)H+WHW)M＝LD
In the formula, the physical meaning of a nonlinear diagonal matrix W is a weight vector of a curvature parameter q, effective reflection is predicted in a tau-p domain, the energy of multiple coherent noise in a tau-p domain horizontal ray parameter nonzero region is inverted through a least square method, the energy of the nonzero region is suppressed in a data-driven automatic weighting mode, and the coherent noise is suppressed in the tau-p domain;
step S53, tau-p domain residual noise suppression: in step S52, the weighted least-squares radon transform noise suppression is performed again on the full-frequency data other than the dominant frequency band by the same method.
And step S54, reconstructing the tau-p domain data processed in the steps back to a time-space domain through the control of a complex conjugate transpose matrix, and finishing the whole processing.
6. The method according to claim 1, wherein the step S6 of performing the variable pulse deconvolution includes:
calculating a deconvolution operator a (t) by using the basic equation of the pulse deconvolution,
calculating a deconvolution operator a (t) by using a fundamental equation of deconvolution,
Figure FDA0002685364750000033
in the formula, rxx(m) is the seismic reflection coefficient, a (t) is the pulse deconvolution operator; utilizing convolution of a '(t) and seismic record x (t) to obtain S (t) ═ a' (t) × (t), wherein S (t) is seismic trace output by pulse deconvolution, performing windowing calculation on input seismic record in the processing process to obtain pulse deconvolution operator a (t), wherein each 2 analysis calculation time windows have an overlapping transition region with a certain length, and the length of the window overlapping region is more than or equal to 4 seismic wavelet lengths; performing pulse deconvolution in a frequency domain by participating all frequency components of an input seismic record in time-varying pulse deconvolution calculation; the deconvolution prediction step size is equal to the sampling interval of the input data.
CN202010974716.7A 2020-09-16 2020-09-16 Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data Active CN112213775B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010974716.7A CN112213775B (en) 2020-09-16 2020-09-16 Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010974716.7A CN112213775B (en) 2020-09-16 2020-09-16 Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data

Publications (2)

Publication Number Publication Date
CN112213775A true CN112213775A (en) 2021-01-12
CN112213775B CN112213775B (en) 2023-01-24

Family

ID=74050044

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010974716.7A Active CN112213775B (en) 2020-09-16 2020-09-16 Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data

Country Status (1)

Country Link
CN (1) CN112213775B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114185095A (en) * 2021-12-02 2022-03-15 中国石油大学(北京) Method for suppressing multiple waves of three-dimensional plane wave domain seismic data

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4964102A (en) * 1988-04-19 1990-10-16 Amoco Corporation Method for enhancing and evaluating seismic data
US20040122596A1 (en) * 2002-12-19 2004-06-24 Core Laboratories, Inc. Method for high frequency restoration of seismic data
US20070064530A1 (en) * 2005-08-26 2007-03-22 Ian Moore Method for processing a record of seismic traces
CN105093291A (en) * 2014-05-14 2015-11-25 中国石油天然气股份有限公司 Method for recovering oil-gas reservoir seismic reflection features
CN106371140A (en) * 2016-08-17 2017-02-01 中国石油化工股份有限公司 Method for raising the resolution of high, middle and deep earthquake data
CN106597532A (en) * 2016-11-14 2017-04-26 中国石油化工股份有限公司 Pre-stack seismic data frequency band expanding method of combining well information and horizon information
WO2017108690A1 (en) * 2015-12-22 2017-06-29 Shell Internationale Research Maatschappij B.V. Method and system for separating blended seismic data
US10191168B2 (en) * 2016-01-12 2019-01-29 Cgg Services Sas AVA compliant pre-stack frequency spectrum enhancement of seismic data
CN109307887A (en) * 2017-07-28 2019-02-05 中国石油化工股份有限公司 The weak reflector recognition methods of earthquake and system
CN110471113A (en) * 2019-08-01 2019-11-19 中国石油大学(北京) Bearing calibration, device and storage medium are moved in inverting based on unstable state seismic data
US20200124752A1 (en) * 2018-10-17 2020-04-23 Anatoly I. Baumstein Accurate Velocity Model Estimation and Imaging in the Presence of Localized Attenuation (Q) Anomalies

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4964102A (en) * 1988-04-19 1990-10-16 Amoco Corporation Method for enhancing and evaluating seismic data
US20040122596A1 (en) * 2002-12-19 2004-06-24 Core Laboratories, Inc. Method for high frequency restoration of seismic data
US20070064530A1 (en) * 2005-08-26 2007-03-22 Ian Moore Method for processing a record of seismic traces
CN105093291A (en) * 2014-05-14 2015-11-25 中国石油天然气股份有限公司 Method for recovering oil-gas reservoir seismic reflection features
WO2017108690A1 (en) * 2015-12-22 2017-06-29 Shell Internationale Research Maatschappij B.V. Method and system for separating blended seismic data
US10191168B2 (en) * 2016-01-12 2019-01-29 Cgg Services Sas AVA compliant pre-stack frequency spectrum enhancement of seismic data
CN106371140A (en) * 2016-08-17 2017-02-01 中国石油化工股份有限公司 Method for raising the resolution of high, middle and deep earthquake data
CN106597532A (en) * 2016-11-14 2017-04-26 中国石油化工股份有限公司 Pre-stack seismic data frequency band expanding method of combining well information and horizon information
CN109307887A (en) * 2017-07-28 2019-02-05 中国石油化工股份有限公司 The weak reflector recognition methods of earthquake and system
US20200124752A1 (en) * 2018-10-17 2020-04-23 Anatoly I. Baumstein Accurate Velocity Model Estimation and Imaging in the Presence of Localized Attenuation (Q) Anomalies
CN110471113A (en) * 2019-08-01 2019-11-19 中国石油大学(北京) Bearing calibration, device and storage medium are moved in inverting based on unstable state seismic data

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
YI BAO ET AL.: "Phase Stabilization Prestack time Q migration of BWH seismic data in northern Songliao Basin", 《2019 SEG INTERNATIONAL EXPOSITION AND 89TH ANNUAL MEETING》 *
YIBAO ET AL.: "Deep Formation Phase Stabilization Q Migration of WBH seismic data on Songliao", 《SEG 2018 WORKSHOP: SEG SEISMIC IMAGING WORKSHOP》 *
丁进杰等: "基于高频补偿方法提高地震资料", 《地球物理学进展》 *
刘志伟等: "薄互层地震成像中高分辨率处理方法", 《地球物理学报》 *
包燚等: "大庆油田" 两宽一高" 地震资料稳相 Q 偏移", 《大庆石油地质与开发》 *
徐泰然等: "深反射地震剖面数据处理的主要技术方法", 《地球物理学进展》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114185095A (en) * 2021-12-02 2022-03-15 中国石油大学(北京) Method for suppressing multiple waves of three-dimensional plane wave domain seismic data

Also Published As

Publication number Publication date
CN112213775B (en) 2023-01-24

Similar Documents

Publication Publication Date Title
EP3193193B1 (en) Ava-compliant enhancement of pre-stack frequency spectrum of seismic data
US5173879A (en) Surface-consistent minimum-phase deconvolution
Yuan et al. Stable inversion-based multitrace deabsorption method for spatial continuity preservation and weak signal compensation
US4884247A (en) Method of processing geophysical data to compensate for earth filter attenuation
US20150168573A1 (en) Geologic quality factor inversion method
CN102043165B (en) Basis tracking algorithm-based surface wave separation and suppression method
CN102288994A (en) Method for regularizing high-dimensional seismic data under constraint of Radon spectrum
CN112305612B (en) High-resolution complex spectrum decomposition time-frequency space domain amplitude variation correction method along with offset distance
Liu et al. Seismic quality factor estimation using frequency-dependent linear fitting
CN104635264B (en) The processing method of earthquake data before superposition and equipment
CN112213775B (en) Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data
Wang et al. Robust singular value decomposition filtering for low signal-to-noise ratio seismic data
Li et al. A generalized seismic attenuation compensation operator optimized by 2-D mathematical morphology filtering
Bai et al. A U-Net based deep learning approach for seismic random noise suppression
CN111257938A (en) Time-lapse seismic virtual source wave field reconstruction method and system based on wavelet cross-correlation
CN104820244A (en) Method for improving signal-to-noise ratio in processing petroleum exploration data
CN107678065B (en) The guarantor for improving seismic resolution constructs well control space the Method of Deconvolution and device
CN114137606A (en) Stable spectrum simulation deconvolution method
Bellezza et al. Multidimensional deconvolution and processing of seismic-interferometry Arctic data
Bakulin et al. Bootstrapping invisible signals: Prestack land data enhancement using nonlinear beamforming with local waveform corrections
US9014985B2 (en) System and method for compensating time and offset varying near-surface effects in seismic data background
Wang et al. Inversion-based non-stationary normal moveout correction along with prestack high-resolution processing
CN117434606B (en) Seismic section denoising method based on improved Laplace filtering reverse time migration
Bao et al. A new technique to support future energy exploration of continental sedimentary basin: BWH full-frequency fidelity and amplitude preserving processing technology
Wang et al. Structure-oriented edge-preserving smoothing based on accurate estimation of orientation and edges

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant