CN104536044A - Interpolation and denoising method and system for seismic data - Google Patents

Abstract

The invention provides an interpolation and denoising method and system for seismic data. The method comprises the steps that a preset sampling matrix is acquired; observation data are collected according to the sampling matrix; an objective functional is constructed; according to the observation data and the objective functional, seismic data volume, after interpolation and denoising are performed, corresponding to the observation data is determined. By studying an interpolation and denoising integration scheme of low-signal-to-noise-ratio seismic data, denoising is performed through a threshold value after each time of interpolation, and finally the complete and regular seismic data volume high in signal to noise ratio is obtained. Accordingly, the interpolation and denoising method and system lay a data foundation for multiple processes such as multiple attenuation based on wave equation surface correlation, migration imaging and wave equation inversion.

Description

A kind of interpolation denoising method of geological data and system

Technical field

The present invention, about technical field of geophysical exploration, particularly about the Seismic Data Processing Technique in exploration of oil and gas field field, is a kind of interpolation denoising method and system of geological data concretely.

Background technology

In exploration of oil and gas field field, due to acquisition cost restriction and prohibit exploiting field, barrier, the drift of marine streamer pinniform impact, make gather geological data be irregular at direction in space.Secondly, the rejecting that in processing procedure, blown-out shot gives up also is cause one of irregular reason of data.And in follow-up process of seismic data processing, based on the compacting of wave equation surface-related multiple, migration imaging and wave equation inversion, need based on complete rule and the higher seismic data volume of signal to noise ratio (S/N ratio).Therefore, geological data interpolation is an important step in seism processing.

In prior art, numerous disposal routes is generally that the interpolation of geological data is separated with denoising and performed, but to realize geological data interpolation prerequisite be that the signal to noise ratio (S/N ratio) of geological data is higher, and geological data disappearance affects denoising effect, and both exist contradiction.Therefore, seismic data with low signal-to-noise ratio interpolation is adopted usually to the mode of weighting, its treatment effect is unsatisfactory.

Therefore, how to develop a kind of interpolation denoising scheme of new geological data, it farthest can rebuild missing data, removes the random noise in section simultaneously, and the complete regular geological data finally obtaining high s/n ratio is this area technical barrier urgently to be resolved hurrily.

Summary of the invention

In order to overcome the above-mentioned technical barrier that prior art exists, the invention provides a kind of interpolation denoising method and system of geological data, by studying interpolation and the denoising integrated programme of seismic data with low signal-to-noise ratio, after each interpolation processing, denoising is carried out by threshold value, finally obtain the complete regular seismic data volume of high s/n ratio, for establishing data basis based on multiple tracks process such as the compacting of wave equation surface-related multiple, migration imaging and wave equation inversions.

An object of the present invention is, provides a kind of interpolation denoising method of geological data, comprising: obtain the sampling matrix preset; Observation data is gathered according to described sampling matrix; Establishing target functional; The seismic data volume after the interpolation denoising that described observation data is corresponding is determined according to described observation data and described cost functional.

In a preferred embodiment of the invention, the cost functional of structure is:

$\mathrm{\Φ}\left(x\right)={\left|\right|d_{\mathrm{obs}}-{\mathrm{RC}}^{T}x\left|\right|}_2^2+\mathrm{\λP}\left(x\right)$

Wherein, Φ (x) is cost functional, and x is bent wave system number vector, C ^tfor bent ripple inverse transformation, C is warp wavelet, and λ is regularization factors, and P (x) is sparse constraint, d _obsfor observation data, R is sampling matrix.

In a preferred embodiment of the invention, the seismic data volume after determining according to described observation data and described cost functional the interpolation denoising that described observation data is corresponding comprises: according to described observation data definite threshold; According to described threshold value and described cost functional determination hard threshold function; According to described threshold value and described cost functional determination soft-threshold function; Build interpolation denoising iterative equation; Solve described cost functional according to described observation data, hard threshold function, soft-threshold function and interpolation denoising iterative equation, obtain the seismic data volume after interpolation denoising.

In a preferred embodiment of the invention, comprise according to described observation data definite threshold: adopt sparse transformation warp wavelet that described observation data is transformed to bent wave zone; Observation data described in acquisition transforms to the maximum bent wave system number of the amplitude after bent wave zone; Obtain the weight limit factor, the minimal weight factor and the maximum iteration time that preset; According to the described weight limit factor and the maximum bent wave system number determination max-thresholds of amplitude; According to the described minimal weight factor and the maximum bent wave system number determination minimum threshold of amplitude; According to described minimum threshold, max-thresholds, maximum iteration time and exponential function definite threshold.

In a preferred embodiment of the invention, described exponential function is:

τ _k＝τ _maxe ^c(k-1)/(N-1)

c＝ln(τ _min/τ _max)

Wherein, N is maximum iteration time, k=1,2 ..., N, τ _kfor the threshold value of kth time iteration, τ _maxfor max-thresholds, τ _minfor minimum threshold, c is intermediate variable.

In a preferred embodiment of the invention, solve described cost functional according to described observation data, hard-threshold, soft-threshold and interpolation denoising iterative equation, the seismic data volume obtained after interpolation denoising comprises: iterative initial value and the weight factor of determining data space territory; According to described weight factor and described observation data determining section plane of vision; According to the solution that the iterative initial value in data space territory, weight factor and interpolation denoising iterative equation determination data space territory iteration produce; The solution that described data space territory iteration produces is projected on described part plane of vision; Adopt threshold strategies to remove noise respective components at bent wave zone, obtain the solution in model space territory; Solution in model space territory is projected to data space territory, obtains the seismic data volume after interpolation denoising.

An object of the present invention is, provides a kind of system of interpolation denoising of geological data, comprising: sampling matrix acquisition device, for obtaining the sampling matrix preset; Observation data harvester, for the observation data that the sampling matrix acquiring seismic data according to described is corresponding; Cost functional construction device, for establishing target functional; Seismic data volume determining device, for determining the seismic data volume after the interpolation denoising that described observation data is corresponding according to described observation data and described cost functional.

In a preferred embodiment of the invention, described seismic data volume determining device comprises: threshold determination module, for according to described observation data definite threshold; Hard threshold function determination module, for according to described threshold value and described cost functional determination hard threshold function; Soft-threshold function determination module, for according to described threshold value and described cost functional determination soft-threshold function; Iterative equation builds module, for building interpolation denoising iterative equation; Seismic data volume determination module, for solving described cost functional according to described observation data, hard threshold function, soft-threshold function and interpolation denoising iterative equation, obtains the seismic data volume after interpolation denoising.

In a preferred embodiment of the invention, described threshold determination module comprises: converter unit, for adopting sparse transformation warp wavelet, described observation data is transformed to bent wave zone; Bent wave system number acquiring unit, transforms to the maximum bent wave system number of the amplitude after bent wave zone for obtaining described observation data; Acquiring unit, for obtaining the weight limit factor, the minimal weight factor and the maximum iteration time that preset; Max-thresholds determining unit, for according to the described weight limit factor and the maximum bent wave system number determination max-thresholds of amplitude; Minimum threshold determining unit, for according to the described minimal weight factor and the maximum bent wave system number determination minimum threshold of amplitude; Threshold value determination unit, for according to described minimum threshold, max-thresholds, maximum iteration time and exponential function definite threshold.

In a preferred embodiment of the invention, described seismic data volume determination module comprises: iterative initial value determining unit, for determining iterative initial value and the weight factor in data space territory; Part plane of vision determining unit, for according to described weight factor and described observation data determining section plane of vision; Determining unit is separated in data space territory, for the solution that the iterative initial value according to data space territory, weight factor and interpolation denoising iterative equation determination data space territory iteration produce; Projecting cell, projects to described part plane of vision for the solution described data space territory iteration produced; Determining unit is separated in model space territory, for adopting threshold strategies to remove noise respective components at bent wave zone, obtains the solution in model space territory; Seismic data volume determining unit, for the solution in model space territory is projected to data space territory, obtains the seismic data volume after interpolation denoising.

Beneficial effect of the present invention is, provides a kind of interpolation denoising method and system of geological data, farthest rebuilds missing data, removes the random noise in section simultaneously, finally obtains the complete regular geological data of high s/n ratio.The present invention is after have studied geological data interpolation problem and seismic data noise attenuation problem, and the basis in conjunction with interpolation, Denoising Problems common feature proposes, and interpolation problem only considers data reconstruction, often affected by noise, makes reconstruction effect undesirable; Shortage of data generally can affect the effect of denoising, therefore the present invention is directed to the disappearance geological data of low signal-to-noise ratio, propose the method for interpolation, denoising integration, denoising is carried out by threshold strategies after each interpolation processing, obtain the complete regular data body of high s/n ratio after iteration terminates, relatively conventional interpolation method, the method has stronger noise immunity, the signal to noise ratio (S/N ratio) of reconstructed results is high, reaches the advantage that interpolation denoising processes simultaneously.

For above-mentioned purpose of the present invention, feature and advantage can be become apparent, preferred embodiment cited below particularly, and coordinate institute's accompanying drawings, be described in detail below.

Accompanying drawing explanation

In order to be illustrated more clearly in the embodiment of the present invention or technical scheme of the prior art, be briefly described to the accompanying drawing used required in embodiment or description of the prior art below, apparently, accompanying drawing in the following describes is only some embodiments of the present invention, for those of ordinary skill in the art, under the prerequisite not paying creative work, other accompanying drawing can also be obtained according to these accompanying drawings.

The process flow diagram of the interpolation denoising method of a kind of geological data that Fig. 1 provides for the embodiment of the present invention;

Fig. 2 is the particular flow sheet of the step S104 in Fig. 1;

Fig. 3 is the particular flow sheet of the step S201 in Fig. 2;

Fig. 4 is the particular flow sheet of the step S205 in Fig. 2;

The structured flowchart of the interpolation denoising system of a kind of geological data that Fig. 5 provides for the embodiment of the present invention;

The structured flowchart of seismic data volume determining device in the interpolation denoising system of a kind of geological data that Fig. 6 provides for the embodiment of the present invention;

The structured flowchart of the threshold determination module in the interpolation denoising system of a kind of geological data that Fig. 7 provides for the embodiment of the present invention;

The structured flowchart of the seismic data volume determination module in the interpolation denoising system of a kind of geological data that Fig. 8 provides for the embodiment of the present invention;

Fig. 9 is the process flow diagram of geological data interpolation denoising integral method in specific embodiment provided by the invention;

Figure 10 is the schematic diagram of simulated earthquake data in specific embodiment provided by the invention;

Figure 11 is the schematic diagram of missing data in specific embodiment provided by the invention;

Figure 12 is the schematic diagram of noisy simulated data in specific embodiment provided by the invention;

Figure 13 is the schematic diagram of noisy missing data corresponding in specific embodiment provided by the invention;

Figure 14 is the schematic diagram of actual oceanographic data in specific embodiment provided by the invention;

Figure 15 is disappearance oceanographic data schematic diagram noisy in specific embodiment provided by the invention;

Figure 16 is the signal to noise ratio (S/N ratio) curve synoptic diagram that weighting convex set projection P OCS method is rebuild;

Figure 17 is the signal to noise ratio (S/N ratio) curve synoptic diagram rebuild in specific embodiment provided by the invention;

Figure 18 is convex set projection POCS method noisy simulated data reconstructed results schematic diagram;

Figure 19 is residual error schematic diagram corresponding to convex set projection POCS method;

Figure 20 is weighting convex set projection P OCS method noisy simulated data reconstructed results schematic diagram;

Figure 21 is residual error schematic diagram corresponding to weighting convex set projection P OCS method;

Figure 22 is noisy simulated data reconstructed results schematic diagram in specific embodiment provided by the invention;

Figure 23 is residual error schematic diagram corresponding in specific embodiment provided by the invention;

Figure 24 is convex set projection POCS method noisy real data reconstructed results schematic diagram;

Figure 25 is residual error schematic diagram corresponding to convex set projection POCS method;

Figure 26 is weighting convex set projection P OCS method noisy real data reconstructed results schematic diagram;

Figure 27 is residual error schematic diagram corresponding to weighting convex set projection P OCS method;

Figure 28 is noisy real data reconstructed results schematic diagram in specific embodiment provided by the invention;

Figure 29 is residual error schematic diagram corresponding in specific embodiment provided by the invention.

Embodiment

Below in conjunction with the accompanying drawing in the embodiment of the present invention, be clearly and completely described the technical scheme in the embodiment of the present invention, obviously, described embodiment is only the present invention's part embodiment, instead of whole embodiments.Based on the embodiment in the present invention, those of ordinary skill in the art, not making the every other embodiment obtained under creative work prerequisite, belong to the scope of protection of the invention.

The present invention relates to exploration of oil and gas field field For Processing Seismic Data, particularly a kind of geological data interpolation denoising integral treatment method, the rule finally obtaining high s/n ratio rebuilds seismic data volume.

The process flow diagram of the interpolation denoising method of a kind of geological data that Fig. 1 provides for the embodiment of the present invention, as shown in Figure 1, described method comprises:

S101: obtain the sampling matrix preset.In the specific embodiment of the present invention, sampling matrix R represents.

S102: gather observation data according to described sampling matrix.Fundamental purpose of the present invention is that the observation data gathered according to this step constructs complete geological data, the seismic data volume namely after interpolation denoising.

In the specific embodiment of the present invention, observation data d _obsrepresent, complete geological data d ₀represent.Observation data d _obswith partial data d ₀and the relation of random noise n is as shown in formula (1):

d _obs＝Rd ₀+n (1)

Figure 10 is the schematic diagram of simulated earthquake data in specific embodiment provided by the invention, Figure 11 is the schematic diagram of missing data in specific embodiment provided by the invention, Figure 12 is the schematic diagram of noisy simulated data in specific embodiment provided by the invention, Figure 13 is the schematic diagram of missing data corresponding in specific embodiment provided by the invention, Figure 14 is the schematic diagram of actual oceanographic data in specific embodiment provided by the invention, and Figure 15 is disappearance oceanographic data schematic diagram noisy in specific embodiment provided by the invention.

S103: establishing target functional.

Due to the impact of shortage of data and noise, solving of equation (1) is ill posed, and consider openness in warp wavelet territory of geological data, adopt sparse promotion strategy, establishing target functional is as follows,

$\mathrm{\Φ}\left(x\right)={\left|\right|d_{\mathrm{obs}}-{\mathrm{RC}}^{T}x\left|\right|}_2^2+\mathrm{\λP}\left(x\right)---\left(2\right)$

Wherein, Φ (x) is cost functional, and x is bent wave system number vector, C ^tfor bent ripple inverse transformation, C is warp wavelet, and λ is regularization factors, P (x) for sparse constraint (|| x|| ₀or || x|| ₁), d _obsfor observation data, R is sampling matrix.

S104: determine the seismic data volume after the interpolation denoising that described observation data is corresponding according to described observation data and described cost functional.Fig. 2 is the particular flow sheet of step S104, and as shown in Figure 2, this step specifically comprises:

S201: according to described observation data definite threshold.Fig. 3 is the particular flow sheet of step S201, and as shown in Figure 3, step S201 specifically comprises:

S301: adopt sparse transformation warp wavelet that described observation data is transformed to bent wave zone, namely utilize sparse transformation warp wavelet, observation data is transformed to bent wave zone.

S302: the observation data described in acquisition transforms to the maximum bent wave system number of the amplitude after bent wave zone, in a particular embodiment, selects the bent wave system number C that amplitude is maximum _maxas a reference.

S303: obtain the weight limit factor, the minimal weight factor and the maximum iteration time that preset, in a particular embodiment, the weight limit factor is p ₁, the minimal weight factor is p ₂, maximum iteration time is N.

S304: according to the described weight limit factor and the maximum bent wave system number determination max-thresholds of amplitude, i.e. max-thresholds τ _max=p ₁c _max.

S305: according to the described minimal weight factor and the maximum bent wave system number determination minimum threshold of amplitude, i.e. minimum threshold τ _min=p ₂c _max.

S306: according to described minimum threshold, max-thresholds, maximum iteration time and exponential function definite threshold.

In a particular embodiment, given maximum iteration time N, then threshold value can be determined by exponential function, concrete formula as shown in formula (3),

τ _k＝τ _maxe ^c(k-1)/(N-1)(3)

c＝ln(τ _min/τ _max)

Wherein, N is maximum iteration time, k=1,2 ..., N, τ _kfor the threshold value of kth time iteration, τ _maxfor max-thresholds, τ _minfor minimum threshold, c is intermediate variable.

As k=N, then τ _n=τ _min=p ₂c _max.

As k=1, then τ ₁=τ _max=p ₁c _max.

As shown in Figure 2, step S104 also comprises:

S202: according to described threshold value and described cost functional determination hard threshold function.

After threshold value is determined in a particular embodiment, when sparse constraint P (x) in cost functional=|| x|| ₀namely sparse constraint P (x) is L ₀during norm, hard threshold function T _λx () is adopted, it meets

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x,& \left|x\right|\&GreaterEqual;\mathrm{\&tau;}\\ 0,& \left|x\right|<\mathrm{\&tau;}\end{array}\right.,---\left(4\right)$

Wherein, for threshold value.

S203: according to described threshold value and described cost functional determination soft-threshold function.

In a particular embodiment, when sparse constraint P (x) in cost functional=|| x|| ₁namely sparse constraint P (x) is L ₁during norm, soft-threshold function T _λx () meets

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&tau;},& x>\mathrm{\&tau;}\\ 0,& \left|x\right|\&le;\mathrm{\&tau;}\\ x+\mathrm{\&tau;},& x<-\mathrm{\&tau;}\end{array}\right.---\left(5\right)$

Wherein, τ=0.5 λ is threshold value.

S204: build interpolation denoising iterative equation.Interpolation denoising iterative equation is as shown in formula (6):

$\begin{array}{c}{\tilde{d}}_k=\mathrm{\&alpha;}{d}_{\mathrm{obs}}+(I-\mathrm{\&alpha;R})d_{k-1}\\ x_k=T_{{\mathrm{\&lambda;}}_k}\left(C({\tilde{d}}_k+(1-\mathrm{\&alpha;})(d_{\mathrm{obs}}-{\mathrm{Rd}}_k))\right)\end{array}.---\left(6\right)$

d _k＝C ^Tx _k

Wherein, d _k=C ^tx _kthe solution in data space, it is the solution after part observation data is inserted.

S205: solve described cost functional according to described observation data, hard threshold function, soft-threshold function and interpolation denoising iterative equation, obtain the seismic data volume after interpolation denoising.Fig. 4 is the particular flow sheet of step S205, and as shown in Figure 4, this step specifically comprises:

S401: iterative initial value and the weight factor of determining data space territory, in a particular embodiment, make k=1, selects the iterative initial value d in data space territory ₀=d _obs, weight factor α ∈ [0,1].

S402: according to described weight factor and described observation data determining section plane of vision.In a particular embodiment, part plane of vision is α d _obs.

S403: the solution produced according to the iterative initial value in data space territory, weight factor and interpolation denoising iterative equation determination data space territory iteration.

S404: the solution that described data space territory iteration produces is projected on described part plane of vision.

Namely ${\tilde{d}}_k={\mathrm{\&alpha;d}}_{\mathrm{obs}}+(I-\mathrm{\&alpha;R})d_{k-1}.$

S405: adopt threshold strategies to remove noise respective components at bent wave zone, obtain the solution x in model space territory _k.

Namely $x_k=T_{{\mathrm{\&lambda;}}_k}\left(C({\tilde{d}}_k+(1-\mathrm{\&alpha;})(d_{\mathrm{obs}}-{\mathrm{Rd}}_k))\right)$

S406: the solution in model space territory is projected to data space territory, obtains the seismic data volume after interpolation denoising.

In a particular embodiment, the solution in model space territory is projected to data space territory d _k=C ^tx _k, and make k=k+1, if k<N, then turn back to step S404, circulation performs; Otherwise, return reconstructed results d _k, the solution d in data space territory _kbe exactly interpolation, the complete seismic data volume of rule that obtains of denoising integrated treatment.

The relative convex set projection of this invention (POCS) method or weighting convex set projection (POCS) method, which increase constraint and the threshold denoising step of data residual error, both the object of interpolation can have been played, can also random noise be removed, finally obtain the partial data body of high s/n ratio.α is a weight factor, the speed of control convergence; During α=0, this invention deteriorates to iteration hard thresholding method but it upgrades the solution in data space territory; During α=1, it deteriorates to POCS method+Threshold Denoising Method.Therefore, this invention both can process the interpolation problem of data, also can the Denoising Problems of processing seismic data, and most important advantage is exactly that this invention can carry out interpolation, denoising integrated treatment simultaneously, obtains the data reconstruction of high s/n ratio.

As implied above, be the interpolation denoising method of a kind of geological data provided by the invention, owing to taking above innovative solution, it has the following advantages: 1, to without to make an uproar data, this invention has the efficient speed of convergence of POCS method and perfectly rebuilds effect; 2, to noisy data, it, in conjunction with iteration hard-threshold strategy, adds threshold denoising step in data reconstruction processes, and interpolation, denoising hocket, and finally obtains the reconstructed results of high s/n ratio, processes while achieving geological data interpolation, denoising.

The structured flowchart of the interpolation denoising system of a kind of geological data that Fig. 5 provides for the embodiment of the present invention, as shown in Figure 5, described system comprises:

Sampling matrix acquisition device 101, for obtaining the sampling matrix preset.In the specific embodiment of the present invention, sampling matrix R represents.

Observation data harvester 102, for gathering observation data according to described sampling matrix.Fundamental purpose of the present invention is that the observation data gathered according to this step constructs complete geological data, the seismic data volume namely after interpolation denoising.

In the specific embodiment of the present invention, observation data d _obsrepresent, complete geological data d ₀represent.Observation data d _obswith partial data d ₀and the relation of random noise n is as shown in formula (1):

d _obs＝Rd ₀+n (1)

Figure 10 is the schematic diagram of simulated earthquake data in specific embodiment provided by the invention, Figure 11 is the schematic diagram of missing data in specific embodiment provided by the invention, Figure 12 is the schematic diagram of noisy simulated data in specific embodiment provided by the invention, Figure 13 is noisy missing data schematic diagram corresponding in specific embodiment provided by the invention, Figure 14 is the schematic diagram of actual oceanographic data in specific embodiment provided by the invention, and Figure 15 is disappearance oceanographic data schematic diagram noisy in specific embodiment provided by the invention.

Cost functional construction device 103, for establishing target functional.

Due to the impact of shortage of data and noise, solving of equation (1) is ill posed, and consider openness in warp wavelet territory of geological data, adopt sparse promotion strategy, establishing target functional is as follows,

$\mathrm{\&Phi;}\left(x\right)={\left|\right|d_{\mathrm{obs}}-{\mathrm{RC}}^{T}x\left|\right|}_2^2+\mathrm{\&lambda;P}\left(x\right)---\left(2\right)$

Wherein, Φ (x) is cost functional, and x is bent wave system number vector, C ^tfor bent ripple inverse transformation, C is warp wavelet, and λ is regularization factors, P (x) for sparse constraint (|| x|| ₀or || x|| ₁), d _obsfor observation data, R is sampling matrix.

Seismic data volume determining device 104, for determining the seismic data volume after the interpolation denoising that described observation data is corresponding according to described observation data and described cost functional.Fig. 6 is the structured flowchart of seismic data volume determining device 400, and as shown in Figure 6, seismic data volume determining device 104 specifically comprises:

Threshold determination module 201, for according to described observation data definite threshold.Fig. 7 is the structured flowchart of threshold determination module 201, and as shown in Figure 7, threshold determination module 201 specifically comprises:

Converter unit 301, for adopting sparse transformation warp wavelet that described observation data is transformed to bent wave zone, namely utilizes sparse transformation warp wavelet, observation data is transformed to bent wave zone.

Bent wave system number acquiring unit 302, transforms to the maximum bent wave system number of the amplitude after bent wave zone for obtaining described observation data, in a particular embodiment, selects the bent wave system number C that amplitude is maximum _maxas a reference.

Acquiring unit 303, for obtaining the weight limit factor, the minimal weight factor and the maximum iteration time that preset, in a particular embodiment, the weight limit factor is p ₁, the minimal weight factor is p ₂, maximum iteration time is N.

Max-thresholds determining unit 304, for according to the described weight limit factor and the maximum bent wave system number determination max-thresholds of amplitude, i.e. max-thresholds τ _max=p ₁c _max.

Minimum threshold determining unit 305, for according to the described minimal weight factor and the maximum bent wave system number determination minimum threshold of amplitude, i.e. minimum threshold τ _min=p ₂c _max.

Threshold value determination unit 306, for according to described minimum threshold, max-thresholds, maximum iteration time and exponential function definite threshold.

In a particular embodiment, given maximum iteration time N, then threshold value can be determined by exponential function, concrete formula as shown in formula (3),

τ _k＝τ _maxe ^c(k-1)/(N-1)(3)

c＝ln(τ _min/τ _max)

Wherein, N is maximum iteration time, k=1,2 ..., N, τ _kfor the threshold value of kth time iteration, τ _maxfor max-thresholds, τ _minfor minimum threshold, c is intermediate variable.

As k=N, then τ _n=τ _min=p ₂c _max.

As k=1, then τ ₁=τ _max=p ₁c _max.

Fig. 6 is known, and seismic data volume determining device 104 also comprises:

Hard threshold function determination module 202, for according to described threshold value and described cost functional determination hard threshold function.

After threshold value is determined in a particular embodiment, when sparse constraint P (x) in cost functional=|| x|| ₀namely sparse constraint P (x) is L ₀during norm, hard threshold function T _λx () is adopted, it meets

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x,& \left|x\right|\&GreaterEqual;\mathrm{\&tau;}\\ 0,& \left|x\right|<\mathrm{\&tau;}\end{array}\right.,---\left(4\right)$

Wherein, for threshold value.

Soft-threshold function determination module 203, for according to described threshold value and described cost functional determination soft-threshold function.

In a particular embodiment, when sparse constraint P (x) in cost functional=|| x|| ₁namely sparse constraint P (x) is L ₁during norm, soft-threshold function T _λx () meets

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&tau;},& x>\mathrm{\&tau;}\\ 0,& \left|x\right|\&le;\mathrm{\&tau;}\\ x+\mathrm{\&tau;},& x<-\mathrm{\&tau;}\end{array}\right.---\left(5\right)$

Wherein, τ=0.5 λ is threshold value.

Iterative equation builds module 204, for building interpolation denoising iterative equation.Interpolation denoising iterative equation is as shown in formula (6):

$\begin{array}{c}{\tilde{d}}_k=\mathrm{\&alpha;}{d}_{\mathrm{obs}}+(I-\mathrm{\&alpha;R})d_{k-1}\\ x_k=T_{{\mathrm{\&lambda;}}_k}\left(C({\tilde{d}}_k+(1-\mathrm{\&alpha;})(d_{\mathrm{obs}}-{\mathrm{Rd}}_k))\right)\end{array}.---\left(6\right)$

d _k＝C ^Tx _k

Wherein, d _k=C ^tx _kthe solution in data space, it is the solution after part observation data is inserted.

Seismic data volume determination module 205, for solving described cost functional according to described observation data, hard threshold function, soft-threshold function and interpolation denoising iterative equation, obtains the seismic data volume after interpolation denoising.Fig. 8 is the structured flowchart of seismic data volume determination module 205, and as shown in Figure 8, seismic data volume determination module 205 specifically comprises:

Iterative initial value determining unit 401, for determining iterative initial value and the weight factor in data space territory, in a particular embodiment, makes k=1, selects the iterative initial value d in data space territory ₀=d _obs, weight factor α ∈ [0,1].

Part plane of vision determining unit 402, for according to described weight factor and described observation data determining section plane of vision.In a particular embodiment, part plane of vision is α d _obs.

Determining unit 403 is separated in data space territory, for the solution that the iterative initial value according to data space territory, weight factor and interpolation denoising iterative equation determination data space territory iteration produce.

Projecting cell 404, projects to described part plane of vision for the solution described data space territory iteration produced.

Namely ${\tilde{d}}_k={\mathrm{\&alpha;d}}_{\mathrm{obs}}+(I-\mathrm{\&alpha;R})d_{k-1}.$

Determining unit 405 is separated in model space territory, for adopting threshold strategies to remove noise respective components at bent wave zone, obtains the solution x in model space territory _k.

Namely $x_k=T_{{\mathrm{\&lambda;}}_k}\left(C({\tilde{d}}_k+(1-\mathrm{\&alpha;})(d_{\mathrm{obs}}-{\mathrm{Rd}}_k))\right)$

Seismic data volume determining unit 406, for the solution in model space territory is projected to data space territory, obtains the seismic data volume after interpolation denoising.

In a particular embodiment, the solution in model space territory is projected to data space territory d _k=C ^tx _k, and make k=k+1, if k<N, then turn back to projecting cell 404, circulation performs; Otherwise, return reconstructed results d _k, the solution d in data space territory _kbe exactly interpolation, the complete seismic data volume of rule that obtains of denoising integrated treatment.

The relative convex set projection of this invention (POCS) method or weighting convex set projection (POCS) method, which increase constraint and the threshold denoising step of data residual error, both the object of interpolation can have been played, can also random noise be removed, finally obtain the partial data body of high s/n ratio.α is a weight factor, the speed of control convergence; During α=0, this invention deteriorates to iteration hard thresholding method but it upgrades the solution in data space territory; During α=1, it deteriorates to POCS method+Threshold Denoising Method.Therefore, this invention both can process the interpolation problem of data, also can the Denoising Problems of processing seismic data, and most important advantage is exactly that this invention can carry out interpolation, denoising integrated treatment simultaneously, obtains the data reconstruction of high s/n ratio.

As implied above, be the interpolation denoising system of a kind of geological data provided by the invention, owing to taking above innovative solution, it has the following advantages: 1, to without to make an uproar data, this invention has the efficient speed of convergence of POCS method and perfectly rebuilds effect; 2, to noisy data, it, in conjunction with iteration hard-threshold strategy, adds threshold denoising step in data reconstruction processes, and interpolation, denoising hocket, and finally obtains the reconstructed results of high s/n ratio, processes while achieving geological data interpolation, denoising.

Below in conjunction with specific embodiment, introduce technical scheme of the present invention in detail.Fig. 9 is the process flow diagram of geological data interpolation denoising integral method in specific embodiment provided by the invention, and for realizing above object, this method takes technical scheme as shown in Figure 9, specific as follows:

(1) design sampling matrix R, collect observation data d according to sampling matrix R _obs;

(2) according to observation data d _obscorresponding bent wave system number, determines maximum, minimum threshold, according to maximum iteration time N, by index threshold function definite threshold;

(3) determine to carry out interpolation denoising iterative equation, select iterative initial value and the weight factor in data space territory;

(4) solution that data space territory iteration produces is projected on part plane of vision;

(5) remove noise respective components in warp wavelet territory by threshold strategies, obtain the solution in model space territory;

(6) solution in model space territory is projected to data space territory, the above step of iteration (4), (5) and (6), to maximum iteration time, obtain the complete seismic data volume of rule after interpolation denoising.

Utilize this invention and weighting POCS method to carry out interpolation to the missing data in Figure 11, maximum iteration time is 50 times, and sparse constraint P (x) elects L as ₀norm constraint, threshold function table elects hard threshold function as, and the signal to noise ratio (S/N ratio) curve of reconstruct is as shown in Figure 16,17.Figure 16 is the signal to noise ratio (S/N ratio) curve synoptic diagram that weighting POCS method is rebuild, and indicates the drawback of weighting POCS method when weight factor α is less, because it lacks the constraint of data parameters.Figure 17 is the signal to noise ratio (S/N ratio) curve synoptic diagram rebuild in specific embodiment provided by the invention, and demonstrate the validity of this invention, weight factor α only affects speed of convergence, on last reconstruction SNR without impact.

Figure 12 is the schematic diagram of noisy simulated data in specific embodiment provided by the invention; Figure 13 is noisy missing data schematic diagram corresponding in specific embodiment provided by the invention, POCS, weighting POCS (α=0.6) are carried out to noisy simulated data and noisy real data and the present invention (α=0.6) processes, rebuild effect and residual error as shown in Figure 18-29, wherein residual error is defined as reconstruct data and original nothing and makes an uproar the difference of data.Wherein, Figure 18 is POCS method noisy simulated data reconstructed results schematic diagram, Figure 19 is residual error schematic diagram corresponding to POCS method, Figure 20 is weighting POCS method noisy simulated data reconstructed results schematic diagram, Figure 21 is residual error schematic diagram corresponding to weighting POCS method, Figure 22 is noisy simulated data reconstructed results schematic diagram in specific embodiment provided by the invention, Figure 23 is residual error schematic diagram corresponding in specific embodiment provided by the invention, Figure 24 is POCS method noisy real data reconstructed results schematic diagram, Figure 25 is residual error schematic diagram corresponding to POCS method, Figure 26 is weighting POCS method noisy real data reconstructed results schematic diagram, Figure 27 is residual error schematic diagram corresponding to weighting POCS method, Figure 28 is noisy real data reconstructed results schematic diagram in specific embodiment provided by the invention, Figure 29 is residual error schematic diagram corresponding in specific embodiment provided by the invention.To noisy simulated data, POCS, weighting POCS and reconstruction SNR of the present invention are respectively 8.1, and 11.6 and 16.6dB; To noisy real data, POCS, weighting POCS and reconstruction SNR of the present invention are respectively 6.5, and 9.5 and 12.9dB.As can be seen from Figure 18-29 and reconstruction SNR, relative POCS method and weighting POCS method, the present invention is with the obvious advantage, and the signal to noise ratio (S/N ratio) of reconstruction result is higher, serves the integrated effect of interpolation, denoising.

In sum, the invention provides a kind of interpolation denoising method and system of geological data, farthest rebuild missing data, remove the random noise in section simultaneously, finally obtain the complete regular geological data of high s/n ratio.The present invention is after have studied geological data interpolation problem and seismic data noise attenuation problem, and the basis in conjunction with interpolation, Denoising Problems common feature proposes, and interpolation problem only considers data reconstruction, often affected by noise, makes reconstruction effect undesirable; Shortage of data generally can affect the effect of denoising, therefore the present invention is directed to the disappearance geological data of low signal-to-noise ratio, propose the method for interpolation, denoising integration, denoising is carried out by threshold strategies after each interpolation processing, obtain the complete regular data body of high s/n ratio after iteration terminates, relatively conventional interpolation method, the method has stronger noise immunity, the signal to noise ratio (S/N ratio) of reconstructed results is high, reaches the advantage that interpolation denoising processes simultaneously.

One of ordinary skill in the art will appreciate that all or part of flow process realized in above-described embodiment method, the hardware that can carry out instruction relevant by computer program has come, described program can be stored in general computer read/write memory medium, this program, when performing, can comprise the flow process of the embodiment as above-mentioned each side method.Wherein, described storage medium can be magnetic disc, CD, read-only store-memory body (Read-Only Memory, ROM) or random store-memory body (Random AccessMemory, RAM) etc.

Those skilled in the art can also recognize that the various functions that the embodiment of the present invention is listed are the designing requirements realizing depending on specific application and whole system by hardware or software.Those skilled in the art for often kind of specifically application, can use the function described in the realization of various method, but this realization can should not be understood to the scope exceeding embodiment of the present invention protection.

Apply specific embodiment in the present invention to set forth principle of the present invention and embodiment, the explanation of above embodiment just understands method of the present invention and core concept thereof for helping; Meanwhile, for one of ordinary skill in the art, according to thought of the present invention, all will change in specific embodiments and applications, in sum, this description should not be construed as limitation of the present invention.

Claims (12)

1. an interpolation denoising method for geological data, is characterized in that, described method comprises:

Obtain the sampling matrix preset;

Observation data is gathered according to described sampling matrix;

Establishing target functional;

The seismic data volume after the interpolation denoising that described observation data is corresponding is determined according to described observation data and described cost functional.

2. method according to claim 1, is characterized in that, the cost functional of structure is:

$\mathrm{\&Phi;}\left(x\right)={\left|\right|d_{\mathrm{obs}}-{\mathrm{RC}}^{T}x\left|\right|}_2^2+\mathrm{\&lambda;P}\left(x\right)$

Wherein, Φ (x) is cost functional, and x is bent wave system number vector, C ^tfor bent ripple inverse transformation, C is warp wavelet, and λ is regularization factors, and P (x) is sparse constraint, d _obsfor observation data, R is sampling matrix.

3. method according to claim 2, is characterized in that, the seismic data volume after determining according to described observation data and described cost functional the interpolation denoising that described observation data is corresponding comprises:

According to described observation data definite threshold;

According to described threshold value and described cost functional determination hard threshold function;

According to described threshold value and described cost functional determination soft-threshold function;

Build interpolation denoising iterative equation;

Solve described cost functional according to described observation data, hard threshold function, soft-threshold function and interpolation denoising iterative equation, obtain the seismic data volume after interpolation denoising.

4. method according to claim 3, is characterized in that, comprises according to described observation data definite threshold:

Adopt sparse transformation warp wavelet that described observation data is transformed to bent wave zone;

Observation data described in acquisition transforms to the maximum bent wave system number of the amplitude after bent wave zone;

Obtain the weight limit factor, the minimal weight factor and the maximum iteration time that preset;

According to the described weight limit factor and the maximum bent wave system number determination max-thresholds of amplitude;

According to the described minimal weight factor and the maximum bent wave system number determination minimum threshold of amplitude;

According to described minimum threshold, max-thresholds, maximum iteration time and exponential function definite threshold.

5. method according to claim 4, is characterized in that, described exponential function is:

τ _k＝τ _maxe ^c(k-1)/(N-1)

c＝ln(τ _min/τ _max)

Wherein, N is maximum iteration time, k=1,2 ..., N, τ _kfor the threshold value of kth time iteration, τ _maxfor max-thresholds, τ _minfor minimum threshold, c is intermediate variable.

6. method according to claim 4, is characterized in that, solve described cost functional according to described observation data, hard-threshold, soft-threshold and interpolation denoising iterative equation, the seismic data volume obtained after interpolation denoising comprises:

Determine iterative initial value and the weight factor in data space territory;

According to described weight factor and described observation data determining section plane of vision;

According to the solution that the iterative initial value in data space territory, weight factor and interpolation denoising iterative equation determination data space territory iteration produce;

The solution that described data space territory iteration produces is projected on described part plane of vision;

Adopt threshold strategies to remove noise respective components at bent wave zone, obtain the solution in model space territory;

Solution in model space territory is projected to data space territory, obtains the seismic data volume after interpolation denoising.

7. an interpolation denoising system for geological data, is characterized in that, described system comprises:

Sampling matrix acquisition device, for obtaining the sampling matrix preset;

Observation data harvester, for the observation data that the sampling matrix acquiring seismic data according to described is corresponding;

Cost functional construction device, for establishing target functional;

Seismic data volume determining device, for determining the seismic data volume after the interpolation denoising that described observation data is corresponding according to described observation data and described cost functional.

8. system according to claim 7, is characterized in that, the cost functional of structure is:

$\mathrm{\&Phi;}\left(x\right)={\left|\right|d_{\mathrm{obs}}-{\mathrm{RC}}^{T}x\left|\right|}_2^2+\mathrm{\&lambda;P}\left(x\right)$

Wherein, Φ (x) is cost functional, and x is bent wave system number vector, C ^tfor bent ripple inverse transformation, C is warp wavelet, and λ is regularization factors, and P (x) is sparse constraint, d _obsfor observation data, R is sampling matrix.

9. system according to claim 8, is characterized in that, described seismic data volume determining device comprises:

Threshold determination module, for according to described observation data definite threshold;

Hard threshold function determination module, for according to described threshold value and described cost functional determination hard threshold function;

Soft-threshold function determination module, for according to described threshold value and described cost functional determination soft-threshold function;

Iterative equation builds module, for building interpolation denoising iterative equation;

Seismic data volume determination module, for solving described cost functional according to described observation data, hard threshold function, soft-threshold function and interpolation denoising iterative equation, obtains the seismic data volume after interpolation denoising.

10. system according to claim 9, is characterized in that, described threshold determination module comprises:

Converter unit, transforms to bent wave zone for adopting sparse transformation warp wavelet by described observation data;

Bent wave system number acquiring unit, transforms to the maximum bent wave system number of the amplitude after bent wave zone for obtaining described observation data;

Acquiring unit, for obtaining the weight limit factor, the minimal weight factor and the maximum iteration time that preset;

Max-thresholds determining unit, for according to the described weight limit factor and the maximum bent wave system number determination max-thresholds of amplitude;

Minimum threshold determining unit, for according to the described minimal weight factor and the maximum bent wave system number determination minimum threshold of amplitude;

Threshold value determination unit, for according to described minimum threshold, max-thresholds, maximum iteration time and exponential function definite threshold.

11. systems according to claim 10, is characterized in that, described exponential function is:

τ _k＝τ _maxe ^c(k-1)/(N-1)

c＝ln(τ _min/τ _max)

Wherein, N is maximum iteration time, k=1,2 ..., N, τ _kfor the threshold value of kth time iteration, τ _maxfor max-thresholds, τ _minfor minimum threshold, c is intermediate variable.

12. systems according to claim 10, is characterized in that, described seismic data volume determination module comprises:

Iterative initial value determining unit, for determining iterative initial value and the weight factor in data space territory;

Part plane of vision determining unit, for according to described weight factor and described observation data determining section plane of vision;

Determining unit is separated in data space territory, for the solution that the iterative initial value according to data space territory, weight factor and interpolation denoising iterative equation determination data space territory iteration produce;

Projecting cell, projects to described part plane of vision for the solution described data space territory iteration produced;

Determining unit is separated in model space territory, for adopting threshold strategies to remove noise respective components at bent wave zone, obtains the solution in model space territory;

Seismic data volume determining unit, for the solution in model space territory is projected to data space territory, obtains the seismic data volume after interpolation denoising.